[gaza]

Hi Foamers,


I created simple 1D test case for condensation, which is basically the Stefan problem for which analytical solution exists. The test case is a rod filled with a water vapour (alpha.liquid = 0) at saturation temperature of 380.26 K. In the first case (StefanCond_left) I set:

• on the left wall T = 370.26 K, the right wall is adiabatic and for the rest of the surfaces, empty boundary condition is appllied . In this configuration one can observe rapid condensation and after time 0.1 s all vapour is condensed (alpha.liquid = 1)

In the second case (StefanCond_right) I changed only the side of the cold wall, namely:

• the left wall is adiabatic, the right wall is at T = 370.26 K and for the rest of the surfaces, empty boundary condition is appllied . In this configuration one can observe very slow condensation and after time 0.1 s only 0.0254598 of the vapour condensed (alpha.liquid = 0.0254598).

Is it a bug? The solution should be the mirror image. Can you explain why the condensation rate decreases so much while it shouldn't change? Here are the two mentioned test cases (35MB):


http://fluid.itcmp.pwr.wroc.pl/~pbla...fanCond.tar.gz

[gaza]

Hi again,
I did described above the two test cases in Fluent and there I got mirror image solutions when I only changed the side of the cooled wall. So in interCondensingEvaporatingFoam is bug. It is more interesting because I found the same strange behavior in interThermalPhaseChangeFoam, as described here

https://github.com/MahdiNabil/CFD-PC/issues/3

So it can be bug that can be inherent in all interFoam based family solvers.

[Andrea_85]

Hi,

i am also working with evaporation/condensation and i want to implement my own phase change model in openFOAM. I started by having a look at the existing solvers: interPhaseChangeFoam and interCondesingEvaporatingFoam.
I do not understand the energy formulation in interCondesingEvaporatingFoam, do you know why there are no source terms to account for extra energies due to evaporation/condensation?
By starting from single-phase equations and deriving the equation for the mixture i think source terms should appear....or am i missing something?


Best,
Andrea

[gaza]

Originally Posted by Andrea_85

Hi,

i am also working with evaporation/condensation and i want to implement my own phase change model in openFOAM. I started by having a look at the existing solvers: interPhaseChangeFoam and interCondesingEvaporatingFoam.
I do not understand the energy formulation in interCondesingEvaporatingFoam, do you know why there are no source terms to account for extra energies due to evaporation/condensation?
By starting from single-phase equations and deriving the equation for the mixture i think source terms should appear....or am i missing something?


Best,
Andrea

Hi Andrea,
source terms to account for extra energies due to evaporation/condensation are included in the definition of internal energy e
please see the file twoPhaseMixtureThermo.H where they defined
e1 = Cv1(T - Tsat) + Hv1
e2 = Cv2(T - Tsat) + Hv2
e = (alpha1*rho1*e1 + alpha2*rho2*e2)/(alpha1*rho1 + alpha2*rho2)

and then in the transportProperties dict we set Hv1 = 0 for the liquid phase

[Andrea_85]

Thanks for the reply,

however I am not yet convinced of the energy formulation. In my opinion the source term for the energy equation should be proportional to the specific mass that is evaporating/condinsing in the unit time (mDotAlphal() according to constant.C) multiply by the specific latent heat of vaporization/condensation. isn't it? here it does not seem to be the case....

Another question...the heat flux term, div(q), in the energy equation can be rewritten using Fourier law as

div(q) = -div(kappa*grad(T))

where k is the thermal conductivitiy. Then assuming e = Cv*T, where "e" is the internal energy, the term should read

-div(kappa/Cv * grad(e))

but in eEqn.H we have

- fvm::laplacian(kappaEff/cp, e)

So, why kappa is divided by Cp and not by Cv?...then Cv and Cp are not constant in space so i am not sure we can take them out of the gradient of T.

Last thing and i am sorry if i bother you.
The term

fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), e)

Does it have a mathematical derivation? or it is artificially added to improve stability (negative implicit source term)? a similar term is also added in the momentum eq. It is the continuity eq. which should be zero...

Best,

andrea

[gaza]

Originally Posted by Andrea_85

Thanks for the reply,

however I am not yet convinced of the energy formulation. In my opinion the source term for the energy equation should be proportional to the specific mass that is evaporating/condinsing in the unit time (mDotAlphal() according to constant.C) multiply by the specific latent heat of vaporization/condensation. isn't it? here it does not seem to be the case....

Hi Andrea,

I did not track all the equations and I cannot give the answer. The correct equations you can find here

http://www.tandfonline.com/doi/abs/1...07780903423908

Originally Posted by Andrea_85

Another question...the heat flux term, div(q), in the energy equation can be rewritten using Fourier law as

div(q) = -div(kappa*grad(T))

where k is the thermal conductivitiy. Then assuming e = Cv*T, where "e" is the internal energy, the term should read

-div(kappa/Cv * grad(e))

but in eEqn.H we have

- fvm::laplacian(kappaEff/cp, e)

