SIC^2 : Logiciel de Simulation Intégrée des Canaux et de leur Contrôle

Détachement Algual en réponse à une chasse

David Dorchies — 2017-11-07T11:47:02Z

Définition

Le taux de détachement accidentel de biomasse fixée est calculé, en chaque section comme suit :

$ {\frac{\partial B}{\partial t} (x,t) = - \frac{1}{\delta }\left ( \frac{\tau_{0}(x,t) - \tau_{0,cr}(x)}{\tau_{ref}(x)} - s_B \right )^\eta B(x,t) $ , si $ \frac{\tau_{0}(x,t) - \tau_{0,cr}(x)}{\tau_{ref}(x)} > s_B} $

$ {\frac{\partial B}{\partial t} = 0} $ , sinon

Avec :

$ {B(x,t)} $ : la biomasse fixée en kg m^-1
$ {\tau _ 0(x,t)} $ : la contrainte de cisaillement à la paroi en N m^-2
$ {\tau_{0,cr}(x)} $ : la valeur critique de cisaillement qui intègre l'effet mémoire des algues dans leur capacité de résistance.
$ {\tau_{ref}(x)} $ : la valeur de référence du cisaillement
$ {s_B} $ : le seuil de sensibilité qui prend en compte les fluctuations antérieures à la chasse
$ {\delta} $ : une constante de temps en secondes
$ {\eta} $ : un exposant adimensionnel

Caractéristiques

Identifiant de la loi : 201
Nombre de classes intervenant : 1
Nombre de paramètres : 5

Classes intervenantes :

$ {B} $

Paramètres fonction de l'abscisse :

$ {\tau_{0,cr}} $
$ {\tau_{ref}} $

Paramètres éventuellement fonction du temps

$ {s_B} $
$ {\delta} $
$ {\eta} $

Engelund-Hansen (1967)

Louis Poirel — 2012-04-30T13:01:27Z

Loi d'Engelund-Hansen

La Concentration d'équilibre est calculée de la façon suivante :

$$C_{eq}=0.05\rho_S\frac{LU^2}{Q}\frac{(JR)^{\frac{3}{2}}}{\sqrt{g} (\rho_S/\rho-1)^2d}$$

Avec :

$ \rho_S$ : la masse volumique du sédiment en kg/m³
$ L $ : la largeur du cours d'eau en m
$ U $ : la vitesse moyenne en m/s
$ Q $ : le débit en m³/s
$ J $ : la pente en m/m
$ R $ : le rayon hydraulique en m
$ g $ : l'accélération de la gravité en m/s²
$ \rho $ : la masse volumique de l'eau en kg/m³
$ d $ : le diamètre du sédiment

Caractéristiques

Identifiant de la loi : 561
Nombre de classes intervenant : 4
Nombre de paramètres : 6

Classes intervenantes :

