CFD Online Logo CFD Online URL
Home > Forums > OpenFOAM Running, Solving & CFD

Implementtaion of Crank Nicolson scheme in OpenFOAM

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Display Modes
Old   December 9, 2015, 11:15
Default Implementtaion of Crank Nicolson scheme in OpenFOAM
Hao Chen
Join Date: Aug 2014
Posts: 66
Rep Power: 0
hchen is on a distinguished road

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>


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

    if (evaluate (ddt0))

    fvm.source() = 
        + 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?
hchen is offline   Reply With Quote

Old   January 28, 2016, 23:57
New Member
Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 6
stevemoore1981 is on a distinguished road

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):171203

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. 5067, 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!
stevemoore1981 is offline   Reply With Quote

Old   February 12, 2016, 08:51
Exclamation solution to off-centre implementation...
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 9
andyru is on a distinguished road

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

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

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)
+ offCentre_(ddt0.internalField())

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.

andyru is offline   Reply With Quote

Old   February 15, 2016, 05:21
New Member
Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 6
stevemoore1981 is on a distinguished road

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?
stevemoore1981 is offline   Reply With Quote

Old   February 17, 2016, 18:58
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 9
andyru is on a distinguished road
As I understand, the following is true:
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,.
andyru is offline   Reply With Quote


Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 06:36
2D Mesh Generation Tutorial for GMSH aeroslacker Open Source Meshers: Gmsh, Netgen, CGNS, ... 12 January 19, 2012 04:52
Cross-compiling OpenFOAM 1.7.0 on Linux for Windows 32 and 64bits with Mingw-w64 wyldckat OpenFOAM Announcements from Other Sources 3 September 8, 2010 06:25
Adventure of fisrst openfoam installation on Ubuntu 710 jussi OpenFOAM Installation 0 April 24, 2008 14:25
Crank Nicolson. ! Main CFD Forum 0 September 5, 2005 12:52

All times are GMT -4. The time now is 00:59.