Calcul du transport en régime transitoire

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}