Calculation at a section

Discretization for Preissmann Scheme

The advection equation (without diffusion) for drift classes is :

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

Preissmann scheme discretization equations are :

$\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)$

We note : $\Delta (U)_{i}=U_{i}^{n+1}-U_{i}^{n}$

These equations for section j+1 at time k+1 yield :

$\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]$

The transport equation can be rewritten as $F(C_{j+1}^{k+1})=0$ since $C_{j}^{k+1}$ is given by the upstream condition.

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

Resolution by conjugate gradient method

This equation $F(C_{j+1}^{k+1})=0$ is solved using a conjugate gradient method, requiring the computation of a Jacobian matrix.

Let $S(C_{j+1}^{k+1})= ^{t}F(C_{j+1}^{k+1})F(C_{j+1}^{k+1})$. The equation is solved by minimizing $S$, using the conjugate gradient method. One need to calculate the partial derivatives of exchange terms with respect to every quality class in order to calculate $F$’s Jacobian matrix $J$ which gives $S$’s gradient $\nabla S$.

Discretization at a node

Mass balance at a node

The mass balance equation at a node is integrated between $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$

This integral is calculated using the trapezoidal rule formula :

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

The source term E is the only non-linear term in this equation.

Newton’s method for a node with a pond

Let $f\left(C_{\mathit{nd}}^{i}\right)=0$ be the result of the integration at the ith iteration.

Newton’s method gives :
$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}$

If $b^{i}=b_{1}+b_{2}^{i}$ where $b_1$ is constant over iterations et $b_2$ we get :

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

We define
$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}$
and
$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}$

yielding : $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)$

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

The derivative with respect to $C_{nd}$ is :
$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}$

As $C_{nd}$ is a vector, this equation is true for every quality class $C_{nd}(i)$.

The different equations may be interdependent due to the exchange rate E.

As a consequence, matrix a is non-diagonal due to E’. If this term is negligible compared to the others, the Newton’s method can still solve the system.

Calculation at a node without a pond

A node without a pond is a nodes where the volume (or area) of the pod is zero at time k+1 and k.

The equation is then linear and the exact result is found within one Newton iteration. If distribution coefficients are heterogeneous, they are adjusted to ensure mass balance at the node.

The formula for $C_{nd}$ is :
$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}$

Equations for a node with a pond are :
$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}$

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

Given the adjustable coefficient :
$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}$ , we get :

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