趙西增，付英男，張大可

(1. 浙江大學(xué) 海洋學(xué)院，浙江 杭州 310058; 2. 國家海洋局第二海洋研究所 衛(wèi)星海洋環(huán)境動力學(xué)國家重點實驗室,浙江 杭州 310012)

?

柱體繞流的CIP方法模擬

趙西增1，2，付英男1，張大可1

(1. 浙江大學(xué) 海洋學(xué)院，浙江 杭州 310058; 2. 國家海洋局第二海洋研究所 衛(wèi)星海洋環(huán)境動力學(xué)國家重點實驗室,浙江 杭州 310012)

摘要：柱體繞流和流致振動現(xiàn)象是一個復(fù)雜的工程問題，利用自主研發(fā)的CIP-ZJU模型，對低雷諾數(shù)(Re<300)圓柱和方柱繞流問題開展了數(shù)值模擬。模型在直角坐標(biāo)系統(tǒng)下建立，采用緊致插值曲線CIP方法作為流場的基本求解器離散了Navier-Stokes方程，基于多相流的理論實現(xiàn)流-固耦合同步求解，利用浸入邊界方法處理固體邊界。模擬結(jié)果與文獻結(jié)果進行比較，二者吻合情況較好。通過引入壓力阻力項和摩擦阻力項，分析了Strouhal數(shù)、阻力系數(shù)、升力系數(shù)等參數(shù)隨雷諾數(shù)的變化情況。結(jié)果表明，圓柱繞流與方柱繞流在水動力參數(shù)的變化規(guī)律上存在多處差異。

關(guān)鍵詞：柱體繞流；CIP方法；Navier-Stokes方程；浸入邊界法；Strouhal數(shù)

鈍體繞流，特別是柱體繞流問題在工程實際中經(jīng)常遇到，如風(fēng)吹過高層建筑物、海水流過石油鉆井平臺、海流流經(jīng)海洋立管等。其共同的特點是，在一定雷諾數(shù)范圍內(nèi)，粘性流體流過柱體結(jié)構(gòu)物時，會發(fā)生邊界層的分離，進而對結(jié)構(gòu)物產(chǎn)生周期性作用。因此，準(zhǔn)確計算出流體對柱體結(jié)構(gòu)的作用荷載及其在荷載作用下的動力響應(yīng)，具有重要的科學(xué)價值和工程意義。

許多學(xué)者采用試驗方法和數(shù)值模擬方法開展了鈍體繞流問題的研究。Tritton[1]在實驗室內(nèi)通過觀察石英試管的彎曲度得到Re在0.5～100范圍內(nèi)圓柱繞流的阻力，與理論值吻合良好。Kawaguti等[2]對Re在1～100圓柱繞流問題進行了數(shù)值模擬，并考慮了壓力和摩擦力對阻力的貢獻，結(jié)果表明壓力阻力項在阻力中占有主導(dǎo)地位。Rajani等[3]利用二階精度的中心差分格式對對流方程進行離散，對圓柱繞流進行了二維和三維數(shù)值模擬，將兩者的渦脫頻率和升力系數(shù)進行對比，發(fā)現(xiàn)大于某一臨界雷諾數(shù)時，二維結(jié)果偏大，三維效應(yīng)不可忽略。方柱繞流也是鈍體繞流的一種典型形式。Okajima[4]研究了不同類型的矩形柱體的渦旋脫落頻率，并指出Strouhal數(shù)是雷諾數(shù)的函數(shù)。Sohankar等[5]引入二階精度的Van Leer迎風(fēng)格式對對流項進行求解，指出Re>50時會有穩(wěn)定的卡門渦街產(chǎn)生。Breuer等[6]將格子玻爾茲曼方法(lattice Boltzmann method，LBM)和有限體積法 (finite volume method，F(xiàn)VM)分別應(yīng)用到二維方柱繞流數(shù)值模擬中，兩種方法的計算結(jié)果顯示了較好的一致性。Sen等[7]利用有限元方法(finite element method，F(xiàn)EM)對Re<150時的方柱繞流進行二維數(shù)值模擬，發(fā)現(xiàn)最初的分離點是在基點處而不是在后角點處。潘小強等[8]利用Fluent分別模擬了二維圓柱和方柱繞流的流場情況，結(jié)果表明同樣雷諾數(shù)的圓柱和方柱繞流，圓柱繞流更易得到穩(wěn)態(tài)非定常解。周云龍等[9]采用FVM模擬了二維圓柱和方柱的流動分離狀態(tài)，發(fā)現(xiàn)兩者最大的不同點在于分離點的不同，圓柱繞流沒有固定分離點而方柱繞流固定分離點在其前后角點處。不同柱體繞流過程中的水動力特性的差異，仍未明朗，這將是本文工作的重點。

