羅 寅 張阿漫 朱永凱 徐珊珊

哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱 150001

結(jié)合亞格子模型的Lattice－Boltzmann算法

羅 寅 張阿漫 朱永凱 徐珊珊

哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱 150001

介紹一種亞格子模型與LB（Lattic－Boltzmann）方法相結(jié)合的算法，并用其對(duì)亞臨界雷諾數(shù)情況下的圓柱繞流問題進(jìn)行數(shù)值模擬，得到流場(chǎng)的流線圖以及圓柱表面的壓力系數(shù)。通過與實(shí)驗(yàn)值相比較，發(fā)現(xiàn)數(shù)值結(jié)果與實(shí)驗(yàn)值是比較吻合的，證明該算法可以用于計(jì)算雷諾數(shù)較大的流體問題，同時(shí)具有LB方法計(jì)算效率高的優(yōu)點(diǎn)。

格子Boltzmann方法；亞臨界雷諾數(shù)；圓柱繞流

1 引言

近年來，備受人們關(guān)注的LB方法通過利用粒子的簡單且規(guī)則統(tǒng)一的假想運(yùn)動(dòng)模型來代替復(fù)雜的真實(shí)微觀運(yùn)動(dòng)，從而達(dá)到實(shí)現(xiàn)模擬宏觀物理現(xiàn)象的目的。LB方法的前身是格子氣自動(dòng)機(jī)（Lattice Gas Automata，LGA［1］）。 LB 方法使用中觀的實(shí)數(shù)型變量—粒子分布函數(shù)替代了LGA的微觀的布爾型變量，從而不再需要具體計(jì)算微觀量，而是直接按Lattic－Boltzmann方程來模擬宏觀現(xiàn)象。LB方法應(yīng)用最廣泛的模型是單松弛模型［2］，它與分子動(dòng)理論中Boltzmann方程的BGK（Bhatnagar－Gross－Krook）算子類似，故這類模型又稱為“LBGK 模型”［3］。LB 方法與傳統(tǒng)的 CFD 方法相比，在問題的求解規(guī)模和求解速度方面，有著后者無法比擬的優(yōu)勢(shì)。在水力學(xué)問題、多相流、對(duì)流擴(kuò)散問題及求解偏微分方程等許多方面，LB方法都展示了廣闊的應(yīng)用前景，甚至傳統(tǒng)的CFD方法難以模擬的多孔介質(zhì)問題，LB方法也能從容處理。LB方法還能夠模擬一些沒有宏觀控制議程的問題。但是，和傳統(tǒng)方法相比，LB方法也有不足的地方，表現(xiàn)在所能求解的實(shí)際問題雷諾數(shù)范圍太低，一般只有低于103時(shí)，才有滿意的精度，而當(dāng)雷諾數(shù)進(jìn)一步提高時(shí)，計(jì)算容易失穩(wěn)。本文將傳統(tǒng)的亞格子模型成功地運(yùn)用到LB方法中，顯著地提高了LB方法的適用范圍。

2 LB方法簡介

在 LB 方法中，粒子分布函數(shù) fα（x，t）的演化方程可寫成如下形式：

式中，fα（x，t） 為沿 α 方向的粒子分布函數(shù)；c＝Δx／Δt，Δx、Δt分別為空間和時(shí)間步長， 可取Δx＝Δt，即 c＝1。

LB方法應(yīng)用較多的是二維九點(diǎn)模型（D2Q9）。

圖1 D2Q9模型速度張量圖

這里，α ＝0，…，b（b ＝8）；eα為粒子的運(yùn)動(dòng)方向；τ為馳譽(yù)時(shí)間；f（x，t）為沿 α 方向的粒子平衡分布函數(shù)。

將演化方程（1）按照 Chapman－Enskog展開，并進(jìn)行尺度粘合［4］可以得到流場(chǎng)控制方程，即連續(xù)性方程和N－S方程。在節(jié)點(diǎn)上的宏觀量（密度、速度和壓力）與粒子分布函數(shù)滿足方程：

