[Lorenzo210]

Hello everyone,
I've been enjoying openFOAM for more than a year, and I'm still learning so much that I consider myself a very new user.

I'd like to share with you a topic about the solver interCondensatingEvaporatingFoam. I am using openFOAM v2006, and I am trying to model an evaporating two-phase flow inside an horizontal microchannel.

Starting off with the mesh, I used blockMesh to build a 2D axis-symmetric geometry. I can use this because the inner diametre is very small, and even the effect of gravity forces are neglected.
The patches are: topWall (wall), axis (empty) front and back (wedge), inlet and outlet (patch).
I am simulating an annular flow, with liquid near the wall and vapour through the middle of the channel. To do this, I've split the block into two.
I have already a fully-developed pattern at the inlet due to an adiabatic numerical simulation.

Let's talk about the boundary conditions:
I fixed the inlet velocities, the outlet pressures; I use the k-omega SST turbulence model because the vapour Reynolds numbers are very high.
I fixed the inlet temperature of both phases (equal to saturation temperature), and then I used the type externalWallHeatFluxTemperature for the wall, because I'd like to have a fixed heat flux (20000 [kW/m2]).

And finally here comes the problem:
I tried using the Lee evaporation model to implement the phase change model inside the simulation. The problem with this one is that the phase change process occurs uniformly inside the liquid zone, in particular near the hot wall; this seems obvious beacuse the process is activated by a temperature difference ( so it seems normal that the liquid starts evaporating from the hottest region), but I need the condition that the evaporation process is located only near the interface liquid-vapour.