本文采用自主研發(fā)的CIP方法模型[10-11]，通過浸入邊界方法處理流-固耦合問題，對不同雷諾數(shù)(Re<300)的圓柱繞流和方柱繞流問題進行數(shù)值模擬，并通過與文獻結(jié)果的比較驗證模型的有效性。為了得到穩(wěn)定的流場信息，本文采用了與傳統(tǒng)差分格式不同的高階格式CIP方法對N-S方程進行了離散，該方法通過對網(wǎng)格內(nèi)結(jié)點函數(shù)值及其導(dǎo)數(shù)值聯(lián)立，實現(xiàn)了三階精度的差分求解對流方程，具有格式穩(wěn)定和低耗散性的特點；同時，在直角網(wǎng)格系統(tǒng)內(nèi)，通過引入浸入邊界方法，有效處理了流-固耦合問題，可有效避免動態(tài)網(wǎng)格系統(tǒng)中的大量信息交換，保證了模型的計算效率，為進一步模擬運動物體做好準(zhǔn)備。本文將從多個角度分析圓柱繞流與方柱繞流的區(qū)別。

1數(shù)學(xué)模型的建立

1.1控制方程

本流場模型以二維N-S方程為控制方程，其矢量形式為

(1)

(2)

為了把CIP方法應(yīng)用于對流方程，對方程(2)取空間導(dǎo)數(shù)，可得到

(3)

模型在直角笛卡爾系統(tǒng)下建立，采用多相流理論處理固-液的相互作用，滿足下面的控制方程：

(4)

式中：m=1, 2,φ1為液體相，φ2為固體相，在一個網(wǎng)格內(nèi)滿足φ1+φ2=1。固體相的處理基于浸入邊界方法得到，將在下面的部分給出解釋。網(wǎng)格內(nèi)的流體特性可用下式來表示：

(5)

式中：λ為密度ρ或者粘性系數(shù)μ。

1.2數(shù)值方法

1.2.1CIP方法

CIP方法是Takewaki等[12]于20世紀80年代中期提出并發(fā)展起來，用于求解雙曲型偏微分方程的一種有效的數(shù)值計算方法。其基本原理是基于空間網(wǎng)格點的變量值及其空間導(dǎo)數(shù)值，利用三次多項式進行插值近似，反演出網(wǎng)格單元內(nèi)部變量的真實信息，得到三階精度的顯式格式。為了更好地解釋CIP方法，考慮簡單的常系數(shù)一維對流方程：

(6)

式中：u為大于零的常數(shù)。此類偏微分方程可用不同差分方法求解。下面將給出CIP方法的工作原理。圖1(a)～(c)給出了一階迎風(fēng)差分方法的工作原理，對于一階差分格式來說，它利用直線的方式聯(lián)立相鄰兩個節(jié)點的信息來工作，而忽略了網(wǎng)格內(nèi)部的信息，導(dǎo)致較大的數(shù)值耗散。為了真實再現(xiàn)網(wǎng)格內(nèi)部的信息，有必要求助高階差分方法，而對于常規(guī)高階差分的建立需要更多的網(wǎng)格點的信息。CIP采用一種獨特的方式，在一個網(wǎng)格內(nèi)實現(xiàn)了高階差分格式，通過利用空間網(wǎng)格點的變量值及其空間導(dǎo)數(shù)值，來描述該網(wǎng)格內(nèi)的信息，可真實再現(xiàn)網(wǎng)格內(nèi)的信息。

