[chriss85]

I am wondering how additional source terms should be added to the velocity equation in sonicFoam.

Should it be done in

1. The velocity equation
2. The flux correction
3. Both

In the flux correction a term is subtracted from U that is proportional to the pressure gradient. I'm wondering if I need to do this with other forces as well.
I guess I need to, because otherwise the forces won't ever be considered unless I enable the momentumPredictor?
Looking at http://openfoamwiki.net/index.php/IcoFoam seems to indicate that the pressure force term which is added in the PISO loop is very important, but it doesn't talk about other forces which are independent from the pressure.

I'm wondering if this special treatment of the pressure term is due to the p-U coupling, or due to the matrix magic being done in the PIMPLE implementation.

Can anyone shed some light on this?

[chriss85]

I made some tests:

1. Additional source terms only in momentum predictor & momentumPredictor off
2. Additional source terms only in momentum predictor & momentumPredictor on
3. Additional source terms in momentum predictor and flux corrector & momentumPredictor off
4. Additional source terms in momentum predictor and flux corrector & momentumPredictor on

Looking at the results, it appears that the terms are needed in the flux corrector as well. (1 & 2) and (3 & 4) are very similar to each other.
However, the terms only lead to slightly different global pressure, while I would expect an unsteady pressure distribution in a steady case, due to other forces (namely, Lorentz force) acting against the pressure. Is the pressure which is being calculated here the real pressure or is it somewhat different and can't be completely related to the real pressure?

Edit: for completeness, here are the lines of code:

//flux predictor:
if (pimple.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p) + j^B);
    K = 0.5*magSqr(U);
}
//...
//flux corrector:
volScalarField rAU(1.0/UEqn.A());
U = HbyA - rAU*fvc::grad(p) + rAU*(j^B);

[chriss85]

Can anyone shed some light on this topic?
I have tried a few different approaches with different results, but I'm not sure which one is correct.

[chriss85]

Looking at the implementation of the fvOptions framework leads me to believe that the source terms should be added to the RHS of the UEqn before the momentum predictor step. I think that the pressure calculation should take this into account properly then because
a) fvOptions should allow to add source forces
b) Pressure calculation uses the UEqn matrix.

However, I'm still wondering the same thing for the rhoCentralFoam solver.

[danny123]

Hello,

If you want to add a source term to the momentum equation, you may add it to Ueq, but it is mandatory to add it to pEq.

pEq uses the matrix coefficient of Ueq (derived only from left side of Ueq) to iterate a pressure field consistent to both momentum and conservation. The new velocity field is calculated thereafter as is the volumetric flux phi.

The momentum predictor calculates a new velocity field, which respects the momentum equation, but since the pressure is a guess field of a previous iteration, it does not respect mass conservation necessarily. This is why it is called a "predictor".

Regards,

Daniel

[chriss85]

Are you sure it is only using the LHS of UEqn?
Because the fvOptions term is specified on the RHS of UEqn and would not be considered then?

[danny123]

This is my current understanding. Look at http://infofich.unl.edu.ar/upload/3b...7523c8ea52.pdf

This source may not be perfect, but look e.g. into the interfoam code for UEqn:

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U)
    );

At this point, the source term of the matrix eq. if any would be

fvOptions(rho, U)

,no?

I am not quite sure what this is, but for sure the momentum equation is not complete.

This is added later, only if there is a momentum predictor:

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );

Note that UEqn is not updated, even if the momentum predictor is called! So, the source part (RHS within the predictor) in respect to U is not considered in the pressure equation. The pressure equation uses UEqn.A (which anyway does not contain a source part) and UEqn.H (which does, but it is

fvOptions(rho, U)

, not the real source stuff.

By the way, the RHS needs to be considered in pEqn, so it is indeed to be added, but this is done seperatly (as phig).

I am not sure what the fvOptions thing does. But maybe you can find out and post it.

Regards,

Daniel

[chriss85]

The fvOptions framework is meant to allow adding source terms to equations, e.g. mass source, impulse source (force) or energy source (power). The types of sources can be selected by the user, so the type of terms is not known yet at compile time. Examples include things like coriolis force, gravity or a heat source.

Here's a quote from the (interesting) document you linked to regarding the == operator:

This operator is here entirely for stylistic reasons,
since the code automatically rearranges the equation (all implicit terms go into the matrix, and all
explicit terms contribute to the source vector).

This means that the fvOptions terms are included in the UEqn. There are also sources for the pEqn, so maybe I still need to add something there...It would be nice to have an example for adding an explicit impulse source term. I will probably have to do some more research regarding the fvOptions implementation. Here's a good starting point for fvOptions btw: http://www.sourceflux.de/blog/softwa...ign-fvoptions/

Edit: I looked in the fvOptions code and I believe that impulse sources only work on the UEqn but not on the pEqn. I'm not really sure about this as the code is nested through a few classes and I don't understand it fully right now.

[danny123]

This is very helpful, what you found. It basecally means that fvOptions is a file in the systemfolder where the user can add source terms to UEqn without the need of recompiling the solver. In my case, there is no such file. There are a couple of examples in the tutorial, where they are present.

You are correct that the == is meant to generate an fvmatrix object (which is UEqn). The source parts are typically RHS. I do not know if that is mandatory or not. I had plotted the fvmatrix object out and the structure is somewhat like this:

- the matrix coefficients
- psi (the initial guess for the equation A x psi = S)
- the source term
- the patches

The fvmatrix object needs to know what psi is. So, the rule, I think, is that all terms that are psi related (containing U in our case since U = psi) have to be LHS. So, at least the p terms have to be RHS.

Now, since the remaining part of the source is not built into UEqn, it would be interesting to know what happens if we do build it into. I think it boils down to the change that we might see in UEqn.H (since UEqn.A does not contain source parts).

Another point: why does fvOptions does not appear in pEqn? It should show up in phig. If it does not it is kind of either or (in UEqn or phig), but how do you judge where to put it? Does it make a difference?

Regards,

Daniel

[chriss85]

What I meant to say is that it doesn't matter if you put your terms on the LHS or RHS of the equation, because OF will decide where to put it in the discretization itself. Implicit terms (depending on the field being solved for) go on the LHS while explicit terms (sources mostly, but they are allowed to depend on the field from a previous iteration as well) go on the RHS.
I don't think it's possible to treat source terms which don't depend on the field implicitly.

fvOptions does appear in pEqn in sonicFoam, but I think the implementations that add impulse sources don't make use of it, they only modify UEqn as far as I can see.

[danny123]

I think you are correct on this. The difference wheter it is treated as psi is the way how the syntax is formulated e.g. fvm and the like. You have to put the arguments in the right order though.

However, only that stuff, which is defined as UEqn is used in pEqn later on. Within the example I have given, the solve function acts on UEqn + some other source term + the p term in the momentum predictor. I wonder if the "other source term" should not be part of the UEqn formulation. You then will get back to the standard having only the

==-fvc::grad(p)

to be added within solve.

Regards,

Daniel

[Thaw Tar]

Thank you very much you all!!

This discussion is really a useful one.

But I do not understand well yet. Do I need to write a new solver if
I want to add body force (inertia force) to pimpleFoam and interFoam solvers?

I need to simulate the flow around a bunch of circular cylinders (with all vortex shedding and separation phenomena neglected) by using body force. Do I need to add body force for each cylinder or just a body force field for all cylinders?

Regards,
Thaw Tar

[chriss85]

You might be able to use the fvOptions framework to add a force to the solver. You'll have to check if you can find a proper module for it though. Otherwise you can add the force term to the equations like we discussed it here and recompile the solver (wmake).