韓兆龍，周 岱，陳亞楠，桂曉瀾，李俊龍

（上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240）

譜單元法及其在多圓柱繞流分析中的應(yīng)用

韓兆龍，周 岱，陳亞楠，桂曉瀾，李俊龍

（上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240）

系統(tǒng)闡明譜單元方法，基于譜單元方法對(duì)低雷諾數(shù)Re＝200時(shí)不同間距下的順排兩圓柱和Re＝150正方形排列的四圓柱繞流及其阻力系數(shù)、升力系數(shù)等進(jìn)行數(shù)值模擬。研究比較不同間距比L／D（兩圓柱圓心距離與圓柱直徑之比）對(duì)兩圓柱和四圓柱繞流的影響，計(jì)算分析了渦量圖分布、平均阻力系數(shù)和斯托羅哈數(shù)隨間距比的變化。研究表明，間距比對(duì)順排兩圓柱和正方形四圓柱繞流影響顯著；順排兩圓柱繞流存在臨界間距比，在Re＝200時(shí)臨界間距比約為3．6。正方形排列四圓柱存在三種流態(tài)。當(dāng)流場(chǎng)從一種流動(dòng)形態(tài)變成另外一種流動(dòng)形態(tài)時(shí)，力學(xué)參數(shù)發(fā)生顯著變化，在某些間距比區(qū)間內(nèi)出現(xiàn)驟升或驟降現(xiàn)象。

譜單元法；順排兩圓柱；正方形四圓柱；間距比；流動(dòng)形態(tài)

0 引 言

自從1977年Gottlieb和Orszag［1］系統(tǒng)性地從數(shù)學(xué)方面對(duì)譜方法進(jìn)行了理論闡述以來，譜方法被廣泛地應(yīng)用于更多的領(lǐng)域。1984年，Gottlieb和Hussaini開始將譜方法向計(jì)算流體動(dòng)力學(xué)推廣［2］。20世紀(jì)80年代初期，Patera結(jié)合譜方法的精度和有限元的思想提出了所謂的譜單元法［3］。譜單元法是譜方法與有限元方法結(jié)合起來，將求解區(qū)域剖分成單元，采用等參元，在每個(gè)單元內(nèi)用N階多項(xiàng)式展開來近似變量，對(duì)方程采用Galerkin方法，用Gauss點(diǎn)求積分，得到離散化的方程組。當(dāng)單元數(shù)目固定時(shí)，計(jì)算精度隨多項(xiàng)式階數(shù)N的增加而增加。因此，譜單元法是求解偏微分方程的一種高階加權(quán)余量法，它繼承了有限元方法對(duì)復(fù)雜幾何區(qū)域的適應(yīng)性和譜方法的高精度與收斂特性等優(yōu)點(diǎn)。

圓柱繞流的理論研究和工程應(yīng)用意義重大。順排兩圓柱或多圓柱是常見的排列形式，例如發(fā)電廠冷卻塔群等。圓柱尾流中的旋渦脫落對(duì)圓柱繞流的周期性升力、阻力產(chǎn)生影響，引起壓力脈動(dòng)、切應(yīng)力脈動(dòng)和升阻力系數(shù)的脈動(dòng)。相對(duì)于單圓柱的繞流，兩圓柱或多圓柱繞流及其流致受力更為復(fù)雜?，F(xiàn)今，有限體積法、Boltzmann等多種數(shù)值模擬方法可用于分析圓柱繞流，獲取柱體表面上的作用力和關(guān)鍵參數(shù)［4-6］。

本文較系統(tǒng)闡述譜單元方法及其表達(dá)，運(yùn)用譜單元方法數(shù)值模擬低雷諾數(shù)Re＝200下的順排兩圓柱、Re＝150正方形排列四圓柱的繞流及其互擾效應(yīng)，分析不同間距比對(duì)流場(chǎng)特性的影響，揭示平均升阻力系數(shù)和斯托羅哈數(shù)（旋渦脫落頻率）等隨間距比的變化，并解釋其形成機(jī)制。

本文的數(shù)值計(jì)算采用由Monash大學(xué)的Blackburn教授所提供的Semtex程序，對(duì)二維計(jì)算，該程序采用Gauss-Lobatto-Legendre（GLL）多項(xiàng)式作為譜單元基函數(shù)，時(shí)間和空間的數(shù)值精度可由使用人員自由控制。該方法已成功應(yīng)用于多項(xiàng)科學(xué)研究［7-10］。

