3.4.6. Discrete Transport Equations

The discrete form of the linearized equations are presented in this section. The nonlinear solution procedure consists of repeated approximate Newton linearizations and linear solves of the discrete equation,

(3.839)A \delta \phi = b .

The matrix A is based on an approximate linearization of F from (3.838) about a predicted value \phi^\ast,

(3.840)A = \left. {1 \over {\Delta t}} - \tilde{{{\partial F} \over {\partial \phi}}} \right|^{\ast} .

The right-hand side, b, of the linearized equation represents the residual of the nonlinear equation,

(3.841)b = F \left( \phi^\ast, \phi^n \right) - {{\phi^{\ast} - \phi{^n}} \over {\Delta t}}

If the nonlinear equation is converged, the right-hand side will be zero. The linear equations are solved in delta-form. The solution vector consists of the change in the unknown rather than the new value of the unknown.

There are four solution states in the nonlinear solver algorithm. The time level n is the old time level. The state \ast represents predicted values at the new time level before the linear solve. The state \ast\ast represents the values after the linear solve. The time level (n+1) is the new time level. Within the nonlinear iteration cycle, values at the new time level (n+1) are copied to the predicted level \ast before the next iteration.

There are three stages to the assembly of the matrix that result from the linearization. The first stage is the assembly of element contributions. The elements contain control-volume sub-faces that are internal to the mesh. The second state is the assembly of flux boundary conditions. The flux boundary conditions contribute to the control-volume sub-faces on the boundary of the mesh. The flux boundary condition contributions are full element contributions because they may involve both boundary and interior nodes. The third stage is the enforcement of Dirichlet boundary conditions.

The element matrix contributions are processed by first evaluating surface integral fluxes at sub-faces and then evaluating volume integral terms at sub-volumes. The flux is evaluate at a sub-face and then added or subtracted from the two adjacent control-volumes. The sub-face area components are constructed such that the face normal direction points from the left adjacent node to the right adjacent node. Fluxes are subtracted from the left node (L) and added to the right node (R). The left and right adjacent nodes for a give sub-face number within an element are given in Tables Element variable values and differentials at control-volume faces for hexahedral elements. Face-to-edge number mapping., Element variable values and differentials at control-volume faces for tetrahedral elements. Face-to-edge number mapping., and Element variable values and differentials at control-volume faces for wedge elements. Face-to-edge number mapping..

The linearization of each transport equation can be broken into contributions from the time term, convection, diffusion, and sources.

(3.842)A = A^{t} + A^{c} + A^{d} + A^{s}

(3.843)b = b^{t} + b^{c} + b^{d} + b^{s}

The linear system is assembled on an element-by-element basis. Each element contributes and N \times N element matrix where N is the number of nodes in the element. The nodal contribution from node J for the control volume about node I is A_{I,J}. Nodal variables in the following discussion are symbolized by capital letters. Linear averages of variables at face k are

(3.844)\mu_k = \sum_J \left. N_J \right|_k \mu_J

(3.845)\kappa_k = \sum_J \left. N_J \right|_k \kappa_J

The density predictor (see Variable Density) may be used to compute the density at the new time level for the time derivative term.

The convection operator for a face i is C_{i,J} and is described in Upwind Interpolation for Convection.

Gradients of variables at face k are:

(3.846)p_x = \sum_J \left. {{\partial N_J} \over {\partial x}}\right|_k P_J \qquad p_y = \sum_J \left. {{\partial N_J} \over {\partial y}}\right|_k P_J \qquad p_z = \sum_J \left. {{\partial N_J} \over {\partial z}}\right|_k P_J

(3.847)u_x = \sum_J \left. {{\partial N_J} \over {\partial x}}\right|_k U_J \qquad u_y = \sum_J \left. {{\partial N_J} \over {\partial y}}\right|_k U_J \qquad u_z = \sum_J \left. {{\partial N_J} \over {\partial z}}\right|_k U_J

(3.848)v_x = \sum_J \left. {{\partial N_J} \over {\partial x}}\right|_k V_J \qquad v_y = \sum_J \left. {{\partial N_J} \over {\partial y}}\right|_k V_J \qquad v_z = \sum_J \left. {{\partial N_J} \over {\partial z}}\right|_k V_J