I am recently using the latest version of OpenFOAM v2006, and I found that there is another phase change model available: interfaceHeatResistance ( https://www.openfoam.com/documentati...8H_source.html). So I tried writing the following text in PhaseChangeProperties file, inside the constant directory:

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object phaseChangeProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

phaseChangeTwoPhaseModel interfaceHeatResistance;

R 0.01;
Tactivate 278.98;
spread 3;



// ************************************************** *********************** //

R value is an initial guess, while Tactivate is equal to the saturation temperature. I need to define "spread", which value is explained in the documentation I link just above.

I am running openFOAM so that I can find whether something should be added, but I found this error and I can't figure out what the error is about:

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 295eef47-20201012 OPENFOAM=2006 patch=201012
Arch : "LSB;label=32;scalar=64"
Exec : interCondensatingEvaporatingFoam
Date : Dec 09 2020
Time : 11:17:38
PID : 11291
I/O : uncollated
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


PIMPLE: Operating solver in PISO mode

Reading field p_rgh

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
Creating temperaturePhaseChangeTwoPhaseMixture

Selecting phaseChange model interfaceHeatResistance
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in /lib/x86_64-linux-gnu/libpthread.so.0
#3 void Foam::divide<double, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensioned<double> const&) in /usr/lib/openfoam/openfoam2006/platforms/linux64GccDPInt32Opt/bin/interCondensatingEvaporatingFoam
#4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:perator/<double, Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::dimensioned<double> const&) in /usr/lib/openfoam/openfoam2006/platforms/linux64GccDPInt32Opt/bin/interCondensatingEvaporatingFoam
#5 Foam::temperaturePhaseChangeTwoPhaseMixtures::inte rfaceHeatResistance::correct() at ??:?
#6 Foam::temperaturePhaseChangeTwoPhaseMixtures::inte rfaceHeatResistance::interfaceHeatResistance(Foam: :thermoIncompressibleTwoPhaseMixture const&, Foam::fvMesh const&) at ??:?
#7 Foam::temperaturePhaseChangeTwoPhaseMixture::addco mponentsConstructorToTable<Foam::temperaturePhaseC hangeTwoPhaseMixtures::interfaceHeatResistance>::N ew(Foam::thermoIncompressibleTwoPhaseMixture const&, Foam::fvMesh const&) at ??:?
#8 Foam::temperaturePhaseChangeTwoPhaseMixture::New(F oam::thermoIncompressibleTwoPhaseMixture const&, Foam::fvMesh const&) at ??:?
#9 ? in /usr/lib/openfoam/openfoam2006/platforms/linux64GccDPInt32Opt/bin/interCondensatingEvaporatingFoam
#10 __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#11 ? in /usr/lib/openfoam/openfoam2006/platforms/linux64GccDPInt32Opt/bin/interCondensatingEvaporatingFoam
Floating point exception (core dumped)

I know that a Floating point exception (core dumped) error occurs when maybe there is some bad operation, like for example dividing to zero, but I really can't figure out what to change or what I am doing wrong. I tried to read all the documentation, trying to undestand how it works but I am not really confident with C++ so maybe I'm just not keeping attention on right things.

I hope someone may help me with that issue, or just give some information in order to better understand the solver and all the applications.

Thank you very much for the attention.

Lorenzo

PS: I tried to use also icoReactingMultiphaseInterFoam, which can use the Hertz-Knudsen phase change model "kineticGasEvaporation" in order to localize the evaporation process where alpha.liquid is between 0.01-0.99 (for istance), so that liquid becomes vapour only at the interface. It works, but I had other issues linked to the fact that I am not able to use dynamic mesh refinement with that solver ( I'm not expert, and it seems too hard to modify the solver for that), so I believe that interCondensatingEvaporatingFoam could be more useful to me.

Thank you again.

EDIT:
I found a sort of tutorial case, which is the stefanProblem validation case located in OpenFOAM/OpenFOAM-v2006/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam. I copied all the features for phaseChangeProperties etc. but when I run the application I find the same error.

[Lorenzo210]

This is the phaseChangeProperties file:

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object phaseChangeProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

phaseChangeTwoPhaseModel interfaceHeatResistance;//constant;


R 1e6;
maxAlphaRate 1;
spread 3;


coeffC 0;
coeffE 500;

// ************************************************** *********************** //

Here there is the controlDict file:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application interCondensatingEvaporatingFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 1;

deltaT 1e-5;

writeControl adjustable;

writeInterval 0.001;

purgeWrite 0;

writeFormat ascii;

writePrecision 12;

writeCompression off;

timeFormat general;

timePrecision 10;

runTimeModifiable yes;

adjustTimeStep yes;

maxCo 0.01;

maxAlphaCo 0.01;

maxDeltaT 0.01;

/*interfaceHeight1
{
// Mandatory entries
type interfaceHeight;
libs (fieldFunctionObjects);
locations ((0 0.0001 1e-5));

// Optional entries
alpha alpha.liquid;
direction (1 0 0);
liquid true;
interpolationScheme cellPoint;

// Optional (inherited) entries
writePrecision 16;
writeToFile true;
useUserTime true;

region region0;
enabled true;
log true;
timeStart 0;
timeEnd 1000;
executeControl timeStep;
executeInterval 1;
writeControl timeStep;
writeInterval 3;
}

*/
// ************************************************** *********************** //

I changed a little bit, cutting off the functionObject of the tutorial case.

Maybe the error is due to the fact that initially there is not a clear interface, but really I can'figure out why it's not working.

[Lorenzo210]

I finally found the error, I'll explain the way I solve it:
I was looking at some discussions in this forum, and I learnt the way to recognise the nature of the problem: the line #3 begins with "Foam::divide", it means that somewhere in the code there is a division per zero, or something like that. In line #5 we can see where the problem occurs: it is in Foam::temperaturePhaseChangeTwoPhaseMixtures::inte rfaceHeatResistance::correct() , so that means the case has some problem regarding the thermophysical properties.
Then, I compared my case with the tutorial case stefanProblem, and I found that in transportProperties file I put the same value of hf (heat of fusion) for both liquid and vapour species. This probably causes the problem with the correction made by Foam::temperaturePhaseChangeTwoPhaseMixtures::inte rfaceHeatResistance.

Finally, I set hf liquid to zero, and ran successfully the case.

Edit: analyzing interfaceHeatResistance.C I found that the latent heat L used for the calculation of m_dot_evap and m_dot_cond is: L= h_f2 - h_f1. Therefore, this is why the solver gives an errore when h_f2 = h_f1.. there is a division per zero (L).

[rocco96]

Hi Lorenzo,
I have a question regarding interCondensatingEvaporatingFoam.
I'm simulating an Open Channel (so in my domain there are air and water) with phase change. I want that the phase change to happen only at interface. I'm using a phaseChangeTwoPhaseModel constant, but the results don't convince me.
Do you think i should use interfaceHeatResistance model?
Could you also please recommend me some papers explaining how to choose the CoeffC and CoeffE?




Best regards,
Rocco

[Lorenzo210]

Hello Rocco,

Nice to see your interest in this solver.

I think the interfaceHeatResistance phase change model fits better with your application. This is based on Hardt-Wondra phase change model, which you can find described in the article "Evaporation model for interfacial flows based on a continuum-field representation of the source terms" by S. Hardt and F. Wondra.
The Hardt-Wondra model is based on kinetic gas theory, where mass transfer and heat transfer occur at the interface. In particular, they used a "spreading" factor so that numerical instabilities are generally decreased.
There are many articles you can read about numerical modelling of two-phase flow, evaporation/condensation and so on.

I can suggest you a review by Z. Guo, D.F. Fletcher and B.S. Haynes "A Review of Computational Modelling of Flow Boiling in Microchannels".

For what concerns interfaceHeatResistance coefficients, I usually use spread = 1 and maxAlphaRate = 1; the R coefficient is calculated via thermophysical properties of the liquid-vapour mixture (water and vapour in your case).

R = ((2*sigma_gg)/(2-sigma_gg)) * ( (rho_gas*h_fg^2)/(2 * Pi * T^3 * MW)^0.5 )

Where:

sigma_gg = accomodation coefficient [/] (varies from 0.001 to 1)
rho_gas = vapour density [kg/m^3]
h_fg = Latent heat of vaporisation [J/kg]
T = Saturation temperature [K]
MW = Molecular weight [kg/kmol]

I hope I wrote it correctly, please control the paper so you can be sure to put the correct value in your simulation.

You can look for articles where Schrage Model or Tanasawa model are implemented to learn more about what has been done in literature.

Here you find an article where the evaporation model and the meaning of sigma_gg are explained: " Numerical investigation of hydrodynamics and heat transfer of elongated bubbles during flow boiling in a microchannel " by M. Magnini, B. Pulvirenti and J. R. Thome.

Apart from these tips, I suggest you to pay attention at the boundary conditions, discretizations schemes and tolerance.

Wish this is helpful, and that you will find good results.
Let me know how it proceeds.
Best regards,

Lorenzo

PS: Reading again your message I saw you also asked about CoeffC and CoeffE values. Let me do a little recap: these are used exclusively for constant model. Their values are semi-empirical, which means you can evaluate them by iteration.
In other words, you 1) make a first simulation with initial guessed values 2) See the results, if the interface Temperature differs too much from Saturation temperature, you need to take a higher coefficient; if your simulation runs too slow and hardly goes on, you need a smaller CoeffE or CoeffC. 3) repeat the simulation.
Giving a rough explanation, CoeffC and CoeffE express the time that liquid takes to evaporate (or vapour takes to condense). If the CoeffE is high, the evaporation occurs quickly, which could cause numerical instabilities; otherwise, if evaporation runs slowly some energy would overheat the interface.

