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

question about fvm::laplacian(turbulence->alphaEff(), he) in rhoSimpleFoam

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

Like Tree8Likes
  • 3 Post By Tobermory
  • 2 Post By Bloerb
  • 2 Post By mAlletto
  • 1 Post By mAlletto

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   December 3, 2020, 07:49
Default question about fvm::laplacian(turbulence->alphaEff(), he) in rhoSimpleFoam
  #1
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 448
Rep Power: 9
mAlletto is on a distinguished road
In rhoSimpleFoam and other solvers as well the energy equation has the form:


Code:
     volScalarField& he = thermo.he();    

   fvScalarMatrix EEqn     (         fvm::div(phi, he)   

    + (             he.name() == "e"         

  ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))       

     : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))         )       - 

 fvm::laplacian(turbulence->alphaEff(), he)      ==     

    fvOptions(rho, he)     );
The if statement is to use the sensible enthalpy h or sensible internal energy e variable to be solved for.


According to the book of Wilcox "Wilcox, David C. Turbulence modeling for CFD. Vol. 3. La Canada, CA: DCW industries, 2006."


equation 5.54 the turbulent heat flux is defined as follows:


q_{t,i} = -\frac{\mu_t c_p}{Pr_t}\frac{\partial T}{\partial x_i} = -\frac{\mu_t}{Pr_t}\frac{\partial h}{\partial x_i} (1)


which is the equation used to calculate the turbulent heat flux for the case the sensible enthalpy h is solved for. The above equation can also written ( I hop my thermodynamic knowledge is correct):


q_{t,i} = -\frac{\mu_t c_p c_v}{Pr_t c_v}\frac{\partial T}{\partial x_i} =  -\frac{\mu_t c_p}{Pr_tc_v}\frac{\partial e}{\partial x_i} (2)


However if the sensible internal energy is solved for the turbulent heat flux is modeled like this:


q_{t,i}  =  -\frac{\mu_t }{Pr_t}\frac{\partial e}{\partial  x_i} (3)




So equation (2) and (3) differ by a factor of cp/cv.



So in may opinion two slightly different equation are solved whether one uses h or e. Regarding the rest of the terms in the equation they are equivalent regardless if one solves for h or e.



Or maybe I have overseen something.



Did someone else encounter the same doubts?
mAlletto is offline   Reply With Quote

Old   December 4, 2020, 06:31
Default
  #2
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 208
Rep Power: 6
Tobermory is on a distinguished road
Your equations are correct, but perhaps a little misleading. The turbulent Prandtl numbers in eqns 1 & 3 are different, since the definition of the turbulent diffusivity, alpha, is different depending on whether you are solving for h or e:
  • For h, then \alpha_h = \kappa / \rho C_v and Pr_h = \mu / (\rho \alpha_h).
  • Whereas for e we have: \alpha_e = \kappa / \rho C_p and Pr_e = \mu / (\rho \alpha_e).
  • And so: \alpha_h = (C_p / C_v) \alpha_e  = \gamma \alpha_e.
Tobermory is offline   Reply With Quote

Old   December 4, 2020, 06:36
Default
  #3
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 254
Blog Entries: 1
Rep Power: 10
dlahaye is on a distinguished road
cp/cv = gamma = 1.4 (for diatomic gas), and thus the difference matters.

Definition of Prandtl number involves c_p in case that energy equation is formulated in terms of enthalpy. No equivalent in terms of c_v seems to exist.

Not sure.
dlahaye is offline   Reply With Quote

Old   December 4, 2020, 06:42
Default
  #4
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 208
Rep Power: 6
Tobermory is on a distinguished road
Ok, let me try explain again. Ultimately you are trying to model q = -\kappa \nabla T, and you use the chain rule to write this as one of the following (I have left out the "at constant p" etc. limits, out of laziness):

\frac{\partial e}{\partial x} = \frac{\partial e}{\partial T} \frac{\partial T}{\partial x} =  C_v \frac{\partial T}{\partial x}

\frac{\partial h}{\partial x} = \frac{\partial h}{\partial T} \frac{\partial T}{\partial x} =  C_p \frac{\partial T}{\partial x}

From this you get either:

q = -\frac{\kappa}{C_v} \nabla e

or

q = -\frac{\kappa}{C_p} \nabla h

and you can manipulate these to get your expressions since \kappa = \rho\alpha_h C_p = \rho \alpha_e C_v.

