等邊布置三方柱繞流漩渦相互作用機(jī)制的數(shù)值分析

劉 歡

（天津城建大學(xué) 能源與安全工程學(xué)院，天津 300384）

眾所周知，在一定雷諾數(shù)條件下，黏性流體流經(jīng)柱體時(shí)，在逆壓梯度作用下，流體邊界層發(fā)生分離，引發(fā)周期性的漩渦脫落.這種渦旋的產(chǎn)生和脫落現(xiàn)象給工程帶來(lái)不利的影響，甚至發(fā)生工程事故.例如在1940年，美國(guó)華盛頓州建成剛4個(gè)月的塔科馬海峽大橋，在未達(dá)其設(shè)計(jì)風(fēng)速一半的氣候條件下，因渦激振動(dòng)而造成橋梁主體的共振，最終結(jié)構(gòu)潰塌[1]；飛機(jī)與船只等，在流體中運(yùn)動(dòng)時(shí)，由于尾部所產(chǎn)生的渦旋消耗動(dòng)能，形成航行時(shí)的阻力.這些都是不利于生產(chǎn)實(shí)踐的.同時(shí)還應(yīng)該注意到，有時(shí)流體的渦旋運(yùn)動(dòng)也具有積極的工程意義.例如，西班牙的初創(chuàng)公司研發(fā)了一套全新的風(fēng)力渦輪發(fā)電機(jī)——Vortex Bladeless，這是一套沒(méi)有葉片的渦輪機(jī)，它是靠周圍空氣的運(yùn)動(dòng)形成漩渦，然后帶動(dòng)設(shè)備前后晃動(dòng)起來(lái)，再利用塔架底部的兩個(gè)環(huán)形相斥磁鐵作為非電動(dòng)機(jī)馬達(dá)，通過(guò)底部相斥的磁鐵推動(dòng)時(shí)椎體持續(xù)產(chǎn)生最大震蕩，震蕩時(shí)產(chǎn)生的機(jī)械能就會(huì)透過(guò)錐體底部的交流發(fā)電機(jī)，將力學(xué)能轉(zhuǎn)換為電能.因此，研究柱體繞流渦旋生成、演化及脫落規(guī)律，不僅具有基礎(chǔ)的理論意義，也具有重要的實(shí)際工程價(jià)值.

近年來(lái)，柱體的繞流已經(jīng)從圓柱轉(zhuǎn)向方柱，柱體個(gè)數(shù)及排列方式也多種多樣.張偉等[2]采用有限體積法對(duì)雷諾數(shù)Re為100的等邊布置三方柱的二維繞流問(wèn)題進(jìn)行了數(shù)值模擬.考慮左“品”字形和右“品”字形兩種柱體幾何布置，間距比在1.5～5.0變化、右“品”字形布置時(shí)，對(duì)流場(chǎng)特性與結(jié)構(gòu)受間距比的影響較小，在所研究范圍內(nèi)，未出現(xiàn)偏流現(xiàn)象.在較小間距比時(shí)，右“品”布置的流動(dòng)穩(wěn)定性比左“品”字形布置更穩(wěn)定.

Okajima[3]運(yùn)用實(shí)驗(yàn)研究的方法，研究了Re在70～20 000時(shí)，不同截面縱橫比的矩形截面柱體的斯特勞哈爾數(shù)St與雷諾數(shù)Re的關(guān)系，并詳細(xì)分析了在不同的雷諾數(shù)下，矩形柱的長(zhǎng)高比率不同時(shí)，方柱后尾跡的基本頻率和能譜的變化.

趙心廣[4]對(duì)高雷諾數(shù)下右“品”字等邊布置的三方柱繞流進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)上游兩個(gè)方柱的阻力系數(shù)比下游方柱小，上游兩方柱升力系數(shù)比下游方柱大，并且下游方柱的升力系數(shù)在零附近振蕩；郭明旻等[5]以三方柱模型實(shí)驗(yàn)結(jié)果為依據(jù)，討論了建筑物表面風(fēng)壓脈動(dòng)的極值分布情況，并發(fā)現(xiàn)在間距比2.5～3.0間流場(chǎng)可能存在突變；張愛(ài)社等[6]用有限元方法對(duì)雷諾數(shù)Re為200的等邊布置的二維繞流問(wèn)題進(jìn)行數(shù)值模擬，研究表明三圓柱之間的干擾在小間距比情況下比較嚴(yán)重，流動(dòng)偏向下游的某個(gè)圓柱.