So, why kappa is divided by Cp and not by Cv?...then Cv and Cp are not constant in space so i am not sure we can take them out of the gradient of T.

The energy equation that is implemented is almost the same as in the Anderson's book on page 77

ftp://210.212.172.242/Digital_Librar...20anderson.pdf

However in the interCondensingEvaporatingFoam we have e under laplacian
and your question is justified why we devide kappa by cp and not by cv. I do not understand it as well.

Originally Posted by Andrea_85

Last thing and i am sorry if i bother you.
The term

fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), e)

Does it have a mathematical derivation? or it is artificially added to improve stability (negative implicit source term)? a similar term is also added in the momentum eq. It is the continuity eq. which should be zero...

Best,

andrea

Regarding the term

fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), e)

I do not fully understand it but my interpretation is the same as yours, namely it is additional term that enhance convergence and it is zero because it is something like continuity equation. However I have some doubts because in the OpenFOAM, the continuity equation for incompressible fluid should be

fvm::Sp(fvc::ddt(rho) + fvc::div(phi), e)

or am I wrong?

[Andrea_85]

Hi,

i understand the formulation in Kunkelmann paper. There the source term in the energy eq. is calculated as the product of the evaporation/condensation rate and the latent heat, i.e.

J=(T-Tsat)/(Ri DeltaH) * DeltaH
which means that if the cell is not evaporating/condensing the source term is zero for that cell.

In interCondesingEvaporatingFoam the evaporation rate is calculated as

Revap = coeff*alpha1*rho1*max(T-Tsat,0)
where it is assumed alpha1=1 in the liquid.

Now if we have a cell completely full with liquid (alpha1=1) but at a temperature lower than Tsat, the liquid in the cell will not evaporate (as it should be, Revap=0). However, if i've understood the energy formulation in OF (and probably i've not), the energy for that cell will be calculated by adding the contribution of the latent heat to the internal energy even if the liquid is not evaporating,

i.e. by assuming we only have evaporation (Hf2=0) and alpha1=1

e = (T-Tsat)*(Cv+Hf1) < 0

which looks weird to me....but maybe i am missing something.


Regarding the energy equation, i derived it from the total energy equation, neglecting viscous dissipation and gravity term, and i got

ddt(rho*e) + div(rhoPhi*e) + ddt(rho*K) + div(rhoPhi*K) + div(pU) = + div(kappa*grad(T)) + Se

which is not exactly what is implemented is eEqn.H. The term p*div(U) that is in eEqn.H appears if you simplify the mechanical contribution using the momentum eq. Have a look here for example

http://cfd.direct/openfoam/energy-equation/


The continuity equation for the mixture is written as

ddt(rho)+div(rho*U)=0 where rho=alpha1*rho1+alpha2*rho2

In the OF syntax this reads

fvc::ddt(rho)+fvc::div(rhoPhi)

phi= volumetric flux
rhoPhi= mass flux = rho*phi

So the term fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), e) is the continuity eq. multiply by the internal energy.


Best,
Andrea

[gaza]

Hi Andrea,
Thanks a lot for explanation.
Maybe you know why does interCondensingEvaporatingFoam behave strangely as I described in post #1 ??

[Andrea_85]

Honestly, i do not know.
I'll do some tests in the next days, if i have news ill let you know.

andrea

[gaza]

Hi
It is very strange because I did it the same in Fluent and everything works well there.

[jaymeen721]

Originally Posted by gaza

Hi Foamers,


I created simple 1D test case for condensation, which is basically the Stefan problem for which analytical solution exists. The test case is a rod filled with a water vapour (alpha.liquid = 0) at saturation temperature of 380.26 K. In the first case (StefanCond_left) I set:

• on the left wall T = 370.26 K, the right wall is adiabatic and for the rest of the surfaces, empty boundary condition is appllied . In this configuration one can observe rapid condensation and after time 0.1 s all vapour is condensed (alpha.liquid = 1)

In the second case (StefanCond_right) I changed only the side of the cold wall, namely:

• the left wall is adiabatic, the right wall is at T = 370.26 K and for the rest of the surfaces, empty boundary condition is appllied . In this configuration one can observe very slow condensation and after time 0.1 s only 0.0254598 of the vapour condensed (alpha.liquid = 0.0254598).

Is it a bug? The solution should be the mirror image. Can you explain why the condensation rate decreases so much while it shouldn't change? Here are the two mentioned test cases (35MB):


http://fluid.itcmp.pwr.wroc.pl/~pbla...fanCond.tar.gz

Hi, Gaza, can you share your case as it is not available over this link. http://fluid.itcmp.pwr.wroc.pl/~pbla...fanCond.tar.gz

[gaza]

Originally Posted by jaymeen721

Hi, Gaza, can you share your case as it is not available over this link. http://fluid.itcmp.pwr.wroc.pl/~pbla...fanCond.tar.gz

yes you can download it here
https://thermores.pwr.edu.pl/pracown...asiak/download