5.4 Discretisation of the Viscous Fluxes 黏性通量的离散[cfd-5-4]

为了计算式(5.2)中的扩散通量 $\vec{F}_v$ ï¼必须知道控制体面上的流动量及其一阶导数。为了得到一致的空间离散并简化数据结构ï¼黏性通量的控制体宜选取成与对流通量的相同。由于黏性通量的椭圆性质ï¼速度分量 $$(u, v, w)$$ 、动力黏度 $\mu$ 以及热传导系数 $$k$$ 的值——即计算黏性项(2.23)、(2.24)与应力(2.15)所需要的量——直接在面上取平均。于是ï¼对单元中心格式(图5.13a)ï¼控制体面 $$IJ$$ 上的值由下式给出en

In order to evaluate the diffusive fluxes $\vec{F}_v$ in Eq. (5.2), flow quantities and their first derivatives have to be known at the faces of the control volumes. The control volume for the viscous fluxes is conveniently chosen to be the same as for the convective fluxes in order to obtain a consistent spatial discretisation and to simplify the data structure. Because of the elliptic nature of the viscous fluxes, values of the velocity components $$(u, v, w)$$ , the dynamic viscosity $\mu$ , and of the heat conduction coefficient $$k$$ , which are required for the computation of the viscous terms (2.23), (2.24) and of the stresses (2.15), are simply averaged at a face. Thus, in the case of the cell-centred scheme (Fig. 5.13a), the values at the face $$IJ$$ of the control volume result from

U_{IJ} = \frac{1}{2}\left(U_I + U_J\right), \tag{5.70}

其中 $$U$$ 为上述任一流动变量。中点对偶格式在面 $$ij$$ 处也有类似表达式——见图5.13b。en

where $$U$$ is any of the above flow variables. A similar expression holds in the case of the median-dual scheme for the face $$ij$$ - see Fig. 5.13b.

剩下的任务是计算式(2.15)中速度分量的一阶导数(梯度)以及式(2.24)中温度的梯度。这可以用以下两种方法之一完成:

基于单元的梯度;或
梯度平均。

The remaining task is the evaluation of the first derivatives (gradients) of the velocity components in Eq. (2.15) and of temperature in Eq. (2.24). This can be accomplished in one of two ways, i.e., by using

element-based gradients, or
average of gradients.

下面我们将进一步了解这两种方法。en

In the following, we shall learn more about both approaches.

5.4.1 Element-Based Gradients 基于单元的梯度[cfd-5-4-1]

这类梯度计算的共同特点是:必须存储关于网格单元的信息ï¼或存储与单元几何有关的一些系数。因此ï¼我们必须把数据结构扩展到前面为对流通量介绍的基于面/边的表述之外。下面讨论用于单元中心离散与中点对偶离散的三种成熟方法。en

A common feature of this type of gradient computation is the necessity to store either information about the grid elements or some coefficients related to the geometry of the elements. Hence, we have to extend the data structure beyond the face-/edge-based formulation presented earlier for the convective fluxes. Below, we discuss three well-established methods for the cell-centred and the median-dual discretisation.

Face-Centred Control Volume 面心控制体

计算控制体面上梯度的一种可能办法ï¼是定义一个以该面为中心的辅助控制体ï¼并使用Green-Gauss定理。我们已在4.4节的结构网格有限体积离散框架内讨论过这一方法。例如ï¼对中点对偶格式ï¼可以把边中点处的梯度计算为共享该边的所有单元梯度的体积平均[70]。基于单元的梯度按式(5.49)计算:遍历所有网格单元ï¼并把梯度累加到各条边上。单元面上的 $$U$$ 值通过平均节点值得到ï¼方式与式(5.53)类似。这种方法在内存与运算次数方面代价相对较高ï¼但可用于任意单元组合。en

One possible way of evaluating the gradients at a face of the control volume is to define an auxiliary control volume centred at the face and to employ the Green-Gauss theorem. We already discussed this approach in Section 4.4 in the framework of the structured finite volume discretisation. For example, in the case of the median-dual scheme we can compute the gradient at the edge-midpoint as the volume average of gradients for all elements which share the edge [70]. The element-based gradients are evaluated according to Eq. (5.49) by looping over all grid cells and accumulating the gradients at the edges. The values of $$U$$ at the cell-faces are obtained by averaging the nodal values, in a manner similar to Eq. (5.53). This approach is relatively costly in terms of the memory and the number of operations. However, it can be implemented for any mix of grid elements.

