[Deep111090]

Hi All,
I asked this question in another forum also, but I haven't reached any conclusion yet. I apologize if you have to read this again.
I am solving a 2D lid-driven flow using icoFOAM algorithm, one of the tutorial cases.


I decided to experiment with it, the First step of which involved the calculation of the velocity and pressure up to convergence.


In the second step I modified the code of the icoFoam solver by commenting out the pEqn.solve() line (i.e. not solving the pressure) and used the pressure obtained in the first step to calculating the velocity profiles. In this step, the velocity was initialed from zeros with boundary conditions similar to the first step.


My assertion was that the velocity profile obtained in the second step should be exactly similar to that of the first step since they correspond to the same pressure field.


But the results are very much Reynolds number dependent. At low Reynolds number, I do have similar velocity profiles as of both the cases (though not exact match ~differ by order of 1e-4). But at Reynolds's number greater than 200 the velocity profiles are not remotely similar, i.e. the velocity fields obtained for the second step are incorrect and exhibit discontinuities in their profiles.


Can you please guide me in the right direction as why this is happening.

Thanks in advance.
PS: if my question is not clear please feel free to ask for more details.

[geth03]

if you want to know what your steady-state velocity profile looks like when pressure is pre-determined, you do not need to engage with pressure-velocity coupled solvers.

if you know your pressure field, just solve

fvVectorMatrix UEqn
(
fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == -fvc::grad(p));

as many times as you need to reach steady-state.

[Deep111090]

Thanks Geth for the response. I tried that a few times. The problem here is that the "phi" is not conserving the mass locally. Right now the question before me is more fundamental that is it even possible to get a solution by this method.

Essentially this problem summarises to two eqautions (momentum in x and y) and two unknowns(u_x,u_y).

But the given the equations are non linear is it mathematically possible to reach the solution (I may be wrong while stating this line)

[geth03]

mathematically, it is true that the NS-equation is non-linear.

but openfoam does not consider that, it linearizes the convection term in the NS-equation.
fvm::div(phi, U) : look at that expression. here you can see that velocity transports itself.

this way you can actually use an iterative solver, and not need to engage with non-linear systems.

[Deep111090]

yes..now i have question that how shall i calculate phi ,..i am doing it by fvc::interpolate(U)&mesh.Sf().
But by default the openFoam solves it by phi = phiHby-pEqn.flux();

Question 1: is there any difference (as both approaches give slighty different answers for phi).
Question 2: Will this phi be locally mass conserving. as it is not if i don't solve pressure eqaution and just feed pre-determined pressure?

[geth03]

before entering the piso loop, you create phi with fvc::interpolate(U)&mesh.Sf().

in the pressure-velocity coupling part, phiHbyA is the volume flux derived from the UEqn, which does not take into account the pressure part, so it is overguessed.

fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

the pEqn is derived when inserting the momentum eqn in the mass conservation eqn, basically you are looking for a pressure field that conserves your mass flow. also called pressure poisson equation.

phi = phiHbyA - pEqn.flux();

here you correct phiHbyA, bc it was overguessed previously, to solve for a volume flux that considers a pressure field that conserves mass flow.

i hope you understand by now what each line is doing.

so if you do not solve for pEqn, you will always substract the same amount from your previous volume flux, which will change its value, of course.


to sum up, you solve for U with known pressure field, calculate then the volume flux without it and then correct the volume flux again with the known pressure field.

[Deep111090]

Hi I have tried the aforementioned method...my code looks like



fvVectorMatrix UEqn
(
fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == -fvc::grad(p));


phi = fvc::interpolate(U)&mesh.Sf();


I got following results for Reynolds number 250. The phi is not conserving mass in this case. Hence the unrealistic velocity profiles. The results for using pre-determined pressure are attached as images

[Deep111090]

The actual results should look like as attached. I am not sure what is causing this phi to be not conserving mass locally. Even though we are solving from pressure which is predetermined and known to be correct.

[geth03]

how many times did you execute the code?
only once?
did you try to solve it many more times? in a loop for example?

[Deep111090]

Hi.. i looped it many times. the momentum euations had resuidals in range of 1e-6 and the continuty errs were high for local.

[geth03]

i think the main problem here is phi. the volume flux which is not known but is treated as known and then re-calculated again for a known pressure field.

i think as U_wall or U_top becomes larger and therefor Re becomes larger, the difference or discrepancy of this approach becomes more obvious.

the main problem of a non-linear set of equations is that for large systems it is much slower than the iterative approach, given an unknown pressure field. the way it is coded in OF both U and p are determined with the iterative approach or segregated, within a loop, to ensure mass conservation.

your problem is, that you already know your pressure field. so you don't need to the approach of fluid transporting itself, i guess.

i don't know if OF is capable of solving non-linear systems.

[Deep111090]

Thanks for the response.
Right now i am under the impression that it can not be done without solving for continuity..which seems to be tricky as I am unable to figure out the way to enforce continuity on "phi" without solving the Poisson's equation. Can you suggest some way to enforce continuity on each cell locally without solving Poisson's .

[geth03]

unfortunately not.

i think without the pressure-velocity coupling, your solution of the velocity for a given pressure field will be highly dependant of your initial values for velocity if you try to solve for velocity in the manner how it is implemented in openfoam.

if you have time, could you use different initial velocity fields to see how sensitive your end solution is for the initial velocity field?

