OpenFOAM 成长之路 首页关于Timeline
OpenFOAM 速度输运方程中的应力张量项
| code analysis|字数总计: 2.9k|阅读时长: 16 分钟

无摘要

速度输运方程

tmp<fvVectorMatrix> tUEqn
(
    fvm::ddt(rho, U) + fvm::div(phi, U)
  + turbulence->divDevRhoReff(U)
 ==
  - fvc::grad(p)
);

对应的公式：

∂∂t(ρˉui~)+∂ρˉui~uj~∂xj=∂σ‾ij∂xj+∂∂xj(−ρˉui′uj′~)\frac{ \partial }{ \partial t }( \bar\rho\widetilde{u_i} ) + \frac{ \partial\bar\rho\widetilde{u_i}\widetilde{u_j} }{ \partial x_j } = \frac{ \partial \overline \sigma _{ij} }{ \partial x_j } + \frac{ \partial }{ \partial x_j}( -\bar\rho\widetilde{ u'_iu'_j } )

其中公式中的粘性应力
σ‾ij=τ‾ij−p‾δij=2μ(S‾ij−13S‾kkδij)−p‾δij\overline \sigma _{ij} = \overline\tau_{ij} - \overline p\delta_{ij} = 2\mu( \overline S_{ij} - \frac{1}{3}\overline S_{kk}\delta_{ij} ) - \overline p\delta_{ij}。

雷诺应力
Γij=−ρˉui′uj′~\Gamma_{ij}=-\bar\rho \widetilde {u'_iu_j'}，
根据 Boussinesq 假设：Γij−13Γkkδij=2μt(Sij~−13Skk~δij)\Gamma_{ij}-\frac{1}{3}\Gamma_{kk}\delta_{ij}=2\mu_t(\widetilde{S_{ij}}-\frac{1}{3}\widetilde{S_{kk}}\delta_{ij})，这里简化成了不可压版本：Γij=2μt(Sij~−13Skk~δij)\Gamma_{ij}=2\mu_t(\widetilde{S_{ij}}-\frac{1}{3}\widetilde{S_{kk}}\delta_{ij})。

这个 turbulence->divDevRhoReff(U) 就是公式中的粘性应力和雷诺应力移到方程左边了。
下面详细解释了这个过程：

template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
    volVectorField& U
) const
{
    return
    (
      - fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
    );
}

代码中表示的公式是

−∇⋅(ρ∗νeff((∇⊗u)T−23tr(∇⊗u)T))−∇⋅(ρ∗νeff∇⊗u)=−∇⋅(ρ∗νeff(∇⊗u+(∇⊗u)T)−23ρ∗νefftr(∇⊗u)T)=−∇⋅(ρ∗νeff(∇⊗u+(∇⊗u)T)−23ρ∗νeff(∇⋅u)I)=−∇⋅(2ρ∗νeffSij−23ρ∗νeffSkkδij)=−∇⋅(2ρ∗(ν+νt)Sij−23ρ∗(ν+νt)Skkδij)=−∇⋅(τij+Γij)\begin{array}{ll} &-\nabla\cdot\left(\rho*\nu_{eff}\left((\nabla\otimes\mathbf{u})^T-\dfrac{2}{3}tr(\nabla\otimes\mathbf{u})^T\right)\right)\\ &-\nabla\cdot(\rho*\nu_{eff}\nabla\otimes\mathbf{u})\\ &=-\nabla\cdot \left( \rho*\nu_{eff}(\nabla\otimes\mathbf{u}+(\nabla\otimes\mathbf{u})^T)-\dfrac{2}{3}\rho*\nu_{eff} tr(\nabla\otimes\mathbf{u})^T \right)\\ &=-\nabla\cdot \left( \rho*\nu_{eff}(\nabla\otimes\mathbf{u}+(\nabla\otimes\mathbf{u})^T)-\dfrac{2}{3}\rho*\nu_{eff}(\nabla \cdot \mathbf{u})\mathbf{I} \right)\\ &=-\nabla\cdot \left( 2\rho*\nu_{eff}S_{ij}-\dfrac{2}{3}\rho*\nu_{eff}S_{kk}\delta_{ij} \right)\\ &=-\nabla\cdot \left( 2\rho*(\nu+\nu_t)S_{ij}-\dfrac{2}{3}\rho*(\nu+\nu_t)S_{kk}\delta_{ij} \right)\\ &=-\nabla\cdot(\tau_{ij}+\Gamma_{ij}) \end{array}

