一種結合壁面模型模擬湍流流動的間斷界面浸入邊界法

張仕釗，吳杰，楊德武，叢歆雨

南京航空航天大學 航空學院 空氣動力學系，南京 210016

近年來，一種名為浸入邊界法（IBM）的處理方法引起了人們的關注，其主要用于解決復雜外形網格生成困難和動邊界問題［1-2］。作為一種不要求邊界和網格重合的方法，IBM 使用簡單的、包含邊界的笛卡爾網格進行流場模擬。對于動邊界問題，IBM 使控制方程能夠在一個固定的網格上得到求解，同時可以解決與網格更新相關的關鍵問題。

IBM 最初由Peskin［3］在20 世紀70 年代提出，用于模擬心臟中的血液流動。在Peskin 的開創性工作之后，許多科學家為提高IBM 的準確性和效率做出了貢獻。Huang 和Tian［4］介紹了IBM的基礎知識，并評估了最新的進展。Cui 等［5］提出IBM 可以分為基于流體-結構界面表示的擴散界面法和間斷界面法兩類。

間斷界面IBM 是一種離散力法，與擴散界面IBM 不同，它能精確識別物面邊界，避免了界面模糊的情況。其主要特點是區分出邊界內部、邊界外部、以及邊界節點，并采用插值方法計算邊界節點的流動參數。間斷界面IBM 主要包括單元切割法（Cut-cell Method）［6-7］、虛擬單元法（Ghost-cell Method）［8］、笛卡爾浸入邊界法（Car‐tesian Immersed Boundary Method）［5，9］。在這些方法中，單元切割法得到的界面最為清晰，這是因為在單元切割法中，邊界面將網格節點劃分為固相和液相2 個區域，因此，流體的質量、動量、能量守恒定律都能被嚴格地執行。然而，在模擬運動物體的過程中，單元切割法的網格節點涉及復雜的重塑過程，這給模擬帶來了困難。虛擬單元法需要修改的邊界節點位于物體內部，該邊界節點稱為Ghost-cell，虛擬單元法選擇在這些點上重建流場變量。笛卡爾浸入邊界法類似于虛擬單元法，與之不同的是，該方法的邊界節點位于流場內部，邊界節點稱為IB 點或邊界點，重構這些點處的流動變量并加強邊界條件可以保證流場的準確性。

自然界中大部分流動都是湍流，盡管可以使用直接數值模擬（Direct Numerical Simulation，DNS）對湍流進行模擬，但由于該方法需要足夠小的空間和時間分辨率，現階段很難使用DNS 模擬較高雷諾數的湍流流動。一種替代方法是使用湍流模型［10］降低數值模擬對網格和時間步長的要求，以這種近似的方法來描述湍流。湍流模型的使用降低了數值模擬對網格的要求，但對于高雷諾數湍流模擬，在壁面附近仍然需要很高的分辨率。由于IBM 大多采用笛卡爾網格，無法做到像貼體網格那樣只在物面法向細化網格，滿足高雷諾數湍流模擬在物面附近所需要的高分辨率網格的要求，因此大部分IBM 的研究都停留在層流或者無黏流動［11］。為了降低高雷諾數湍流模擬在物面附近的網格要求，可以在IBM 中引入壁面模型［12］。壁面模型提供了一組以壁面剪應力的形式存在的邊界條件，從而允許在近壁區域存在較為粗糙的網格。壁面模型的引入使得IBM 在模擬高雷諾數湍流流動時，能保持本身優點且降低對物面附近網格細化的要求，避免網格細化帶來的計算效率下降的問題，這拓展了IBM在高雷諾數湍流中的應用［13］。

