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

evapVOFHardt discussion! Come join!

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

Like Tree3Likes
  • 1 Post By jameswilson620
  • 2 Post By Swagga5aur

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 24, 2015, 20:29
Default evapVOFHardt discussion! Come join!
  #1
Member
 
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 38
Rep Power: 8
jameswilson620 is on a distinguished road
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:

Code:
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:
Code:
hCSource = fvc::ddt(rhoCp) + fvc::div(rhoCpPhi);
and
Code:
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:

Code:
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!
Attached Files
File Type: pdf stefan_twophase_GUPTA_interface.pdf (56.3 KB, 175 views)
File Type: txt stefan_twophaseDOTm.txt (1.9 KB, 88 views)
Elham likes this.
jameswilson620 is offline   Reply With Quote

Old   May 6, 2015, 12:45
Default EvapVODHArdt
  #2
Member
 
Anastasios
Join Date: Mar 2009
Posts: 34
Rep Power: 13
ageorg is on a distinguished road
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
ageorg is offline   Reply With Quote

Old   August 14, 2015, 13:48
Default
  #3
Senior Member
 
Przemek
Join Date: Jun 2011
Posts: 234
Rep Power: 12
gaza is on a distinguished road
Quote:
Originally Posted by jameswilson620 View Post
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.
__________________
best regards
pblasiak

Last edited by wyldckat; August 14, 2015 at 19:23. Reason: merged posts that were a few minutes apart and that quote the same post
gaza is offline   Reply With Quote

Old   January 13, 2020, 18:30
Default evapVOFHardt OFV6
  #4
Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 97
Rep Power: 7
Swagga5aur is on a distinguished road
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
gaza and FerdiFuchs like this.

Last edited by Swagga5aur; January 17, 2020 at 04:20.
Swagga5aur is offline   Reply With Quote

Old   March 30, 2020, 04:52
Default evapVOFHardt
  #5
New Member
 
Muyiwa
Join Date: Feb 2020
Posts: 11
Rep Power: 3
Muyiwa is on a distinguished road
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.
Muyiwa is offline   Reply With Quote

Old   March 30, 2020, 05:17
Default
  #6
Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 97
Rep Power: 7
Swagga5aur is on a distinguished road
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
Swagga5aur is offline   Reply With Quote

Old   March 30, 2020, 21:58
Default evapVOFHardt
  #7
New Member
 
Muyiwa
Join Date: Feb 2020
Posts: 11
Rep Power: 3
Muyiwa is on a distinguished road
Attached are compilation errors from wmake of evapVOFHardt.
Attached Images
File Type: jpg wmake error_evapVOFHardt1.jpg (194.2 KB, 12 views)
File Type: jpg wmake error_evapVOFHardt2.jpg (122.8 KB, 9 views)
Muyiwa is offline   Reply With Quote

Old   April 1, 2020, 07:26
Default
  #8
Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 97
Rep Power: 7
Swagga5aur is on a distinguished road
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
Swagga5aur is offline   Reply With Quote

Old   April 17, 2020, 00:01
Default evapVOFHardt
  #9
New Member
 
Muyiwa
Join Date: Feb 2020
Posts: 11
Rep Power: 3
Muyiwa is on a distinguished road
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
Muyiwa is offline   Reply With Quote

Old   April 17, 2020, 04:03
Default
  #10
Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 97
Rep Power: 7
Swagga5aur is on a distinguished road
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
Swagga5aur is offline   Reply With Quote

Old   April 18, 2020, 01:57
Default evapVOFHardt
  #11
New Member
 
Muyiwa
Join Date: Feb 2020
Posts: 11
Rep Power: 3
Muyiwa is on a distinguished road
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
Muyiwa is offline   Reply With Quote

Old   April 18, 2020, 02:12
Default
  #12
Senior Member
 
Przemek
Join Date: Jun 2011
Posts: 234
Rep Power: 12
gaza is on a distinguished road
Hi
Just run in terminal
Code:
setFields
to do this you have to have setFieldsDict file in system directory.
Befor runing it edit setFieldsDict accordingly to your case.
__________________
best regards
pblasiak
gaza is offline   Reply With Quote

Old   May 29, 2020, 10:03
Default The last term in the energy equation
  #13
New Member
 
Join Date: Mar 2020
Posts: 11
Rep Power: 2
Huihui Xiao is on a distinguished road
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.
Huihui Xiao is offline   Reply With Quote

Old   June 23, 2020, 07:44
Default
  #14
New Member
 
Yağmur GÜLEÇ
Join Date: May 2014
Posts: 9
Rep Power: 8
yagmur_89 is on a distinguished road
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?
yagmur_89 is offline   Reply With Quote

Old   September 24, 2020, 11:18
Default Modeling dorplet evaporation with evapVOFHardt solver
  #15
New Member
 
Join Date: Mar 2020
Posts: 11
Rep Power: 2
Huihui Xiao is on a distinguished road
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
Huihui Xiao is offline   Reply With Quote

Old   September 29, 2020, 16:19
Default
  #16
New Member
 
Yağmur GÜLEÇ
Join Date: May 2014
Posts: 9
Rep Power: 8
yagmur_89 is on a distinguished road
Quote:
Originally Posted by ageorg View Post
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
yagmur_89 is offline   Reply With Quote

Old   November 29, 2020, 12:54
Default
  #17
New Member
 
M Naarendharan
Join Date: Aug 2020
Posts: 19
Rep Power: 2
Naaren is on a distinguished road
Quote:
Originally Posted by yagmur_89 View Post
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
Naaren is offline   Reply With Quote

Reply

Tags
boiling, interfoam, phase change, source terms, validation

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
Let us make a skype discussion group for online discussion on Openfoam 13msmemusman OpenFOAM 4 April 27, 2014 08:52
Can't Join Connectors AeroCat Pointwise & Gridgen 2 January 25, 2014 12:04
Pointwise join problem travis.moulding Pointwise & Gridgen 1 December 11, 2011 17:42
New Discussion Forum Launched pete Site News & Announcements 0 March 15, 2009 15:55
How to join a news group? maximus Main CFD Forum 1 January 15, 2003 20:05


All times are GMT -4. The time now is 22:47.