在前面介绍的矢量运算和张量表达法的基础上，本文将主要介绍CFD中与张量相关的运算。从本文开始，如无特别说明，提到张量时均指的是二阶张量。
矢量的并矢积
两个矢量的并矢积可以写作 
 或 
 ，它是将两个矢量的各个方向进行并置，从而得到一个二阶张量：
也可以将并矢积看作是一个列向量与一个行向量相乘所得的矩阵：
根据这一运算规则，我们可以得到矢量的梯度，其结果是一个张量：
张量的运算
张量的运算与矩阵运算具有高度的一致性：
张量相加
2. 标量乘以张量
3. 张量点乘矢量
根据这一运算规则，我们可以得到张量的散度，其结果是一个矢量：
4. 两个张量（或并矢）的点乘
按照矩阵乘法进行运算，其结果仍为张量：
上式中 
 均为矢量，
 和 
 则为两个矢量的并矢积，即张量；等式右边 
 是一个标量， 
 为张量。
5. 张量的双点积
两个张量的双点积的运算是基于以下两个矢量并矢积的双点乘运算规则：
上式中 
 均为矢量， 
 则为两个矢量的并矢积，即张量；而 
 和 
 均是两个矢量的点乘，其结果是标量。由此可见，
张量的双点积的结果是一个标量。
将这一规则应用到笛卡尔坐标系的单位矢量上，可以得到：
上面仅列举了部分结果，其他单位矢量的双点积结果也容易得到。
由此可以得到，黏性力张量和速度梯度张量的双点积结果为：
以上用四篇小文章给出了CFD中常见的张量运算法则，在公式推导中还会涉及更多张量运算公式，但基本都可以通过以上这些基本法则来进行推导和证明，在后面文章用到时会进行简单说明。下一篇文章将给出推导N-S方程时会用到的几个积分定理。
参考资料
[1] Moukalled, F., Mangani, L., and Darwish, M. The Finite Volume Method in Computational Fluid Dynamics : An Advanced Introduction with OpenFOAM and Matlab. 2016.

张量的几个物理场景
应力状态
材料内部受力状况就是应力状态，应力在不同方向的大小不同，并且与感应面有关。如图所示材料某点的应力状态，可用矩阵表示

$\sigma=\begin{pmatrix}
\sigma_x&\tau_{xy}&\tau_{xz}\\ \tau_{yx}&\sigma_y&\tau_{yz}\\ \tau_{zx}&\tau_{zy}&\sigma_z
\end{pmatrix}$

这是个对称矩阵，即$\tau_{xy}=\tau_{yx},\tau_{xz}=\tau_{zx},\tau_{yz}=\tau_{zy}$

现在已知某任意平面，那么怎么求该平面上的应力状态呢?如下图所示，以平面示意

红色的箭头是该平面的总的应力，该总应力可以分解成垂直于平面的正应力和平行于平面的切应力。若已知平面的法矢为$n=(n_x,n_y,n_z)$，其总应力计算如下

$\sigma'=n\cdot \sigma=(n_x,n_y,n_z) \begin{pmatrix}
\sigma_x&\tau_{xy}&\tau_{xz}\\ \tau_{yx}&\sigma_y&\tau_{yz}\\ \tau_{zx}&\tau_{zy}&\sigma_z
\end{pmatrix}$

那么正应力为

$\sigma'_n=\sigma'/n\cdot n=\sigma'\cdot n\cdot n\\=(n_x,n_y,n_z) \begin{pmatrix}
\sigma_x&\tau_{xy}&\tau_{xz}\\ \tau_{yx}&\sigma_y&\tau_{yz}\\ \tau_{zx}&\tau_{zy}&\sigma_z
\end{pmatrix}
\begin{pmatrix}
n_x\\n_y\\ n_z
\end{pmatrix}(n_x,n_y,n_z)$

切应力的为

$\sigma'_\tau=\sigma'\div n=\sigma'-\sigma'_n$
转动惯量
一个物体的转动，和平移一样，也存在惯性。转动惯量就是衡量这个惯性的大小。对于规则的物体，其转动比较容易求得，如下图所示

该立方块绕$x$轴的转动惯量为$J_x$，绕$y$轴的转动惯量为$J_y$，绕$z$轴的转动惯量为$J_z$。那么如何求得其绕任意轴$n$的转动惯量$J'$，如图红色所示。转动惯量也是一种张量，其数学形式如下

$J=\begin{pmatrix}
J_x&0&0\\0&J_y&0\\ 0& 0&J_z
\end{pmatrix}$

那么绕$n$的总转动惯量计算如下

$J'=n\cdot J=(n_x,n_y,n_z)\begin{pmatrix}
J_x&0&0\\0&J_y&0\\ 0& 0&J_z
\end{pmatrix}$

总的转动惯量到转轴的投影就是绕这个轴的转动惯量

$J'_n=J'/n\cdot n$