$ C_i $ : Classe dérivante variant sous l'effet de la loi
$ C_j $ : Classe fixée variant sous l'effet de la loi
$ C_k $ : Classe paramètre de la loi ($i=k$ pour un calcul standard)
$ T $ : Température (utilisée pour calculer la viscosité de l'eau)

Paramètres éventuellement fonction du temps

$ d $ : le diamètre du sédiment
$ \rho_S$ : la masse volumique du sédiment
$ p $ : la porosité du sédiment
$ \alpha$
$ i_{ech}$ : formule d'échange choisie :1 pour Han, 2 pour Hazen
$\beta$

Croissance des algues (Thèse O. Fovet, 2010, p.101)

David Dorchies — 2012-04-30T08:07:27Z

Définition

La croissance de la biomasse alguale fixée est calculé au temps t, en chaque section d'abscisse x comme suit :

$\frac{\partial B}{\partial t}(x,t) = \mu(x,t) B(x,t) F_{lim}(B'(x,t))$

avec :

$F_{lim}(B(x,t)) = \left ( 1 - \frac{B'(x,t)}{B_{Max}} \right )$
$ \mu(x,t) = \mu_{0} \theta^{T(t)-T_{0}} \frac{I(x,t)}{I_{opt}} e^{1- \frac{I(x,t)}{I_{opt}}} \textup{min} \left ( \frac{N_{i}(x,t)}{N_{i}(x,t) + K_{N_{I}}} \right )$
$ I(x,t) = I_{s}(x,t) e^{-k_{ext} h(x,t)} $
$ I_{s}(x,t) = (1-C_m(x))(1-a)R_N(t) $

Variables et paramètres

$ B(x,t) $ : la biomasse fixée en kg m^-1
$ B'(x,t) $ : la biomasse fixée en kg m^-1
$ B_{Max} $ : La valeur maximale de la biomasse fixée pour un mètre linéaire de canal en kg m^-1
$ \mu_{0} $ : le taux de croissance de référence en s ^-1
$ \theta $ : un coefficient de croissance
$ T(x,t) $ : la température de l'eau en °C
$ T_{0} $ : la température de référence en °C
$ I(x,t) $ : l'intensité lumineuse solaire au fond du canal en W m^-2
$ I_{s}(x,t) $ : l'intensité lumineuse solaire en W m^-2
$ C_{m}(x)$ : le coefficient de masquage (cf. Simulation de la température)
$ a $ : albédo (cf. Simulation de la température)
$ R_{N}(t) $ le rayonnement net en W m^-2 (cf. Simulation de la température)
$ k_{ext} $ : facteur d'extinction lié à la turbidité
$ h(x,t) $ : hauteur d'eau moyenne ($ h = S / L $) en m
$ I_{opt}(x,t) $ : l'intensité lumineuse optimale en W m^-2
$ N_{i}(x,t) $ : concentration du nutriment i en kg m^-3
$ K_{N_{i}} $ : constante de demi-saturation du nutriment limitant i en kg m^-3

Nutriments limitants

Cette loi d'échange permet de saisir au maximum 3 nutriments au choix intervenant dans la croissance de l'algue. POur ne pas faire intervenir de nutriment limitant, il suffit de laisser les valeurs $ K_{N_{i}} $ à zéro.

Caractéristiques

Identifiant de la loi : 301
Nombre de classes intervenant : 6
Nombre de paramètres météo : 3
Nombre de paramètres utilisateur : 9

Classes intervenantes :

$ B(x,t) $ : Classe variant sous l'effet de la loi
$ B'(x,t) $ : Classe paramètre de la loi
$ T(x,t) $ : Température de l'eau
$ N_{1}(x,t) $ : Nutriment 1
$ N_{2}(x,t) $ : Nutriment 2
$ N_{3}(x,t) $ : Nutriment 3

Paramètres météorologiques :

$ C_{m}(x)$
$ a $
$ R_{N}(t) $

Paramètres éventuellement fonction du temps

$ B_{Max} $ : Biomasse maximale
$ \mu_{0} $ : Taux de croissance de référence
$ \theta $ : Coefficient de croissance
$ T_{0} $ : Température de référence
$ k_{ext} $ : Coefficient d'extinction
$ I_{opt} $ : Intensité lumineuse optimale
$ K_{N_{1}} $ : Constante de demi-saturation du nutriment 1
$ K_{N_{2}} $ : Constante de demi-saturation du nutriment 2
$ K_{N_{3}} $ : Constante de demi-saturation du nutriment 3

Calcul de la dérivée de l'échange

La dérivée de l'échange par rapport à $B(x,t)$ est égale à :
$ E_{B}'(x,t) = \mu (x,t) \left ( 1 - \frac{2B(x,t)}{B_{Max}} \right ) $

Principe des lois d'évolution

Louis Poirel — 2012-04-27T14:17:32Z

Principe des lois d'évolution

Une loi d'évolution fait varier une ou plusieurs classes de qualité, en fonction de paramètres spécifiques à la loi, de variables hydrauliques, et éventuellement des concentrations de certaines classes qualité.

Ainsi, les termes d'échanges relatifs à la loi $i$ sont calculés de la façon suivante :

$$\left(E_{k_1}^i, \dots, E_{k_n}^i\right)=L^i\left(C_{k_1}, \dots, C_{k_n}, p_1, \dots, V, h, \dots\right)$$

Le terme d'échange d'une classe $k$ pour une itération est la somme des termes d'échange pour toutes les lois où cette classe $k$ intervient :

$$ E_k=E_{k}^{i_1}+E_{k}^{i_2}+\cdots+E_{k}^{i_n}$$

Pour la résolution du Newton, $E'_k$ est calculé de la même façon en sommant sur les lois où la classe intervient.

Calcul des classes fixées en régime transitoire

David Dorchies — 2011-01-28T14:33:23Z

L'équation différentielle de variation des classes
fixées est la suivante :

$\frac{\partial C_k}{\partial t}(x,t)=E_k(x,t)$

Avec :

$C_k$ : la masse par unité
de longueur de cours d'eau de la classe au point considéré ;
$E_k$ : l'évolution de $C_k$ en fonction du temps.

La discrétisation pour le temps $j+1=t+\Delta t$ de cette équation est
faite par une méthode d'Euler semi-implicite :

$C_k^{j+1}-C_k^{j}-\frac{\Delta t}{2}\left(E_k^{j}+E_k^{j+1}\right)=0$

Sa résolution est ensuite réalisée dans les mêmes itérations de Newton
que les classes de qualité dérivantes.

La formule de l'itération de Newton pour une itération
i est :

$C_k^{j+1,i}=C_k^{j+1,i-1}-\frac{f\left(C_k^{j+1,i-1}\right)}{f'\left(C_k^{j+1,i-1}\right)}=C_k^{j+1,i-1}-\frac{b^{i}}{a^{i}}$

avec $b_{i}=C_k^{j+1,i-1}-C_k^{j}-\frac{\Delta t}{2}\left(E_k^{j}+E_k^{j+1,i}\right)$ et $a_{i}=1-\frac{\Delta t}{2}E'_k^{j+1,i}$

$E_k^{j+1,i}$ et $E'_k^{j+1,i}$ sont fonction de $C_k^{j+1,i-1}$ mais
aussi des classes de qualité dérivantes $C_{k'}^{j+1,i-1}$

Calcul du transport en régime transitoire

David Dorchies — 2011-01-28T10:39:14Z

Résolution à la section de calcul

Discrétisation avec le schéma de Preissmann

L'équation de transport (sans la diffusion) des classes dérivantes est :

$\frac{\partial \mathit{CS}}{\partial t}+\frac{\partial \mathit{CQ}}{\partial x}=E(x,t,C)S$

Les équations de discrétisation selon le schéma de Preissmann sont :

$\frac{\partial U}{\partial t}\approx (1-\psi )\frac{U_{j}^{n+1}-U_{j}^{n}}{\Delta t}+\psi \frac{U_{j+1}^{n+1}-U_{j+1}^{n}}{\Delta t}$

$\frac{\partial U}{\partial x}\approx (1-\theta )\frac{U_{j+1}^{n}-U_{j}^{n}}{\Delta x_{j+1/2}}+\theta \frac{U_{j+1}^{n+1}-U_{j}^{n+1}}{\Delta x_{j+1/2}}$
$U(x,t)=\theta \left(\psi U_{j+1}^{n+1}+(1-\psi )U_{j}^{n+1}\right)+(1-\theta )\left(\psi U_{j+1}^{n}+(1-\psi )U_{j}^{n}\right)$

Si on pose : $\Delta (U)_{i}=U_{i}^{n+1}-U_{i}^{n}$

L'application du schéma de Preissmann à
l'équation de transport pour la section
j+1 au temps k+1 donne :

$\frac{(1-\psi )\Delta (\mathit{CS})_{j}+\psi .\Delta (\mathit{CS})_{j+1}}{\Delta t}+\frac{(\mathit{CQ})_{j+1}^{k}-(\mathit{CQ})_{j}^{k}+\theta \left(\Delta (\mathit{CQ})_{j+1}-\Delta (\mathit{CQ})_{j}\right)}{\Delta x}$
$=\left[(1-\psi )(\mathit{ES})_{j}^{k}+\psi (\mathit{ES})_{j+1}^{k}\right]+\theta \left[(1-\psi )\Delta (\mathit{ES})_{j}+\psi \Delta (\mathit{ES})_{j+1}\right]$

L'équation de transport peut s'écrire
sous la forme $F(C_{j+1}^{k+1})=0$ car $C_{j}^{k+1}$ est connu par
descente de la condition amont.

$\delta x\left[(1-\psi )\Delta (\mathit{CS})_{j}+\psi .\Delta (\mathit{CS})_{j+1}\right]+\Delta t\left[(\mathit{CQ})_{j+1}^{k}-(\mathit{CQ})_{j}^{k}+\theta \left(\Delta (\mathit{CQ})_{j+1}-\Delta (\mathit{CQ})_{j}\right)\right]$
$-\delta x\Delta t\left[(1-\psi )(\mathit{ES})_{j}^{k}+\psi (\mathit{ES})_{j+1}^{k}\right]-\delta x\Delta t\theta \left[(1-\psi )\Delta (\mathit{ES})_{j}+\psi \Delta (\mathit{ES})_{j+1}\right]=0$

$\delta x\left[(1-\psi )\left((\mathit{CS})_{j}^{k+1}-(\mathit{CS})_{j}^{k}\right)+\psi .\left((\mathit{CS})_{j+1}^{k+1}-(\mathit{CS})_{j+1}^{k}\right)\right]$
$+\Delta t\left[(\mathit{CQ})_{j+1}^{k}-(\mathit{CQ})_{j}^{k}+\theta \left(\left((\mathit{CQ})_{j+1}^{k+1}-(\mathit{CQ})_{j+1}^{k}\right)-\left((\mathit{CQ})_{j}^{k+1}-(\mathit{CQ})_{j}^{k}\right)\right)\right]$
$-\delta x\Delta t\left[(1-\psi )(\mathit{ES})_{j}^{k}+\psi (\mathit{ES})_{j+1}^{k}\right]$
$-\delta x\Delta t\theta \left[(1-\psi )\left((\mathit{ES})_{j}^{k+1}-(\mathit{ES})_{j}^{k}\right)+\psi \left((\mathit{ES})_{j+1}^{k+1}-(\mathit{ES})_{j+1}^{k}\right)\right]=0$

Résolution par la méthode du gradient conjugué

L'équation $F(C_{j+1}^{k+1})=0$ est ensuite résolue par la méthode du gradient conjugué : on pose $S(C_{j+1}^{k+1})= ^{t}F(C_{j+1}^{k+1})F(C_{j+1}^{k+1})$ le carré de la norme de $F(C_{j+1}^{k+1})$. La résolution de l'équation est obtenue par minimisation de $S$. On a pour cela besoin du gradient $\nabla S$ de $S$, et donc du jacobien $J$ de $F$ qui prend en compte les dérivées partielles des termes d'échanges en fonction de chaque classe de qualité.

Discrétisation au nœud

Équation de continuité au Nœud

On intègre l'équation de conservation de la masse au nœud entre $t$ et $t+\Delta t$ :

$\int _{t}^{t+\Delta t}{\left(\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right)\mathit{dt}}-\int _{C_{\mathit{nd}}}^{C_{\mathit{nd}}+\Delta C_{\mathit{nd}}}{V_{\mathit{cas}}\mathit{dC}_{\mathit{nd}}}$
$+\int _{t}^{t+\Delta t}{\left(C_{\mathit{nd}}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]\right)\mathit{dt}}-\int _{V_{\mathit{cas}}}^{V_{\mathit{cas}}+\Delta V_{\mathit{cas}}}{C_{\mathit{nd}}\mathit{dV}_{\mathit{cas}}}=0$

On calcule cette intégrale avec la méthode des trapèzes :

$\frac{\Delta t}{2}\left(\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right)^{k}$
$+{\frac{\Delta t}{2}}\left(\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right)^{k+1}$
$-{\frac{V_{\mathit{cas}}^{k+1}-V_{\mathit{cas}}^{k}}{2}}\left[C_{\mathit{nd}}^{k}+C_{\mathit{nd}}^{k+1}\right]-\frac{C_{\mathit{nd}}^{k+1}-C_{\mathit{nd}}^{k}}{2}\left[V_{\mathit{Cas}}^{k}+V_{\mathit{Cas}}^{k+1}\right]$
$+{\frac{\Delta t}{2}}C_{\mathit{nd}}^{k}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]^{k}$
$+{\frac{\Delta t}{2}}C_{\mathit{nd}}^{k+1}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]^{k+1}=0$

L'équation est non linéaire uniquement du fait du terme source E.

Calcul pour un nœud avec casier : résolution par itération de Newton

Soit $f\left(C_{\mathit{nd}}^{i}\right)=0$,
l'intégration de l'équation de
continuité au Nœud par la méthode des trapèzes à
l'itération i.

La formule de l'itération de Newton est :
$C_{\mathit{nd}}^{i}=C_{\mathit{nd}}^{i-1}-\frac{f\left(C_{\mathit{nd}}^{i-1}\right)}{f'\left(C_{\mathit{nd}}^{i-1}\right)}=C_{\mathit{nd}}^{i-1}-\frac{b^{i}}{a^{i}}$

$b^{i}=\frac{\Delta t}{2}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})+C_{\mathit{nd}}^{k}\left(\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k}$
$+{\frac{\Delta t}{2}}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}}^{i-1})+C_{\mathit{nd}}^{i-1}\left(\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k+1}$
$-V_{\mathit{cas}}^{k+1}C_{\mathit{nd}}^{k}-V_{\mathit{cas}}^{k}C_{\mathit{nd}}^{k+1,i-1}$

