Difference between revisions of "User:Hoogs/buoyantPisoFoam"

From OpenFOAMWiki
(Created page with '=Overview= The <tt>buoyantPisoFoam</tt> solver is suitable for unsteady, compressible flow and includes the [http://en.wikipedia.org/wiki/Buoyancy gravitational body force] on t...')
 
(Fixed fluxes replaced dot with \cdot, removed OpenCFD [M[U]] notation)
Line 22: Line 22:
 
| enthalpy
 
| enthalpy
 
|-
 
|-
| <math> \underline{ \phi } \,</math>
+
| <math> \mathbf{\Phi} \,</math>
 +
|
 +
| mass flux, <math> = \rho \mathbf{ U }\,</math>
 +
|-
 +
| <math> \phi \,</math>
 
| <tt>phi</tt>
 
| <tt>phi</tt>
| mass flux through face, <math> = \rho \mathbf{ U } . \mathbf{ S }_{ f } \,</math>
+
| mass flux through faces, <math> = \rho \mathbf{ U } . \mathbf{ S }_{ f } \,</math>
 
|-
 
|-
 
| <math> \alpha_{ t } \,</math>
 
| <math> \alpha_{ t } \,</math>
 
| <tt>alphaEff</tt>
 
| <tt>alphaEff</tt>
 
| effective turbulent thermal diffusivity
 
| effective turbulent thermal diffusivity
|-
 
| <math> S_{ hd } \,</math>
 
| <tt>sourceEnthalpy</tt>
 
| enthalpy source in continuous phase resulting from particles
 
|-
 
| <math> S_{ ud } \,</math>
 
| <tt>sourceMomentum</tt>
 
| momentum source in continuous phase resulting from particles
 
 
|-
 
|-
 
| <math> \psi \,</math>
 
| <math> \psi \,</math>
Line 57: Line 53:
 
| <tt>divRhoR</tt>
 
| <tt>divRhoR</tt>
 
| Reynolds stresses <math> = \rho \overline{ U_{ i }' U_{ j }' } \,</math>
 
| Reynolds stresses <math> = \rho \overline{ U_{ i }' U_{ j }' } \,</math>
|-
 
| <math> N_{ d } \,</math>
 
| <tt>nParticle</tt>
 
| Number of particles in dispersed phase parcel
 
|-
 
| <math> m_{ d } \,</math>
 
| <tt>mass0</tt>
 
| Mass of particles in dispersed phase parcel
 
|-
 
| <math> U_{ d0 } \,</math>
 
| <tt>U0</tt>
 
| Velocity of dispersed phase parcel at beginning of time step
 
|-
 
| <math> U_{ d1 } \,</math>
 
| <tt>U1</tt>
 
| Velocity of dispersed phase parcel at end of time step
 
|-
 
| <math> T_{ d0 } \,</math>
 
| <tt>T0</tt>
 
| Temperature of dispersed phase parcel at beginning of time step
 
|-
 
| <math> T_{ d1 } \,</math>
 
| <tt>T1</tt>
 
| Temperature of dispersed phase parcel at end of time step
 
|-
 
| <math> c_{ pd } \,</math>
 
| <tt>cp0</tt>
 
| Specific heat capacity of dispersed phase parcel
 
|-
 
| <math> \Delta_{ 01 } \,</math>
 
|
 
| Difference operator acting on a parameter that changes during a timestep, e.g. <math> \Delta_{ 01 } U_{ d } = U_{ d0 } - U_{ d1 } \,</math>
 
|-
 
| <math> \Delta_{ cd } \,</math>
 
|
 
| Difference operator acting on a parameter common to the continuous (c) and dispersed (d) phases, e.g. <math> \Delta_{ cd } U = U_{ c } - U_{ d } \,</math>
 
|-
 
 
|}
 
|}
  
