Multi phase flow capabilities in the modules

[lindsayad]

Reason

Multi phase flow is very important physics for some applications and users including @joe61vette

Design

One possible design is adding additional terms to the Rhie-Chow interpolation as outlined in the references from https://github.com/joe61vette/ASP/pull/10#issuecomment-1087772360. Other designs can also be considered, given our ownership of relap-7 which I understand works well for some inputs and given our new common layer collaboration with SAM.

Impact

Powerful new physics capability in the moose repository

[joe61vette]

@joshuahansel @andrsd @GiudGiud As I understand it, rDG forms the basis of the THM and, as such, is used by both RELAP-7 and Sockeye. It would appear to be a good choice for a two-fluid type code as the face values of the fluxes are guaranteed to be bounded by the cell-centered values on either side of the face.

At present, this is not the case for the NS FVM approach using RC interpolation. This works fine as long as the two fluids are not strongly coupled. However, when the interfacial terms dominate, the face velocity can be seriously in error.

One quick question: does the rDG scheme as implemented in the THM have the pressure checkerboard problem when used for a nearly incompressible fluid like water?

[joshuahansel]

Hi Joe, I have never observed any pressure checkerboarding in any THM-based application. The rDG scheme (which in our case, is exactly equivalent to a FV scheme) has options for slope reconstruction. Sometimes for second-order slope reconstruction, you can get some oscillations at discontinuities, but I don't think that this is the checkerboard phenomena that you refer to.

[lindsayad]

I'm really surprised that a "naive" (e.g. doesn't have something like Rhie-Chow) collocated FV implementation would not show pressure checkerboarding for incompressible flows. Maybe it's the HLLC implementation?

[lindsayad]

@tanoret is a key stakeholder and developer in the two-phase world so pinging him here. Other key developers of two-phase flow capability include @lingzou and Rui (will mention once I have user-name)

[lindsayad]

@snschune is another important stakeholder

[joshuahansel]

I'm really surprised that a "naive" (e.g. doesn't have something like Rhie-Chow) collocated FV implementation would not show pressure checkerboarding for incompressible flows. Maybe it's the HLLC implementation?

Yes, I believe this would be an advantage of upwind methods over centered methods, since the upwind methods shouldn't have the cancellation of pressure terms from the adjacent pressure gradient terms.

[lindsayad]

So this is the phy.velocity_t_3eqn.i matrix for 3 elements for the first timestep, where liquid water is the fluid. The mass equation does indeed have a huge diagonal, which does make sense given the variable set choice (rhoA, rhoEA, rhouA) and the default absence of slope reconstruction (resulting in upwind). However, now it is the momentum equation that has a diagonal that is a couple of orders of magnitude less than the off-diagonal ... but it does still have a diagonal

Mat Object: () 1 MPI process
  type: seqaij
row 0: (0, 494.851)  (1, 0.000420325)  (2, 0.000203622)  (3, -492.464)  (4, -0.000421564)  (5, 0.500109) 
row 1: (0, 3587.32)  (1, 0.0031173)  (2, 0.017478)  (3, -3585.78)  (4, -0.00306953)  (5, 3.64145) 
row 2: (0, -788344.)  (1, -0.674162)  (2, 1604.19)  (3, 787541.)  (4, 0.674159)  (5, -799.767) 
row 3: (0, -492.464)  (1, -0.000421563)  (2, -0.499892)  (3, 988.262)  (4, 0.000843127)  (5, -0.000217456)  (6, -492.464)  (7, -0.000421564)  (8, 0.500109) 
row 4: (0, -3594.23)  (1, -0.00309301)  (2, -3.65583)  (3, 7180.)  (4, 0.00619587)  (5, 0.0143805)  (6, -3585.78)  (7, -0.00306953)  (8, 3.64144) 
row 5: (0, -787908.)  (1, -0.675841)  (2, -800.416)  (3, 366.996)  (4, 0.00168301)  (5, 1603.52)  (6, 787541.)  (7, 0.674158)  (8, -799.766) 
row 6: (3, -492.464)  (4, -0.000421563)  (5, -0.499891)  (6, 987.763)  (7, 0.000843127)  (8, 0.49947) 
row 7: (3, -3594.22)  (4, -0.003093)  (5, -3.65583)  (6, 7176.35)  (7, 0.00619585)  (8, 3.65273) 
row 8: (3, -787908.)  (4, -0.675842)  (5, -800.416)  (6, 1165.69)  (7, 0.00168357)  (8, 804.423) 