Please refer to literature to better understand CoeffC and CoeffE meaning: in this article "Numerical Simulation of Laminar Liquid Film Condensation in a Horizontal Circular Minichannel" by E. Da Riva and D. Del Col they used a very high coefficient, and they mention a work where low value are used. It depends on the situation simulated.

[rocco96]

Hi Lorenzo,
thank you very much for the advice and explanation, it was all very useful to me. In the end I chose to use the constant model, as by setting a saturation temperature close to that of the interface, the phase change occurs right at the interface. In particular, I added the equation of the transport of a scalar (humidity), and I created a source term for the humidity so that the evaporated and condensed mass also modified the variation of humidity. I'm finally ready to launch the simulation, but I still have two doubts that I hope you can solve.
1) How do you advise me to set the latent heat for the two phases? For now I have set only the air side, and the water side I left it at 0. I am considering the steam with the thermophysical properties of air, as my domain must represent a flat channel containing water and air.
2) My case takes into account a fixed wind source, therefore a speed that is set via fvOption which represents a source term for the momentum equation. In the solver folder, however, there is no Ueqn equation, ie the one that solves the velocity field. This equation is taken by calling the solver interPhaseChangeFoam in the option file?
For now I have launched a simulation, but the deltaT becomes of the order of 10 ^ -6 to guarantee a courant of 0.5, the deltaT I have set is 0.02. The problem I think is precisely when calculating p_rgh, as the number of iterations to calculate a single value of p_rgh is 406.
Sorry for the many questions, but I'm having a hard time finding detailed documentation for this solver.
Thank you again for your availability and for your reply,
Sincerely, Rocco

