Non-Reflecting Boundary Conditions in OpenFOAM

What we want to achieve

When we simulate fluid flow, we have to cut a finite computational domain out of an entire flow region. For accurate simulation, we need to let fluid and sound wave flow smoothly out of the domain through the boundary.

The reflection at the boundary has a larger effect on the solution especially when we perform a compressible flow simulation as can be seen in the following movie (Upper: With reflection, Lower: Without reflection of sound wave).

In OpenFOAM, we can use two approximate non-reflecting boundary conditions:

They determine the boundary value by solving the following equation

\frac{D \phi}{D t} = \frac{\partial \phi}{\partial t} + \boldsymbol{U} \cdot \nabla \phi = 0, \tag{1} \label{eq:advection}

where \(D/Dt\) is the material derivative and \(\boldsymbol{U}(\boldsymbol{x}, t)\) is the advection velocity.

We assume that the advection velocity \(\boldsymbol{U}\) is parallel to the boundary (face) normal direction and rewrite the eqn. \eqref{eq:advection} as

\frac{D \phi}{Dt} \approx \frac{\partial \phi}{\partial t} + U_{n} \cdot \frac{\partial \phi}{\partial \boldsymbol{n}}= 0, \tag{2} \label{eq:advection2}

where \(\boldsymbol{n}\) is the outward-pointing unit normal vector.

These boundary conditions are different in how the advection speed (scalar quantity) \(U_{n}\) is calculated and it is calculated in advectionSpeed() member function.

advective B.C.

The advection speed is the component of the velocity normal to the boundary

U_n = u_n. \tag{3} \label{eq:advectiveUn}

waveTransmissive B.C.

The advection speed is the sum of the component of the velocity normal to the boundary and the speed of sound \(c\)

U_n = u_n + c = u_n + \sqrt{\gamma/\psi}, \tag{4} \label{eq:waveTransmissiveUn}

where \(\gamma\) is the ratio of specific heats \(C_p/C_v\) and \(\psi\) is compressibility.

What do lInf and fieldInf mean?

Coming soon.

Thermal Boundary Conditions in OpenFOAM

Many thermal boundary conditions are available in OpenFOAM.
I will upload some basic cases that explain the usage of these boundary conditions.

Source Code

  • convectiveHeatTransfer

    It calculates the heat transfer coefficients from the following empirical correlations for forced convection heat transfer:
    Nu = 0.664 Re^{\frac{1}{2}} Pr^{\frac{1}{3}} \left( Re \lt 5 \times 10^5 \right) \\
    Nu = 0.037 Re^{\frac{4}{5}} Pr^{\frac{1}{3}} \left( Re \ge 5 \times 10^5 \right) \tag{1} \label{eq:NuPlate}
    where \(Nu\) is the Nusselt number, \(Re\) is the Reynolds number and \(Pr\) is the Prandtl number.

  • externalCoupledTemperature
  • externalWallHeatFluxTemperature

    This boundary condition can operate in the following two modes:
    Mode#1 Specify the heat flux \(q\)
    -k \frac{T_p – T_b}{\vert \boldsymbol{d} \vert} = q + q_r \tag{2} \label{eq:fixedHeatFlux}
    * \(k\): thermal conductivity
    * \(q_r\): radiative heat flux
    * \(T_b\): temperature on the boundary

    Mode#2 Specify the heat transfer coefficient \(h\) and the ambient temperature \(T_a\) (Fig. 1)
    -k \frac{T_p – T_b}{\vert \boldsymbol{d} \vert} = \frac{T_a – T_b}{R_{th}} + q_r \tag{3} \label{eq:fixedHeatTransferCoeff}
    * \(R_{th}\): total thermal resistance of convective and conductive heat transfer
    R_{th} = \frac{1}{h} + \sum_{i=1}^{n} \frac{l_i}{k_i} \tag{4} \label{eq:Rth}

    Figure 1
  • fixedIncidentRadiation
  • lumpedMassWallTemperature

    There is a dimensionless quantity called the Biot number, which is defined as
    Bi = \frac{l/k}{1/h} = \frac{hl}{k}, \tag{5} \label{eq:Biot}
    where \(h\) is the heat transfer coefficient, \(k\) is the thermal conductivity of a solid and \(l\) is the characteristic length of the solid. As the definition in Eq. \eqref{eq:Biot} indicates, it represents the ratio of the internal conduction resistance \(l/k\) and the external convection resistance \(1/h\). If the Biot number is small (\(Bi \ll 1\)), the solid may be treated as a simple lumped mass system of an uniform temperature. This boundary condition calculates the uniform temperature variation \(\Delta T\) on the boundary from the following equation:
    m c_p \Delta T = Q \Delta t. \tag{6} \label{eq:lumpedmass}
    * \(m\): total mass [kg]
    * \(c_p\): specific heat capacity [J/(kg.K)]
    * \(Q\): net heat flux on the boundary [W]
    * \(\Delta t\): time step [s]

  • outletMappedUniformInletHeatAddition
  • totalFlowRateAdvectiveDiffusive
  • wallHeatTransfer
  • compressible::thermalBaffle1D
  • compressible::turbulentHeatFluxTemperature
  • compressible::turbulentTemperatureCoupledBaffleMixed
  • compressible::turbulentTemperatureRadCoupledMixed
  • compressible::alphatJayatillekeWallFunction
  • compressible::alphatPhaseChangeWallFunction
  • compressible::alphatWallFunction