(3.849)w_x = \sum_J \left. {{\partial N_J} \over {\partial x}}\right|_k W_J \qquad w_y = \sum_J \left. {{\partial N_J} \over {\partial y}}\right|_k W_J \qquad w_z = \sum_J \left. {{\partial N_J} \over {\partial z}}\right|_k W_J

(3.850)t_x = \sum_J \left. {{\partial N_J} \over {\partial x}}\right|_k T_J \qquad t_y = \sum_J \left. {{\partial N_J} \over {\partial y}}\right|_k T_J \qquad t_z = \sum_J \left. {{\partial N_J} \over {\partial z}}\right|_k T_J

3.4.6.1. Positive-Flow Convention and Integration Quadrature

The sign on a flux integral is defined such that flow into a control volume is positive and flow out of a control volume is negative. The equations are assembled into the implicit matrix and right-hand side such that the time derivative contribution of an unknown is positive. In reference to the model differential equation, (3.838), any implicit terms that contribute to the control volume balance, F(\phi), in a positive sense must be moved to the implicit left-hand side, switching signs.

The control volume balance is assembled on an element-by-element basis. Each element contributes terms from fluxes over its internal sub-control volume faces and volumetric terms from its internal sub-control volumes. A flux is computed for each sub-control volume face. The flux contribution is then summed into the two adjacent control volumes, adjusting the sign according to whether the flux is in or out of the control volume. The convention is that the sub-face normal direction between two adjacent control volumes is positive from the lower local sub-volume number to the higher sub-volume number in a local node numbering sense. The consistent treatment of fluxes is a requirement for conservation. Each sub-control volume face is numbered the same as the element edge number. The two adjacent control volumes for each edge number are given in Tables Element variable values and differentials at control-volume faces for hexahedral elements. Face-to-edge number mapping., Element variable values and differentials at control-volume faces for tetrahedral elements. Face-to-edge number mapping., and Element variable values and differentials at control-volume faces for wedge elements. Face-to-edge number mapping. for different element types.

The elemental flux contributions are assembled into a global control volume matrix. Each control volume balance is written in terms of coefficients multiplying the surrounding nodal values. In terms of matrix terminology for two-dimensional elements, the matrix coefficient for Node 5 of Figure 3.28, associated with the control volume center, is the diagonal term and should be positive. All other nodal coefficients for the control volume balance are the off-diagonal terms and complete one row of a global flux-balance matrix.

The control volume flux integrals are evaluated using numerical quadrature. The integral term for each control volume sub-face and sub-volume is evaluated using a single quadrature point. The number of quadrature points for the surface fluxes in an element is equivalent to the number of sub-faces. For example, a quadrilateral element will have four sub-face quadratures and four sub-volume quadratures. A hexahedral element will have twelve sub-face quadratures and eight sub-volume quadratures.

In three-dimensional elements, the control-volume sub-faces may not be planar. Care must be taken to conserve surface area over a control-volume to prevent non-physical sources and sinks. The sub-faces in a three-dimensional element are defined by bilinear surfaces and the discrete surface area differential is also a bilinear function. Since the quadrature for a bilinear function is exact if evaluated at the mid-point, the current quadrature strategy will ensure surface area conservation.

The quadrature coefficients are customarily derived such that the integration ranges from -1 to 1, so a mapping is required to quadrature space.

(3.851)\int_a^b F \left( \xi \right) {\rm d}\xi = {{b - a} \over 2} \int_{-1}^1 f \left( \bar{\xi} \right) {\rm d}\bar{\xi}

(3.852)\xi = {{b + a} \over 2} + {{b - a} \over 2} \bar{\xi}

The integrand is evaluated at discrete points, called Gauss points, and summed using weighting functions.

(3.853)\int_{-1}^1 F \left( \bar{\xi} \right) {\rm d}\bar{\xi} = w_i F \left( \bar{\xi}_i \right)

For a one-point quadrature, \bar{\xi}_1 = 0 and w_1 = 2.

3.4.6.2. X-Momentum, 3D Laminar Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.854)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.855)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} U_{I}^{\ast} - \rho_{I}^{n} U_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.856)A^c_{IL,J} += C^{\ast}_{k,J}

(3.857)A^c_{IR,J} -= C^{\ast}_{k,J}

(3.858)b^c_{IL} -= \sum_J C^{\ast}_{k,J} U^{\ast}_J

