how to add source-term in momentum equation for interFoam?
 

hi everybody

I'm trying to modify the interFoam equations by adding a source term to the momentum equation, i.e. UEqn.H.
Whether it's better to change UEqn.H or pEqn.H (i.e. left or right side of equation (?)), I dont know, but since my source term is dependent on U, I assumed it would make sense to add to UEqn.H.

the source-term that I want to add is s.th. like this (in the code I had the actual function +.. * ... - etc) (see attachment)

f(alpha1) * g(rho1, rho2) * h(g, ddt(U))

my understanding is that the solver calculates the variables and their values at the cell faces rather than at the cell-centers?! If that's true, how do I incorporate that? (I tried ..::interpolate(...))

Every help, idea, hint... will be greatly apprecitated!

***EDIT:
- so far I couldnt find an answer in this forum (and others), but maybe there is already a similar thread? a link to such a thread would be equally good for an answer/help ;)

maybe someone has some experience with other solver/code-modifications?

If you are adding a source term to the momentum equation, you will need to modify the momentum equation, which is the equation being solved in UEqn.H. If your source term depends on U, you will need to think carefully about whether to treat it explicitly or implicitly. See Patankar's "Numerical Heat Transfer and Fluid Flow" for implicit/explicit handling of source terms. If your source term is proportional to ddt(U) as you suggest in your post, then it seems like you should be modifying the existing time derivative.

1st: thx for replying...

that's what i thought, thx 4 confirming

it seems to me that in OpenFOAM the implicit/explicit handling is addressed by the fvc/fvm-functions, but what do i know :D. (i would've tryied both anyway to see which one to pick)

i read through patankar, but i didnt find any passages, that would (explicitly) cover source-term handling as implicit/explicit. i dont expect it, but maybe you remember which chapter you were thinkin about?

also, my problem seems to be more of a - programming/implementing one: when i add my source term and try to compile it, i get all sorts of error messages. i hoped someone could give me general advice on how to implement a term.

i thought about that, but that wouldnt do me any good, since i still have to deal with the gravitational vector g that i want to include in my source term,...
...and, it would be nice to be able to add/remove a source term without editing the original equation.

any other ideas/suggestions/advice? help's much appriciated

I found it in the appendix, under "source-term linearization."

Try looking at the distributed solvers. There are 40 different momentum equations there, and probably a hundred more equations with various source terms to see for examples.

What kind of a source term is it? Can you post the equation you want to code in?
If the term should be treated implicitly, I would use fvm, otherwise fvc. They are not interchangeable, as the return types are different.
fvm returns the contribution to the coefficient matrix A(the Ax=b system yiou get after discretizing your equation) , while fvc returns a vector that you can just add to the right hand side of the equation.
My understanding for implicit and explicit handling is that if by implicit handling the matrix A becomes easier to invert, then do it. Else, use explicit source term.

i have to admit that this is - at least for a novice - quite difficult. finding a solution that way is sort of very time intensive and not guaranteed as i have trouble with implementing the code in the first place. but i'm on to it and if i find something in the existing equations that resembles my problem, i'll post where to find it, so that others find it as well.

But maybe i should start giving more details so that my problem becomes clearer. ;)

mathematically speaking: S = (1-alpha1)*(rho2-rho1)*(g - d(U)/dt)

i figuered that my problems, i.e. the error messages upon compiling, arise as soon as a put the 'g' in there. the rest is implemented as is, except fvm::ddt(U) ^^

another question on the last post:
maybe that's a stupid question, i apologize for my non-knowledge, but how do i know if by implementing a term implicitely, the matrix will become easier to invert? (And also, why is this important?) I should explain that i started with openfoam from a "use as tool" kind of standpoint, i.e. right now i'm more interested in the cases that i'm running than in how OF achieves this...

I added the line you describe to the UEqn in UEqn.H and got no errors for both OpenFOAM 1.6-ext and OpenFOAM 2.0.x. Maybe you need to attach your modified UEqn.H, tell us what version you are using, and most importantly, tell us what the error is that the compiler gives?

you are right ;)

here's why it didnt work for me before: as i said i thought the equations where/are solved at the cell-faces... but the variables, being volScalarFields/-VectorFields, represent -as far as i know, please correct me if i'm wrong - the whole cell value, stored at the center of the cell. therefore i thought it would be necessary to project/interpolate the values to the surfaces. i tried this...:

surfaceScalarField rhoF = fvc::interpolate(rho)

i saw this in other parts of the code and thought it would be reasonable to implement it that way. but again, what do i know, i'm just at the first step of a high ladder when it comes to openfoam knowlegde. i would be curious though if my attempt was that bad at all and where the mistake in my "reasoning" exactly was. if somebody has the answer for this, please feel free to share ;)

anyway, when using my "created variables" (rhoF etc.) in the equation, compilation failed.

You're correct (now) that the dependent variables are usually stored as cell-centered values (volScalarFields and volVectorFields). But there are also many situations where the face-centered values (e.g., surfaceScalarFields) are required. I'd suggest finding Hrv Jasak's PhD dissertation (available on-line, check Google) for the explanation of how the equations are discretized to understand much of that reasoning.

In the future, I'd also suggest posting more detailed explanations of your problems. Saying that 'it failed' is not nearly as helpful as posting the exact error message and the portion of the code that the error refers to. But from what you've said, it sounds like the compiler might have been telling you that you couldn't mix and match your surfaceScalarFields and volScalarFields. For example, you wouldn't be allowed to add one to the other.