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}$