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

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

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

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

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

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

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

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

1 數(shù)值方法

1.1 控制方程

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

式中：U為速度矢量；p為壓強(qiáng)；ρ為密度；μ為黏度；t為時(shí)間。

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

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

1.2 雷諾平均N-S 方程

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

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

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

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

1.3 湍流模型

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

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

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

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

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

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

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

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

模型常數(shù)為

1.4 壁面模型

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

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

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

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

FRei，m使用Reichardt 的壁面定律

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

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

2 浸入邊界法

2.1 邊界處理

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

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

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

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

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

2.2 IB 修正

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

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

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

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

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

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

2）根據(jù)（1）得到的uτ計(jì)算出B點(diǎn)處的。

4）由式（27）可得B點(diǎn)的切向速度UB，tangential。

B點(diǎn)的壓強(qiáng)可由給定的Neumann 邊界條件和P點(diǎn)的壓強(qiáng)求得，計(jì)算公式為

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

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

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

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

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

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

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

由μ、ω可得IB 的湍動(dòng)能k［19］值為

2.3 數(shù)值模擬流程

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

1）通過式（29）計(jì)算出IP 處 的速度U和壓強(qiáng)p。

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

3）通過OpenFOAM 求解流場(chǎng)控制方程離散矩陣。

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

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

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

3 數(shù)值算例

3.1 平板繞流

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

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

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

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

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

圖4 展示了Re=7×105時(shí)平板表面Cf分布，它與冪律公式的計(jì)算結(jié)果基本一致，體現(xiàn)了本文方法對(duì)于一定雷諾數(shù)范圍的湍流平板問題的模擬較為準(zhǔn)確。

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

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

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

3.2 NACA0012 翼型繞流

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

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

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

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

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

3.3 Buice 二維擴(kuò)散管

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

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

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

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

3.4 二維圓柱繞流

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

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

圖12、13 展示了在Re=63 100，126 000 下圓柱繞流的阻力系數(shù)CD和升力系數(shù)CL隨時(shí)間變化的曲線，表1 為平均阻力系數(shù)、平均升力系數(shù)、升力系數(shù)均方根CL，rms、斯特勞哈爾數(shù)St實(shí)驗(yàn)值和數(shù)值模擬結(jié)果［28-31］對(duì)比?？梢钥闯?，本文的數(shù)值模擬結(jié)果與文獻(xiàn)結(jié)果基本一致。但是，與Yeon 等［29］使用LES 得到的結(jié)果對(duì)比，CL、St稍大，這可能是由于LES 和RANS 對(duì)湍流?；呗圆煌?，以及本文只對(duì)近壁面附近的網(wǎng)格加密而沒有對(duì)尾流區(qū)加密導(dǎo)致誤差產(chǎn)生。CD、CL、St的計(jì)算公式為

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

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

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

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

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

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

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

3.5 二維方柱繞流

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

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

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

圖18 Re=22 000 時(shí)方柱繞流計(jì)算結(jié)果Fig.18 Calculation results of flow around square cylin‐der at Re=22 000

4 結(jié)論

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

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

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