1 基本理論和計(jì)算公式

二維粘性不可壓縮非定常流動(dòng)的控制方程為Navier-Stokes方程：

采用基于混合剛性穩(wěn)定格式的高階時(shí)間分裂算法［7］，可將Navier-Stokes方程的求解轉(zhuǎn)化為如下三部分：

流體壓力滿足泊松方程：

對(duì)流體域的空間離散，采用譜單元法，且選取Legendre多項(xiàng)式［7］的極值點(diǎn)作為插值點(diǎn)。

Φ（Ω）（Ω上任意連續(xù)函數(shù)構(gòu)成的空間）中的任何連續(xù)實(shí)函數(shù)u（x）可表示為：

插值多項(xiàng)式

（1）在［-1，1］以外為0；（2）在［-1，1］內(nèi)，φi（xj）＝δij。其中LN（x）為L(zhǎng)egendre多項(xiàng)式。

若將式（4）、（5a）、（5b）、（6a）寫成變分形式，則有：

對(duì)計(jì)算區(qū)域Ω采取等參變換，且取φ＝φj，Navier-Stokes方程的離散格式變?yōu)椋?/p>

其中

選取Legendre多項(xiàng)式，并利用Gauss積分，各系數(shù)計(jì)算公式為

其中

其中，［φb］表示相應(yīng)變量（壓力或速度）在計(jì)算區(qū)域邊界點(diǎn)上的未知值，［φi］表示相應(yīng)變量是在計(jì)算區(qū)域的內(nèi)點(diǎn)上的未知值。

先由下式解出邊界上的未知物理量值：

再由下式解出計(jì)算區(qū)域內(nèi)點(diǎn)的未知物理量值：

若求解的方程組維數(shù)較高，可采用對(duì)單元（而不是對(duì)計(jì)算區(qū)域）求解的策略，即先采用類似于式（19）的方程求得各個(gè)單元邊界上的未知量值，而后再對(duì)各個(gè)單元分別求解其內(nèi)點(diǎn)上的未知物理量值。

2 單圓柱繞流模擬與數(shù)值驗(yàn)證

為驗(yàn)證計(jì)算的準(zhǔn)確性和可靠性，茲選用圖1示的單圓柱繞流問題進(jìn)行數(shù)值方法驗(yàn)證?？紤]直徑為D的圓柱受未經(jīng)擾動(dòng)的均勻來流作用，基于圓柱直徑和來流流速的雷諾數(shù)取Re＝200。選定計(jì)算域?yàn)?0D×40D，圓柱位于坐標(biāo)系原點(diǎn)（0，0）。入口邊界和出口邊界分別位于圓柱中心上游20D和下游30D處，流域頂部和底部離圓柱中心20D。相應(yīng)的邊界條件為：進(jìn)口處自由來流速度為特征速度，即u＝U∞＝1，v＝0．0；上下邊界條件與進(jìn)口邊界條件相同；出口邊界處壓力均為0．0，即p＝0．0；圓柱表面處為不可滑移邊界條件，即u＝0．0，v＝0．0。計(jì)算域網(wǎng)格劃分采用四邊形非結(jié)構(gòu)化譜單元網(wǎng)格，網(wǎng)格單元數(shù)為354，如圖1（a）所示。在靠近圓柱壁面的地方進(jìn)行幾層細(xì)化網(wǎng)格加密，離圓柱壁面最近的一層網(wǎng)格厚度為0．1D，如圖1（b）示。同時(shí)，在圓柱尾流區(qū)域亦進(jìn)行加密處理。

圖1 單圓柱繞流的譜單元網(wǎng)格示意Fig.1 Mesh system of spectral element for cross-flow around a single cylinder

研究圓柱繞流流場(chǎng)常用的無量綱化系數(shù)包括阻力系數(shù)CD、升力系數(shù)CL和斯托羅哈數(shù)St，定義為：

其中，F(xiàn)D為阻力，與來流方向一致，主要由流體圓柱前后壓力差和繞圓柱柱表面摩擦阻力造成；FL為升力，與來流方向垂直，主要由渦交替從圓柱上下表面脫落產(chǎn)生上下表面壓力脈動(dòng)造成；fs為旋渦脫落頻率，D為圓柱直徑。

