Transient calculation
Calculation at a section
Discretization for Preissmann Scheme
The advection equation (without diffusion) for drift classes is :
∂CS∂t+∂CQ∂x=E(x,t,C)S
Preissmann scheme discretization equations are :
∂U∂t≈(1−ψ)Un+1j−UnjΔt+ψUn+1j+1−Unj+1Δt
∂U∂x≈(1−θ)Unj+1−UnjΔxj+1/2+θUn+1j+1−Un+1jΔxj+1/2
U(x,t)=θ(ψUn+1j+1+(1−ψ)Un+1j)+(1−θ)(ψUnj+1+(1−ψ)Unj)
We note : Δ(U)i=Un+1i−Uni
These equations for section j+1 at time k+1 yield :
(1−ψ)Δ(CS)j+ψ.Δ(CS)j+1Δt+(CQ)kj+1−(CQ)kj+θ(Δ(CQ)j+1−Δ(CQ)j)Δx
=[(1−ψ)(ES)kj+ψ(ES)kj+1]+θ[(1−ψ)Δ(ES)j+ψΔ(ES)j+1]
The transport equation can be rewritten as F(Ck+1j+1)=0 since Ck+1j is given by the upstream condition.
δx[(1−ψ)Δ(CS)j+ψ.Δ(CS)j+1]+Δt[(CQ)kj+1−(CQ)kj+θ(Δ(CQ)j+1−Δ(CQ)j)]
−δxΔt[(1−ψ)(ES)kj+ψ(ES)kj+1]−δxΔtθ[(1−ψ)Δ(ES)j+ψΔ(ES)j+1]=0
δx[(1−ψ)((CS)k+1j−(CS)kj)+ψ.((CS)k+1j+1−(CS)kj+1)]
+Δt[(CQ)kj+1−(CQ)kj+θ(((CQ)k+1j+1−(CQ)kj+1)−((CQ)k+1j−(CQ)kj))]
−δxΔt[(1−ψ)(ES)kj+ψ(ES)kj+1]
−δxΔtθ[(1−ψ)((ES)k+1j−(ES)kj)+ψ((ES)k+1j+1−(ES)kj+1)]=0
Resolution by conjugate gradient method
This equation F(Ck+1j+1)=0 is solved using a conjugate gradient method, requiring the computation of a Jacobian matrix.
Let S(Ck+1j+1)=tF(Ck+1j+1)F(Ck+1j+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 ∇S.
Discretization at a node
Mass balance at a node
The mass balance equation at a node is integrated between t et t+Δt:
∫t+Δtt(∑sbQb>0CbsbQb+∑Qp>0CpQp+VCasE(t,Cnd))dt−∫Cnd+ΔCndCndVcasdCnd
+∫t+Δtt(Cnd[∑sbQb<0kbsbQb+∑Qp<0kpQp−kinfScasvinf])dt−∫Vcas+ΔVcasVcasCnddVcas=0
This integral is calculated using the trapezoidal rule formula :
Δt2(∑sbQb>0CbsbQb+∑Qp>0CpQp+VCasE(t,Cnd))k
+Δt2(∑sbQb>0CbsbQb+∑Qp>0CpQp+VCasE(t,Cnd))k+1
−Vk+1cas−Vkcas2[Cknd+Ck+1nd]−Ck+1nd−Cknd2[VkCas+Vk+1Cas]
+Δt2Cknd[∑sbQb<0kbsbQb+∑Qp<0kpQp−kinfScasvinf]k
+Δt2Ck+1nd[∑sbQb<0kbsbQb+∑Qp<0kpQp−kinfScasvinf]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(Cind)=0 be the result of the integration at the ith iteration.
Newton’s method gives :
Cind=Ci−1nd−f(Ci−1nd)f′(Ci−1nd)=Ci−1nd−biai
bi=Δt2[∑sbQb>0CbsbQb+∑Qp>0CpQp+VCasE(t,Cnd)+Cknd(∑sbQb<0kbsbQb+∑Qp<0kpQp−kinfScasvinf)]k
+Δt2[∑sbQb>0CbsbQb+∑Qp>0CpQp+VCasE(t,Ci−1nd)+Ci−1nd(∑sbQb<0kbsbQb+∑Qp<0kpQp−kinfScasvinf)]k+1
−Vk+1casCknd−VkcasCk+1,i−1nd
If bi=b1+bi2 where b1 is constant over iterations et b2 we get :
b1=Δt2[∑sbQb>0CbsbQb+∑Qp>0CpQp+VCasE(t,Cnd)+Cknd(∑sbQb<0kbsbQb+∑Qp<0kpQp−kinfScasvinf)]k
+Δt2[∑sbQb>0CbsbQb+∑Qp>0CpQp]k+1−Vk+1casCknd
We define
bk11=[∑sbQb>0CbsbQb+∑Qp>0CpQp]k
and
ak1=[∑sbQb<0kbsbQb+∑Qp<0kpQp]k
yielding : b1=Δt2([b11+VCasE(t,Cnd)]k+bk+111)+Cknd(Δt2[a1−kinfScasvinf]k−Vk+1cas)
and bi2=Δt2[VCasE(t,Ci−1nd)+Ci−1nd(a1−kinfScasvinf)]k+1−VkcasCk+1,i−1nd
The derivative with respect to Cnd is :
ai=Δt2[VCasE′(t,Ci−1nd)+(a1−kinfScasvinf)]k+1−Vkcas
As Cnd is a vector, this equation is true for every quality class Cnd(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 Cnd is :
Ck+1nd=Cknd−f(Cknd)f′(Cknd)=Cknd−ba
Equations for a node with a pond are :
b=Δt2[bk11+bk+111+Cknd(ak1+ak+11)]
et a=Δt2ak+11
And :
Ck+1nd=Cknd−bk11+bk+111+Cknd(ak1+ak+11)ak+11
Given the adjustable coefficient :
ak1=[∑sbQb<0kbsbQb+∑Qp<0kpQp]k , we get :
ak1=[∑sbQb<0,bi=0kbsbQb+∑Qp<0,bi=0kpQp+ka(∑sbQb<0,bi=1kbsbQb+∑Qp<0,bi=1kpQp)]k