[Lorenzo210]

Hi Rocco,
I'm glad to help you, and I wish you can succeed in you work.

I'll give you some information based on my experience, hoping that it will make this solver easier to understand. Sorry if this reply will be a little longer.
1) As you can see in in the phase change model "constant.C", the latent heat is calculated as " L = mixture_.Hf2() - mixture_.Hf1() " (line 204).
When you write down the properties of liquid and vapour in "transportProperties" you have to define hf of each phase (for example, in condensatingVessel liquid is phase1 and vapour is phase2).
Then, in order to have positive L values you should set hf_vapour equals to the latent heat of vaporisation, and hf_liquid to zero (like you already did).

2) Yeah, the files you cannot file into the solver folder are taken from other solvers. interCondensating is based on interPhaseChangeFoam, so they have a lot in common (for example also alphaEqn.H is in common). You can see that also interFoam and VoF folders are attached in "options", so if there's something you don't find in interPhaseChangeFoam it is probably located in interFoam or in VoF.

It's interesting the way you apply the wind source term, let me know if this works. I'm struggling to fully understand how mass source terms are implemented in alphaEqn.H and how MULES works in it. It's kinda different from interFoam, and I want to be really sure that all the things are set in a proper way.

For the last part of your message, it's a bit hard to say something without more precise informations about the settings and the conditions of your simulations.
Starting with the time step, the deltaT depends on your mesh size, on the velocity field and on the Courant number you can set in controlDict file. So basically you cannot set deltaT and Courant at the same time. Without any other info, I can't say if the time step is strange or not, for istance my simulations run with orders of 10^-7 s. For your first simulations, you can try a rough geometry and see how it behaves, and then keep on refining the mesh or reducing Courant number in order to obtain more accurate results (this is pretty much what I've done with my simulations).

