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

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

南京航空航天大學(xué) 航空學(xué)院 空氣動力學(xué)系，南京 210016

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

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

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

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

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

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

1 數(shù)值方法

1.1 控制方程

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

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

不可壓縮流體的連續(xù)性方程為

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

1.2 雷諾平均N-S 方程

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

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

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

式中：μt為湍流黏度。式（11）與式（1）的唯一區(qū)別在于耗散項系數(shù)由μ轉(zhuǎn)變?yōu)棣?μt。

1.3 湍流模型

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

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

式中：Φ為k-ωSST 模型經(jīng)驗常數(shù)；Φ1、Φ2為k-ε模型經(jīng)驗常數(shù)。

k-ω模型經(jīng)驗常數(shù)為

k-ε模型經(jīng)驗常數(shù)為

輔助函數(shù)F1表示為

式中：d為與壁面最近的網(wǎng)格單元中心與壁面的距離。

輔助函數(shù)F2表示為

模型常數(shù)為

1.4 壁面模型

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

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

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

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

FRei，m使用Reichardt 的壁面定律

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

聯(lián)立式（22）～式（27）即可求得壁面切應(yīng)力uτ。

2 浸入邊界法

2.1 邊界處理

由于本文采用間斷界面IBM，為了實施邊界條件，需要對網(wǎng)格節(jié)點進(jìn)行區(qū)分。將固體邊界內(nèi)部節(jié)點標(biāo)記為固體點（Solid Points，SP），固體邊界外部為流場，在流場中根據(jù)與邊界的距離劃分出需要進(jìn)行處理的邊界點（IB Points，IB），以及不需要特殊處理流場點（Fluid Points，F(xiàn)P）。

本文根據(jù)網(wǎng)格點與壁面的距離計算出SDF，對于SDF ≤0 的點識別為SP（圖1 的紅色區(qū)域），對于0

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

本文研究發(fā)現(xiàn)，h1、h2的選取對模擬結(jié)果稍有影響。h1偏小會影響穩(wěn)定性，h1偏大會影響精度；而h2偏小會導(dǎo)致采樣點分布不均，進(jìn)而影響IP 的計算精度，h2偏大會因為距離IB 距離太遠(yuǎn)而導(dǎo)致影響IB 的計算精度。因此，對于不同的流動問題需要不同的參數(shù)，本文參數(shù)的選取標(biāo)準(zhǔn)為

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

2.2 IB 修正

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

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

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

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

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

1）根據(jù)P點的切向速度計算出壁面切應(yīng)力uτ，詳細(xì)的計算流程可參見文獻(xiàn)［22］。

2）根據(jù)（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］給出了黏性層和對數(shù)層中k-ω湍流模型的近似解

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

采用如下混合公式和漸近關(guān)系

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

圖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 數(shù)值模擬流程

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

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

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

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

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

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

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

3 數(shù)值算例

3.1 平板繞流

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

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

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

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

圖3 Re=1×106時使用不同y+網(wǎng)格得到的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分布，它與冪律公式的計算結(jié)果基本一致，體現(xiàn)了本文方法對于一定雷諾數(shù)范圍的湍流平板問題的模擬較為準(zhǔn)確。

圖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方向的速度型，可以看出基本復(fù)合對數(shù)律。對數(shù)律公式為

圖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 翼型繞流

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

圖6 NACA0012 翼型周圍網(wǎng)格Fig.6 Grid around NACA0012 airfoil

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

圖7 Re=1×106、α=0°時的壓力系數(shù)云圖及無量綱化的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+網(wǎng)格下用IBM 模擬得到的Cp、Cf在翼型的表面分布情況，以及與貼體網(wǎng)格的結(jié)果對比?？梢钥吹剑琁BM 的Cp與貼體網(wǎng)格的結(jié)果幾乎完全吻合。而Cf在前緣附近的結(jié)果有些許偏差，這可能是由于前緣流場變化劇烈，現(xiàn)有網(wǎng)格密度不足以完全解析流動導(dǎo)致的數(shù)值模擬誤差，而本算例只在物體附近加密而非全局網(wǎng)格加密，對數(shù)值模擬結(jié)果的準(zhǔn)確性也會造成一定影響。對比2 套不同網(wǎng)格得到的Cf，可以看出前緣偏差隨著網(wǎng)格的加密而減小，證明本方法模擬不規(guī)則物體繞流問題的準(zhǔn)確性。

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

3.3 Buice 二維擴(kuò)散管

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

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

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

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

3.4 二維圓柱繞流

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

圖11 圓柱繞流計算網(wǎng)格Fig.11 Computational grid for cylinder

圖12、13 展示了在Re=63 100，126 000 下圓柱繞流的阻力系數(shù)CD和升力系數(shù)CL隨時間變化的曲線，表1 為平均阻力系數(shù)、平均升力系數(shù)、升力系數(shù)均方根CL，rms、斯特勞哈爾數(shù)St實驗值和數(shù)值模擬結(jié)果［28-31］對比?？梢钥闯?，本文的數(shù)值模擬結(jié)果與文獻(xiàn)結(jié)果基本一致。但是，與Yeon 等［29］使用LES 得到的結(jié)果對比，CL、St稍大，這可能是由于LES 和RANS 對湍流模化策略不同，以及本文只對近壁面附近的網(wǎng)格加密而沒有對尾流區(qū)加密導(dǎo)致誤差產(chǎn)生。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 時的無量綱速度、無量綱渦量和壓強系數(shù)云圖，其中速度云圖用入口速度無量綱化，渦量用入口速度Uin與圓柱直徑d無量綱化，即Wzd/Uin，其中Wz為z方向的渦量。IBM 的結(jié)果與文獻(xiàn)［30］的結(jié)果吻合較好，證明了本文的IBM 可以很好地捕獲臨界流狀態(tài)下重要的全局流特征。圖15 用流線展示了Re=126 000 時渦在圓柱表面產(chǎn)生、發(fā)展、壯大、脫落的全過程。

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

3.5 二維方柱繞流

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

圖16 方柱繞流計算網(wǎng)格Fig.16 Computational grid for square cylinder

圖17 展示了方柱繞流CD、CL隨時間變化的曲線，表2 為用本文方法計算得到的平均阻力系數(shù)、升力系數(shù)均方根CL，rms、施特魯哈爾數(shù)St與其他實驗和數(shù)值模擬結(jié)果［32-34］的對比，可以看到結(jié)果基本一致。圖18 給出了無量綱化的速度、無量綱渦量和壓力系數(shù)云圖，以及對應(yīng)時刻的流線圖，其中速度云圖用入口速度Uin無量綱化，渦量用入口速度與圓柱直徑d無量綱化?？梢钥闯霰疚牡姆椒芎芎玫夭蹲降奖诿娓浇拔矞u的流場細(xì)節(jié)。

表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 時方柱繞流計算結(jié)果Fig.18 Calculation results of flow around square cylin‐der at Re=22 000

4 結(jié)論

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

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

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