∇⊗u=(∂u∂x∂u∂y∂u∂z∂v∂x∂v∂y∂v∂z∂w∂x∂w∂y∂w∂z)\nabla\otimes\mathbf{u}= \left( \begin{matrix} \dfrac{\partial u}{\partial x}&\dfrac{\partial u}{\partial y}&\dfrac{\partial u}{\partial z}\\ \dfrac{\partial v}{\partial x}&\dfrac{\partial v}{\partial y}&\dfrac{\partial v}{\partial z}\\ \dfrac{\partial w}{\partial x}&\dfrac{\partial w}{\partial y}&\dfrac{\partial w}{\partial z} \end{matrix} \right)

(∇⊗u)T=(∂u∂x∂v∂x∂w∂x∂u∂y∂v∂y∂w∂y∂u∂z∂v∂z∂w∂z)(\nabla\otimes\mathbf{u})^T= \left( \begin{matrix} \dfrac{\partial u}{\partial x}&\dfrac{\partial v}{\partial x}&\dfrac{\partial w}{\partial x}\\ \dfrac{\partial u}{\partial y}&\dfrac{\partial v}{\partial y}&\dfrac{\partial w}{\partial y}\\ \dfrac{\partial u}{\partial z}&\dfrac{\partial v}{\partial z}&\dfrac{\partial w}{\partial z} \end{matrix} \right)

tr(∇⊗u)=tr((∇⊗u)T)=∂uk∂xk=∇⋅u=tr(Sij)tr(\nabla\otimes\mathbf{u})=tr\left((\nabla\otimes\mathbf{u})^T\right)=\dfrac{\partial u_k}{\partial x_k}=\nabla \cdot \mathbf{u}=tr(S_{ij})

strain rateSijS_{ij} :

Sij=12(∂ui∂xj+∂uj∂xi)=12(∇⊗u+(∇⊗u)T)S_{ij}=\frac{1}{2}( \frac{ \partial u_i }{ \partial x_j } + \frac {\partial u_j }{ \partial x_i } )=\dfrac{1}{2}\left(\nabla\otimes\mathbf{u}+(\nabla\otimes\mathbf{u})^T\right)

Sij=(∂u∂x12(∂u∂y+∂v∂x)12(∂u∂z+∂w∂x)S21∂v∂y12(∂v∂z+∂w∂y)S31S32∂w∂z)S_{ij}= \left( \begin{matrix} \dfrac{\partial u}{\partial x}&\dfrac{1}{2}(\dfrac{\partial u}{\partial y}+\dfrac{\partial v}{\partial x})&\dfrac{1}{2}(\dfrac{\partial u}{\partial z}+\dfrac{\partial w}{\partial x})\\ S_{21}&\dfrac{\partial v}{\partial y}&\dfrac{1}{2}(\dfrac{\partial v}{\partial z}+\dfrac{\partial w}{\partial y})\\ S_{31}&S_{32}&\dfrac{\partial w}{\partial z} \end{matrix} \right)

SijS_{ij} 是对称矩阵，S12=S21, S13=S31, S23=S32S_{12}=S_{21},\ S_{13}=S_{31},\ S_{23} = S_{32}。

tr(Sij)=Skk=∂uk∂xk=∇⋅utr(S_{ij})=S_{kk}=\dfrac{\partial u_k}{\partial x_k}=\nabla\cdot \mathbf{u}

SijS_{ij} 在代码中是：symm(fvc::grad(U))，其中 symm(tensor) = 1/2*(tensor + T(tensor))。

deformation rate tensorDijD_{ij}:

Dij=Sij−13tr(Sij)=(S11−13(S11+S22+S33)S12S13S21S22−13(S11+S22+S33)S23S31S32S33−13(S11+S22+S33))=Sij−13∇⋅uδij\begin{array}{ll} &D_{ij}= S_{ij}-\frac{1}{3}tr(S_{ij})\\ &= \left( \begin{matrix} S_{11}-\frac{1}{3}(S_{11}+S_{22}+S_{33}) & S_{12} & S_{13} \\ S_{21} & S_{22}-\frac{1}{3}(S_{11}+S_{22}+S_{33}) & S_{23} \\ S_{31} & S_{32} & S_{33}-\frac{1}{3}(S_{11}+S_{22}+S_{33}) \end{matrix} \right) \\ &=S_{ij}-\frac{1}{3}\nabla\cdot \mathbf{u} \delta_{ij} \end{array}

DijD_{ij} 在代码中是：dev(symm(fvc::grad(U)))，也就是 dev(S_{ij})。

template <class Cmpt>
inline Tensor <Cmpt> dev(const Tensor <Cmpt >& t)
{
return t - SphericalTensor <Cmpt >::oneThirdsI*tr(t);
}

template <class Cmpt>
inline Tensor <Cmpt> dev2(const Tensor <Cmpt >& t)
{
return t - SphericalTensor <Cmpt >::twoThirdsI*tr(t);
}

参考文献：
Tobias Holzmann. Mathematics, Numerics, Derivations and OpenFOAM®, Holzmann CFD,
URL www.holzmann-cfd.de, DOI: 10.13140/RG.2.2.27193.36960.

以下公式适用于可压缩流动：

Cauchy stress tensorσ=σ=−pI+[σ−13tr(σ)I]\mathbf{\sigma}=\mathbf{\sigma}=-p\mathbf{I}+[\mathbf{\sigma}-\frac{1}{3}tr(\mathbf{\sigma})\mathbf{I}]

shear−rate tensorτ=σ−13tr(σ)I=σdev\mathbf{\tau}=\mathbf{\sigma}-\frac{1}{3}tr(\mathbf{\sigma})\mathbf{I}=\mathbf{\sigma}^{dev}

即：σ=−pI+τ\mathbf{\sigma}=-p\mathbf{I}+\mathbf{\tau}

Reynolds-Stress tensor
σ‾t=−ρˉui′uj′~\overline \mathbf{\sigma}_t=-\bar\rho \widetilde {u'_iu_j'}

Boussinesq 假设：

σ‾t−13tr(σ‾t)=2μtDˉ−23tr(Dˉ)I\overline \mathbf{\sigma}_t-\frac{1}{3}tr(\overline \mathbf{\sigma}_t)=2\mu_t \bar \mathbf{D}-\frac{2}{3}tr(\bar \mathbf{D})\mathbf{I}

这里的Dˉ\bar \mathbf{D} 是变形速率张量，deformation rate tensor，一般也称为SijS_{ij}。

tr(σ‾t)=tr(−ρˉui′uj′~)=(−ρˉu′u′~)+(−ρˉv′v′~)+(−ρˉw′w′~)=−2ρˉIktr(\overline \mathbf{\sigma}_t)=tr(-\bar\rho \widetilde {u'_iu_j'}) =(-\bar\rho \widetilde {u'u'})+(-\bar\rho \widetilde {v'v'})+(-\bar\rho \widetilde {w'w'})=-2\bar\rho \mathbf{I}k

故雷诺应力可以写为：

σ‾t=2μtDˉ−23tr(Dˉ)I+13tr(σ‾t)=2μtDˉ−23tr(Dˉ)I−23ρˉIk\overline \mathbf{\sigma}_t=2\mu_t \bar \mathbf{D}-\frac{2}{3}tr(\bar \mathbf{D})\mathbf{I}+\frac{1}{3}tr(\overline \mathbf{\sigma}_t) =2\mu_t \bar \mathbf{D}-\frac{2}{3}tr(\bar \mathbf{D})\mathbf{I}-\frac{2}{3}\bar\rho \mathbf{I}k

