dynLES

动态 SGS 模型

输运方程

可压缩流动的N-S方程如下:

ρt+ρuixi=0ρuit+ρuiujxj=pxi+σijxj+ρfi\begin{array}{l} \frac{\partial \rho }{\partial t} + \frac{\partial \rho {u_i}}{\partial {x_i}} = 0\\ \frac{\partial \rho u_i} {\partial t} + \frac{\partial \rho u_iu_j} {\partial x_j} = - \frac{\partial p}{\partial x_i} + \frac{\partial \sigma _{ij}} {\partial x_j} + \rho f_i \end{array}

其中,σij=2μSij23μ(ukxk)δij ,Sij=12(uixj+ujxi)\sigma _{ij} = 2\mu S_{ij} - \frac{2}{3}\mu (\frac{\partial u_k} {\partial x_k}) \delta _{ij}\ ,S_{ij} = \frac{1}{2}(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})

对N-S方程进行雷诺平均,并引入Favre平均 :fi~=ρfiρˉ\widetilde {f_i} = \frac {\overline {\rho f_i} } {\bar \rho },可以得到:

ρˉt+ρˉui~xi=0t(ρˉui~)+ρˉui~uj~xj=pˉxi+σij~xj+ρˉfi~+xj(ρˉuiuj~)\begin{array}{l} \frac{\partial \bar \rho }{\partial t} + \frac{\partial \bar \rho \widetilde {u_i}}{\partial x_i} = 0\\ \frac{\partial }{\partial t}(\bar \rho \widetilde {u_i}) + \frac{\partial \bar \rho \widetilde {u_i}\widetilde {u_j}}{\partial x_j} = - \frac{\partial \bar p}{\partial x_i} + \frac{\partial \widetilde {\sigma _{ij}}}{\partial x_j} + \bar \rho \widetilde {f_i} + \frac{\partial }{\partial x_j}( - \bar \rho \widetilde {u'_iu'_j}) \end{array}

其中,

σij~=2μSij~23μSkk~δij,Sij~=12(ui~xj+uj~xi)\widetilde {\sigma _{ij}} = 2\mu \widetilde {S_{ij}} - \frac{2}{3}\mu \widetilde {S_{kk}}{\delta _{ij}}{\rm{ , }}\widetilde {S_{ij}} = \frac{1}{2}(\frac{\partial \widetilde {u_i}}{\partial {x_j}} + \frac{\partial \widetilde {u_j}}{\partial x_i})

uiuj~=ui~uj~(ui~uj~uiuj~)\widetilde {u_iu_j} = \widetilde {u_i}\widetilde {u_j} - (\widetilde {u_i}\widetilde {u_j} - \widetilde {u_iu_j})

LES的动量方程为:

t(ρˉui~)+ρˉui~uj~xj=pˉxi+σij~xj+ρˉfi~xjρˉ(uiuj~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 \bar p}{\partial x_i} + \frac{\partial \widetilde {\sigma _{ij}}}{\partial x_j} + \bar \rho \widetilde {f_i} - \frac{\partial }{\partial x_j}\bar \rho (\widetilde {u_iu_j} - \widetilde {u_i}\widetilde {u_j})

定义亚格子应力为τij=ρˉ(uiuj~ui~uj~)\tau _{ij} = \bar \rho (\widetilde {u_iu_j} - \widetilde {u_i}\widetilde {u_j})

这样,动量方程就变为:

t(ρˉui~)+ρˉui~uj~xj=pˉxi+σij~xj+ρˉfi~τijxj\frac{\partial }{\partial t}(\bar \rho \widetilde {u_i}) + \frac{\partial \bar \rho \widetilde {u_i}\widetilde {u_j}}{\partial x_j} = - \frac{\partial \bar p}{\partial x_i} + \frac{\partial \widetilde {\sigma _{ij}}}{\partial x_j} + \bar \rho \widetilde {f_i} - \frac{\partial \tau _{ij}}{\partial x_j}