為分析插值函數(shù)的階數(shù)對(duì)計(jì)算結(jié)果的影響，對(duì)單圓柱繞流的三種計(jì)算工況進(jìn)行基于三種不同階數(shù)的插值函數(shù)的數(shù)值模擬。計(jì)算工況1中，譜函數(shù)插值采用N＝5階Gauss-Lobatto-Legendre多項(xiàng)式形函數(shù)；計(jì)算工況2中，N＝7；工況3中，N＝9。計(jì)算時(shí)間步長(zhǎng)Δi＝0．005。從表1中看出，采取高階譜插值的譜單元法的計(jì)算精度和效率較高，即使對(duì)較粗糙網(wǎng)格劃分，當(dāng)N＝5、7、9時(shí)，亦可得到滿意的結(jié)果，且與既有文獻(xiàn)中的實(shí)驗(yàn)結(jié)果和數(shù)值模擬結(jié)果吻合。為取得計(jì)算精度和計(jì)算效率的平衡，本文取N＝7。這樣，計(jì)算的空間精度為7階，時(shí)間上的精度為2階。

需要說明的是，盡管在本單圓柱擾流的算例中Re＝200有弱三維效應(yīng)，但是本研究采用二維計(jì)算，仍能得到與其他文獻(xiàn)接近的結(jié)果。

表1 平均阻力系數(shù)和斯托羅哈數(shù)（雷諾數(shù)Re＝200）Table 1 Average drag coefficients and Strouhal numbers（Re＝200）

3 順排兩圓柱繞流數(shù)值模擬

將兩圓柱順排排列放置在流體中，通過改變圓柱軸線之間的距離L，分析流場(chǎng)的變化。取Re＝200，計(jì)算時(shí)間步長(zhǎng)Δi＝0．005，采用N＝7階插值函數(shù)，數(shù)值模擬間距比L／D＝1．2、1．5、2．0、3．6、4．0、6．0、8．0和10．0情形下的流場(chǎng)。

3.1 計(jì)算域和邊界條件

流場(chǎng)計(jì)算域和邊界條件與上述單個(gè)圓柱繞流數(shù)值模擬情形相同，圖2為計(jì)算域和邊界條件示意圖。圖3（a）為兩圓柱圓心間距為2．5D時(shí)的網(wǎng)格劃分示意圖，共劃分432個(gè)網(wǎng)格，在靠近圓柱壁面的幾層進(jìn)行細(xì)化網(wǎng)格加密處理，離圓柱壁面最近的一層網(wǎng)格厚度為0．1D，如圖3（b）所示。同時(shí)，在圓柱尾流區(qū)域進(jìn)行加密處理。

圖2 雙圓柱順排排列的計(jì)算域和邊界條件示意圖Fig.2 Schematic diagram of the computational domain and boundary conditions

圖3 兩圓柱繞流的譜單元網(wǎng)格示意圖Fig.3 Mesh system of spectral element for cross-flow around two cylinders

3.2 不同間距比下的流態(tài)

圖4為不同間距比下的順排兩圓柱繞流渦量圖，從圖4中看出，當(dāng)兩圓柱串列時(shí)存在一個(gè)臨界間距，小于臨界間距時(shí)，上游圓柱后不存在旋渦脫落；大于臨界間距時(shí)，上游圓柱后發(fā)生周期性旋渦脫落。據(jù)圖4可定性判定出此臨界間距比大約為L(zhǎng)／D≈3．6；即當(dāng)L／D＜3．6時(shí)，上游圓柱后沒有產(chǎn)生渦脫落，上游圓柱繞流的分離分界層附著在下游圓柱上，而下游圓柱的后側(cè)產(chǎn)生渦脫落，尾流模式類似于單圓柱繞流尾流模式，當(dāng)L／D＝3．6時(shí)，上游圓柱后開始形成渦脫落，當(dāng)L／D＝4時(shí)，上游圓柱的后側(cè)已形成Kármán渦街，且旋渦脫落后撞擊下游圓柱前表面，干擾下圓柱的旋渦脫落模式，隨著間距比的繼續(xù)增大，上游圓柱后的脫落渦繼續(xù)影響下游圓柱的尾流。Re＝200時(shí)，圓柱繞流發(fā)生三維效應(yīng)；可以看到上游圓柱后的渦脫落穩(wěn)定，但下游圓柱后的尾流渦脫落形態(tài)變得復(fù)雜，因此二維數(shù)值模擬不能精確確定臨界間距比周圍的互擾區(qū)域。

