[karman]

Hi all,

I need to add an electric body force term (let's call it Fci) in the Navier-Stokes momentum equation and I chose the simpleFoam solver because I have an incompressible and steady-state flow...so I first modified the createFields.H file in this solver by adding the new entries highlighted in red:

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
);

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

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh
);

#include "createPhi.H"


label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());


singlePhaseTransportModel laminarTransport(U, phi);

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

#include "createMRF.H"
 #include "createFvOptions.H"

(Guess I need the fluid density rho so that dimensions match...)


Afterwards I added the bodyforce term Fci divided by rho (also highlighted in red) in the UEqn.H file:

    // Momentum predictor

    MRF.correctBoundaryVelocity(U);

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p)+Fci/rho);

        fvOptions.correct(U);
     }

But now I'm a little bit confused because (in at least one former post here in this forum) I've also seen a different approach of implementing a source term into the momentum equation, accordingly:

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
      + Fci/rho
       );

Anyone knows which of these two approaches is the right way to add a source term??
Guess underrelaxation and the command fvOptions.constrain(UEqn) have the effect that I get different results...


Many thanks in advance and stay healthy in these difficult times!

[HPE]

I would create a new fvOptions instead of modifying simpleFoam.

[Mahmoud Abbaszadeh]

Originally Posted by HPE

I would create a new fvOptions instead of modifying simpleFoam.

Hi,

I am simulating an open channel flow with simpleFoam.
I used fvOptions to add gravity as a momentum source. What I'd like to know is how to set the outlet BC? I can set it as zero fixedValue as we have hydrostatic Pressure there. Do you have any idea how to set HP pressure at the outlet?

Cheers