筆者對(duì)Re為100和200的低雷諾數(shù)條件下，等邊左“品”字形布置的三方柱二維繞流問(wèn)題進(jìn)行數(shù)值模擬研究.考慮到間距過(guò)小三方柱會(huì)出現(xiàn)重合，而間距過(guò)大三方柱間又不再有相互影響，所以取間距比為1.2～6.0.通過(guò)局部區(qū)間進(jìn)行間距比的加密，對(duì)流場(chǎng)流線、渦量圖以及阻力系數(shù)、升力系數(shù)、斯特勞哈爾數(shù)的分析，研究流場(chǎng)流動(dòng)特性、柱體阻力系數(shù)、升力系數(shù)及斯特勞哈爾數(shù)的變化規(guī)律，特別關(guān)注了對(duì)流場(chǎng)流動(dòng)和動(dòng)力特性具有重要影響的臨界間距比.

1 數(shù)值研究方法

1.1 控制方程和相關(guān)參數(shù)

采用有限體方法，求解二維黏性不可壓非定常流場(chǎng)；壓力采用二階離散格式，動(dòng)量采用二階迎風(fēng)格式，壓力與速度耦合采用SIMPLE算法.

二維黏性不可壓非定常流動(dòng)的控制方程為連續(xù)性方程和不可壓Navier-Stokes方程

方程中的下標(biāo)采用了求和約定.vi、p、t分別為速度分量、壓力和時(shí)間；ρ為流體密度；μ為流體動(dòng)力黏性系數(shù).

阻力系數(shù)、升力系數(shù)、間距比及St數(shù)等無(wú)量綱數(shù)的定義為

式中：CD為阻力系數(shù)，正方向與來(lái)流方向相反；CL為升力系數(shù)，正方向與來(lái)流方向垂直；FD為阻力；FL為升力；ρ為流體密度；U為來(lái)流速度；D為方柱邊長(zhǎng)；g為間距比；L為方柱中心之間的距離；St為斯特勞哈爾數(shù)；f為渦脫落頻率.

1.2 計(jì)算域及邊界條件

三方柱布置形式及計(jì)算區(qū)域如圖1所示.計(jì)算域?yàn)殚L(zhǎng)方形，三方柱呈等邊左“品”字型布置，上游為單方柱A，下游為雙方柱B、C并列.方柱邊長(zhǎng)為無(wú)量綱邊長(zhǎng)D=1，A方柱形心至入流邊界距離為20D，B、C方柱形心至出流邊界的距離為40D，上下邊界距離最近方柱中心為10D.入口為速度邊界條件，來(lái)流速度為無(wú)量綱數(shù)U=1，出口取自由出流邊界條件，上下壁面和方柱表面為Wall邊界條件，參見(jiàn)圖1.

圖1 等邊三方柱布置形式及計(jì)算域示意

1.3 網(wǎng)格劃分

在柱體繞流中，柱體周圍靠近柱體的地方需要更高的求解精度.為了保證本研究的求解精度，驗(yàn)證算例的網(wǎng)格劃分采用以下形式：在方柱外圍劃分一層越靠近方柱越密的矩形網(wǎng)格，目的是在數(shù)值模擬過(guò)程中，保證流體參數(shù)變化最為激烈的區(qū)域的求解精度；再在這層矩形網(wǎng)格的外圍用非結(jié)構(gòu)三角形網(wǎng)格對(duì)求解區(qū)域進(jìn)行劃分，并通過(guò)確定邊界網(wǎng)格節(jié)點(diǎn)的方式，控制網(wǎng)格的疏密程度，越靠近柱體的地方網(wǎng)格越致密，如圖2所示.為了加快收斂過(guò)程，時(shí)間導(dǎo)數(shù)項(xiàng)采用了二階隱式離散格式.

圖2 方柱網(wǎng)格劃分示意

2 結(jié)果與討論

2.1 驗(yàn)證算例

