CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

add a phase diffusion term in both continuity and monmentum equations for twoPhEuler

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 6, 2020, 10:33
Default
  #21
Senior Member
 
kimy
Join Date: Mar 2019
Posts: 150
Rep Power: 4
qi.yang@polimi.it is on a distinguished road
It is so weird. When I used the mesh with larger boundary layer thickness, the result is reasonable and accurate. After I decreased the boundary layer thickness, it was crashed cause the bounding k and epsilon were infinite. Who can give me some suggestions? Thanks!
qi.yang@polimi.it is offline   Reply With Quote

Old   February 7, 2020, 03:22
Default
  #22
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 15
Cyp is on a distinguished road
Quote:
Originally Posted by qi.yang@polimi.it View Post
I think I solved the problem with the following code:
- fvc::div(fvc::grad(alpha1)*rho1*phase2.turbulence( ).nut()/sigma*U1)
it was compiled successfully and I ran the test case. The results seem resonable.
However, I doubt that whether I need write the code like:
- fvc::Sp(div(fvc::grad(alpha1)*rho1*phase2.turbulen ce().nut()/sigma, U1))

Because I found the following two codes will lead different results.
solve( fvm::laplacian(k, T) + fvc::Sp(sp, T) );
solve( fvm::laplacian(k, T) + sp*T );

However, when I wrote - fvm::Sp(div(fvc::grad(alpha1)*rho1*phase2.turbulen ce().nut()/sigma, U1)), it failed to be compiled.
it failed to compile because the operator fvc::div() needs fluxes at the face center (so, a surfaceScalarField and not a volVectorField). When you do fvc::grad(alpha1), you create a volVectorField. when you write phase2.turbulence().nut(), you generate a volScalarField, where your values are computed at the cell centers. This is why you have to interpolate.
Cyp is offline   Reply With Quote

Old   February 7, 2020, 03:40
Default
  #23
Senior Member
 
kimy
Join Date: Mar 2019
Posts: 150
Rep Power: 4
qi.yang@polimi.it is on a distinguished road
Quote:
Originally Posted by Cyp View Post
it failed to compile because the operator fvc::div() needs fluxes at the face center (so, a surfaceScalarField and not a volVectorField). When you do fvc::grad(alpha1), you create a volVectorField. when you write phase2.turbulence().nut(), you generate a volScalarField, where your values are computed at the cell centers. This is why you have to interpolate.
Thanks. In your opinion,
- fvc::div(fvc::grad(alpha1)*rho1*phase2.turbulence( ).nut()/sigma*U1) this code is reasonable or not? After I added this code in the equation, it can be compiled and run well for the corse mesh. However, when I use a lower thickness of boundary layer, the simulation crashed.
qi.yang@polimi.it is offline   Reply With Quote

Old   February 7, 2020, 03:59
Default
  #24
Senior Member
 
kimy
Join Date: Mar 2019
Posts: 150
Rep Power: 4
qi.yang@polimi.it is on a distinguished road
-fvm::div(fvc::snGrad(alpha1)*mesh.magSf()*fvc::int erpolate(rho1*phase2.turbulence().nut())/sigma, U1)

I changed it to be like this, it compiled. But the simulation is also crashed for the fine boundary layer mesh.
qi.yang@polimi.it is offline   Reply With Quote

Old   February 7, 2020, 04:23
Default
  #25
Senior Member
 
kimy
Join Date: Mar 2019
Posts: 150
Rep Power: 4
qi.yang@polimi.it is on a distinguished road
Quote:
Originally Posted by qi.yang@polimi.it View Post
-fvm::div(fvc::snGrad(alpha1)*mesh.magSf()*fvc::int erpolate(rho1*phase2.turbulence().nut())/sigma, U1)

I changed it to be like this, it compiled. But the simulation is also crashed for the fine boundary layer mesh.
My obejective is to add the phase dispersion term in both U and p equations
Please help me to verify the codes I wrote are right or not?
In Uequation I wrote like this:
U1Eqn =
(
fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
+ MRF.DDt(alpha1*rho1, U1)
+ - fvm::laplacian(alpha1*rho1*phase2.turbulence().nuE ff(), U1)
- fvc::div(fvc::grad(alpha1)*rho1*phase2.turbulence( ).nut()/sigma*U1)

While in p equation I added one code in p equacomp1 and equacomp2:
pEqnComp1 =
(
contErr1
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)) - fvc::laplacian(phase2.turbulence().nut()/sigma,alpha1);
Attached Images
File Type: jpg Picture2.jpg (17.6 KB, 21 views)
qi.yang@polimi.it is offline   Reply With Quote

Old   February 8, 2020, 23:30
Default
  #26
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 783
Rep Power: 14
sharonyue is on a distinguished road


The equation looks not correct. Perhaps the first term should be a \nabla\cdot instead of \nabla. Meanwhile,
1) to add a diffusion term, you better to add that after MULES to ensure boundedness. See how it does in driftFluxFoam.
2) There is already a turbulent dispersion force in E-E model and it was implemented already, see "turbulent dispersion force".
3) The momentum interfacial exchange term needs to be addressed in UEqn, not pEqn.
4) div(grad()) employs an extended stencil instead of a compact stencil, which means it may introduce possible oscillations.
__________________
My OpenFOAM algorithm website: http://dyfluid.com
By far the largest Chinese CFD-based forum: http://www.cfd-china.com/category/6/openfoam
By far the largest WECHAT CFD media ID: cfdresearch
sharonyue is offline   Reply With Quote

Old   February 10, 2020, 03:37
Default
  #27
Senior Member
 
kimy
Join Date: Mar 2019
Posts: 150
Rep Power: 4
qi.yang@polimi.it is on a distinguished road
Thanks a lot Dongyue. The equations that I wrote here are not accurate and you are right.
We will not use the dispersion force defined in the openfoam because our group provided a new two fluid model that add two dispersion terms in both continuity and momentum equations which makes numerical stability and fast calculating.
In terms of "div(grad()) employs an extended stencil instead of a compact stencil", how to solve the problem you mentioned. However we have to add the div(grad(alpha))
qi.yang@polimi.it is offline   Reply With Quote

Reply

Tags
multiphaseflow, twophaseeulerfoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 16:04.