(3.859)b^c_{IR} += \sum_J C^{\ast}_{k,J} U^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix. The stress term may or may not include the molecular viscosity, depending on the user specified model.

(3.860)F_{k,J} = - \mu_k \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.861)A^d_{IL,J} += F_{k,J}

(3.862)A^d_{IR,J} -= F_{k,J}

(3.863)\tau_{xx} = \mu_k \left( u_x^{\ast} + u_x^{\ast} \right)

(3.864)\tau_{xy} = \mu_k \left( u_y^{\ast} + v_x^{\ast} \right)

(3.865)\tau_{xz} = \mu_k \left( u_z^{\ast} + w_x^{\ast} \right)

(3.866)f_k = - \left( \tau_{xx} A_x + \tau_{xy} A_y + \tau_{xz} A_z \right)

(3.867)b^d_{IL} -= f_k

(3.868)b^d_{IR} += f_k

The pressure is assembled in the form of a volume integral. The pressure gradients have been pre-computed at nodes use a surface-integral approximation.

(3.869)b^s_{I} -= \left. {{\partial P} \over {\partial x}}\right|_I^\ast \Delta V_I

3.4.6.3. Y-Momentum, 3D Laminar Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.870)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.871)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} V_{I}^{\ast} - \rho_{I}^{n} V_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.872)A^c_{IL,J} += C^{\ast}_{k,J}

(3.873)A^c_{IR,J} -= C^{\ast}_{k,J}

(3.874)b^c_{IL} -= \sum_J C^{\ast}_{k,J} V^{\ast}_J

(3.875)b^c_{IR} += \sum_J C^{\ast}_{k,J} V^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix.

(3.876)F_{k,J} = - \mu_k \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.877)A^d_{IL,J} += F_{k,J}

(3.878)A^d_{IR,J} -= F_{k,J}

(3.879)\tau_{yx} = \mu_k \left( v_x^{\ast} + u_y^{\ast} \right)

(3.880)\tau_{yy} = \mu_k \left( v_y^{\ast} + v_y^{\ast} \right)

(3.881)\tau_{yz} = \mu_k \left( v_z^{\ast} + w_y^{\ast} \right)

(3.882)f_k = - \left( \tau_{yx} A_x + \tau_{yy} A_y + \tau_{yz} A_z \right)

(3.883)b^d_{IL} -= f_k

(3.884)b^d_{IR} += f_k

The pressure is assembled in the form of a volume integral. The pressure gradients have been pre-computed at nodes use a surface-integral approximation.

(3.885)b^s_{I} -= \left. {{\partial P} \over {\partial y}}\right|_I^\ast \Delta V_I

3.4.6.4. Z-Momentum, 3D Laminar Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.886)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.887)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} W_{I}^{\ast} - \rho_{I}^{n} W_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.888)A^c_{IL,J} += C^{\ast}_{k,J}

(3.889)A^c_{IR,J} -= C^{\ast}_{k,J}

(3.890)b^c_{IL} -= \sum_J C^{\ast}_{k,J} W^{\ast}_J

(3.891)b^c_{IR} += \sum_J C^{\ast}_{k,J} W^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix.

(3.892)F_{k,J} = - \mu_k \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.893)A^d_{IL,J} += F_{k,J}

(3.894)A^d_{IR,J} -= F_{k,J}

(3.895)\tau_{zx} = \mu_k \left( w_x^{\ast} + u_z^{\ast} \right)

(3.896)\tau_{zy} = \mu_k \left( w_y^{\ast} + v_z^{\ast} \right)

(3.897)\tau_{zz} = \mu_k \left( w_z^{\ast} + w_z^{\ast} \right)

(3.898)f_k = - \left( \tau_{zx} A_x + \tau_{zy} A_y + \tau_{zz} A_z \right)

(3.899)b^d_{IL} -= f_k

(3.900)b^d_{IR} += f_k

The pressure is assembled in the form of a volume integral. The pressure gradients have been pre-computed at nodes use a surface-integral approximation.

(3.901)b^s_{I} -= \left. {{\partial P} \over {\partial z}}\right|_I^\ast \Delta V_I

3.4.6.5. Buoyancy, Momentum Transport

The body force imposed by the buoyancy term can be constructed in one of three ways.

3.4.6.5.1. Boussinesq Form

For the Boussinesq approximation, the body force is evaluated at the sub-volume centroid, k, for sub-volume I.

