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 x2k_1,x_{j+\frac12}\right)=E\left(x_{j+\frac12},y_j+\frac\Delta x2k_1\right)S_{j+\frac12}=E\left(x_{j+\frac12},\frac{y_j+\frac\Delta x2k_1}{Q_{j+\frac12}}\right)S_{j+\frac12}\) where \(S_{j+\frac12}=\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 x2k_2,x_{j+\frac12}\right)=E\left(x_{j+\frac12},y_j+\frac\Delta x2k_2\right)S_{j+\frac12}=E\left(x_{j+\frac12},\frac{y_j+\frac\Delta x2k_2}{Q_{j+\frac12}}\right)S_{j+\frac12}\)
• \(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_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}}\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]=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_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]\)

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_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_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_bQ_b<0}k_bs_bQ_b+\sum _{Q_p<0}k_pQ_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_bQ_b>0}C_bs_bQ_b+\sum _{Q_p>0}C_pQ_p}{\sum _{s_bQ_b>0}s_bQ_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”.