Let's focus on p_rgh: the resolution of the PIMPLE loop shouldn't affect the time step, which is stricly linked to the CFL condition (see https://www.cfd-online.com/Wiki/Cour...Lewy_condition ). However, a large number of iterations can increase the computational timing cost. This means that it takes longer to calculate a "deltaT" timestep of simulation. 400 iterations for pressure are quite a lot in my opinion, but I can't say for sure what's the problem if there is any.
If you can share some output from your calculations it will help understanding what's going on. It would be useful to know also the fvSolution, fvSchemes and boundary conditions.
Which solver are you using for p_rgh? And what is the tolerance? These two are generally the main things to take care in my opinion, and also the usage of the right boundary conditions.

I look forward to hearing from you,

Lorenzo

[rocco96]

Hi Lorenzo, thanks for the immediate reply and for the clarifications.
As for the deltaT, I certainly explained myself badly. In my simulation I fix the courant, so it is obvious that the deltaT is variable. Such a low deltaT worries me because the same simulation was launched with the solver "interFoam", which I modified in "interTempFoam" in order to solve the thermal balance equation and therefore to solve the temperature field. The case is identical (obviously evaporation and condensation are not calculated here), therefore the same mesh, boundary conditions, etc. (obviously without some thermophysical properties present in the new interCondensatingEvaporatingFoam solver). So summarizing with interTempFoam my simulation has a deltaT of about 0.02, while with interCondensatingEvaporatingFoam modified as I showed you in the previous post, it has a deltaT of the order of 10 ^ -6, with the same Courant number equal to 0.5. To tell the truth after a day of simulation the situation seems to have improved, as the deltaT is about 0.0025.
However, I leave you attached the files relating to the controlDict, fvSchemes, fvSolution and a log file of the simulation. The boundary conditions instead are these:
1) U top = slip, bottom = noSlip, sides = cyclic, inlet / outlet = cyclic
2) p_rgh top = Zero Gradient, bottom = Zero Gradient, sides = cyclic, inlet / outlet = cyclic
3) alpha top = Zero Gradient, bottom = Zero Gradient, sides = cyclic, inlet / outlet = cyclic
4) T top = Zero Gradient, bottom = Zero Gradient, sides = cyclic, inlet / outlet = cyclic
5) q (humidity) top = Zero Gradient, bottom = ZeroGradient, sides = cyclic, inlet / outlet = cyclic.
Finally, the mesh is made up of approximately 10 million cells. I'm curious to see if what I'm doing makes sense, because not having a very powerful PC, I risk waiting weeks for something that doesn't actually work.
Thanks a lot again for your availability.
Yours sincerely,
Rocco.

[Lorenzo210]

Hello Rocco,
Ok I see what you mean. As far as I know, phase change affects the time step; the faster vapour is generated, the lower the time step will be. This is linked to the fact that convergence is harder to reach with higher CoeffE/CoeffC.
It's not surprising to me that delta T is reduced.

Boundary conditions: Help me undestand better your case: the geometry is an open channel 3D, where the bottom is the ground/wall channel/whatever, the top is a patch where outside there is atmosphere, is this correct?
Furthermore, I never used type cyclic, so I'm not really sure what that stands for. It means that the sides patch's variables are connected each other, doesn't it?
Said that, I think your BC are correct. For p_rgh I would use dirrefent BCs, for example for walls (bottom?) I would use fixedFluxPressure, and for top (atmosphere?) I would use totalPressure. Please refer to damBreak tutorial for details about these BCs.

10 Million cells are quite a lot!!

Let me get back to the schemes:
fvSchemes:

Never used crankNicolson, always Euler. So I'm sorry but I can't help with it.
I can say one thing about it: comparing alphaEqn.H of interFoam and interCondensatingEvaporatingFoam I noticed that in interFoam there is a specific calculation of volumetric flux for CrankNicolson option (phiCN), whereas interCondensating doesn't have it.
I'm recently learning these things, so actually I don't know much about it. I just want to warn you about that.

divSchemes are quite the same I use, I think they're okay.
Nothing to say about the others.

fvSolution:

In my experience, GAMG solver gives faster resolution than PCG. I'll give a try with:

"pcorr.*"
{
solver PCG;
preconditioner
{
preconditioner GAMG;

smoother GaussSiedel;

tolerance 1e-5;
relTol 0;
}
}

p_rgh
{
solver GAMG;
preconditioner DIC;
smoother DIC;
tolerance 1e-07;
relTol 0.01;
}

Last thing on the Courant number, usually courant number should be lower. This is because you can get more accurate results, and limiting time step can also help with convergence problems.

To sum up, I think you can first try changing p.corr and p_rgh solvers, and see what happens in first time steps. Maybe you can change controlDict writeInterval so you can see if something is wrong with the solution without waiting too much.
Then I would also try with BC, depending on what happened in the first modification. This is because I think it's not a big trouble having zeroGradient for p_rgh BC.

Hope this helps, let me know what happens.
Cheers,
Lorenzo

[ajaymanjunath]

Hi Lorenzo,

I am new to OpenFOAM and I am trying to simulate a droplet evaporation case. Can I use interCondensatingevaporating Foam to simulate the case??

[rocco96]

