[syavash]

Dear All,

How should I add an explicit source term to momentum equation in the incompressible explicit 4th-order runge-kutta method? I have already implemented the method but I can't be sure about the details of adding a source term. At the moment, I have added the explicit source term at the end of the 4th step by adding the source to the velocity quantity, right before solving the Poisson equation for pressure. However, I wonder if it is the correct place.
Please help me to find the answer, possibly by introducing an appropriate reference addressing the issue.

Kind Regards,
Syavash

[FMDenaro]

Originally Posted by syavash

Dear All,

How should I add an explicit source term to momentum equation in the incompressible explicit 4th-order runge-kutta method? I have already implemented the method but I can't be sure about the details of adding a source term. At the moment, I have added the explicit source term at the end of the 4th step by adding the source to the velocity quantity, right before solving the Poisson equation for pressure. However, I wonder if it is the correct place.
Please help me to find the answer, possibly by introducing an appropriate reference addressing the issue.

Kind Regards,
Syavash

That depends on the linearity or non-linearity of the source term, if it is stiff, etc

[syavash]

Originally Posted by FMDenaro

That depends on the linearity or non-linearity of the source term, if it is stiff, etc

Thanks Prof. Denaro for your quick response. It is a simple constant source which is supposed to drive the flow inside a planar channel for LES.

The implementation is like the following:

//4th step
        U = U + deltaP * dt;          //deltaP is the source vector
        phi = (fvc::interpolate(U) & mesh.Sf());
        dU = dt*(fvc::laplacian(turbulence->nuEff(), U) - fvc::div(phi,U)); //nuEff is the summation of nuLam and nuSgs
        Uc = Uc +b4*dU; U = Uc;
        #include "pressureCorrection.H"    //pressure correction

The file pressureCorrection.H contains the following lines:
Please note that pEqn is the laplacian of pressure.

U.correctBoundaryConditions();
solve(pEqn == fvc::div(U)/runTime.deltaT());
U = U - (fvc::grad(p) * runTime.deltaT());
U.correctBoundaryConditions();

Is this the correct implementation?

Kind Regards,
Syavash

[FMDenaro]

if this is the constante pressure gradient that drives the main flow you just add at the end the term dt*grad p, this an exact integration.

[syavash]

Originally Posted by FMDenaro

if this is the constante pressure gradient that drives the main flow you just add at the end the term dt*grad p, this an exact integration.

Thanks. Let me ask if I got it right. You mean something like the following in the 4th step:

4th step
        #include "pressureCorrection.H"    //pressure correction

         U = U + deltaP * dt;                 //Ending line of the code

In other words, should it not be included in the pressure correction equation?

Kind Regards,
Syavash

[FMDenaro]

Originally Posted by syavash

Thanks. Let me ask if I got it right. You mean something like the following in the 4th step:

4th step
        #include "pressureCorrection.H"    //pressure correction

         U = U + deltaP * dt;                 //Ending line of the code

In other words, should it not be included in the pressure correction equation?

Kind Regards,
Syavash

yes, you can also add it out of the RK4 steps