趙體豪, 趙 欣

(1.北京理工大學 機械與車輛學院,北京 100081；2.北京理工大學 宇航學院,北京 100081)

計算成本及計算能力的限制使得我們在做計算時常常要將計算對象假定為不可壓縮無粘流體. 這一物理假設雖然可以在一定程度上減小計算成本，但同時也會在結果中引入一定程度的誤差. 為解決上述問題，本文提出了一種適用于求解弱可壓縮粘性流動問題的數(shù)值解法. 該方法的數(shù)學推導基于邊界數(shù)據(jù)浸入法(Boundary Data Immersion Method, BDIM)，并且該方法充分考慮了運動固體對流場計算結果的影響.

對于計算域內含有動邊界的問題，浸入邊界法(Immersed Boundary Method, IBM)是一種準確有效的求解方法. 該法由Peskin[1]于1972年首次提出，用于對人類心臟中血液流動的模擬. 其主要思想是將流場的復雜邊界等效為一個體積力并施加于N-S動量方程中. 此后的幾十年中，該法得到了不斷發(fā)展且被廣泛使用. 為改善體積不守恒問題，Peskin[2]引入新的投影算子替代原有方程中的哈密頓算子；自適應加密算法與IBM的結合[3-5]使得計算效率得到了提高;針對流體與運動剛體之間相互作用的研究，這些年來逐漸發(fā)展出了離散直接力法[6-8]、虛擬邊界法[9-10]、反饋力法[11-13]等不同方法. 此外，如何提高邊界浸入法的計算精度也引起了不少研究學者的關注. Yusof[14]將B樣條(B-spline)與浸入邊界法結合以獲得高階精度，且消除了時間步長的限制；Tseng[10]通過在邊界上靠近固體一側的區(qū)域進行虛擬網(wǎng)格重建的方法來提高計算精度；Yang[15]通過對離散delta函數(shù)的光滑處理，有效抑制了體積力中的非物理振蕩；宮兆新等[16]通過虛擬解法驗證了正則化δ函數(shù)的改變對于計算精度的提高沒有任何幫助；李秋實等[17]通過改變力源的構造以及邊界速度插值的方式獲得二階精度結果；及春寧等[18-20]將力源項的求解內嵌到壓力泊松方程的迭代過程，得到了一種高精度的嵌入式迭代浸入邊界法.

然而與原始浸入邊界法相同，上述研究對象大多為不可壓縮無粘流體，并未成功地將邊界浸入法推廣到可壓縮流動中來. 對于不可壓縮流動來說，使用邊界浸入法只需要給出速度邊界條件(即無滑移速度邊界)即可；但是對于可壓縮流動來說，因流體的一些變量，如密度、溫度等都會在運動中發(fā)生較大的變化，所以需要除了速度邊界條件以外的其他邊界條件才可以順利完成流場的求解. 為了將浸入邊界法應用到可壓縮流動中來，研究人員進行了大量工作[21-24]. Palma等[21-22]將IBM與求解預置可壓縮N-S方程的解法相結合用于求解流場中復雜幾何結構的問題；Qiu等[23]依次使用無滑移速度邊界條件對流體的密度、速度、溫度以及壓力依次進行修正，推導出了一種邊界條件浸入法(Boundary condition-enforced)；Schlanderer[25]提出一種適用于求解可壓縮流動的浸入邊界法——邊界數(shù)據(jù)浸入法(Boundary Data Immersion Method). 該方法使用廣義積分核公式將計算域內不同子域的控制方程進行重組，得到具有統(tǒng)一形式的方程，并且使得重組后的方程不僅在各子域內有效，且在邊界處也可以實現(xiàn)光滑過渡. 此外，現(xiàn)在對于弱可壓縮粘性流動問題的求解方法主要是predictor-corrector方法[26-27]或者人工可壓縮的方法[28]. 在這里我們將借鑒Caltagirone[29]對流體弱可壓縮粘性的處理方式，并結合Schlanderer[25]的邊界數(shù)據(jù)浸入法推導出一種可以求解水下弱可壓縮粘性流動問題的簡單算法.

1 算法推導