D2Q9模型的馳譽(yù)時(shí)間與粘性系數(shù)滿足下面的關(guān)系：

3 格子Boltzmann亞格子模型

從數(shù)值計(jì)算角度看，LB方法可以看作是求解N－S方程的一種數(shù)值方法，同傳統(tǒng)的流體力學(xué)方法一樣，要直接計(jì)算雷諾數(shù)較高流體問題就意味需要降低網(wǎng)格尺度，這樣對(duì)計(jì)算機(jī)的硬件要求是很高的，相應(yīng)的計(jì)算時(shí)間也成倍增長，并且在LB方法中對(duì)馬赫數(shù)是有要求的，因此，流場(chǎng)速度也是受限制的。傳統(tǒng)的方法在求解湍流問題時(shí)將其分為兩種尺度，大尺度運(yùn)動(dòng)由N－S方程描述，小尺度運(yùn)動(dòng)的耗散作用及其對(duì)大尺度運(yùn)動(dòng)的反饋，則通過引入湍流模型來描述。亞格子模型便是其中一種。亞格子模型沿襲了湍流模式化方法，其中最早最常用的是氣象學(xué)家 Smagorinsky［5，6］在 1963 年提出的渦粘性模型。該模型假定用各向同性濾波器濾掉的小尺度脈動(dòng)是局部平衡的，假定渦粘性正比于亞格子尺度的特征長度和特征湍流速度，以此求出亞格子應(yīng)力。根據(jù)Smagorinsky提出的渦粘性模型可知：

為了將亞格子模型引入LB方法，假設(shè)粒子的碰撞過程僅與局部的粒子信息相關(guān)，保持模型中的平衡分布函數(shù)形式不變，其中的宏觀物理量均用過濾后的平均量替代。D2Q9模型中弛豫時(shí)間與粘性的關(guān)系變?yōu)椋?/p>

式中，粘性系νtotal分為物理動(dòng)力粘性ν0和渦粘性νt兩個(gè)部分。將亞格子模型渦粘性系數(shù)應(yīng)用到式（6）中可以得到：

本文中濾波尺度Δ＝δx。根據(jù)節(jié)點(diǎn)的分布函數(shù)可以很容易的確定應(yīng)變率張量Sij：

結(jié)合了亞格子模型的D2Q9模型演化方程最終形式為：

公式（9）與公式（1）有相同的形式，不同之處在于原模型中弛豫時(shí)間是固定值，本文通過引入亞格子模型對(duì)弛豫時(shí)間進(jìn)行了修正。在計(jì)算過程中，弛豫時(shí)間根據(jù)每步計(jì)算結(jié)果的變化而變化。這個(gè)特點(diǎn)說明，使用本文所引入的方法同使用原LB方法具有相同的計(jì)算流程，只是在每一個(gè)迭代時(shí)間步開始前，利用公式（7）完成對(duì)弛豫時(shí)間的修正即可。

4 圓柱繞流算例

為了證明本文方法的有效性，本文使用D2Q9模型，對(duì)亞雷諾數(shù)條件Re＝1.6×104圓柱繞流進(jìn)行數(shù)值模擬 流場(chǎng)區(qū)域結(jié)構(gòu)如圖2所示

眾所周知，對(duì)圓柱繞流問題而言，當(dāng)雷諾數(shù)超過一定數(shù)值時(shí)，在圓柱后面會(huì)產(chǎn)生周期性的旋渦，俗稱“卡門渦街”。圖3所示為計(jì)算結(jié)果穩(wěn)定后，截取圓柱附近的流場(chǎng)向量圖。圖4是對(duì)應(yīng)的流線圖，截取間隔時(shí)間為2 s。由該圖可以明顯地看到，在圓柱后面出現(xiàn)了旋渦，并且旋渦在柱體左右兩側(cè)交替地產(chǎn)生、發(fā)放。圓柱繞流的這種旋渦發(fā)放是周期性的。