(3.902)b^s_{I} -= {{\rho g} \over {T_\circ}} \left( \sum_J \left. N_J \right|_k T_J - T_\circ \right) \Delta V_I

3.4.6.5.2. Differential Form

For the “differential” form, the hydrostatic component of pressure has been removed. The body force is evaluated at the control-volume centroid, for sub-volume I.

(3.903)b^s_{I} += \left( \rho_I^\ast - \rho_\circ \right) g \Delta V_I

3.4.6.5.3. Full Form

The body force is evaluated at the control-volume centroid, for sub-volume I.

(3.904)b^s_{I} += \rho_I^\ast g \Delta V_I

3.4.6.6. Mass Transport – 3D Continuity

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.905)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} - \rho_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes using the Rhie/Chow scheme from Flow Solver.

(3.906)F_{k,J} = - f \Delta t \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.907)A^d_{IL,J} += F_{k,J}

(3.908)A^d_{IR,J} -= F_{k,J}

(3.909)u_k^{\ast} = \sum_J \left. N_J \right|_k U_J^{\ast} + f {{\Delta t} \over {\rho}} \left( \sum_J \left. {{\partial P} \over {\partial x}} \right|_J^{\ast} - p^{\ast}_x \right) + f \left( u_k^{n} - \sum_J \left. N_J \right|_J U_J^{n} \right)

(3.910)v_k^{\ast} = \sum_J \left. N_J \right|_k V_J^{\ast} + f {{\Delta t} \over {\rho}} \left( \sum_J \left. {{\partial P} \over {\partial y}} \right|_J^{\ast} - p^{\ast}_y \right) + f \left( v_k^{n} - \sum_J \left. N_J \right|_J V_J^{n} \right)

(3.911)w_k^{\ast} = \sum_J \left. N_J \right|_k W_J^{\ast} + f {{\Delta t} \over {\rho}} \left( \sum_J \left. {{\partial P} \over {\partial z}} \right|_J^{\ast} - p^{\ast}_z \right) + f \left( w_k^{n} - \sum_J \left. N_J \right|_J W_J^{n} \right)

(3.912)\dot{m}_k = \rho \left( u_k^{\ast} A_x + v_k^{\ast} A_y + w_k^{\ast} A_z \right)

(3.913)b^c_{IL} -= \dot{m}_k

(3.914)b^c_{IR} += \dot{m}_k

Velocity correction and new mass flow rate…..

3.4.6.7. Energy, 3D Laminar Transport

The laminar energy equation is linearized with respect to the temperature. The time term is lumped. The time-term contribution is evaluated for each sub-volume. The density must also be linearized for stability.

(3.915)A^t_{I,I}\ +=\ \left( \rho_{I}^{\ast} C_{p,I}^{\ast} - \rho_{I}^{\ast} {{H^{\ast}_I} \over {T^{\ast}_I}} \right) {{\Delta V_I} \over {\Delta t}}

(3.916)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} H_{I}^{\ast} - \rho_{I}^{n} H_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.917)A^c_{IL,J} += C^{n+1}_{k,J} C_{p,J}^{\ast}

(3.918)A^c_{IR,J} -= C^{n+1}_{k,J} C_{p,J}^{\ast}

(3.919)b^c_{IL} -= \sum_J C^{n+1}_{k,J} H^{\ast}_J

(3.920)b^c_{IR} += \sum_J C^{n+1}_{k,J} H^{\ast}_J

The heat conduction term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.921)F_{k,J} = - \kappa_k \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.922)A^d_{IL,J} += F_{k,J}

(3.923)A^d_{IR,J} -= F_{k,J}

(3.924)q_k = - \kappa_k \left( t_x^{\ast} A_x + t_y^{\ast} A_y + t_z^{\ast} A_z \right)

(3.925)b^d_{IL} -= q_k

(3.926)b^d_{IR} += q_k

3.4.6.8. Temperature, 3D Laminar Transport

The laminar temperature equation is linearized with respect to the temperature. The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.927)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.928)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} T_{I}^{\ast} - \rho_{I}^{n} T_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.929)A^c_{IL,J} += C^{n+1}_{k,J}

(3.930)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.931)b^c_{IL} -= \sum_J C^{n+1}_{k,J} T^{\ast}_J