[lingzou]

I'm really surprised that a "naive" (e.g. doesn't have something like Rhie-Chow) collocated FV implementation would not show pressure checkerboarding for incompressible flows. Maybe it's the HLLC implementation?

Edit: I may misread your comments. You asked about incompressible flow, but the following I was answering compressible flow. My understanding is that HLLC is strictly for compressible flow, for which you could construct a Riemann problem on the interface of finite volume, and which is the core algorithm of the so-called 'Godunov' methods.

@lindsayad The Rhie-Chow is more for incompressible NS. It is an extension of SIMPLE (staggered-grid) to the collocated FV. The collocated FV for compressible flow is based on a different family of schemes (the Godunov method). The flux terms on the edge (interface between finite volumes) are only depending on the two neighboring finite volumes by solving an approximated Riemann problem. For first-order implementation, the 'monotonicity' is guaranteed; while if you go to higher-order, such as second-order, the monotonicity is preserved by carefully choosing the 'slope limiter' or 'flux limiter' as @joshuahansel explained.

Here is a good reference paper you may find useful A. Kurganov and E. Tadmor, New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations, Journal of Computational Physics, 160, 241–282 (2000)

[lingzou]

The link in the OP does not work by the way.

[lindsayad]

Yea @joe61vette's ASP repository is currently private. I may repaste some of the information there here as it's important for the discussion.

Yea we currently have both HLLC and Kurganov-Tadmor for doing compressible finite volume simulations in the navier_stokes module. And we use Rhie-Chow for incompressible and weakly compressible finite volume simulations in navier_stokes.

However, the thermal_hydraulics module uses HLLC with finite volumes for all classes of flow conditions, e.g. liquid water, ideal gas, etc. And when applying HLLC finite volume to incompressible flows there is not any checkerboarding as @joshuahansel stated and as I have found out it in my own testing today.

[joe61vette]

Hello everyone. If you would like access to the ASP repo, please let me know. I’ll let you know which branch to look at tomorrow as I am out of town w/o my laptop. Some background first.

ASP was a first attempt to build a two-fluid three-field (liquid, gas, drops) model for two-phase flow. It is a feasibility study for the possible future migration of NRC codes to the MOOSE framework. I decided to go with the new FVM model and so it extends classes from WCNSFV. Thanks to help from Alex and Guillaume, I made a good bit of progress with a relatively low level of effort. But….

Some problems ran well, other simple problems didn’t converge at all. This was traced to my explicit usage of RC interpolation. Basically one UO for each field. This meant that the interfacial terms (ie, drag as I was only modeling mass & energy) contributed to the residuals used in RC but did not have an implicit connection to the velocity of the other fields. So, there was no guarantee that the face velocities would be bounded by the cell-centered values. Indeed, both cell values could be positive while RC returned a negative value for the face. The result was that I halted work on ASP about 6 months ago while waiting on a possible two-phase implementation of RC (that is too deep in the framework for me to attempt at present). ASP is thus in a very unstable state as it was halted in the middle of some new models being implemented.

[joe61vette]

I had hoped to resume working on ASP this fall. So, I need to choose between WCNSFV, rDG or CFEM. Each have pros/cons. I like the idea of rDG due the face values of the fluxes always being bounded. But…

I am using the primitive variable set: (P, V, T) for single-phase flow. Add phasic volume fractions for 2 of the other 3 fields along with their V & T. For single-phase flow of a nearly incompressible fluid like water or liquid sodium, the derivative of the mass eq residual with respect to pressure is essentially zero. So, unless the scheme used to provide the flux values at the faces brings in a dependency upon pressure, I don’t know how this would work.

Would a two-phase version of HLLC need to be developed?

Should I wait for a segregated solver? Then the momentum eqs could be solved separately while mass & energy are still solved together (kind of like TRACE and RELAP5). Then velocity can be treated as a function of pressure as far as mass/energy are concerned.

Your advice is greatly appreciated. Thanks.

[joshuahansel]

In short, I'm not sure what the right answer is. My feeling is that it's going to take quite a bit of theoretical work/experimentation just to determine the right scheme for your equations/application.

I can speak a little to rDG and CFEM. Sockeye and RELAP-7 use a two-phase rDG (FV) scheme with a two-phase version of HLLC, but this corresponds to the 7-equation model. I'm guessing your equation set would need some theory development to derive an HLLC flux formulation. I don't know your equations, but I think they would need to be fully hyperbolic to apply any HLLC. Plus, it may not be the right method to begin with, just because it's geared toward fully compressible flows.