圖4 雙圓柱不同間距比下的渦量云圖Fig.4 Vorticity contours with different spacing ratios for the two cylinders

3.3 圓柱表面受力特性隨間距比的變化

圖5為兩圓柱繞流平均阻力系數(shù)隨間距比的變化。從圖5可以看到，兩圓柱的平均阻力系數(shù)均小于單個(gè)圓柱的平均阻力系數(shù)，圓柱1（上游圓柱）的平均阻力系數(shù)要大于圓柱2（下游圓柱）的平均阻力系數(shù)。對(duì)于小間距比下圓柱1，先觀察到其平均阻力系數(shù)隨著間距比的增大而減小，間距比L／D＝3．6的時(shí)取得最小值；然后隨著間距比增大，平均阻尼增大，在L／D＝3．6和L／D＝4之間劇增后，其值緩慢增大，逐漸接近于單個(gè)圓柱的平均阻力系數(shù)。對(duì)于圓柱2，間距比L／D＜3．6時(shí)，其值隨著間距比增大而增大，在L／D＝3．6和L／D＝4之間劇增，從負(fù)值變到相當(dāng)大的正值，之后又緩慢減小，趨于穩(wěn)定。從中易見臨界間距比的存在（3．6＜L／D＜4．0）。當(dāng)間距比小于臨界值時(shí)，上游圓柱的分離剪切層附著在下游圓柱上，渦脫落現(xiàn)象只在下游圓柱出現(xiàn)，此時(shí)下游圓柱受到較小的阻力，且為負(fù)值。當(dāng)間距比大于臨界值時(shí)，上下游圓柱均產(chǎn)生渦脫落，且升力振幅均驟然增大（圖4g），平均阻力系數(shù)突然變大，下游圓柱變化更劇烈。之后，隨著間距比繼續(xù)增大，上下游圓柱的渦脫落均呈現(xiàn)出穩(wěn)定性，變化緩慢。

圖5 雙圓柱平均阻力系數(shù)隨間距比的變化Fig.5 Variation of mean drag coefficients with different spacing ratios for the two cylinders

3.4 斯托羅哈數(shù)隨間距比的變化

本文所有變量均經(jīng)過無量綱化處理，因此斯托羅哈數(shù)即等于渦脫落頻率。渦脫落頻率可通過對(duì)升力曲線作FFT分析求得。

圖6為斯托羅哈數(shù)隨間距比的變化?？傮w上，兩圓柱的斯托羅哈數(shù)基本保持相等，均隨著間距比的增大先減小再增大，并趨向于單個(gè)圓柱繞流時(shí)的斯托羅哈數(shù)。但間距比在某些區(qū)間內(nèi)對(duì)斯托羅哈數(shù)的影響顯著，即間距比從3．6增大到4．0時(shí)，斯托羅哈數(shù)劇烈增大。當(dāng)間距比小于臨界值時(shí)，上游圓柱后無渦脫落，僅下游圓柱后產(chǎn)生渦脫落；當(dāng)間距比大于臨界值時(shí)，上游圓柱和下游圓柱均存在周期性的渦脫落，從而使得斯托羅哈數(shù)在臨界間距之間發(fā)生顯著跳躍。同樣，在雙圓柱擾流的算例中Re＝200有弱三維效應(yīng)，二維計(jì)算，仍能得到與其他文獻(xiàn)較好擬合的結(jié)果。

圖6 斯托羅哈數(shù)隨間距比的變化Fig.6 Variation of the Strouhal numbers，St，with different spacing ratios，L／D for the two tandem cylinders

4 正方形順排排列四圓柱繞流數(shù)值模擬