近年來，國內對IBM 的研究和應用與日俱增。李旭等［14］對Goldstein 提出的反饋力IBM 進行了新的思考，改進了其對力源項的計算，拓展了該方法的使用范圍。胡國暾等［15］使用IBM 求解振蕩轉子葉片的快速計算模型，從而避免了傳統方法中由于需要不斷重構貼體網格造成的數值模擬的復雜性。陳浩等［16］發展了基于IBM 思想的虛擬單元重構技術，構造了高保真的非貼體笛卡爾網格邊界條件，開展了該技術在小展弦比飛翼布局低速流動問題中的應用研究。唐志共等［17］從IBM 出發，結合虛擬鏡像對稱方法和曲率修正技術進行黏性物面邊界條件的處理，同時發展了多值點問題的處理技術，提出了一種在笛卡爾網格下可有效模擬黏性物面邊界條件的方法。

本文根據笛卡爾浸入邊界法的思想，將邊界節點放置在流場中，通過計算代表物面距離的符號向量場（Signed Distance Field，SDF）來區分出固體節點、流場節點以及邊界節點，SDF 的結果也可以用于k-ωSST（Shear Stress Transfer）湍流模型。對于Xu 和Liu［18］提到的由于邊界節點與壁面距離不規律導致的壁面剪切力計算不光滑的問題，本文采用壁面距離一致的鏡像點（Im‐age Point，IP）計算切應力uτ和湍動能耗散率ω，并通過Knopp 等［19］提出的壁面模型計算出邊界節點的切向速度，避免了這個問題。

1 數值方法

1.1 控制方程

非定常不可壓Navier-Stokes（N-S）方程為

式中：U為速度矢量；p為壓強；ρ為密度；μ為黏度；t為時間。

不可壓縮流體的連續性方程為

對式（1）取散度并結合式（2），可以得到壓力泊松方程

1.2 雷諾平均N-S 方程

采用雷諾平均N-S 方程（RANS）求解湍流問題，RANS 的思路是將流場變量分解為時均分量和脈動分量，即

將（4）～式（7）式代入式（1）、式（2），可得

由于速度分量u′、v′、w′在相關性中可互換，雷諾應力張量只有6 個獨立分量，其中u′，v′，w′分別為速度在x、y、z這3 個方向上的脈動分量。為了封閉方程，需要使用湍流模型計算雷諾應力張量，基于渦黏度法［10］，可以將式（8）轉化為

式中：μt為湍流黏度。式（11）與式（1）的唯一區別在于耗散項系數由μ轉變為μ+μt。

1.3 湍流模型

選擇兩方程模型k-ωSST 湍流模型計算μt，k-ωSST 是目前工業界較為流行的湍流模型。湍動能k和湍動能耗散率ω的輸運方程以微分形式［20］給出

式中：Pk為Faver 平均湍流應力張量和應變律張量的積；σk、σω、γ、β、β*為模型相關的常數；νt為湍流運動黏度；curl（U）為速度的旋度。由于k-ωSST 模型結合了k-ε和k-ω湍流模型，因此k-ωSST 湍流模型所用常數σk、σω、γ、β，通過k-ε和kω湍流模型所用的常數計算得到，計算方式為

式中：Φ為k-ωSST 模型經驗常數；Φ1、Φ2為k-ε模型經驗常數。

k-ω模型經驗常數為

k-ε模型經驗常數為

輔助函數F1表示為

式中：d為與壁面最近的網格單元中心與壁面的距離。

輔助函數F2表示為

模型常數為

1.4 壁面模型

使用的壁面模型是基于壁面和對數層外邊緣之間的流動由湍流邊界層（Turbulent Bound‐ary Layer，TBL）控制的假設得出，其中，二維TBL 方程為

式中：y為壁面法向方向；u為切向速度；S為非定常、對流及壓力梯度項之和。

采用Knopp 等［19］提出的一種與k-ωSST 湍流模型一致性的壁面函數，一致性代表壁面函數的解與TBL 方程指定的外邊界的位置無關，壁面函數被描述為

使用Spalding［21］提出的擬合公式來描述FSp

FRei，m使用Reichardt 的壁面定律