隨著渦街的周期性發(fā)放， 在柱體表面會(huì)受到周期性的升力FL和阻力FD。在圓柱繞流計(jì)算中，參數(shù)Strouhal（簡稱St）數(shù)是表明旋渦脫落特性的相似準(zhǔn)則數(shù)，是升力頻率fs的無因次表達(dá)。在本文算例中，經(jīng)過計(jì)算得到St的值為0.21。而根據(jù)實(shí)驗(yàn)研究得到，在亞臨界雷諾數(shù)區(qū)范圍時(shí)，即30≤Re＜3 ×105時(shí)，St約等于 0.2， 兩者相對(duì)誤差是5%。流體流過圓柱的時(shí)候，圓柱會(huì)受到升力和阻力是因?yàn)樵趫A柱表面受對(duì)壓力的結(jié)果。圖6是圓柱表面無量綱穩(wěn)態(tài)壓力系數(shù)Cp計(jì)算值與實(shí)驗(yàn)值對(duì)比圖。

圖5 圓柱表面的Cp 值變化圖

從圖5可以看出，Cp的計(jì)算值與實(shí)驗(yàn)值是非常吻合的。另外，由于隨時(shí)間變化的渦激勵(lì)是一個(gè)平穩(wěn)隨機(jī)過程，因此通常采用無因次脈動(dòng)阻力系數(shù) CD′，脈動(dòng)升力系數(shù) CL′來表達(dá)脈動(dòng)阻力值、脈動(dòng)升力值。 本文中 CL′＝ 0.55，CD′＝ 0.1，這是符合實(shí)驗(yàn)值［6］的。

5 結(jié)論

本文將傳統(tǒng)的亞格子模型應(yīng)用于LB方法，用該算法模擬了亞臨界雷諾數(shù)條件下的圓柱繞流問題，并通過將數(shù)值模擬的結(jié)果與實(shí)驗(yàn)值對(duì)比，證明了這種方法的有效性。

［1］FRISCH U，HASSIACHER B，POMEAU Y.Lattice gas automata for the Navier－Stokes equation［J］.Physical Review Letters， 1986，56（14）：1505－1508.

［2］QIAN Y H，HUMIERES D，LALLEMAND P.Lattice BGK models for the Navier－Stokes Equation［J］.Europhys.Lett，1992，（6）：479－484.

［3］GUO Z L，SHI B C，WANG N C.Lattice BGK model for incompressible Navier－Stokes equation ［J］.J Stat Phys，2000（165）：288－306.

［4］SUCCI S， HUMIRES D D＇， QIAN Y H， el at.On the small－scale dynamical behavior of lattice BGK and lattice Bolzmann schemes ［J］.J.Sci.Computing，1993，8 （3）：219－230.

［5］SMAGORINSKY J.General circulation experiments with primitive equations［J］.Mon Weather Rev，1963（91）：99－165.

［6］SMAGORINSKY J.General circulation experiments with the primitive equations part I：the basic experiment［J］.Mon.Weather Rev.，1963， 91（1）：99－165.

Combination of the Lattice-Boltzmann Method and the Sub-grid Model

Luo Yin Zhang A－man Zhu Yong－kaiXu Shan－shan

College of Shipbuilding Engineering， Harbin Engineering University， Harbin 150001，China

An algorithm to combine the Lattice－Boltzmann （LB） method and the sub-grid model is introduced and is used to simulate the flow around a cylinder under subcritical Reynolds number.The obtained streamline of the flow field and the pressure coefficient around the cylinder are found to agree well with the data from tests.The agreement indicates that the discussed method suits those problems with relatively high Reynolds number and inherits the computation efficiency of the LB method.

Lattice Boltzmann method； subcritical Raynolds number； flow around cylinder

U661.73

A

1673－3185（2010）01－10－04

2009 - 06 - 26

羅寅（1985 － ） ，男，碩士研究生。研究方向：船舶與海洋結(jié)構(gòu)物性能與安全性。E-mail：luoyinlsq@ 163. com

張阿漫（1981 - ） ，男，副教授。研究方向：船舶與海洋工程結(jié)構(gòu)動(dòng)力學(xué)