Line 99: Line 58:
  
 
The primary [http://en.wikipedia.org/wiki/Navier_Stokes#Compressible_flow_of_Newtonian_fluids Navier Stokes] equation is for the transport of momentum with sources and sinks.
 
The primary [http://en.wikipedia.org/wiki/Navier_Stokes#Compressible_flow_of_Newtonian_fluids Navier Stokes] equation is for the transport of momentum with sources and sinks.
 +
 +
<center>
 +
<math>
 +
\frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla \cdot ( \rho \mathbf{U} \mathbf{U} ) + \nabla \cdot ( \rho \mathbf{R} ) = -\nabla p + \rho \mathbf{g}
 +
</math>
 +
</center>
 +
 +
or in terms of the flux which must be discretised at cell faces (unlike most other quantities which are discretised at cell centres),
  
 
<center>
 
<center>
 
{| cellpadding="2" style="border:2px solid #ccccff"
 
{| cellpadding="2" style="border:2px solid #ccccff"
 
|<math>
 
|<math>
\frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla . ( \underline{\phi} \mathbf{U} ) + \nabla . ( \rho \mathbf{R} ) = -\nabla p + \rho \mathbf{g}
+
\frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla \cdot ( \mathbf{\Phi} \mathbf{U} ) + \nabla \cdot ( \rho \mathbf{R} ) = -\nabla p + \rho \mathbf{g}
 
</math>
 
</math>
 
|}
 
|}
Line 148: Line 115:
 
{| cellpadding="2" style="border:2px solid #ccccff"
 
{| cellpadding="2" style="border:2px solid #ccccff"
 
|<math>
 
|<math>
\frac{\partial \rho}{\partial t} + \nabla . \underline{\phi} = 0
+
\frac{\partial \rho}{\partial t} + \nabla \cdot \phi = 0
 
</math>
 
</math>
 
|}
 
|}
Line 168: Line 135:
 
{| cellpadding="2" style="border:2px solid #ccccff"
 
{| cellpadding="2" style="border:2px solid #ccccff"
 
|<math>
 
|<math>
\frac{\partial ( \rho h ) }{\partial t} + \nabla . ( \underline{\phi} h ) + \nabla . ( \alpha_{t} \nabla h ) = \frac{Dp}{Dt}
+
\frac{\partial ( \rho h ) }{\partial t} + \nabla \cdot ( \phi h ) + \nabla \cdot ( \alpha_{t} \nabla h ) = \frac{Dp}{Dt}
 
</math>
 
</math>
 
|}
 
|}
Line 275: Line 242:
 
=PISO=
 
=PISO=
  
The non linearity of the convective term of the momentum equation requires an iterative solution technique if we do not wish to be limited to small time steps.  But while the momentum equation provides a way to find the velocity in terms of the pressure gradient, we need some way to determine the pressure as we iterate.  In the Pressure-Implicit Split-Operator (PISO) approach, we formulate a linear relationship between the pressure gradient and the velocity in order to obtain a prediction of the velocity (which violates mass conservation), then apply the mass conservation requirement in order to obtain a pressure field that is mass conservative so that it can be fed back to obtain a corrected velocity field.  At the same time, we correct for face skewness by updating the mass fluxes at cell faces with the new pressure field.
+
The non linearity of the convection term of the momentum equation requires an iterative solution technique if we do not wish to be limited to small time steps.  But while the momentum equation provides a way to find the velocity in terms of the pressure gradient, we need some way to determine the pressure as we iterate.  In the Pressure-Implicit Split-Operator (PISO) approach, we formulate a linear relationship between the pressure gradient and the velocity in order to obtain a prediction of the velocity (which violates mass conservation), then apply the mass conservation requirement in order to obtain a pressure field that is mass conservative so that it can be fed back to obtain a corrected velocity field.  At the same time, we correct for face skewness by updating the mass fluxes at cell faces with the new pressure field.
  
The momentum equation is discretised to obtain
+
The momentum equation is discretised at <math>N</math> discrete locations to obtain
  
 
::<math>
 
::<math>
\mathbf{[\mathcal{M}}[\mathbf{U}]\mathbf{]} = -\nabla p
+
A \mathbf{U} - \mathbf{H}
 
</math>
 
</math>
  
where <math>\mathbf{[\mathcal{M}}[\mathbf{U}]\mathbf{]}</math> represents the square matrix system of coefficients and a source vector over <math>N</math> discrete locations,
+
or
  
 
::<math>
 
::<math>
Line 306: Line 273:
 
</math>
 
</math>
  
The <math>A</math> and <math>H</math> are then modified to decompose <math>\mathbf{[\mathcal{M}}[\mathbf{U}]\mathbf{]}</math> into a diagonal matrix and new square source matrix  
+
<math>A</math> and <math>\mathbf{H}</math> are then reconstituted to make <math>A</math> into a diagonal matrix with the "offcuts" placed into the square matrix <math>\mathbf{H}</math>  
 
+
::<math>
+
\mathbf{[\mathcal{M}}[\mathbf{U}]\mathbf{]} \rightarrow A \mathbf{U} - \mathbf{H}
+
</math>
+
  
 
Multiplying both sides by the inverse of <math>A</math> yields the ''momentum corrector equation''
 
Multiplying both sides by the inverse of <math>A</math> yields the ''momentum corrector equation''
Line 322: Line 285:
 
::<math>
 
::<math>
 
\begin{align}
 
\begin{align}
& \frac{\partial \rho}{\partial t} + \nabla . (\rho \mathbf{U}) = 0 \\
+
& \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{U}) = 0 \\
& \frac{\partial \rho}{\partial t} + \mathbf{U} . \nabla \rho + \rho \nabla . \mathbf{U} = 0 \\
+
& \frac{\partial \rho}{\partial t} + \mathbf{U} \cdot \nabla \rho + \rho \nabla \cdot \mathbf{U} = 0 \\
& \frac{\partial \rho}{\partial t} + \frac{\mathbf{H}}{A} . \nabla \rho - \frac{1}{A} \nabla p . \nabla \rho + \rho \nabla . \left( \frac{\mathbf{H}}{A} \right) - \rho \nabla . \left( \frac{1}{A} \nabla p \right) = 0 \\
+
& \frac{\partial \rho}{\partial t} + \frac{\mathbf{H}}{A} \cdot \nabla \rho - \frac{1}{A} \nabla p \cdot \nabla \rho + \rho \nabla \cdot \left( \frac{\mathbf{H}}{A} \right) - \rho \nabla \cdot \left( \frac{1}{A} \nabla p \right) = 0 \\
 
\end{align}
 
\end{align}
 
</math>
 
</math>
Line 332: Line 295:
 
::<math>
 
::<math>
 
\begin{align}
 
\begin{align}
\underline{\phi} & = \mathbf{S}_{f} . \mathbf{U}_{f} \\
+
\phi & = \mathbf{S}_{f} . \mathbf{U}_{f} \\
& = \mathbf{S}_{f} . \left[ \left( \frac{\mathbf{H}}{A} \right)_{f} - \left( \frac{1}{A} \right)_{f} ( \nabla p )_{f} \right] \\
+
& = \mathbf{S}_{f} \cdot \left[ \left( \frac{\mathbf{H}}{A} \right)_{f} - \left( \frac{1}{A} \right)_{f} ( \nabla p )_{f} \right] \\
 
\end{align}
 
\end{align}
 
</math>
 
</math>
Line 346: Line 309:
 
|-
 
|-
 
| align="center" | 1
 
| align="center" | 1
| Accept <math>p, \mathbf{U}, \underline{\phi}</math>, and require <math>\rho</math>
+
| Accept <math>p, \mathbf{U}, \phi</math>, and require <math>\rho</math>
 
| <pre>
 
| <pre>
 
while (runTime.loop())
 
while (runTime.loop())
Line 356: Line 319:
 
|-
 
|-
 
| align="center" | 2
 
| align="center" | 2
| <math>\mathbf{[\mathcal{M}}[\mathbf{U}]\mathbf{]} = \frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla . ( \underline{\phi} \mathbf{U} ) + \nabla . ( \rho \mathbf{R} )</math>
+
| <math>\mathrm{LHS} = \frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla . ( \phi \mathbf{U} ) + \nabla . ( \rho \mathbf{R} )</math>
 
| <pre>
 
| <pre>
 
     fvVectorMatrix UEqn
 
     fvVectorMatrix UEqn
Line 368: Line 331:
 
|-
 
|-
 
| align="center" | 3
 
| align="center" | 3
| <math>\mathbf{[\mathcal{M}}[\mathbf{U}]\mathbf{]} = -\nabla p + \rho \mathbf{g}</math>
+
| <math>\mathrm{RHS} = -\nabla p + \rho \mathbf{g}</math>
 
| <pre>
 
| <pre>
 
     UEqn.relax();
 
     UEqn.relax();

Revision as of 05:35, 18 September 2009

1 Overview

The buoyantPisoFoam solver is suitable for unsteady, compressible flow and includes the gravitational body force on the fluid. The Richardson number is a dimensionless measure of the importance of buoyancy in a flow.

2 Nomenclature

 p \, p pressure
 \rho \, rho density
 h \, h enthalpy
 \mathbf{\Phi} \, mass flux,  = \rho \mathbf{ U }\,
 \phi \, phi mass flux through faces,  = \rho \mathbf{ U } . \mathbf{ S }_{ f } \,
 \alpha_{ t } \, alphaEff effective turbulent thermal diffusivity
 \psi \, psi compressibility  = \rho / p
 \mu_{ eff } \, muEff effective viscosity  = \mu + \mu_{ t } \,
 \mu_{ t } \, turbulent viscosity
 \mu \, mu laminar molecular viscosity
 \rho \mathbf{ R } \, divRhoR Reynolds stresses  = \rho \overline{ U_{ i }' U_{ j }' } \,

3 Momentum

The primary Navier Stokes equation is for the transport of momentum with sources and sinks.


\frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla \cdot ( \rho \mathbf{U} \mathbf{U} ) + \nabla \cdot ( \rho \mathbf{R} ) = -\nabla p + \rho \mathbf{g}

or in terms of the flux which must be discretised at cell faces (unlike most other quantities which are discretised at cell centres),


\frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla \cdot ( \mathbf{\Phi} \mathbf{U} ) + \nabla \cdot ( \rho \mathbf{R} ) = -\nabla p + \rho \mathbf{g}


File: OpenFOAM-1.6/applications/solvers/heatTransfer/buoyantPisoFoam/UEqn.H
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(phi, U)
      + turbulence->divDevRhoR(U)
    );

    UEqn.relax();

    solve(UEqn == -fvc::grad(p) - fvc::grad(rho)*gh);