選取雷諾數(shù)為100、200時(shí)的單方柱繞流算例，與經(jīng)典的計(jì)算和實(shí)驗(yàn)結(jié)果進(jìn)行比較，驗(yàn)證本文所采用的數(shù)值方法和求解參數(shù)的有效性.表1為本文計(jì)算的時(shí)間平均阻力系數(shù)和斯特勞哈爾數(shù)與相關(guān)文獻(xiàn)中的數(shù)值計(jì)算和實(shí)驗(yàn)結(jié)果的比較.從表1可以看到，本文的計(jì)算結(jié)果與經(jīng)典結(jié)果比較吻合，證明本文的數(shù)值方法是合理有效的.

表1 單方柱繞流模擬結(jié)果比較

2.2 流動(dòng)特性

對(duì)雷諾數(shù)為100、200的左品字三方柱繞流流場(chǎng)進(jìn)行模擬，并進(jìn)行間距比局部加密模擬.圖3-4為t=300 s、Re為 100和 200時(shí)間距比 g分別為 1.2、1.5、2.0、2.5、3.0、3.5、4.0、4.5、5.0、5.5、6.0 下的渦量圖.

圖3 Re為100時(shí)不同間距比的渦量圖

由圖3可知：當(dāng)g≤1.5時(shí)，三方柱非常接近，此時(shí)流場(chǎng)形態(tài)類似于一個(gè)大的孤立方柱，在尾流區(qū)形成一組漩渦；當(dāng)g＞1.5時(shí)，上流柱體脫落的剪切層附著于下流方柱，出現(xiàn)偏流現(xiàn)象，方柱B、C后產(chǎn)生了各自的漩渦，漩渦脫落發(fā)生在一個(gè)方柱后；兩方柱之間的間隙流呈現(xiàn)出雙穩(wěn)態(tài)偏流現(xiàn)象；在g為2.0～3.0時(shí)，流場(chǎng)十分不穩(wěn)定，在方柱A后出現(xiàn)了回流區(qū)，漩渦發(fā)生脫落；隨著g進(jìn)一步增大，偏流現(xiàn)象減弱；在g為3.5～4.5之間時(shí)，漩渦發(fā)生相位的相對(duì)偏移，上述現(xiàn)象與文獻(xiàn)[1]的現(xiàn)象大致相同；當(dāng)g≥5.0時(shí)，方柱B和C后的漩渦卻呈現(xiàn)出反相位脫落；當(dāng)g=6.0時(shí)，三方柱后產(chǎn)生各自漩渦，相互之間幾乎不再有影響.

圖4 Re為200時(shí)不同間距比的渦量圖

由圖4可知：當(dāng)g=1.2時(shí)，已經(jīng)出現(xiàn)偏流現(xiàn)象，間隙流對(duì)方柱的影響較大，方柱B、C后差生了各自的漩渦，漩渦脫落發(fā)生在一個(gè)方柱后，兩方柱之間的間隙流呈現(xiàn)出雙穩(wěn)態(tài)偏流現(xiàn)象；在g為2.0～3.0時(shí)，流場(chǎng)十分不穩(wěn)定，在方柱A后出現(xiàn)了回流區(qū)，漩渦發(fā)生脫落，偏流現(xiàn)象逐漸增強(qiáng)，這是因?yàn)殡S著g進(jìn)一步增大，間隙流對(duì)后方對(duì)稱的方柱B、C影響趨于相同，同時(shí)前方柱A與后方柱B、C之間流場(chǎng)的相互干擾逐漸減小，偏流現(xiàn)象慢慢減弱；當(dāng)g=3.5時(shí)，下游方柱后的漩渦出現(xiàn)反相位脫落，流場(chǎng)為同步對(duì)稱結(jié)構(gòu)，偏流現(xiàn)象消失；在g為4.5～5.5之間時(shí)，漩渦始終處于反相位同步脫落；當(dāng)g=6.0時(shí)，三方柱后產(chǎn)生各自漩渦，相互之間干擾現(xiàn)象并未消失.

2.3 Re為100時(shí)阻力系數(shù)和升力系數(shù)特性

