TEqn in compressibleInterFoam

sabrahams · August 9, 2023, 11:24

Hi everyone!

I wonder if anyone can help me to understand the formulation of the TEqn in compressibleInterFoam in OF9. The TEqn is:

Code:

fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T)
      - fvm::laplacian(turbulence.alphaEff(), T)
      + (
            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
          - (U()&(fvModels.source(rho, U)&U)()) - contErr*K 
        )
       *(
           alpha1()/mixture.thermo1().Cv()()
         + alpha2()/mixture.thermo2().Cv()()
        )
        
     ==
           fvModels.source(rho, T)

I believe

should be inside the derivatives since

is averaged over the two phases and can therefore vary in time and space.

I believe the TEqn should be equivalent to
$\frac{\partial}{\partial t}(\rho C_v T)+\nabla\cdot(\rho C_v U T) + \frac{\partial}{\partial t}(\rho k)+\nabla\cdot(\rho U k)+\nabla \cdot (pU)-\nabla\cdot (K \nabla T)=S,$

but the formulation above seems to be equivalent to
$\frac{C_{1}C_{2}}{{\alpha_1 C_{2}+\alpha_2 C_{1}}} \left[\frac{\partial}{\partial t}(\rho T)+\nabla\cdot(\rho U T) \right]+ \frac{\partial}{\partial t}(\rho k)+\nabla\cdot(\rho U k)+\nabla \cdot (pU)-\nabla\cdot (K \nabla T)=S.$

Can anyone help me to understand how

has been taken out in this way?

Many thanks,
SAbrahams

sabrahams · August 15, 2023, 10:34

Does no-one know the answer to this?

Tobermory · August 15, 2023, 11:55

I think the problem lies with your first Teqn equation ... Start from basics - the underlying energy equation is stated in terms of the internal energy, i.e. in conservative form:

$\frac{\partial \rho e}{\partial t} + \frac{\partial \rho U_i e}{\partial x_i} + ...$

which we can expand out using continuity to the non-conservative form:

$\rho \frac{\partial e}{\partial t} + \rho U_i \frac{\partial e}{\partial x_i} + ...$
Now using the definition of the specific heat, we write

, so substitute and you have:

$\rho C_v \frac{\partial T}{\partial t} + \rho C_v U_i \frac{\partial T}{\partial x_i} + ...$

or in terms of the material derivative:
$\rho C_v \frac{DT}{Dt} + ...$

with the specific heat capacity on the outside. Hopefully that helps.

sabrahams · August 15, 2023, 12:22

Hi Tobermory,

Thanks for your reply. However, since

is not a constant in this case, I don't believe we can say

. Instead, it should be

.

In a two phase system, as we have in compressibleInterFoam,

is averaged over the two phases, that is
$C_v = \alpha_1 C_{v,1}+\alpha_2 C_{v,2}.$

While $C_{v,1}$ and $C_{v,2}$ may be constant,

is not since $\alpha_1$ and $\alpha_2$ vary in time and space.

Furthermore, TEqn in compressibleInterFoam (as shown in my original post) is not
$\rho C_v \frac{\partial T}{\partial t}+ ...$
it seems to be (if you divide through by $\alpha_1/C_{v,1}+\alpha_2/C_{v,2}$ , which is multiplying the second term in TEqn)
$\frac{C_{v,1} C_{v,2}}{\alpha_1 C_{v,2}+\alpha_2 C_{v,1}} \frac{\partial \rho T}{\partial t}+ ...$

So I'm trying to understand how the energy equation has been rearranged to get the form in the TEqn as in my original post. Do you know how this is derived?

Thanks again,
SAbrahams

Tobermory · August 16, 2023, 12:32

Okay - let me break this down. First up, the following is not correct I am afraid:

Quote:

Originally Posted by sabrahams

However, since

is not a constant in this case, I don't believe we can say

. Instead, it should be

.

You need to start with the definition of the specific heat capacity, which is $C_v=(\frac{\partial e}{\partial T})_V$ ... and now can you see where I got my original

from? We have said nothing so far about the variability of Cv in space, just that the process is assumed to be a constant volume process (which makes sense for two immiscible fluids). If you're still not convinced, then think of the chain rule:

$\frac{\partial e}{\partial t} = \frac{\partial e}{\partial T} \frac{\partial T}{\partial t} = C_v \frac{\partial T}{\partial t}$

