[jameswilson620]

Hello all!

You will find some fantastic work submitted: http://www.ttd.tu-darmstadt.de/forsc...vofhard.en.jsp

that makes reference to the work of Dr. Kumklemann (http://tuprints.ulb.tu-darmstadt.de/...Kunkelmann.pdf) based on the theory presented by Hardt and Wondra in 2008.

I've been developing my own phase change model for some time and dealing with an unpredictable temperature field. The conservative form of my equation yielded great results in the absence of phase change but would barf due to the presence of the divergence in the velocity field.

The formulation of the temperature equation provided in the above link is as follows:

volScalarField k = alpha1*k1+(scalar(1)-alpha1)*k2;
// turbulent heat transfer is to be added some time

volScalarField rhoCp = alpha1*rho1*cp1+(scalar(1)-alpha1)*rho2*cp2;
rhoCp.oldTime() = alpha1.oldTime()*rho1*cp1+(scalar(1)-alpha1.oldTime())*rho2*cp2;

surfaceScalarField alphaPhi = (rhoPhi-phi*rho2)/(rho1-rho2);
surfaceScalarField rhoCpPhi = alphaPhi*(rho1*cp1-rho2*cp2)+phi*rho2*cp2;

// source terms due to mass sources (due to non divergence-free velocity field)
hCSource = fvc::ddt(rhoCp) + fvc::div(rhoCpPhi);

// energy equation
fvScalarMatrix TEqn
(
	fvm::ddt(rhoCp,T) + fvm::div(rhoCpPhi,T) - fvm::laplacian(k,T) - fvm::Sp(hESource,T) - fvm::Sp(hCSource,T) == -hESource*Tsat
);

TEqn.solve();

The source term: hCSource = fvc::ddt(rhoCp) + fvc::div(rhoCpPhi); accounts for the "enthalpy collapse" at the interface. This is discussed by Hardt and Wondra on page 6 of their journal article "Evaporation model for interfacial flows based on a continuum-field representation of the source terms". (http://www.sciencedirect.com/science...21999108001228)

I know there are some OF buffs out there who love to educate and this is why im posting here. Ive struggled with finding the approporiate form of this source term and here it is; however, I hope some of you will help in delineating what is happening here so im not blindly following a functioning piece of code.

Lets look at:

hCSource = fvc::ddt(rhoCp) + fvc::div(rhoCpPhi);

and

rhoCp.oldTime() = alpha1.oldTime()*rho1*cp1+(scalar(1)-alpha1.oldTime())*rho2*cp2;

I assume using alpha1.oldTime() etc. in rhoCP.oldTime() gives us access to the location of the quantity rhoCp prior to advecting the interface. I assume that this is necessary to ensure that fvc::ddt(rhoCp) calculates the difference in rhoCp between time steps explicitly. Now, I also assume that this change in rhoCp results from

1) convection of the quantity by phi and
2) the regression of the interface due to phase change

This isnt very clear to me since the change seems to undermine the purpose of including these terms to begin with but clearly I am not correct in assuming this..

Now for the term:

fvc::div(rhoCpPhi)

where rhoCpPhi = alphaPhi*(rho1*cp1-rho2*cp2)+phi*rho2*cp2; (this is just a creative and effective way of writing phi*(interpolate(alpha1)*rho1*cp1 + (1-interpolate(alpha1))*rho1*cp1) made available after solving the alpha equation). Here, fvc::div(rhoCpPhi) is using the current values of the interface location and consequently the new location of the material properties rho and cp.

The combination of these two terms form hCSource and is the interfacial boundary condition accounting for the destruction of liquid and creation of vapor in the presence of phase change. This condition appears as a correction in this solver but is more realistically a physical interfacial boundary condition and can be seen in the analytical solution of 1D Stefan type freezing/melting problems.

I have attached a plot verifying the validity of the source terms. NOTE that this plot is my solver with the addition of the source term hcSource. Without this source term, my results had ~3% error where I now have <<1% error where grid refinement/and temporal resolution only improves the results, i.e. im converging on the solution and not past it.

The theory supporting my claim that this hcSource term is a coupling boundary condition at the interface can be seen in eq. 1.3.6 of :

