Commit ad6808bb authored by Dorchies David's avatar Dorchies David
Browse files

feat: Add parameter rQ

- parametrisation of Cx, maxCd, rV and rQ
- parametrisation of the use of min(6, Cd) or Cd=Cx * min(3, fh*)

Refs #295
parent 456916dd
function [f] = calcfFr(Frg)
if abs(Frg)>1.3
f=(Frg.^(-2/3)).^2;
function [f] = calcfFr(Frg, rQ)
if ~exists("rQ","local") then
if(Cd0 > 1.99) then
rQ = 1.25
else
rQ = 1
end
end
if abs(Frg) < 1
f=(min(rQ/(1-(Frg.^2)/4),abs(Frg).^(-2/3))).^2;
else
f=(min(1./(1-(Frg.^2)/4),abs(Frg).^(-2/3))).^2;
f = 1
end
endfunction
......@@ -19,7 +19,7 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
fFr = calcfFr(Frg)
//Cd=Cd0.*(1+0.4./(pf./D).^2);
Cd=Cd0.*min((1+1./(pf./D).^2), 3);
Cd=Cd0.*min(3,(1+0.8./(pf./D).^1.5));
//Cd=Cd0.*(0.8-2*C).*(1+0.4./(pf./D).^2).*fFr;
R=(1-sigma*C);//%.*(1-C.^0.5).^2;
if(bDbg) then
......@@ -88,8 +88,6 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
else
alpha=1-(1.*C).^0.5-1/2*sigma.*C;
if bDbg then print_r('alpha'); end;
Cd2=Cd.*fFr;
if bDbg then print_r('Cd2'); end;
......@@ -104,15 +102,17 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
end
if bDbg then print_r('Cf'); end;
N= (1 .*Cf)./(pf./D.*Cd2.*C);
alpha=1-(1.*C).^0.5-1/2*sigma.*C;
if bDbg then print_r('alpha'); end;
N= (alpha .*Cf)./(pf./D.*Cd2.*C);
if bDbg then print_r('N'); end;
U0a=(2*g.*S.*D.*(R)./(Cd2.*C.*(1+N))).^0.5;
if bDbg then print_r("U0a");end;
U0complet = find_U0_complet(U0a, pf,C,D,sigma,Cd0,Cf,S, bDbg)
U0complet = find_U0_complet(U0a, pf,C,D,sigma,Cd,Cf,S, bDbg)
if bDbg then print_r("U0complet");end;
[u fval] = fminsearch(list(find_U0_complet, pf,C,D,sigma,Cd0,Cf,S), U0a, opt)
[u fval] = fminsearch(list(find_U0_complet, pf,C,D,sigma,Cd,Cf,S), U0a, opt)
res=abs(U0-u);
......
function [res]=find_U0_complet(U01,h,C,D,sigma,Cd0,cfmean,S, bDbgU0)
function [res]=find_U0_complet(U01,h,C,D,sigma,Cd,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(1./(1-(Frexp.^2)/4),Frexp.^(-2/3))).^2;
if Frexp>=1.5
fFr=(Frexp.^(-2/3)).^2;
end
fFr = calcfFr(Frexp)
alpha=1-(1.*C).^0.5-1/2*sigma.*C;
//Cd=Cd0.*(0.8-2*C).*(1+0.4./(h./D).^2).*fFr;
Cd=Cd0.*(1+1./(h./D).^2).*fFr;
Cd2=Cd.*fFr;
N= (alpha.*cfmean)./(h./D.*Cd.*C);
N= (alpha.*cfmean)./(h./D.*Cd2.*C);
if bDbgU0 then
print_r("Frexp")
......@@ -25,7 +19,7 @@ if bDbgU0 then
print_r("N")
end
res=abs(U01-(2*g.*S.*D.*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5);
res=abs(U01-(2*g.*S.*D.*(1-(sigma*C))./(Cd2.*C.*(1+N))).^0.5);
if bDbgU0 then
print_r("res")
end
......
......@@ -35,12 +35,12 @@ function macrorugo_resultComp(z_amont, S, long, Q, L, pf, C, Cd0, h, D)
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
if Cd0 > 2
if Cd0 > 1.99
coeff_contraction=1.5;
else
coeff_contraction=1.1;
coeff_contraction=1.2;
end
Vmax=Vg.*coeff_contraction*calcfFr(Fr);
Vmax=max(Vg, Vg.*coeff_contraction*sqrt(calcfFr(Fr, 1)));
print_r("Vmax");
print_r("V_technique");
end
......
......@@ -16,34 +16,33 @@ z_amont = 12.5; // Cote amont (m)
long = 6; // Longueur rampe (m)
// *****************************************************************************
printf("\n*** Emergent conditions Cd=1.1003***\n")
printf("\n*** Emergent conditions Cd=1.15***\n")
// *****************************************************************************
h = 0.6
k = 0.7
Cd0 = 1.1003 // Forme ronde
h = 0.4
k = 0.6
Cd0 = 1.15 // Forme ronde
macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg)
// *****************************************************************************
printf("\n*** Emergent conditions Cd=2.592***\n")
printf("\n*** Emergent conditions Cd=2***\n")
// *****************************************************************************
h = 0.6
k = 0.7
Cd0 = 2.592 // Forme plane
h = 0.4
k = 0.6
Cd0 = 2.07// Forme plane
macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg)
// *****************************************************************************
printf("\n*** Submerged conditions Cd=1.1003***\n")
printf("\n*** Submerged conditions Cd=1.15***\n")
// *****************************************************************************
k = 0.7
k = 0.6
h = 0.8
Cd0 = 1.1003 // Forme ronde
C = 0.13
Cd0 = 1.15 // Forme ronde
macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg)
// *****************************************************************************
printf("\n*** JalHyd #85 ***\n")
// *****************************************************************************
Cd0 = 1.1003
Cd0 = 1.15
k = 0.8
C = 0.2; // Concentration
Q = []
......
......@@ -4,29 +4,29 @@ import { MacrorugoParams } from "../../src/macrorugo/macrorugo_params";
import { checkResult } from "../test_func";
/*
** Emergent conditions Cd=1.1003***
*** Emergent conditions Cd=1.15***
*** INPUT ***
ks: 0.010000,
D: 0.500000,
k = 0.700000
Cd0: 1.100300,
k = 0.600000
Cd0: 1.150000,
S: 0.050000,
B: 1.000000,
h: 0.600000,
h: 0.400000,
C: 0.130000,
RESULTS:
find_Q_nat(1.052827)=0.600769
Q=0.692366 fVal=0.000000
find_Q_nat(0.686550)=0.719273
Q=0.398841 fVal=0.000000
ZF2: 12.200000,
Vdeb: 1.153943,
Vg: 1.804601,
Fr: 0.743826,
PV: 566.009124,
Vdeb: 0.997103,
Vg: 1.559326,
Fr: 0.787177,
PV: 489.078970,
flowcond: emergent,
Vmax: 2.673506,
V_technique: 1.991299,
q_technique: 0.561860,
Strickler: 7.254351,
Vmax: 2.194842,
V_technique: 1.748990,
q_technique: 0.312101,
Strickler: 8.213879,
*/
function macroRugoInstanceEmergentCd1(): MacroRugo {
......@@ -36,12 +36,12 @@ function macroRugoInstanceEmergentCd1(): MacroRugo {
6, // L
1, // B
0.05, // If
0.692366, // Q
0.6, // h
0.398841, // Q
0.4, // h
0.01, // Ks
0.13, // C
0.5, // D
0.7, // k
0.6, // k
1 // Cd0
)
);
......@@ -50,95 +50,95 @@ function macroRugoInstanceEmergentCd1(): MacroRugo {
const macroRugoExtraResultEmergentCd1: { [key: string]: number | MacroRugoFlowType } = {
ENUM_MacroRugoFlowType: MacroRugoFlowType.EMERGENT,
ZF2: 12.200000,
Vdeb: 1.153943,
Vg: 1.804601,
Fr: 0.743826,
PV: 566.009124,
Vmax: 2.673506,
Strickler: 7.254351
Vdeb: 0.997103,
Vg: 1.559326,
Fr: 0.787177,
PV: 489.078970,
Vmax: 2.194842,
Strickler: 8.213879
};
/*
** Emergent conditions Cd=2.592***
*** Emergent conditions Cd=2***
*** INPUT ***
ks: 0.010000,
D: 0.500000,
k = 0.700000
Cd0: 2.592000,
k = 0.600000
Cd0: 2.070000,
S: 0.050000,
B: 1.000000,
h: 0.600000,
h: 0.400000,
C: 0.130000,
RESULTS:
find_Q_nat(0.675214)=0.321449
Q=0.482344 fVal=0.000000
find_Q_nat(0.503712)=0.614201
Q=0.258032 fVal=0.000000
ZF2: 12.200000,
Vdeb: 0.803907,
Vg: 1.257195,
Fr: 0.518194,
PV: 394.316317,
Vdeb: 0.645079,
Vg: 1.008811,
Fr: 0.509267,
PV: 316.411237,
flowcond: emergent,
Vmax: 2.166970,
V_technique: 1.991299,
q_technique: 0.561860,
Strickler: 5.053822
*/
Vmax: 1.618133,
V_technique: 1.748990,
q_technique: 0.312101,
Strickler: 5.313996
*/
function macroRugoInstanceEmergentCd2(): MacroRugo {
const nub: MacroRugo = macroRugoInstanceEmergentCd1();
nub.prms.Cd0.singleValue = 2;
nub.prms.Q.singleValue = 0.481698;
nub.prms.Q.singleValue = 0.258032;
return nub;
}
const macroRugoExtraResultEmergentCd2: { [key: string]: number | MacroRugoFlowType } = {
ENUM_MacroRugoFlowType: MacroRugoFlowType.EMERGENT,
ZF2: 12.200000,
Vdeb: 0.803907,
Vg: 1.257195,
Fr: 0.518194,
PV: 394.316317,
Vmax: 2.166970,
Strickler: 5.053822
Vdeb: 0.645079,
Vg: 1.008811,
Fr: 0.509267,
PV: 316.411237,
Vmax: 1.618133,
Strickler: 5.313996
};
/*
** Submerged conditions Cd=1.1003***
*** Submerged conditions Cd=1.15***
*** INPUT ***
ks: 0.010000,
D: 0.500000,
k = 0.700000
Cd0: 1.100300,
k = 0.600000
Cd0: 1.150000,
S: 0.050000,
B: 1.000000,
h: 0.800000,
C: 0.130000,
RESULTS:
find_Q_nat(1.403770)=0.143414
Q=1.289038 fVal=0.000000
find_Q_nat(1.373101)=0.102117
Q=1.454795 fVal=0.000000
ZF2: 12.200000,
Vdeb: 1.611298,
Vg: 2.519839,
Fr: 0.899484,
PV: 790.341606,
Vdeb: 1.818493,
Vg: 2.843863,
Fr: 1.015147,
PV: 891.970987,
flowcond: immerge,
q_technique: 0.940450,
Strickler: 8.361756
q_technique: 1.060933,
Strickler: 9.436988,
*/
function macroRugoInstanceSubmerged(): MacroRugo {
const nub: MacroRugo = macroRugoInstanceEmergentCd1();
nub.prms.Y.singleValue = 0.8;
nub.prms.Q.singleValue = 1.289038;
nub.prms.Q.singleValue = 1.454795;
return nub;
}
const macroRugoExtraResultSubmerged: { [key: string]: number | MacroRugoFlowType } = {
ENUM_MacroRugoFlowType: MacroRugoFlowType.SUBMERGED,
ZF2: 12.200000,
Vdeb: 1.611298,
PV: 790.341606,
Strickler: 8.361756
Vdeb: 1.818493,
PV: 891.970987,
Strickler: 9.436988
};
function MacroRugoFactory(sInstance: string): MacroRugo {
......@@ -186,7 +186,7 @@ function testMacroRugoConfig(sInstance: string, Q0: number, fVal0: number, mrExt
// tslint:disable-next-line:no-string-literal
nubit["setFlowType"]();
// tslint:disable-next-line:no-string-literal
expect(nubit["resolveQ"]()).toBeCloseTo(fVal0, 5);
expect(Math.abs(nubit["resolveQ"]())).toBeCloseTo(fVal0, 5);
});
const nub = MacroRugoFactory(sInstance);
for (const prm of nub.prms) {
......@@ -215,8 +215,8 @@ function macroRugoInstanceJalHyd85(): MacroRugo {
describe("Class MacroRugo: ", () => {
testMacroRugoConfig("EmergentCd1", 1.052827, 0.600769, macroRugoExtraResultEmergentCd1);
testMacroRugoConfig("EmergentCd2", 0.675214, 0.321449, macroRugoExtraResultEmergentCd2);
testMacroRugoConfig("EmergentCd1", 0.686550, 0.719273, macroRugoExtraResultEmergentCd1);
testMacroRugoConfig("EmergentCd2", 0.503712, 0.614201, macroRugoExtraResultEmergentCd2);
// Le test passe en debug mais pas sous Jasmine !??
// describe(`Condition Submerged` , () => {
// it(`resolveAlpha_t(0.07) should be 0.074752`, () => {
......@@ -224,7 +224,7 @@ describe("Class MacroRugo: ", () => {
// expect(macroRugoInstanceSubmerged()["resolveAlpha_t"](0.07)).toBeCloseTo(0.074752, 5);
// });
// });
testMacroRugoConfig("Submerged", 1.403770, 0.143414, macroRugoExtraResultSubmerged);
testMacroRugoConfig("Submerged", 1.373101, 0.102117, macroRugoExtraResultSubmerged);
describe("JalHyd #85", () => {
it("CalcSerie Q should return the good result :)", () => {
......@@ -232,7 +232,7 @@ describe("Class MacroRugo: ", () => {
const nubMR2 = macroRugoInstanceJalHyd85();
nubMR.prms.Y.setValues(0.7, 1.2, 0.1);
const aQ: number[] = [0.688476, 0.824387, 1.166239, 1.475853, 1.850393, 2.288123];
const aQ: number[] = [0.681442, 0.808861, 1.132263, 1.429118, 1.790897, 2.215953];
nubMR.CalcSerie();
for (let i = 0; i < aQ.length; i++) {
......
......@@ -874,7 +874,8 @@ describe("vérificateur de franchissement −", () => {
mr.prms.B.singleValue = 1.789; // 2 patterns
mr.prms.C.singleValue = 0.2; // limite la vitesse max
mr.prms.PBH.singleValue = 0.6; // limite la vitesse max
mr.prms.Cd0.singleValue = 2.6; // limite la vitesse max
mr.prms.Cd0.singleValue = 2; // limite la vitesse max
mr.prms.If.singleValue = 0.04;
// vérificateur
const v = new Verificateur();
v.nubToVerify = mr;
......
......@@ -39,10 +39,23 @@ export class MacroRugo extends FishPass {
private _cache: { [key: string]: number };
public paramFhStar: [number, number, number] = [1, 0.8, 1.5];
/** Coefficients used in f_h*(h*) */
private paramFhStar: [number, number, number] = [1, 0.8, 1.5];
/** Rapport pour le calcul de fFr */
public rQ: number = 1.;
/** Coefficient used in Cx */
private paramCx: [number, number] = [0.92, 0.23];
/** Coefficient used in rQ */
private paramRQ: [number, number] = [0.25, 0.75];
/** Coefficient used in rQ */
private paramRV: [number, number] = [0.3, 0.9];
/** Maximum value for Cd */
private paramMaxCd: number = 6;
/** true: Cd0 * min(3, fh), false : min(6, Cd0 * fh) */
private paramCdNewVersion: boolean = true;
constructor(prms: MacrorugoParams, dbg: boolean = false) {
super(prms, dbg);
......@@ -306,11 +319,11 @@ export class MacroRugo extends FishPass {
/**
* Interpolation of Cd0 for Cd from calibrated Cd0 (See #291)
* Cx = 1.1 for Cd0 = 1 and Cx = 2.592 for Cd0 = 2
* Cx = 1.15 for Cd0 = 1 and Cx = 2 for Cd0 = 2
*/
private get Cx(): number {
if (this._cache.Cx === undefined) {
this._cache.Cx = this.prms.Cd0.v //* 1.4917 -0.3914;
this._cache.Cx = this.paramCx[0] * this.prms.Cd0.v + this.paramCx[1];
}
return this._cache.Cx;
}
......@@ -320,7 +333,11 @@ export class MacroRugo extends FishPass {
*/
private get Cd(): number {
if (this._cache.Cd === undefined) {
this._cache.Cd = Math.min(6, this.Cx * (this.paramFhStar[0] + this.paramFhStar[1] / Math.pow(this.prms.Y.v / this.prms.PBD.v, this.paramFhStar[2])));
if(this.paramCdNewVersion) {
this._cache.Cd = this.Cx * Math.min(this.paramMaxCd, (this.paramFhStar[0] + this.paramFhStar[1] / Math.pow(this.prms.Y.v / this.prms.PBD.v, this.paramFhStar[2])));
} else {
this._cache.Cd = Math.min(this.paramMaxCd, this.Cx * (this.paramFhStar[0] + this.paramFhStar[1] / Math.pow(this.prms.Y.v / this.prms.PBD.v, this.paramFhStar[2])));
}
}
return this._cache.Cd;
}
......@@ -421,31 +438,22 @@ export class MacroRugo extends FishPass {
}
private resolveQEmergent(): number {
// tslint:disable-next-line: variable-name
const Cd = this.Cd * this.fFr;
/** N from Cassan 2016 eq(2) et Cassan 2014 eq(12) */
const N = (1 * this.calcCf(this.prms.Y.v)) / (this.prms.Y.v / this.prms.PBD.v * Cd * this.prms.C.v);
// const U0i = Math.sqrt(
// 2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v *
// (1 - this.sigma * this.prms.C.v) / (Cd * this.prms.C.v * (1 + N))
// );
return this.U0 - uniroot(this.resolveU0Complete, this, 1E-6, 100);
}
private resolveU0Complete(U0: number): number {
const fFr = this.CalcfFr(U0);
const alpha = 1 - Math.pow(1 * this.prms.C.v, 0.5) - 0.5 * this.sigma * this.prms.C.v;
// tslint:disable-next-line: variable-name
const Cd = this.Cd * this.CalcfFr(U0);
const N = (alpha * this.calcCf(this.prms.Y.v)) /
(this.prms.Y.v / this.prms.PBD.v * this.Cd * fFr * this.prms.C.v);
(this.prms.Y.v / this.prms.PBD.v * Cd * this.prms.C.v);
return U0 - Math.sqrt(
2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v *
(1 - this.sigma * this.prms.C.v) / (this.Cd * fFr * this.prms.C.v * (1 + N))
(1 - this.sigma * this.prms.C.v) / (Cd * this.prms.C.v * (1 + N))
);
}
......@@ -455,11 +463,22 @@ export class MacroRugo extends FishPass {
*/
private get rV(): number {
if (this._cache.rV === undefined) {
this._cache.rV = 0.4 * this.prms.Cd0.v + 0.7;
this._cache.rV = this.paramRV[0] * this.prms.Cd0.v + this.paramRV[1];
}
return this._cache.rV;
}
/**
* Perte de charge supplémentaire due à la forme (voir fFr)
* r = 1 pour un plot circulaire Cd0​=1 et r = 1.25 pour un plot à face plane Cd0​=2
*/
private get rQ(): number {
if (this._cache.rQ === undefined) {
this._cache.rQ = this.paramRQ[0] * this.prms.Cd0.v + this.paramRQ[1];
}
return this._cache.rQ;
}
/**
* Froude correction function (Cassan et al. 2014, Eq. 19)
*/
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment