3.4.7. Discrete Boundary Conditions

The Dirichlet boundary conditions are applied directly in the linear solver. The flux boundary conditions are linearized and then assembled to the linear system. The flux boundary conditions are processed on a face-by-face basis. The data available with each face includes all the data on the parent element.

3.4.7.1. Symmetry, 3D Momentum

The viscous stresses can only impart a normal force at a symmetry boundary. The only other force contribution is from the pressure. The pressure is integrated over the boundary using the boundary nodal values.

The normal viscous force component is assembled to the right hand side only for the laminar equations.

The viscous stress and sub-face normal are computed at each sub-face on the element face. The integrated sub-face force is assembled to its adjacent node.

(3.1087) $F_{wi} = \mu \left( {{\partial u_i} \over {\partial x_j}}+ {{\partial u_j} \over {\partial x_i}} \right) n_j A_w$

where $n_j$ is the unit sub-face normal vector and $A_w$ is the area of the sub-face.

3.4.7.2. Outflow, 3D Mass

The mass flux at a pressure-specified outflow boundary is given by (3.837). The pressure at the face in the equation is $P_{fc}$ and is the specified value (see Figure 3.31). The interior sub-face pressure is $P_{ss}$ and is an average of nodal pressures. The fully assembled Poisson equation for pressure will have positive diagonal coefficients. Note that the form of (3.837) will contribute a positive diagonal value. The nodal pressure gradient, $G_{ij} p^{ast}_j$ , contains the influence of the specified pressure. The difference of the nodal pressure gradient and the boundary pressure gradient cancels the influence of the specified pressure in the outflow boundary condition. The specified pressure at the boundary only directly influences the momentum balance.

3.4.7.3. Outflow, 3D Momentum

The outflow boundary condition is applied to boundaries with either pressure-specified inflow or pressure-specified outflow. The viscous stresses are integrated over the boundary, but the viscous force normal to the boundary is neglected.

If the flow is entering the domain, the convected velocity is a combination of a specified tangential velocity (coflow) and a normal velocity. The normal velocity is constructed from the local nodal values.

If the flow exits the domain, the convected velocity values are interpolated from nodal velocities in the element adjacent to the boundary, similar to the interior scheme discussed in Upwind Interpolation for Convection. The convected velocities are blended from an upwind interpolation (nearest boundary node) and centered interpolation. The shape functions for the centered interpolation are taken from the interior sub-face that is directly opposite the boundary sub-face. The upwind scheme will extrapolate from the nearest node and the linear profile skew upwind scheme will interpolate to the boundary sub-face centroid.

3.4.7.4. Outflow, 3D Energy and Temperature

The outflow boundary condition is applied to boundaries with either pressure-specified inflow or pressure-specified outflow. The heat conduction is integrated over the boundary. The transport of enthalpy by mass diffusion for a multicomponent system is not yet implemented (cdm – 9/26/10).

If the flow is entering the domain, the convected enthalpy is set to a far-field reference value.

The convected enthalpy values are interpolated from nodal enthalpies in the element adjacent to the boundary, similar to the interior scheme discussed in Upwind Interpolation for Convection. The convected enthalpies are blended from an upwind interpolation (nearest boundary node) and centered interpolation. The shape functions for the centered interpolation are taken from the interior sub-face that is directly opposite the boundary sub-face. The upwind scheme will extrapolate from the nearest node and the linear profile skew upwind scheme will interpolate to the boundary sub-face centroid.

3.4.7.5. Outflow, 3D Species and Soot

The outflow boundary condition is applied to boundaries with either pressure-specified inflow or pressure-specified outflow. The mass diffusion is integrated over the boundary.

If the flow is entering the domain, the convected mass fractions are set to far-field reference values.

The convected species mass fraction values are interpolated from nodal mass fractions in the element adjacent to the boundary, similar to the interior scheme discussed in Upwind Interpolation for Convection. The convected mass fractions are blended from an upwind interpolation (nearest boundary node) and centered interpolation. The shape functions for the centered interpolation are taken from the interior sub-face that is directly opposite the boundary sub-face. The upwind scheme will extrapolate from the nearest node and the linear profile skew upwind scheme will interpolate to the boundary sub-face centroid.

3.4.7.6. Outflow, 3D Turbulent Kinetic Energy

The outflow boundary condition is applied to boundaries with either pressure-specified inflow or pressure-specified outflow. If the flow is entering the domain, the convected turbulent kinetic energy is set by one of two ways:

user specified value for turbulent kinetic energy, e.g. 0.0.,

calculated entrainment value based on user specified turbulence intensity, $T_{in}$ , and the relationship