Approximate Galerkin Finite Element Approach 近似Galerkin有限元方法

另一种适用于中点对偶格式的方法由Galerkin有限元方法导出[31]。大体上说ï¼该方法把梯度在控制体表面上的积分转化为在中心节点处计算Hessian矩阵(二阶导数)。此时黏性项遵循笛卡尔坐标系下Navier-Stokes方程的微分形式(式(A.4)ï¼其中 $\xi = x$ 、 $\eta = y$ 、 $\zeta = z$ 且 $J^{-1} = 1$ )ï¼其中包含如下形式的项en

Another methodology, which is applicable to the median-dual scheme, was derived from the Galerkin finite element method [31]. Basically speaking, the approach transforms the integration of gradients over the surface of the control volume into an evaluation of the Hessian matrix (second derivatives) at the central node. The viscous terms then follow the differential form of the Navier-Stokes equations in Cartesian coordinates (Eq. (A.4) with $\xi = x$ , $\eta = y$ , $\zeta = z$ and $J^{-1} = 1$ ), which contains terms such as

\partial_x\left(\mu\,\partial_x U\right)\,, \text{ etc.} \tag{1}

其中 $\partial_m(\cdot) = \partial(\cdot)/\partial m$ 。这样ï¼黏性通量就无需再沿控制体的面进行积分。en

with $\partial_m(\cdot) = \partial(\cdot)/\partial m$ . Hence, no further integration of the viscous fluxes over the faces of the control volume is required.

原始格式是针对纯三角形/四面体网格建立的。它采用包含相应节点的所有单元的并集。为了简化实现ï¼动力黏度系数由节点值平均得到ï¼这是与Galerkin方法的一个不同之处。于是ï¼二阶导数可在节点 $$i$$ 处按如下方式计算[71]en

The original scheme was formulated for purely triangular/tetrahedral grids. It employs a union of all elements that contain the particular node. In order to simplify the implementation, the dynamic viscosity coefficient is averaged from the nodal values, which is a difference to the Galerkin method. Then, the second derivatives can be evaluated at node $$i$$ as follows [71]