(3.932)b^c_{IR} += \sum_J C^{n+1}_{k,J} T^{\ast}_J

The heat conduction term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.933)F_{k,J} = - {{\kappa_k} \over {C_{p,k}}} \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.934)A^d_{IL,J} += F_{k,J}

(3.935)A^d_{IR,J} -= F_{k,J}

(3.936)q_k = - {{\kappa_k} \over {C_{p,k}}} \left( t_x^{\ast} A_x + t_y^{\ast} A_y + t_z^{\ast} A_z \right)

(3.937)b^d_{IL} -= q_k

(3.938)b^d_{IR} += q_k

A correction for variable specific heat is applied as a volume term. The correction is computed at the centroid of the sub-volume, k, for control volume I.

(3.939)b^d_{I} += {{\kappa} \over {C_p^2}} \left( t_x C_{p,x} + t_y C_{p,y} + t_z C_{p,z} \right) \Delta V_I

3.4.6.9. Species, 3D Laminar Transport

There is a species equations for each species. The mass fraction is Y_s, where s is the species number. The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.940)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.941)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} Y_{s,I}^{\ast} - \rho_{I}^{n} Y_{s,I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.942)A^c_{IL,J} += C^{n+1}_{k,J}

(3.943)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.944)b^c_{IL} -= \sum_J C^{n+1}_{k,J} Y^{\ast}_{s,J}

(3.945)b^c_{IR} += \sum_J C^{n+1}_{k,J} Y^{\ast}_{s,J}

The mass diffusion term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.946)F_{k,J} = - \rho_k D_{s,k} \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.947)A^d_{IL,J} += F_{k,J}

(3.948)A^d_{IR,J} -= F_{k,J}

(3.949)f_k = - \rho_k D_{s,k} \left( ys_x^{\ast} A_x + ys_y^{\ast} A_y + ys_z^{\ast} A_z \right)

(3.950)b^d_{IL} -= f_k

(3.951)b^d_{IR} += f_k

3.4.6.10. X-Momentum, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.952)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.953)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} U_{I}^{\ast} - \rho_{I}^{n} U_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.954)A^c_{IL,J} += C^{\ast}_{k,J}

(3.955)A^c_{IR,J} -= C^{\ast}_{k,J}

(3.956)b^c_{IL} -= \sum_J C^{\ast}_{k,J} U^{\ast}_J

(3.957)b^c_{IR} += \sum_J C^{\ast}_{k,J} U^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix.

(3.958)F_{k,J} = - \left( \mu_k + \mu_{T,k} \right) \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.959)A^d_{IL,J} += F_{k,J}

(3.960)A^d_{IR,J} -= F_{k,J}

(3.961)\tau_{xx} = \left( \mu_k + \mu_{T,k} \right) \left( u_x^{\ast} + u_x^{\ast} \right)

(3.962)\tau_{xy} = \left( \mu_k + \mu_{T,k} \right) \left( u_y^{\ast} + v_x^{\ast} \right)

(3.963)\tau_{xz} = \left( \mu_k + \mu_{T,k} \right) \left( u_z^{\ast} + w_x^{\ast} \right)

(3.964)f_k = - \left( \tau_{xx} A_x + \tau_{xy} A_y + \tau_{xz} A_z \right)

(3.965)b^d_{IL} -= f_k

(3.966)b^d_{IR} += f_k

The pressure is assembled in the form of a volume integral. The pressure gradients have been pre-computed at nodes use a surface-integral approximation.

(3.967)b^s_{I} -= \left. {{\partial P} \over {\partial x}}\right|_I^\ast \Delta V_I

3.4.6.11. Y-Momentum, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.968)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.969)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} V_{I}^{\ast} - \rho_{I}^{n} V_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.970)A^c_{IL,J} += C^{\ast}_{k,J}

(3.971)A^c_{IR,J} -= C^{\ast}_{k,J}

(3.972)b^c_{IL} -= \sum_J C^{\ast}_{k,J} V^{\ast}_J

(3.973)b^c_{IR} += \sum_J C^{\ast}_{k,J} V^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix.

(3.974)F_{k,J} = - \left( \mu_k + \mu_{T,k} \right) \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.975)A^d_{IL,J} += F_{k,J}

(3.976)A^d_{IR,J} -= F_{k,J}