Si on pose $b^{i}=b_{1}+b_{2}^{i}$ avec
$b_1$ invariant dans les itérations et $b_2$ dépendant de
l'itération i, on a :

$b_{1}=\frac{\Delta t}{2}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})+C_{\mathit{nd}}^{k}\left(\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k}$
$+{\frac{\Delta t}{2}}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}\right]^{k+1}-V_{\mathit{cas}}^{k+1}C_{\mathit{nd}}^{k}$

Si on pose :
$b_{11}^{k}=\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}\right]^{k}$
et
$a_{1}^{k}=\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}\right]^{k}$

On obtient : $b_{1}=\frac{\Delta t}{2}\left(\left[b_{11}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right]^{k}+b_{11}^{k+1}\right)+C_{\mathit{nd}}^{k}\left(\frac{\Delta t}{2}\left[a_{1}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]^{k}-V_{\mathit{cas}}^{k+1}\right)$

et $b_{2}^{i}=\frac{\Delta t}{2}\left[V_{\mathit{Cas}}E(t,C_{\mathit{nd}}^{i-1})+C_{\mathit{nd}}^{i-1}\left(a_{1}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k+1}-V_{\mathit{cas}}^{k}C_{\mathit{nd}}^{k+1,i-1}$

Pour la dérivée en fonction de
$C_{nd}$, on a :
$a^{i}=\frac{\Delta t}{2}\left[V_{\mathit{Cas}}E'(t,C_{\mathit{nd}}^{i-1})+\left(a_{1}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k+1}-V_{\mathit{cas}}^{k}$

Dans le cas où
$C_{nd}$ est un
vecteur de concentration cette équation est vérifiée pour chaque
concentration
$C_{nd}(i)$.
Les équations sont couplées uniquement par le terme source
E.

La matrice a est non diagonale uniquement du fait du
terme E'. Si ce
terme est négligeable devant les apports au nœud ou si seulement la
dérivée par rapport à la concentration $C_{nd}(i)$ est
prépondérante alors la matrice peut être supposée diagonale et on peut
résoudre équation par équation à chaque itération.

Calcul pour un nœud sans casier

Un nœud sans casier est un nœud ayant un volume (ou une surface)
de casier nulle au temps k+1 et au temps
k.

L'équation est linéaire et le résultat exact est trouvé
en une seule itération du Newton. Si les coefficients de répartition
sont hétérogènes, les
coefficients de répartition doivent être ajustés afin
d'assurer la conservation de la matière au nœud.

La formule de calcul de $C_{nd}$ devient :
$C_{\mathit{nd}}^{k+1}=C_{\mathit{nd}}^{k}-\frac{f\left(C_{\mathit{nd}}^{k}\right)}{f'\left(C_{\mathit{nd}}^{k}\right)}=C_{\mathit{nd}}^{k}-\frac{b}{a}$

Les équations du calcul pour un nœud avec casier deviennent :
$b=\frac{\Delta t}{2}\left[b_{11}^{k}+b_{11}^{k+1}+C_{\mathit{nd}}^{k}\left(a_{1}^{k}+a_{1}^{k+1}\right)\right]$
et $a=\frac{\Delta t}{2}a_{1}^{k+1}$

Et donc :
$C_{\mathit{nd}}^{k+1}=C_{\mathit{nd}}^{k}-\frac{b_{11}^{k}+b_{11}^{k+1}+C_{\mathit{nd}}^{k}\left(a_{1}^{k}+a_{1}^{k+1}\right)}{a_{1}^{k+1}}$

Si on prend en compte le coefficient ajustable dans la formule $a_{1}^{k}=\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}\right]^{k}$ , on obtient :

$a_{1}^{k}=\left[\sum _{s_{b}Q_{b}<0,b_{i}=0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{i}=0}k_{p}Q_{p}+k_{a}\left(\sum _{s_{b}Q_{b}<0,b_{i}=1}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{i}=1}k_{p}Q_{p}\right)\right]^{k}$