圖5給出Re=100時(shí)平均阻力系數(shù)CD隨g的變化曲線.由圖5可知：①在g為1.2～1.8區(qū)間內(nèi)，上游方柱A所受的阻力略大于下游方柱B、C所受的阻力；當(dāng)g=1.5時(shí)，方柱B、C的CD差值最大；②當(dāng)g=2.0時(shí)，方柱A的CD值出現(xiàn)了最小值CD，A=1.10（這里下標(biāo)表示具體的方柱，下同）；當(dāng)g＞2.0時(shí)，三方柱的CD隨 g 的增加而增大，之后隨著 g 的增大，CD，B、CD，C始終大于 CD，A；g=2.5 時(shí)，CD，B=1.63，CD，C=1.60，接近于單方柱 CD，CD，A=1.22 依然小于單方柱 CD；③在 g＞3.0 時(shí)，B、C柱CD幾乎相等，平緩變化，這與流線圖中的變化相吻合；當(dāng)g=4.0時(shí)，下游方柱B、C的CD約為單柱繞流時(shí)的1.15倍，與文獻(xiàn)[1]的1.17倍大致相同；當(dāng)g達(dá)到最大值5.5時(shí)，各柱的CD也分別達(dá)到最大值，CD，A=1.46，CD，B=1.89，CD，C=1.89，上游柱 A 的 CD，A依然小于單方柱繞流時(shí)的CD=1.61，下游B、C柱的CD比單方柱繞流的情形大.

圖5 Re為100時(shí)CD-g的變化曲線

圖5中的現(xiàn)象出現(xiàn)的原因主要是：①g過(guò)小時(shí)，A的自由剪切層會(huì)將下游的方柱B、C完全包圍，來(lái)流阻力將主要作用在柱A上；②隨著g的增加，下游方柱B、C已不再全部被上游柱A的自由剪切層吞沒(méi)，來(lái)流會(huì)部分撞擊到方柱B、C上；③g繼續(xù)增大，B、C柱的來(lái)流阻力減小，而上游柱A的來(lái)流阻力增大，并未超過(guò)單柱繞流的數(shù)值；即使g較大，方柱A與方柱B、C之間區(qū)域的壓力也遠(yuǎn)大于單柱繞流的情形，是由于下游并列方柱B、C的“夾心”效應(yīng).

圖6給出Re=100時(shí)平均升力系數(shù)CL隨g的變化曲線.

圖6 Re=100時(shí)CL-g的變化曲線

由圖6可以看出：上游方柱A的CL，A始終接近于常數(shù)0，與單方柱的CL一致（CL=-0.001），這表明下游方柱B、C周圍的流動(dòng)對(duì)上游柱A周圍的流場(chǎng)已基本無(wú)干擾；下游方柱 B、C 的 CL，B和 CL，C大小幾乎相等，方向相反，這表明柱體間具有更加強(qiáng)烈的橫向干擾效應(yīng)；其中g(shù)=1.5時(shí)，方柱A的CL，A突然增大，達(dá)到最大值；隨著 g 的增加，CL，B和 CL，C對(duì)稱地趨近于 0.綜合分析繞流流場(chǎng)和CL變化規(guī)律可以推測(cè)，臨界間距比應(yīng)該在3.5附近.

2.4 Re為200時(shí)阻力系數(shù)和升力系數(shù)特性

圖7給出Re=200時(shí)CD隨g的變化曲線.

圖7 Re=200時(shí)CD-g的變化曲線

由圖7可以看出，Re為200時(shí)各方柱CD隨g的變化比100時(shí)的劇烈：①在g為1.2～1.8區(qū)間內(nèi)，上游方柱A所受的阻力略大于下游方柱B、C所受的阻力，兩種雷諾數(shù)下的變化趨勢(shì)一致；當(dāng)g=1.5時(shí)，柱B、C的 CD差值最大；當(dāng) g=2.0 時(shí)，CD，A=1.08（這與 Re 為100時(shí)，當(dāng)g=1.9時(shí)，柱A的CD值出現(xiàn)了一個(gè)下跌CD，A=1.10 相似），之后隨著 g 的增大，CD，B、CD，C始終大于 CD，A；②隨 g的增加，三方柱的 CD增大（g=1.9 除外）；當(dāng) g=2.2 左右時(shí)，CD，B、CD，C接近于單方柱阻力系數(shù)，CD，A依然小于單方柱阻力系數(shù)，這與Re為100時(shí)g=2.5處的現(xiàn)象相似；隨后隨著g的增加，三方柱的CD均增大，但A方柱的CD依然小于單方柱阻力系數(shù)；③在g=3.5時(shí)，CD突然增加，隨后降低，平緩變化，這與Re為100時(shí)，在g=3.7處的變化相似，這與流線圖中的突變相吻合；當(dāng)g=3.5時(shí)，B、C柱的CD也分別達(dá)到最大值，CD，B=2.15，CD，C=2.19，顯然，CD，A＞CD，B≈CD，C＞CD=1.63.