(3.977)\tau_{yx} = \left( \mu_k + \mu_{T,k} \right) \left( v_x^{\ast} + u_y^{\ast} \right)

(3.978)\tau_{yy} = \left( \mu_k + \mu_{T,k} \right) \left( v_y^{\ast} + v_y^{\ast} \right)

(3.979)\tau_{yz} = \left( \mu_k + \mu_{T,k} \right) \left( v_z^{\ast} + w_y^{\ast} \right)

(3.980)f_k = - \left( \tau_{yx} A_x + \tau_{yy} A_y + \tau_{yz} A_z \right)

(3.981)b^d_{IL} -= f_k

(3.982)b^d_{IR} += f_k

The pressure is assembled in the form of a volume integral. The pressure gradients have been pre-computed at nodes use a surface-integral approximation.

(3.983)b^s_{I} -= \left. {{\partial P} \over {\partial y}}\right|_I^\ast \Delta V_I

3.4.6.12. Z-Momentum, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.984)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.985)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} W_{I}^{\ast} - \rho_{I}^{n} W_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.986)A^c_{IL,J} += C^{\ast}_{k,J}

(3.987)A^c_{IR,J} -= C^{\ast}_{k,J}

(3.988)b^c_{IL} -= \sum_J C^{\ast}_{k,J} W^{\ast}_J

(3.989)b^c_{IR} += \sum_J C^{\ast}_{k,J} W^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix.

(3.990)F_{k,J} = - \left( \mu_k + \mu_{T,k} \right) \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.991)A^d_{IL,J} += F_{k,J}

(3.992)A^d_{IR,J} -= F_{k,J}

(3.993)\tau_{zx} = \left( \mu_k + \mu_{T,k} \right) \left( w_x^{\ast} + u_z^{\ast} \right)

(3.994)\tau_{zy} = \left( \mu_k + \mu_{T,k} \right) \left( w_y^{\ast} + v_z^{\ast} \right)

(3.995)\tau_{zz} = \left( \mu_k + \mu_{T,k} \right) \left( w_z^{\ast} + w_z^{\ast} \right)

(3.996)f_k = - \left( \tau_{zx} A_x + \tau_{zy} A_y + \tau_{zz} A_z \right)

(3.997)b^d_{IL} -= f_k

(3.998)b^d_{IR} += f_k

The pressure is assembled in the form of a volume integral. The pressure gradients have been pre-computed at nodes use a surface-integral approximation.

(3.999)b^s_{I} -= \left. {{\partial P} \over {\partial z}}\right|_I^\ast \Delta V_I

3.4.6.13. Turbulent Kinetic Energy, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.1000)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.1001)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} K_{I}^{\ast} - \rho_{I}^{n} K_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1002)A^c_{IL,J} += C^{n+1}_{k,J}

(3.1003)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.1004)b^c_{IL} -= \sum_J C^{n+1}_{k,J} K^{\ast}_J

(3.1005)b^c_{IR} += \sum_J C^{n+1}_{k,J} K^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix.

(3.1006)F_{k,J} = - {{\mu_{T,k}} \over {\sigma_k}} \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.1007)A^d_{IL,J} += F_{k,J}

(3.1008)A^d_{IR,J} -= F_{k,J}

(3.1009)f_k = - {{\mu_{T,k}} \over {\sigma_k}} \left( k_x^{\ast} A_x + k_y^{\ast} A_y + k_z^{\ast} A_z \right)

(3.1010)b^d_{IL} -= f_k

(3.1011)b^d_{IR} += f_k

The turbulence production is assembled in the form of a volume integral. The velocity derivatives are computed at the sub-volume centroids.

(3.1012)\Phi = 2 \left( u_x^2 + v_y^2 + w_z^2 \right) - {2 \over 3} \left( u_x + v_y + w_z \right)^2

(3.1013)+ \left( u_y + v_x \right)^2 + \left( v_z + w_y \right)^2 + \left( w_x + u_z \right)^2

(3.1014)b^s_{I} += \mu_{T} \Phi \Delta V_I

The turbulence dissipation is assembled in the form of a volume integral. The terms are evaluated at the node associated with the control volume.

(3.1015)A^s_{I,I}\ +=\ \rho_I {{E_I^{\ast}} \over {K_I^{\ast}}} \Delta V_I

(3.1016)b^s_{I}\ -=\ \rho_I E_I^{\ast} \Delta V_I

