diff --git a/doc/matlab_macro_rugo/find_Q_nat.m b/doc/matlab_macro_rugo/find_Q_nat.m
index feea83b7aab92e56295666f2a3ff89304f473fbd..9e786d2e949ac83842fbb2fc44539af96288bcd8 100755
--- a/doc/matlab_macro_rugo/find_Q_nat.m
+++ b/doc/matlab_macro_rugo/find_Q_nat.m
@@ -14,13 +14,10 @@ Frg=Fr/(1-C^0.5);
if ks==0;ks=1e-34;end
-
-coeff_contraction=0.4*Cd0+0.7;
-
%if Cd0==2
% fFr=(min(2.5,Frg^(-4/3)));%
%else
- fFr=(min(coeff_contraction./(1-(Frg.^2)/4),Frg.^(-2/3))).^2;
+ fFr=(min(1/(1-(Frg.^2)/4),Frg.^(-2/3))).^2;
%end
if Frg>1.5
@@ -37,76 +34,76 @@ R=(1-sigma*C);%.*(1-C.^0.5).^2;
if pf/h>1; %fFr=1;
-
+
choixturb=1;
htilde=h./D;
hstar=pf./h;
Rh=h.*(hstar-1);
ustar=(g.*S.*Rh).^0.5;
-
+
Cd1=Cd;
CdCh=Cd1.*C.*htilde;
-
+
Cf=2./(5.1.*log10(h./ks)+6).^2;
U0b=(2*g.*S.*R./(Cd1.*C.*h./D+Cf.*R).*h).^0.5;
-
+
% U0=(2*g.*S.*D.*R./(Cd1.*C)).^0.5;
-
+
[P]=fminbnd(@(alphai) resolve_alphaR(alphai,CdCh,R,U0b,hstar,h,C,D,Cd1,ustar,choixturb),1e-5*h,h,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx));
-
+
alpha=P(1);
-
+
beta2=h.*CdCh./alpha./R;
beta=(beta2).^0.5;
a1=beta*(hstar-1)/(cosh(beta));
c=1;
-
+
UhU0=(a1*sinh(beta)+c)^0.5;
Uh=UhU0*U0b;
-
+
dhp=1-1/kappa*alpha./h.*Uh./ustar;
z0hp=(1-dhp).*exp(-1*(kappa*Uh./ustar));
-
+
qsup=ustar./kappa.*h.*((hstar-dhp).*(log((hstar-dhp)./z0hp) - 1)-((1-dhp).*(log((1-dhp)./z0hp) - 1)));
-
- %calcul intégrale dans la canopée----
+
+ %calcul int�grale dans la canop�e----
dzinf=0.01;
Zinf=(0:dzinf:1);
Uinf=U0b .*(beta.*Rh./h.*sinh(beta*Zinf)./cosh(beta)+1).^0.5;
Ub=zeros(size(Uinf));
- Ub(1:end-1)=Uinf(2:end);
+ Ub(1:end-1)=Uinf(2:end);
qinf=sum((Uinf(1:end-1)+Ub(1:end-1))/2*dzinf.*h);
-
+
qtot=qinf+qsup;
-
+
PI=0;
delta=1;
-
+
Umoy=qtot./pf;
res=abs(U0-Umoy);
-
+
else
-
+
hstar=pf/D;
Re=U0.*pf/1e-6;
-
+
Cd=Cd.*fFr;
-
+
if ks==0
- Cf=0.3164/4.*Re.^(-0.25);
+ Cf=0.3164/4.*Re.^(-0.25);
else
- Cf=2/(5.1*log10(pf/ks)+6)^2;
- end
-
+ Cf=2/(5.1*log10(pf/ks)+6)^2;
+ end
+
N= (1.*Cf)./(pf./D.*Cd.*C);
U0a=(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5;
-
-
-
-[u]=fminsearch(@(U0i) find_U0_complet(U0i,pf,C,D,sigma,Cd0,Cf,coeff_contraction,S),U0a,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx));
+
+
+
+[u]=fminsearch(@(U0i) find_U0_complet(U0i,pf,C,D,sigma,Cd0,Cf,S),U0a,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx));
res=abs(U0-u);
-
+
end
diff --git a/doc/matlab_macro_rugo/find_U0_complet.m b/doc/matlab_macro_rugo/find_U0_complet.m
index 14c033afee2934155fa70bda7f7765aba003e2fb..932f559fec327c12c689abb20657c84d14c45b3a 100755
--- a/doc/matlab_macro_rugo/find_U0_complet.m
+++ b/doc/matlab_macro_rugo/find_U0_complet.m
@@ -1,12 +1,12 @@
-function [res]=find_U0_complet(U01,h,C,D,sigma,Cd0,cfmean,r,S)
+function [res]=find_U0_complet(U01,h,C,D,sigma,Cd0,cfmean,S)
g=9.81;
Frexp=U01./(1-(1.*C).^0.5)./(9.81*h)^0.5;
%Frexp=U01./(9.81*h)^0.5;
- fFr=(min(r./(1-(Frexp.^2)/4),Frexp.^(-2/3))).^2;
-
+ fFr=(min(1./(1-(Frexp.^2)/4),Frexp.^(-2/3))).^2;
+
+
-
if Frexp>=1.5
fFr=(Frexp.^(-2/3)).^2;
end
diff --git a/doc/matlab_macro_rugo/rampe_rugueuse.m b/doc/matlab_macro_rugo/rampe_rugueuse.m
index 051092efb3f98a9153d0e13b0531832750a0ceca..b3828a76fe3b2fd0dc9dab08d54fc86ba5d036ac 100755
--- a/doc/matlab_macro_rugo/rampe_rugueuse.m
+++ b/doc/matlab_macro_rugo/rampe_rugueuse.m
@@ -595,14 +595,10 @@ set(handles.edit11,'string',num2str(cote_bas))
Vdeb=Q/L/pf;
set(handles.edit12,'string',num2str(Vdeb))
-coeff_contraction=0.4*Cd+0.7;
Vg=Q/L/pf/(1-C^0.5);
Fr=Vg./(g.*pf).^0.5;
-if Cd==2
- Vmax=Vg.*(min(coeff_contraction./(1-(Fr.^2)/4),Fr.^(-2/3)));%
-else
- Vmax=Vg.*(min(coeff_contraction./(1-(Fr.^2)/4),Fr.^(-2/3)));
-end
+coeff_contraction=0.2*Cd+0.9;
+Vmax=Vg.*(min(coeff_contraction./(1-(Fr.^2)/4),Fr.^(-2/3)));
set(handles.edit16,'string',num2str(Vmax))
diff --git a/doc/scilab_tests/macrorugo_find_Q_nat.sci b/doc/scilab_tests/macrorugo_find_Q_nat.sci
index a80e5f88b850202a8dfedd0e6d88fd4982a99e38..a3d9462839ca73255ccae8a7adbbd470590f7958 100755
--- a/doc/scilab_tests/macrorugo_find_Q_nat.sci
+++ b/doc/scilab_tests/macrorugo_find_Q_nat.sci
@@ -18,12 +18,10 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
if ks==0;ks=1e-34;end
- coeff_contraction=0.4*Cd0+0.7;
-
//if Cd0==2
// fFr=(min(2.5,Frg^(-4/3)));//
//else
- fFr=(min(coeff_contraction./(1-(Frg.^2)/4),abs(Frg).^(-2/3))).^2;
+ fFr=(min(1./(1-(Frg.^2)/4),abs(Frg).^(-2/3))).^2;
//end
if abs(Frg)>1.5
@@ -121,10 +119,10 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
U0a=(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5;
if bDbg then print_r("U0a");end;
- U0complet = find_U0_complet(U0a, pf,C,D,sigma,Cd0,Cf,coeff_contraction,S, bDbg)
+ U0complet = find_U0_complet(U0a, pf,C,D,sigma,Cd0,Cf,S, bDbg)
if bDbg then print_r("U0complet");end;
- [u fval] = fminsearch(list(find_U0_complet, pf,C,D,sigma,Cd0,Cf,coeff_contraction,S), U0a, opt)
+ [u fval] = fminsearch(list(find_U0_complet, pf,C,D,sigma,Cd0,Cf,S), U0a, opt)
res=abs(U0-u);
diff --git a/doc/scilab_tests/macrorugo_find_U0_complet.sci b/doc/scilab_tests/macrorugo_find_U0_complet.sci
index bd01b5b87cca7d912d0a4989090e13c78c836e67..d819d5b68ed935e1f44e84ef153d8653d816b636 100644
--- a/doc/scilab_tests/macrorugo_find_U0_complet.sci
+++ b/doc/scilab_tests/macrorugo_find_U0_complet.sci
@@ -1,11 +1,11 @@
-function [res]=find_U0_complet(U01,h,C,D,sigma,Cd0,cfmean,r,S, bDbgU0)
+function [res]=find_U0_complet(U01,h,C,D,sigma,Cd0,cfmean,S, bDbgU0)
if ~exists("bDbgU0","local"); bDbgU0 = %f; end;
g=9.81;
Frexp=U01./(1-(1.*C).^0.5)./(9.81*h)^0.5;
//Frexp=U01./(9.81*h)^0.5;
- fFr=(min(r./(1-(Frexp.^2)/4),Frexp.^(-2/3))).^2;
+ fFr=(min(1./(1-(Frexp.^2)/4),Frexp.^(-2/3))).^2;
if Frexp>=1.5
fFr=(Frexp.^(-2/3)).^2;
diff --git a/src/macrorugo/macrorugo.ts b/src/macrorugo/macrorugo.ts
index 1bde85f2050be3ba35b1825709c74fe05be7fb22..ff93fd5c35a9d1a5b0a3910bcf3d773d6fb26b88 100644
--- a/src/macrorugo/macrorugo.ts
+++ b/src/macrorugo/macrorugo.ts
@@ -116,15 +116,7 @@ export class MacroRugo extends FishPass {
}
r.resultElement.values.Fr = resFr;
// Vitesse maximale
- const cc = 0.4 * this.prms.Cd0.v + 0.7;
- let resVmax = vg * Math.min(
- cc / (1 - Math.pow(resFr, 2) / 4),
- Math.pow(resFr, -2 / 3)
- );
- if (isNaN(resVmax)) { // if Y == 0
- resVmax = 0;
- }
- r.resultElement.values.Vmax = resVmax;
+ r.resultElement.values.Vmax = this.r * vg * this.fFr;
// Puissance dissipée
r.resultElement.values.PV = 1000 * MacroRugo.g * r.resultElement.values.Vdeb * this.prms.If.v;
// Type d'écoulement
@@ -440,6 +432,28 @@ export class MacroRugo extends FishPass {
);
}
+ /**
+ * Ratio entre la vitesse moyenne à l'aval d'un block et la vitesse maximale
+ */
+ private get r(): number {
+ if (this._cache.r === undefined) {
+ this._cache.r = this.Calc_r();
+ }
+ return this._cache.r;
+ }
+
+ /**
+ * Calcul du ratio entre la vitesse moyenne à l'aval d'un block et la vitesse maximale
+ * r = 1.1 pour un plot circulaire Cd0​=1 et r = 1.5 pour un plot à face plane Cd0​=2.5
+ * Voir #283
+ */
+ private Calc_r(): number {
+ return 0.2 * this.prms.Cd0.v + 0.9;
+ }
+
+ /**
+ * Froude correction function (Cassan et al. 2014, Eq. 19)
+ */
private get fFr(): number {
if (this._cache.fFr === undefined) {
this._cache.fFr = this.CalcfFr(this.U0);
@@ -448,7 +462,7 @@ export class MacroRugo extends FishPass {
}
/**
- * Froude correction function
+ * Calculation of Froude correction function (Cassan et al. 2014, Eq. 19)
*/
private CalcfFr(U0: number): number {
// tslint:disable-next-line:variable-name
@@ -456,11 +470,8 @@ export class MacroRugo extends FishPass {
(1 - Math.sqrt(MacroRugo.fracAxAy * this.prms.C.v)) /
Math.sqrt(MacroRugo.g * this.prms.Y.v);
- /** Interpolation linéaire entre le bloc rond (Cd0=1) et le carré (Cd0=2) */
- const r = 0.4 * this.prms.Cd0.v + 0.7;
-
if (Fr < 1.3) {
- return Math.pow(Math.min(r / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2);
+ return Math.pow(Math.min(1 / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2);
} else {
return Math.pow(Fr, -4 / 3);
}