通過對方程(6)求關(guān)于x的偏導(dǎo)，得到如下的空間導(dǎo)數(shù)方程：

(7)

(8)

在n+1時刻的單元格剖面函數(shù)fn+1可以通過將n時刻的剖面函數(shù)fn平移-u△t得到，函數(shù)f和g的時間演變可以通過下面的拉格朗日變換得到

(9)

圖1　CIP方法的基本原理Fig. 1　Principle of CIP method

圖2　半拉格朗日方法的CIP格式Fig. 2　CIP method of semi-Lagrangian format

(10)

CIP方法采用一種獨特的方式，在一個網(wǎng)格內(nèi)實現(xiàn)了高階差分格式，使得本文的數(shù)值模型可以應(yīng)用于復(fù)雜流動問題的模擬，這也是本文所采用的差分方法與其他方法的不同之處。

1.2.2時間積分

采用分步算法對動量方程進行時間積分。首先忽略擴散項和壓力項，只考慮對流項的中間速度；其次求解擴散項；然后求解壓力方程，計算下一時間步的壓力；最后考慮壓力梯度項，計算速度的最后值。設(shè)△t為時間步長，在t=n△t到t=(n+1)△t時刻的計算時間內(nèi)，具體的時間積分過程如下：

1) 對流項(I)：

(11)

通過CIP方法[13]可得到方程的解：

(12)

式中：‘*’為對流項計算結(jié)束后的中間時間標(biāo)志。

2)非對流項(I):

擴散項的計算：

(13)

擴散項的時間離散采用顯示格式，中間速度可表示為下面的形式：

(14)

其中,方程(14)的空間離散采用中心差分格式。

3)非對流項(II):

壓力和速度的匹配:

(15)

通過對方程 (15)第一行公式取散度，并引入連續(xù)方程，可得到如下形式的泊松方程：

(16)

泊松方程的求解通過SOR迭代得到。

考慮動量方程的壓力梯度項，計算速度的最終值，計算式為

(17)

根據(jù)式(4)、(5)和(18)更新界面和網(wǎng)格內(nèi)流體信息，然后返回到步驟1，這樣就完成了整個計算過程，一直到設(shè)定的步驟結(jié)束。

1.2.3固體邊界處理

采用浸入邊界法(immersedboundarymethod)[14]處理固體對流場的影響:

(18)

式中：ufn+1為在歐拉坐標(biāo)體系中通過式(17)計算得到的“局部”流場速度；Ubn+1為通過拉格朗日方法得到的固體速度[15]，本文中柱體靜止，固體的速度設(shè)定為零；un+1為二者耦合作用的“整體”速度，這樣就得到了流-固耦合流場的求解。上述方法也可稱為固體體積流-固耦合處理方法?；诮脒吔绶椒ǖ牧?固耦合技術(shù)，采用固定網(wǎng)格，在計算過程中無需更新網(wǎng)格，可保證模型的高效率，使得模型在處理大位移的流固耦合計算中具有強大優(yōu)越性。本文采用的方法具有較易實施的優(yōu)勢，方便模型向三維拓展。

1.3參數(shù)定義

升力Fl為

(19)

式中：Re為雷諾數(shù),Cd為阻力系數(shù),Fd為阻力,Cl為升力系數(shù),St為Strouhal數(shù),上角標(biāo)p和v分別表示由于壓力和摩擦力所產(chǎn)生的阻力或升力，Sij=(?ui/

?xj+?uj/?xi)/2。

1.4數(shù)值模型

計算區(qū)域及網(wǎng)格劃分如圖3和圖4所示。圓柱與方柱的特征長度均為D=1m，上游邊界距離柱體中心7.5D，下游邊界距離柱體中心22.5D，兩側(cè)邊界距離柱體中心均為7.5D。為了捕捉柱體周圍的復(fù)雜流動，柱體近壁面采用細網(wǎng)格。