(3.1088) $k_{ip} = {{3} \over {2}}\left(U_{ref} T_{in} \right)^2 .$

The reference velocity at the integration point, $U_{ref}$ , is determined by the current integration point mass flow rate divided by a characteristic area divided by the integration point density.

The convected turbulent kinetic energy is blended from an upwind interpolation (nearest boundary node) and centered interpolation. The shape functions for the centered interpolation are taken from the interior sub-face that is directly opposite the boundary sub-face. The upwind scheme will extrapolate from the nearest node and the linear profile skew upwind scheme will interpolate to the boundary sub-face centroid.

3.4.7.7. Outflow, 3D Turbulence Dissipation

The outflow boundary condition is applied to boundaries with either pressure-specified inflow or pressure-specified outflow. If the flow is entering the domain, the convected turbulence dissipation rate is set by one of two ways:

user specified value for turbulent dissipation rate, e.g. 0.0.,

calculated entrainment value based on user specified turbulence intensity, characteristic length and the relationship

(3.1089) $\epsilon_{ip} = {C_{\mu}}^{3/4} {{k_{ip}^{3/2}} \over {l}} ,$

where $l = 0.07L$ ; $L$ represents the user-specified characteristic length of large turbulent structures. The integration point turbulent kinetic energy is again based on the user specified turbulence intensity in conjunction with Equation (3.1088).

The convected turbulent dissipation rate is blended from an upwind interpolation (nearest boundary node) and centered interpolation. The shape functions for the centered interpolation are taken from the interior sub-face that is directly opposite the boundary sub-face. The upwind scheme will extrapolate from the nearest node and the linear profile skew upwind scheme will interpolate to the boundary sub-face centroid.

3.4.7.8. Wall, 3D Turbulent Momentum

The effect of the wall force imparted by the wall on the fluid, as outlined in Wall Functions; Momentum, is handled by the standard law of the wall formulation. To explain this procedure, consider a two dimensional element with two faces that consist of a wall boundary side set, Figure 3.32.

Fig. 3.32 Integration locations for a wall boundary.

The resulting discretization of the $i^{th}$ -component of velocity, for the boundary face that is a wall can be expressed as follows,

(3.1090) $F_{wi} = - \int \tau_{ij} n_j dS = \lambda_w A_w u_{i\|} ,$

where $A_w$ is the area, $n_j$ is the unit normal to the wall, and $\lambda_w$ is the wall shear stress factor from law of the wall,

(3.1091) $\lambda_w = { {\rho \kappa u_{\tau}} \over {ln\left(E y^+ \right)} } .$

The parallel velocity component in (3.1090) is determined by the projection of the nodal velocity onto each of the four (hex) or three (tet) subcontrol boundary faces (see (3.206)). In many respects, this procedure is very much like that of a cell-centered scheme in that the nodal velocity is assumed to act over all boundary faces. The paramount difference is the ability of one nodal velocity to be applied to a multitude of faces of potentially different orientation.

As indicated in Wall Functions; Momentum, the friction velocity at the centroid of the boundary face is determined by a nonlinear solution procedure that will now be described. The procedure begins by use of (3.195), rearranged to form the function $F$ ,

(3.1092) $F \left( u_{\tau} \right) = u_\| - {u_{\tau} \over \kappa} ln \left( {{E \rho Y_p u_{\tau}} \over {\mu}} \right) .$

The objective is to determine the value of the friction velocity such that the function, $F$ , is minimized. A Newton solve is therefore constructed that has the following standard iteration form,