etc. Now consider the single phase energy equation, which in OpenFOAM is (refer to buoyantPimpleFoam for example):

$\frac{\partial \rho e}{\partial t} + \frac{\partial \rho U_j e}{\partial x_j} + \frac{\partial \rho K}{\partial t} + \frac{\partial \rho U_j K}{\partial x_j} + \frac{\partial U_j p}{\partial x_j} - \frac{\partial}{\partial x_j}(\rho\alpha_{eff} \frac{\partial e}{\partial x_j})$

expand out the first two terms into the nonconservative form using continuity, substitute for $\partial e$ and convert back to conservative form and you get:

$\frac{\partial \rho T}{\partial t} + \frac{\partial \rho U_j T}{\partial x_j} + \frac{1}{C_v}[\frac{\partial \rho K}{\partial t} + \frac{\partial \rho U_j K}{\partial x_j} + \frac{\partial U_j p}{\partial x_j}] - \frac{1}{C_v}\frac{\partial}{\partial x_j}(\rho C_v \alpha_{eff} \frac{\partial T}{\partial x_j})$

NOW a minor approximation is made in the laplacian, where it assumed that local variations in Cv are small, so that the terms cancel leaving:

$\frac{\partial \rho T}{\partial t} + \frac{\partial \rho U_j T}{\partial x_j} + \frac{1}{C_v}[\frac{\partial \rho K}{\partial t} + \frac{\partial \rho U_j K}{\partial x_j} + \frac{\partial U_j p}{\partial x_j}] - \frac{\partial}{\partial x_j}(\rho \alpha_{eff} \frac{\partial T}{\partial x_j})$

and we are almost there. All that is left is to realise that the finite volume discretised form of the energy (or temperature) budget equation is the sum of alpha1 times the equation for phase 1 plus alpha2 times that for phase two, which results in:

$(\alpha_1 + \alpha_2)[\frac{\partial \rho T}{\partial t} + \frac{\partial \rho U_j T}{\partial x_j}]$ $+ (\frac{\alpha_1}{C_{v1}}+\frac{\alpha_2}{C_{v2}})[\frac{\partial \rho K}{\partial t} + \frac{\partial \rho U_j K}{\partial x_j} + \frac{\partial U_j p}{\partial x_j}] - (\alpha_1 + \alpha_2)\frac{\partial}{\partial x_j}(\rho \alpha_{eff} \frac{\partial T}{\partial x_j})$

or on noting that $\alpha_1 + \alpha_2=1$ , we have:

$\frac{\partial \rho T}{\partial t} + \frac{\partial \rho U_j T}{\partial x_j} + (\frac{\alpha_1}{C_{v1}}+\frac{\alpha_2}{C_{v2}})[\frac{\partial \rho K}{\partial t} + \frac{\partial \rho U_j K}{\partial x_j} + \frac{\partial U_j p}{\partial x_j}] - \frac{\partial}{\partial x_j}(\rho \alpha_{eff} \frac{\partial T}{\partial x_j})$

