ChtMultiRegionFoam

ChtMultiRegionFoam

   Solver for steady or transient fluid flow and solid heat conduction, with
   conjugate heat transfer between regions, buoyancy effects, turbulence,
   reactions and radiation modelling.

1 Equations

For each region defined as fluid, the according equation for the fluid is solved and the same is done for each solid region. The regions are coupled by a thermal boundary condition. A short description of the solver can be found also in

1.1 Equations Fluid

For each fluid region the compressible Navier Stokes equation are solved.

1.1.1 Mass conservation

The variable-density continuity equation is

$\frac{\partial \rho}{\partial t} + \frac{\partial {\rho u}_j}{\partial x_j} = 0$

(1)

The source code can be found in src/finiteVolume/cfdTools/compressible/rhoEqn.H:

 
 
{
    fvScalarMatrix rhoEqn
    (
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
        fvOptions(rho)
    );
 
    fvOptions.constrain(rhoEqn);
 
    rhoEqn.solve();
 
    fvOptions.correct(rho);
}
 

1.1.2 Momentum conservation

$\frac{ \partial (\rho {u}_i)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j u_i \right) = - \frac{\partial p_{rgh}} {\partial{x_i}} - \frac{\partial \rho g_j x_j}{\partial x_i} + \frac{\partial}{\partial x_j} \left( \tau_{ij} + \tau_{t_{ij}} \right)$

(2)

$u$ represent the velocity, $g_i$ the gravitational acceleration, $p_{rgh} = p - \rho g_j x_j$ the pressure minus the hydrostatic pressure and $\tau_{ij}$ and $\tau_{t_{ij}}$ are the viscose and turbulent stresses.

The source code can be found in Ueqn.H:

 
 
    // Solve the Momentum equation
 
    MRF.correctBoundaryVelocity(U);
 
    tmp<fvVectorMatrix> tUEqn
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence.divDevRhoReff(U)
     ==
        fvOptions(rho, U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();
 
    UEqn.relax();
 
    fvOptions.constrain(UEqn);
 
    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                )*mesh.magSf()
            )
        );
 
        fvOptions.correct(U);
        K = 0.5*magSqr(U);
    }
 
    fvOptions.correct(U);
 
 

1.1.3 Energy conservation

The energy equation can be found in: https://cfd.direct/openfoam/energy-equation/

The total energy of a fluid element can be seen as the sum of kinetic energy $k = 0.5 u_i u_i$ and internal energy $e$ . The rate of change of the kinetic energy within a fluid element is the work done on this fluid element by the viscous forces, the pressure and eternal volume forces like the gravity:

$\frac{ \partial (\rho k)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j k \right) = - \frac{\partial p u_j} {\partial{x_i} }- \rho g_j u_j + \frac{\partial}{\partial x_j}\left( \tau_{ij} u_i \right)$

(3)

The rate of change of the internal energy $e$ of a fluid element is the heat transferred to this fluid element by diffusion and turbulence $q_i + q_{ti}$ plus the heat source term $r$ plus the heat source by radiation $Rad$  :

$\frac{ \partial (\rho e)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j e \right) = - \frac{\partial q_i} {\partial{x_i} } + \rho r + Rad$

(4)

The change rate of the total energy is the sum of the above two equations:

$\frac{ \partial (\rho e)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j e \right) + \frac{ \partial (\rho k)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j k \right) = - \frac{\partial ( q_i + q_{ti}) } {\partial{x_i} } + \rho r + Rad - \frac{\partial p u_j} {\partial{x_i} }- \rho g_j u_j + \frac{\partial}{\partial x_j}\left( \tau_{ij} u_i \right)$

(5)

Instead of the internal energy $e$ there is also the option to solve the equation for the enthalpy $h = e + p/\rho$:

$\frac{ \partial (\rho h)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j h \right) + \frac{ \partial (\rho k)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j k \right) = - \frac{\partial ( q_i + q_{ti})} {\partial{x_i} } + \rho r + Rad + \frac{\partial p} {\partial{t} }- \rho g_j u_j + \frac{\partial}{\partial x_j}\left( \tau_{ij} u_i \right)$

(5)

The source code can be found in EEqn.H:

 
{
    volScalarField& he = thermo.he();
 
    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + fvm::div(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(turbulence.alphaEff(), he)
     ==
        rho*(U&g)
      + rad.Sh(thermo, he)
      + Qdot
      + fvOptions(rho, he)
    );
 
    EEqn.relax();
 
    fvOptions.constrain(EEqn);
 
    EEqn.solve();
 
    fvOptions.correct(he);
 
    thermo.correct();
    rad.correct();
 
    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
        << max(thermo.T()).value() << endl;
}

1.1.4 Species conservation

In order to account for the chemical reactions occurring between different chemical species a conservation equation for each species k has to be solved:

$\frac{ \partial (\rho Y_k)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j Y_k \right) = \frac{\partial } {\partial{x_j} }\mu_{eff}\frac{\partial Y_k}{\partial x_j} + R_k$

(6)

$R_k$ is the reaction rate of the species k.

The source code can be found in YEqn.H:

 
tmp<fv::convectionScheme<scalar>> mvConvection(nullptr);
 
if (Y.size())
{
    mvConvection = tmp<fv::convectionScheme<scalar>>
    (
        fv::convectionScheme<scalar>::New
        (
            mesh,
            fields,
            phi,
            mesh.divScheme("div(phi,Yi_h)")
        )
    );
}
 
{
    reaction.correct();
    Qdot = reaction.Qdot();
    volScalarField Yt
    (
        IOobject("Yt", runTime.timeName(), mesh),
        mesh,
        dimensionedScalar("Yt", dimless, 0)
    );
 
    forAll(Y, i)
    {
        if (i != inertIndex && composition.active(i))
        {
            volScalarField& Yi = Y[i];
 
            fvScalarMatrix YiEqn
            (
                fvm::ddt(rho, Yi)
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence.muEff(), Yi)
             ==
                reaction.R(Yi)
              + fvOptions(rho, Yi)
            );
 
            YiEqn.relax();
 
            fvOptions.constrain(YiEqn);
 
            YiEqn.solve(mesh.solver("Yi"));
 
            fvOptions.correct(Yi);
 
            Yi.max(0.0);
            Yt += Yi;
        }
    }
 
    if (Y.size())
    {
        Y[inertIndex] = scalar(1) - Yt;
        Y[inertIndex].max(0.0);
    }
}