Calcul du transport en régime permanent

David Dorchies — 2011-01-28T09:52:01Z

Calcul du transport dans un bief de section à section

Équation de transport

En régime permanent, l'équation de transport d'une classe dérivante devient :

$\frac{\mathit{dCQ}}{\mathit{dx}}=E(x,C)S$

Intégration par la méthode de Runge-Kutta d'ordre 4

L'intégration de l'équation précédente se
fait par la méthode de Runge-Kutta d'ordre 4 (RK4).

Si on pose $y=\mathit{CQ}$ , on a
$\frac{\mathit{dy}}{\mathit{dx}}=f\left(y,x\right)=E\left(x,y\right)S$
, la méthode de Runge-Kutta donne pour la section
j+1 : $y_{j+1}=y_{j}+\frac{\Delta x}{6}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)$ avec :

$k_{1}=f\left(y_{j},x_{j}\right)=E\left(x_{j},y_{j}\right)S_{j}=E\left(x_{j},C_{j}\right)S_{j}$
$k_{2}=f\left(y_{j}+\frac{\Delta x}{2}k_{1,}x_{j+\frac{1}{2}}\right)=E\left(x_{j+\frac{1}{2}},y_{j}+\frac{\Delta x}{2}k_{1}\right)S_{j+\frac{1}{2}}=E\left(x_{j+\frac{1}{2}},\frac{y_{j}+\frac{\Delta x}{2}k_{1}}{Q_{j+\frac{1}{2}}}\right)S_{j+\frac{1}{2}}$ avec $S_{j+\frac{1}{2}}=\frac{S_{j}+S_{j+1}}{2}$ et
$Q_{j+\frac{1}{2}}=\frac{Q_{j}+Q_{j+1}}{2}$
$k_{3}=f\left(y_{j}+\frac{\Delta x}{2}k_{2,}x_{j+\frac{1}{2}}\right)=E\left(x_{j+\frac{1}{2}},y_{j}+\frac{\Delta x}{2}k_{2}\right)S_{j+\frac{1}{2}}=E\left(x_{j+\frac{1}{2}},\frac{y_{j}+\frac{\Delta x}{2}k_{2}}{Q_{j+\frac{1}{2}}}\right)S_{j+\frac{1}{2}}$
$k_{4}=f\left(y_{j}+\Delta xk_{3,}x_{j+1}\right)=E\left(x_{j+1},y_{j}+\Delta xk_{3}\right)S_{j+1}=E\left(x_{j+1},\frac{y_{j}+\Delta xk_{3}}{Q_{j+1}}\right)S_{j+1}$

