Processing math: 100%

Transient calculation

Calculation at a section

Discretization for Preissmann Scheme

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

CSt+CQx=E(x,t,C)S

Preissmann scheme discretization equations are :

Ut(1ψ)Un+1jUnjΔt+ψUn+1j+1Unj+1Δt

Ux(1θ)Unj+1UnjΔxj+1/2+θUn+1j+1Un+1jΔxj+1/2
U(x,t)=θ(ψUn+1j+1+(1ψ)Un+1j)+(1θ)(ψUnj+1+(1ψ)Unj)

We note : Δ(U)i=Un+1iUni

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))dtCnd+ΔCndCndVcasdCnd
+t+Δtt(Cnd[sbQb<0kbsbQb+Qp<0kpQpkinfScasvinf])dtVcas+Δ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+1casVkcas2[Cknd+Ck+1nd]Ck+1ndCknd2[VkCas+Vk+1Cas]
+Δt2Cknd[sbQb<0kbsbQb+Qp<0kpQpkinfScasvinf]k
+Δt2Ck+1nd[sbQb<0kbsbQb+Qp<0kpQpkinfScasvinf]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=Ci1ndf(Ci1nd)f(Ci1nd)=Ci1ndbiai

bi=Δt2[sbQb>0CbsbQb+Qp>0CpQp+VCasE(t,Cnd)+Cknd(sbQb<0kbsbQb+Qp<0kpQpkinfScasvinf)]k
+Δt2[sbQb>0CbsbQb+Qp>0CpQp+VCasE(t,Ci1nd)+Ci1nd(sbQb<0kbsbQb+Qp<0kpQpkinfScasvinf)]k+1
Vk+1casCkndVkcasCk+1,i1nd

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<0kpQpkinfScasvinf)]k
+Δt2[sbQb>0CbsbQb+Qp>0CpQp]k+1Vk+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[a1kinfScasvinf]kVk+1cas)

and bi2=Δt2[VCasE(t,Ci1nd)+Ci1nd(a1kinfScasvinf)]k+1VkcasCk+1,i1nd

The derivative with respect to Cnd is :
ai=Δt2[VCasE(t,Ci1nd)+(a1kinfScasvinf)]k+1Vkcas

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=Ckndf(Cknd)f(Cknd)=Ckndba

Equations for a node with a pond are :
b=Δt2[bk11+bk+111+Cknd(ak1+ak+11)]
et a=Δt2ak+11

And :
Ck+1nd=Ckndbk11+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