feat(macrorugo): new formulas for fh*(h*), fF(F) and r

- All the tests crash due to different calculation results!

Refs #283

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 intgrale dans la canope----
+
+        %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
 

doc/matlab_macro_rugo/find_U0_complet.m

-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

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))
 

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);
 

doc/scilab_tests/macrorugo_find_U0_complet.sci

-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;

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);
         }