Gravitational flux in multiphase equations

ELwardi · April 24, 2019, 17:55

Hello,

I'm stuck with a good-for-nothing-but-working-fine solver with no gravity capabilities for around a week now. I can't find where the problem is, so, please, please help me out

.

I'm solving two-phase compressible flow in porous media, so, my first attempt to calculate the (water) flux was:

Code:

// Kwf is a surfaceScalarField as it's important to do this
// because of very specific "interpolation" schemes that 
// have to be used
phiwG = (Kwf * rhowf * g) & mesh.Sf();

// Then, in the equation, I do
fvScalarMatrix wEqn
(
   Coeff1*fvm::ddt(water.alpha())
   + fvc::div(phiwG)
   + sources
);

The handling of the other phase is similar, only I have fvm::laplacian(-Ko, p) instead of the ::ddt term, and both equations are in "dimless/dimTime" dimensions.

So, when I set g = (0 0 0), I get fine results. But when g = (0 0 -9.81), everything gets messed up, ie. bounded variables get out of their bounds.

In my many attempts to debug the problem, I guessed either:

- My calculation of phiwG is wrong (You can't mess up g, so Kwf or rhowf, but they are fine).
- My calculation of other fluxes is wrong, and after inspection

Code:

// This is not valid OpenfOAM code
// units missmatch
fvc::laplacian(-Kw, p)/(fvc::div(phiwG)+VSMALL)

was ranging between 40~70 which is not that bad. It just means my pressure forces are ~50 times are stronger than gravity.
- rhow strongly depends on the pressure (which is a main variable), and I don't have enough loops

. But in an incompressible case, this ceases to be true and I can't even get that case working

Most of the other approaches, that are used by standard solver, ex. taking out rho (with snGrad) are not suitable in this case because even if the flow is reduced to be incompressible, difference of rho of both fluids still matters.

If anyone has even a hint of a suggestion, please don't hold back.
Best wishes, and thanks in advance.

wyldckat · April 28, 2019, 16:07

Quick answer: I came here after receiving the PM you sent me, but I'm not understanding what's going on either.

Because:

I don't know on which solver(s) you've based yourself on.
I don't know what the field "Kwf" actually stands for, so it could be anything from potatoes to the meaning of life...
I don't know what the original equation is meant to be.
I didn't understand if without "rhowf" the solver even runs...
I don't know how "rhowf" was created in the first place...

To me, there are too many unknowns... it could something as simple as having defined or updated the variables incorrectly, if they are based on tmp variables.

The other thing that comes to mind is that the use o "mesh.Sf()" might be wrong, given that it's the product of face area with face normal, so perhaps you were not meant to have the face area?

Perhaps the "div(phiwG)" term is utterly incorrect, in the sense that maybe it should be something like "div(phi,g)".

Or perhaps the wrong discretization scheme is being used... from the description you've given, nothing other than "upwind" should perhaps be used?

Last but not least, perhaps you're not familiar with this page: https://openfoamwiki.net/index.php/HowTo_debugging ?

ELwardi · April 29, 2019, 13:52

Thanks @wyldckat for your informative input; I will try to clarify things a bit; I know I've been over protective about my solver

.

1. The problematic solver is a coupled solver, similar, to some degree, to pUcoupledFoam in foam-extend.

2. "Kwf" is the face interpolation of what you, CFD guys, may call diffusion. It's closer to a "transmissiblity" really. The important thing that it has some factors that must be interpolated using a specific scheme (upwind & harmonic).

3. Well, I got ambitious and tried to solve the most complex form of equations there is. Then I went back and simplified them to two-phase, incompressible BlackOil equations:

$\frac{\partial }{\partial t}(\phi S_w) +\nabla(-K_w \nabla p) +\nabla(K_w \rho_w g \nabla z) +\nabla(K_w \nabla p_{cow}) + q_w = 0$

$\nabla(-K_t \nabla p) +\nabla(K_w \rho_w g \nabla z) +\nabla(K_o \rho_o g \nabla z) +\nabla(K_w \nabla p_{cow}) + q_w + q_o = 0$

where

, and one may treat $\phi$ as constant porosity. The main variables are

and

. $p_{cow}$ is the capillary pressure (treated explicitly).

4. I tried neutralizing "rho" effect before by setting both densities to 1 but the solver fails anyway.

5. 'rhowf' is the face interpolation of water density (which I keep constant through the domain for the purpose of testing).

... And of course I know about GDB, and gdbFoam; I don't think I would get this far if I didn't use them

.

Your suggestion that "div(phiwG)" term should be expressed in an other form got me thinking, so I went back and tested it on a "segregated" version of the solver (But the timeStep is so small It's painful to watch) and it worked fine. That's why I assumed It would work in a coupled solver in the first place.

I'm very careful with discretization schemes and all my face fluxes need to be multiplied by mesh.magSf() somehow if I want to "div" them as the end result should be in dimless/dimTime dimensions

.

Thanks for reading this long thread, but I'm starting to think my coupling of these equations is not entirely correct so I'm re-building on top of the working segregated solver; I'll let you know how that goes.

wyldckat · April 29, 2019, 21:29

Quick answer: Unfortunately this really goes outside of my current level of experience.

One of the things that popped into my mind was that perhaps "phiwG" might be accounting as many times as faces per cell, for the actual divergence value... although it doesn't make much sense, given that there are other discretization terms handled this way...

Take a look into the "nextRelease" branch in foam-extend 4.0, because they have additional solvers there for coupling...

And I took a second look to the equations you provided and I'm confused as to which variables you're actually solving for here... and I got the feeling that you need to restructure how the coupling is being done.
From what I can briefly see in the coupled solver you mentioned, there are terms that have to be coupled manually because they are present in more than one equation, so they have to be handled in a special way; others had to be constructed based on the other equation.

This may be too late for you to attent, but just in case it's still on time, you might want to attend this, since the people that worked on this kind of problem-solving will likely be there or nearby: NUMAP FOAM Summer School 2019

The other group of people I believe have worked with this kind of thing are the ones who develop Caelus-CML.

ELwardi · April 30, 2019, 17:53

Thanks a lot Mr. @wyldckat for your insightful assistance on this subject.

As we suspected, the coupling of the

Code:

fvc::div(phiwG)

term was not handled well; Apparently I forgot to "insert the coupling to the block matrix".

Without your assistance, I know I would never caught the bug

. Thanks again.

Now, things work as I expected, well, I still have some issues, but nothing that GDB can't handle.

And by the way,

Quote:

The other group of people I believe have worked with this kind of thing are the ones who develop Caelus-CML.

I haven't got such on-the-spot orientation to a group of developers in years

.

wyldckat · April 30, 2019, 20:24

Quote:

Originally Posted by ELwardi

Thanks a lot Mr. @wyldckat for your insightful assistance on this subject.

As we suspected, the coupling of the

Code:

fvc::div(phiwG)

term was not handled well; Apparently I forgot to "insert the coupling to the block matrix".

Without your assistance, I know I would never caught the bug

. Thanks again.

You're welcome and I'm glad I was able to help!
I'm not entirely certain how at 1AM I was able to pinpoint that being one of the likely issues... we both got lucky

April 24, 2019, 17:55	Gravitational flux in multiphase equations	#1
ELwardi Member Elwardi Fadeli Join Date: Dec 2016 Location: Boumerdes, Algeria Posts: 41 Rep Power: 9	Hello, I'm stuck with a good-for-nothing-but-working-fine solver with no gravity capabilities for around a week now. I can't find where the problem is, so, please, please help me out . I'm solving two-phase compressible flow in porous media, so, my first attempt to calculate the (water) flux was: Code: // Kwf is a surfaceScalarField as it's important to do this // because of very specific "interpolation" schemes that // have to be used phiwG = (Kwf * rhowf * g) & mesh.Sf(); // Then, in the equation, I do fvScalarMatrix wEqn ( Coeff1*fvm::ddt(water.alpha()) + fvc::div(phiwG) + sources ); The handling of the other phase is similar, only I have fvm::laplacian(-Ko, p) instead of the ::ddt term, and both equations are in "dimless/dimTime" dimensions. So, when I set g = (0 0 0), I get fine results. But when g = (0 0 -9.81), everything gets messed up, ie. bounded variables get out of their bounds. In my many attempts to debug the problem, I guessed either: - My calculation of phiwG is wrong (You can't mess up g, so Kwf or rhowf, but they are fine). - My calculation of other fluxes is wrong, and after inspection Code: // This is not valid OpenfOAM code // units missmatch fvc::laplacian(-Kw, p)/(fvc::div(phiwG)+VSMALL) was ranging between 40~70 which is not that bad. It just means my pressure forces are ~50 times are stronger than gravity. - rhow strongly depends on the pressure (which is a main variable), and I don't have enough loops . But in an incompressible case, this ceases to be true and I can't even get that case working Most of the other approaches, that are used by standard solver, ex. taking out rho (with snGrad) are not suitable in this case because even if the flow is reduced to be incompressible, difference of rho of both fluids still matters. If anyone has even a hint of a suggestion, please don't hold back. Best wishes, and thanks in advance.

April 28, 2019, 16:07		#2
wyldckat Retired Super Moderator Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,982 Blog Entries: 45 Rep Power: 128	Quick answer: I came here after receiving the PM you sent me, but I'm not understanding what's going on either. Because: I don't know on which solver(s) you've based yourself on. I don't know what the field "Kwf" actually stands for, so it could be anything from potatoes to the meaning of life... I don't know what the original equation is meant to be. I didn't understand if without "rhowf" the solver even runs... I don't know how "rhowf" was created in the first place... To me, there are too many unknowns... it could something as simple as having defined or updated the variables incorrectly, if they are based on tmp variables. The other thing that comes to mind is that the use o "mesh.Sf()" might be wrong, given that it's the product of face area with face normal, so perhaps you were not meant to have the face area? Perhaps the "div(phiwG)" term is utterly incorrect, in the sense that maybe it should be something like "div(phi,g)". Or perhaps the wrong discretization scheme is being used... from the description you've given, nothing other than "upwind" should perhaps be used? Last but not least, perhaps you're not familiar with this page: https://openfoamwiki.net/index.php/HowTo_debugging ? __________________ OpenFOAM: FAQ \| Getting started Forum: How to get help, to post code/output and forum guide Read this before sending me PM

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
NS Equations for multiphase flow in porous media	quir	OpenFOAM Running, Solving & CFD	3	December 26, 2020 12:25
Multiphase Homogeneous Model equations	marcoymarc	CFX	3	July 3, 2013 20:11
Question about heat transfer coefficient setting for CFX	Anna Tian	CFX	1	June 16, 2013 07:28
Flux splitting	Dr B.M. Smith (Smith)	OpenFOAM Running, Solving & CFD	19	January 9, 2013 05:07
Replace periodic by inlet-outlet pair	lego	CFX	3	November 5, 2002 21:09

April 29, 2019, 13:52		#3
ELwardi Member Elwardi Fadeli Join Date: Dec 2016 Location: Boumerdes, Algeria Posts: 41 Rep Power: 9	Thanks @wyldckat for your informative input; I will try to clarify things a bit; I know I've been over protective about my solver . 1. The problematic solver is a coupled solver, similar, to some degree, to pUcoupledFoam in foam-extend. 2. "Kwf" is the face interpolation of what you, CFD guys, may call diffusion. It's closer to a "transmissiblity" really. The important thing that it has some factors that must be interpolated using a specific scheme (upwind & harmonic). 3. Well, I got ambitious and tried to solve the most complex form of equations there is. Then I went back and simplified them to two-phase, incompressible BlackOil equations: $\frac{\partial }{\partial t}(\phi S_w) +\nabla(-K_w \nabla p) +\nabla(K_w \rho_w g \nabla z) +\nabla(K_w \nabla p_{cow}) + q_w = 0$ $\nabla(-K_t \nabla p) +\nabla(K_w \rho_w g \nabla z) +\nabla(K_o \rho_o g \nabla z) +\nabla(K_w \nabla p_{cow}) + q_w + q_o = 0$ where , and one may treat $\phi$ as constant porosity. The main variables are and . $p_{cow}$ is the capillary pressure (treated explicitly). 4. I tried neutralizing "rho" effect before by setting both densities to 1 but the solver fails anyway. 5. 'rhowf' is the face interpolation of water density (which I keep constant through the domain for the purpose of testing). ... And of course I know about GDB, and gdbFoam; I don't think I would get this far if I didn't use them . Your suggestion that "div(phiwG)" term should be expressed in an other form got me thinking, so I went back and tested it on a "segregated" version of the solver (But the timeStep is so small It's painful to watch) and it worked fine. That's why I assumed It would work in a coupled solver in the first place. I'm very careful with discretization schemes and all my face fluxes need to be multiplied by mesh.magSf() somehow if I want to "div" them as the end result should be in dimless/dimTime dimensions . Thanks for reading this long thread, but I'm starting to think my coupling of these equations is not entirely correct so I'm re-building on top of the working segregated solver; I'll let you know how that goes.

April 29, 2019, 21:29		#4
wyldckat Retired Super Moderator Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,982 Blog Entries: 45 Rep Power: 128	Quick answer: Unfortunately this really goes outside of my current level of experience. One of the things that popped into my mind was that perhaps "phiwG" might be accounting as many times as faces per cell, for the actual divergence value... although it doesn't make much sense, given that there are other discretization terms handled this way... Take a look into the "nextRelease" branch in foam-extend 4.0, because they have additional solvers there for coupling... And I took a second look to the equations you provided and I'm confused as to which variables you're actually solving for here... and I got the feeling that you need to restructure how the coupling is being done. From what I can briefly see in the coupled solver you mentioned, there are terms that have to be coupled manually because they are present in more than one equation, so they have to be handled in a special way; others had to be constructed based on the other equation. This may be too late for you to attent, but just in case it's still on time, you might want to attend this, since the people that worked on this kind of problem-solving will likely be there or nearby: NUMAP FOAM Summer School 2019 The other group of people I believe have worked with this kind of thing are the ones who develop Caelus-CML.