圖3　圓柱、方柱繞流計算區(qū)域示意圖Fig. 3　Sketch of the computational domain

圖4　圓柱、方柱繞流網(wǎng)格示意圖Fig. 4　Sketch of the computational mesh

1.5邊界條件和初始條件設(shè)置

計算流場的邊界條件以及初始條件設(shè)置如下：初始條件：初始設(shè)定一個從左向右的速度場，u=U0，v=0；其中，u、v分別為x、y方向的速度；初始壓力場都設(shè)定為零。邊界條件：來流邊界保持為均勻來流速度U0；出口邊界條件，開邊界；柱體表面邊界條件：無滑移邊界，即un=0；兩側(cè)壁面：自由滑移邊界。

2模型驗證及數(shù)據(jù)分析

2.1網(wǎng)格收斂性分析

為驗證模型的有效性及網(wǎng)格收斂性，采用三套網(wǎng)格對Re=100時圓柱繞流進行模擬。三種網(wǎng)格從網(wǎng)格1到網(wǎng)格3，最小網(wǎng)格分別為0.04、0.02和0.01m，模擬結(jié)果如表1所示。由于網(wǎng)格1的網(wǎng)格尺寸較大，Cd值與文獻結(jié)果相差較為明顯，而網(wǎng)格2與網(wǎng)格3誤差較小。考慮到計算效率和精度，選擇網(wǎng)格2對下文其他工況進行模擬。

2.2圓柱繞流參數(shù)分析

下面對Re=1～300的圓柱繞流進行模擬，結(jié)果如表2所示。從表中可看出，Re<50時，沒有漩渦脫落，St數(shù)為零。

表1　網(wǎng)格數(shù)據(jù)及主要參數(shù)對比

表2不同Re圓柱繞流主要參數(shù)結(jié)果

Table 2Main parameters of flow around a cylinder at differentRe

ReCdSt113.4—54.53—103.11—202.18—301.81—401.59—501.470.131601.440.139801.370.1571001.330.1661201.320.1751401.3250.1821601.3150.1852001.3250.1922501.330.1993001.350.206

圖5給出了圓柱繞流主要水動力參數(shù)隨Re的變化情況。在Re<40時，沒有漩渦脫落為層流狀態(tài)，如圖5(a)所示在開始段Cd驟減，主要是因為4

圖5　主要參數(shù)隨Re的變化Fig. 5　The variation of main parameters vs Re

圖6給出了50100時由于考慮到實際情況的三維效應(yīng)，數(shù)模的二維結(jié)果會高估St數(shù)。主波長λ為卡門渦街一側(cè)前后2個渦之間的距離，如圖7所示，λ與St數(shù)的關(guān)系可近似表示為λ/D=1/St[22]。由圖可知隨著Re的增大，λ呈單調(diào)下降趨勢，卡門渦街漩渦之間的距離越來越小。

2.3方柱繞流參數(shù)分析

下面將對Re<300情況下的方柱繞流進行模擬，結(jié)果如表3所示。與圓柱繞流結(jié)果類似，本文模擬出的結(jié)果在Re<50時顯示沒有漩渦脫落，依舊為層流，St數(shù)值為零；而在Re=50時，不穩(wěn)定狀態(tài)首次顯現(xiàn)，伴隨有漩渦脫落，流動轉(zhuǎn)化為非層流，說明在Re=40～50，漩渦開始脫落，與Sohankar等[23]通過數(shù)值模擬測得的Re=52較為接近。

圖6　Strouhal數(shù)和主波長λ隨Re的變化曲線Fig. 6　The variation of Strouhal number and primary wavelength vs Re

圖7　主波長λ示意圖Fig. 7　Sketch of primary wavelength

Table 3Main parameters of flow around a cylinder at differentRe

ReCdSt114.1—54.85—103.28—202.33—301.95—401.74—501.6250.116601.5650.126801.510.1361001.490.1431201.480.1471401.4850.151601.490.152001.510.1432501.630.1363001.720.136