Smagorinsky模型:

τijSmag13τkkSmagδij=2μSGSSij~\tau _{ij}^{Smag} - \frac{1}{3}\tau _{kk}^{Smag}{\delta _{ij}} = - 2{\mu _{SGS}}\widetilde {S_{ij}}

其中,μSGS=ρˉ(CsΔˉ)2S~{\mu _{SGS}} = \bar \rho {({C_s}\bar \Delta )^2}|\widetilde S|,S~=2Sij~Sij~|\widetilde S| = \sqrt {2\widetilde {S_{ij}}\widetilde {S_{ij}}}

Yoshizawa.1986对其进行了可压缩修正,我们称之为可压 Smagorinsky 模型:

τij13τkkδij=2ρˉ(CsΔˉ)2S~(Sij~13Skk~δij)=Cs2αijτkk=CI2ρˉΔˉ2S~2=CIα\begin{array}{l} \tau _{ij} - \frac{1}{3}\tau _{kk}\delta _{ij} = - 2\textcolor{red}{\bar \rho (C_s{\bar \Delta })^2 |\widetilde S|}(\widetilde {S_{ij}} - \frac{1}{3}\widetilde {S_{kk}}\delta _{ij}) = C_s^2\alpha _{ij}\\ \tau _{kk} = C_I 2\bar \rho {\bar \Delta }^2|\widetilde S{|^2} = {C_I}\alpha \end{array}

其中,Cs=0.16,CI=0.09,S~=2Sij~Sij~{C_s} = 0.16{\rm{ }},{\rm{ }}{C_I} = 0.09{\rm{ }},{\rm{ }}|\widetilde S| = \sqrt {2\widetilde {S_{ij}}\widetilde {S_{ij}}}
1991 年,Germano 提出了动态 Smagorinsky 模型,该模型通过局部自相似假设,将两个不同尺度的亚格子应力模型联系起来,得到在计算过程中动态决定模型系数的方法。动态大涡模型的过程中,引入了二次滤波,称为检验滤波(test filter),检验滤波的滤波宽度为Δ^\widehat \Delta,比一次滤波的滤波宽度要大,一般为Δ^=2Δˉ\widehat \Delta = 2\bar \Delta 。但是 Germano 提出的动态模型是基于不可压流体的,这里我们采用Moin.1991提出的可压动态 SGS 模型。对一次滤波之后的动量方程进行检验滤波,得到如下方程:

t(ρˉui~^)+ρˉui~uj~^xj=pˉ^xi+σij~^xj+ρˉfi~^xjρˉ(uiuj~ui~uj~)^\frac{\partial }{\partial t}(\widehat {\bar \rho \widetilde {u_i}}) + \frac{\partial \widehat {\bar \rho \widetilde {u_i}\widetilde {u_j}}}{\partial x_j} = - \frac{\partial \widehat {\bar p}}{\partial x_i} + \frac{\partial \widehat {\widetilde {\sigma _{ij}}}}{\partial x_j} + \widehat {\bar \rho \widetilde {f_i}} - \frac{\partial }{\partial x_j}\widehat {\bar \rho (\widetilde {u_iu_j} - \widetilde {u_i}\widetilde {u_j})}

根据ρˉui~=ρui\bar \rho \widetilde {u_i} = \overline {\rho {u_i}}ρˉuiuj~=ρuiuj\bar \rho \widetilde {u_iu_j} = \overline {\rho {u_i}{u_j}} ,我们得到:

t(ρui^)+xj(ρui^ρui^ρˉ^)=pˉ^xi+σij~^xj+ρˉfi~^xjρuiuj^+xjρˉui~uj~^+xj(ρuiρuj^^ρˉ^)xjρˉui~uj~^\frac{\partial } {\partial t} (\widehat{\overline{\rho {u_i}}}) +\frac{\partial } {\partial {x_j}} (\frac{\widehat{\overline{\rho u_i}}\widehat{\overline{\rho u_i}}}{\widehat{\bar\rho}}) = - \frac{\partial \widehat{\bar p}} {\partial x_i} + \frac{\partial \widehat{\widetilde{\sigma_{ij}}}} {\partial x_j} + \widehat{\bar\rho \widetilde{f_i}} - \frac{\partial } {\partial x_j} \widehat{\overline{\rho u_iu_j}} \bcancel{+\frac{\partial }{\partial {x_j}}\widehat{\bar\rho \widetilde{u_i}\widetilde{u_j}}} + \frac{\partial } {\partial x_j} (\frac{\widehat{\overline{\rho u_i}\widehat{\overline{\rho u_j}}}} {\widehat{\bar\rho }}) \bcancel{-\frac{\partial }{\partial x_j}\widehat{\bar{\rho }\widetilde{u_i}\widetilde{u_j}}}

t(ρui^)+xj(ρui^ρui^ρˉ^)=pˉ^xi+σij~^xj+ρˉfi~^Tijxj\frac{\partial }{\partial t} (\widehat{\overline{\rho {u_i}}})+\frac{\partial }{\partial {x_j}}(\frac{\widehat{\overline{\rho u_i}}\widehat{\overline{\rho {u_i}}}}{\widehat{\bar{\rho }}})=-\frac{\partial \widehat{\bar{p}}}{\partial x_i}+\frac{\partial \widehat{\widetilde{\sigma_{ij}}}}{\partial x_j}+\widehat{\bar{\rho }\widetilde{f_i}}-\frac{\partial T_{ij}}{\partial x_j}

Tij=ρuiuj^ρui^ρuj^ρˉ^T_{ij}= \widehat{\overline{\rho u_iu_j}} - \frac{\widehat{\overline{\rho u_i}}\widehat{\overline{\rho u_j}}} {\widehat{\bar\rho }}

TijT_{ij} 是两次滤波之后的亚格子应力。

Lij=Tijτij^=ρuiuj^ρui^ρuj^ρˉ^(ρˉ^uiuj~ρˉ^ui~uj~) =ρuiuj^ρui^ρuj^ρˉ^ρuiuj^+ρuiρuj^ρˉ^ =ρuiρuj^ρˉ^ρui^ρuj^ρˉ^ =ρˉui~uj~^ρˉui~^ρˉuj~^ρˉ^\begin{array}{l} L_{ij}=T_{ij}-\widehat {\tau _{ij}}=\widehat{\overline{\rho {u_i}{u_j}}}-\frac{\widehat{\overline{\rho {u_i}}} \widehat{\overline{\rho u_j}}}{\widehat{\bar{\rho }}} -(\widehat{\bar\rho }\widetilde{u_iu_j} -\widehat{\bar\rho } \tilde{u_i}\tilde{u_j})\\ \text{ }=\bcancel{\widehat{\overline{\rho {u_i}u_j}}}-\frac{\widehat{\overline{\rho {u_i}}}\widehat{\overline{\rho u_j}}}{\widehat{\bar{\rho }}}\bcancel{-\widehat{\overline{\rho {u_{i}}{u_j}}}}+\frac{\widehat{\overline{\rho u_i}\overline{\rho {u_j}}}}{\hat{\bar{\rho }}}\\\text{ }=\frac{\widehat{\overline{\rho {u_{i}}}\overline{\rho u_j}}}{\hat{\bar{\rho }}}-\frac{\widehat{\overline{\rho {u_{i}}}}\widehat{\overline{\rho {u_{j}}}}}{\widehat{\bar{\rho }}}\\\text{ }=\widehat{\bar{\rho }\widetilde{u_{i}}\widetilde{u_{j}}}-\frac{\widehat{\bar{\rho }\widetilde{u_i}}\widehat{\bar{\rho }\widetilde{u_j}}}{\widehat{\bar{\rho }}} \end{array}

这就是所谓的 Leonard 应力,上次称为 Germano 恒等式?

