OpenFOAM 速度输运方程中的应力张量项

无摘要

速度输运方程

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

对应的公式:

t(ρˉui~)+ρˉui~uj~xj=σijxj+xj(ρˉuiuj~)\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=τijpδij=2μ(Sij13Skkδ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=ρˉuiuj~\Gamma_{ij}=-\bar\rho \widetilde {u'_iu_j'}
根据 Boussinesq 假设:Γij13Γ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)T23tr(u)T))(ρνeffu)=(ρνeff(u+(u)T)23ρνefftr(u)T)=(ρνeff(u+(u)T)23ρνeff(u)I)=(2ρνeffSij23ρνeffSkkδij)=(2ρ(ν+νt)Sij23ρ(ν+ν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=(uxuyuzvxvyvzwxwywz)\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=(uxvxwxuyvywyuzvzwz)(\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)=ukxk=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(uixj+ujxi)=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=(ux12(uy+vx)12(uz+wx)S21vy12(vz+wy)S31S32wz)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=ukxk=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=Sij13tr(Sij)=(S1113(S11+S22+S33)S12S13S21S2213(S11+S22+S33)S23S31S32S3313(S11+S22+S33))=Sij13uδ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=ρˉuiuj~\overline \mathbf{\sigma}_t=-\bar\rho \widetilde {u'_iu_j'}

Boussinesq 假设:

σt13tr(σ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(ρˉuiuj~)=(ρˉuu~)+(ρˉvv~)+(ρˉww~)=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ˉ)I23ρˉ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(τt23ρˉIk)=2ρˉIktr(σt)=tr(2μtDˉ23μttr(D)I23ρˉ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=σijxj+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=pxixj(23ρˉδijk)+xj[2(μ+μt)(Sijukxkδ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=σijxjxj(Γ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 假设:

Γij13Γ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=pxi+xj(23ρˉδijk)+xj[2(μ+μt)(Sijukxkδ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 并不是真正的压力,而是p23ρ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 成长之路
微信打赏给博主更多动力吧~