RELAP-7 used to use CFEM before we abandoned it for FV. The tricky part with CFEM is obviously stabilization. We implemented a variety of schemes using artificial viscosity. I think entropy viscosity was the most successful, but the problem was that it needed tuning some stabilization coefficients to achieve good results. The first-order scheme was reliable of course, but just like with Godunov in FV, artificial dissipation was sometimes quite excessive.

[lingzou]

@joe61vette A Riemann-solver based first-order 3D FVM for two-fluid six-equation model is available, see the paper published by Guojun: https://www.sciencedirect.com/science/article/pii/S0306454921007179 [edited: I said is available, but I should really say could. Guojun's paper deals with 1D, which can be extended to 3D]

I would argue that it is possible to have it work for the two-fluid three-field eight-equation two-phase flow model, but clearly, as @joshuahansel pointed out, much theoretical work is needed here

• first, prove that the equations system is hyperbolic, or to fix the equation system (like CATHARE did for six-equation) to make it hyperbolic,
• second, find the eigenvalues and eigenvectors for the 8by8 matrix
• develop the Riemann solver which produces the correct fluxes As soon as this is done, a first-order FV can be easily achieved as shown in the Guojun's paper.

Obviously, compared to way how mass/energy fluxes are treated with the staggered-grid FVM using upwind method (Ok, I never hide my love to staggered-grid method), the collocated FVM is more expensive. The trouble is the implementation in 3D since as we have discussed (@lindsayad), a dual mesh for the momentum equation is needed. [This is also the reason that staggered-grid FVM is PERFECT for 1D system code] Cost related to DG is another thing, for a second-order implementation, the # of DOF is pretty much doubled, compared to the staggered-grid method. For 1D system code, that may not be a big issue, for 3D, that is, and you have to seriously think about the cost vs. benefit.

rDG is a different story. [PS: we are not using the rDG word right here, it is misleading. THM uses collocated FVM, it is not rDG] rDG is a way to achieve higher-order, which requires the extension of element stencil to do the 'r' (reconstructed) part, the idea is in some sense how FV achieve higher-order, which I can argue that it will be challenging for existing MOOSE framework.

[lindsayad]

I like this paper for discussion of using rDG for incompressible flow problems.

@lingzou THM is using FV ... but it is a sub-class of rDG, e.g. rDG(P0P1), right? It uses piecewise constant basis functions (the P0 portion) and then reconstructs a first order polynomial using the surrounding neighbor elements (the P1 portion).

So that being said @joe61vette to say you want to use rDG vs. (collocated) FV is kind of a misnomer I think. I think what you are unhappy with in the collocated FV (sub-class of rDG) scheme is the Rhie-Chow piece in the advection discretization.

[lindsayad]

THM does not use Rhie-Chow in its advection discretization; it uses HLLC. That is the key distinction between navier_stokes WCNSFV vs. THM.

[lingzou]

@lindsayad Thanks for sharing the paper.

• Figure 1 of the paper is exactly what I did for some early version of Pronghorn, P0P1 (P0 for velocities and P1 for pressure), which is super stable but it later went to other directions.
• It is argued that FV could be categorized as a sub-class of FEM, Not trying to be picky, but it helps if we are more specific when referring to the methods.
• FV a sub-class of rDG, maybe. Because of my limit knowledge, I don't see much such statement, but the similarity between the two methods are obvious.
• I don't think THM FV is the typical P0P1. Typically, when you say P1, the P1 variable does have two nodal values plus the P1 shape functions. For THM FV, although we do reconstruct the linear function, but it is not the typical P1 variable.
• I have limit or practically zero knowledge on Rhie-Chow two-phase flow solving, so I can't help much here.

[joshuahansel]

we are not using the rDG word right here, it is misleading

Yes, I agree with this - this is why I try to call it FV instead of rDG. I think that the scheme could be described as P0P1 rDG, but there could be some subtle differences that I'm not aware of. In any case, "FV" is a better description, since as you said, there is no extension to higher order.

[lindsayad]

• I don't think THM FV is the typical P0P1. Typically, when you say P1, the P1 variable does have two nodal values plus the P1 shape functions. For THM FV, although we do reconstruct the linear function, but it is not the typical P1 variable.