把四圓柱繞正方形順排排列放置在流體中，通過改變圓柱軸線之間的距離L，分析流場(chǎng)的變化。取Re＝150，計(jì)算時(shí)間步長(zhǎng)Δi＝0．005，采用N＝7階插值函數(shù)，數(shù)值模擬L＝1．4D、2．5D、3．5D和4．0D時(shí)的流場(chǎng)。在Re＝150這樣的雷諾數(shù)下，流體無三維效應(yīng)。同樣，計(jì)算的空間精度為7階，時(shí)間上的精度為2階。

4.1 計(jì)算域和邊界條件

圖7 四圓柱排列的計(jì)算域和邊界條件示意圖Fig.7 Schematic diagram of the computational domain and boundary conditions for four cylinders

流場(chǎng)計(jì)算域和邊界條件與上述單圓柱繞流數(shù)值模擬情形相同，圖7為計(jì)算域和邊界條件示意圖。圖8（a）為兩圓柱圓心間距為2．5D時(shí)的網(wǎng)格劃分示意圖，共劃分778個(gè)網(wǎng)格，在靠近圓柱壁面的地方進(jìn)行幾層網(wǎng)格加密，離圓柱壁面最近的一層網(wǎng)格厚度為0．1D，如圖8（b）所示。同時(shí)，在圓柱尾流區(qū)域進(jìn)行加密處理。

圖8 四圓柱繞流的譜單元網(wǎng)格示意圖Fig.8 Mesh system of spectral element for cross-flow around four cylinders

4.2 不同間距比下的流態(tài)

圖9為不同間距比下的四圓柱繞流的渦量云圖。觀察發(fā)現(xiàn)三種不同的流動(dòng)形態(tài)：間距比較?。↙／D＝1．4）時(shí)，圓柱互擾以臨近效應(yīng)為主，上游兩個(gè)圓柱后面靠?jī)?nèi)側(cè)的自由剪切層附著在下游兩個(gè)圓柱表面，但是靠外側(cè)的兩個(gè)自由剪切層并未附著在下游兩個(gè)圓柱表面，而是覆蓋下游兩個(gè)圓柱，并從下游兩個(gè)圓柱后面流出，內(nèi)外兩側(cè)的自由剪切層未出現(xiàn)擺動(dòng)，圓柱尾流呈現(xiàn)為單鈍體繞流流動(dòng)形態(tài)；在中等間距比（L／D＝2．5）條件下，結(jié)構(gòu)互擾以臨近效應(yīng)和剪切層干擾為主，上游兩個(gè)圓柱后面靠?jī)?nèi)側(cè)的自由剪切層再附著在下游兩個(gè)圓柱表面，但是靠外側(cè)的兩個(gè)自由剪切層并未再附著在下游兩個(gè)圓柱表面，而是在下游兩個(gè)圓柱表面交替出現(xiàn)擺動(dòng)，此時(shí)下游圓柱尾部形成旋渦脫落流態(tài)，且干擾上游圓柱尾部旋渦發(fā)展并抵制其脫落；圓柱間距比較大（L／D＝3．5／4．0）時(shí)，結(jié)構(gòu)互擾以尾流效應(yīng)為主，此時(shí)所有圓柱尾流均充分發(fā)展，上游圓柱的自由剪切層形成成熟的渦脫落，然后撞擊到下游圓柱，且剪切層不斷向下游泄送Kármán旋渦，圓柱互擾表現(xiàn)為相鄰剪切層之間的相互干擾。由上可見，流動(dòng)形態(tài)明顯影響圓柱表面的壓力分布，尤其是當(dāng)從一種流動(dòng)形態(tài)轉(zhuǎn)變到另一種流動(dòng)形態(tài)，壓力分布變化顯著。

圖9 不同間距比下的渦量云圖Fig.9 Vorticity contours with different spacing ratios for the four cylinders

4.3 圓柱表面受力特性隨間距比的變化

圖10為四圓柱繞流情形下的柱體平均阻力系數(shù)隨間距比的變化，同時(shí)與相關(guān)文獻(xiàn)［11-12］在不同雷諾數(shù)下的結(jié)果對(duì)比。從圖10中看出，圓柱1和2的平均阻力系數(shù)CD1和CD2基本上相等，圓柱3和4CD3和CD4也呈現(xiàn)相同的特征，這是對(duì)稱性造成的結(jié)果，但是在間距比較小時(shí)，則有細(xì)微差異。另外，比較不同雷諾數(shù)的結(jié)果，發(fā)現(xiàn)雷諾數(shù)對(duì)阻力系數(shù)的影響較大，但是數(shù)值在同一數(shù)量級(jí)范圍內(nèi)。