We note that the mass flux is on hard-wired linear interpolation from cell centres to faces,

File: OpenFOAM-1.6.x/src/finiteVolume/cfdTools/compressible/compressibleCreatePhi.H
surfaceScalarField phi
(
    IOobject
    (
    "phi",
    runTime.timeName(),
    mesh,
    IOobject::READ_IF_PRESENT,
    IOobject::AUTO_WRITE
    ),
    linearInterpolate(rho*U) & mesh.Sf()
);

4 Mass

Conservation of mass provides a relation between the local rate of change of density and the degree to which the local mass flux is acting as a source or sink.


\frac{\partial \rho}{\partial t} + \nabla \cdot \phi = 0


File: OpenFOAM-1.6/src/finiteVolume/cfdTools/compressible/rhoEqn.H
{
    solve(fvm::ddt(rho) + fvc::div(phi));
}

5 Energy

Actually the more general concept of enthalpy is used. This is more relevant for gases and chemical reactions because it accounts for internal energy and the work done by the system, or in other words kinetic and potential energy.


\frac{\partial ( \rho h ) }{\partial t} + \nabla \cdot ( \phi h ) + \nabla \cdot ( \alpha_{t} \nabla h ) = \frac{Dp}{Dt}


File: applications/solvers/heatTransfer/buoyantPisoFoam/hEqn.H
{
    solve
    (
        fvm::ddt(rho, h)
      + fvm::div(phi, h)
      - fvm::laplacian(turbulence->alphaEff(), h) 
     ==
        DpDt
    );

    hEqn.relax();
    hEqn.solve();

    thermo.correct();
}

