CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   how to add source-term in momentum equation for interFoam? (http://www.cfd-online.com/Forums/openfoam/93304-how-add-source-term-momentum-equation-interfoam.html)

 anon_g October 11, 2011 09:10

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

1 Attachment(s)
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 ;)

 anon_g October 12, 2011 09:09

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

 boger October 13, 2011 09:37

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.

 anon_g October 14, 2011 09:43

Quote:
 Originally Posted by boger (Post 327818) 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.
that's what i thought, thx 4 confirming

Quote:
 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.
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.

Quote:
 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.
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

 boger October 14, 2011 14:59

Quote:
 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?
I found it in the appendix, under "source-term linearization."

Quote:
 i hoped someone could give me general advice on how to implement a term.
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.

 anon_g October 17, 2011 07:56

Quote:
 Originally Posted by boger (Post 327978) 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.
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. ;)

Quote:
 Originally Posted by adhiraj (Post 327980) Can you post the equation you want to code in?
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:
Quote:
 ...if by implicit handling the matrix A becomes easier to invert...
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...

 boger October 17, 2011 09:27

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?

 anon_g October 18, 2011 07:54

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.

 boger October 18, 2011 12:47

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.

 All times are GMT -4. The time now is 02:14.