CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Time-marching scheme in dbnsFoam

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

LinkBack Thread Tools Search this Thread Display Modes
Old   August 20, 2018, 10:15
Default Time-marching scheme in dbnsFoam
New Member
Sam Maloney
Join Date: Aug 2018
Location: Zürich, Switzerland
Posts: 2
Rep Power: 0
Aitrus is on a distinguished road
I am trying to understand the time-marching scheme that is implemented in the dbnsFoam (and dbnsTurbFoam) solvers in foam-extend (I am using version 4.0).

The solver source file specifies it only as a "explicit time-marching" scheme and in my googling around, I found on page 11 of this presentation by Jemcov et al. [1] a reference to this paper by Arnone et al. [2] which says that the scheme should look like the following:

Q^{(i+1)}=Q^{(0)}+\beta_i \Delta t R(Q^{(i)})

where Q^{(0)}=Q^n and Q^{(n+1)}=Q^{(i_{max})} and i specifies the stage of the Runge-kutta scheme. One difference, however, is that the coefficients used in the code appear to be \beta_i=[0.11, 0.2766, 0.5, 1], while the paper [2] appears to specify \beta_i=[1/4, 1/3, 1/2, 1]. If anyone knows a reference for why these coefficients were chosen, that would be great?

However, my bigger question comes when looking at the code itself in dbnsFoam.C, it looks like this:

        forAll (beta, i)
            // Solve the approximate Riemann problem for this time step

            // Time integration
              + fvc::div(dbnsFlux.rhoFlux())

              + fvc::div(dbnsFlux.rhoUFlux())

              + fvc::div(dbnsFlux.rhoEFlux())

#           include "updateFields.H"
My understanding of the "solve" operation is that it will solve the discrete equation created using an Euler discretization of the time-derivative, which in general looks like:

\frac{Q^{(n+1)}-Q^n}{\Delta t}+R(Q^n)=0

so in the above code it seems to me that it would produce update equations at each iteration of the for loop that look like:

Q^{(i+1)}=Q^{(i)}+\beta_i \Delta t R(Q^{(i)})

which differs from the original form in [1] and [2] in that the first RHS term is Q^{(i)} instead of Q^{(0)}.

My question then is whether I am missing something in my understanding of how the solve operation works? It seems like on the first iteration of the for loop, the convervative variables rho, rhoU, and rhoE should be updated to the (n+0.11) \Delta t timestep and then the code in updateFields.H will update the remaining primitive variables. Thus on the next for loop iteration, the original Q^{(0)} was not saved anywhere, so it will have to use Q^{(1)}=Q^{(n+0.11)} instead.

If anyone can shed light on what I might be missing here it would be most appreciated!

[1] Density Based Navier Stokes Solver for Transonic Flows,

[2] Multigrid time-accurate integration of Navier-Stokes equations,
Aitrus is offline   Reply With Quote

Old   August 26, 2018, 04:25
New Member
Sam Maloney
Join Date: Aug 2018
Location: Zürich, Switzerland
Posts: 2
Rep Power: 0
Aitrus is on a distinguished road
Alright, so after more digging through code and documentation, I have learned that the "GeometricField" class has a member "oldTime()" (and some other associated methods) which allow it to automatically store old values in the background during the solve operation. The Euler ddt scheme uses this oldTime() value to build its discretization, and since the time index is only incremented after all four Runge-Kutta stages have completed, this means it is indeed Q^{(0)} and not Q^{(i)} which appears on the RHS of the update equation, giving the desired form shown in the references.

Therefore my only remaining question is the motivation for the choice of the \beta_i=[0.11, 0.2766, 0.5, 1] coefficients? If anyone knows a reference for these values, or has any other information on their properties (order of accuracy, is it SSP?, ...) that would be very helpful!
Aitrus is offline   Reply With Quote

Old   November 15, 2019, 14:28
New Member
Join Date: Mar 2016
Posts: 3
Rep Power: 10
TommyGun is on a distinguished road
Hi, I've found only one paper with these coefficients so far. Try this, page 14:

It seems like the first two coefficients 0.11 and 0.2766 are result of trial and error around the classical coefficients 1/4 and 1/3 which cause bigger stability region of four-stage Runge-Kutta scheme.

Hope this helps.
TommyGun is offline   Reply With Quote


dbnsfoam, extend 4.0, foam-extend, foam-extend 4.0, runge-kutta

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
[Other] Contribution a new utility: refine wall layer mesh based on yPlus field lakeat OpenFOAM Community Contributions 58 December 23, 2021 02:36
pimpleDyMFoam computation randomly stops babapeti OpenFOAM Running, Solving & CFD 5 January 24, 2018 05:28
High Courant Number @ icoFoam Artex85 OpenFOAM Running, Solving & CFD 11 February 16, 2017 13:40
[solidMechanics] solidMechanics gear contact in rotation nlc OpenFOAM CC Toolkits for Fluid-Structure Interaction 3 January 11, 2015 06:41
Star cd es-ice solver error ernarasimman STAR-CD 2 September 12, 2014 00:01

All times are GMT -4. The time now is 18:50.