6 Turbulence

6.1 Momentum source

The Reynolds stresses are estimated [1] via the tensor


\begin{align}
\mathbf{R} & = \frac{2}{3} \mathbf{I} k - 2 \frac{\mu_{t} }{\rho} \mathrm{dev} \left[ \mathrm{symm} \left( \nabla \mathbf{U} \right) \right ] \\
             & = \frac{2}{3} \mathbf{I} k - 2 \frac{\mu_{t}}{\rho} \left( \mathrm{symm} \left( \nabla \mathbf{U} \right) - \frac{1}{3} \mathrm{ tr } \left( \mathrm{symm} \left( \nabla \mathbf{U} \right) \right) \mathbf{I} \right) \\
             & = \frac{2}{3} \mathbf{I} k - 2 \frac{\mu_{t}}{\rho} \left( \frac{1}{2} \left( \nabla \mathbf{U} + ( \nabla \mathbf{U} )^{T} \right) - \frac{1}{3} \mathrm{tr} \left( \nabla \mathbf{U} \right) \mathbf{I} \right) \\
\end{align}

because the trace of the symmetric component of a tensor is just the trace of the original tensor,


\mathrm{tr}( \mathrm{symm} \mathbf{T} ) = \mathrm{tr} ( T )


File: OpenFOAM-1.6/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
273 tmp<volSymmTensorField> kOmegaSST::R() const
274 {
275     return tmp<volSymmTensorField>
276     (
277         new volSymmTensorField
278         (
279             IOobject
280             (
281                 "R",
282                 runTime_.timeName(),
283                 mesh_,
284                 IOobject::NO_READ,
285                 IOobject::NO_WRITE
286             ),
287             ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
288             k_.boundaryField().types()
289         )
290     );
291 }

