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

New solver for two-phase flows with phase-change heat transfer

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

Like Tree17Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 16, 2016, 08:52
Smile New solver for two-phase flows with phase-change heat transfer
  #1
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Hi all,

Our group recently released a two-phase flow solver for phase-change heat transfer based on the interFoam code. We have posted the code, explanation of the solver algorithm, sample cases, and installation instructions on gitHub: https://github.com/MahdiNabil/CFD-PC.

A number of phase-change heat transfer solvers have been described in the literature. The standard approach is to incorporate closure models that apply heating/alpha1/dilatation source terms at the interface. However, almost all studies reimplemented these formulations from scratch, which adds a lot of software development and validation effort. In this solver, we have added support for runtime selectable phase-change closure models (much like turbulence models, viscosity models, etc.). The base code provides 5 different closure models to choose from, and it is relatively simple to implement new ones. We hope that this will help end users more quickly investigate problems of interest.

The current version of the code includes tutorial cases for the horizontal film condensation (Stefan problem), smooth falling film condensation (Nusselt problem), wavy film condensation, bubble nucleation on a wall (picture below), and rising bubble condensation.

We also add support to select between the default surface tension force model in OpenFOAM, and the Sharp Surface Tension force model of Raeini et al. (often more accurate). We also add a switch so that you can select between the standard split hydrostatic/dynamic pressure formulation, or combined pressure. The second case is useful for specifying pressure jump conditions in cyclic domains (used for the wavy film condensation case here).

We hope this code is helpful to you all. Please let us know your thoughts so that we can improve/extend the code and refine the documentation.

liquidspoon is offline   Reply With Quote

Old   April 7, 2016, 21:28
Unhappy
  #2
New Member
 
Julia
Join Date: Jan 2016
Posts: 6
Rep Power: 10
Wangzhaoting is on a distinguished road
hi,Alex

I'm interested in your new slover. I try to simulate temperature field of two phase flow by your solver,but I'm failed.errors are listed:

FOAM FATAL IO ERROR:
keyword phase2 is undefined in dictionary "/home/wzt/OpenFOAM/wzt-3.0.0/run/re1-0.003L1g3/constant/transportProperties"

file: /home/wzt/OpenFOAM/wzt-3.0.0/run/re1-0.003L1g3/constant/transportProperties from line 23 to line 41.

From function dictionary::subDict(const word& keyword)
in file db/dictionary/dictionary.C at line 666.

I don't know where the problem is.please help me!
Wangzhaoting is offline   Reply With Quote

Old   April 15, 2016, 11:32
Default
  #3
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Thank you for the message - sorry for my delayed response.

We built the code for OpenFOAM 2.4.0, and haven't updated for 3.0 yet. One of the changes from the versions is that the phases are called phase1 (liquid) and phase2 (gas, or other liquid) in transportProperties file. This version doesn't support other names for the two phases. It looks like that may be the error you are getting.

Could you please try copying one of the constant/transportProperties files from a tutorial case (maybe NusseltSmooth) and use that as a template.
liquidspoon is offline   Reply With Quote

Old   April 17, 2016, 09:50
Default
  #4
New Member
 
Julia
Join Date: Jan 2016
Posts: 6
Rep Power: 10
Wangzhaoting is on a distinguished road
Hi ,Alex

Thank you for your help. It works now! But there is another problem .The object of the simulation is about 2.6-meters long pipeline with radius of 0.013 meter.Temperature of air and water are both 313 K.And the temperature of tube wall is 270 K.But some of final value of temperature is much higher than anticipated results(which is physically impossible).However,when I try to siumlate temperature field of damBreak by the solver,it works. I don't know what the problem is.

Julia
Wangzhaoting is offline   Reply With Quote

Old   April 17, 2016, 11:37
Default
  #5
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Hi Julia,

I am happy to take a look at the case. Would you be able to share the case setup files?
liquidspoon is offline   Reply With Quote

Old   April 18, 2016, 02:39
Default
  #6
Member
 
Zhiheng Wang
Join Date: Mar 2016
Posts: 72
Rep Power: 10
Zhiheng Wang is on a distinguished road
Can we solved liquidFilm evaporation cases with this olver ????
Zhiheng Wang is offline   Reply With Quote

Old   April 18, 2016, 06:51
Default
  #7
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Yes. If you check out the latest release, and open the MeshSensitivity branch there is an example case for smooth falling-film evaporation: https://github.com/MahdiNabil/CFD-PC...eltSmoothsEvap .
liquidspoon is offline   Reply With Quote

Old   April 18, 2016, 09:36
Default regarding the source
  #8
Member
 
Zhiheng Wang
Join Date: Mar 2016
Posts: 72
Rep Power: 10
Zhiheng Wang is on a distinguished road
Thanks for the quick reply,
I would lie to impliment code with spices diffusion and combustion but I could not Find any spices diffussion equation specially for two phase reacting solver. ??
Zhiheng Wang is offline   Reply With Quote