圖8　主要參數(shù)隨Re的變化Fig. 8　The variation of main parameters vs Re

ReCdCd(p)Cd(v)1601.491.5-0.012001.511.54-0.032501.631.67-0.043001.721.75-0.03

與圓柱繞流類似，在Re<40時，是層流狀態(tài)。圖8(a)給出了Cd、Cd(p)、Cd(v)的變化，可發(fā)現(xiàn)整體趨勢相同，最后近似發(fā)展為3條平行的直線，與Sen等[7]在2100時，與畢繼紅等[25]的結(jié)果相差較小，與其他文獻相差較大，其中的影響因素較多，像離散方法不同、計算域的大小、邊界條件設(shè)置的不同均可導(dǎo)致計算結(jié)果的誤差。在Re≈140時，阻力系數(shù)達到極小值，此時分離點由后角點處轉(zhuǎn)移到前角點處(本文將此時定義為“轉(zhuǎn)折點”，方便后文敘述)，使得方柱側(cè)面受力發(fā)生變化，導(dǎo)致Cd(v)的合力方向發(fā)生變化，在Re<140時，Cd(v)數(shù)值為正；而如表4所示，在Re>140時，Cd(v)數(shù)值為負，方向發(fā)生了變化，這也是方柱繞流與圓柱繞流最大的區(qū)別之一。由圖8(c)可看出方柱繞流中，壓力阻力項仍然占有主導(dǎo)地位，且隨著雷諾數(shù)的增大，摩擦阻力項相比于壓力阻力項基本可以忽略。從圖8(d)中可以看出在Re≈140時，升力系數(shù)的增長率逐漸減小，這也是分離點變化導(dǎo)致的結(jié)果，過了轉(zhuǎn)折點之后，升力系數(shù)的增長率變大，結(jié)果與現(xiàn)有文獻吻合良好。從圖8(e)可看出St數(shù)的變化趨勢與現(xiàn)有文獻基本一致，在轉(zhuǎn)折點處達到最大值，而數(shù)值方面的誤差可能由多種原因?qū)е拢珺reuer等[6]指出阻塞比的增加會導(dǎo)致St數(shù)的增大。主波長λ的變化與St數(shù)相反，先減小再增大，在轉(zhuǎn)折點處的主波長達到最小。

2.4圓柱與方柱繞流參數(shù)對比

為了清晰說明圓柱繞流與方柱繞流的區(qū)別，將結(jié)果中的重要參數(shù)進行對比。

圖9給出了圓柱、方柱繞流主要水動力參數(shù)隨Re的變化曲線。由圖9(a)可以看出，Cd mean隨Re的變化趨勢大體一致，在Re<10時，Cd mean隨Re的變大急劇減小，然后趨于平穩(wěn)，且方柱的阻力系數(shù)大于圓柱，這是兩者形狀不同所導(dǎo)致的。兩者在Re≈140～150達到極小值，此時方柱的分離點由后角點轉(zhuǎn)移到前角點，方柱摩擦阻力的方向發(fā)生變化，而圓柱沒有此現(xiàn)象發(fā)生，其分離點只圍繞一點往復(fù)運動，且摩擦阻力保持為正值，方向不變。圖9(b)給出了圓柱和方柱Cl max隨Re的變化曲線。在Re<150時，兩者的大小基本一致，經(jīng)過轉(zhuǎn)折點之后，方柱的升力系數(shù)顯著增大，大體呈線性增長，分離點在前角點時，方柱兩側(cè)面將受到附面渦的作用，使得升力增大，而圓柱不存在此現(xiàn)象，升力系數(shù)變化相對平穩(wěn)。圖9(c)給出了圓柱和方柱St隨Re的變化曲線，圓柱的St隨著Re的增大呈單調(diào)上升趨勢，而圓柱的St先增加，在轉(zhuǎn)折點處到達極大值，后隨著Re的增大而減小。方柱的St值一直小于圓柱的St值，這也說明方柱的漩渦脫落頻率比圓柱的漩渦脫落頻率小。圖9(d)給出了圓柱和方柱λ/D隨Re的變化曲線?？煽闯?，圓柱繞流的主波長隨著Re的增大而單調(diào)下降，而方柱繞流的主波長先隨著Re的增大先下降再上升，在轉(zhuǎn)折點附近到達極小值，且方柱繞流的主波長一直大于圓柱繞流的主波長。

