CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Difference between UEqn in sonicFoam & Theory (http://www.cfd-online.com/Forums/openfoam-programming-development/119333-difference-between-ueqn-sonicfoam-theory.html)

sasanghomi June 14, 2013 18:13

Difference between UEqn in sonicFoam & Theory
 
Hi Foamers ,

I found a difference between UEqn in sonicFoam and Ferziger's book (Computational Methods for fluid Dynamics ).
In the code for UEqn in sonicFoam we have :
Code:

fvVectorMatrix UEqn
(
    fvm::ddt(rho, U)
    + fvm::div(phi, U)
    + turbulence->divDevRhoReff(U)
);

solve(UEqn == -fvc::grad(p));

This code ensures below equation :

\frac{\partial \rho u_{i}}{\partial t} + \frac{\partial \rho u_{i}u_{j}}{\partial x_{j}} = -\frac{\partial p}{\partial x_{i}} + \frac{\partial }{\partial x_{j}} (\mu _{eff}(\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}})) - \frac{2}{3}\mu _{eff}\frac{\partial }{\partial x_{j}} (\frac{\partial u_{k}}{\partial x_{k}}\delta _{ij})

But in the Ferziger's book (generally in theory) there is another equation :

\frac{\partial \rho u_{i}}{\partial t} + \frac{\partial \rho u_{i}u_{j}}{\partial x_{j}} = -\frac{\partial p}{\partial x_{i}} + \frac{\partial }{\partial x_{j}} (\mu _{eff}(\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}})) - \frac{2}{3}\frac{\partial }{\partial x_{j}} (\rho k \delta _{ij})-\frac{2}{3}\frac{\partial }{\partial x_{j}}(\mu _{eff}\frac{\partial u_{k}}{\partial x_{k}}\delta _{ij})

The last term in two equations are different . Please guide me that what is the origin of this difference ? And please correct me if I am wrong.

I appreciate any help from you.
Thanks and best regards,
Sasan.

vwibaut June 26, 2013 09:20

Your ferziger's equation comes from incompressible theory, no? You haven't any term
dui/dxi

immortality June 26, 2013 11:25

I think so ;)

lfgmarc June 27, 2013 01:34

Dear sasanghomi,

I found the same incongruence some time ago. I think that is a bug, but if you have some time you can review the
Code:

divDevRhoReff
member function to try to figure out this problem. I try to do this, but my conclusion was that it is an error, because the term

\frac{2}{3}\bar{\rho}\widetilde{k}\mathcal{I}

is missing :confused:

Best regards.

Felipe G.

vwibaut June 27, 2013 02:54

Quote:

Originally Posted by lfgmarc (Post 436240)
Dear sasanghomi,

I found the same incongruence some time ago. I think that is a bug, but if you have some time you can review the
Code:

divDevRhoReff
member function to try to figure out this problem. I try to do this, but my conclusion was that it is an error, because the term

\frac{2}{3}\bar{\rho}\widetilde{k}\mathcal{I}

is missing :confused:

Best regards.

Felipe G.


The term with 2/3 rho k is present in the equation. It is taken into account in the pressure (in the equation p = ps + 2/3 rho k). The results you have, after simulations, give you a pressure wich contains the real pressure and turbulent kinetic energy.

Chris Lucas June 27, 2013 03:09

Hi,

have a closer look at the function:

divDevRhoReff
and you will find your "missing" term

Regards,
Christian

lfgmarc June 27, 2013 03:48

Hi to all,

Can you say me the OF version ? because at least in OF1.7.0 I can't find this term.

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

Maybe I looked in a wrong location in the code.

sonicFoam pressure equation.

Code:

rho = thermo.rho();

volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();

surfaceScalarField phid
(
    "phid",
    fvc::interpolate(psi)
  *(
        (fvc::interpolate(U) & mesh.Sf())
      + fvc::ddtPhiCorr(rUA, rho, U, phi)
    )
);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
    fvScalarMatrix pEqn
    (
        fvm::ddt(psi, p)
      + fvm::div(phid, p)
      - fvm::laplacian(rho*rUA, p)
    );

    pEqn.solve();

    if (nonOrth == nNonOrthCorr)
    {
        phi = pEqn.flux();
    }
}

#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"

divDevRhoReff kEpsilon :

Code:


tmp<fvVectorMatrix> kEpsilon::divDevRhoReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(muEff(), U)
      - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
    );
}

Thanks in advance.


Best regards.

sasanghomi June 27, 2013 03:50

Hi Guys ,
I corrected the momentum equation in first post (Theory equation). But I think below term is missed.

\frac{2}{3}\frac{\partial }{\partial x_{j}} (\rho k \delta _{ij})
Best regards
sasan.

sasanghomi June 27, 2013 03:56

Quote:

Originally Posted by vwibaut (Post 436257)
The term with 2/3 rho k is present in the equation. It is taken into account in the pressure (in the equation p = ps + 2/3 rho k). The results you have, after simulations, give you a pressure wich contains the real pressure and turbulent kinetic energy.

Hi Valentin ,
Which file have this equation for pressure ?
Can you explain more and clarify the issue ?

Thanks and best regards,
Sasan.

vwibaut June 27, 2013 04:12

You will not find this term because it is into the pressure term. More clearly, the pressure in your equation isn't the real pressure. It is the real pressure+ 2/3 rho k.

But it is close to the real pressure because k << p


All times are GMT -4. The time now is 15:49.