you could
1) initialize your new solver with the end solution of the pressure and velocity field of icoFoam to see if things change
2) initiliaze with fields from a time you get with icoFoam, lets say endTime/2 or 0.75 x endTime.

[Deep111090]

Thanks Geth. You are right.

Its all my hunch and i would appreciate to be corrected. The problem here is coupling. It is way the openfoam is and use rhie chow interpolating in spirit of it not explicitly. I may be wrong but from what i can see is that openfoam sort of spread the elements of the rhie chow in entrie solution algorithm which causes problem if we comment some sections out.
The reason why it should not work with phi = fvc::interpolate(U)& mesh is because there is no coupling between phi and pressure due to this interpolation.
The reason it is not working with phi = phiHbyA-pEqn.flux() was al ready stated by you.

[Tobermory]

I came here from your thread on Rhie Chow, and can understand now why you were asking about that. To reiterate:

how shall i calculate phi ,..i am doing it by fvc::interpolate(U)&mesh.Sf().
But by default the openFoam solves it by phi = phiHbyA - pEqn.flux();

These 2 approaches generally do not give the same answer; the simple interpolation approach is incorrect, and will not conserve mass. The OpenFOAM solvers updates their face fluxes, phi, based on the momentum interpolation (also called Rhie Chow) approach which is the second method you mentioned above.

Coming back to your original experiment (which I approve, by the way - the best way to learn is to tinker - so well done!), you have only commented out the p.solve() line in the corrector loop ... the code is still calculating the momentumless face flux and then trying to correct this with the pressure gradient - but you have broken the link between the coefficients in the pEqn and the pressure values since you're not solving the pressure equation. This means that there's no pressure velocity coupling, as was observed above, and so I can see why the solution might go off the rails.

For example, in the normal solver, if you perturb the velocities in a few cells, then these will generate a pressure field that will ultimately dampen them down again since there is nothing there to maintain the pressure perturbation. In your solver you can perturb the velocities and yet no pressure perturbation occurs, and so no restoring force occurs. The only stabilising effect is viscosity ... which explains the Re dependence.

[Deep111090]

Thanks for your response Tobermory. It seems that now we can reach some conclusion from ideas given by yourself and Geth. I am not sure which coefficients have have links broken. As the rAU remains the same through out all the iterations. And judging from the discussion on pEqn.flux() in this thread ( Description of flux() method ). The pEq.flux() depends upon (1/Ap) which is rAU and other thing it depends on is p.

[Tobermory]

but you have broken the link between the coefficients in the pEqn and the pressure values since you're not solving the pressure equation

To understand this, start by looking at the corrector loop. The pEqn is set up as follows:

fvScalarMatrix pEqn
   (
      fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
   );

This builds the fvMatrix for the pressure equation, populating the coefficients using phi, H and A - ie these depend on the velocity and flux fields, as well as the geometric constants (if mesh is stationary). The pEqn().solve would then solve the pressure field to give a p field that is "in equilibrium" with the matrix coefficients. You have commented out the solve, so the p-field and the pEqn matrix coefficients are no longer in balance - that is the broken link that I am referring to.

Later on in the corrector loop, the face flux is corrected with

   phi = phiHbyA - pEqn.flux();

where pEqn.flux() pulls the data directly from the pEqn fvMatrix (to avoid continuity errors) ... the problem is that the p values are not correct - they are not in balance with the face fluxes, and so will not preserve continuity.

Have a think of the analogy that I mentioned in an earlier post: if you were to perturb the velocities in a few cells, then you would expect to see some impact on the pressure field, which itself would generate some kind of feedback correction to the velocity field, smoothing out the perturbation. In your setup, you have disabled that feedback, so the velocity is free to wander.

[Deep111090]

Thanks I now have a cleared idea abt this.

[geth03]

hi again

i once more looked at this problem and i think i could write the code for that problem. the solution strategy is the same i described in post #2 with the difference that the updated flux is updated with a slight difference.

if you know the pressure somehow (like from measurement or from another simulation) the main problem you need to deal with is the non-linearity of the divergence of velocity bc the flux is the product of velocity and the velocity is carried by the flux!

i did run the icoFoam-cavity tutorial and got the pressure field p as result. afterwards i inserted p in my solver and solved the UEqn and updated phi again and again until no change.

Here is the code:
pKnownFoam.C

#include "fvCFD.H"
#include "simpleControl.H"

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

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

	simpleControl simple(mesh);

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
	
	while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.userTimeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {
            fvVectorMatrix UEqn
			(
				fvm::div(phi, U)
			  - fvm::laplacian(nu, U)
			  ==
			  -fvc::grad(p)
			);

            UEqn.solve();
			
			phi = 0.5*phi + 0.5*fvc::flux(U);
        }
		
        runTime.write();
    }

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

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

    return 0;
}

createFields.H

Info<< "Reading physicalProperties\n" << endl;

IOdictionary physicalProperties
(
    IOobject
    (
        "physicalProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar nu
(
    "nu",
    dimViscosity,
    physicalProperties.lookup("nu")
);

Info<< "Reading field p\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


#include "createPhi.H"


label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solution().dict().subDict("PISO"), pRefCell, pRefValue);
mesh.schemes().setFluxRequired(p.name());

Here is the result of both simulations. in the icoFoam simulation velocity is calculated with pressure-velocity coupling, the pKnownFoam simulation has the pressure field as known input from the icoFoam simulation and velocity is calculated according to given code.