圖9　圓柱繞流與方柱繞流參數(shù)對比Fig. 9　Comparison between flow around a cylinder and a square cylinder

圖10(a)、(b)給出了圓柱、方柱阻力和升力系數(shù)振幅隨Re的變化?？煽闯觯ο禂?shù)的振幅比其對應(yīng)的阻力系數(shù)的振幅大一個數(shù)量級，兩者均隨著Re的增大而增大，且在方柱繞流的轉(zhuǎn)折點之前，圓柱和方柱的阻力升力系數(shù)振幅大體一致，轉(zhuǎn)折點之后，方柱的阻力升力系數(shù)振幅顯著增大，圓柱的變化依舊平緩，且方柱的振幅大于圓柱的振幅。圖10(c)給出了Cd(p)/Cd(v)隨Re的變化曲線，可看出Cd(p)相比于Cd(v)一直處于主導(dǎo)地位，且隨著Re的增大越明顯。圓柱繞流Cd(p)/Cd(v)隨Re呈線性增長，而方柱繞流Cd(p)/Cd(v)隨Re呈非線性增長，隨著Re的增加，變化率越來越大，且在方柱繞流中，Cd(p)的主導(dǎo)地位相對圓柱繞流要明顯。

圖10　主要參數(shù)隨Re編號曲線Fig. 10　The variation of main parameters vs Re

3結(jié)論

本文利用自主研發(fā)的CIP-ZJU模型，對低雷諾數(shù)(Re<300)圓柱和方柱繞流問題開展了數(shù)值模擬。模型在直角坐標(biāo)系統(tǒng)下建立，采用緊致插值曲線CIP方法作為流場的基本求解器，離散了Navier-Stokes方程,利用浸入邊界方法處理固體邊界，基于多相流的理論實現(xiàn)流-固耦合同步求解。分別分析了圓柱和方柱繞流的繞流水動力特性，并將兩者的主要水動力參數(shù)進行比較，得出如下結(jié)論：

1)圓柱和方柱繞流的Cd在Re≈140～150達到極小值，此時方柱的分離點由后角點轉(zhuǎn)移到前角點，方柱的摩擦阻力的方向發(fā)生變化，而圓柱沒有此現(xiàn)象發(fā)生，其分離點只圍繞一點往復(fù)運動，且摩擦阻力保持為正值。

2)圓柱繞流的Cd(p)、Cd(v)均為正值，而方柱繞流在轉(zhuǎn)折點之后，由于分離點的變化，導(dǎo)致摩擦阻力的合力方向發(fā)生改變，Cd(v)變?yōu)樨撝怠?/p>

3)方柱繞流的St值一直小于圓柱繞流的St值，說明方柱的漩渦脫落頻率比圓柱的漩渦脫落頻率小。相反方柱繞流的主波長大于圓柱繞流的主波長。

4)圓柱繞流與方柱繞流升力系數(shù)的振幅比其對應(yīng)的阻力系數(shù)的振幅大一個數(shù)量級，兩者均隨著Re的增大而增大，且在方柱繞流的轉(zhuǎn)折點之前，圓柱和方柱的阻力升力系數(shù)振幅大體一致，轉(zhuǎn)折點之后，方柱的阻力升力系數(shù)振幅顯著增大，圓柱的變化依舊平緩，且方柱的振幅大于圓柱的振幅。

本文的計算結(jié)果與文獻結(jié)果吻合良好，驗證了此方法的可行性，為進一步探討柱體形狀對繞流阻力的影響，開展流致振動等研究打下良好基礎(chǔ)。