對于單相弱可壓縮粘性的流體，從連續(xù)性方程與流體狀態(tài)方程出發(fā)，經(jīng)過推導可以得到如下關于壓力項的方程[29]：

(1)

式中:γ為比熱容比系數(shù)，這里取6.1；B為水的剛度系數(shù)，其具有與壓力相同的量綱， 取值345 MPa. 將上式代入式(1)，并對其進行離散化處理，有

(2)

式中，上標n表示物理量的當前時刻，n+1表示下一時刻.

對動量方程進行離散可以得到速度項的離散方程[25]如下：

(3)

式中:ρn為n時刻的流體密度；μ為流體的動力黏度系數(shù)，且上述物理量單位均采用國際單位制.

在得到流體的壓力與速度的離散方程之后，可以對流場進行求解. 但在實際的流場計算中，常常會遇到整個計算域中存在多個不同物理性質的子域，因此需要得到求解整個流場的計算方法. 為解決上述問題，Schlanderer等[25]提出了基于廣義積分核公式的邊界數(shù)據(jù)浸入法. 這種方法通過積分核公式，將不同子域物理變量的控制方程和界面條件結合起來，實現(xiàn)了對整個計算域的模擬. 使用這種方法得到的具有統(tǒng)一形式的控制方程不僅保留了原系統(tǒng)中的物理變量特性，還做到了不同子域間的光滑過渡，使得我們可以對包含運動物體的可壓縮粘性流動進行計算. 本文在這里直接使用可將各子域控制方程連結起來的控制方程[25]：

(4)

(5)

對于剛體內部，其壓力方程具有如下形式：

將該式改寫后可得

(6)

為滿足剛體與流體邊界處上壓力連續(xù)變化的條件，剛體邊界上壓力方程修正后有如下形式：

(7)

其中式(7)中右邊第二項恒為0. 流體的壓力方程有以下形式:

(8)

將式(7)和式(8)進行離散化處理，并且結合式(4)整理可得適用于整個計算域中壓強p′的方程為

(9)

至此，我們在弱可壓縮粘性流動與BDIM的基礎上，得到了可用于求解包含運動固體的可壓縮粘性流動的完整算法.

2 算法驗證

為對算法的有效性及準確性進行充分說明，我們在本章計算了3個二維圓柱的經(jīng)典算例，并將計算結果與前人在相同工況下的計算結果進行了對比.

2.1 靜止圓柱繞流

在本節(jié)，我們首先計算了靜止圓柱繞流問題，驗證了算法在計算包含靜止固體的可壓縮粘性流動問題時的有效性及準確性.

基于入流速度u與圓柱直徑d，設置流動的雷諾數(shù)Re=ud/υ=100；設置計算域的大小為30d×8d，并在整個計算域上均布尺寸大小為1/30d的計算網(wǎng)格；時間步長為1×10-4. 計算采用入口速度條件，出口壓力條件，其余邊界采用無滑移邊界條件.

表1Re=100時靜止圓柱繞流的阻力、升力系數(shù)與斯特勞哈爾數(shù)

Tab.1 Drag and lift coefficients and Strouhal number for flow around a stationary cylinder atRe=100

方法C-DCL'St本文方法1.3600.3340.183文獻[30]1.5010.3490.172 文獻 [30]*1.4530.3390.169文獻[31]1.3500.3390.165文獻[32]1.3300.332—文獻[33]1.3420.2970.166

注：*為增大計算域后所得結果.