圖10 四圓柱平均阻力系數(shù)隨間距比的變化Fig.10 Variation of mean drag coefficients with different spacing ratios for the four cylinders

圖11 四圓柱平均升力系數(shù)隨間距比的變化Fig.11 Variation of mean lift coefficients with different spacing ratios for the four cylinders

圖11為四圓柱繞流情形下的柱體平均升力系數(shù)隨間距比的變化，同時(shí)與相關(guān)文獻(xiàn)［14-15］在不同雷諾數(shù)下的結(jié)果對(duì)比。從圖中可見，圓柱1和2平均升力系數(shù)CL1和CL2相反，圓柱3和4平均升力系數(shù)CL3和CL4也相反，這反映出兩對(duì)圓柱間互相排斥。隨著間距比的增大，4個(gè)圓柱的平均升力系數(shù)的絕對(duì)值均減小，且趨于0。上游兩個(gè)圓柱的平均升力系數(shù)絕對(duì)值大于下游兩個(gè)圓柱。在層流范圍內(nèi)不同雷諾數(shù)的情況下，柱體的平均升力系數(shù)接近。但是在湍流情況下，間距越小，平均升力系數(shù)會(huì)更大。

4.4 斯托羅哈數(shù)隨間距比的變化

圖12為四圓柱繞流情形下的斯托羅哈數(shù)St1-St4隨間距比的變化。因本文變量均經(jīng)無量綱化處理，所以斯托羅哈數(shù)即等于旋渦脫落頻率。旋渦脫落頻率可通過對(duì)升力曲線作傅立葉變換（FFT）分析求得。從圖9中看出，四個(gè)圓柱旋渦脫落頻率近乎相同。從L／D＝1．4到L／D＝2．5是從一種尾流形態(tài)轉(zhuǎn)變成另一種完全不同尾流形態(tài)，旋渦脫落頻率有所減小，從L／D＝2．5到L／D＝3．5尾流形態(tài)再發(fā)生轉(zhuǎn)變，旋渦脫落頻率增大。L／D＝1．4時(shí)，尾流形態(tài)類似于單鈍體繞流時(shí)的尾流，尾流只有一個(gè)占主要地位的旋渦脫落；L／D＝2．5時(shí)，上下側(cè)兩對(duì)圓柱尾流中分別形成旋渦脫落，類似于兩對(duì)反對(duì)稱的串列圓柱，但是每一對(duì)串列圓柱后面的尾流又不同于兩串列圓柱時(shí)的尾流，尾流中有一個(gè)占主要地位的旋渦，還有一個(gè)小渦伴隨，而并非兩排相同大小的旋渦；之后隨著間距比繼續(xù)增大，旋渦脫落頻率隨之增大，并逐漸趨近于單圓柱繞流時(shí)的旋渦脫落頻率，上下兩側(cè)尾流形態(tài)也類似于兩對(duì)串列圓柱尾流形態(tài)。

圖12 四圓柱斯托羅哈數(shù)隨間距比的變化Fig.12 Variation of the Strouhal numbers with different spacing ratios for the four cylinders

比較雷諾數(shù)的影響，發(fā)現(xiàn)從層流到湍流范圍內(nèi)，旋渦脫落頻率在0．05～0．2之間，在同一數(shù)量級(jí)。

5 結(jié) 論

本文系統(tǒng)闡明了譜單元方法的基本原理，基于譜單元方法對(duì)低雷諾數(shù)和不同間距比下的順排雙圓柱繞流和正方形排列的四圓柱擾流及其阻力系數(shù)、升力系數(shù)等關(guān)鍵參數(shù)進(jìn)行數(shù)值模擬和比較分析。