noting that these should now not be differentials etc., but should instead be volume integrals of the differentials yielding volume average values or surface fluxes (I didn't have the time to tidy up the terminology). This is the form that appears in TEqn.h for compressibleInterFoam.

There are clearly a bunch of other implicit assumptions here ... but this is my reverse engineering of the code. I couldn't find any explicit reference to the source of the model equations, and suspect it was authored by one of the OpenFOAM team. If you do find a reference - share it with the forum please.

sabrahams · August 17, 2023, 11:38

That makes perfect sense. Thank you for the clear explanation!

I haven't found any references for the equations used in compressibleInterFoam in OF9. If I find anything I'll certainly share. UEqn and alphaEqn appear to be in the same form as equations (7) and (25) of Shi et al.. Although, this paper was published later than OF9 so isn't a reference but just a useful resource to understand the equations (for anyone else reading this).

Thanks again Tobermory.

sabrahams · August 23, 2023, 11:36

Hi again Tobermory,

I've just been working through this and I'm stuck on this step:

Quote:

Originally Posted by Tobermory

All that is left is to realise that the finite volume discretised form of the energy (or temperature) budget equation is the sum of alpha1 times the equation for phase 1 plus alpha2 times that for phase two, which results in:

$(\alpha_1 + \alpha_2)[\frac{\partial \rho T}{\partial t} + \frac{\partial \rho U_j T}{\partial x_j}]$ $+ (\frac{\alpha_1}{C_{v1}}+\frac{\alpha_2}{C_{v2}})[\frac{\partial \rho K}{\partial t} + \frac{\partial \rho U_j K}{\partial x_j} + \frac{\partial U_j p}{\partial x_j}] - (\alpha_1 + \alpha_2)\frac{\partial}{\partial x_j}(\rho \alpha_{eff} \frac{\partial T}{\partial x_j})$

I believe that the equation for phase 1 would have $\rho_1$ and the equation for phase 2 would have $\rho_2$ so I'm struggling to see how the two equations could be summed in the way that you've shown above. I'm getting an equation of the form:

$\alpha_1 [\frac{\partial \rho_1 T}{\partial t} + \frac{\partial \rho_1 U_j T}{\partial x_j}]$ $+ \alpha_2[\frac{\partial \rho_2 T}{\partial t} + \frac{\partial \rho_2 U_j T}{\partial x_j}] + ...$

Am I missing a step or an assumption somewhere?

Tobermory · August 24, 2023, 11:34

The way to think of it, I reckon, is as follows: we need to integrate the differential equations over the cell volume to get the finite volume expressions. Focus just on the first two terms, for simplicity:

$\frac{\partial \rho T}{\partial t} + \frac{\partial \rho U_j T}{\partial x_j} + ...$

integrate these to get:

$\int_V( \frac{\partial \rho T}{\partial t}).dV + \int_V(\frac{\partial \rho U_j T}{\partial x_j}).dV + ...$

$\frac{\partial }{\partial t}(\int_V \rho T . dV) + \int_S(\rho U_j T).dA + ...$

$\frac{\partial }{\partial t}(\rho_P T_P V_P) + \sum(\phi_f T_f) + ...$

where it's assumed that the mesh is not moving or morphing, $\phi_f$ are the face fluxes, and the P subscript denotes a volume-averaged value, e.g.

$\int_V T . dV = T_P V_p$

Note though that we are integrating over volume $\alpha_1 V_P$ for phase 1, and $\alpha_2 V_P$ for phase 2. This is the origin of the $\alpha_1$ and $\alpha_2$ coefficients in each of the terms. This is what I meant rather laconically in my earlier post by "noting that these should now not be differentials etc.". Hope this helps.

August 9, 2023, 11:24	TEqn in compressibleInterFoam	#1
sabrahams New Member S Abrahams Join Date: Mar 2022 Location: UK Posts: 15 Rep Power: 4	Hi everyone! I wonder if anyone can help me to understand the formulation of the TEqn in compressibleInterFoam in OF9. The TEqn is: Code: fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T) - fvm::laplacian(turbulence.alphaEff(), T) + ( fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rhop + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - (U()&(fvModels.source(rho, U)&U)()) - contErrK ) ( alpha1()/mixture.thermo1().Cv()() + alpha2()/mixture.thermo2().Cv()() ) == fvModels.source(rho, T) I believe should be inside the derivatives since is averaged over the two phases and can therefore vary in time and space. I believe the TEqn should be equivalent to $\frac{\partial}{\partial t}(\rho C_v T)+\nabla\cdot(\rho C_v U T) + \frac{\partial}{\partial t}(\rho k)+\nabla\cdot(\rho U k)+\nabla \cdot (pU)-\nabla\cdot (K \nabla T)=S,$ but the formulation above seems to be equivalent to $\frac{C_{1}C_{2}}{{\alpha_1 C_{2}+\alpha_2 C_{1}}} \left[\frac{\partial}{\partial t}(\rho T)+\nabla\cdot(\rho U T) \right]+ \frac{\partial}{\partial t}(\rho k)+\nabla\cdot(\rho U k)+\nabla \cdot (pU)-\nabla\cdot (K \nabla T)=S.$ Can anyone help me to understand how has been taken out in this way? Many thanks, SAbrahams Last edited by sabrahams; August 9, 2023 at 14:09. Reason: typos*