3.4.6.14. Turbulence Dissipation, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.1017)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.1018)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} E_{I}^{\ast} - \rho_{I}^{n} E_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1019)A^c_{IL,J} += C^{n+1}_{k,J}

(3.1020)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.1021)b^c_{IL} -= \sum_J C^{n+1}_{k,J} E^{\ast}_J

(3.1022)b^c_{IR} += \sum_J C^{n+1}_{k,J} E^{\ast}_J

The viscous stress term is computed at each face k and assembled to the left (IL) and right (IR) control volumes. Only the solenoidal part of the stress term is used for the matrix. As with the turbulent kinetic energy transport equation, the molecular viscosity may augment the effective diffusivity.

(3.1023)F_{k,J} = - {{\mu_{T,k}} \over {\sigma_\epsilon}} \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.1024)A^d_{IL,J} += F_{k,J}

(3.1025)A^d_{IR,J} -= F_{k,J}

(3.1026)f_k = - {{\mu_{T,k}} \over {\sigma_\epsilon}} \left( \epsilon_x^{\ast} A_x + \epsilon_y^{\ast} A_y + \epsilon_z^{\ast} A_z \right)

(3.1027)b^d_{IL} -= f_k

(3.1028)b^d_{IR} += f_k

The velocity derivatives are computed at the sub-volume centroids using velocities at the new time level (n+1).

(3.1029)\Phi = 2 \left( u_x^2 + v_y^2 + w_z^2 \right) - {2 \over 3} \left( u_x + v_y + w_z \right)^2

(3.1030)+ \left( u_y + v_x \right)^2 + \left( v_z + w_y \right)^2 + \left( w_x + u_z \right)^2

(3.1031)b^s_{I} += \mu_{T} C_{\epsilon_1} \Phi {{E_I^{\ast}} \over {K_I^{\ast}}} \Delta V_I

The turbulence dissipation is assembled in the form of a volume integral. The terms are evaluated at the node associated with the control volume.

(3.1032)A^s_{I,I}\ +=\ \rho_I C_{\epsilon_2} {{E_I^{\ast}} \over {K_I^{\ast}}} \Delta V_I

(3.1033)b^s_{I}\ -=\ \rho_I C_{\epsilon_2} {{E_I^{\ast}} \over {K_I^{\ast}}} E_I^{\ast} \Delta V_I

3.4.6.15. Energy, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.1034)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.1035)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} H_{I}^{\ast} - \rho_{I}^{n} H_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1036)A^c_{IL,J} += C^{n+1}_{k,J}

(3.1037)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.1038)b^c_{IL} -= \sum_J C^{n+1}_{k,J} H^{\ast}_J

(3.1039)b^c_{IR} += \sum_J C^{n+1}_{k,J} H^{\ast}_J

The heat conduction term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1040)F_{k,J} = - \left( {{\mu_k} \over {\rm Pr}} + {{\mu_{T,k}} \over {{\rm Pr}_T}} \right) \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.1041)A^d_{IL,J} += F_{k,J}

(3.1042)A^d_{IR,J} -= F_{k,J}

(3.1043)q_k = - \left( {{\mu_k} \over {\rm Pr}} + {{\mu_{T,k}} \over {{\rm Pr}_T}} \right) \left( h_x^{\ast} A_x + h_y^{\ast} A_y + h_z^{\ast} A_z \right)

(3.1044)b^d_{IL} -= q_k

(3.1045)b^d_{IR} += q_k

3.4.6.16. Species, 3D Turbulent Transport

There is a species equations for each species. The mass fraction is Y_s, where s is the species number. The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.1046)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.1047)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} Y_{s,I}^{\ast} - \rho_{I}^{n} Y_{s,I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1048)A^c_{IL,J} += C^{n+1}_{k,J}

(3.1049)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.1050)b^c_{IL} -= \sum_J C^{n+1}_{k,J} Y^{\ast}_{s,J}

(3.1051)b^c_{IR} += \sum_J C^{n+1}_{k,J} Y^{\ast}_{s,J}

The mass diffusion term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1052)F_{k,J} = - \left( {{\mu_k} \over {\rm Sc}} + {{\mu_{T,k}} \over {{\rm Sc}_T}} \right) \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.1053)A^d_{IL,J} += F_{k,J}

