本人南航CFD研究所學生，歡迎加qq：1019003721互相學習讨論！
本文将介紹OpenFOAM裡面的大渦模拟相關代碼。起因是最近學習OpenFOAM中的大渦模拟，在一篇OpenFOAM中的LES湍流模型及其植入的專欄中看到，但是裡面的公式不知道是什麼格式，雖然知道一些變量，但一時間也無法弄清楚。是以一邊學習張兆順教授的湍流大渦數值模拟的理論與應用及相關文獻，一邊看代碼校對。（相關文獻資料下載下傳）在這期間，看到CFD中文網的李東嶽教授發的CFD中的LES湍流模型。這裡面講得非常詳細。我将結合上面的内容和OpenFOAM對應的Smagorinsky.C中的代碼盡可能仔細地講解。
LES控制方程
OpenFOAM内含有非常豐富的操作符和算子等。隻要确定了控制方程，就能在求解器裡自定義要程式解的方程了。在大渦模拟中，省去推導過程，控制方程一般如下形式（張兆順，2007）：
如果寫成矩陣形式，引用李東嶽的網站：
先看上面的分量形式。在這個控制方程中，要想在OpenFOAM裡求解流場變量，仍需要知道 τ k k / 3 \tau_{kk}/3 τkk/3和 ν t \nu_t νt才能封閉方程組。在這裡， ν t \nu_t νt即為亞格子渦粘系數，它在OpenFOAM中是nuSgs。而正應力 τ k k \tau_{kk} τkk則可以用亞格子動能代替： k s g s = τ k k / 2 k_{sgs}=\tau_{kk}/2 ksgs=τkk/2。 k s g s k_{sgs} ksgs在OpenFOAM中是k。
Boussinesq假定
具體推導可參照李東嶽的博文，這裡給出不可壓時的結果：
該方程描述了亞格子應力 τ i j \tau_{ij} τij和 ν s g s \nu_{sgs} νsgs、 k s g s k_{sgs} ksgs的關系。
Smagorinsky模型
在1967年Lilly給出了 τ k k / 3 \tau_{kk}/3 τkk/3的表達式（相關文獻）：
τ k k 3 = 2 3 ν s g s 2 / ( C k △ ) 2 \frac{\tau_{kk}}{3}=\frac{2}{3}\nu_{sgs}^2/(C_k\triangle)^2 3τkk=32νsgs2/(Ck△)2
而正應力 τ k k \tau_{kk} τkk又可以通過 k s g s = τ k k / 2 k_{sgs}=\tau_{kk}/2 ksgs=τkk/2來替換，最後得到 ν s g s \nu_{sgs} νsgs和 k s g s k_{sgs} ksgs的關系：
ν s g s = C k △ k s g s \nu_{sgs}=C_k\triangle\sqrt{k_sgs} νsgs=Ck△ksgs
由此，方程的封閉還差最後一個變量 k s g s k_{sgs} ksgs。Smagorinsky提出模型：
S ˉ i j : τ i j + C e k s g s 3 2 / △ = 0 \bar{S}_{ij}:\tau_{ij}+C_ek_{sgs}^\frac{3}{2}/\triangle=0 Sˉij:τij+Ceksgs23/△=0
通過代入上述各式，即可求解 k s g s k_{sgs} ksgs。過程可參考CFD中文網中一二發的文章。
Ps:運算符“：”是雙重内積，即兩張量先作矩陣乘法運算（點積），再進行縮并（Contraction）。因為參與運算的兩個張量都是二階，點積後仍為二階，但縮并後降二階，最後是零階，即一個數。最後這個方程是一個一進制二次方程，按公式法求解即可。
在此引用一二發的文章中的結果：
OpenFOAM中代碼
在OpenFOAM中，公式與上文所述一緻，隻是所用的變量名字不同：
\verbatim
B = 2/3*k*I - 2*nuSgs*dev(D)
where
D = symm(grad(U));
k from D:B + Ce*k^3/2/delta = 0
nuSgs = Ck*sqrt(k)*delta
\endverbatim
B就是 τ i j \tau_{ij} τij，D就是 S i j S_{ij} Sij。
在Smagorinsky.C的代碼中，先計算了 k s g s k_{sgs} ksgs：
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都與上面的結果一一對應。
再計算 ν s g s \nu_{sgs} ν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 C_k Ck在Lilly（1967）的文章中是0.094，和OpenFOAM中的對應。而 C e C_e Ce的來源尚不清楚，博文說 C e = 1.048 C_e=1.048 Ce=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可以在我上傳的檔案中下載下傳。