式中：κ為模型常數；y為第1 層網格高度；u為切向速度；uτ為壁面切應力；ν為運動黏度；y+、u+為表征近壁面流場信息的無量綱量，計算公式為

聯立式（22）～式（27）即可求得壁面切應力uτ。

2 浸入邊界法

2.1 邊界處理

由于本文采用間斷界面IBM，為了實施邊界條件，需要對網格節點進行區分。將固體邊界內部節點標記為固體點（Solid Points，SP），固體邊界外部為流場，在流場中根據與邊界的距離劃分出需要進行處理的邊界點（IB Points，IB），以及不需要特殊處理流場點（Fluid Points，FP）。

本文根據網格點與壁面的距離計算出SDF，對于SDF ≤0 的點識別為SP（圖1 的紅色區域），對于0

圖1 根據SDF 的單元分區及IP、IB 各點的相對位置Fig.1 Distinguish mesh cell according to SDF and rela‐tive locations of IP and IB points

本文研究發現，h1、h2的選取對模擬結果稍有影響。h1偏小會影響穩定性，h1偏大會影響精度；而h2偏小會導致采樣點分布不均，進而影響IP 的計算精度，h2偏大會因為距離IB 距離太遠而導致影響IB 的計算精度。因此，對于不同的流動問題需要不同的參數，本文參數的選取標準為

式中：Hmin為最小網格單元長度。本研究均采用直角網格。

2.2 IB 修正

P、B、W點的位置如圖1 所示，使用反距離加權插值算法，使用已知的流場變量計算出P點處流場變量的值

式中：R為以P點為圓心的插值模板的半徑，本文取R=2Hmin；ri為插值模板中插值點與P點的距離；?i為插值點的流場變量。

B點的法向速度UB，normal可以通過P點和W點法向速度線性插值求得

式中：下標normal、tangential 分別表示法向、切向分量；UP為P點速度，由式（29）計算得到；UW，normal為給定的Dirichlet 邊界條件的法向分量，hB為B點到壁面的距離。

若物面附近網格足夠密，即y+≤1，那么B點的切向速度可以認為同法向速度一樣呈線性分布。為了減少對網格的需求，本文使用壁面函數計算切向速度，計算流程如下：

1）根據P點的切向速度計算出壁面切應力uτ，詳細的計算流程可參見文獻［22］。

2）根據（1）得到的uτ計算出B點處的。

4）由式（27）可得B點的切向速度UB，tangential。

B點的壓強可由給定的Neumann 邊界條件和P點的壓強求得，計算公式為

式中：pB、pP分別為B、P點處的壓強。給定0，所以式（31）可簡化為pB=pP。式（30）、式（31）對于Dirichlet 和Neumann 邊界條件其他流場變量的修改也是通用的。

由于本文使用的是k-ωSST 湍流模型，還需要確定IB 處湍流變量的值。Kalitzin 等［23］給出了黏性層和對數層中k-ω湍流模型的近似解

式中：β=0.075；κ=0.41；Cμ=0.09；y為壁面距離。通常來講，位于過渡層的ω通過插值近似，本文采用Knopp 等［19］提出的方法計算

采用如下混合公式和漸近關系

使用式（32）～式（34）計算IB 的湍動能耗散率ωIB時，令y=h2而不是直接使用IB 的壁面距離hIB。如圖2 所示，一般情況下IB 與壁面的距離hIB不是平滑變化的，這會導致相鄰的IB 計算出的ωIB相差較大，從而在一定程度上導致計算難以收斂。

圖2 壁面附近hIB的分布情況Fig.2 Distribution of hIB near the wall

對于零壓力梯度邊界條件，式（21）中的S近似為0，可得式（35），無量綱化得到式（36）。式（37）為Spalding 提出的y+與u+的擬合公式［20］，將式（37）代入式（36）可得IB 的μt。

由μ、ω可得IB 的湍動能k［19］值為

2.3 數值模擬流程

