Commit dbae67af by Dorchies David

### #22 Passe tous les tests en émergent sauf Vmax

parent de7a95a7
 function [res]=find_Q_nat2(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg) //fonctin pour ax=ay maxfun=5000; maxiter=5000; tolfun=1e-16; tolx=1e-16; g=9.81; kappa=0.41; U0=Q./L./pf; if(bDbg) then print_r("U0"); end; Fr=U0./(9.81*pf).^0.5; Frg=Fr/(1-C^0.5); 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; //end if Frg>1.5 fFr=(Frg.^(-2/3)).^2; end if(bDbg) then print_r("fFr"); end; alpha=1-(1.*C).^0.5-1/2*sigma.*C; //Cd=Cd0.*(1+1./(pf./D).^2).*fFr; Cd=Cd0.*(0.8-2*C).*(1+0.4./(pf./D).^2).*fFr; R=(1-sigma*C).*(1-C.^0.5).^2; if pf/h>1.1; fFr=1; choixturb=1; htilde=h./D; hstar=pf./h; Rh=h.*(hstar-1); ustar=(g.*S.*Rh).^0.5; CdCh=Cd0.*(1+0.4./(pf./D).^2).*C.*htilde; u0=(2*g.*S.*D.*R./(Cd0.*C)).^0.5; // [P]=fminbnd(@(alphai) resolve_alpha(alphai,CdCh,R,u0,hstar,h,C,D,Cd0,ustar,choixturb),1e-5*h,h,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); //[P]=fminbnd(@(alphai) resolve_alpha(alphai,CdCh,R,u0,hstar,h,C,D,Cd0,ustar,choixturb),1e-5*h,h); [alpha fval] = fminsearch(list(resolve_alpha, CdCh,R,u0,hstar,h,C,D,Cd0,ustar), 1e-5*h) beta2=h.*CdCh./alpha./R; beta2=(beta2).^0.5; a1=beta2*(hstar-1)/(cosh(beta2)); c=1; UhU0=(a1*sinh(beta2)+c)^0.5; Uh=UhU0*u0; 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---- U(1)=u0; U(2)=u0.*(beta2.*Rh./h.*sinh(beta2*0.1)./cosh(beta2)+c).^0.5; U(3)=u0.*(beta2.*Rh./h.*sinh(beta2*0.2)./cosh(beta2)+c).^0.5; U(4)=u0.*(beta2.*Rh./h.*sinh(beta2*0.3)./cosh(beta2)+c).^0.5; U(5)=u0.*(beta2.*Rh./h.*sinh(beta2*0.4)./cosh(beta2)+c).^0.5; U(6)=u0.*(beta2.*Rh./h.*sinh(beta2*0.5)./cosh(beta2)+c).^0.5; U(7)=u0.*(beta2.*Rh./h.*sinh(beta2*0.6)./cosh(beta2)+c).^0.5; U(8)=u0.*(beta2.*Rh./h.*sinh(beta2*0.7)./cosh(beta2)+c).^0.5; U(9)=u0.*(beta2.*Rh./h.*sinh(beta2*0.8)./cosh(beta2)+c).^0.5; U(10)=u0.*(beta2.*Rh./h.*sinh(beta2*0.9)./cosh(beta2)+c).^0.5; U(11)=Uh; Ub=zeros(U); Ub(1:\$-1)=U(2:\$); qinf=sum((U(1:\$-1)+Ub(1:\$-1))/2*0.1.*h); qtot=qinf+qsup; PI=0.2; delta=1; Umax=ustar./kappa*(log((delta*(h-h)-dhp*h)/(z0hp*h))+2*PI); Umoy=qtot./pf; res=abs(U0-Umoy); else hstar=pf/D; Re=U0.*pf/1e-6; if ks==0 Cf=0.3164/4.*Re.^(-0.25); else Cf=2/(5.1*log10(pf/ks-1)+6)^2; end if(bDbg) then print_r("Cf"); end; //[u res]=fminsearch(@(U0i) find_U0_complet(U0i,pf,C,D,sigma,Cd0,Cf,coeff_contraction,S),U0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); [u res] = fminsearch ( list(find_U0_R0,pf,C,D,sigma,Cd0,Cf,coeff_contraction,S) , U0 ) // N= (alpha.*Cf)./(pf./D.*Cd.*C); // res=abs(U0-(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5); res=abs(U0-u); end endfunction
 function [res]=find_U0_R0(U01,h,C,D,sigma,Cd0,cfmean,r,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; if Frexp>1.3 fFr=(Frexp.^(-2/3)).^2; end alpha=1-(1.*C); Cd=Cd0.*(1+0.4./(h./D).^2).*fFr; //Cd=Cd0.*1.*(1+1./(h./D).^2).*fFr; N= (alpha.*cfmean)./(h./D.*Cd.*C); hstar=h/D; //res=(U01-(2*g.*S.*D.*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5); res=abs(U01-(2*g.*S.*D*(1-sigma.*C)./(Cd.*C.*(1+N))).^0.5); endfunction
 function [res]= resolve_alpha(alpha,CdCh,R,U0,hstar,hp,C,D,Cd,ustar) g=9.81; kappa=0.41; L=D*(1/C^0.5-1); beta2=hp.*CdCh./alpha./R; beta2=(beta2).^0.5; a1=beta2*(hstar-1)/(cosh(beta2)); c=1; UhU0=(a1*sinh(beta2)+c)^0.5; Uh=UhU0*U0; //choix du modele de turbulence //switch choixturb // // case 1 L1=min(L,0.15*hp); //L1=L;//0.15*hp; // case 2 // // L1=L.*hp./(1/0.41*L+hp); // // case 3 // L1=0.15*min(hp,D./C./Cd);//; // case 4 // if L/hp<1 // L1=(1-L/hp).*L; // if L1<0.05*hp && L/hp>0.8;L1=0.05*hp;end // else // L1=0.41*((1-hp./L))*hp; // if L1<0.05*hp && L/hp<1.2;L1=0.05*hp;end // end //end res=abs(alpha*Uh-L1*ustar); endfunction
 function macrorugo_resultComp(z_amont, S, long, Q, L, pf, C, Cd, h, D) g = 9.81 cote_bas=z_amont-S*long; print_r("cote_bas"); Vdeb=Q/L/pf; print_r("Vdeb"); Vg=Q/L/pf/(1-C^0.5); Fr=Vg./(g.*pf).^0.5; print_r("Fr"); if Cd==2 Vmax=Vg.*(min(1.6./(1-(Fr.^2)/4),Fr.^(-2/3)));% else Vmax=Vg.*(min(1.2./(1-(Fr.^2)/4),Fr.^(-2/3))); end print_r("Vmax"); P=1000*g*Q/L*S; print_r("P"); if pf/h<1 flowcond = 'emergent' elseif pf/h<1.1 & pf/h>=1 flowcond = 'quasi emergent' elseif pf/h>1.2 flowcond = 'immerge' end print_r("flowcond"); if pf/h>1.1 q_technique=0.955*(pf/h)^2.282*S^0.466*C^(-0.23)*(9.81*D)^0.5.*h*L; else if Cd==2 q_technique= 0.648*(pf/D)^1.084*S^0.56*C^(-0.456)*(9.81*D)^0.5.*D*L; V_technique=3.35*(pf/D)^0.27*S^0.53*(9.81*D)^0.5; else q_technique=0.815*(pf/D)^1.45*S^0.557*C^(-0.456)*(9.81*D)^0.5.*D*L; V_technique=4.54*(pf/D)^0.32*S^0.56*(9.81*D)^0.5; end print_r("V_technique"); end print_r("q_technique"); endfunction
 // Test des passes à macro rugosité pour JaLHyd clear sCurrentPath = get_absolute_file_path("main_macrorugo.sce"); getd(sCurrentPath); // Tests parameters ks = 0.01 Cd0=1.5 S = 0.05 C = 0.05; D = 0.5; sigma = 1; B = 1 z_amont = 12.5; long = 6; // ***************************************************************************** printf("\n*** Emergent conditions ***\n") // ***************************************************************************** h = 0.6 k = 0.8 // Test de find_U0_R0 Cf=2/(5.1*log10(h/ks-1)+6)^2 coeff_contraction=0.4*Cd0+0.7 res = find_U0_R0(2.58, h, C, D, sigma, Cd0, Cf, coeff_contraction, S) printf("find_U0_R0=%f\n",res); // Test de find_Q_nat2 [Q fVal] = fminsearch(list(find_Q_nat2, ks,D,k,Cd0,S,B,h,C,sigma,%f), 0.1); printf("Qemerg=%f fVal=%f\n",Q, fVal); res = find_Q_nat2(Q,ks,D,k,Cd0,S,B,h,C,sigma,%t) printf("find_Q_nat2=%f\n",res); macrorugo_resultComp(z_amont, S, long, Q, B, h, C, Cd0, k, D) // ***************************************************************************** printf("\n*** Submerged conditions ***\n") // ***************************************************************************** k = 0.6 h = 0.8 [Q fVal] = fminsearch(list(find_Q_nat2, ks,D,k,Cd0,S,B,h,C,sigma,%f), 0.1); printf("Qsubmerg=%f fVal=%f\n",Q, fVal); res = find_Q_nat2(Q,ks,D,k,Cd0,S,B,h,C,sigma,%t) printf("find_Q_nat2=%f\n",res); macrorugo_resultComp(z_amont, S, long, Q, B, h, C, Cd0, k, D)
 function print_r(s) e = eval(s) if typeof(e) == "string" then pat = "s" else pat = "f" end if size(e,1) > 1 | size(e,2) > 1 then printf("%s = [",s); if size(e,2) > 1 then e = e' end printf("%f ",e) printf("]\n") else printf("%s = %"+pat+"\n",s,e); end endfunction
 ... ... @@ -7,17 +7,18 @@ // import { describe, expect, it, xdescribe, xit } from "../mock_jasmine"; import { ParamCalculability, ParamValueMode } from "../../src"; import { MacroRugo, MacrorugoParams } from "../../src/macrorugo/macrorugo"; import { checkResult } from "../test_func"; import { MacroRugo, MacroRugoFlowType, MacrorugoParams } from "../../src/macrorugo/macrorugo"; import { Result } from "../../src/util/result"; import { checkResult, checkPercent } from "../test_func"; function macroRugoInstance(): MacroRugo { function macroRugoInstanceEmergent(): MacroRugo { return new MacroRugo( new MacrorugoParams( 12.5, // ZF1 6, // L 1, // B 0.05, // If 9.418, // Q 1.57, // Q 0.6, // h 0.01, // Ks 0.05, // C ... ... @@ -28,10 +29,18 @@ function macroRugoInstance(): MacroRugo { ); } function testMacroRugo(varTest: string, valRef: number) { function macroRugoInstanceSubmerged(): MacroRugo { const nub: MacroRugo = macroRugoInstanceEmergent(); nub.prms.Y.v = 0.8; nub.prms.PBH.v = 0.6; nub.prms.Q.v = 4.737; return nub; } function testMacroRugo(nubFactory: () => MacroRugo, varTest: string, valRef: number) { describe("Calc(): ", () => { it(`\${varTest} should be \${valRef}`, () => { const nub = macroRugoInstance(); const nub = nubFactory(); nub.prms.Q.v = nub.Calc("Q").vCalc; const p = nub.getParameter(varTest); p.v = undefined; ... ... @@ -43,12 +52,68 @@ function testMacroRugo(varTest: string, valRef: number) { describe("Class MacroRugo: ", () => { const nub = macroRugoInstance(); const vCalc = nub.Calc("Q").vCalc; describe("Emerged conditions", () => { it(`resolveU0(2.58) should be around 0`, () => { const nubit = macroRugoInstanceEmergent(); // tslint:disable-next-line:no-string-literal expect(nubit["resolveU0"](2.58)).toBeCloseTo(0, 1); }); it(`resolveQ(1.547827) should be around 0`, () => { const nubit = macroRugoInstanceEmergent(); // tslint:disable-next-line:no-string-literal expect(nubit["resolveQ"](1.547827)).toBeCloseTo(0, 1); }); it(`Calc("Q") should be around 1.548`, () => { checkResult(macroRugoInstanceEmergent().Calc("Q", 0.1), 1.548, 1); }); it(`Calc("Q", 0.1).extraResults.ZF2 should be around 12.2`, () => { expect(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.ZF2).toBeCloseTo(12.2, 3); }); it(`Calc("Q", 0.1).extraResults.Vdeb should be around 2.579818`, () => { checkPercent(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.Vdeb, 2.579818, 0.03); }); for (const prm of nub.prms) { if ([ParamCalculability.DICHO, ParamCalculability.EQUATION].includes(prm.calculability)) { testMacroRugo(prm.symbol, prm.v); it(`Calc("Q", 0.1).extraResults.V should be around 2.694279`, () => { checkPercent(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.V, 2.694279, 0.03); }); it(`Calc("Q", 0.1).extraResults.Fr should be around 1.369611`, () => { checkPercent(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.Fr, 1.369611, 0.03); }); it(`Calc("Q", 0.1).extraResults.P should be around 759.240352`, () => { checkPercent(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.P, 759.240352, 0.03); }); it(`Calc("Q", 0.1).extraResults.FlowType should be MacroRugoFlowType.EMERGENT`, () => { expect(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.FlowType).toBe(MacroRugoFlowType.EMERGENT); }); it(`Calc("Q", 0.1).extraResults.Q2 should be around 0.868672`, () => { checkPercent(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.Q2, 0.868672, 0.03); }); it(`Calc("Q", 0.1).extraResults.V2 should be around 1.991299`, () => { checkPercent(macroRugoInstanceEmergent().Calc("Q", 0.1).extraResults.V2, 1.991299, 0.03); }); const nub = macroRugoInstanceEmergent(); for (const prm of nub.prms) { if ([ParamCalculability.DICHO, ParamCalculability.EQUATION].includes(prm.calculability)) { testMacroRugo(macroRugoInstanceEmergent, prm.symbol, prm.v); } } } }); describe("Submerged conditions", () => { it(`resolveQ(4.737) should be around 0`, () => { const nubit = macroRugoInstanceSubmerged(); // tslint:disable-next-line:no-string-literal expect(nubit["resolveQ"](4.737)).toBeCloseTo(0, 1); }); }); });
 ... ... @@ -21,6 +21,12 @@ export function equalEpsilon(val1: number, val2: number, epsilon: number = precD return Math.abs(val1 - val2) < epsilon; } export function checkPercent(valTest: number, valRef: number, rPercent: number) { const r: number = Math.abs(valTest - valRef) / valRef; const bounds = valRef * rPercent; expect(r < rPercent).toBeTruthy(`\${valTest} should be between \${valRef - bounds} and \${valRef + bounds}`); } /** * compare 2 objets * @param s message ... ...
 ... ... @@ -2,6 +2,7 @@ import { Nub } from "../nub"; import { ParamCalculability } from "../param/param-definition"; import { Result } from "../util/result"; import { MacrorugoParams } from "./macrorugo_params"; import { ParamValueMode } from "../param/param-value-mode"; export { MacrorugoParams }; ... ... @@ -28,13 +29,18 @@ export class MacroRugo extends Nub { /** Rugosité de fond (m) */ private ks: number; /** Averaged velocity at the bed (m.s-1) */ /** Averaged velocity (m.s-1) */ private U0: number; /** Velocity at the bed (m.s-1) */ private u0: number; private _cache: { [key: string]: number }; constructor(prms: MacrorugoParams, dbg: boolean = false) { super(prms, dbg); this._cache = {}; } /** ... ... @@ -52,12 +58,17 @@ export class MacroRugo extends Nub { * @param rPrec Précision attendue */ public Calc(sVarCalc: string, rInit?: number, rPrec: number = 0.001): Result { /** @todo Voir pour déclarer le paramètre en calcul dans nub */ this.getParameter(sVarCalc).valueMode = ParamValueMode.CALCUL; const r: Result = super.Calc(sVarCalc, rInit, rPrec); // Ajout des résultats complémentaires // Cote de fond aval r.extraResults.ZF2 = this.prms.ZF1.v - this.prms.If.v * this.prms.L.v; // Vitesse débitante r.extraResults.Vdeb = this.V(this.prms.Q) / this.prms.B.v / this.prms.Y.v; // Froude r.extraResults.Fr = r.extraResults.Vdeb / (1 - Math.sqrt(MacroRugo.fracAyAx * this.prms.C.v)) / Math.sqrt(MacroRugo.g * this.prms.Y.v); // Vitesse maximale r.extraResults.V = r.extraResults.Vdeb * this.calc_fFr(r.extraResults.Vdeb); // Puissance dissipée ... ... @@ -65,7 +76,7 @@ export class MacroRugo extends Nub { // Type d'écoulement if (this.prms.Y.v / this.prms.PBH.v < 1) { r.extraResults.FlowType = MacroRugoFlowType.EMERGENT; } else if (this.prms.Y.v / this.prms.PBH.v < 1.1) { } else if (this.prms.Y.v / this.prms.PBH.v < MacroRugo.limitSubmerg) { r.extraResults.FlowType = MacroRugoFlowType.QUASI_EMERGENT; } else { r.extraResults.FlowType = MacroRugoFlowType.IMMERGE; ... ... @@ -73,9 +84,12 @@ export class MacroRugo extends Nub { // Vitesse et débit du guide technique let cQ: [number, number, number, number]; let cV: [number, number, number]; if (this.prms.Y.v / this.prms.PBH.v > 1.1) { let hdk: number; if (this.prms.Y.v / this.prms.PBH.v > MacroRugo.limitSubmerg) { cQ = [0.955, 2.282, 0.466, -0.23]; hdk = this.prms.PBH.v; } else { hdk = this.prms.PBD.v; if (Math.abs(this.prms.Cd0.v - 2) < 0.05) { cQ = [0.648, 1.084, 0.56, -0.456]; cV = [3.35, 0.27, 0.53]; ... ... @@ -84,10 +98,10 @@ export class MacroRugo extends Nub { cV = [4.54, 0.32, 0.56] } } r.extraResults.Q2 = cQ[0] * Math.pow(this.prms.Y.v / this.prms.PBH.v, cQ[1]) * r.extraResults.Q2 = cQ[0] * Math.pow(this.prms.Y.v / hdk, cQ[1]) * Math.pow(this.prms.If.v, cQ[2]) * Math.pow(this.prms.C.v, cQ[3]) * Math.sqrt(MacroRugo.g * this.prms.PBD.v) * this.prms.PBH.v * this.prms.B.v; if (this.prms.Y.v / this.prms.PBH.v <= 1.1) { Math.sqrt(MacroRugo.g * this.prms.PBD.v) * this.prms.PBD.v * this.prms.B.v; if (this.prms.Y.v / this.prms.PBH.v <= MacroRugo.limitSubmerg) { r.extraResults.V2 = cV[0] * Math.pow(this.prms.Y.v / this.prms.PBD.v, cV[1]) * Math.pow(this.prms.If.v, cQ[2]) * Math.sqrt(MacroRugo.g * this.prms.PBD.v); } ... ... @@ -95,7 +109,7 @@ export class MacroRugo extends Nub { } public Equation(sVarCalc: string): Result { const Q = uniroot(this.calcQ, this, 0, 1E7) * this.prms.B.v; const Q = uniroot(this.resolveQ, this, 0, 1E7) * this.prms.B.v; return new Result(Q); } ... ... @@ -121,7 +135,7 @@ export class MacroRugo extends Nub { * Knowledge & Management of Aquatic Ecosystems 45. * @param sVarCalc Variable à calculer */ private calcQ(this: MacroRugo, Q: number): number { private resolveQ(this: MacroRugo, Q: number): number { // Reset cached variables depending on Q (or not...) this._cache = {}; /** Longueur (m) */ ... ... @@ -146,19 +160,18 @@ export class MacroRugo extends Nub { const g = MacroRugo.g; const kappa = 0.41; // von Karman constant // u0 = Averaged velocity at the bed (m.s-1) this.u0 = Q / this.prms.L.v / this.prms.Y.v; // U0 = Averaged velocity (m.s-1) this.U0 = Q / this.prms.B.v / this.prms.Y.v; /** Calulated average velocity */ let u: number; let uMoy: number; if (h / k > MacroRugo.limitSubmerg) { // Submerged conditions /** Velocity at the bed §2.3.2 Cassan et al., 2016 */ /** @todo Faut-il un point fixe ici vu que Cd dépend de u0 ? */ this.u0 = Math.sqrt(2 * g * S * D * this.R / (this.calcCd(this.calc_fFr(this.u0)) * C)); this.u0 = Math.sqrt(2 * g * S * D * this.R / (this.calcCd(this.calc_fFr(this.U0)) * C)); /** turbulent length scale (m) within the blocks layer (alpha_t) */ const alpha = uniroot(this.calcAlpha_t, this, 0, 100); const alpha = uniroot(this.resolveAlpha_t, this, 0, 100); /** averaged velocity at the top of blocks (m.s-1) */ const uk = this.calcUz(alpha); /** Equation (13) Cassan et al., 2016 */ ... ... @@ -175,7 +188,7 @@ export class MacroRugo extends Nub { // calcul intégrale dans la canopée---- // tslint:disable-next-line:variable-name let Qinf: number = this.u0; u = this.u0; let u = this.u0; let uOld: number; const step = 0.1; for (let z = step; z <= 1; z += step) { ... ... @@ -186,15 +199,15 @@ export class MacroRugo extends Nub { Qinf = Qinf * step * k; // Calcul de u moyen u = (Qinf + Qsup) / k; uMoy = (Qinf + Qsup) / k; } else { // Emergent conditions // Resolve equation (4) Cassan et al., 2016 u = uniroot(this.resolveU0, this, 0, 1E7); uMoy = uniroot(this.resolveU0, this, 0, 1E7); } return this.u0 - u; return this.U0 - uMoy; } private get CdChD(): number { ... ... @@ -281,7 +294,7 @@ export class MacroRugo extends Nub { return Math.sqrt(MacroRugo.g * this.prms.If.v * (this.prms.Y.v - this.prms.PBH.v)); } private calcAlpha_t(alpha: number): number { private resolveAlpha_t(alpha: number): number { /** s: minimum distance between blocks */ const s = this.prms.PBD.v * ( 1 / Math.sqrt(this.prms.C.v) - 1); /** Equation(11) Cassan et al., 2016 */ ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!