August 15, 2023, 12:22		#4
sabrahams New Member S Abrahams Join Date: Mar 2022 Location: UK Posts: 15 Rep Power: 4	Hi Tobermory, Thanks for your reply. However, since is not a constant in this case, I don't believe we can say . Instead, it should be . In a two phase system, as we have in compressibleInterFoam, is averaged over the two phases, that is $C_v = \alpha_1 C_{v,1}+\alpha_2 C_{v,2}.$ While $C_{v,1}$ and $C_{v,2}$ may be constant, is not since $\alpha_1$ and $\alpha_2$ vary in time and space. Furthermore, TEqn in compressibleInterFoam (as shown in my original post) is not $\rho C_v \frac{\partial T}{\partial t}+ ...$ it seems to be (if you divide through by $\alpha_1/C_{v,1}+\alpha_2/C_{v,2}$ , which is multiplying the second term in TEqn) $\frac{C_{v,1} C_{v,2}}{\alpha_1 C_{v,2}+\alpha_2 C_{v,1}} \frac{\partial \rho T}{\partial t}+ ...$ So I'm trying to understand how the energy equation has been rearranged to get the form in the TEqn as in my original post. Do you know how this is derived? Thanks again, SAbrahams Last edited by sabrahams; August 15, 2023 at 12:28. Reason: correcting mistake in equation

August 17, 2023, 11:38		#6
sabrahams New Member S Abrahams Join Date: Mar 2022 Location: UK Posts: 15 Rep Power: 4	That makes perfect sense. Thank you for the clear explanation! I haven't found any references for the equations used in compressibleInterFoam in OF9. If I find anything I'll certainly share. UEqn and alphaEqn appear to be in the same form as equations (7) and (25) of Shi et al.. Although, this paper was published later than OF9 so isn't a reference but just a useful resource to understand the equations (for anyone else reading this). Thanks again Tobermory. Tobermory likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Correction of T after TEqn in compressibleInterFoam	krikre	OpenFOAM Programming & Development	3	June 26, 2022 23:07
TEqn (passive scalar) diverges all of a sudden	backscatter	OpenFOAM Running, Solving & CFD	11	October 21, 2018 08:56
TEqn. in PorousSimpleFoam	svramana	Main CFD Forum	0	February 14, 2018 09:53
EEqn to TEqn	jishnuhari25	OpenFOAM Programming & Development	1	June 4, 2015 12:06
Replace H equation with Teqn in chtMultiregionFoam OF 2.3	Moncef	OpenFOAM Running, Solving & CFD	0	May 22, 2014 21:37

August 15, 2023, 10:34		#2
sabrahams New Member S Abrahams Join Date: Mar 2022 Location: UK Posts: 15 Rep Power: 4	Does no-one know the answer to this?

August 15, 2023, 11:55		#3
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 745 Rep Power: 14	I think the problem lies with your first Teqn equation ... Start from basics - the underlying energy equation is stated in terms of the internal energy, i.e. in conservative form: $\frac{\partial \rho e}{\partial t} + \frac{\partial \rho U_i e}{\partial x_i} + ...$ which we can expand out using continuity to the non-conservative form: $\rho \frac{\partial e}{\partial t} + \rho U_i \frac{\partial e}{\partial x_i} + ...$ Now using the definition of the specific heat, we write , so substitute and you have: $\rho C_v \frac{\partial T}{\partial t} + \rho C_v U_i \frac{\partial T}{\partial x_i} + ...$ or in terms of the material derivative: $\rho C_v \frac{DT}{Dt} + ...$ with the specific heat capacity on the outside. Hopefully that helps.

August 24, 2023, 11:34		#8
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 745 Rep Power: 14	The way to think of it, I reckon, is as follows: we need to integrate the differential equations over the cell volume to get the finite volume expressions. Focus just on the first two terms, for simplicity: $\frac{\partial \rho T}{\partial t} + \frac{\partial \rho U_j T}{\partial x_j} + ...$ integrate these to get: $\int_V( \frac{\partial \rho T}{\partial t}).dV + \int_V(\frac{\partial \rho U_j T}{\partial x_j}).dV + ...$ $\frac{\partial }{\partial t}(\int_V \rho T . dV) + \int_S(\rho U_j T).dA + ...$ $\frac{\partial }{\partial t}(\rho_P T_P V_P) + \sum(\phi_f T_f) + ...$ where it's assumed that the mesh is not moving or morphing, $\phi_f$ are the face fluxes, and the P subscript denotes a volume-averaged value, e.g. $\int_V T . dV = T_P V_p$ Note though that we are integrating over volume $\alpha_1 V_P$ for phase 1, and $\alpha_2 V_P$ for phase 2. This is the origin of the $\alpha_1$ and $\alpha_2$ coefficients in each of the terms. This is what I meant rather laconically in my earlier post by "noting that these should now not be differentials etc.". Hope this helps.