(3.1054)A^d_{IR,J} -= F_{k,J}

(3.1055)f_k = - \left( {{\mu_k} \over {\rm Sc}} + {{\mu_{T,k}} \over {{\rm Sc}_T}} \right) \left( ys_x^{\ast} A_x + ys_y^{\ast} A_y + ys_z^{\ast} A_z \right)

(3.1056)b^d_{IL} -= f_k

(3.1057)b^d_{IR} += f_k

The chemical production source terms from the EDC model are applied at the centroid of the control volume. The production term is constructed from the rate, the fine structure mass fractions, and the average mass fractions.

(3.1058)A^s_{I,I}\ +=\ \dot{r}_{s,I} \Delta V_I

(3.1059)\dot{\omega}_{s,I} = \dot{r}_{s,I} \left( Y_{s,I}^{fs} - Y_{s,I} \right)

(3.1060)b^s_{I} += \dot{\omega}_{s,I} \Delta V_I

3.4.6.17. Soot Transport, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.1061)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.1062)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} S_{I}^{\ast} - \rho_{I}^{n} S_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1063)A^c_{IL,J} += C^{n+1}_{k,J}

(3.1064)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.1065)b^c_{IL} -= \sum_J C^{n+1}_{k,J} S^{\ast}_{J}

(3.1066)b^c_{IR} += \sum_J C^{n+1}_{k,J} S^{\ast}_{J}

The diffusion term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1067)F_{k,J} = - \left( {{\mu_k} \over {\rm Sc}} + {{\mu_{T,k}} \over {{\rm Sc}_T}} \right) \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.1068)A^d_{IL,J} += F_{k,J}

(3.1069)A^d_{IR,J} -= F_{k,J}

(3.1070)f_k = - \left( {{\mu_k} \over {\rm Sc}} + {{\mu_{T,k}} \over {{\rm Sc}_T}} \right) \left( s_x^{\ast} A_x + s_y^{\ast} A_y + s_z^{\ast} A_z \right)

(3.1071)b^d_{IL} -= f_k

(3.1072)b^d_{IR} += f_k

The soot production source term from the EDC model is applied at the centroid of the control volume.

(3.1073)b^s_{I} += \dot{\omega}_{soot,I} \Delta V_I

3.4.6.18. Soot Nuclei Transport, 3D Turbulent Transport

The time term is lumped. The time-term contribution is evaluated for each sub-volume.

(3.1074)A^t_{I,I}\ +=\ \rho_{I}^{\ast} {{\Delta V_I} \over {\Delta t}}

(3.1075)b^t_{I}\ -=\ \left( \rho_{I}^{\ast} N_{I}^{\ast} - \rho_{I}^{n} N_{I}^{n} \right) {{\Delta V_I} \over {\Delta t}}

The convection term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1076)A^c_{IL,J} += C^{n+1}_{k,J}

(3.1077)A^c_{IR,J} -= C^{n+1}_{k,J}

(3.1078)b^c_{IL} -= \sum_J C^{n+1}_{k,J} N^{\ast}_{J}

(3.1079)b^c_{IR} += \sum_J C^{n+1}_{k,J} N^{\ast}_{J}

The diffusion term is computed at each face k and assembled to the left (IL) and right (IR) control volumes.

(3.1080)F_{k,J} = - \left( {{\mu_k} \over {\rm Sc}} + {{\mu_{T,k}} \over {{\rm Sc}_T}} \right) \left( \left. {{\partial N_J} \over {\partial x}}\right|_k A_x + \left. {{\partial N_J} \over {\partial y}}\right|_k A_y + \left. {{\partial N_J} \over {\partial z}}\right|_k A_z \right)

(3.1081)A^d_{IL,J} += F_{k,J}

(3.1082)A^d_{IR,J} -= F_{k,J}

(3.1083)f_k = - \left( {{\mu_k} \over {\rm Sc}} + {{\mu_{T,k}} \over {{\rm Sc}_T}} \right) \left( n_x^{\ast} A_x + n_y^{\ast} A_y + n_z^{\ast} A_z \right)

(3.1084)b^d_{IL} -= f_k

(3.1085)b^d_{IR} += f_k

The soot nuclei production source term from the EDC model is applied at the centroid of the control volume.

(3.1086)b^s_{I} += \dot{\omega}_{nucl,I} \Delta V_I