圖8給出Re為200時(shí)CL隨g的變化曲線.

圖8 Re=200時(shí)CL-g的變化曲線

由圖8可以看出：Re為200時(shí)CL隨g的變化比Re為100時(shí)的變化劇烈；上游柱A的CL，A始終與單方柱的平均升力系數(shù)（CL=0.004）一致，接近于常數(shù)0，這表明下游方柱B、C周圍的流動(dòng)對(duì)上游柱A的流場(chǎng)幾乎沒(méi)有影響；下游并列的柱 B、C 的 CL，B和 CL，C大小幾乎相等，方向相反，這表明柱體間具有更加強(qiáng)烈的橫向干擾效應(yīng)，這與Re為100時(shí)相似；其中g(shù)=3.5時(shí)，CL，B和CL，C呈現(xiàn)大小相等、方向相反的對(duì)稱趨勢(shì).綜合繞流流場(chǎng)和CL變化規(guī)律可以推測(cè)，臨界間距比在4.0附近.

2.5 臨界間距比的驗(yàn)證

為準(zhǔn)確捕捉到臨界間距比的數(shù)值，分別對(duì)g=3.5（Re=100）和 g=4（Re=200）附近進(jìn)行加密.圖 9 給出了間距比加密區(qū)間內(nèi)CD和CL隨g的變化曲線.

圖9 不同Re時(shí)CD-g與CL-g的變化曲線

由圖9a-9b可以看出：Re為100、g為3.0～4.5時(shí)，變化平緩，方柱A與方柱B、C的CD的變化趨勢(shì)一致，CD，B=CD，C，CL，B=|CL，C|；由圖 9c 可知，Re 為 200、g 為3.0～4.4，方柱A與方柱B、C的CD的變化趨勢(shì)一致，在g為3.5與3.9時(shí)，三方柱CD突然增大.

圖10 不同Re時(shí)St-g的變化曲線

圖10給出斯特勞哈爾數(shù)St隨g的變化曲線.由圖 10a-10b中可以看出：Re為 100、g為2.0～3.0 之間時(shí)，St隨g的增大變化幅度較大；隨后，St趨于一致，隨g的變化趨勢(shì)與文獻(xiàn)[2，13]基本相同，與單方柱繞流的 St=0.145接近；g≥3.5時(shí)，方柱 A、B、C 的 St與單方柱情形基本相等，在g=3.9處出現(xiàn)突然增大，方柱A、B、C的St幾乎相等，這與CL變化趨勢(shì)一致.說(shuō)明此時(shí)三方柱之間的干擾最小，也就是柱B和柱C間偏流消失的原因；當(dāng)g≥5.0時(shí)，三方柱之間的相互干擾現(xiàn)象消失.

由圖10c-10d可以看出：Re為200、g=2.0～3.0之間，St隨g的增大變化幅度較大；隨后，St趨于一致，趨近單方柱繞流的St=0.163，與Re為100時(shí)的變化規(guī)律相似；g≥3.7時(shí)，方柱A、B、C的St幾乎相等；g=4.1處出現(xiàn)突然增大；當(dāng)g=6.0時(shí)，三方柱之間始終存在相互干擾現(xiàn)象.

圖11 Re為100、g=3.7時(shí)流場(chǎng)渦量驗(yàn)證圖

圖11-12為臨界間距比的流場(chǎng)渦量驗(yàn)證圖，即Re為100、g=3.7、不同時(shí)間的渦量變化圖11和Re為200、g=4.1、不同時(shí)間的渦量變化圖12.

圖12 Re為200、g=4.1時(shí)流場(chǎng)渦量驗(yàn)證圖

綜合分析圖9-12可以確定，Re為100時(shí)，偏流現(xiàn)象消失的臨界間距比g=3.7；Re為200時(shí)，偏流現(xiàn)象消失的臨界間距比g=4.1.

3 結(jié)論