https://www.dropbox.com/s/p424otswg5...GUPTA.pdf?dl=0

Also note the corrections i made to the solution provided. I made a matlab file for those interested in my validation (See attached and please forgive the programming etiquette. my openfoam code is much cleaner : ) )

This is also the source for the analytical solution seen in the attachment. One important note, this validation is for constant density. If the density is constant and you use the solver provided in the links, the term:
surfaceScalarField alphaPhi = (rhoPhi-phi*rho2)/(rho1-rho2);
will result in divide by zero error.

For those who would like to use the solver provided, wyldcat has ported the referenced solver to work with various versions of OF: https://github.com/wyldckat/evapVOFHardt

i look forward to your responses!

[ageorg]

Dear James

I am working in a user defined code using almost the same approach as Kunkelmann in EvapVOFHard. The only difference is that i am using a slightly different approach for smoothing the alpha field and dampen out spurious currents (from the hydrodynamic adiabatic point of view). Then I added heat transfer and phase change following the same procedure as evapVOFHard from the work of Hardt and Wondra. I have so far validated the solver with experimental results on pool boiling. I am trying now to validate for cases of convective boiling but I notice very fast evaporation with respect to experimental results. What is your field of application? Maybe we can cooperate on this.

Thanks in advance

ageorg

[gaza]

Originally Posted by jameswilson620

Hello all!

You will find some fantastic work submitted: http://www.ttd.tu-darmstadt.de/forsc...vofhard.en.jsp
[...]

Hi jameswilson620,

Can you upload your solver?
I checked evapVOFHardt solver and it seems it does not work well.
I created thread on it here:

http://www.cfd-online.com/Forums/ope...pvofhardt.html


and see also here

http://www.cfd-online.com/Forums/ope...tml#post559592

Maybe you know the answer? Any help is really appreciated.

[Swagga5aur]

Hello everyone interested,

I have updated the evapVOFHardt solver to openFOAM v6 and inserted a link it in this post. Let me know if you have any issues or comments.

I will finish testing and comparing the solver with the results of the 2.3.x variant to validate the update in the following days.

Best regards, Lasse

https://github.com/Swagga5aur/evapVOFHardt

[Muyiwa]

I'm a new user of openfoam and I'm working on boiling. I'm using OpenFoam v1812. I downloaded your evapVOFHardt from https://github.com/Swagga5aur/evapVOFHardt but I don't know how to go about running it in my openfoam version. Can you explain to me how to run it.

[Swagga5aur]

I haven't used openFOAM by esi in quite some time, but I believe the differences are only minor. What compilation errors do you get?

Regards,
Lasse

[Muyiwa]

Attached are compilation errors from wmake of evapVOFHardt.

[Swagga5aur]

Hello Muyiwa,

I'm currently quite busy in life outside of openFOAM but I'll try to assist you when I have the time, no need to PM me.

I would prefer all correspondence and support regarded material, that aren't confidential, to be here in the forum for helping future users having similar issues.

I will try to get time to look at it tonight, but I would suggest you try compiling openFOAM v6 from the foundation (for comparison) or just look in the source code to see where the differences are between the standard solver the evapVOFHardt is based upon, I believe it was interDyMFoam.

Regards,
Lasse

[Muyiwa]

Hi Swagga5aur

The evapVOFHardt code works fine on ofv6. I use setFields dict for my new case with 2D rectangular domain instead of using initFieldVOFHardt but the simulation stopped after running for a while. Is it possible to introduce setFields dict into the evapVOFHardt code and if yes, how?

Regards

Muyiwa

[Swagga5aur]

What do you means by introducing setFields dict into the evapVOFHardt code? Are you currently unable to use the setFields command in the case?

Regards,
Lasse

[Muyiwa]

Hi Swagga5aur

Actually, I'm able to use the "initFieldVOFHardt" setField command for the example case provided with the solver. But, I want use the solver to simulate a film boiling according to Welch and Wilson, 2000 in "A Volume of Fluid Based Method for Fluid Flows with Phase Change" and I don't know how to incorporate the setField since the initialization of alpha.phase1 is different from the example case.

Regards,
Muyiwa