There is a very unfortunate overlap of terminology here. When describing rDG methodologies it is very typical to use a PNPM notation in which 'N' refers to the order of the underlying DG basis functions and 'M' refers to the polynomial order that is reconstructed using nonlocal element information (e.g. neighbors).

And then there is notation like P2P1 for describing Taylor-Hood continuous finite elements on simplices where in this notation '2' means that the velocity bases are quadratic and the '1' means the pressure basis is linear.

[lindsayad]

@rui-hu now that I have your username I wanted to make you aware of this discussion

[lingzou]

@lindsayad thanks, so the P0P1 naming does make sense for THM FV. How about lets keep the FVM naming for THM to avoid future confusion (as @joshuahansel also suggested so). If needed, we can call it FV+HLLC when in the context of comparing it to FV+Rhie-Chow.

[lindsayad]

:+1: I agree that will help minimize confusion

[tanoret]

Hi all!

We've been doing some experimenting on two-phase flow since we need to accurately relatively-accurately model bubble transport in MSRs this FY. Here are some partial conclusions from our experiments that could help this discussion. Please tell me if you disagree with anything or if you think that we are missing something so that we can re-adjust.

• When approaching the incompressible limit, standard Reimann solvers have issues due to the pressure-independent density (basically, the two sound-speed-related eigenvalues in the Riemann problem go to infinity). We solved this issue by switching locally the way in which advective fluxes are computed to the cell faces, i.e., using Limited HLLC with Einfeldt signal speed estimates for high Mach numbers and standard arithmetic means interpolation (with RC momentum fluxes) for low Mach numbers. Currently, I set the cutoff manually at Ma=0.3 (there are probably smarter ways of doing this). The resulting set of equations is hyperbolic by borrowing CATHARE trick.
• In the few test cases that we have performed, there is not much difference in terms of accuracy/stability of the collocated vs. staggered grids. However, to stabilize face fluxes in the collocated formulation and preserve monotonicity we needed flux-splitting and weighted ENO schemes. We obtained relatively good results with WENO9-JS and TENO5
• In our problems, the standard reconstruction of the interface velocity via rDG or mixture models led to too diffusive results. For this purpose, we had to implement a cell-cutting algorithm for interface reconstruction, which allows us to model the effect of surface tension. This method also improved the results in 1D.

Here are some pictures of some preliminary results: Sod's Tube one-phase (left results with the code + 1 exact results):

Bubble undergoing a two-phase shock problem (left experiment - right simulations):

Merging droplets in a gas bed:

This testing is currently on our prototype code. We should be implementing it in MOOSE shortly once we are a bit more satisfied with the model and results. Our current main issue to solve in the short term is that we cannot appropriately model bubble break up...

[lingzou]

@tanoret That's fantastic progress!

I have some comments/questions as follows:

• My understanding is that there are a couple of issues related to Riemann problem associated with two-phase flow. As you mentioned the Riemann solver becomes poorly performed when approaching incompressible limit, and it becomes 'stiff' for the implicit solver because of the very large (tend to infinite) sound speed. The other issue is related to its performance when one phase goes to zero.
• For the 2d problem you are showing here, is it structure mesh or unstructured mesh you are using? Is WENO/TENO easy to be implemented MOOSE for unstructured mesh?

[tanoret]

• For the 2d problem you are showing here, is it structure mesh or unstructured mesh you are using? Is WENO/TENO easy to be implemented MOOSE for unstructured mesh?

These results are in a structured mesh but I have similar ones for unstructured ones. I am focusing currently mostly on structured meshes because my model for MCRE is currently a structured mesh :) Should not be impossible to implement WENO/TENO in MOOSE but I am unsure about how much will be the actual work for doing it in the 'right way'. I'd ask @lindsayad about this one.

[lingzou]

Thank you @tanoret

[lindsayad]

These do look like really exciting results @tanoret. You're rolling out awesome stuff as always!

[lindsayad]

@grmnptr I realized you hadn't been tagged here

[joe61vette]

@lindsayad A couple of things related to your post from yesterday. My concern about RC for my application is simply that developing the capability for two-phase multi-field flow requires a much deeper understanding of the new FV framework than I am ever likely to learn. Basically, I would need up to 3 RC UO's that are solved simultaneously so the the interfacial terms would be implicitly treated.

The paper (Pandere & Lou, 2016) is quite interesting. For their hybrid method, there seems to be two important parts:

