An error occurred while loading the file. Please try again.
-
Dorchies David authoredb90e0bfe
import { Nub } from "../nub";
import { ParamCalculability } from "../param/param-definition";
import { ParamValueMode } from "../param/param-value-mode";
import { Result } from "../util/result";
import { MacrorugoParams } from "./macrorugo_params";
export { MacrorugoParams };
export enum MacroRugoFlowType {
EMERGENT,
QUASI_EMERGENT,
IMMERGE
}
export class MacroRugo extends Nub {
private static readonly g = 9.81;
/** nu: water kinematic viscosity */
private static readonly nu = 1E-6;
// Water at 20 °C has a kinematic viscosity of about 10−6 m2·s−1
// (https://en.wikipedia.org/wiki/Viscosity#Kinematic_viscosity,_%CE%BD)
/** Ratio between the width (perpendicular to flow) and the lenght (parallel to flow) of a cell (-) */
private static readonly fracAxAy = 1;
/** Limit between emergent and submerged flow */
private static readonly limitSubmerg = 1.1;
/** Rugosité de fond (m) */
private ks: number;
/** 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 = {};
}
/**
* paramètres castés au bon type
*/
get prms(): MacrorugoParams {
return this._prms as MacrorugoParams;
}
/**
* Calcul du débit total, de la cote amont ou aval ou d'un paramètre d'une structure
* @param sVarCalc Nom du paramètre à calculer :
* "Q", "Z1", "Z2" ou "n.X" avec "n" l'index de l'ouvrage et "X" son paramètre
* @param rInit Valeur initiale
* @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.fracAxAy * this.prms.C.v))
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
/ 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
r.extraResults.PV = 1000 * MacroRugo.g * this.V(this.prms.Q) / this.prms.B.v * this.prms.If.v;
// Type d'écoulement
if (this.prms.Y.v / this.prms.PBH.v < 1) {
r.extraResults.ENUM_MacroRugoFlowType = MacroRugoFlowType.EMERGENT;
} else if (this.prms.Y.v / this.prms.PBH.v < MacroRugo.limitSubmerg) {
r.extraResults.ENUM_MacroRugoFlowType = MacroRugoFlowType.QUASI_EMERGENT;
} else {
r.extraResults.ENUM_MacroRugoFlowType = MacroRugoFlowType.IMMERGE;
}
// Vitesse et débit du guide technique
let cQ: [number, number, number, number];
let cV: [number, number, number];
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];
} else {
cQ = [0.815, 1.45, 0.557, -0.456];
cV = [4.54, 0.32, 0.56];
}
}
r.extraResults.Q_GuideTech = 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.PBD.v * this.prms.B.v;
if (this.prms.Y.v / this.prms.PBH.v <= MacroRugo.limitSubmerg) {
r.extraResults.V_GuideTech = 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);
}
return r;
}
public Equation(sVarCalc: string): Result {
const Q = uniroot(this.resolveQ, this, 0, 1E7) * this.prms.B.v;
return new Result(Q);
}
/**
* paramétrage de la calculabilité des paramètres
*/
protected setParametersCalculability() {
this.prms.ZF1.calculability = ParamCalculability.FREE;
this.prms.L.calculability = ParamCalculability.FREE;
this.prms.Ks.calculability = ParamCalculability.FREE;
this.prms.B.calculability = ParamCalculability.DICHO;
this.prms.If.calculability = ParamCalculability.DICHO;
this.prms.Q.calculability = ParamCalculability.EQUATION;
this.prms.Y.calculability = ParamCalculability.DICHO;
this.prms.C.calculability = ParamCalculability.DICHO;
this.prms.PBD.calculability = ParamCalculability.FREE;
this.prms.PBH.calculability = ParamCalculability.FREE;
this.prms.Cd0.calculability = ParamCalculability.FREE;
}
/**
* Equation from Cassan, L., Laurens, P., 2016. Design of emergent and submerged rock-ramp fish passes.
* Knowledge & Management of Aquatic Ecosystems 45.
* @param sVarCalc Variable à calculer
*/
private resolveQ(this: MacroRugo, Q: number): number {
// Reset cached variables depending on Q (or not...)
this._cache = {};
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
/** Longueur (m) */
const L: number = this.prms.L.v;
/** Tirant d'eau (m) */
const h: number = this.prms.Y.v;
/** Paramètre de bloc : Forme (1 pour rond, 2 pour carré)
* drag coefficient of a block considering a single block
* infinitely high with F ≪ 1;
*/
// tslint:disable-next-line:variable-name
const Cd0: number = this.prms.Cd0.v;
/** Concentration de blocs (-) */
const C: number = this.prms.C.v;
/** Paramètre de bloc : Diamètre (m) */
const D: number = this.prms.PBD.v;
/** Paramètre de bloc : Hauteur (m) */
const k: number = this.prms.PBH.v;
/** Pente (m/m) */
const S: number = this.prms.If.v;
const g = MacroRugo.g;
const kappa = 0.41; // von Karman constant
// U0 = Averaged velocity (m.s-1)
this.U0 = Q / this.prms.B.v / h;
/** Calulated average velocity */
let uMoy: number;
if (h / k > MacroRugo.limitSubmerg) {
// Submerged conditions
/** Velocity at the bed §2.3.2 Cassan et al., 2016 */
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.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 */
const d = k - 1 / kappa * alpha * uk / this.ustar;
/** Equation (14) Cassan et al., 2016 */
const z0 = (k - d) * Math.exp(- kappa * uk / this.ustar);
/** Integral of Equation (12) Cassan et al., 2016 */
// tslint:disable-next-line:variable-name
const Qsup = this.ustar / kappa * (
(h - d) * (Math.log((h - d) / z0) - 1)
- ((k - d) * (Math.log((k - d) / z0) - 1))
);
// calcul intégrale dans la canopée----
// tslint:disable-next-line:variable-name
let Qinf: number = this.u0;
let u = this.u0;
let uOld: number;
const step = 0.1;
for (let z = step; z <= 1; z += step) {
uOld = u;
u = this.calcUz(alpha, z);
Qinf += (uOld + u) / 2;
}
Qinf = Qinf * step * k;
// Calcul de u moyen
uMoy = (Qinf + Qsup) / k;
} else {
// Emergent conditions
// Resolve equation (4) Cassan et al., 2016
uMoy = uniroot(this.resolveU0, this, 0, 1E7);
}
return this.U0 - uMoy;
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
}
private get CdChD(): number {
if (this._cache.CdChD !== undefined) {
return this._cache.CdChD;
}
return this.calcCd(1) * this.prms.C.v * this.prms.PBH.v / this.prms.PBD.v;
}
/**
* sigma ratio between the block area in the x, y plane and D2
*/
private get sigma(): number {
if (this._cache.sigma !== undefined) {
return this._cache.sigma;
}
if (this.prms.Cd0.v === 2) {
return 1;
} else {
return Math.PI / 4;
}
}
private get R(): number {
if (this._cache.R !== undefined) {
return this._cache.R;
}
return (1 - this.sigma * this.prms.C.v);
}
/**
* Bed friction coefficient Equation (3) (Cassan et al., 2016)
*/
private calcCf(U0: number): number {
if (this.prms.Ks.v < 1E-6) {
// Between Eq (8) and (9) (Cassan et al., 2016)
const reynolds = U0 * this.prms.Y.v / MacroRugo.nu;
return 0.3164 / 4. * Math.pow(reynolds, -0.25);
} else {
// Equation (3) (Cassan et al., 2016)
return 2 / Math.pow(5.1 * Math.log10(this.prms.Y.v / this.prms.Ks.v - 1) + 6, 2);
}
}
/**
* Calculation of Cd : drag coefficient of a block under the actual flow conditions
* @param fFr
*/
private calcCd(fFr: number): number {
return this.prms.Cd0.v * (1 + 0.4 / Math.pow(this.prms.Y.v / this.prms.PBD.v, 2)) * fFr;
}
/**
* Calcul de Beta force ratio between drag and turbulent stress (Cassan et al. 2016 eq(8))
* \Beta = (k / alpha_t) (C_d C k / D) / (1 - \sigma C)
* @param alpha \alpha_t turbulent length scale (m) within the blocks layer
*/
private calcBeta(alpha: number): number {
return Math.sqrt(this.prms.PBH.v * this.CdChD / alpha / this.R);
}
/**
* Averaged velocity at a given vertical position (m.s-1)
* @param alpha turbulent length scale (m) within the blocks layer
* @param z dimensionless vertical position z / k
*/
private calcUz(alpha: number, z: number = 1): number {
const beta = this.calcBeta(alpha);
// Equation (9) Cassan et al., 2016
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
return this.u0 * Math.sqrt(
beta * (this.prms.Y.v / this.prms.PBH.v - 1) * Math.sinh(beta * z) / Math.cosh(beta) + 1
);
}
private get ustar(): number {
if (this._cache.ustar !== undefined) {
return this._cache.ustar;
}
return Math.sqrt(MacroRugo.g * this.prms.If.v * (this.prms.Y.v - this.prms.PBH.v));
}
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 */
const l0 = Math.min(s, 0.15 * this.prms.PBH.v);
// Equation(10) Cassan et al., 2016
return alpha * this.calcUz(alpha) - l0 * this.ustar;
}
private resolveU0(U0: number): number {
const g = MacroRugo.g;
const alpha = 1 - (1 / MacroRugo.fracAxAy * this.prms.C.v);
// tslint:disable-next-line:variable-name
const Cd = this.calcCd(this.calc_fFr(U0));
/** N from Cassan 2016 eq(2) et Cassan 2014 eq(12) */
const N = (alpha * this.calcCf(U0)) / (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) / (Cd * this.prms.C.v * (1 + N))
);
}
/**
* Froude correction function
* @param u0 Average velocity
*/
private calc_fFr(u0: number): number {
// tslint:disable-next-line:variable-name
const Fr = u0 / (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);
} else {
return Math.pow(Fr, -4 / 3);
}
}
}
/**
* Searches the interval from <tt>lowerLimit</tt> to <tt>upperLimit</tt>
* for a root (i.e., zero) of the function <tt>func</tt> with respect to
* its first argument using Brent's method root-finding algorithm.
*
* Translated from zeroin.c in http://www.netlib.org/c/brent.shar.
*
* Copyright (c) 2012 Borgar Thorsteinsson <borgar@borgar.net>
* MIT License, http://www.opensource.org/licenses/mit-license.php
*
* @param {function} func function for which the root is sought.
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
* @param {number} lowerlimit the lower point of the interval to be searched.
* @param {number} upperlimit the upper point of the interval to be searched.
* @param {number} errorTol the desired accuracy (convergence tolerance).
* @param {number} maxIter the maximum number of iterations.
* @returns an estimate for the root within accuracy.
*
*/
function uniroot<T>(func: (param: number) => number, thisArg: T, lowerLimit: number, upperLimit: number,
errorTol: number = 0, maxIter: number = 1000
) {
let a = lowerLimit;
let b = upperLimit;
let c = a;
let fa = func.call(thisArg, a);
let fb = func.call(thisArg, b);
let fc = fa;
let tolAct; // Actual tolerance
let newStep; // Step at this iteration
let prevStep; // Distance from the last but one to the last approximation
let p; // Interpolation step is calculated in the form p/q; division is delayed until the last moment
let q;
while (maxIter-- > 0) {
prevStep = b - a;
if (Math.abs(fc) < Math.abs(fb)) {
// Swap data for b to be the best approximation
a = b, b = c, c = a;
fa = fb, fb = fc, fc = fa;
}
tolAct = 1e-15 * Math.abs(b) + errorTol / 2;
newStep = (c - b) / 2;
if (Math.abs(newStep) <= tolAct || fb === 0) {
return b; // Acceptable approx. is found
}
// Decide if the interpolation can be tried
if (Math.abs(prevStep) >= tolAct && Math.abs(fa) > Math.abs(fb)) {
// If prev_step was large enough and was in true direction, Interpolatiom may be tried
let t1;
let cb;
let t2;
cb = c - b;
if (a === c) { // If we have only two distinct points linear interpolation can only be applied
t1 = fb / fa;
p = cb * t1;
q = 1.0 - t1;
} else { // Quadric inverse interpolation
q = fa / fc, t1 = fb / fc, t2 = fb / fa;
p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1));
q = (q - 1) * (t1 - 1) * (t2 - 1);
}
if (p > 0) {
q = -q; // p was calculated with the opposite sign; make p positive
} else {
p = -p; // and assign possible minus to q
}
if (p < (0.75 * cb * q - Math.abs(tolAct * q) / 2) &&
p < Math.abs(prevStep * q / 2)) {
// If (b + p / q) falls in [b,c] and isn't too large it is accepted
newStep = p / q;
}
// If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
}
421422423424425426427428429430431432433434435436
if (Math.abs(newStep) < tolAct) { // Adjust the step to be not less than tolerance
newStep = (newStep > 0) ? tolAct : -tolAct;
}
a = b, fa = fb; // Save the previous approx.
b += newStep, fb = func.call(thisArg, b); // Do step to a new approxim.
if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
c = a, fc = fa; // Adjust c for it to have a sign opposite to that of b
}
}
return undefined;
}