綜上所述，當前算法從時間步n到n+1 的求解過程如下：

1）通過式（29）計算出IP 處 的速度U和壓強p。

2）通過式（30）～式（34）計算IB 處的U和p，并代入控制方程離散矩陣。

3）通過OpenFOAM 求解流場控制方程離散矩陣。

4）重復步驟1～3，直到U和p收斂。

5）計算IB 點的湍流變量ω、k、μt，作為求解k-ωSST 湍流模型的邊界條件，并代入湍流模型離散矩陣。

6）通過OpenFOAM 求解k-ωSST 湍流模型。

3 數值算例

3.1 平板繞流

由于湍流平板邊界層具有理論上的近似解，可以用來驗證求解器解析近壁面流動的準確性。采用近壁面網格密度不同的2 套網格，分別模擬雷諾數Re=1×106，7×105情況下的二維平板湍流邊界層問題，并與理論近似解進行對比。網格的y+基于IP 與物面距離計算，即y+=h2uτ/ν。

結果表明，隨著y+的減少，表面摩擦力系數Cf逐漸貼近理論值。在Re=1×106、y+=72 時，Cf的數值模擬結果和理論值基本一致。Cf的理論值采用冪律公式計算

式中：x為到平板前緣的距離；c為平板長度；u∞為遠前方來流速度。

由圖3 可以看出，除了平板前緣附近稍有偏差之外，其他位置的表面摩擦力均與理論值吻合。前緣的誤差也出現在其他的模擬中，例如Kalitizin 和Iaccarino［24］、Lee 和Ruffin［25］、Pu 和Zhu［11］等的模擬分析，這可能是因為前緣流場變化劇烈，網格分辨率不足以解析前緣流動導致的。

圖3 Re=1×106時使用不同y+網格得到的Cf與理論值對比Fig.3 Comparison of Cf between simulated and theoreti‐cal values using different y+ mesh at Re=1×106

圖4 展示了Re=7×105時平板表面Cf分布，它與冪律公式的計算結果基本一致，體現了本文方法對于一定雷諾數范圍的湍流平板問題的模擬較為準確。

圖4 Re=7×105時的Cf與理論值對比Fig.4 Comparison of Cf between simulated and theoreti‐cal values at Re=7×105

圖5 展示了距離前緣x/c=0.9 處x方向的速度型，可以看出基本復合對數律。對數律公式為

圖5 Re=1×106時在x/c=0.9 處的速度型Fig.5 Velocity profiles at position of x/c=0.9 for Re=1×106

3.2 NACA0012 翼型繞流

為了驗證本文方法的魯棒性，模擬一個不規則物體的繞流問題。選用NACA0012翼型，自由來流雷諾數Re=1×106，攻角α=0°，計算域為20cairfoil×20cairfoil，其中cairfoil為翼型弦長（見圖6）。采用網格密度不同的2 套網格，網格最小尺寸分別是3.8×10?4cairfoil、7.6×10?4cairfoil。這里采用Nguyen等［24］發展的二維網格加密方法對壁面附近的區域進行加密，以確保在整體網格量沒有明顯增加的前提下使得壁面附近的網格加密，以達到較小的y+值。

圖6 NACA0012 翼型周圍網格Fig.6 Grid around NACA0012 airfoil

速度、壓力系數Cp云圖如圖7 所示，在攻角為0°的情況下，壓強和速度呈對稱分布，在翼型前緣可以看到明顯的加速，而在1/4 翼型弦長之后速度變化逐漸不明顯，對應Cf在翼型表面的分布。壓強除了在翼型前端和末端為正值，上下表面附近的壓強都為負值。

圖7 Re=1×106、α=0°時的壓力系數云圖及無量綱化的x 方向速度云圖Fig.7 Pressure coefficient contour and dimensionless x-direction velocity contour according to inlet ve‐locity for Re=1×106 and α=0°