The point is - the turbulent Prandtl number (an approximation) is different for the two approaches, and of course is not related to the moelcular properties of the gas. Do remember - this is a modelled diffusion term for the turbulent energy transport. And again, my point is that alphaT is different depending on whether you solve for h or e; the difference being the factor gamma.
oswald, dlahaye and mAlletto like this.
Tobermory is offline   Reply With Quote

Old   December 5, 2020, 03:28
Default
  #5
Senior Member
 
Join Date: Sep 2013
Posts: 269
Rep Power: 16
Bloerb will become famous soon enough
turbulence->alphaEff() returns lambda/c_p or lambda/c_v depending on the solution variable e or h
dlahaye and mAlletto like this.
Bloerb is offline   Reply With Quote

Old   December 5, 2020, 05:13
Default
  #6
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 448
Rep Power: 9
mAlletto is on a distinguished road
Oh I see.


The compressible turbulence model take also the thermodynamic model as template parameter.


alphaEff is calculated as follwos in the file EddyDiffusivity.H:


Code:
        //- Return the effective turbulent thermal diffusivity for enthalpy 
        //  [kg/m/s] 
        virtual tmp<volScalarField> alphaEff() const
        {
            return this->transport_.alphaEff(alphat());
        }
the function
Code:
alphaEff(alphat())
is called int the file heThermo.C:



Code:
template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::volScalarField>
Foam::heThermo<BasicThermo, MixtureType>::alphaEff
(
    const volScalarField& alphat
) const
{
    tmp<Foam::volScalarField> alphaEff(this->CpByCpv()*(this->alpha_ + alphat));
    alphaEff.ref().rename("alphaEff");
    return alphaEff;
}
and the function
Code:
 CpByCpv()
is defined in the same file as above:


Code:

template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::volScalarField>
Foam::heThermo<BasicThermo, MixtureType>::CpByCpv() const
{
    const fvMesh& mesh = this->T_.mesh();

    tmp<volScalarField> tCpByCpv
    (
        new volScalarField
        (
            IOobject
            (
                "CpByCpv",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh,
            dimless
        )
    );

    volScalarField& CpByCpv = tCpByCpv.ref();

    forAll(this->T_, celli)
    {
        CpByCpv[celli] = this->cellMixture(celli).CpByCpv
        (
            this->p_[celli],
            this->T_[celli]
        );
    }
and the function
Code:
CpByCpv(p,T)
is defined differently whether one takes the sensibleEnthalpy or internalEnergy.


For the sensibleEntalpy we find the function in sensibleEnthalpy.H (it returns 1)



Code:

            //- Cp/Cp []
            scalar CpByCpv
            (
                const Thermo& thermo,
                const scalar p,
                const scalar T
            ) const
            {
                return 1;
            }
and for sensibleInternalEnergy in sensibleInternalEnergy.H (it returns gamma)



Code:

            //- Ratio of specific heats Cp/Cv []
            scalar CpByCpv
            (
                const Thermo& thermo,
                const scalar p,
                const scalar T
            ) const
            {
                #ifdef __clang__
                // Using volatile to prevent compiler optimisations leading to
                // a sigfpe
                volatile const scalar gamma = thermo.gamma(p, T);
                return gamma;
                #else
                return thermo.gamma(p, T);
                #endif
            }
dlahaye and Tobermory like this.
mAlletto is offline   Reply With Quote

Old   December 5, 2020, 05:21
Default
  #7
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 448
Rep Power: 9
mAlletto is on a distinguished road
Ok so we can say the the same equations are solve regardless one take h or e as quantity you solve for.



The difference in both equation may lay in the numerical stability since different explicit source terms appear in the equations.


For this see also When choose sensibleInternalEnergy or sensibleEnthalpy in thermophysicalProperties
dlahaye likes this.
mAlletto is offline   Reply With Quote

Reply

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
Question Re Engineering Data Source imnull ANSYS 0 March 5, 2012 13:51
Transonic rhoSimpleFoam Equations eric.m.tridas OpenFOAM 3 January 25, 2012 10:52
internal field question - PitzDaily Case atareen64 OpenFOAM Running, Solving & CFD 2 January 26, 2011 15:26
Question about rhoSimpleFoam "if (transonic)" universez OpenFOAM 4 April 17, 2010 10:21
Poisson Solver question Suresh Main CFD Forum 3 August 12, 2005 04:37


All times are GMT -4. The time now is 00:01.