Hi Lorenzo,
Thank you very much for your help. I ran the simulation with your advices, so using GAMG and the velicity of simulation has improved!
I should ask you another question, because this is my eighth month in which I'm working on this solver. I modified the contant.C file in order to limit the evaporation only in the interface. The way in which I make this, is this: "

const volScalarField Vcoeff
(
coeffE_*mixture_.rho1()*limitedAlpha1*L*pos(T - TSat)*pos(scalar(0.5)-limitedAlpha1)
);


"
For example this is the modify in the temperature source term. The problem is that when there is evaporation, there is a discontinuity in the interface (I see a line with TSat value). So I would ask you, if you know how works the source term of temperature. I think that the source term has the purpose to decrease the temperature (when there is evaporation, because T>Tsat), and to increase the temperature otherwise, in order to make the phase change take place at a constant temperature (Tsat). Have you some advise, or some documentation about this?
Best Regards,
Rocco.

[Pavithra]

Hello Everyone,

I am trying to simulate a heat pipe using interCondensatingEvaporatingFoam. As it is a heat pipe, it has evaporation taking place at the lower end and condensation taking place at the top end.

1) Can interCondensatingEvaporatingFoam handle this ?
2) I am planning to use the constant model. In that case, can I enter same values for CoEffC = CoEffE = r/Tsat.
3) In the transportProperties, can I enter same value of hf_liquid and hf_vapour ?

I already tried these. But, I am facing two issues
1) Evaporation is taking place as expected. But condensation is not seen.
2) The temperature goes beyond the saturation temperature (which is unphysical).

Someone, pls kindly clarify these.

Thanks in advance.

With Thanks,
Pavithra.

[Lorenzo210]

Hi all, I'm sorrry, but I've been quite busy in the last months.

@ajaymanjunath you can try it ;) Personally, I had some problem in runnning the solver with Hardt-Wondra phase change model, so I moved on another customized solver.

@rocco96 Glad to hear the simulation has improved. The discontinuity is probably due to the limiter you used. pos(scalar(0.5)-limitedAlpha1) is equal to 1
in any cell that has alpha1 > 0.5, which is also inside all the liquid. I guess you want the interface to be located in all cells where 0 < alpha1 < 1, so all the cells
with 0<alpha1<0.5 are being left out from the calculation. I would try with pos(alpha1)*pos(1-alpha1), so that you exclude liquid and vapour.
The temperature source term has the purpose you already said: the evaporation consumes energy with constant temperature. Right now I don't remember any useful documentation, but I'll
come back when I find anything useful.

@Pavithra I think you can simulate both phase change processes in a simulation. The solver calculates mass flow rate evaporated and condensated in relation to the equations you find
in Constant.C. In order to have both, you should have CoeffC!=0 and CoeffE!=0.
As I wrote above some time ago, you should look at the code, in line 204 in Contant.C the latent heat is calculated as: mixture_.Hf2() - mixture_.Hf1(), where
2 is phase2 (usually vapour phase) and 1 is phase1 (liquid). If you put the same value for each phase, then L=0, and it gives error.
For what concerns the issues:
1) Please check you boundary conditions, transport properties and phaseCHangeproperties and be sure that everything is set correctly. If you'd like to share your file, we can help you looking
for the issue, you didn't tell enough informations.
2) If CoeffE or CoeffC are too low, the interface temperature will differ from the saturation temperature.

[Pavithra]

Originally Posted by Lorenzo210

Hi all, I'm sorrry, but I've been quite busy in the last months.

@ajaymanjunath you can try it Personally, I had some problem in runnning the solver with Hardt-Wondra phase change model, so I moved on another customized solver.

@rocco96 Glad to hear the simulation has improved. The discontinuity is probably due to the limiter you used. pos(scalar(0.5)-limitedAlpha1) is equal to 1
in any cell that has alpha1 > 0.5, which is also inside all the liquid. I guess you want the interface to be located in all cells where 0 < alpha1 < 1, so all the cells
with 0<alpha1<0.5 are being left out from the calculation. I would try with pos(alpha1)*pos(1-alpha1), so that you exclude liquid and vapour.
The temperature source term has the purpose you already said: the evaporation consumes energy with constant temperature. Right now I don't remember any useful documentation, but I'll
come back when I find anything useful.

