
[Sponsors] 
December 9, 2015, 11:15 
Implementtaion of Crank Nicolson scheme in OpenFOAM

#1 
Member
Hao Chen
Join Date: Aug 2014
Posts: 66
Rep Power: 0 
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: Code:
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; } 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 nonorthoganity correction?) 3. why the code use ddt0 for CN scheme? 

January 28, 2016, 23:57 

#2 
New Member
Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 7 
Hi
I've recently come across this same problem, namely trying to understand how CrankNicolson 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^n1 + 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 freesurface 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 HeatConduction 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! 

February 12, 2016, 08:51 
solution to offcentre implementation...

#3 
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 10 
hello,
I faced the same problem. I wanted to know, how the implementation of CN with offcentre=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 offcentre 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 tt° =implict with t°  t°° You keep with this trick the generic matrix addition approach for fvm calculus. 

February 15, 2016, 05:21 

#4 
New Member
Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 7 
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? 

February 17, 2016, 18:58 

#5 
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 10 
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 semiimplicit approach in my opinion,. 

Thread Tools  
Display Modes  


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 
Crosscompiling OpenFOAM 1.7.0 on Linux for Windows 32 and 64bits with Mingww64  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 