[s9a9m92001]

Hi there,

This is my first post. I read a academic paper "A Partial Average Based Study of Compressible Turbulent Flows"
http://www.jomse.org/paperInfo.aspx?ID=88

For incompressible flow, the governing equations are shown in the figure.

Gao_yong equations.PNG

Equation (1) and (2) are averaged continuity and momentum equations.
Equation (3) and (4) are fluctuation continuity and momentum equations.
Equation (5) is the turbulent energy dissipation equation.

The equations describe a new turbulence model without any empirical coefficients, which can simulate statistical mean behaviors and coherent structures of various benchmark turbulent flows.

The chain of momentum transference contains two stages.The first stage is from the mean flow to the fluctuation flow as the viscous dissipation term of the mean flow equals the source term of the fluctuation flow. The second stage is from the fluctuation flow to molecular motion as the viscous dissipation term of the fluctuation flow plays its role.

I have followed the tutorial that added temperature to icoFoam, now I plan to modify pisoFoam solver to solve the equations in OpenFoam.

Any suggestion is really appreciated.

Sam Chio

[s9a9m92001]

As I known, the equations used by pisoFOAM solver are similar to equation (1) and (2). So I can keep the velocity predictor UEqn using the following code:

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

UEqn.relax();

if (momentumPredictor)
{
   solve(UEqn == -fvc::grad(p));
}

The main difference of the new model is the calculation of Reynolds stress divDevReff(U). How to obtain Reynolds stress by solving equation (3) to (5) in OpenFoam? That is my question that I wish to solve.

[s9a9m92001]

Dear all,

I consider to build up the new turbulence solver by modifying the RASModel. As I mentioned in foregoing thread, OpenFOAM calculates the turbulent shear stresses in incompressible RASModel using the divDevReff term:

tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
    );
}

Noting that 'nuEff' is the so-called effective viscosity that is the sum of the laminar molecular viscosity and turbulent eddy viscosity. The expression is in the form:

virtual tmp<volScalarField> nuEff() const
{
    return tmp<volScalarField>
    (
      new volScalarField("nuEff", nut() + nu())
    );
}

Investigating the new turbulent model without empirical coefficients, I prefer to implement equations (3), (4) and (5) to calculate the value of "nut".

Sam

[s9a9m92001]

Next it is required to solve drift velocity udrift, drift pressure pdrift and turbulence length scale l. Header file createField.h may take the form like:

    udrift_
    (
        IOobject
        (
            "udrift",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        autoCreateUdrift("udrift", mesh_)
    ),

    pdrift_
    (
        IOobject
        (
            "pdrift",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        autoCreatePdrift("pdrift", mesh_)
    ),

    l_
    (
        IOobject
        (
            "l",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        autoCreateL("l", mesh_)
    ),

    nut_
    (
        IOobject
        (
            "nut",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        autoCreateNut("nut", mesh_)
    )

[s9a9m92001]

Considering the C++ code of equation (4), we have the expression in the form

fvVectorMatrix UdriftEqn
(
   fvm::ddt(udrift)
 + fvm::div(phi, udrift)
 + fvm::div(phidrift, U)
 + turbulence->divDevReff(udrift)
 - turbulence->divDevReff(U)
);

UdriftEqn.relax();

if (momentumPredictor)
{
   solve(UdriftEqn == -fvc::grad(pdrift));
}

[sharonyue]

It has been a while. how is it going dude?
AFAIK, Prof. Gao's student is doing the same thing and their works had been published in Chinese Journal.

[sharonyue]

Originally Posted by s9a9m92001

Considering the C++ code of equation (4), we have the expression in the form

fvVectorMatrix UdriftEqn
(
   fvm::ddt(udrift)
 + fvm::div(phi, udrift)
 + fvm::div(phidrift, U)
 + turbulence->divDevReff(udrift)
 - turbulence->divDevReff(U)
);

UdriftEqn.relax();

if (momentumPredictor)
{
   solve(UdriftEqn == -fvc::grad(pdrift));
}

I think it should be:

+ fvc::div(phidrift, U)

rite?