趙體豪, 趙 欣

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

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

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

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

1 算法推導(dǎo)

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

(1)

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

(2)

式中，上標(biāo)n表示物理量的當(dāng)前時(shí)刻，n+1表示下一時(shí)刻.

對動(dòng)量方程進(jìn)行離散可以得到速度項(xiàng)的離散方程[25]如下：

(3)

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

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

(4)

(5)

對于剛體內(nèi)部，其壓力方程具有如下形式：

將該式改寫后可得

(6)

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

(7)

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

(8)

將式(7)和式(8)進(jìn)行離散化處理，并且結(jié)合式(4)整理可得適用于整個(gè)計(jì)算域中壓強(qiáng)p′的方程為

(9)

至此，我們在弱可壓縮粘性流動(dòng)與BDIM的基礎(chǔ)上，得到了可用于求解包含運(yùn)動(dòng)固體的可壓縮粘性流動(dòng)的完整算法.

2 算法驗(yàn)證

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

2.1 靜止圓柱繞流

在本節(jié)，我們首先計(jì)算了靜止圓柱繞流問題，驗(yàn)證了算法在計(jì)算包含靜止固體的可壓縮粘性流動(dòng)問題時(shí)的有效性及準(zhǔn)確性.

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

表1Re=100時(shí)靜止圓柱繞流的阻力、升力系數(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文獻(xiàn)[30]1.5010.3490.172 文獻(xiàn) [30]*1.4530.3390.169文獻(xiàn)[31]1.3500.3390.165文獻(xiàn)[32]1.3300.332—文獻(xiàn)[33]1.3420.2970.166

注：*為增大計(jì)算域后所得結(jié)果.

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

2.2 振動(dòng)圓柱繞流

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

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

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

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

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

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

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

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

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

表2Re=185時(shí)振動(dòng)圓柱繞流的時(shí)均阻力系數(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文獻(xiàn)[30,37]1.2820.0620.223文獻(xiàn)[36]1.25—0.18

表3 圓柱振動(dòng)頻率為0.156時(shí)，不同振幅對應(yīng)的時(shí)均阻力系數(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 圓柱振動(dòng)幅度為0.2d時(shí)，不同振動(dòng)頻率對應(yīng)的時(shí)均阻力系數(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)點(diǎn)呈中心對稱. 此外，在圖5中也給出了流動(dòng)穩(wěn)定時(shí)的渦量圖.

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

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

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

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

2.3 橫向振動(dòng)圓柱

上述2個(gè)算例對算法的有效性及準(zhǔn)確性做了驗(yàn)證，本節(jié)通過驗(yàn)證橫向振動(dòng)圓柱，對算法在聲學(xué)層面的應(yīng)用進(jìn)行驗(yàn)證.

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

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

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

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

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

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

3 結(jié) 論

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

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