@Pavithra I think you can simulate both phase change processes in a simulation. The solver calculates mass flow rate evaporated and condensated in relation to the equations you find
in Constant.C. In order to have both, you should have CoeffC!=0 and CoeffE!=0.
As I wrote above some time ago, you should look at the code, in line 204 in Contant.C the latent heat is calculated as: mixture_.Hf2() - mixture_.Hf1(), where
2 is phase2 (usually vapour phase) and 1 is phase1 (liquid). If you put the same value for each phase, then L=0, and it gives error.
For what concerns the issues:
1) Please check you boundary conditions, transport properties and phaseCHangeproperties and be sure that everything is set correctly. If you'd like to share your file, we can help you looking
for the issue, you didn't tell enough informations.
2) If CoeffE or CoeffC are too low, the interface temperature will differ from the saturation temperature.

Hi Lorenzo,

Thanks for the reply. I was able to get condensation.

The problem is with the choice of CoeffC and CoeffE. Smaller CoeffC leads to nonphysical temperatures. Larger CoeffC leads to numerical instability and requires very very small time steps. I understand that this is a limitation of Lee's model.

I am trying to handle this by implementing Schrage or Energy balance mass transfer models into the framework for interCondensatingEvaporatingFoam.

Thank You.

[SURENDAR KUMAR V]

Originally Posted by Lorenzo210

I finally found the error, I'll explain the way I solve it:
I was looking at some discussions in this forum, and I learnt the way to recognise the nature of the problem: the line #3 begins with "Foam::divide", it means that somewhere in the code there is a division per zero, or something like that. In line #5 we can see where the problem occurs: it is in Foam::temperaturePhaseChangeTwoPhaseMixtures::inte rfaceHeatResistance::correct() , so that means the case has some problem regarding the thermophysical properties.
Then, I compared my case with the tutorial case stefanProblem, and I found that in transportProperties file I put the same value of hf (heat of fusion) for both liquid and vapour species. This probably causes the problem with the correction made by Foam::temperaturePhaseChangeTwoPhaseMixtures::inte rfaceHeatResistance.

Finally, I set hf liquid to zero, and ran successfully the case.

Edit: analyzing interfaceHeatResistance.C I found that the latent heat L used for the calculation of m_dot_evap and m_dot_cond is: L= h_f2 - h_f1. Therefore, this is why the solver gives an errore when h_f2 = h_f1.. there is a division per zero (L).

Hi sir,
I couldnt able to use Lee model in interCondensatingEvaporatingFoam. The solver allows only constant model and InterfaceHeatResistanceModel. May I know why? and I need to know where i need the modifications to use LEE Model in InterCondensatingEvaporatingFoam. I use OpenFoam 2106 version. I am new to openFoam, so question a may feel silly.



Thanks in Advance !!

[Henning86]

FYI: This may be useful



recently Maganini et [1] al published a paper about microchannel that is partially based on the library twophaseflow

https://github.com/DLR-RY/twophaseflow


[1] Magnini, M., Municchi, F., El Mellas, I., & Icardi, M. (2022). Liquid film distribution around long gas bubbles propagating in rectangular capillaries. International Journal of Multiphase Flow, 103939.

[saicharan662000@gmail.com]

Originally Posted by Lorenzo210

Hello everyone,
I've been enjoying openFOAM for more than a year, and I'm still learning so much that I consider myself a very new user.

I'd like to share with you a topic about the solver interCondensatingEvaporatingFoam. I am using openFOAM v2006, and I am trying to model an evaporating two-phase flow inside an horizontal microchannel.