對(duì)順排雙圓柱擾流Re＝200，間距比選擇包括1．2、1．5、2．0、3．6、4．0、6．0、8．0和10．0等八種。分析顯示，間距比對(duì)兩圓柱繞流流場(chǎng)特性影響明顯，當(dāng)流場(chǎng)從一種流動(dòng)形態(tài)變成另外一種流動(dòng)形態(tài)時(shí)，關(guān)鍵力學(xué)參數(shù)變化顯著，在某些間距比區(qū)間內(nèi)出現(xiàn)驟升或驟降的現(xiàn)象。順排雙圓柱繞流存在臨界間距比，在Re＝200時(shí)臨界間距比約為3．6；當(dāng)間距比較小時(shí)（L／D≤3．6），上游圓柱后不發(fā)生渦脫落，而僅在下游圓柱后發(fā)生渦脫落；間距比增大時(shí)（L／D≥4．0），兩圓柱之間產(chǎn)生渦脫落，且上游圓柱后的渦脫落對(duì)下游圓柱產(chǎn)生較大的影響。在臨界間距比區(qū)間（3．6≤L／D≤4．0），平均阻力系數(shù)、平均阻力系數(shù)均方根、升力系數(shù)均方根和斯托羅哈數(shù)變化顯著。

對(duì)于正方形排列四圓柱繞流問題在Re＝150，發(fā)現(xiàn)了三種不同的流動(dòng)形態(tài)：間距比較小（L／D＝1．4）時(shí)，圓柱互擾以臨近效應(yīng)為主，圓柱尾流呈現(xiàn)為單鈍體繞流流動(dòng)形態(tài)；中等間距比（L／D＝2．5）時(shí)，結(jié)構(gòu)互擾以臨近效應(yīng)和剪切層干擾為主，此時(shí)下游圓柱尾部形成旋渦脫落流態(tài)，且干擾上游圓柱尾部旋渦發(fā)展并抑制其脫落；圓柱間距較大（L／D＝3．5／4．0）時(shí)，結(jié)構(gòu)互擾以尾流效應(yīng)為主，此時(shí)所有圓柱尾流均充分發(fā)展，且剪切層不斷向下游泄送Kármán旋渦，圓柱互擾表現(xiàn)為相鄰剪切層之間的相互干擾。圓柱1和圓柱2的平均阻力系數(shù)在L／D＝2．5處發(fā)生跳躍，四個(gè)圓柱的斯托羅哈數(shù)均在L／D＝2．5處有一個(gè)彎折，之后隨著間距比的增大而增大。

［1］GOTTLIEB D，ORGSZAG S A．Numerical analysis of spectral method：theory and application［A］．CB Ms-NFS Monograph No 26．Society for Industry and Applied Mathematics［C］．Philadelphia，1977．

［2］GOTTLIEB D，HUSSAINI M Y，ORGSZAG S A．Theory and application of spectral method［A］．In：Voigt RG．Gottlieb D．Hussaini M Y．Spectral Method for Partial Differential Equation［C］．SIAM Philadelphia，1984：1-54．

［3］PATERA A T．A spectral element method for fluid dynamics：laminar flow in a channel expansion［J］．Journal of computattonal Phystcs，1984，154：468-488．

［4］TAN Lingyan．Numerical simulation of flow past a main cylinder with an affiliated cylinder［J］．Journal of Jtltn Untverstty（Sctence Edttton），2012，50（1）：69-72．（in Chinese）

譚玲燕．?dāng)?shù)值模擬放置附屬圓柱的主圓柱繞流［J］．吉林大學(xué)學(xué)報(bào)（理學(xué)版），2012，50（1）：69-72．

［5］WANG Tong，CAO Shuyang，ZHOU Qiang．Investigation of aerodynamic lift of a circular cylinder in linear shear flows at subcritical Reynolds number［J］．Chtnese Quarterly of Mechantcs，2011，32（4）：473-479．（in Chinese）

王通，曹曙陽，周強(qiáng)．亞臨界雷諾數(shù)線性剪切流場(chǎng)中圓柱的氣動(dòng)升力特性研究［J］．力學(xué)季刊，2011，32（4）：473-479．

［6］GONG Shuai，GUO Zhaoli．Lattice Boltzmann simulation of flow over a transversly oscillation circular cylinder［J］．Chtnese Journal of Theorettcal and Applted Mechantcs，2011，43（5）：809-818．（in Chinese）

龔帥，郭照立．橫向振蕩圓柱繞流的格子Boltzmann方法模擬［J］．力學(xué)學(xué)報(bào)，2011，43（5）：809-818．