分别对τij^\widehat{\tau_{ij}}TijT_{ij} 采用可压 Smagorinsky 模型:

τij^13τkk^=Cs2αij^Tij13Tkkδij=Cs22ρˉ^Δ^2S¨(S¨ij13S¨kkδij)=Cs2βijwhere:f¨=ρf^ρˉ^αij=2ρˉΔˉ2S~(Sij~13Skk~δij)βij=2ρˉ^Δ^2S¨(S¨ij13S¨kkδij)\widehat {\tau_{ij}}-\frac{1}{3}\widehat{\tau_{kk}}=C_s^2\widehat{\alpha_{ij}} \\ T_{ij}-\frac{1}{3} {T_{kk}}\delta_{ij}=-C_s^22\widehat{\bar\rho} \widehat \Delta^2|\ddot S|(\ddot S_{ij}-\frac{1}{3}\ddot S_{kk}\delta_{ij})=C_s^2\beta_{ij}\\ where: \ddot f=\frac{\widehat{\overline{\rho f}}}{\widehat{\bar\rho}}\\ \alpha_{ij} =- 2\textcolor{red}{\bar \rho {\bar \Delta }^2 |\widetilde S|}(\widetilde {S_{ij}} - \frac{1}{3}\widetilde {S_{kk}}\delta _{ij})\\ \beta_{ij}=-2\widehat{\bar\rho} \widehat \Delta^2|\ddot S|(\ddot S_{ij}-\frac{1}{3}\ddot S_{kk}\delta_{ij})

Lij=Tijτij^L_{ij} = T_{ij} - \widehat {\tau _{ij}},我们可以得到:

Lij13Lkk^δij=Cs2βijCs2αij^=Cs2Mijwhere:Mij=βijαij^{L_{ij}} - \frac{1}{3}\widehat {L_{kk}}{\delta _{ij}} = C_s^2{\beta _{ij}} - C_s^2\widehat {\alpha _{ij}} = C_s^2{M_{ij}}\\ where: M_{ij}=\beta_{ij}-\widehat {\alpha _{ij}}

方程组是超定的,因此上式应该看作是在某种平均意义上成立。

Martín.2000采用Lilly.1992提出的最小二乘法,使式(15)的误差最小化,同时采用对均匀方向平均的形式:

Cs2=LijMijMklMkl , CI=Lkkβα^ , β=2Δ^2ρ^S¨2C_s^2 = \frac{\left\langle {L_{ij}{M_{ij}}} \right\rangle }{\left\langle {M_{kl}{M_{kl}}} \right\rangle }{\text{ }},{\text{ }}{C_I} = \frac{\left\langle {L_{kk}} \right\rangle }{\left\langle {\beta - \widehat \alpha } \right\rangle }{\text{ }},{\text{ }}\beta = 2{\widehat \Delta ^2}\widehat {\overline \rho }|\ddot S|^2

这就是可压 Smagorinsky 模型中需要的。

τij13τkkδij=2ρˉ(CsΔˉ)2S~(Sij~13Skk~δij)=Cs2αijτkk=CI2ρˉΔˉ2S~2=CIα\begin{array}{l} \tau _{ij} - \frac{1}{3}\tau _{kk}\delta _{ij} = - 2\textcolor{red}{\bar \rho (C_s{\bar \Delta })^2 |\widetilde S|}(\widetilde {S_{ij}} - \frac{1}{3}\widetilde {S_{kk}}\delta _{ij}) = C_s^2\alpha _{ij}\\ \tau _{kk} = C_I 2\bar \rho {\bar \Delta }^2|\widetilde S{|^2} = {C_I}\alpha \end{array}

滤波宽度比Δ^/Δ\widehat \Delta /\overline \Delta 是动态模型中唯一需要给定的参数。

文章作者: Yan Zhang
文章链接: https://openfoam.top/dynLES/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 OpenFOAM 成长之路
微信打赏给博主更多动力吧~