(3.1093) $u_{\tau}^{k+1} = u_{\tau}^k - {{F^k} \over {F'^k}} ,$

where $F^k$ is defined by (3.1094) evaluated at the $k^{th}$ iteration level, and $F'^{k}$ is defined by

(3.1094) $F'^{k} = - { 1 \over \kappa} \left[ 1 + ln \left( { {E\rho Y_p {u_{\tau}}^k } \over {\mu} } \right) \right] .$

The procedure by which the normal distance to the wall is determined is based on the method outlined by the vertex-centered CFD code TASCflow [140]. In the procedure, the normal distance to the wall is linked to the grid by the evaluation of the normal distance from the subcontrol volume center to the boundary face. Therefore, the normal distance to the wall can be determined by the following steps:

Determination of the coordinates of the subcontrol volume center by a shape function loop over all nodes. This step in the procedure mandates a SIERRA heterogeneous (face-element) workset algorithm.

The determination of a vector, $\boldsymbol{x_i}$ , from the subcontrol volume center to the respective nodal location.

The use of the perpendicular projection operator, $P_{\perp}$ , which is defined by,

(3.1095) $P_{\perp} = - n_i n_j ,$

The determination of the normal distance by

(3.1096) $Y_p = \sqrt {{x_{\perp,1}}^2 + {x_{\perp,2}}^2 + {x_{\perp,3}}^2} .$

For convenience, the density and viscosity used in all of the above equations are nodal quantities. In other words, the physical properties are not interpolated to the centroid of the boundary face.

Once the wall shear stress factor is evaluated, it is required that the appropriate component of the velocity parallel to the boundary face is used appropriately within the respective momentum equations. As was discussed in the section on non-orthogonal momentum math models, Wall Functions; Momentum, the parallel velocity can be written in component form (see (3.207)).

3.4.7.8.1. X-Momentum

The x-momentum wall force, $F_{w1}$ is expressed as

(3.1097) $F_{w1} = -\lambda_w A_w u_{1\|} ,$

where $u_{1\|}$ is defined as

(3.1098) $u_{1\|} = \left(1-n_1^2 \right) u_{1,nd} - \left(1-n_1 n_2\right)u_{2,nd} - \left(1-n_1 n_2\right)u_{3,nd} .$

Note that the form of (3.1098) allows for an implicit treatment of the force imparted by the wall on the fluid by the factor

(3.1099) $\lambda_w A_w \left(1-n_1^2 \right) .$

3.4.7.8.2. Y-Momentum

The y-momentum wall force, $F_{w2}$ is expressed as

(3.1100) $F_{w2} = -\lambda_w A_w u_{2\|} ,$

where $u_{2\|}$ is defined as

(3.1101) $u_{2\|} = \left(1-n_2^2 \right) u_{2,nd} - \left(1-n_2 n_1\right)u_{1,nd} - \left(1-n_2 n_3\right)u_{3,nd} .$

Note that the form of (3.1101) allows for an implicit treatment of the force imparted by the wall on the fluid by the factor

(3.1102) $\lambda_w A_w \left(1-{n_2}^2 \right) .$

3.4.7.8.3. Z-Momentum

The z-momentum wall force, $F_{w3}$ is expressed as

(3.1103) $F_{w3} = -\lambda_w A_w u_{3\|} ,$

where $u_{3\|}$ is defined as

(3.1104) $u_{3\|} = \left(1-n_3^2 \right) u_{3,nd} - \left(1-n_3 n_1\right)u_{1,nd} - \left(1-n_3 n_2\right)u_{2,nd} .$

Note that the form of (3.1104) allows for an implicit treatment of the force imparted by the wall on the fluid by the factor

(3.1105) $\lambda_w A_w \left(1-n_3^2 \right) .$

3.4.7.9. Wall, 3D Turbulent Kinetic Energy

As described in Wall Functions; Turbulent Kinetic Energy, the wall boundary condition for turbulent kinetic energy can be applied in a variety of ways. In general, there are two supported methods.

The first method is specify the near-wall turbulent kinetic energy as a Dirichlet condition whose value is determined by the assumption of local equilibrium between production and dissipation of turbulent kinetic energy.

The second method is to solve a transport equation for the near wall turbulent kinetic energy whose form utilizes a modified production and dissipation term based on the assumption of local equilibrium between production and dissipation of turbulent kinetic energy. The use of a full control volume equation for the near wall turbulent kinetic energy in the presence of non-zero convection and diffusion coefficients is a violation of the very tenants of the law of the wall formulation which implicitly assumes pure shear flow behavior. Nevertheless, this method is frequently used.

The Dirichlet method consists of the determination of each integration point turbulent kinetic energy by use of the following equation,

(3.1106) $k_{ip} = { {u_{\tau}^2} \over C_{\mu}^{1/2}} ,$

The value of $u_{\tau}$ is determined by a nonlinear iteration solve of the law of the wall formulation. The integration point values are area weighted and assembled into the nodal location. The nodal value of the turbulent kinetic energy is given by the accumulated area weighed integration point turbulent kinetic energy divided by the total face area.

3.4.7.10. Wall, 3D Turbulence Dissipation

Consistent with all of literature, the near-wall value of turbulent dissipation is determined from iteration-lagged values of friction velocity,

(3.1107) $\epsilon_{ip}^{n+1} = { { {u_{\tau}}^3 } \over { \kappa Y_p } } .$

As with the implementation of the turbulent kinetic energy, the value computed in (3.1107) is area weighted and assembled to the nodal location. The Dirichlet condition is determined by the assembled quantity divided by the entire area of the boundary faces that are “owned” by the node.