\begin{aligned} &\begin{bmatrix} \partial_x(\mu\,\partial_x U) & \partial_y(\mu\,\partial_x U) & \partial_z(\mu\,\partial_x U)\\ \partial_x(\mu\,\partial_y U) & \partial_y(\mu\,\partial_y U) & \partial_z(\mu\,\partial_y U)\\ \partial_x(\mu\,\partial_z U) & \partial_y(\mu\,\partial_z U) & \partial_z(\mu\,\partial_z U) \end{bmatrix}_{i}\\ &= \frac{1}{\Omega'}\sum_{j=1}^{N_A}\left\{\begin{bmatrix} \alpha_{xx} & \alpha_{xy} & \alpha_{xz}\\ \alpha_{yx} & \alpha_{yy} & \alpha_{yz}\\ \alpha_{zx} & \alpha_{zy} & \alpha_{zz} \end{bmatrix}_{ij}\frac{\mu_i + \mu_j}{2}\left(U_i - U_j\right)\right\}. \end{aligned} \tag{5.71}

体积 $\Omega'$ 包含共享节点 $$i$$ 的所有四面体。式(5.71)中的系数矩阵 $\alpha$ 关于对角线对称[71]ï¼即 $\alpha_{xy} = \alpha_{yx}$ 、 $\alpha_{xz} = \alpha_{zx}$ 、 $\alpha_{yz} = \alpha_{zy}$ 。因此ï¼每条边只需存储六个系数。系数由下式给出[71]en

The volume $\Omega'$ contains all tetrahedra which share the node $$i$$ . The coefficient matrix $\alpha$ in Eq. (5.71) is symmetric about the diagonal [71], i.e., $\alpha_{xy} = \alpha_{yx}$ , $\alpha_{xz} = \alpha_{zx}$ , and $\alpha_{yz} = \alpha_{zy}$ , respectively. Thus, it is necessary to store only six coefficients for each edge. The coefficients are given by [71]

\alpha_{nk} = \sum_{e}\frac{\left(\vec{S}_e^{i}\right)_n\left(\vec{S}_e^{j}\right)_k}{\Omega_e}\,, \tag{5.72}

其中 $$n$$ 、 $$k$$ 表示 $$x$$ 、 $$y$$ 、 $$z$$ 下标, $(\vec{S}_e^{i})_k$ 与 $(\vec{S}_e^{j})_n$ 表示图5.18所示外侧面向量 $\vec{S}_e^{i}$ 、 $\vec{S}_e^{j}$ 的分量。求和遍及共享边 $$ij$$ 的所有四面体(各自体积为 $\Omega_e$ )。若网格静止ï¼系数可在预处理阶段计算。一个可取之处是:可以使用与对流通量相同的基于边的数据结构。en

where $$n$$ , $$k$$ denote the $$x$$ , $$y$$ , $$z$$ subscripts, and $(\vec{S}_e^{i})_k$ , $(\vec{S}_e^{j})_n$ represent components of the outer face vectors $\vec{S}_e^{i}$ , $\vec{S}_e^{j}$ displayed in Fig. 5.18. The summation is carried out over all tetrahedra (with particular volumes $\Omega_e$ ) which share the edge $$ij$$ . If the grid is stationary, the coefficients can be computed in a pre-processing step. A desirable feature is that the same edge-based data structure can be employed as for the convective fluxes.

不过这种方法的缺点在于:完全的黏性项只在三角形或四面体网格上得以保留。对棱柱或六面体之类的单元ï¼该技术简化为三个坐标方向上的TSL型Navier-Stokes方程近似[8]。文献[72]提出了一种保留完全黏性项的非单纯形单元推广。但此时便无法再使用高效的基于边的数据结构ï¼因为扩展的模板涉及与点 $$i$$ 没有边相连的节点。en

The disadvantage of this approach is however that the full viscous terms are retained only on triangular or tetrahedral grids. For elements like prisms or hexahedra, this technique simplifies to a TSL-like approximation of the Navier-Stokes equations in all three coordinate directions [8]. An extension to non-simplex elements which conserves the full viscous terms was presented in [72]. But the efficient edge-based data structure can then no longer be used since the extended stencil involves nodes not connected by an edge to point $$i$$ .

Average of Nodal Values 节点值平均

该方案面向单元中心型控制体与纯四面体网格。它采用文献[61]所引入模板的一个修改版本[4]来计算单元面上的梯度。该方法基于定义单元面的三个节点上数值的平均ï¼并结合单元形心处的量。单元面上的一阶导数由下列线性方程组的解得到[4]en

This scheme is intended for the cell-centred type of control volume and purely tetrahedral grids. It employs a modified version [4] of the stencil introduced in [61] to evaluate the gradients at the cell faces. The approach is based on an average of the values at the three nodes which define the cell face, combined with the quantities at the cell centroids. The first derivatives at a cell face result from the solution of the linear system of equations [4]

\begin{aligned} &\begin{bmatrix} x_J - x_I & y_J - y_I & z_J - z_I\\ \frac{1}{2}(x_2 + x_3) - x_1 & \frac{1}{2}(y_2 + y_3) - y_1 & \frac{1}{2}(z_2 + z_3) - z_1\\ \frac{1}{2}(x_1 + x_3) - x_2 & \frac{1}{2}(y_1 + y_3) - y_2 & \frac{1}{2}(z_1 + z_3) - z_2 \end{bmatrix} \begin{bmatrix} \partial_x U\\ \partial_y U\\ \partial_z U \end{bmatrix}\\ &= \begin{bmatrix} U_J - U_I\\ \frac{1}{2}(U_2 + U_3) - U_1\\ \frac{1}{2}(U_1 + U_3) - U_2 \end{bmatrix}, \end{aligned} \tag{5.73}

参照图5.19ï¼下标 $$I$$ 、 $$J$$ 表示单元形心ï¼下标 $$1$$ 、 $$2$$ 、 $$3$$ 分别表示节点 $$P_1$$ 、 $$P_2$$ 与 $$P_3$$ 。节点处的流动变量可由逆距离加权(式(5.45))或伪拉普拉斯加权[2](类似于式(5.24))确定。en

Referring to Fig. 5.19, the subscripts $$I$$ , $$J$$ denote the cell-centroids, and the subscripts $$1$$ , $$2$$ , $$3$$ stand for the nodes $$P_1$$ , $$P_2$$ and $$P_3$$ , respectively. Flow variables at the nodes can be determined either from the inverse-distance weighting (Eq. (5.45)) or by the pseudo-Laplacian weighting [2] (similar to Eq. (5.24)).

图5.18:节点i处的黏性项:体积为Ωe的四面体及计算与边ij相关系数所涉及的三角面

图5.18:节点 $$i$$ 处的黏性项:体积为 $\Omega_e$ 的四面体ï¼以及计算与边 $$ij$$ 相关系数所涉及的三角面[71]。图例: $$i$$ 、 $$j$$ ——节点; $\Omega_e$ ——共享边 $$ij$$ 的四面体的体积; $\vec{S}_e^{i}$ 、 $\vec{S}_e^{j}$ ——外侧面向量。

图5.19:单元中心格式:四面体网格上梯度计算的模板

图5.19:单元中心格式:四面体网格上梯度计算的模板[4]。图例: $$I$$ 、 $$J$$ ——单元形心; $$P_1$$ 、 $$P_2$$ 、 $$P_3$$ ——定义单元面的节点;叉号表示为计算黏性通量而求梯度的位置(面 $$P_1P_2P_3$$ 的中点);虚线为两形心连线。

5.4.2 Average of Gradients 梯度平均[cfd-5-4-2]

既然我们已经在每个控制体内部算出了梯度(例如用分段线性重构ï¼式(5.41)或(5.42))ï¼自然会想用简单平均来计算面中点处的梯度[73]en

Since we already computed the gradients inside each control volume (e.g., using the piecewise linear reconstruction, Eq. (5.41) or (5.42)), it would be tempting to evaluate the gradient at the face-midpoint by the simple average [73]

\overline{\nabla U}_{IJ} = \frac{1}{2}\left[\nabla U_I + \nabla U_J\right]. \tag{5.74}

这种方法特别有吸引力ï¼因为它只需要基本的基于面或边的数据结构ï¼不需要额外存储。然而ï¼正如文献[71]等指出的ï¼它导致一个权重分布不利的宽模板[49]。此外,[49]中证明ï¼该模板使解在四边形或六面体网格上可能发生失联(decoupling)。en

This approach is particularly attractive, because it requires only the basic face- or edge-based data structure and no additional storage. However, as it was pointed out, e.g., in Ref. [71], it leads to a wide stencil with an unfavourable distribution of the weights [49]. Furthermore, it was demonstrated in [49] that the stencil allows for the decoupling of the solution on quadrilateral or hexahedral grids.

利用沿单元形心连线方向的方向导数(对单元中心格式)ï¼可以改善该方法的性质ï¼尤其可以防止失联ï¼即en

The properties of the method can be improved, and particularly the decoupling can be prevented, by using the directional derivative along the connection between the cell-centroids (in the case of the cell-centred scheme), i.e.,

\left(\frac{\partial U}{\partial\ell}\right)_{IJ} \approx \frac{U_J - U_I}{\ell_{IJ}}\,, \tag{5.75}

其中 $\ell_{IJ}$ 表示两个单元形心 $$I$$ 与 $$J$$ 之间的距离(图5.19中的虚线)。对中点对偶格式也有类似表达式ï¼其中 $\vec{r}_{ij}$ 按式(5.43)定义。定义沿 $$I$$ 与 $$J$$ 连线的单位向量 $\vec{t}_{IJ}$ 为en

where $\ell_{IJ}$ represents the distance between the both cell-centroids $$I$$ and $$J$$ (dashed line in Fig. 5.19). A similar expression holds also for the median-dual scheme with $\vec{r}_{ij}$ according to Eq. (5.43). With the definition of the unit vector $\vec{t}_{IJ}$ along the line connecting $$I$$ and $$J$$ ,

\vec{t}_{IJ} = \frac{\vec{r}_{IJ}}{\ell_{IJ}}\,, \tag{5.76}

修正的平均可写为[74]、[75]en

the modified average may be written as [74], [75]

\nabla U_{IJ} = \overline{\nabla U}_{IJ} - \left[\overline{\nabla U}_{IJ}\cdot\vec{t}_{IJ} - \left(\frac{\partial U}{\partial\ell}\right)_{IJ}\right]\vec{t}_{IJ} \tag{5.77}

其中 $\overline{\nabla U}_{IJ}$ 由式(5.74)给出。这一修正在四面体以及棱柱或六面体网格上都导致强耦合的模板[49]。修正后的方法仍与基于面/边的数据结构相容ï¼且不需要额外存储。因此ï¼只要控制体内部的梯度反正要用于对流通量的计算ï¼它就比基于单元的方法更有吸引力。en

where $\overline{\nabla U}_{IJ}$ is given by Eq. (5.74). The modification leads to strongly coupled stencils on tetrahedral as well as on prismatic or hexahedral grids [49]. The modified approach is also still compatible with the face-/edge-based data structure and requires no additional storage. It is therefore more attractive than the element-based methodology, provided the gradients inside control volumes are utilised for the convective fluxes anyway.