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

Mass transfer model in interCondensatingEvaporatingFoam

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

Like Tree9Likes
  • 3 Post By Lorenzo210
  • 3 Post By Lorenzo210
  • 1 Post By Lorenzo210
  • 1 Post By Lorenzo210
  • 1 Post By Lorenzo210

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   December 9, 2020, 09:06
Default Mass transfer model in interCondensatingEvaporatingFoam
  #1
Member
 
Lorenzo
Join Date: Apr 2020
Location: Italy
Posts: 39
Rep Power: 6
Lorenzo210 is on a distinguished road
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.

Last edited by Lorenzo210; December 11, 2020 at 08:44.
Lorenzo210 is offline   Reply With Quote

Old   December 11, 2020, 08:50
Default
  #2
Member
 
Lorenzo
Join Date: Apr 2020
Location: Italy
Posts: 39
Rep Power: 6
Lorenzo210 is on a distinguished road
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 is offline   Reply With Quote

Old   December 11, 2020, 11:36
Default
  #3
Member
 
Lorenzo
Join Date: Apr 2020
Location: Italy
Posts: 39
Rep Power: 6
Lorenzo210 is on a distinguished road
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).
TeresaT, Pavithra and mani4be like this.

Last edited by Lorenzo210; December 17, 2020 at 10:51.
Lorenzo210 is offline   Reply With Quote

Old   March 29, 2021, 10:28
Default
  #4
New Member
 
rocco
Join Date: Nov 2020
Posts: 13
Rep Power: 5
rocco96 is on a distinguished road
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
rocco96 is offline   Reply With Quote

Old   March 29, 2021, 15:52
Default
  #5
Member
 
Lorenzo
Join Date: Apr 2020
Location: Italy
Posts: 39
Rep Power: 6
Lorenzo210 is on a distinguished road
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.
Pavithra, ZZW and Roham..Seif like this.
Lorenzo210 is offline   Reply With Quote

Old   April 26, 2021, 10:34
Default
  #6
New Member
 
rocco
Join Date: Nov 2020
Posts: 13
Rep Power: 5
rocco96 is on a distinguished road
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
rocco96 is offline   Reply With Quote

Old   April 26, 2021, 17:28
Default
  #7
Member
 
Lorenzo
Join Date: Apr 2020
Location: Italy
Posts: 39
Rep Power: 6
Lorenzo210 is on a distinguished road
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
Pavithra likes this.
Lorenzo210 is offline   Reply With Quote

Old   April 27, 2021, 10:30
Default
  #8
New Member
 
rocco
Join Date: Nov 2020
Posts: 13
Rep Power: 5
rocco96 is on a distinguished road
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.
Attached Files
File Type: txt fvSchemes.txt (2.0 KB, 16 views)
File Type: txt controlDict.txt (1.6 KB, 15 views)
File Type: txt fvSolution.txt (2.3 KB, 15 views)
rocco96 is offline   Reply With Quote

Old   May 1, 2021, 06:59
Default
  #9
Member
 
Lorenzo
Join Date: Apr 2020
Location: Italy
Posts: 39
Rep Power: 6
Lorenzo210 is on a distinguished road
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
Pavithra likes this.
Lorenzo210 is offline   Reply With Quote

Old   May 25, 2021, 10:39
Default
  #10
New Member
 
Ajay
Join Date: May 2021
Location: India
Posts: 1
Rep Power: 0
ajaymanjunath is on a distinguished road
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??
ajaymanjunath is offline   Reply With Quote

Old   June 12, 2021, 07:31
Default
  #11
New Member
 
rocco
Join Date: Nov 2020
Posts: 13
Rep Power: 5
rocco96 is on a distinguished road
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.
rocco96 is offline   Reply With Quote

Old   September 8, 2021, 22:42
Default interCondensatingFoam for Heat Pipe
  #12
Member
 
Join Date: Apr 2019
Location: India
Posts: 81
Rep Power: 7
Pavithra is on a distinguished road
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.
Pavithra is offline   Reply With Quote

Old   September 14, 2021, 15:52
Default
  #13
Member
 
Lorenzo
Join Date: Apr 2020
Location: Italy
Posts: 39
Rep Power: 6
Lorenzo210 is on a distinguished road
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 likes this.
Lorenzo210 is offline   Reply With Quote

Old   September 15, 2021, 07:17
Default
  #14
Member
 
Join Date: Apr 2019
Location: India
Posts: 81
Rep Power: 7
Pavithra is on a distinguished road
Quote:
Originally Posted by Lorenzo210 View Post
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.
Pavithra is offline   Reply With Quote

Old   November 30, 2021, 20:56
Default
  #15
New Member
 
SURENDAR KUMAR V
Join Date: Nov 2021
Posts: 5
Rep Power: 4
SURENDAR KUMAR V is on a distinguished road
Quote:
Originally Posted by Lorenzo210 View Post
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 !!

Last edited by SURENDAR KUMAR V; December 7, 2021 at 02:46.
SURENDAR KUMAR V is offline   Reply With Quote

Old   January 17, 2022, 04:30
Default Micro Changel
  #16
Member
 
Join Date: Nov 2012
Posts: 83
Rep Power: 13
Henning86 is on a distinguished road
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.
Henning86 is offline   Reply With Quote

Old   September 6, 2022, 01:12
Default solution is not converging
  #17
Member
 
hari charan
Join Date: Sep 2021
Location: India,hyderabad
Posts: 96
Rep Power: 4
saicharan662000@gmail.com is on a distinguished road
Quote:
Originally Posted by Lorenzo210 View Post
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
Code:
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.

Code:
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

Last edited by saicharan662000@gmail.com; September 7, 2022 at 02:40.
saicharan662000@gmail.com is offline   Reply With Quote

Old   September 12, 2022, 07:39
Default
  #18
Member
 
hari charan
Join Date: Sep 2021
Location: India,hyderabad
Posts: 96
Rep Power: 4
saicharan662000@gmail.com is on a distinguished road
Hey lorenzo,
Can you attach your case file?
Thanks in advance
saicharan662000@gmail.com is offline   Reply With Quote

Old   December 14, 2023, 20:07
Question
  #19
New Member
 
Bruno B W
Join Date: Aug 2023
Posts: 5
Rep Power: 2
brunobw is on a distinguished road
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.
brunobw is offline   Reply With Quote

Reply

Tags
flow boiling, intercondensatingevaporat, phase change, vof

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
VoF Model, miscible fluids, mass transfer rate yulu FLUENT 1 March 24, 2018 10:59
Eulerian Multiphase Interfacial Species Mass Transfer thomasduffy001_ Fluent Multiphase 1 October 3, 2017 07:21
Interphase mass transfer of a reaction cfx_ws1992 Main CFD Forum 0 May 15, 2017 21:42
Difficulty In Setting Boundary Conditions Moinul Haque CFX 4 November 25, 2014 17:30
Mass and het transfer using eulerian model Pablo FLUENT 0 February 16, 2007 08:36


All times are GMT -4. The time now is 23:28.