Les concentrations à la section j+1 sont alors
égales à : $C_{j+1}=\frac{y_{j+1}}{Q_{j+1}}$

Calcul de la conservation au nœud

Équation de continuité au nœud

En régime permanent, l'équation de continuité au nœud devient :

$\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})+C_{\mathit{nd}}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]=0$

Résolution par la méthode de Newton

Soit $f\left(C_{\mathit{nd}}^{i}\right)$ , l'équation
de continuité au nœud à l'itération i.

La formule de l'itération de Newton est :

$C_{\mathit{nd}}^{i}=C_{\mathit{nd}}^{i-1}-\frac{f\left(C_{\mathit{nd}}^{i-1}\right)}{f'\left(C_{\mathit{nd}}^{i-1}\right)}=C_{\mathit{nd}}^{i-1}-\frac{b^{i}}{a^{i}}$

$b^{i}=\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}}^{i-1})$
$+C_{\mathit{nd}}^{i-1}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]$

Dans le cas d'un nœud sans casier, il
n'y a pas de terme source et le problème est donc
linéaire. Avec des coefficients de répartition hétérogènes, les coefficients
$k_{b}$, $k_{p}$ et $k_{inf}$ sont
ajustés avant l'unique itération de calcul.

Dans le cas d'un nœud avec casier, le problème
n'est plus linéaire du fait du terme source et
nécessite plusieurs itérations cependant les coefficients
$k_{b}$, $k_{p}$ et $k_{inf}$ ne sont
pas ajustés. Ce qui nous permet d'avoir pour le cas
général :

