#22 Ajout du nub MacroRugo

src/macrorugo/macrorugo.ts

+import { Nub } from "../nub";
+import { ParamCalculability } from "../param/param-definition";
+import { Result } from "../util/result";
+import { MacrorugoParams } from "./macrorugo_params";
+
+export class MacroRugo extends Nub {
+
+    /** Rugosité de fond (m) */
+    private ks: number;
+    /** Averaged velocity at the bed (m.s-1) */
+    private u0: number;
+
+    private _cache: { [key: string]: number }
+
+    static readonly g = 9.81;
+
+    /** nu: water kinematic viscosity  */
+    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 lenght (parallel to flow) and the width (perpendicular to flow) of a cell (-) */
+    static readonly frac_ay_ax = 1;
+
+    constructor(prms: MacrorugoParams, dbg: boolean = false) {
+        super(prms, dbg);
+    }
+
+    /**
+     * paramètres castés au bon type
+     */
+    get prms(): MacrorugoParams {
+        return this._prms as MacrorugoParams;
+    }
+
+    public Equation(sVarCalc: string): Result {
+
+
+        const Q = uniroot(this.calcQ, 0, 1E7);
+        return new Result(Q);
+    }
+
+    /**
+     * 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
+     */
+    public calcQ(Q: number): number {
+        // Reset cached variables depending on Q (or not...)
+        this._cache = {};
+        /** 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;
+         */
+        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
+
+
+        /** Calulated average velocity */
+        let u: number;
+        if (h / k > 1.1) {
+            // Submerged conditions
+
+            /** Velocity at the bed §2.3.2 Cassan et al., 2016 */
+            this.u0 = Math.sqrt(2 * g * S * D * this.R / (Cd0 * C));
+            /** turbulent length scale (m) within the blocks layer (alpha_t) */
+            const alpha = uniroot(this.calcAlpha_t,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 */
+            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----
+            let Qinf: number = this.u0;
+            u = this.u0;
+            let u_old: number;
+            const step = 0.1;
+            for (let z = step; z <= 1; z += step) {
+                u_old = u;
+                u = this.calcUz(alpha, z);
+                Qinf += (u_old + u) / 2
+            }
+            Qinf = Qinf * step * k;
+
+            // Calcul de u moyen
+            u = (Qinf + Qsup) / k;
+
+         } else {
+            // Emergent conditions
+
+            // u0 = Averaged velocity at the bed (m.s-1)
+            this.u0 = Q / this.prms.L.v / this.prms.Y.v;
+
+            //
+            u = uniroot(this.calcU0, 0, 1E7)
+        }
+        return this.u0 - u;
+    }
+
+    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) * Math.pow(1 - Math.sqrt(this.prms.C.v), 2);
+    }
+
+    /**
+     * Bed friction coefficient Equation (3) (Cassan et al., 2016)
+     */
+    private get Cf(): number {
+        // Between Eq (8) and (9) (Cassan et al., 2016)
+        const Re = this.u0 * this.prms.Y.v / MacroRugo.nu
+
+        if (this.prms.Ks.v < 1E-6) {
+            return 0.3164/4.* Math.pow(Re, -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);
+        }
+    }
+
+    /**
+     * 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;
+    }
+
+    /**
+     * Calculation of Cd : drag coefficient of a block under the actual flow conditions
+     * @param Cd0
+     * @param pf
+     * @param D
+     * @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
+        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 calcAlpha_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 calcU0(U01: number, C: number, h: number, r: number, Cd0: number, D: number, cfmean: number, S: number, sigma: number) {
+        const g = MacroRugo.g;
+
+        const alpha = 1 - (1 * C);
+        const Cd = this.calcCd(this.calc_fFr(U01));
+
+        /** N from Cassan 2016 eq(2) et Cassan 2014 eq(12) */
+        const N = (alpha * this.Cf) / (h / D * Cd * C);
+
+        return U01 - Math.sqrt(2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v * (1 - this.sigma * this.prms.C.v) / (Cd * C * (1 + N)));
+    }
+
+    /**
+     * Froude correction function
+     * @param u0 Average velocity
+     */
+    private calc_fFr(u0: number): number {
+
+        const Fr = u0 / (1 - Math.sqrt(MacroRugo.frac_ay_ax * 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.
+ * @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(func: Function, lowerLimit: number, upperLimit: number,
+                 errorTol: number = 0, maxIter: number = 1000
+) {
+    let a = lowerLimit
+        , b = upperLimit
+        , c = a
+        , fa = func(a)
+        , fb = func(b)
+        , fc = fa
+        , s = 0
+        , fs = 0
+        , tol_act   // Actual tolerance
+        , new_step  // Step at this iteration
+        , prev_step // Distance from the last but one to the last approximation
+        , p         // Interpolation step is calculated in the form p/q; division is delayed until the last moment
+        , q
+        ;
+
+    while (maxIter-- > 0) {
+
+        prev_step = 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;
+        }
+
+        tol_act = 1e-15 * Math.abs(b) + errorTol / 2;
+        new_step = (c - b) / 2;
+
+        if (Math.abs(new_step) <= tol_act || fb === 0) {
+            return b; // Acceptable approx. is found
+        }
+
+        // Decide if the interpolation can be tried
+        if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
+            // If prev_step was large enough and was in true direction, Interpolatiom may be tried
+            var t1, cb, 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(tol_act * q) / 2) &&
+                p < Math.abs(prev_step * q / 2)) {
+                // If (b + p / q) falls in [b,c] and isn't too large it is accepted
+                new_step = p / q;
+            }
+
+            // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
+        }
+
+        if (Math.abs(new_step) < tol_act) { // Adjust the step to be not less than tolerance
+            new_step = (new_step > 0) ? tol_act : -tol_act;
+        }
+
+        a = b, fa = fb;     // Save the previous approx.
+        b += new_step, fb = func(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;
+
+}
\ No newline at end of file

src/macrorugo/macrorugo_params.ts

+import { ParamDefinition } from "../param/param-definition";
+import { ParamDomainValue } from "../param/param-domain";
+import { ParamsEquation } from "../param/params-equation";
+
+export class MacrorugoParams extends ParamsEquation {
+
+    /** Cote de fond amont (m) */
+    private _ZF1: ParamDefinition;
+
+    /** Longueur (m) */
+    private _L: ParamDefinition;
+
+    /** Largeur (m) */
+    private _B: ParamDefinition;
+
+    /** Pente (m/m) */
+    private _If: ParamDefinition;
+
+    /** Débit (m3/s) */
+    private _Q: ParamDefinition;
+
+    /** Tirant d'eau (m) */
+    private _Y: ParamDefinition;
+
+    /** Rugosité de fond (m) */
+    private _Ks: ParamDefinition;
+
+    /** Concentration de blocs (-) */
+    private _C: ParamDefinition;
+
+    /** Paramètre de bloc : Diamètre (m) */
+    private _PBD: ParamDefinition;
+
+    /** Paramètre de bloc : Hauteur (m) */
+    private _PBH: ParamDefinition;
+
+    /** Paramètre de bloc : Forme (1 pour rond, 2 pour carré) */
+    private _Cd0: ParamDefinition;
+
+    /**
+     *
+     * @param rZF1 Cote de fond amont (m)
+     * @param rL Longueur (m)
+     * @param rB Largeur (m)
+     * @param rIf Pente (m/m)
+     * @param rQ Débit (m3/s)
+     * @param rY Tirant d'eau (m)
+     * @param rRF Rugosité de fond (m)
+     * @param rCB Concentration de blocs (m)
+     * @param rPBD Paramètre de bloc : Diamètre (m)
+     * @param rPBH Paramètre de bloc : Hauteur (m)
+     * @param rCd0 Paramètre de bloc : Forme (1 pour rond, 2 pour carré)
+     */
+    constructor(rZF1: number, rL: number, rB: number, rIf: number, rQ: number,
+                rY: number, rRF: number, rCB: number, rPBD: number, rPBH: number, rCd0: number) {
+        super();
+
+        this._ZF1 = new ParamDefinition("ZF1", ParamDomainValue.POS, rZF1);
+        this.addParamDefinition(this._ZF1);
+
+        this._L = new ParamDefinition("L", ParamDomainValue.POS, rL);
+        this.addParamDefinition(this._L);
+
+        this._B = new ParamDefinition("B", ParamDomainValue.POS, rB);
+        this.addParamDefinition(this._B);
+
+        this._If = new ParamDefinition("If", ParamDomainValue.POS, rIf);
+        this.addParamDefinition(this._If);
+
+        this._Q = new ParamDefinition("Q", ParamDomainValue.POS, rQ);
+        this.addParamDefinition(this._Q);
+
+        this._Y = new ParamDefinition("Y", ParamDomainValue.POS, rY);
+        this.addParamDefinition(this._Y);
+
+        this._Ks = new ParamDefinition("Ks", ParamDomainValue.POS, rRF);
+        this.addParamDefinition(this._Ks);
+
+        this._C = new ParamDefinition("C", ParamDomainValue.POS, rCB);
+        this.addParamDefinition(this._C);
+
+        this._PBD = new ParamDefinition("PBD", ParamDomainValue.POS, rPBD);
+        this.addParamDefinition(this._PBD);
+
+        this._PBH = new ParamDefinition("PBH", ParamDomainValue.POS, rPBH);
+        this.addParamDefinition(this._PBH);
+
+        this._Cd0 = new ParamDefinition("Cd0", ParamDomainValue.POS, rCd0);
+        this.addParamDefinition(this._Cd0);
+
+    }
+
+    /**
+     * Cote de fond amont (m)
+     * @return {ParamDefinition}
+     */
+    public get ZF1(): ParamDefinition {
+        return this._ZF1;
+    }
+
+    /**
+     * Longueur (m)
+     * @return {ParamDefinition}
+     */
+    public get L(): ParamDefinition {
+        return this._L;
+    }
+
+    /**
+     * Getter B
+     * @return {ParamDefinition}
+     */
+    public get B(): ParamDefinition {
+        return this._B;
+    }
+
+    /**
+     * Pente (m/m)
+     * @return {ParamDefinition}
+     */
+    public get If(): ParamDefinition {
+        return this._If;
+    }
+
+    /**
+     * Débit (m3/s)
+     * @return {ParamDefinition}
+     */
+    public get Q(): ParamDefinition {
+        return this._Q;
+    }
+
+    /**
+     * Tirant d'eau (m)
+     * @return {ParamDefinition}
+     */
+    public get Y(): ParamDefinition {
+        return this._Y;
+    }
+
+    /**
+     * Rugosité de fond (m)
+     * @return {ParamDefinition}
+     */
+    public get Ks(): ParamDefinition {
+        return this._Ks;
+    }
+
+    /**
+     * Concentration de blocs (-)
+     * @return {ParamDefinition}
+     */
+    public get C(): ParamDefinition {
+        return this._C;
+    }
+
+    /**
+     * Paramètre de bloc : Diamètre (m)
+     * @return {ParamDefinition}
+     */
+    public get PBD(): ParamDefinition {
+        return this._PBD;
+    }
+
+    /**
+     * Paramètre de bloc : Hauteur (m)
+     * @return {ParamDefinition}
+     */
+    public get PBH(): ParamDefinition {
+        return this._PBH;
+    }
+
+    /**
+     * Paramètre de bloc : Forme (1 pour rond, 2 pour carré)
+     * @return {ParamDefinition}
+     */
+    public get Cd0(): ParamDefinition {
+        return this._Cd0;
+    }
+
+}