总的转动惯量垂直于这个轴的投影为

$J'_\tau=J'\div n$

有个疑问，为什么会有这个垂直的转动惯量。这是因为，在物体转动的时候，两端会甩开远离轴，或者聚拢越来越接近轴，而这个转动的转动惯量就是总转动惯量垂直于这个轴的分量。
各向异性的感应电场
我们知道，一个电荷在真空中形成的电场是向各个方向发散，并且每个方向等强度，即等势面是球形。但如果在物质内部呢?这个电荷在物质内部形成的是感应电场(也常称为电位移矢量)，并且由于很多物质都是各向异性的，也就是说电荷在不同方向的感应效果不一样，所以等势面试是椭球，如下图所示



其实所有三维张量都可以用椭球表示，二维张量可以用椭圆表示。

电荷在真空发出的电场为

$E=\frac{Q}{4\pi r^2}\begin{pmatrix}
\epsilon&0&0\\0&\epsilon&0\\0&0&\epsilon
\end{pmatrix}$

电荷在各向异性物质内感应的电场为

$D=\frac{Q}{4\pi r^2}\begin{pmatrix}
\epsilon_{xx}&\epsilon_{xy}&\epsilon_{xz}\\\epsilon_{yx}&\epsilon_{yy}&\epsilon_{yz}\\\epsilon_{zx}&\epsilon_{zy}&\epsilon_{zz}
\end{pmatrix}$

该式中的矩阵为对称矩阵，即$\epsilon_{xy}=\epsilon_{yx},\epsilon_{xz}=\epsilon_{zx},\epsilon_{yz}=\epsilon_{zy}$。矩阵中$\epsilon_{ij}=\epsilon^r_{ij}\epsilon$，其中$\epsilon^r_{ij}$为方向$ij$的相对介电系数，所以任意方向的感应电场强度为

$D'=n\cdot D=(n_x,n_y,n_z)\begin{pmatrix}
\epsilon_{xx}&\epsilon_{xy}&\epsilon_{xz}\\\epsilon_{yx}&\epsilon_{yy}&\epsilon_{yz}\\\epsilon_{zx}&\epsilon_{zy}&\epsilon_{zz}
\end{pmatrix}$

其大小为

$|D'|=D'/n=D'\cdot n=(n_x,n_y,n_z)\begin{pmatrix}
\epsilon_{xx}&\epsilon_{xy}&\epsilon_{xz}\\\epsilon_{yx}&\epsilon_{yy}&\epsilon_{yz}\\\epsilon_{zx}&\epsilon_{zy}&\epsilon_{zz}
\end{pmatrix}
\begin{pmatrix}
n_x\\n_y\\n_z
\end{pmatrix}$

这就是一个二次型，其图像就是椭球，所以说张量可以用椭球表示。
张量的散度
我们知道张量可以表示感应电场，其实他可以表示任意种不同方向不同大小和方向的场，也就是张量场。既然有场，那就有散度的概念。已知下面的张量

$M=\begin{pmatrix}
a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}
\end{pmatrix}$

那么该场为

$f=(x,y,z)\begin{pmatrix}
a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}
\end{pmatrix}$

散度有两种计算方式，一种是极限的计算方式

$div f=\lim_{v \to 0}\oint_{s}fds/v$

还有一种方式为$div f=\nabla \cdot f=\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}+\frac{\partial f}{\partial z}$

采用第二种比较容易计算，那么张量场的散度为

$div f=\nabla \cdot f=a_{11}+a_{22}+a_{33}$

也就是说张量的散度就是张量的迹，可知一个方阵无论经过什么样的相似变换，它的迹不变；也就是说场无论经过什么样的坐标变换，它的散度不变。

网站、

https://wenku.baidu.com/view/ab76fe504b73f242336c5f81.html
（介绍二阶张量的梯度和散度）

在OpenFOAM中，我们用函数divDevReff 与 divDevRhoReff来计算剪切应力。

文章目录
1.函数名称含义
2.不可压缩流
3.可压缩流

1.函数名称含义
下面，解释这两个函数名称的由来：

$\nabla \cdot \tau=\nabla \cdot \sigma^{dev}$
关于这个公示的推导可查看Cauchy应力张量、剪切力张量与压力的关系。

div:散度（divergence）；
dev:每个二阶张量都可分解为 dev 部分（deviatoric）和 hyd 部分（hydrostatic）；
R:雷诺时均(Reynolds-Average Simulation)；
eff:有效的(effective)流动方式，包括 层流（laminar）或湍流(turbulent)；
rho:取决于流体密度是否发生变化。

由此，不难判断divDevReff 用于不可压缩流体，divDevRhoReff用于可压缩流体 。
2.不可压缩流
2.1代码解读
对于turbulence->divDevReff(U)，即操作对象turbulence通过指针调用名为divDevReff(U)的函数。该函数声明如下：