Old   April 18, 2016, 11:30
Default
  #9
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Hi Zhiheng. That sounds like an interesting problem. It is a bit beyond the scope of this code. I imagine that you could use this as a starting point, and add a species diffusion equation and appropriate source terms. Best wishes on the project.
liquidspoon is offline   Reply With Quote

Old   April 21, 2016, 10:28
Default
  #10
New Member
 
Julia
Join Date: Jan 2016
Posts: 6
Rep Power: 10
Wangzhaoting is on a distinguished road
Hi,Alex
sorry for my delayed response.Phasechange isn't included in my object,so I modified your solver,and I get rid of the part of phasechange and create a new solver.At first,I try to test the solver by damBreak case,but its temperature fileds is physically impossible.temperature range is from 280K to 313K.Here are the picture and setup files.
Attached Images
File Type: jpg IMG_1.jpg (160.5 KB, 190 views)
Attached Files
File Type: zip dambreak.zip (107.5 KB, 37 views)
Wangzhaoting is offline   Reply With Quote

Old   April 21, 2016, 10:49
Default
  #11
New Member
 
Julia
Join Date: Jan 2016
Posts: 6
Rep Power: 10
Wangzhaoting is on a distinguished road
http://www.wolfdynamics.com/images/c...erTempFoam.pdf
I also used the wolfdynamics page to make modifications to incorporate temperature into OpenFOAM (3.0).The damBreak case works well and the temperature field is physically.When I try to simulate the object which is a simulation of slug flow temperature field about 2.6-meters long pipeline with radius of 0.013 meter.Temperature of air and water are both 313 K.And the temperature of tube wall is 277 K.But some of final value of temperature is much higher than anticipated results(which is physically impossible).Some values are physical. The damBreak works well,but the pipeline can't. I really don't know what the problem is. Here are my setup files.(I draw a grid by Gambit,so ployMesh is big,I didn't upload it)
Attached Files
File Type: zip repipeline.zip (20.0 KB, 45 views)
Wangzhaoting is offline   Reply With Quote

Old   April 22, 2016, 09:45
Default
  #12
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Hi Julia,

I ran your case yesterday and saw the same issue. I really appreciate you bringing this to our attention. We hadn't tried a case where there was such a rapid mixing of liquid and gas of different temperatures before. This turns out to be a really good validation case. I plan to dig into this next week.

I did two things to fix the issue. First, I used gass upwind divSchemes for the velocity, alpha1, and enthalpy fields in fvSchemes. I used gauss interfaceCompression for phirb.

Number 2: interFoam uses a compressive velocity field to counteract numerical diffusion of the interface. This modified velocity field is not conservative! As a result, calculation of temperature from energy in cells diverges with a high cAlpha in this case. The temperature field stays correctly bounded if cAlpha is set to 0 in fvSolution (see attached images with cAlpha = 0 and cAlpha = 1).

I hope to put together a solution that will be more robust soon, but hopefully this points you in the right direction.

cAlpha = 0 CAlpha=0.PNG
cAlpha = 1 CAlpha=1.PNG
liquidspoon is offline   Reply With Quote

Old   April 22, 2016, 09:55
Default
  #13
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 29
akidess will become famous soon enough
In both of your cases your interface is highly smeared though, could that be impacting the solutions?

On a side note, the equations and the code in the wolfdynamics PDF are not consistent.
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   April 22, 2016, 10:40
Default
  #14
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Yes, I agree that the interface is highly smeared in this case. It is good that you point this out. The goal here was to try to track down the source of the inconsistency leading to diverging temperature fields.

The smeared interface would definitely affect the solution. The case Wangzhaoting is looking into isn't really physical to begin with, so it is hard to say what a correct solution really means. It is a 2D dambreak problem. The lengthscales and velocities are such that we would expect turbulence and sprays (not accounted for here).

I'm not familiar with the wolfdynamics pdf. This code is based on our JHT paper: http://heattransfer.asmedigitalcolle...icleid=1829850 (with some modifications since then).
liquidspoon is offline   Reply With Quote

Old   April 28, 2016, 14:17
Default
  #15
Member
 
Atul Kumar
Join Date: Dec 2015
Location: National Centre for Combustion Research and Development
Posts: 48
Rep Power: 10
atulkjoy is on a distinguished road
Hi Alex,
Its difficult to understand your code , can you send commentd file for main source code if possible, i have few questions
1. what equation you have used in heat transfer.
2. Is there any radiation model implemented.
atulkjoy is offline   Reply With Quote

Old   April 28, 2016, 14:33
Default
  #16
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Hi atulkjoy,

Thank you for the message. We have added a pretty good amount of comments to the code, although the comments are pretty sparse in the sections from the original interFoam code. If there are particular sections you are finding unclear, we are happy to clarify those.

The thermal energy transport equation in EEqn.H is:

Code:
label nEnergyLoops(readLabel(pimple.dict().lookup("nEnergyLoops")));
for (int EEqnCount=0; EEqnCount < nEnergyLoops; EEqnCount++)
{

    //Form and solve the energy equation
    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, H)
        + fvm::div(rhoPhi, H)
        - fvc::laplacian(kEff, T)
        - ChillaxFac*( fvm::laplacian(alphaEffRho, H) - fvc::laplacian(alphaEffRho, H) )
        + phaseChangeModel->Q_pc()
    );
    EEqn.solve();
    
    //Now reevaluate T for the updated enthalpy fields
    T = T_0 + rho*H/( limAlpha1*rho1*cp1 + (1-limAlpha1)*rho2*cp2 );
}
Here, we are tracking thermal energy transport, neglecting pressure work and viscous dissipation. H is enthalpy and T is temperature. Heat is diffused by the laplacian of effective thermal conductivity and temperature (line 3). The fourth line with the relaxation factor applies numerical diffusivity to stabilize the thermal energy matrix equation. The last line applies the phase change heat source from the user-specified thermalPhaseChangeModel. We apply solution of this equation and update of the temperature field over multiple passes (kind of like the alternating correction of pressure and velocity in PIMPLE). Please let me know if you have any specific questions about this section.