$a^{i}=V_{\mathit{Cas}}E'(C_{\mathit{nd}}^{i-1})+\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}$
avec $E'(C_{\mathit{nd}}^{i-1})$ la dérivée de
l'échange à chaque itération.

Calcul pour un nœud avec casier

On sépare les termes invariants dans l'itération en
posant :
$b_{1}=\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}$
et $a_{1}=\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}$

Ce qui permet d'obtenir dans le cas
d'un nœud avec casier :

$b^{i}=b_{1}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}}^{i-1})+a_{1}C_{\mathit{nd}}^{i-1}$
et $a^{i}=V_{\mathit{Cas}}E'(C_{\mathit{nd}}^{i-1})+a_{1}$

Calcul pour un nœud sans casier

Pour un nœud sans casier, la concentration au nœud
s'obtient directement en moyennant les concentrations
d'apport au nœud :

$C_{\mathit{nd}}=\frac{\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}}{\sum _{s_{b}Q_{b}>0}s_{b}Q_{b}+\sum _{Q_{p}>0}Q_{p}}$

Dans le cas d'un nœud sans casier avec départs hétérogènes, les coefficients de répartition doivent être ajustés afin d'assurer la conservation de la matière au nœud. Il s'agit de calculer Ka avant de calculer la concentration dans la première section du bief aval comme c'est décrit au chapitre "Calcul du réseau".

Calcul du réseau

David Dorchies — 2011-01-27T17:28:13Z

Dans le cas où il n'y a pas de débits négatifs, on calcule le réseau dans l'ordre amont-aval « topologique » calculé par Talweg.

Arrivé à un nœud on connait toutes les concentrations
$C_{p}$ et $C_{b}$ des prises et des biefs confluents, et avec
l'équation au nœud on est capable de calculer
$C_{nd}$ et donc $C_{b}$ pour tous
les biefs avec débit défluents.

La formule de calcul de $C_{b}$ pour les
biefs défluents diffère si le coefficient de répartition est ajustable
ou non. Pour un coefficient non ajustable, on aura :
$C_{b}=C_{\mathit{nd}}k_{b}$ et pour un coefficient ajustable on aura :
$C_{b}=C_{\mathit{nd}}k_{b}k_{a}$.

Dans le cas où il existe des débits négatifs dans les biefs, le calcul nécessite un reclassement amont-aval dans le sens « hydraulique » des biefs et des sections. Cette fonctionnalité n'est pas gérée actuellement.

Concepts et équations au nœud

David Dorchies — 2011-01-27T16:51:01Z

On traite le cas du nœud dans sa plus grande généralité. Un nœud peut être constitué :