（1）當(dāng)雷諾數(shù)Re為100、間距比g≥1.5時(shí)，出現(xiàn)雙穩(wěn)態(tài)間隙偏流；g≥2.5時(shí)，偏流現(xiàn)象減弱；當(dāng)g=3.7時(shí)，方柱B、C后的漩渦出現(xiàn)同相位脫落現(xiàn)象，流場(chǎng)為同步對(duì)稱結(jié)構(gòu)，偏流現(xiàn)象消失；當(dāng)g≥4.9時(shí)，方柱B、C后的漩渦出現(xiàn)反相位脫落現(xiàn)象；當(dāng)g=5.0時(shí)，流場(chǎng)呈對(duì)稱分布；當(dāng)g≥5.0時(shí)，三方柱之間的相互干擾現(xiàn)象消失；偏流現(xiàn)象消失的臨界間距比為3.7，且g=3.7時(shí)，方柱A、B、C的斯特勞哈爾數(shù)St出現(xiàn)突然增大，且?guī)缀跸嗟?

（2）當(dāng)Re為200、g≥1.2時(shí)，出現(xiàn)雙穩(wěn)態(tài)間隙偏流；g≥2.5時(shí)，偏流現(xiàn)象減弱；當(dāng)g=3.5時(shí)，方柱B、C后的漩渦出現(xiàn)反相位脫落現(xiàn)象，流場(chǎng)為同步對(duì)稱結(jié)構(gòu)，偏流現(xiàn)象消失；當(dāng)g≥4.5時(shí)，方柱B、C后的漩渦卻呈現(xiàn)出反相位脫落；在g為4.5～5.5之間，漩渦始終處于反相位同步脫落；當(dāng)g=6.0時(shí)，三方柱之間始終存在相互干擾現(xiàn)象；偏流現(xiàn)象消失的臨界間距比為4.1，且g≥4.1，方柱 A、B、C 的 St幾乎相等.

（3）Re為200時(shí)的流場(chǎng)比100時(shí)的更不穩(wěn)定，隨著雷諾數(shù)的增加，偏流現(xiàn)象消失的臨界間距比增大.

參考文獻(xiàn)：

[1]余化軍.圓柱和方柱繞流及矩形柱渦激振動(dòng)的二維數(shù)值分析[D].天津大學(xué)，2011：41-45.

[2]張 偉，戴玉滿.低雷諾數(shù)下等邊布置三方柱繞流的數(shù)值研究[J].機(jī)械工程學(xué)報(bào)，2015（12）：185-191.

[3]OKAJIMA A.Strouhal numbers of rectangular cylinders[J].Journal of Fluid Mechanics，1982，123（123）：379-398.

[4]趙心廣.三方柱繞流的大渦模擬及頻譜分析[J].水科學(xué)與工程技術(shù)，2011（4）：28-30.

[5]郭明旻，黃東群，林 晨，等.三方柱表面風(fēng)壓脈動(dòng)的時(shí)域特性[J].力學(xué)季刊，2001，22（4）：439-446.

[6]張愛(ài)社，張 陵.等邊布置三圓柱繞流的數(shù)值分析[J].應(yīng)用力學(xué)學(xué)報(bào)，2003，20（1）：31-36.

[7]王掩剛，陳俊旭，先松川.基于POD方法的二維方柱低雷諾數(shù)繞流流場(chǎng)分析研究[J].西北工業(yè)大學(xué)學(xué)報(bào)，2014（4）：612-617.

[8]DAVIS R W，MOORE E F，PURTELL L P.A numericalexperimental study of confined flow around rectangular cylinders[J].Physics of Fluids，1984，27（1）：46-59.

[9]SOHANKAR A，DAVIDSON L，NORBERG C.Numerical simulation of unsteady flow around a square two-dimensional cylinder[C]//The University of Sydney.Australasian Fluid Mechanics Conference.Sydney：The University of Sydney，1995：517-520.

[10]ISLAM S U，ZHOU C Y.Numerical simulation of flow around a row of circular cylinders using the lattice Boltzmann method[J].Information Technology Journal，2009，8（4）：513-520.

[11]SAHA A K，BISWAS G，MURALIDHAR K.Three-dimensional study of flow past a square cylinder at low Reynolds numbers[J].International Journal of Heat&Fluid Flow，2003，24（1）：54-66.

[12]WAN D C，PATNAIK B S V，WEI G W.Discrete singular convolution-finite subdomain method for the solution of incompressible viscous flows[J].Journal of Computational Physics，2002，180（1）：229-255.

[13]沈立龍，劉明維，呂啟兵，等.雙排并列三方柱繞流數(shù)值模擬[J].科學(xué)技術(shù)與工程，2014，14（23）：135-139.