Right now we don't account for radiation heat transfer. It would certainly be something interesting to explore in the future.
liquidspoon is offline   Reply With Quote

Old   April 29, 2016, 01:30
Default
  #17
Member
 
Zhiheng Wang
Join Date: Mar 2016
Posts: 72
Rep Power: 10
Zhiheng Wang is on a distinguished road
Hi Alex,
I am facing problem to write groovy BC while the Boundary Condition is coupled of 0/T and 0/U , can I write a coupled or mapped boundary.I have tried by Code stream but it was a disaster.
My velocity eqyuations is
U=(diff*(gradYi)/(1-Yi))


and
Temperature equation is

gradT= rho*U*Hf/kappa .
I want to write the expression at inlet using groovy BC.
Zhiheng Wang is offline   Reply With Quote

Old   April 29, 2016, 08:12
Default
  #18
Member
 
Alex
Join Date: Jun 2011
Posts: 33
Rep Power: 14
liquidspoon is on a distinguished road
Hi Zhiheng,

I am not sure I understand the case you are trying to setup. Could you please describe the problem you are modeling? I would be happy to suggest some possible boundary condition sets.
liquidspoon is offline   Reply With Quote

Old   May 1, 2016, 05:12
Default
  #19
Member
 
Zhiheng Wang
Join Date: Mar 2016
Posts: 72
Rep Power: 10
Zhiheng Wang is on a distinguished road
Hi Alex,
At frist thanks a lot for giving your time
I am writting an interface condition for inlet where liquid converts to gas Method is very simple
1. first i have assume inlet temp 10k below boiling temp say Ts=290 and Tb=300
2.Then calculation XFs=exp((-hfg*Mw/Ru)((1/Ts)-(1/Tb)) Basic Clausious Claeperon hfg latent heat , Ru gas constant, Mw molecular weight of liquid.
3.Using Xfs we can find Yfs=(Xfs*Mw)/Mw,mix , Mw,mix is molecular weight of liquid
4.By using Yfs we can calculate velocity Us=diff*(gradYf)/(1-Yf).
5.We check that our assumed Ts value is right or not using gradTs = rho*Us*hfg/k.
If assumed Ts is write, Ts will be the boundary condition temp of inlet else it will go to re-iterate. until Ts=0.98Tb.
i can also send you code of c++ i tried and code of codestram can you plz help regarding this. Can we use for loop in groovy BC???
Attached Files
File Type: c my_boundary best.c (9.3 KB, 38 views)
Kummi likes this.
Zhiheng Wang is offline   Reply With Quote

Old   June 19, 2016, 11:21
Default compile the code
  #20
New Member
 
ali
Join Date: Jun 2016
Posts: 10
Rep Power: 9
abdoli is on a distinguished road
i want to compile the code CFD-PC at https://github.com/MahdiNabil/CFD-PC in foam-extend-3.1. but i received errors that are shown in the file. could anyone help me.
Attached Files
File Type: docx geometriczerofield.docx (16.2 KB, 39 views)
raj kumar saini likes this.
abdoli is offline   Reply With Quote

Reply

Tags
interfoam, phase change

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


Similar Threads
Thread Thread Starter Forum Replies Last Post
Effective heat capacity - Phase change does not work!!! papteo CFX 8 October 31, 2013 05:15
Constant velocity of the material Sas CFX 15 July 13, 2010 08:56
Which Code is Best for Two Phase, Heat Transfer earlc Main CFD Forum 0 September 23, 2005 20:26
compressible two phase flow in CFX4.4 youngan CFX 0 July 1, 2003 23:32
help needed about phase change Yanhu Guo Main CFD Forum 4 January 23, 2001 23:16


All times are GMT -4. The time now is 21:44.