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_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_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_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_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 t2\left(\sum _{s_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_p+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right)^k\) \(+{\frac\Delta t2}\left(\sum _{s_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_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 t2}C_{\mathit{nd}}^k\left[\sum _{s_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_p-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]^k\) \(+{\frac\Delta t2}C_{\mathit{nd}}^{k+1}\left[\sum _{s_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_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 t2\left[\sum _{s_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_p+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})+C_{\mathit{nd}}^k\left(\sum _{s_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_p-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^k\) \(+{\frac\Delta t2}\left[\sum _{s_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_p+V_{\mathit{Cas}}E(t,C_{\mathit{nd}}^{i-1})+C_{\mathit{nd}}^{i-1}\left(\sum _{s_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_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 t2\left[\sum _{s_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_p+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})+C_{\mathit{nd}}^k\left(\sum _{s_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_p-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^k\) \(+{\frac\Delta t2}\left[\sum _{s_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_p\right]^{k+1}-V_{\mathit{cas}}^{k+1}C_{\mathit{nd}}^k\)

Si on pose : \(b_11^k=\left[\sum _{s_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_p\right]^k\) et \(a_1^k=\left[\sum _{s_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_p\right]^k\)

On obtient : \(b_1=\frac\Delta t2\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 t2\left[a_1-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]_{\mathit{cas}}^{k+1,k-V}\right)\)

et \(b_2^i=\frac\Delta t2\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}}^kC_{\mathit{nd}}^{k+1},i-1\)

Pour la dérivée en fonction de \(C_{nd}\), on a : \(a^i=\frac\Delta t2\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 t2\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_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_p\right]^k\) , on obtient :

\(a_1^k=\left[\sum _{s_bQ_b<0,b_i=0}k_bs_bQ_b+\sum _{Q_p<0,b_i=0}k_pQ_p+k_a\left(\sum _{s_bQ_b<0,b_i=1}k_bs_bQ_b+\sum _{Q_p<0,b_i=1}k_pQ_p\right)\right]^k\)