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".