圖1中給出了當流動穩(wěn)定時，圓柱表面的時均壓力系數(shù)CP=(P-P)/u2的分布情況. 這里，提取數(shù)據(jù)的壓力測點均為圓柱外側最近的網(wǎng)格節(jié)點. 圖中θ=0°與θ=180°分別對應圓柱上距離來流最近與最遠的兩點. 在圖1中，我們將分別使用網(wǎng)格尺寸大小為1/20d、1/30d及1/40d得到的計算結果與實驗結果[34-35]進行比對. 近些年技術水平的提高使得學者們將關注點放在高雷諾數(shù)的工況上，較少有人對低雷諾數(shù)的工況進行實驗，故本文引用的參考文獻時間較早. 該對比結果表明：使用本算法得到的計算結果與實驗數(shù)據(jù)[34-35]整體趨勢吻合較好. 將兩個實驗結果對比可以發(fā)現(xiàn)，較大的雷諾數(shù)會使得壓力系數(shù)CP在50-180°小于雷諾數(shù)較小時的壓力系數(shù). 該現(xiàn)象可以解釋為什么使用數(shù)值計算得到的結果在θ角度較大時會略小于實驗數(shù)據(jù). 此外，對比發(fā)現(xiàn)使用3種不同網(wǎng)格尺寸大小得到的計算結果相差不大，由此可得出計算結果與網(wǎng)格尺寸大小無關性的結論. 為使結果更加直觀，在圖2中提供了流動穩(wěn)定后的渦量圖. 至此，驗證了算法在計算包含靜止物體的可壓縮粘性流動問題時的有效性及準確性.

2.2 振動圓柱繞流

上節(jié)使用算法計算了包含靜止固體的可壓縮粘性流動的算例，本節(jié)更進一步，驗證當物體發(fā)生運動時算法的有效性及準確性. 本節(jié)考慮另一個二維算例——橫向振動圓柱繞流的算例.

圖1Re=100時，流動穩(wěn)定后靜止圓柱時均表面壓力系數(shù)CP分布圖

Fig.1 Distribution of time-averaged surface pressure coefficient on stationary cylinder when flow is stable atRe=100

圖2 Re=100時靜止圓柱擾流穩(wěn)定時渦量圖

Fig.2 Vortex diagram for flow past a stationary cylinder when the flow is stable atRe=100

基于入流速度u與圓柱直徑d，設置流動的雷諾數(shù)Re=ud/υ=185. 設置計算域的大小為30d×10d，并在整個計算域上均布尺寸為1/30d的計算網(wǎng)格；時間步長為2×10-4. 計算采用入口速度條件，出口壓力條件，其余邊界采用無滑移邊界條件. 圓柱在垂直于來流方向上以正弦規(guī)律進行振蕩，其振蕩幅度A=0.2d，振動頻率為fo=0.156. 圓柱位移隨時間的變化可用如下數(shù)學函數(shù)關系進行表達：

y(t)=Asin(2πfot).

表2對比了采用不同數(shù)值方法計算得到時均阻力系數(shù)以及阻力系數(shù)和升力系數(shù)的均方根. 通過對比，可以發(fā)現(xiàn)：使用我們的算法計算得到的3個數(shù)值在量級上與其他算法所得結果相同，具有較好的一致性；具體而言，計算得到的時均阻力系數(shù)值比其他算法得到的數(shù)值小1.1%～7.8%，阻力系數(shù)的均方根同文獻[30,37]中結果幾乎相同，而升力系數(shù)的均方根則略大于文獻[30,36-37]中的計算結果. 產(chǎn)生這一差異化現(xiàn)象的原因是參考文獻中采用數(shù)值計算方法均為不可壓縮流體，其計算域的壓力出口邊界條件與本文采用的弱可壓縮流體的壓力邊界條件不同.

此外，本文引入不同的圓柱振動特性進行計算，探究其變化對計算結果的影響. 文中采用改變單一變量的方法，計算了振動頻率為0.156時，圓柱振動幅度A分別為0.2d、0.3d、0.4d、0.5d的4種情況以及振動幅度A為0.2d時，圓柱振動頻率分別為0.156、0.3、0.45的3種情況. 對比結果見表3及表4. 從表3可以看出，振動幅度的增大，使得3種系數(shù)均呈現(xiàn)不斷增大的趨勢；而從表4可以看出，隨著圓柱振動頻率的增大，3個系數(shù)均有不同程度的增幅，其中，升力系數(shù)均方根的增幅遠大于時均阻力系數(shù)與阻力系數(shù)均方根，達到了十倍以上.

表2Re=185時振動圓柱繞流的時均阻力系數(shù)、阻力及升力系數(shù)的均方根

Tab.2 Time-averaged drag coefficient and root mean square of drag and lift coefficient for flow around an oscillating cylinder atRe=185

方法C-D(CD)rms(CL)rms本文方法1.1820.0620.231文獻[30,37]1.2820.0620.223文獻[36]1.25—0.18

