[hchen]

Hi,

I just checked the codes for Crank Nicolson temporal schemes in OpenFOAM 2.4.0, and I can not understand how it is implemented. I attach the code below:

template<class Type>
tmp<fvmatrix<Type>>
CrankNicolsonDdtScheme<Type>::fvmDdt
(
....
)

...

fvm.diag () = rDtCoef*rho.internalField()*mesh().V();
vf.oldTime().oldTime();
rho.oldTime().oldTime();

if(mesh().moving())
{
...
}
else
{
    if (evaluate (ddt0))
    {
    ...
    }

    fvm.source() = 
    (
        rDtCoef*rho.oldTime().internalField()*vf.oldTime().internalField()
        + offCentre_(ddt0.internalField())
    ) * mesh().V();

    return tfvm;
}

I can see that the main difference with Euler scheme is the last term: offCentre_(ddt0.internalField).

My question is:

1. Is this ddt0 the time derivative of (rho, U) for previous time step?

2. My understanding on CN scheme is that the spacial discretization should be half implicit and half explicit, i.e. 0.5*fvm::div(rho,U) + 0.5 * fvc::div(rho,U).
(Maybe the diffusion term should be always evaluated implicitly with explicit non-orthoganity correction?)

3. why the code use ddt0 for CN scheme?

[stevemoore1981]

Hi

I've recently come across this same problem, namely trying to understand how Crank-Nicolson is implemented in OpenFOAM. I don't have an answer, but thought I might add what I have found.. maybe someone has figured this out by now?

While googling around I came across the book "The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab", F. Moukalled, L. Mangani, M.Darwish, which describes the implementation as an explicit procedure, which for an ODE of the form dphi/dt=f(phi) would look like:

phi^n+1 = phi^n-1 + 2*Delta_t f(phi^n)

Which was quite different to my understanding of the method, which would look like:

phi^n+1 = phi^n + Delta_t/2 ( f(phi^n+1) + f(phi^n) )

The book claims that this is how the CN method is implemented in OpenFOAM, and references the paper:

Moukalled F, Darwish M (2012) Transient schemes for capturing interfaces of free-surface flows. Numer Heat Transf Part B Fundam 61(3):171–203

which describes the CN method in the same way, but references the paper by Crank and Nicolson:

J. Crank and P. Nicolson, A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat-Conduction Type, Proc. Cambr. Phil.
Soc., vol. 43, pp. 50–67, 1947.

which as far as I can see presents the method in the way that I was used to. This book mentions that the CN method is implemented in OpenFOAM as a two step procedure, with an implicit solution step, followed by an explicit update step... but I found this confusing as well because my understanding was that the fvm::ddt operator was just going to basically add the contribution of the unsteady term into a linear system of equations (along with the div, laplacian, etc operators) in a generic way that works in the same way for the other time marching schemes, and so I can't see how a two step procedure could be implemented.

Anyway, if anyone knows this answer to this question, I'd love to know... it's been bugging me for a while now!

[andyru]

hello,

I faced the same problem. I wanted to know, how the implementation of CN with off-centre=psi value is working.

Well, I took paper and pencil and the implementation is quit straight forward.

1.) Assume you have a convection diffusion problem with unsteady term for temperature and you blend between phi=0 (EE=explicit euler), phi=0.5 (CN) and phi=1.0 (IE)
you gain the well known relation:

(ac' + phi * ac) * Phi_c + phi*sum(af * Phi_F) = ac'° * Phi_c° + (1 - phi) * (-ac) * Phi_c° + (1 - phi) * sum(af * Phi_F + b

Ok.
Phi_c is derivate at cell
Phi_F are neighbour values
ac' is central coefficient for cell
ac'° old time step central coefficient for cell
ac is central coefficient for central cell of contribution by convection/diffusion
af are coefficients for face flux of neighbouring values
phi is blending factor like found in literature (explanation above)
b is right side source
Phi_c° old time value of cell

Trick:
Now you multiplicate this equation with (1+psi) where psi is off-centre coefficient of openfoam (0=IE, 1=CN)
Additionally, you can replace phi with relation phi=1/(1+psi)

You gain the following equation:
(1 + psi ) * ac' * Phi_c + ac * Phi_c + sum(af * Phi_F) = b + (1 + psi) * ac'° * Phi_c° ...
This is the matrix summed up for fvm. With CN: Two times the contribution of time derivate matrix + convection = right side...
The rest is (I splitted the (1 + psi) * b !):
... = ... + psi * b - psi * ac * Psi_c° - psi * sum(af * Psi_F°)
This is the source (explicit on right hand side)
This is equal to explicit euler:
psi * (Phi_c - Phi_c°)/(dt) = - psi * L(Phi_c°)
L(Phi_c°) is spatial derivation

Ok. For making it easier now: CN is IE+EE -> psi = 1
So the following relation is true (and this is the most important step):
1 * (Phi_c - Phi_c°)/(dt) = - 1 * L(Phi_c°) = 1 * (Phi_c° - Phi_c°°)/(dt)

So you have standard matrix addition for solving equals the implicit part of the explicit source:
fvm.source() = (1 + psi) * ac'° * Phi_c° + (psi * (Phi_c° - Phi_c°°)/(dt)*Vc)
==rDtCoef*vf.oldTime().internalField()
+ offCentre_(ddt0.internalField())
)*mesh().V();

In summary: multiplicate equation with (1+psi) and replace phi with 1/(1+psi)
Right side is replaced with explicit with t-t° =implict with t° - t°°
You keep with this trick the generic matrix addition approach for fvm calculus.

[stevemoore1981]

Hi

Thanks a lot for your detailed response, it was very helpful. The only bit that I don't quite yet understand fully is the step that you mentioned:

"So the following relation is true (and this is the most important step):
1 * (Phi_c - Phi_c°)/(dt) = - 1 * L(Phi_c°) = 1 * (Phi_c° - Phi_c°°)/(dt)"

To me it seems that

1 * (Phi_c - Phi_c°)/(dt) = - 1 * L(Phi_c°)

is explicit Euler, whereas

1 * (Phi_c° - Phi_c°°)/(dt) = - 1 * L(Phi_c°)

is implicit Euler. Both of those equations are only first order approximations and have numerical error that affects the two methods in different ways (e.g. in terms of stability). I'm wondering if I'm just thinking about this in the wrong way, but I can't why one can say that the above equation is valid and would still retain the second order accuracy of the CN method?

[andyru]

As I understand, the following is true:
CN= IE+EE
Hence the treatment of right side is allowed.
You can check that easily with paper and pencil.

And yes, for me, CN is not implicit.
It is an semi-implicit approach in my opinion,.

[rajibroy]

Please see the derivation and assumption that goes into OpenFOAM CrankNicolson implementation at a separate forum post


Crank Nicolson scheme implemented wrong?