本人南航CFD研究生，歡迎加qq：1019003721互相學(xué)習(xí)討論！
本文將介紹OpenFOAM里面的大渦模擬相關(guān)代碼。起因是最近學(xué)習(xí)OpenFOAM中的大渦模擬，在一篇OpenFOAM中的LES湍流模型及其植入的專欄中看到，但是里面的公式不知道是什么格式，雖然知道一些變量，但一時(shí)間也無(wú)法弄清楚。因此一邊學(xué)習(xí)張兆順教授的湍流大渦數(shù)值模擬的理論與應(yīng)用及相關(guān)文獻(xiàn)，一邊看代碼校對(duì)。（相關(guān)文獻(xiàn)資料下載）在這期間，看到CFD中文網(wǎng)的李東岳教授發(fā)的CFD中的LES湍流模型。這里面講得非常詳細(xì)。我將結(jié)合上面的內(nèi)容和OpenFOAM對(duì)應(yīng)的Smagorinsky.C中的代碼盡可能仔細(xì)地講解。

LES控制方程

OpenFOAM內(nèi)含有非常豐富的操作符和算子等。只要確定了控制方程，就能在求解器里自定義要程序解的方程了。在大渦模擬中，省去推導(dǎo)過(guò)程，控制方程一般如下形式（張兆順，2007）：

如果寫成矩陣形式，引用李東岳的網(wǎng)站：


先看上面的分量形式。在這個(gè)控制方程中，要想在OpenFOAM里求解流場(chǎng)變量，仍需要知道${\tau }_{kk}/3$和${\nu }_t$才能封閉方程組。在這里，${\nu }_t$即為亞格子渦粘系數(shù)，它在OpenFOAM中是nuSgs。而正應(yīng)力${\tau }_{kk}$則可以用亞格子動(dòng)能代替：$k_{sgs}={\tau }_{kk}/2$。$k_{sgs}$在OpenFOAM中是k。

Boussinesq假定

具體推導(dǎo)可參照李東岳的博文，這里給出不可壓時(shí)的結(jié)果：

該方程描述了亞格子應(yīng)力${\tau }_{ij}$和${\nu }_{sgs}$、$k_{sgs}$的關(guān)系。

Smagorinsky模型

在1967年Lilly給出了${\tau }_{kk}/3$的表達(dá)式（相關(guān)文獻(xiàn)）：
$$\frac{{\tau }_{kk}}{3}=\frac{2}{3}{\nu }_{sgs}^2/(C_k△{)}^2$$
而正應(yīng)力${\tau }_{kk}$又可以通過(guò)$k_{sgs}={\tau }_{kk}/2$來(lái)替換，最后得到${\nu }_{sgs}$和$k_{sgs}$的關(guān)系：
$${\nu }_{sgs}=C_k△\sqrt{k_{s}gs}$$
由此，方程的封閉還差最后一個(gè)變量$k_{sgs}$。Smagorinsky提出模型：
$${\bar{S}}_{ij}:{\tau }_{ij}+C_e{k}_{sgs}^{\frac{3}{2}}/△=0$$
通過(guò)代入上述各式，即可求解$k_{sgs}$。過(guò)程可參考CFD中文網(wǎng)中一二發(fā)的帖子。
Ps:運(yùn)算符“：”是雙重內(nèi)積，即兩張量先作矩陣乘法運(yùn)算（點(diǎn)積），再進(jìn)行縮并（Contraction）。因?yàn)閰⑴c運(yùn)算的兩個(gè)張量都是二階，點(diǎn)積后仍為二階，但縮并后降二階，最后是零階，即一個(gè)數(shù)。最后這個(gè)方程是一個(gè)一元二次方程，按公式法求解即可。
在此引用一二發(fā)的帖子中的結(jié)果：

OpenFOAM中代碼

在OpenFOAM中，公式與上文所述一致，只是所用的變量名字不同：

\verbatimB = 2/3*k*I - 2*nuSgs*dev(D)whereD = symm(grad(U));k from D:B + Ce*k^3/2/delta = 0nuSgs = Ck*sqrt(k)*delta\endverbatim

B就是${\tau }_{ij}$，D就是$S_{ij}$。
在Smagorinsky.C的代碼中，先計(jì)算了$k_{sgs}$：

template<class BasicTurbulenceModel> tmp<volScalarField> Smagorinsky<BasicTurbulenceModel>::k (const tmp<volTensorField>& gradU ) const {volSymmTensorField D(symm(gradU));volScalarField a(this->Ce_/this->delta());volScalarField b((2.0/3.0)*tr(D));volScalarField c(2*Ck_*this->delta()*(dev(D) && D));return volScalarField::New(IOobject::groupName("k", this->alphaRhoPhi_.group()),sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a))); }

代碼中abc都與上面的結(jié)果一一對(duì)應(yīng)。
再計(jì)算${\nu }_{sgs}$并修正控制方程中的nut：

template<class BasicTurbulenceModel> void Smagorinsky<BasicTurbulenceModel>::correctNut() {volScalarField k(this->k(fvc::grad(this->U_)));this->nut_ = Ck_*this->delta()*sqrt(k);this->nut_.correctBoundaryConditions();fv::options::New(this->mesh_).correct(this->nut_);BasicTurbulenceModel::correctNut(); }

correctNut()在各求解器中都會(huì)有的，在計(jì)算完這一時(shí)間步之后，修正湍流粘度系數(shù)，作下一時(shí)間步所用。這就是OpenFOAM植入湍流模型的邏輯：在不改變其控制方程的情況下，只修改湍流粘度系數(shù)，就能達(dá)到計(jì)算湍流的目的。
參數(shù)$C_k$在Lilly（1967）的文章中是0.094，和OpenFOAM中的對(duì)應(yīng)。而$C_e$的來(lái)源尚不清楚，博文說(shuō)$C_e=1.048$是因?yàn)橥牧髂Ｐ突贕enEddyVisc model（一般渦粘性假設(shè)？），如有誤，請(qǐng)加qq1019003721進(jìn)行勘誤，謝謝！

參考文獻(xiàn)

Smagorinsky J. General circulation experiments with the primitive
equations: I. the basic experiment[J]. Monthly weather review, 1963,
91(3): 99-164.

Smagorinsky, J. , Manabe, S. , & Holloway, J. L. . (1963).
“numerical results from a nine-layer general circulation model,”.
monthly weather review, 91(12), 263-270.

Deardorff, J. W. . (1970). A numerical study of three-dimensional
turbulent channel flow at large reynolds number. Journal of Fluid
Mechanics, 41.

LILLY, D. K. 1967 Proceedings of the IBM Scientific Computing
Symposium on Environmental Sciencer, IBM Form no. 320-1951, 195.

張兆順, 崔桂香, 許春曉. 湍流大渦數(shù)值模擬的理論與應(yīng)用[M]. 清華大學(xué)出版社, 2008.
其中1235可以在我上傳的文件中下載。

總結(jié)

以上是生活随笔為你收集整理的OpenFOAM大涡模拟湍流模型之Smagorinsky模型代码详解的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。