d'un casier ;
de prises avec débits entrants ou sortants ;
de biefs confluents ou défluents.

Les concepts et variables intervenant au nœud sont résumés dans le schéma suivant (seules les classes de qualité dérivantes sont considérées) :

Un casier est connu sous la forme d'une suite monotone
de couples (surface, cote). Ce casier peut être soumis à
l'évaporation connue sous la forme
d'une vitesse d'évaporation
$V_{ep}(t)$ et à une vitesse d'infiltration sous la forme $v_{inf}(Z_{nd})$,
$Z_{nd}$ étant la cote dans le casier.

Le Nœud reçoit ou perd du débit via les prises
($Q_p$) et les biefs ($Q_b$). Les débits entrants dans le nœud sont positifs, et les sortants négatifs.

Les biefs sont identifiés par un sens d'écoulement lié
à leur topologie. Ce sens « topologique » d'écoulement
est invariable dans le temps :

$s_{b}=1$ pour les biefs confluents
$s_{b}=-1$ pour les biefs défluents

Le débit $Q_b$ est
positif s'il coule dans le sens
d'écoulement du bief et négatif s'il
coule dans le sens inverse. On obtient le débit entrant ou sortant du
nœud avec le bon signe par la formule : $s_{b}Q_{b}$

Débits entrants dans le nœud (i.e. débits positifs)

Pour les biefs confluents ou défluents avec , $s_{b}Q_{b}>0$ ou les
prises avec $Q_{p}>0$ , on a une condition limite
$C$ qui est la concentration d'apport dans le nœud.

Débits sortants du nœud (i.e. débits négatifs)

Biefs

Pour les biefs confluents ou défluents avec , $s_{b}Q_{b}<0$ , on a
$C_{b}=k_{b}C_{nd}$ .

Prises

Pour les prises avec $Q_{p}<0$ , on a $C_{p}=k_{p}C_{nd}$.

Infiltration au nœud

On a $C_{\mathit{inf}}=k_{\mathit{inf}}C_{\mathit{nd}}$ et
$Q_{\mathit{inf}}=S_{\mathit{cas}}v_{\mathit{inf}}$ avec :
$S_{cas}$ la surface
du casier, $v_{inf}$
la vitesse d'infiltration,
$C_{\mathit{nd}}$ la
concentration au Nœud.

Débits de stockage au casier

On considère le casier comme un bassin à connexion directe au noeud,
la concentration du noeud est donc la même que celle du casier. La
quantité de matière présente dans le noeud-casier est :
$m_{\mathit{Cas}}=V_{\mathit{Cas}}C_{\mathit{nd}}$ et le débit entrant
ou sortant du noeud s'obtient avec :
$Q_{\mathit{cas}}=S_{\mathit{cas}}\frac{\mathit{dZ}_{\mathit{nd}}}{\mathit{dt}}$.

Équation du bilan de masse au nœud

Écrivons l'équation de conservation de la masse
d'une classe de qualité au nœud $m_{\mathit{nd}}$
à l'instant t du calcul en tenant
compte d'un terme source au nœud.

Au temps t on connaît :

la surface $S_{cas}$ du casier, son volume $V_{cas}$ ainsi que
sa cote $Z_{nd}$,
les débits entrants et sortants du nœud,
les concentrations d'apport liées aux débits entrants dans le nœud,
les coefficients de répartition liés aux débits sortants du nœud.

$\frac{\mathit{dm}_{\mathit{nd}}}{\mathit{dt}}=\sum_{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum_{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}C_{\mathit{nd}}$
$+\sum_{Q_{p}>0}C_{p}Q_{p}+\sum_{Q_{p}<0}k_{p}Q_{p}C_{\mathit{nd}}-k_{\mathit{inf}}S_{\mathit{cas}}C_{\mathit{nd}}v_{\mathit{inf}}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})$

Avec $\frac{\mathit{dm}_{\mathit{nd}}}{\mathit{dt}}=\frac{d\left(C_{\mathit{nd}}V_{\mathit{cas}}\right)}{\mathit{dt}}=\frac{dV_{\mathit{cas}}}{\mathit{dt}}C_{\mathit{nd}}+\frac{dC_{\mathit{nd}}}{\mathit{dt}}V_{\mathit{cas}}$,
l'équation de conservation de la masse au nœud devient :

$\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}\left(E(t,C_{\mathit{nd}})-\frac{\mathit{dC}_{\mathit{nd}}}{\mathit{dt}}\right)$
$+C_{\mathit{nd}}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}-\frac{\mathit{dV}_{\mathit{cas}}}{\mathit{dt}}\right]=0$

Contraintes au nœud pour les k

3 cas sont envisageables pour les coefficients de répartition de concentration sur les départ de matière du nœud :

Tous les départs sont homogènes

