CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Source terms in sonicFoam / PIMPLE

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   November 10, 2014, 07:53
Default Source terms in sonicFoam / PIMPLE
  #1
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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 is offline   Reply With Quote

Old   November 11, 2014, 04:24
Default
  #2
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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:
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);

Last edited by chriss85; November 11, 2014 at 11:45.
chriss85 is offline   Reply With Quote

Old   June 18, 2015, 08:26
Default
  #3
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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 is offline   Reply With Quote

Old   June 22, 2015, 02:56
Default
  #4
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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.

Last edited by chriss85; June 22, 2015 at 10:48.
chriss85 is offline   Reply With Quote

Old   June 24, 2015, 09:47
Default
  #5
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5
danny123 is on a distinguished road
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
danny123 is offline   Reply With Quote

Old   June 25, 2015, 04:56
Default
  #6
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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?
chriss85 is offline   Reply With Quote

Old   June 25, 2015, 06:04
Default
  #7
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5
danny123 is on a distinguished road
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:

Code:
    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

Code:
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:

Code:
    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
Code:
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
danny123 is offline   Reply With Quote

Old   June 25, 2015, 07:08
Default
  #8
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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:
Quote:
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.
chriss85 is offline   Reply With Quote

Old   June 25, 2015, 07:56
Default
  #9
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5
danny123 is on a distinguished road
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
danny123 is offline   Reply With Quote

Old   June 25, 2015, 11:35
Default
  #10
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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.
chriss85 is offline   Reply With Quote

Old   June 25, 2015, 12:20
Default
  #11
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5
danny123 is on a distinguished road
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
Code:
==-fvc::grad(p)
to be added within solve.

Regards,

Daniel
danny123 is offline   Reply With Quote

Old   July 5, 2015, 23:18
Default
  #12
New Member
 
Thaw Tar's Avatar
 
Thaw Tar
Join Date: Apr 2013
Location: Yangon, Myanmar
Posts: 19
Rep Power: 4
Thaw Tar is on a distinguished road
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
Thaw Tar is offline   Reply With Quote

Old   July 6, 2015, 04:59
Default
  #13
Senior Member
 
Join Date: Oct 2013
Posts: 278
Rep Power: 4
chriss85 is on a distinguished road
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).
chriss85 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem compiling a custom Lagrangian library brbbhatti OpenFOAM Programming & Development 2 July 7, 2014 11:32
swak4Foam-groovyBC build problem zxj160 OpenFOAM 18 July 30, 2013 13:14
Source Terms in momentum Balance vidyaraja FLUENT 0 May 25, 2009 13:16
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51
UDFs for Scalar Eqn - Fluid/Solid HT Greg Perkins FLUENT 0 October 11, 2000 03:43


All times are GMT -4. The time now is 06:47.