圖8 展示了Re=1×106、α=0°時，2 套不同y+網格下用IBM 模擬得到的Cp、Cf在翼型的表面分布情況，以及與貼體網格的結果對比。可以看到，IBM 的Cp與貼體網格的結果幾乎完全吻合。而Cf在前緣附近的結果有些許偏差，這可能是由于前緣流場變化劇烈，現有網格密度不足以完全解析流動導致的數值模擬誤差，而本算例只在物體附近加密而非全局網格加密，對數值模擬結果的準確性也會造成一定影響。對比2 套不同網格得到的Cf，可以看出前緣偏差隨著網格的加密而減小，證明本方法模擬不規則物體繞流問題的準確性。

圖8 表面摩擦力系數和壓力系數Fig.8 Skin friction coefficient and pressure coefficient

3.3 Buice 二維擴散管

為了驗證本文方法求解分離流問題的有效性，模擬了一個二維非對稱擴散器流動問題，外形為美國航空航天局（NASA）提供的模型［26］。根據入口高度和速度計算出的雷諾數為Re=2×104，入口高度H=1 m，入口速度大小Uin=0.2 m/s。

不同位置的水平速度分布與實驗數據［27］和貼體網格結果對比如圖9 所示，圖中的y坐標根據入口高度無量綱化，x坐標根據入口速度和剖面位置無量綱化。IBM 使用了y+=63，112 這2 套不同網格，可以看出得到的數值模擬結果和實驗以及貼體網格結果基本一致。y+=112 時得到的分離點較為靠后；隨著網格加密，y+=63 時的分離點位置和分離區尺寸與實驗結果、貼體網格計算結果基本一致。圖10 展示了y+=63 使用入口速度無量綱的速度云圖及流線圖，其中上半部分是IBM 的結果，下半部分是貼體網格的結果，可以看出兩者基本一致。

圖9 不同位置的水平速度分布Fig.9 Horizontal velocity profiles at different locations

圖10 速度云圖和流線圖Fig.10 Velocity contours and streamlines

3.4 二維圓柱繞流

圓柱繞流問題一直是測試求解器精度的標準算例，因此本小節模擬了二維圓柱繞流問題。其中，圓柱直徑d=1 m，基于圓柱直徑的雷諾數分別為Re=63 100，126 000，計算域45d×30d，速度入口條件為均勻來流，壓力入口條件為零梯度，速度出口條件為零梯度邊界條件，壓力出口給定值設置為0，上下邊界設置為自由流邊界條件。網格單元最小尺寸為0.003d。由于本節主要關注圓柱的受力情況，為了節省計算資源，只對圓柱附近的網格進行加密，網格如圖11 所示。

圖11 圓柱繞流計算網格Fig.11 Computational grid for cylinder

圖12、13 展示了在Re=63 100，126 000 下圓柱繞流的阻力系數CD和升力系數CL隨時間變化的曲線，表1 為平均阻力系數、平均升力系數、升力系數均方根CL，rms、斯特勞哈爾數St實驗值和數值模擬結果［28-31］對比。可以看出，本文的數值模擬結果與文獻結果基本一致。但是，與Yeon 等［29］使用LES 得到的結果對比，CL、St稍大，這可能是由于LES 和RANS 對湍流模化策略不同，以及本文只對近壁面附近的網格加密而沒有對尾流區加密導致誤差產生。CD、CL、St的計算公式為

表1 圓柱繞流Re=63 100，126 000 時的、CL，rms、StTable 1 ，CL，rms and St for flow around cylinder at Re=63 100，126 000

表1 圓柱繞流Re=63 100，126 000 時的、CL，rms、StTable 1 ，CL，rms and St for flow around cylinder at Re=63 100，126 000

圖12 圓柱繞流CD隨時間變化曲線Fig.12 Time histories of drag coefficients of flow around cylinder

圖13 圓柱繞流CL隨時間變化曲線Fig.13 Time histories of lift coefficients of flow around cylinder