(文件位置：)
tmp<fvVectorMatrix>laminar::divDevReff(volVectorField &U) const 
{
	Return
	(
		-fvm::laplacian(nuEff(), U)
		-fvc::div(nuEff()*dev(T(fvc::grad(U))))
	)
}

注意：计算过程中，采用的是运动粘度。方程一部分用了隐式(fvm)，而另一部分用了显式(fvc)。此外，还调用了一个名为dev()的函数:
//- Return the deviatoric part of a tensor
template<class Cmpt>
inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
{
    return t - SphericalTensor<Cmpt>::oneThirdI*tr(t);
}

（文件位置：$FOAM_SRC/turbulenceModel/incompressible/RAS/laminar/laminar.C）
2.2函数方程
关于函数dev()，可以简单地理解为一个矩阵的dev(deviatoric part)部分：

$\mathbf A^{dev}=\mathbf A-\mathbf A^{hyd}=\mathbf A-\frac13 tr(\mathbf A)\mathbf I$
函数dev()的参数为速度场的转置梯度矩阵:

$\mathbf A	=[\nabla \mathbf U]^T$
我们对代码进行重写：
-fvm::laplacian(nuEff(), U)=$-\nabla\cdot\left(\nu^{eff}(\nabla \mathbf U)\right)$
-fvc::div(nuEff()*dev(T(fvc::grad(U))))=$-\nabla\cdot\left(\nu^{eff}*dev((\nabla \mathbf U)^T)\right)$
$dev((\nabla \mathbf U)^T)=(\nabla \mathbf U)^T-\frac 13 tr((\nabla \mathbf U)^T)\mathbf I$
将上述代码整合在一起，有：

$-\nabla\cdot \mathbf R^{eff}=-\nabla\cdot\left(\nu^{eff}(\nabla \mathbf U)\right)-\nabla\cdot
\left(\nu^{eff}
\left[
(\nabla \mathbf U)^T-\frac 13 tr((\nabla \mathbf U)^T)\mathbf I
\right]
\right)$

整理，可得：

$-\nabla\cdot \mathbf R^{eff}=-\nabla\cdot
\left(
\nu^{eff}
\left[
\nabla \mathbf U+(\nabla \mathbf U)^T-\frac 13 tr((\nabla \mathbf U)^T)\mathbf I
\right]
\right)$

根据CFD张量公式

$\mathrm{tr}\left(\nabla\mathbf{U}\right)\mathbf I=\mathrm{tr}\left((\nabla\mathbf{U})^{\mathrm{T}}\right)\mathbf I=\left(\nabla\cdot\mathbf{U}\right)\bf I$

对于不可压缩流连续性方程 $\nabla\cdot \mathbf U=0$，则：

$-\nabla\cdot \mathbf R^{eff}=-\nabla\cdot
\left(
\nu^{eff}
\left[
\nabla \mathbf U+(\nabla \mathbf U)^T
\right]
\right)=-\nabla\cdot
\left(
2\nu^{eff}
\underbrace{\left[\frac12
\left(
\nabla \mathbf U+(\nabla \mathbf U)^T
\right)
\right] }_{`形变率 \mathbf D`}
\right)=\nabla\cdot \tau$
对于不可压缩流，$\frac 13 tr((\nabla \mathbf U)^T)$ 应该恒等于零，是可以忽略的，但考虑到连续性不可能百分百为零，这里我们将其作为误差引入进去。
3.可压缩流
对于可压缩流，调用turbulence->divDevRhoReff(U)函数。该函数包含关键字Rho，这说明在计算可压缩流的剪切应力时，需要考虑由密度变化引起的膨胀粘度项。
3.1代码解读
函数divDevRhoReff(U)
tmp<fvVectorMatrix>laminar::divDevRhoReff(volVectorField &U) const
{
	 Return
	(
		 -fvm::laplacian(muEff(), U) 
		 -fvc::div(muEff()*dev2(T(fvc::grad(U)))) 
	)
}

注意，这里我们应用的是dev2()函数，这种新函数的代码如下：
//- Return the deviatoric part of a tensor
template<class Cmpt>
inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
{
    return t - SphericalTensor<Cmpt>::twoThirdI*tr(t);
}

3.2 函数方程
关于函数dev2()，可以理解为一个矩阵减去两次hyd部分
$\mathbf A^{dev2}=\mathbf A-2\mathbf A^{hyd}=\mathbf A-\frac23 tr(\mathbf A)\mathbf I$

经整理（方法同上），其方程为：

$-\nabla\cdot \mathbf R^{eff}=-\nabla\cdot
\left(
2\nu^{eff}
\underbrace{\left[\frac12
\left(
\nabla \mathbf U+(\nabla \mathbf U)^T
\right)
\right] }_{`形变率 \mathbf D`}-\frac23\mu^{eff}(\nabla \cdot \mathbf U)\mathbf I
\right)=\nabla\cdot \tau$