[liquidspoon]

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.

[Wangzhaoting]

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!

[liquidspoon]

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.

[Wangzhaoting]

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

[liquidspoon]

Hi Julia,

I am happy to take a look at the case. Would you be able to share the case setup files?

[Zhiheng Wang]

Can we solved liquidFilm evaporation cases with this olver ????

[liquidspoon]

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 .

[Zhiheng Wang]

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. ??

[liquidspoon]

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.

[Wangzhaoting]

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.

[Wangzhaoting]

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)

[liquidspoon]

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

[akidess]

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.

[liquidspoon]

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).

[atulkjoy]

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.

[liquidspoon]

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:

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.

[Zhiheng Wang]

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.

[liquidspoon]

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.

[Zhiheng Wang]

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???

[abdoli]

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.