式中：D、L分別為圓柱的阻力、升力；U為給定的入口來流速度大小；fs為渦脫落頻率，周期T=1/fs。

圖14 給出了Re=126 000 時的無量綱速度、無量綱渦量和壓強系數云圖，其中速度云圖用入口速度無量綱化，渦量用入口速度Uin與圓柱直徑d無量綱化，即Wzd/Uin，其中Wz為z方向的渦量。IBM 的結果與文獻［30］的結果吻合較好，證明了本文的IBM 可以很好地捕獲臨界流狀態下重要的全局流特征。圖15 用流線展示了Re=126 000 時渦在圓柱表面產生、發展、壯大、脫落的全過程。

圖14 Re=126 000 時圓柱繞流計算結果Fig.14 Calculation results of flow around a cylinder at Re=126 000

3.5 二維方柱繞流

為了進一步測試本求解器求解鈍體繞流問題的精度，模擬了二維方柱的繞流問題。方柱切面為一正方形，邊長d=1 m，基于邊長的雷諾數為Re=22 000，計算域大小為45d×30d，速度入口條件為均勻來流，壓力入口條件為零梯度，速度出口條件為零梯度邊界條件，壓力出口給定值設置為0，上下邊界設置為自由流邊界條件。網格單元最小尺寸為0.006d。與3.5 節類似，這里也只對方柱附近的網格進行加密，網格如圖16 所示。

圖16 方柱繞流計算網格Fig.16 Computational grid for square cylinder

圖17 展示了方柱繞流CD、CL隨時間變化的曲線，表2 為用本文方法計算得到的平均阻力系數、升力系數均方根CL，rms、施特魯哈爾數St與其他實驗和數值模擬結果［32-34］的對比，可以看到結果基本一致。圖18 給出了無量綱化的速度、無量綱渦量和壓力系數云圖，以及對應時刻的流線圖，其中速度云圖用入口速度Uin無量綱化，渦量用入口速度與圓柱直徑d無量綱化。可以看出本文的方法能很好地捕捉到壁面附近及尾渦的流場細節。

表2 方柱繞流的、CL，rms、StTable 2 ，CL，rms and St for flow around square cylinder

表2 方柱繞流的、CL，rms、StTable 2 ，CL，rms and St for flow around square cylinder

圖17 方柱繞流CD、CL隨時間變化曲線Fig.17 Time histories of drag and lift coefficients of flow around square cylinder

圖18 Re=22 000 時方柱繞流計算結果Fig.18 Calculation results of flow around square cylin‐der at Re=22 000

4 結論

提出了一種可用于模擬不可壓黏性流動的間斷界面浸入邊界法，并將該方法植入開源軟件OpenFOAM 中，該方法采用k-ωSST 湍流模型和壁面函數技術降低了湍流模擬時的網格分辨率和時間分辨率要求，為了解決浸入邊界法實施壁面函數時出現的物面距離變化劇烈問題，本文采用了IP 與物面的距離參與計算，這有利于減少因為IB 分布變化劇烈導致的壁面附近ω分布變化劇烈的情況，從而提高計算的穩定性。

在應用本文方法時，IB 層的厚度h1、IP 與壁面的距離h2的選取對數值模擬結果稍有影響，選取的值太小不利于計算的穩定性，太大會影響精度。這里給出的參考值是h1=1.5Hmin，h2=3Hmin，對于h1、h2如何影響數值模擬結果還要進一步研究。

為了驗證提出方法的有效性，本文進行了相關問題的數值模擬，包括高雷諾數湍流平板繞流、NACA0012 表面壓強和切應力分布、Buice 非對稱擴散管、圓柱繞流和方柱繞流。計算結果與文獻結果相比有較好的一致性，從而驗證了本文方法模擬湍流流動的有效性和準確性。未來工作可以在本文方法的基礎上，進一步拓展到高雷諾數可壓縮流動。