Hi lorenzo,
I am developing the solver similar to yours. I am following the source terms implemented in satos paper. A sharp-interface phase change model for a mass-conservative interface tracking method Yohei Sato , Bojan Niceno.
My mass transfer source term is like this

Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::myPhaseChange::mDotAlphal() const
{
const volScalarField limalpha1(min(max((alpha1()), scalar(0)), scalar(1)));
const volVectorField& U = db().lookupObject<volVectorField>("U");
const fvMesh & mesh = U.mesh();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
const dimensionedScalar secondValue("secondValue",dimensionSet(1,-2,-1,0,0,0,0),0.0 );
const volVectorField gradAlpha(fvc::grad(alpha1()));
const dimensionedScalar deltaN = 1e-8/pow(average(alpha1().mesh().V()), 1.0/3.0);
volVectorField nHatfv(gradAlpha/(mag(gradAlpha) + deltaN));
const volScalarField Tgrad = fvc::grad(T) & nHatfv;
Info << "grad : " << max(fvc::grad(T)).value() << " m" <<endl;
Info << "normal : " << max(nHatfv).value() << " m" <<endl;

return Pair<tmp<volScalarField>>
    ( 
      (limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * (mag(Tgrad))/h_lv_,limalpha1 * 0.0 * secondValue
    );
}

affter few iterations my solution is not converging. The value of grad(T) is increasing to 1e10 and so on.

grad : (3.76772e+82 0 0) m
normal : (1 0 0) m
mass transfer: 4.98472e+73 kg/m2/s
Max interFaceArea: 0.02 m2
Courant Number mean: 0.0807874 max: 0.95333
Interface Courant Number mean: 0.000176765 max: 0.0456387
deltaT = 2.00643e-79
Time = 891.993

grad : (3.76772e+82 0 0) m
normal : (1 0 0) m
MULES: Solving for alpha.water
Liquid phase volume fraction = 7.44071e+08  Min(alpha.water) = -3.1905e-06  Max(alpha.water) = 3.72035e+11
grad : (3.76772e+82 0 0) m
normal : (1 0 0) m
MULES: Solving for alpha.water
Liquid phase volume fraction = 8.43143e+08  Min(alpha.water) = -2.8265e-06  Max(alpha.water) = 4.21571e+11
grad : (3.76772e+82 0 0) m
normal : (1 0 0) m
DICPCG:  Solving for p_rgh, Initial residual = 0.805597, Final residual = 8.2184e-13, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.56469, Final residual = 1.41949e-13, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.615492, Final residual = 3.83556e-13, No Iterations 2
grad : (3.76772e+82 0 0) m
normal : (1 0 0) m
DICPCG:  Solving for p_rgh, Initial residual = 0.799576, Final residual = 6.62343e-09, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.447132, Final residual = 1.45708e-13, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.529525, Final residual = 1.00829e-13, No Iterations 2
grad : (3.76772e+82 0 0) m
normal : (1 0 0) m
DICPCG:  Solving for p_rgh, Initial residual = 0.549217, Final residual = 1.27283e-13, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.603145, Final residual = 5.07359e-14, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.628802, Final residual = 4.96658e-12, No Iterations 2
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PBiCGStab::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#6  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/hari/OpenFOAM/hari-8/platforms/linux64GccDPInt32Opt/bin/interPhaseChangeFoam_12"
#7  Foam::fvMatrix<double>::solve() in "/home/hari/OpenFOAM/hari-8/platforms/linux64GccDPInt32Opt/bin/interPhaseChangeFoam_12"
#8  ? in "/home/hari/OpenFOAM/hari-8/platforms/linux64GccDPInt32Opt/bin/interPhaseChangeFoam_12"
#9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10  ? in "/home/hari/OpenFOAM/hari-8/platforms/linux64GccDPInt32Opt/bin/interPhaseChangeFoam_12"
Floating point exception (core dumped)

Do you have any idea about this problem?
Thanks in advance

[saicharan662000@gmail.com]

Hey lorenzo,
Can you attach your case file?
Thanks in advance

[brunobw]

Dear Fomers,
I'm trying to model the evaporation of a small pool of ethanol on a room submited to natural convection. Wondering to track the cloud to assess its concentration and dispersion over the room.
Read some nice works about gas cloud dispersion using fireFoam with combustion disabled.
Wonder if some new solver are more suitable for my case.
Any directions would be very nice.
Thanks.