The divergence of the product appearing in the momentum equations is estimated as


\nabla . ( \rho \mathbf{R} ) = - \nabla . ( \mu_{eff} \nabla \mathbf{U} ) - \nabla . \left( \mu_{eff} \nabla \left( \nabla ( \mathbf{U}^{T} ) - \frac{2}{3} \mathrm{tr} ( \nabla ( \mathbf{U}^{T} ) ) \mathbf{I} \right) \right)


File: OpenFOAM-1.6/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
314 tmp<fvVectorMatrix> kOmegaSST::divDevRhoReff(volVectorField& U) const
315 {
316     return
317     (
318       - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
319     );
320 }

For reference,


\mathrm{ grad \ } \mathbf{ U } = \nabla \mathbf{ U } = \left[ \begin{array}{ccc}
    \dfrac{\partial U_{x}}{\partial x} & \dfrac{\partial U_{y}}{\partial x} & \dfrac{\partial U_{z}}{\partial x} \\
    \\
    \dfrac{\partial U_{x}}{\partial y} & \dfrac{\partial U_{y}}{\partial y} & \dfrac{\partial U_{z}}{\partial y} \\
    \\
    \dfrac{\partial U_{x}}{\partial z} & \dfrac{\partial U_{y}}{\partial z} & \dfrac{\partial U_{z}}{\partial z} \\
\end{array} \right]

7 PISO

The non linearity of the convection term of the momentum equation requires an iterative solution technique if we do not wish to be limited to small time steps. But while the momentum equation provides a way to find the velocity in terms of the pressure gradient, we need some way to determine the pressure as we iterate. In the Pressure-Implicit Split-Operator (PISO) approach, we formulate a linear relationship between the pressure gradient and the velocity in order to obtain a prediction of the velocity (which violates mass conservation), then apply the mass conservation requirement in order to obtain a pressure field that is mass conservative so that it can be fed back to obtain a corrected velocity field. At the same time, we correct for face skewness by updating the mass fluxes at cell faces with the new pressure field.

