Steady state calculation

Calculation at a section

Transport equation

The steady-state advection equation is:

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

Integration using fourth-order Runge-Kutta méthod

Previous equation is rewriten as $\frac{\mathit{dy}}{\mathit{dx}}=f\left(y,x\right)$ (where $y=\mathit{CQ}$ and $f\left(y,x\right)=E\left(x,y\right)S$) and integrated using a fourth-order Runge-Kutta méthod(RK4) :

$y_{j+1}=y_{j}+\frac{\Delta x}{6}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)$ where :

  • $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}}$ where $S_{j+\frac{1}{2}}=\frac{S_{j}+S_{j+1}}{2}$ and
    $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}$

Concentrations in section j+1 are therefore :

$C_{j+1}=\frac{y_{j+1}}{Q_{j+1}}$

Calculation at a node

Continuity equation at a node

The steady state mass balance at a node is :

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

Newton’s method

Let $f\left(C_{\mathit{nd}}^{i}\right)=0$ be the mass balance equation for the node on iteration i.

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

Is the node does not contain a pond, there is no source term and the equation is linear. If distribution coefficients $k_{b}$, $k_{p}$ and $k_{inf}$ are heterogeneous, they are adjusted before the first and only iteration.

If there is a pond, the equation is no longer linear and several iterations are needed but distribution coefficients are not adjusted.

The general equation for $a^{i}$ is :

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

where $E’(C_{\mathit{nd}}^{i-1})$ is the derivative of the exchange.

Calculation for a node with a pond

Some terms in the previous equation are constant and do not need to be recomputed on every iteration :

$b_{1}=\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}$
$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}}$

Which yields :

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

Calculation for a node without a pond

For a node without a pond, the node concentration is directly calculated as a pondered mean of inflow concentrations.

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

If the distribution coefficients are adjustable, $K_a$ is computed before the iterations, as described in "Network computation".