[haghgoo_reza]

Hello All

I'm trying to implement an explicit Runge-Kutta 4 method to an incompressible fluid inside a 2-D cavity. I have already modified icoFoam.C but it does not work.

could anybody advise me how to do it or refer me to a suitable tutorial?

my modified icoFoam gives out the following error:
[Make/linux64GccDPOpt/myicoFoam.o] Error 1


modified icoFoam.C


int main(int argc, char *argv[])
{
#include "setRootCase.H"

#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

volVectorField H("H", U);
H = -fvc::div(phi, U)+fvc::laplacian(nu, U);
surfaceScalarField phitilta
(
"phitilta",
(fvc::interpolate(H) & mesh.Sf())
);

adjustPhi(phitilta, U, p);


fvScalarMatrix pEqn
(
fvc::laplacian(p) == fvc::div(phitilta)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();



#include "continuityErrs.H"

volVectorField k1("k1", U);
k1 = runTime.deltaT()*[-fvc::div(phi, U)-fvc::grad(p)+fvc::laplacian(nu, U)];

volVectorField k2("k2", U);
k2 = runTime.deltaT()*[-fvc::div(phi, U+0.5*k1)-fvc::grad(p)+fvc::laplacian(nu, U+0.5*k1)];

volVectorField k3("k3", U);
k3 = runTime.deltaT()*[-fvc::div(phi, U+0.5*k2)-fvc::grad(p)+fvc::laplacian(nu, U+0.5*k2)];

volVectorField k4("k4", U);
k4 = runTime.deltaT()*[-fvc::div(phi, U+k3)-fvc::grad(p)+fvc::laplacian(nu, U+k3)];

U = U+(k1+2.0*k2+2.0*k3+k4)/6.0;
U.correctBoundaryConditions();

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
}

[thomasArk47]

Hi,

many misunderstanding I think.

First, your equation on the laplacian(p) cannot work since you use fvc namespace which is the namespace for explicit discretisation. For building a matrix as in a implicit scheme, you need to use fvm namespace.

Second, and much more drastic, your algorithm for solving NS cannot work.
Recall that incompressible NS system is constraint by div(u)=0 and this constraint is reflected at numerical level by the fact that one cannot use just a "consistant" algorithm to solve equation. Well defined predictor/corrector is the simplest scheme that can work. Roughly speaking, 3 phases:
a- solve equation on U (with/without p)
b- solve laplacian on p
c- correct U
The order and the phase 3 is fundamental. Note also that since you are working with collocated shceme (velocity and pressure are defined both at the cell centers), you need a well defined stabilisation on the velocity flux field (so-called Rhie-Chow stabilisation).

I think you should review your algorithm before of all.

Good luck.