
[Sponsors] 
November 10, 2014, 07:53 
Source terms in sonicFoam / PIMPLE

#1 
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
I am wondering how additional source terms should be added to the velocity equation in sonicFoam.
Should it be done in
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 pU coupling, or due to the matrix magic being done in the PIMPLE implementation. Can anyone shed some light on this? 

November 11, 2014, 04:24 

#2 
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
I made some tests:
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. 

June 18, 2015, 08:26 

#3 
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
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. 

June 22, 2015, 02:56 

#4 
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
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. 

June 24, 2015, 09:47 

#5 
Senior Member
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5 
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 

June 25, 2015, 04:56 

#6 
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
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? 

June 25, 2015, 06:04 

#7 
Senior Member
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5 
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) ); Code:
fvOptions(rho, U) 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() ) ); Code:
fvOptions(rho, U) 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 

June 25, 2015, 07:08 

#8  
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
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:
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. 

June 25, 2015, 07:56 

#9 
Senior Member
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5 
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 

June 25, 2015, 11:35 

#10 
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
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. 

June 25, 2015, 12:20 

#11 
Senior Member
Daniel Witte
Join Date: Nov 2011
Posts: 127
Rep Power: 5 
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) Regards, Daniel 

July 5, 2015, 23:18 

#12 
New Member
Thaw Tar
Join Date: Apr 2013
Location: Yangon, Myanmar
Posts: 19
Rep Power: 4 
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 

July 6, 2015, 04:59 

#13 
Senior Member
Join Date: Oct 2013
Posts: 278
Rep Power: 4 
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).


Thread Tools  
Display Modes  


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 
swak4FoamgroovyBC 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 