1) Using rDG(PnPm) & CG(Pn) so that the pressure is defined at the nodes while the velocity is cell-centered.

2) Using a projection method so that an estimate for the velocity field (V*) is computed using the time lagged pressure. The divergence (or equivalently mass eq) is then solved for the necessary pressure correction. This appear to be very similar to what is done in TRACE & RELAP5.

The step that I am not sure of is that leading to eq. 3.2:

Specifically, it appears that the advection and shear stress terms are included in the calculation of V* but that the corrected velocity (Vn+1) only considers the difference in the pressure gradient. It seems like this could be implemented using the MOOSE rDG module. Would this be a reasonable approach? Would it require a segregated solver?

Thanks to all for the great discussion here.

[joe61vette]

After a bit more reflection, the pressure correction in the above hybrid scheme is too simple. Basically, it would ignore the coupling of the phasic velocities. In COBRA-TF and other codes I have worked with, we solved a 3x3 matrix eq to get estimates of the derivative of the phasic velocities wrt pressure that we then used in the mass/energy eqs along with the tentative velocities. So, very much a segregated approach. In essence:

• solve the 3 phasic momentum eqs for tentative velocities
• solve the 3x3 for every face (a material property in MOOSE?)
• solve all 6 mass/energy eqs together

I believe that I either need to go back to the WCNSFV approach or switch to CFEM.

For my applications, flows in reactor vessels for PWRs, BWRs and their SMR versions, I only need a weakly compressible form. The only place where a sonic velocity occurs is at a break and that is in a 1D component and usually treated by an empirical model. Water hammer is another matter but that is way down the list of priorities.

[lindsayad]

