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

LES控制方程

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

如果寫成矩陣形式，引用李東岳的網站：


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

Boussinesq假定

具體推導可參照李東岳的博文，這里給出不可壓時的結果：

該方程描述了亞格子應力${\tau }_{ij}$和${\nu }_{sgs}$、$k_{sgs}$的關系。

Smagorinsky模型

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

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的代碼中，先計算了$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都與上面的結果一一對應。
再計算${\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()在各求解器中都會有的，在計算完這一時間步之后，修正湍流粘度系數，作下一時間步所用。這就是OpenFOAM植入湍流模型的邏輯：在不改變其控制方程的情況下，只修改湍流粘度系數，就能達到計算湍流的目的。
參數$C_k$在Lilly（1967）的文章中是0.094，和OpenFOAM中的對應。而$C_e$的來源尚不清楚，博文說$C_e=1.048$是因為湍流模型基于GenEddyVisc model（一般渦粘性假設？），如有誤，請加qq1019003721進行勘誤，謝謝！

參考文獻

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.

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。