表3 圓柱振動頻率為0.156時，不同振幅對應的時均阻力系數(shù)、阻力及升力系數(shù)的均方根

Tab.3 Time-averaged drag coefficient and root mean square of drag and lift coefficient associated with different oscillating amplitudes when the oscillating frequency is 0.156

圓柱振幅/dC-D(CD)rms(CL)rms0.21.1820.0620.2310.31.2650.1230.2750.41.3470.1870.2880.51.4260.2570.306

表4 圓柱振動幅度為0.2d時，不同振動頻率對應的時均阻力系數(shù)、阻力及升力系數(shù)的均方根

Tab.4 Time-averageddrag coefficient and root mean square of drag and lift coefficient associated with different oscillating frequencies when the oscillating amplitude is 0.2d

圓柱振頻C-D(CD)rms(CL)rms0.1561.1820.0620.2310.301.3790.1451.3090.451.5060.2432.494

圖3和圖4給出了阻力系數(shù)及升力系數(shù)隨圓柱位移變化的圖像. 從圖中可以看出，無論是阻力還是升力，其數(shù)值大小均為對稱分布，其中，前者是以靜止位置(yc=0)軸對稱，而后者則以(0,0)點呈中心對稱. 此外，在圖5中也給出了流動穩(wěn)定時的渦量圖.

圖3 Re=185時振動圓柱擾流的阻力系數(shù)隨位移變化圖

Fig.3 Variation of drag coefficient with displacement for flow past an oscillating cylinder atRe=185

圖4 Re=185時振動圓柱擾流的升力系數(shù)隨位移變化圖

Fig.4 Variation of lift coefficient with displacement for flow past an oscillating cylinder atRe=185

2.3 橫向振動圓柱

上述2個算例對算法的有效性及準確性做了驗證，本節(jié)通過驗證橫向振動圓柱，對算法在聲學層面的應用進行驗證.

計算域的大小為30d×30d，整個計算域的網(wǎng)格尺寸為0.05d，時間步長為0.01. 圓柱在橫向上的振動規(guī)律同2.2節(jié)，即圓柱位移與時間的函數(shù)關系式為

y(t)=Asin(2πfot).

圖5 Re=185時振蕩圓柱流動穩(wěn)定時渦量圖

Fig.5 Vortex diagram for flow past an oscillating cylinder when flow is stable atRe=185

圖6給出了圓柱在振動穩(wěn)定時，計算域內壓力云圖變化的整個周期. 通過圖6可以看出，圓柱在發(fā)生振動時所產(chǎn)生的壓力脈動在聲學上等同于一個偶極子，這與理論分析相一致. 此外，從圖7中可以看出，對方向譜而言，計算結果同偶極子振動的理論解有很好的一致性. 從而驗證了算法在計算包含運動物體的可壓縮粘性流動時聲學層面的有效性.

(a)振動周期階段一 (b)振動周期階段二圖6 圓柱振動壓力云圖Fig.6 Pressure contour of oscillating cylinder

3 結 論

相比于傳統(tǒng)的浸入邊界法，本文所采用的邊界數(shù)據(jù)浸入法不僅可對包含運動物體的流域進行求解，更可與弱可壓縮粘性流結合起來以考慮流體的可壓縮性及粘性對計算結果的影響. 本文從可壓縮粘性流的連續(xù)性方程及水的狀態(tài)方程出發(fā)，使用邊界數(shù)據(jù)浸入法將運動物體的固體子域與流體子域的變量方程進行重組，嚴格推導了適用于整個流場的壓力與速度方程. 該方程不僅有效保留了在各子域內的所有物理性質，而且成功地在各子域之間交界面上實現(xiàn)了光滑過渡.

此外，通過3個二維經(jīng)典算例，在力學與聲學層面上驗證了算法在計算包含靜止物體及運動物體的流場時的有效性. 通過將計算所得結果與其他算法所得結果進行對比，對算法的準確性進行了驗證說明. 相較于目前已有的可壓縮粘性流動問題的求解方法，本算法推導簡單，無需過多假設，且計算結果與實驗數(shù)據(jù)具有良好的吻合程度.