Tous les k valent alors 1 : la conservation de la masse au nœud est assurée par l'équation de conservation de la masse au nœud quelque soit les variations de débits et de concentration.

Les départs sont hétérogènes et il y a un casier au nœud

La masse de polluant entrant ou sortant du casier assure la conservation de la masse au nœud en permanent et en transitoire.

Les départs sont hétérogènes et il n'y a pas de casier au nœud

Il faut ajuster la valeur d'au moins un coefficient de répartition pour assurer la continuité au nœud avec la contrainte que ce coefficient reste positif.

L'utilisateur défini les coefficients de répartition aux prises, au bief et à l'infiltration. Pour chacun de ces coefficients, il définit aussi si ceux-ci sont fixes ou ajustable pour assurer la continuité au nœud (Sauf pour l'infiltration car elle n'a lieu que dans un nœud avec casier et il n'y a donc pas d'ajustement dans ce cas).

La continuité s'effectue en introduisant un coefficient
d'ajustement $k_{a}$ tel que chaque coefficient de répartition puisse s'écrire
ainsi : ${k_{a}}^{b_{i}}k_{i}$ où $b_{i}$ est le choix pour le
coefficient $k_{i}$ entre rester fixe $b_{i}=0$ et être ajustable
$b_{i}=1$ pour assurer la conservation de la masse. De ce fait, les
deux possibilités sont représentées de cette façon :

Coefficient de répartition fixe : $b_{i}=0$ et
${k_{a}}^{b_{i}}k_{i}=k_{i}$
Coefficient de répartition ajustable : $b_{i}=1$ et
${k_{a}}^{b_{i}}k_{i}=k_{a}k_{i}$

Si on simplifie l'écriture de
l'équation de conservation de la masse au nœud en assimilant les débits aux biefs et aux prises sous forme de débits entrants $Q_{i}>0$ et de débits
sortants $Q_{i}<0$ , on obtient :

$\sum _{Q_{i}>0}C_{i}Q_{i}+C_{\mathit{nd}}\left[\sum _{Q_{i}<0,b_{i}=0}k_{i}Q_{i}+k_{a}\sum _{Q_{i}<0,b_{i}=1}k_{i}Q_{i}\right]=0$ et donc
$k_{a}=-{\frac{\frac{1}{C_{\mathit{nd}}}\sum _{Q_{i}>0}C_{i}Q_{i}+\sum _{Q_{i}<0,b_{i}=0}k_{i}Q_{i}}{\sum _{Q_{i}<0,b_{i}=1}k_{i}Q_{i}}}$

Or pour un nœud sans casier, la concentration au nœud peut
s'écrire :

$C_{\mathit{nd}}=\frac{\sum _{Q_{i}>0}C_{i}Q_{i}}{\sum _{Q_{i}>0}Q_{i}}$

L'équation de $k_{a}$ devient :
$k_{a}=-{\frac{\sum _{Q_{i}>0}Q_{i}+\sum _{Q_{i}<0,b_{i}=0}k_{i}Q_{i}}{\sum _{Q_{i}<0,b_{i}=1}k_{i}Q_{i}}}$

L'équation complète du calcul de $k_{a}$ donne :

$k_{a}=-{\frac{\left(\sum _{s_{b}Q_{b}>0}s_{b}Q_{b}+\sum _{Q_{p}>0}Q_{p}\right)+\left(\sum _{s_{b}Q_{b}<0,b_{b}=0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{p}=0}k_{p}Q_{p}\right)}{\sum _{s_{b}Q_{b}<0,b_{b}=1}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{p}=1}k_{p}Q_{p}}}$

L'équation de conservation de la masse au nœud devient alors :

$\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+C_{\mathit{nd}}\left(\sum _{s_{b}Q_{b}<0}{k_{a}}^{b_{b}}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}{k_{a}}^{b_{p}}k_{p}Q_{p}\right)=0$

Afin d'éviter de placer trop
d'exponentielles qui pourraient pénaliser lourdement
les calculs, on pourra aussi l'écrire
cette équation de cette façon :

$\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}$
$+C_{\mathit{nd}}\left[\sum _{s_{b}Q_{b}<0,b_{b}=0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{p}=0}k_{p}Q_{p}+k_{a}\left(\sum _{s_{b}Q_{b}<0,b_{b}=1}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{p}=1}k_{p}Q_{p}\right)\right]=0$

Évolution de la température

David Dorchies — 2011-01-27T08:47:55Z

Définition

Cette loi fait évoluer la température en fonction de la température de l'eau et des données météorologiques selon les formules développées dans l'article « Simulation de la température ».

Caractéristiques

Identifiant de la loi : 101
Nombre de classes intervenant : 1
Nombre de paramètres : 0

Classes intervenantes :

Température de l'eau