參考文獻：

[1]TRITTON D J. Experiments on the flow past a circular cylinder at low Reynolds numbers[J]. Journal of fluid mechanics, 1959, 6(4): 547-567.

[2]KAWAGUTI M, JAIN P. Numerical study of a viscous fluid flow past a circular cylinder[J]. Journal of the physical society of Japan, 1966, 21(10): 2055-2062.

[3]RAJANI B N, KANDASAMY A, MAJUMDAR S. Numerical simulation of laminar flow past a circular cylinder[J]. Applied mathematical modelling, 2009, 33(3): 1228-1247.

[4]OKAJIMA A. Strouhal numbers of rectangular cylinders[J]. Journal of fluid mechanics, 1982, 123: 379-398.

[5]SOHANKAR A, DAVIDSON L, NORBERG C. Numerical simulation of unsteady flow around a square two-dimensional cylinder[C]//Proceedings of the Twelfth Australasian Fluid Mechanics Conference，1995.

[6]BREUER M, BERNSDORF J, ZEISER T, et al. Accurate computations of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume[J]. International journal of heat and fluid flow, 2000, 21(2): 186-196.

[7]SEN S, MITTAL S, BISWAS G. Flow past a square cylinder at low Reynolds numbers[J]. International journal for numerical methods in fluids, 2011, 67(9): 1160-1174.

[8]潘小強, 袁璟. CFD軟件在工程流體數(shù)值模擬中的應(yīng)用[J]. 南京工程學(xué)院學(xué)報: 自然科學(xué)版, 2004, 2(1): 62-66.

PAN Xiaoqiang, YUAN Jing. Numerical simulation of engineering fluid by using CFD software[J]. Journal of Nanjing institute of technology: Natural Science Edition, 2004, 2(1): 62-66.

[9]周云龍, 鄧冬, 刁成東,等. 方柱繞流和圓柱繞流旋渦脫落的數(shù)值模擬[C]//中國工程熱物理學(xué)會2008多相流學(xué)術(shù)會議論文集. 青島, 2008: 086020.

ZHOU Yunlong, DENG Dong, DIAO Chengdong et al. Numerical simulation of vortex shedding from a square and circular cylinders[C]// Conferences of multiphase flow of Chinese Society of Engineering Thermophysics. Qingdao, 2008: 086020.

[10]趙西增. 一種自由面大變形流動的CFD模型[J]. 浙江大學(xué)學(xué):工學(xué)版, 2013, 47(8): 1384-1392.

ZHAO Xizeng. CFD modeling of distorted free surface flow[J]. Journal of Zhejiang university: engineering science, 2013, 47(8): 1384-1392.

[11]ZHAO Xizeng, YE Zhouteng, FU Yingnan, et al. A CIP-based numerical simulation of freak wave impact on a floating body[J]. Ocean engineering, 2014, 87: 50-63.

[12]TAKEWAKI H, NISHIGUCHI A, YABE T. Cubic interpolated pseudo-particle method (CIP) for solving hyperbolic-type equations[J]. Journal of computational physics, 1985, 61(2): 261-268.

[13]YABE T, XIAO Feng, UTSUMI T. The constrained interpolation profile method for multiphase analysis[J]. Journal of computational physics, 2001, 169(2): 556-593.

[14]PESKIN C S. Flow patterns around heart valves: A numerical method [J]. Journal of computational physics, 1972, 10(2): 252-271.

[15]HU Changhong, KASHIWAGI M. Two-dimensional numerical simulation and experiment on strongly nonlinear wave-body interactions[J]. Journal of marine science and technology, 2009, 14(2): 200-213.

[16]WILLIAMSON C H K. Defining a universal and continuous Strouhal-Reynolds number relationship for the laminar vortex shedding of a circular cylinder[J]. Physics of fluids, 1988, 31(10): 2742-2744.

[17]SCHLICHTING H, GERSTEN K, KRAUSE E, et al. Boundary-layer theory[M]. Berlin: Springer, 1968.