[gaza]

Hi
Just run in terminal

setFields

to do this you have to have setFieldsDict file in system directory.
Befor runing it edit setFieldsDict accordingly to your case.

[Huihui Xiao]

Hello,

I am interested in this model as well, and I read the source codes of the solver. Concerning the implementations of the latent heat source term in the energy equation, I have the following questions.

In the evapVOFHardt solver (TEqn.H), they calculated the latent heat term by decomposing -1*(mass source term)*(Enthalpy of vaporization) into hESource*T - hESource*Tsat, and calculate them as implicit and explicit terms respectively.

Actually, I know this method works, however, I tried to calculate the latent heat with a direct way, i.e. using the single explicit term ''-1*(mass source term)*(Enthalpy of vaporization) ''. But I got something strange there when I validate the solver with 1D stefan problem, and the problem is that I got the minimum temperature around 369.677 K (for the BCs in my case (Left wall temperature : 383.15 K, and the rest is 373.15K )), and the simulation was really slow although we can see the interface was moving.

Could you please tell me why I got the strange results? Thank you so much.

[yagmur_89]

I solved the problem using instead of using

dimensionedScalar k1(phase1.lookup("k"));
dimensionedScalar cp1(phase1.lookup("cp"));
Instead of dimensionedScalar k1=phase1.lookup("k");


I set up a case with a vapor bubble attached to a wall with contact angle boundary condition, but I have problems with higher temperature gradients, which results in a significantly higher growth rates? Do you have any idea that causes the problem like initial boundary layer surrounding the bubble? Should I prescribe temperature field around the interface of the bubble? Or should it work without this process?

[Huihui Xiao]

Hello everyone,

I compiled the evapVOFHardt solver, and try to model 2D static droplet evaporation. Actually, I successfully validated the solver with 1D Stefan problem, however, I failed to validate the solver with 2D static droplet evaporation case. As reported in some literature, this is probably because the non-divergence-free one-filed velocity U in the solver failed to correctly advect the free surface, and in some paper, someone tried to reconstruct the divergence-free liquid velocity field to advect the free surface (e.g. https://arxiv.org/abs/2001.03477).

Additionally, I tried to only set the source terms on the interface cells, and rather than removing the source terms from the interface cells implemented in the solver, and I found that the 2D droplet seems expand rather than shrink due to the evaporation.

I would like to know does anyone encounter the same problem? and any suggestions will be appreciated.

Thanks,
Huihui

[yagmur_89]

Originally Posted by ageorg

Dear James

I am working in a user defined code using almost the same approach as Kunkelmann in EvapVOFHard. The only difference is that i am using a slightly different approach for smoothing the alpha field and dampen out spurious currents (from the hydrodynamic adiabatic point of view). Then I added heat transfer and phase change following the same procedure as evapVOFHard from the work of Hardt and Wondra. I have so far validated the solver with experimental results on pool boiling. I am trying now to validate for cases of convective boiling but I notice very fast evaporation with respect to experimental results. What is your field of application? Maybe we can cooperate on this.

Thanks in advance

ageorg

Dear Anastasios,
I use the same code for growing bubbles attached to a surface. I have the same problem as you had, especially with conjugate heat transfer. Did you solve the problem? If so, could you explain what improvements you did?

Thanks

[Naaren]

Originally Posted by yagmur_89

I solved the problem using instead of using

dimensionedScalar k1(phase1.lookup("k"));
dimensionedScalar cp1(phase1.lookup("cp"));
Instead of dimensionedScalar k1=phase1.lookup("k");


I set up a case with a vapor bubble attached to a wall with contact angle boundary condition, but I have problems with higher temperature gradients, which results in a significantly higher growth rates? Do you have any idea that causes the problem like initial boundary layer surrounding the bubble? Should I prescribe temperature field around the interface of the bubble? Or should it work without this process?

Hello everyone,
I am working on growth of a single bubble. I would like to know if anyone has validated the solver "evapVOFHardt" for simulation of growth and departure of a single bubble attached to the base. I would also like to know whether the solver takes into account the changes in contact angle when the contact angle boundary condition is given for base.


Thanking you,
M Naarendharan