tr(σ‾t)=tr(τ‾t−23ρˉIk)=−2ρˉIktr(σ‾t)=tr(2μtDˉ−23μttr(D)I−23ρˉIk)=−2ρˉIktr(\overline \mathbf{\sigma}_t)=tr(\overline \mathbf{\tau}_t-\frac{2}{3}\bar\rho \mathbf{I}k)=-2\bar\rho \mathbf{I}k\\ tr(\overline \mathbf{\sigma}_t)=tr(2\mu_t\bar \mathbf{D}-\frac{2}{3}\mu_ttr(\mathbf{D})\mathbf{I}-\frac{2}{3}\bar\rho \mathbf{I}k)=-2\bar\rho \mathbf{I}k

对于 RANS

速度输运方程：

∂∂t(ρˉui~)+∂ρˉui~uj~∂xj=∂σ‾ij∂xj+∂∂xj(Γij)\frac{ \partial }{ \partial t }( \bar\rho\widetilde{u_i} ) + \frac{ \partial\bar\rho\widetilde{u_i}\widetilde{u_j} }{ \partial x_j } = \frac{ \partial \overline \sigma _{ij} }{ \partial x_j } + \frac{ \partial }{ \partial x_j}( \Gamma_{ij} )

这里的Γij\Gamma_{ij} 对应上文中的σt\sigma_t。

最后变形为：

∂∂t(ρˉui~)+∂ρˉui~uj~∂xj=−∂p∂xi−∂∂xj(23ρˉδijk)+∂∂xj[2(μ+μt)(Sij−∂uk∂xkδij)]\frac{ \partial }{ \partial t }( \bar\rho\widetilde{u_i} ) + \frac{ \partial\bar\rho\widetilde{u_i}\widetilde{u_j} }{ \partial x_j } = -\frac{\partial p}{\partial x_i} -\frac{\partial}{\partial x_j}\left(\frac{2}{3}\bar\rho \delta_{ij}k\right) +\frac{\partial}{\partial x_j}\left[2(\mu+\mu_t)(S_{ij}-\frac{\partial u_k}{\partial x_k}\delta_{ij})\right]

OF 速度输运方程代码里边并没有包含kk 这一项，因此代码里的pp 并不是真正的压力，而是p+23ρkp+\frac{2}{3}\rho k。Page 76 in the book of Tobias.

对于 LES

速度输运方程：

∂∂t(ρˉui~)+∂ρˉui~uj~∂xj=∂σ‾ij∂xj−∂∂xj(Γij)\frac{ \partial }{ \partial t }( \bar\rho\widetilde{u_i} ) + \frac{ \partial\bar\rho\widetilde{u_i}\widetilde{u_j} }{ \partial x_j } = \frac{ \partial \overline \sigma _{ij} }{ \partial x_j } - \frac{ \partial }{ \partial x_j}( \Gamma_{ij} )

Boussinesq 假设：

Γij−13Γkkδij=−2μt(Sij~−13Skk~δij)\Gamma_{ij}-\frac{1}{3}\Gamma_{kk}\delta_{ij}=-2\mu_t(\widetilde{S_{ij}}-\frac{1}{3}\widetilde{S_{kk}}\delta_{ij})

LES 中有tr(Γij)=Γkk=2ρˉktr(\Gamma_{ij})=\Gamma_{kk}=2\bar\rho k

因此：

Γij=−2μt(Sij~−13Skk~δij)−23ρˉkδij\Gamma_{ij}=-2\mu_t(\widetilde{S_{ij}}-\frac{1}{3}\widetilde{S_{kk}}\delta_{ij})-\frac{2}{3}\bar\rho k \delta_{ij}

速度输运方程可以变形为：