The momentum equation is discretised at N discrete locations to obtain


A \mathbf{U} - \mathbf{H}

or


\left[ \begin{array}{cccc}
    A_{11} & A_{12} & \ldots & A_{1N} \\
    A_{21} & A_{22} & \ldots & A_{2N} \\
    \vdots & \vdots & \ddots & \vdots \\
    A_{N1} & A_{N2} & \ldots & A_{NN} \\
\end{array} \right]
\left[ \begin{array}{c}
    \mathbf{U}_{1} \\
    \mathbf{U}_{2} \\
    \vdots \\
    \mathbf{U}_{N} \\
\end{array} \right] = 
\left[ \begin{array}{c}
    H_{1} \\
    H_{2} \\
    \vdots \\
    H_{N} \\
\end{array} \right]

A and \mathbf{H} are then reconstituted to make A into a diagonal matrix with the "offcuts" placed into the square matrix \mathbf{H}

Multiplying both sides by the inverse of A yields the momentum corrector equation


\mathbf{U} = \frac{\mathbf{H}}{A} - \frac{1}{A} \nabla p

Applying mass conservation, we obtain the pressure corrector equation


\begin{align}
& \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{U}) = 0 \\
& \frac{\partial \rho}{\partial t} + \mathbf{U} \cdot \nabla \rho + \rho \nabla \cdot \mathbf{U} = 0 \\
& \frac{\partial \rho}{\partial t} + \frac{\mathbf{H}}{A} \cdot \nabla \rho - \frac{1}{A} \nabla p \cdot \nabla \rho + \rho \nabla \cdot \left( \frac{\mathbf{H}}{A} \right) - \rho \nabla \cdot \left( \frac{1}{A} \nabla p \right) = 0 \\
\end{align}

This is solved to provide a pressure field inserted into a flux corrector equation


\begin{align}
\phi & = \mathbf{S}_{f} . \mathbf{U}_{f} \\
& = \mathbf{S}_{f} \cdot \left[ \left( \frac{\mathbf{H}}{A} \right)_{f} - \left( \frac{1}{A} \right)_{f} ( \nabla p )_{f} \right] \\
\end{align}

7.1 Algorithm

Step Equations Code Comments
1 Accept p, \mathbf{U}, \phi, and require \rho
while (runTime.loop())
{
...
}
The "old" values go into the loop, either initialised in createFields.H or provided by previous time step.
2 \mathrm{LHS} = \frac{\partial ( \rho \mathbf{U} ) }{\partial t} + \nabla . ( \phi \mathbf{U} ) + \nabla . ( \rho \mathbf{R} )
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(phi, U)
      + turbulence->divDevRhoR(U)
    );
Assemble the LHS of the momentum equation and discretise it into matrix form using the vector calculus schemes in fvSchemes.
3 \mathrm{RHS} = -\nabla p + \rho \mathbf{g}
    UEqn.relax();
        
    if (momentumPredictor)
    {       
        solve
        (   
            UEqn
         == 
            fvc::reconstruct
            (
                fvc::interpolate(rho)*(g & mesh.Sf())
              - fvc::snGrad(p)*mesh.magSf()
            )
        );
    }

This predictor step uses the old pressure and so while momentum conservation is satisfied, mass conservation is normally not satisfied. In some cases (e.g. low Re flow), switching off this momentum predictor via the switch in fvSolution can improve solution speed. Notice that the LHS is formulated in implicit (matrix) terms (i.e. fvm:: while the RHS terms are explicit (i.e. fvc::). The equation is therefore solved implicitely.


7.2 References

  1. Ferziger, J. H. and Peric, M., Computational Methods for Fluid Dynamics, 3rd Edition, Springer, 2002, Eqn. 9.40