After going to war with the Water97FluidProperties class for the last three days to add support for running it with conservative variables, the matrix I get out on the THM test ends up looking very similar to that when using a StiffenedGas :laughing: (see https://github.com/idaholab/moose/issues/22637#issuecomment-1314370338 for that matrix)

      SVD: condition number 2.441843322518e+11, 0 of 9 singular values are (nearly) zero
      SVD: singular values:
             1.501529298631e+06   1.038973912264e+06   3.701408259763e+05   4.668696900287e+01   1.192502051978e+01  
             2.869699441978e+00   8.516552976210e-06   7.377570511006e-06   6.149163153854e-06  
      0 Linear |R| = 2.824152e-04
      1 Linear |R| = 4.153350e-13
Mat Object: () 1 MPI process
  type: seqaij
row 0: (0, 583.966)  (1, 0.000186027)  (2, 0.000510781)  (3, -581.157)  (4, -0.000186357)  (5, 0.499837) 
row 1: (0, 4221.37)  (1, 0.00139813)  (2, 0.0180354)  (3, -4224.11)  (4, -0.00135452)  (5, 3.63304) 
row 2: (0, -833398.)  (1, -0.267009)  (2, 1437.37)  (3, 832664.)  (4, 0.267006)  (5, -716.151) 
row 3: (0, -581.162)  (1, -0.000186357)  (2, -0.500163)  (3, 1165.65)  (4, 0.000372714)  (5, 0.00032579)  (6, -581.154)  (7, -0.000186357)  (8, 0.499837) 
row 4: (0, -4233.51)  (1, -0.00136985)  (2, -3.64972)  (3, 8457.57)  (4, 0.00275771)  (5, 0.0166714)  (6, -4224.07)  (7, -0.00135452)  (8, 3.63302) 
row 5: (0, -833558.)  (1, -0.26775)  (2, -717.614)  (3, 883.188)  (4, 0.000743158)  (5, 1437.09)  (6, 832654.)  (7, 0.267005)  (8, -716.147) 
row 6: (3, -581.157)  (4, -0.000186357)  (5, -0.500163)  (6, 1165.14)  (7, 0.000372715)  (8, 0.499976) 
row 7: (3, -4233.46)  (4, -0.00136985)  (5, -3.64971)  (6, 8453.85)  (7, 0.0027577)  (8, 3.64832) 
row 8: (3, -833547.)  (4, -0.267749)  (5, -717.609)  (6, 1600.23)  (7, 0.000743599)  (8, 721.21) 

[lindsayad]

Specifically, it appears that the advection and shear stress terms are included in the calculation of V* but that the corrected velocity (Vn+1) only considers the difference in the pressure gradient. It seems like this could be implemented using the MOOSE rDG module. Would this be a reasonable approach? Would it require a segregated solver?

Using a predictor-corrector scheme is not necessary. We could explore it now that we have the capability to solve multiple (nonlinear) systems on the same mesh and after @grmnptr develops Executioners/Executors for the classic segregated CFD schemes. (See issue #22356 and associated pull request #22699). But I see the use of a predictor-corrector solve scheme and the choice of discretization as two separate concepts. I personally like the idea of a discretization that is innately LBB stable, such as:

• Taylor-Hood continuous finite elements
• hybrid rDG-CFEM discretization proposed by Pandare and Luo
• staggered grid finite volume

However, it seems that @tanoret is getting really good results with collocated finite volume for multi-phase and he is using classic arithmetic mean-based interpolations and Rhie-Chow for velocity for Mach numbers below 0.3 (HLLC-like scheme for Mach numbers above 0.3).

In @joe61vette's ASP issue he shared several papers including this 2021 Hanimann article and this 2015 article by Moukalled. They outline a couple of modifications to Rhie-Chow that we would have to make in navier_stokes in order to run multi-phase (these are shown clearly in equation 25 of the Moukalled paper):

• consideration of the phase fractions
• including a contribution to RC dependent on the other phase's velocities, e.g. an "off-diagonal" velocity-to-velocity coupling

@tanoret do you have modifications like this for RC for your low Mach number simulations?

[joe61vette]

Would the hybrid rDG-CG scheme w/o the predictor-corrector step work if primitive variables were used? For single-phase nearly incompressible (ie, solid water), the derivative of the mass residual with respect to pressure (the non-linear variable) would essentially be zero. The off-diagonal terms would then be derivatives with respect to velocity which in turn is dependent upon the pressure gradient.

Thanks.

[lindsayad]

Yes it could work. You can always introduce the same kind of on-diagonal ellipticity for the pressure equation in a monolithic scheme that you can in a segregated scheme; but you could also choose not to. If you do add the ellipticity then hopefully with a good choice of field split preconditioner for the monolithic system you can reproduce the benefits of the segregated solve behavior.

If you do not add the ellipticity (or sometimes even if you do), you will need to either go with LU or field split with a Schur complement or sometimes ASM with a fair amount of overlap will work too.

Taylor-Hood finite elements with no pressure-stabilized Petrov-Galerkin have no on-diagonal for the pressure equation.

[lindsayad]

Before making any decisions @joe61vette I would wait for @tanoret to respond about whether he is doing any RC modifications in his multi-phase prototype. We will eventually want to incorporate the work he is doing in his prototype into navier_stokes and it's possible the translation won't be that difficult and that it could be applied to your use cases

[tanoret]

@lindsayad I am doing standard RC for low Mach and that is working fine for my toy problems. However, I like the mods in Hanimann. Should be worth trying. I can do a rapid test on my side and report if you want. What do you think?

[lindsayad]

That sounds great. Maybe @joe61vette can outline a problematic input

[joe61vette]

@lindsayad @tanoret Sorry this response is so late. Some of the really difficult test cases were fill/drain involving phase appearance and disappearance. However, I also had failures where the face velocity was not bounded by the cell-center values for some more normal two-phase cases (don't remember which right now).

First, however, I might be able to explain why Mauricio hasn't seen a problem like those I encountered. Basically while both problems are two-phase, they are actually very different. In the problem Mauricio showed (Mauricio please correct me if I'm wrong here):

• there was one vapor structure (a large Taylor cap type bubble)
• the bubble diameter was much larger then the element size
• Level set was used to find the liquid/vapor interface
• At the liquid/vapor interface, a no-slip condition was applied between the two velocity fields
• There is no equivalent to the interfacial drag force term that depends on the phasic relative velocity

Whereas in my test problems:

• there can be hundreds of bubbles in a single 3D element
• the bubble diameter is much smaller than the element size
• The vapor-liquid interface is not defined
• Rather a volumetric force term due to interfacial drag is applied. This is typically of the same magnitude as the gravitational force term and has the form Fi = Ci Ai (Vgas - Vliq). Where Ci is a drag coefficient (value ~0.5), Ai is the interfacial area density that depends on flow regime and phase volume fraction (value ~100 m^2/m^3).

This is similar to the way two-phase flow is treated in codes like RELAP5, TRACE & CTF. Many of my test problems were basically a balance between the gravity term and the interfacial drag term. The problem I had with the simple-minded way I used RC interpolation (one UO for each field) was that while the interfacial drag term contributed to the residual part of the RC interpolation, the derivative part would only have a dependence on that phase velocity and not the other. So, I would get situations where the face velocity for a phase was negative while both of the cell-centered values (on each side of that face) were positive. This caused a death spiral where the time step got smaller and the residual larger (if I remember correctly).

For a robust two-phase model for my type of problems that uses the new FVM capability in Moose, I believe that the interpolated velocities for all of the fields (I had 3, liquid, vapor & drops) would need to be solved for together so that the interfacial drag term is implicitly handled in the RC interpolation.

I can invite Mauricio to be an ASP collaborator if he would like to see my code & test problems. The master branch is where I stopped with the FVM implementation last March.

[lingzou]

@joe61vette Thanks for sharing more details.

I'd like to contribute some discussions toward the two-phase flow problems. I am more familiar with the classical two-phase flow models (two-fluid two(three)-field two-phase flow model used in RELAP, TRACE, CTF) you are working on, while less on the model @tanoret works on.

• In the classic two-phase flow model, the interfacial drag term in the momentum equation is one of the largest source term (sometime the single largest term), so intuitively and based on our past experience [1,2], it has to be treated implicitly to make them behave, especially when one phase is (dis)appearing.
• The question is: isn't it being implicitly treated in the cell center velocities? If the phasic velocities behave correctly in cell centers, naively I'd expect that they should also behave on interfaces between two cells given that they are somewhat a weighted averaged value. (Ok, I might be too naive here).
• Phase appearing and disappearing are quite tricky to deal with base on our past experience [1,2]. The velocity and temperature of the (dis)appearing phase needs special treatments for them to behave if solved implicitly. As a matter of fact, we rely on the interfacial drag/heat transfer to make them behave during phase (dis)appearing. Again, cell centered method is a bit beyond my knowledge base, but I would expect that some special treatments are also needed dealing with phase (dis)appearing. I would suggest to get the 'normal' (void fraction not too close to 0 or 1) two-phase flow working before tackling the phase (dis)appearing issue.

[1] Computers and Fluids 129 (2016) 179–188 [2] Progress in Nuclear Energy 88 (2016) 198-210

[joe61vette]

@lingzou Hi Ling. To learn about RC interpolation see Moukaled's book (ref in Moose/NS module/FVM). But in layman's terms, it basically solves a pseudo-momentum equation to get the face velocity based on the residuals and their derivatives from the adjacent cell-centers. But in my simple implementation I had a separate RC UO for each field that was not solved simultaneously. It was kind of a cross your fingers and hope effort as actually solving two UO's together with caching etc. would require a knowledge of the framework that is far beyond me.

[lingzou]

@joe61vette Thanks. I see now. So the face velocity is not simply a weighted average, while it is indeed the solution of the pseudo-momentum equation using neighboring cell velocities and provisional pressure. The separate RC UO missing the important pieces of interfacial drag term (as discussed, a very important piece), and thus it generates not-well-behaved values. I guess, besides the moose infrastructure development needed for simultaneously solving the face velocities, it is quite possible that there need some sort of theory development on the coupled phasic pseudo-momentum equation involving the interfacial drag terms. This could be an interesting R&D topic and I wonder if such R&D has already been done by someone.

[lindsayad]

The Hanniman article has off-diagonal contribution to one phase's Rhie-Chow from other coupled phases. I think this is the critical piece we are missing

[lingzou]

Ahh...neat, I missed that paper in the long discussion.

[alimostafavi24]

Is simulating laminar multiphase flows still in demand at MOOSE? I have done this by coupling phase field and INS, and tested it on a number of benchmark problems. It works well. I can send a report on that.

[lindsayad]

I think any multiphase flow capability is great!

[GiudGiud]

I think as long as it’s well documented as to what it can and cannot do we should be good to include this in moose. Is there a particular name for this coupled NS phase field problem/equations? To decide where this goes

[alimostafavi24]

I am using coupled Cahn-Hilliard Navier-Stokes equations. I got reasonable agreement in the problems that I tried to simulate.

[GiudGiud]

I think this will fit in the combined module. NS doesn't depend on PF, and PF does not depend on NS.

Either that or we start a multiphase_flow module. @lindsayad ?

[lindsayad]

I think combined makes the most sense for now

[alimostafavi24]

That's good! I am PhD student at U of Arizona. How does code addition to the MOOSE environment works here? Is there any way that our research group can collaborate on this with MOOSE?