|
[Sponsors] |
Implementtaion of Crank Nicolson scheme in OpenFOAM |
![]() |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
![]() |
![]() |
#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 non-orthoganity correction?) 3. why the code use ddt0 for CN scheme? |
|
![]() |
![]() |
![]() |
![]() |
#2 |
New Member
Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 13 ![]() |
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! |
|
![]() |
![]() |
![]() |
![]() |
#3 |
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 17 ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#4 |
New Member
Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 13 ![]() |
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? |
|
![]() |
![]() |
![]() |
![]() |
#5 |
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 17 ![]() |
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,. |
|
![]() |
![]() |
![]() |
![]() |
#6 |
New Member
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 12 ![]() |
Please see the derivation and assumption that goes into OpenFOAM CrankNicolson implementation at a separate forum post
Crank Nicolson scheme implemented wrong? |
|
![]() |
![]() |
![]() |
|
|
![]() |
||||
Thread | Thread Starter | Forum | Replies | Last Post |
Superlinear speedup in OpenFOAM 13 | msrinath80 | OpenFOAM Running, Solving & CFD | 18 | March 3, 2015 05:36 |
[Gmsh] 2D Mesh Generation Tutorial for GMSH | aeroslacker | OpenFOAM Meshing & Mesh Conversion | 12 | January 19, 2012 03: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 |