[18]余化軍. 圓柱和方柱繞流及矩形柱渦激振動的二維數(shù)值模擬分析[D]. 天津: 天津大學(xué), 2011: 17.

YU Huajun. Two dimensional numerical analysis of flow over a circular and square cylinder and vortex-induced vibration of rectangular cylinder[D]. Tianjin: Tianjin University, 2011: 17.

[19]呂林. 海洋工程中小尺度物體的相關(guān)水動力數(shù)值計算[D]. 大連: 大連理工大學(xué), 2006: 66-67.

LYU Lin. Numerical simulations for the hydrodynamics of small scale objects in ocean engineering[D]. Dalian: Dalian University of Technology, 2006: 66-67.

[20]ZDRAVKOVICH M M. Flow around circular cylinder[M]. Oxford: Oxford Science Publications, 1997.

[21]HE J W, GLOWINSKI R, METCALFE R, et al. Active control and drag optimization for flow past a circular cylinder: I. oscillatory cylinder rotation[J]. Journal of computational physics, 2000, 163(1): 83-117.

[22]VOROBIEFF P, GEORGIEV D, INGBER M S. Onset of the second wake: dependence on the Reynolds number[J]. Physics of fluids, 2002, 14(7): 53-56.

[23]SOHANKAR A, NORBERG C, DAVIDSON L. Low-Reynolds-number flow around a square cylinder at incidence: study of blockage, onset of vortex shedding and outlet boundary condition[J]. International journal for numerical methods in fluids, 1998, 26(1): 39-56.

[24]GERA B, SHARMA P K, SINGH R K. CFD analysis of 2D unsteady flow around a square cylinder[J]. International journal of applied engineering research, 2010, 1(3): 602-610.

[25]畢繼紅, 余化軍, 任洪鵬. 靜止方柱和圓柱繞流的二維數(shù)值分析[J]. 三峽大學(xué)學(xué)報: 自然科學(xué)版, 2012, 34(1): 41-45.

BI Jihong, YU Huajun, REN Hongpeng. Two dimensional numerical simulation of flow over a static square cylinder and a static circular cylinder[J]. Journal of China Three Gorges university: natural sciences, 2012, 34(1): 41-45.

Numerical simulation of flow past a cylinder using a CIP-based model

ZHAO Xizeng1，2, FU Yingnan1, ZHANG Dake1

(1. Ocean College, Zhejiang University, Hangzhou 310058, China; 2. Second Institute of Oceanography, SOA State Key Laboratory of Satellite Ocean Environment Dynamics, Hangzhou 310012, China)

Abstract：Flow past a cylinder and the flow-induced vibration phenomenon are complex engineering problems. We developed a CIP-ZJU (constrained interpolation profile-Zhejiang University) model to study the flow past circular and square cylinders for Reynolds numbers Re< 300. The model was established in the Cartesian coordinate system, using the CIP method as the base flow solver to discretise the Navier-Stokes equations. The fluid-structure interaction was treated as multiphase flow, with liquid and solid phases solved simultaneously. We used an immersed boundary method to deal with the boundary of the solid body. We then compared our computations with available results and obtained good agreements. By introducing pressure and viscous drags, we analyzed the variations in the Strouhal number, drag coefficient, and lift coefficient with the Reynolds numbers. Results show that the variations of the dynamic characteristics differ between the flows past circular and square cylinders.

Keywords：flow past cylinder; CIP method; Navier-Stokes equation; immersed boundary method; Strouhal number

中圖分類號：O352

文獻標(biāo)志碼：A

文章編號：1006-7043(2016)03-297-09

doi:10.11990/jheu.201411003

作者簡介：趙西增(1979-), 男, 副教授, 博士生導(dǎo)師；付英男(1990-),男，碩士研究生.通信作者：付英男, E-mail:fyn_zju@163.com.

基金項目：國家自然科學(xué)基金資助項目(51209184, 51479175).

收稿日期：2014-11-02.

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20160104.1427.002.html

網(wǎng)絡(luò)出版日期：2016-01-04.