［7］BLACKBURN H M，SHERWIN S J．Formulation of a Galerkin spectral element-Fourier method for three-dimensional incompressible flows in cylindrical geometries［J］．Journal of Computattonal Phystcs，2004，197：759-778．

［8］BLACKBURN H M，LOPEZ J M．Modulated waves in a periodically driven annular cavity［J］．Journal of Flutd Mechantcs，2011，667：336-357．

［9］BLACKBURN H M，SHEARD G J．On quasi-periodic and subharmonic Floquet wake instabilities［J］．Phystcs of Flutds，2010，22：031701-1-4．

［10］BLACKBURN H M，BARKLEY D，SHERWIN S J．Convective instability and transient growth in flow over a backward-facing step［J］．Journal of Flutd Mechantcs，2008，603：271-304．

［11］FRANKE R，RODI W，SCHONUNG B．Numerical calculation of laminar vortex shedding flow past cylinders［J］．Journal of Wtnd Engtneertng and Industrtal Aerodynamtcs，1990，35：237-257．

［12］MENEGHINI J R，SALTARA F，SIQUEIRA C L R，et al．Numerical simulation of flow interference between two circular cylinders in tandem and side-by-side arrangements［J］．Journal of Flutds and Structures，2001，15：327-350．

［13］FARRANT T，TAN M，PRICE W G．A cell boundary element method applied to laminar vortex-shedding from arrays of cylinders in various arrangements［J］．Journal of Flutds and Structures，2000，14：375-402．

［14］LAM K，LI J Y，SO R M C．Force coefficient and strouhal numbers of four cylinders in cross flow［J］．Journal of Flutds and Structures，2003，18：305-324．

［15］LAM K，GONG W Q，SO R M C．Numerical simulation of cross-flow around four cylinders in an in-line square configuration［J］．Journal of Flutds and Structures，2008，24：34-57．

Numerical simulation of cross-flow around multiple circular cylinders by spectral element method

HAN Zhaolong，ZHOU Dai，CHEN Yanan，GUI Xiaolan，LI Junlong
（School of Naval architecture，Ocean and Civil Engineering，Shanghai Jiao Tong University，Shanghai 200240，China）

The details of the spectral element method are presented．Flow past two tandem circular cylinders at the Reynolds numberRe＝200 and four square-arranged circular cylinders with various spacing ratios atRe＝150 are numerically investigated．The effects of spacing ratiosL／D（ratio of center-to-center distance to diameter of cylinder）for the two cases above are main concerned．The relationship between the spacing ratio and the flow patterns，mean drag and lift coefficients as well as the Strouhal numbers are studied．The results reveal that the spacing ratio has a great impact on the flow characteristics of both the two-cylinder andthe four-cylinder cases．There is a critical spacing ratio of around 3．6 atRe＝200 for the case of two tandem cylinders．Three different flow patterns are observed in the case of four-cylinders．Among different flow patterns，a transformation from one pattern to another will change the drag／lift forces and vortex shedding frequencies．The key mechanical parameters jump suddenly in some critical spacing ratio intervals．

spectral element method；two tandem cylinders；four square-arranged cylinders；spacing ratio；flow patterns

O353.4

A doi:10.7638／kqdlxxb-2012.0070

0258-1825（2014）01-0021-10

2012-08-07；

2012-11-22

國(guó)家自然科學(xué)基金（51078230，11172174，51278297）、上海市優(yōu)秀學(xué)術(shù)帶頭人計(jì)劃（13XD1402100）項(xiàng)目資助

韓兆龍（1983-），男，工學(xué)博士，專業(yè)：結(jié)構(gòu)工程，研究方向：風(fēng)工程與流固耦合．E-mail：han．a(chǎn)rkey＠gmail．com

周岱，男（1963-），工學(xué)博士，教授，研究方向：大跨空間結(jié)構(gòu)與流固耦合．E-mail：zhoudai＠sjtu．edu．cn

韓兆龍，周 岱，陳亞楠，等．譜單元法及其在多圓柱繞流分析中的應(yīng)用［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32（1）：21-29．

10．7638／kqdlxxb-2012．0070．HAN Z L，ZHOU D，CHEN Y N，et al．Numerical simulation of cross-flow around multiple circular cylinders by spectral element method［J］．ACTA Aerodynamica Sinica，2014，32（1）：21-29．