∂∂t(ρˉui~)+∂ρˉui~uj~∂xj=−∂p∂xi+∂∂xj(23ρˉδijk)+∂∂xj[2(μ+μt)(Sij−∂uk∂xkδij)]\frac{ \partial }{ \partial t }( \bar\rho\widetilde{u_i} ) + \frac{ \partial\bar\rho\widetilde{u_i}\widetilde{u_j} }{ \partial x_j } = -\frac{\partial p}{\partial x_i} +\frac{\partial}{\partial x_j}\left(\frac{2}{3}\bar\rho \delta_{ij}k\right) +\frac{\partial}{\partial x_j}\left[2(\mu+\mu_t)(S_{ij}-\frac{\partial u_k}{\partial x_k}\delta_{ij})\right]

那么按照 Tobias 的说法，OF 速度输运方程代码里边并没有包含kk 这一项，因此代码里的pp 并不是真正的压力，而是p−23ρkp-\frac{2}{3}\rho k？

按照上边的观点，OF代码中的 p 其实并不是真实的压力，那么它怎么能用于气体状态方程呢？这个问题是论坛里网友提出来的：cfd-online。

其实最简单的就是把缺失的那个 k 补上，为什么大家不这么做呢？看看网上的讨论。
https://www.cfd-online.com/Forums/openfoam-solving/58214-calculating-divdevreff.html

一些矢量和张量的操作

LES 中经常出现的∣S~∣=2Sij~Sij~|\tilde S|=\sqrt{2\widetilde{S_{ij}}\widetilde{S_{ij}}}，用代码表示是：

volSymmTensorField Sij = symm(fvc::grad(U));
volScalarField magSij = sqrt(2*(Sij && Sij));

这里的 && 的操作是缩并，也叫双内积，两个二阶张量的缩并是一个标量。

//- Double-dot-product between a symmetric tensor and a symmetric tensor
template<class Cmpt>
inline Cmpt
operator&&(const SymmTensor<Cmpt>& st1, const SymmTensor<Cmpt>& st2)
{
    return
    (
        st1.xx()*st2.xx() + 2*st1.xy()*st2.xy() + 2*st1.xz()*st2.xz()
                          +   st1.yy()*st2.yy() + 2*st1.yz()*st2.yz()
                                                +   st1.zz()*st2.zz()
    );
}

& 是内积，也叫数量积，也叫点乘，两个矢量的内积是一个标量。
^ 是叉乘，两个矢量的叉乘，仍旧是一个矢量。在三维几何中，向量a和向量b的叉乘结果是一个向量，更为熟知的叫法是法向量，该向量垂直于a和b向量构成的平面。

template<class Cmpt>
inline Vector<Cmpt> operator^(const Vector<Cmpt>& v1, const Vector<Cmpt>& v2)
{
    return Vector<Cmpt>
    (
        (v1.y()*v2.z() - v1.z()*v2.y()),
        (v1.z()*v2.x() - v1.x()*v2.z()),
        (v1.x()*v2.y() - v1.y()*v2.x())
    );
}

外积 a*b，也叫量积（tensor-product），两个矢量可以形成一个二阶张量。

template<class Cmpt>
inline typename outerProduct<Vector<Cmpt>, Vector<Cmpt>>::type
operator*(const Vector<Cmpt>& v1, const Vector<Cmpt>& v2)
{
    return Tensor<Cmpt>
    (
        v1.x()*v2.x(), v1.x()*v2.y(), v1.x()*v2.z(),
        v1.y()*v2.x(), v1.y()*v2.y(), v1.y()*v2.z(),
        v1.z()*v2.x(), v1.z()*v2.y(), v1.z()*v2.z()
    );
}

sqr(U)=uiujsqr(U)=u_iu_j

这个全局函数没看懂什么意思：

template<class Cmpt>
inline Vector<Cmpt> perpendicular(const Vector<Cmpt>& v)
{
    Vector<Cmpt> u(Zero);
    u[findMin(cmptMag(v))] = 1;
    return u ^ v;
}
文章作者: Yan Zhang
文章链接: https://openfoam.top/shearRateTensor/
版权声明: 本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 OpenFOAM 成长之路！
TurbulenceModels
您的肯定会给我更大动力~