適用于渦激振蕩問(wèn)題研究的并行高精度方法

邱滋華，徐敏, *，張斌，梁春雷

1. 西北工業(yè)大學(xué) 航天學(xué)院，西安 710072 2. 喬治華盛頓大學(xué) 機(jī)械與航空工程系，華盛頓特區(qū) 20052

因材料特性和流體作用力形式的不同，流體流過(guò)固體結(jié)構(gòu)時(shí)可以引發(fā)多種形式的振動(dòng)[1]。例如，橋梁或機(jī)翼的顫振和馳振等不穩(wěn)定振動(dòng)；高樓或電線的抖振和渦激振蕩(Vortex-Induced Vibration, VIV)等共振振動(dòng)；多邊形截面和直升機(jī)旋翼的自旋轉(zhuǎn)動(dòng)。在上述諸多形式中，鈍體的渦脫落及其VIV現(xiàn)象由于具有廣泛的工程應(yīng)用，且較為基礎(chǔ)和易于實(shí)驗(yàn)，近年來(lái)持續(xù)受到廣泛關(guān)注[2-5]。VIV常發(fā)生在鈍體繞流伴隨渦脫落的情況下，不穩(wěn)定的渦施加振蕩力于物體表面，從而引起振動(dòng)。當(dāng)渦脫落的頻率接近結(jié)構(gòu)的自然頻率時(shí)，會(huì)發(fā)生共振，振動(dòng)的幅值也會(huì)大大增加。

Feng在1968年進(jìn)行的實(shí)驗(yàn)被認(rèn)為是第一批現(xiàn)代VIV研究之一[6]，在其研究中，圓柱只能垂直于來(lái)流振動(dòng)。Khalak和Williamson在實(shí)驗(yàn)研究了質(zhì)量和阻尼對(duì)單個(gè)圓柱振動(dòng)的影響后，發(fā)現(xiàn)相比于經(jīng)典的Feng式響應(yīng)，參數(shù)的改變會(huì)使得響應(yīng)模式更為豐富[7-8]。Mittal和Kumar[9]、Jeon和Gharib[10]、Jauvtis和Williamson[11]、Prasanth等[12]研究了圓柱的二自由度渦激振蕩，發(fā)現(xiàn)沿來(lái)流方向的振動(dòng)幅值比垂直來(lái)流方向的振動(dòng)幅值小若干個(gè)量級(jí)。Yao和Jaiman[13-14]提出了適用于渦激振蕩問(wèn)題研究的降階模型，并研究了針對(duì)VIV的反饋控制算法。以上研究發(fā)現(xiàn)，VIV對(duì)于系統(tǒng)的物理參數(shù)非常敏感，表現(xiàn)出了極強(qiáng)的非線性，其響應(yīng)機(jī)理還有待進(jìn)一步研究。

VIV問(wèn)題的數(shù)值模擬主要存在以下幾個(gè)難點(diǎn)。首先，VIV對(duì)于系統(tǒng)的物理參數(shù)非常敏感，較小的參數(shù)變化也有可能引起響應(yīng)的截然不同[12]，因此為了準(zhǔn)確地模擬系統(tǒng)響應(yīng)，需要盡可能精確地刻畫(huà)流場(chǎng)細(xì)節(jié)，比如渦，得到流體作用于固體上的力，從而對(duì)求解器的高精度以及方法的低耗散有要求。其次，為了適應(yīng)結(jié)構(gòu)的振動(dòng)，通常需要網(wǎng)格變形。但當(dāng)結(jié)構(gòu)振動(dòng)的幅值很大、或者多個(gè)振動(dòng)物體之間間距很小時(shí)，變形網(wǎng)格通常會(huì)很扭曲、甚至出現(xiàn)負(fù)體積，從而引入過(guò)大的幾何結(jié)構(gòu)誤差，或者直接使得仿真發(fā)散。

相對(duì)于低階精度方法而言，高精度方法(三階及以上)只引入很小的數(shù)值耗散。譜方法是最經(jīng)典的高精度方法之一，在相同的條件下具有最高的精度[15]。但由于譜方法固有的全域性，其應(yīng)用范圍受到了很大限制，難以計(jì)算復(fù)雜區(qū)域的流場(chǎng)。國(guó)內(nèi)學(xué)者在譜方法的推進(jìn)與改良中做出了很大貢獻(xiàn)[16-17]?；诰W(wǎng)格單元連續(xù)的譜差分(Spectral Difference, SD)方法是一種高效的非結(jié)構(gòu)網(wǎng)格高精度算法[18-24]，本文引入了SD方法在三角形/四邊形混合網(wǎng)格上對(duì)Navier-Stokes方程直接進(jìn)行空間離散，對(duì)復(fù)雜流場(chǎng)的適應(yīng)性非常好?；赟D方法，Zhang等[25-26]提出了滑移網(wǎng)格譜差分(Sliding-mesh Spectral Difference, SSD)方法，用于求解自由轉(zhuǎn)動(dòng)物體周圍的流場(chǎng)，相比于傳統(tǒng)的網(wǎng)格協(xié)調(diào)技術(shù)，該方法通過(guò)引入滑移界面可以大大減小網(wǎng)格的扭曲程度，從而實(shí)現(xiàn)高精度和高效率計(jì)算。本文將SSD方法進(jìn)一步推廣，使其可以處理非均勻滑移網(wǎng)格界面，從而減小VIV仿真時(shí)網(wǎng)格的扭曲程度，特別是大幅振動(dòng)或者多個(gè)物體振動(dòng)的情況。最后，研究了并行策略，對(duì)求解器進(jìn)行消息傳遞接口(Message Passing Interface, MPI)并行，大大提高了仿真效率，增加了求解器的實(shí)用性。

1 控制方程

1.1 流體方程

對(duì)于流體而言，考慮如下守恒形式的二維非定常Navier-Stokes方程

(1)

式中：Qf為守恒變量；Ff和Gf分別為沿x和y方向的通量，下標(biāo)f代表流體。為了模擬動(dòng)網(wǎng)格上的流動(dòng)，采用任意拉格朗日-歐拉(Arbitrary-Lagrangian-Eulerian, ALE)算法。在此方法中，將真實(shí)的時(shí)間和空間(t,x(t),y(t))投影到計(jì)算空間上(τ,ξ,η)，其中(ξ,η)表示時(shí)間無(wú)關(guān)的計(jì)算空間，時(shí)間τ=t。由此，Navier-Stokes方程在計(jì)算空間中可表示為

(2)

式中：

(3)

其中:|J|為變換雅克比矩陣J的秩；J-1為J的逆矩陣。

1.2 固體振動(dòng)方程

研究鈍體VIV問(wèn)題所廣泛采用的模型是，將一個(gè)圓柱彈性連接，置于來(lái)流之中，如圖1所示。若將振動(dòng)限制在縱向(比如y向)，則固體的振動(dòng)就由如下阻尼諧振子方程控制

(4)

式中：m和y分別為固體的質(zhì)量和縱向位移；c為阻尼系數(shù)；k為彈簧的剛度系數(shù)；FL為流體作用在固體上的升力(縱向力)。

圖1 圓柱的質(zhì)量彈簧阻尼系統(tǒng)示意圖Fig.1 Diagram of a mass-spring-damper system for a circular cylinder

1.3 流固耦合方程

為了耦合固體方程和流體方程(對(duì)時(shí)間具有一階導(dǎo)數(shù))，需要將式(4)變換為兩個(gè)時(shí)間一階導(dǎo)的方程組

(5)

式(5)進(jìn)一步可以寫(xiě)成如下殘差形式

(6)

式中：

(7)

(8)

Qs和Rs分別為固體的自變量和殘差的向量,下標(biāo)s代表固體。

類似地，將流體方程也寫(xiě)成如下的殘差形式

(9)

式中：

(10)

(11)

聯(lián)立式(6)和式(9)可得到全耦合的守恒系統(tǒng)

(12)

式中：

(13)

Q和R(Q)分別為耦合的自變量和殘差的向量。VIV數(shù)值模擬中，將對(duì)式(12)進(jìn)行時(shí)間推進(jìn)。

2 數(shù)值方法

2.1 高階網(wǎng)格下的SD方法

對(duì)流體方程進(jìn)行空間離散首先要通過(guò)式(14)的映射方程，將物理空間內(nèi)的動(dòng)網(wǎng)格單元映射到計(jì)算空間的靜止網(wǎng)格單元上

(14)

式中：(x,y)表示物理空間；(ξ,η)表示計(jì)算空間；K為一個(gè)網(wǎng)格單元的節(jié)點(diǎn)數(shù)；(xi,yi)為隨時(shí)間變化的第i個(gè)節(jié)點(diǎn)坐標(biāo)；Mi為和節(jié)點(diǎn)對(duì)應(yīng)的形函數(shù)。

圖2顯示了用不同精度的網(wǎng)格單元去擬合真實(shí)物理單元的情況。比如，對(duì)于一個(gè)二階三角形和四邊形網(wǎng)格單元，K=6,8。顯然，相比于低精度網(wǎng)格，高精度曲邊網(wǎng)格能夠把具有曲面(邊)幾何形狀的物體擬合得更精確?；谶@個(gè)原因，本文中所有曲線邊界都用三階精度網(wǎng)格進(jìn)行擬合。

SD方法通過(guò)在網(wǎng)格單元內(nèi)部分布的解點(diǎn)(Solution Points，SPs)與通量點(diǎn)(Flux Points，F(xiàn)Ps)來(lái)構(gòu)造單元內(nèi)連續(xù)的解多項(xiàng)式和通量多項(xiàng)式。關(guān)于SD方法在三角形/四邊形混合網(wǎng)格上解點(diǎn)與通量點(diǎn)分布的詳細(xì)描述，可參考文獻(xiàn)[27]。

圖2 不同階數(shù)的網(wǎng)格單元(細(xì)黑線)與真實(shí)的物理單元(粗灰線)Fig.2 Grid elements (thin black lines) in different orders for approximating curved physical cells (thick gray lines)

對(duì)于三階精度SD格式而言，三角形和四邊形單元上的解點(diǎn)與通量點(diǎn)分布如圖3所示，灰色箭頭表示通量點(diǎn)處的自由度方向。

由于SD方法所構(gòu)造的解多項(xiàng)式和通量多項(xiàng)式只在網(wǎng)格單元內(nèi)部連續(xù)，在網(wǎng)格單元邊界處是間斷的。為了保證方法的守恒以及穩(wěn)定，本文采用Rusanov黎曼求解器[28]在單元邊界處計(jì)算相鄰網(wǎng)格之間的共同無(wú)黏通量。

圖3 三階SD格式四邊形和三角形單元上的解點(diǎn)(藍(lán)色圓圈)和通量點(diǎn)(橘色方塊)分布示意圖Fig.3 Schematic of distribution of solution points (blue circles) and flux points (orange squares) for a third-order SD scheme on a quadrilateral element and a triangular element

2.2 動(dòng)網(wǎng)格策略

對(duì)于VIV而言，當(dāng)固體結(jié)構(gòu)位置變化時(shí)，流場(chǎng)區(qū)域也隨之改變。本文將調(diào)配函數(shù)[29]與滑移網(wǎng)格相結(jié)合，使之適用于極端情況的網(wǎng)格變形問(wèn)題。

以圖4中單個(gè)縱向振動(dòng)圓柱的例子，詳細(xì)解釋一下該方法。通過(guò)兩個(gè)滑移界面，將整個(gè)流場(chǎng)區(qū)域分為3個(gè)子區(qū)域。兩邊的兩個(gè)子區(qū)域是靜止的，中間的區(qū)域可以變形以適應(yīng)圓柱的運(yùn)動(dòng)。中間的區(qū)域進(jìn)一步以距圓心d1和d2為界,分為3個(gè)虛擬的區(qū)域：在區(qū)域Ⅰ中，網(wǎng)格是剛性的，隨著圓柱一起運(yùn)動(dòng)；區(qū)域Ⅱ中，網(wǎng)格可以變形；區(qū)域Ⅲ中，網(wǎng)格保持靜止。

圖5用兩個(gè)振動(dòng)圓柱網(wǎng)格變形的例子來(lái)說(shuō)明滑移網(wǎng)格方法的優(yōu)勢(shì)。圖5(a)是初始網(wǎng)格，兩個(gè)圓柱相隔兩倍直徑放置。圖5(b)是當(dāng)兩個(gè)圓柱在縱向相對(duì)運(yùn)動(dòng)了一個(gè)直徑的距離時(shí)，用連續(xù)協(xié)調(diào)變形方法得到的網(wǎng)格。顯然，網(wǎng)格在第1個(gè)圓柱的上方和第2個(gè)圓柱的下方變得非常扭曲。圖5(c)顯示了滑移網(wǎng)格方法的結(jié)果，充分說(shuō)明了通過(guò)引入滑移界面，可以很大程度地提高網(wǎng)格質(zhì)量。

圖4 用兩個(gè)滑移界面將單個(gè)振動(dòng)圓柱計(jì)算域分割的示意圖Fig.4 Schematic of a computational domain with two sliding interfaces for a vibrating cylinder

圖5 圍繞兩個(gè)振動(dòng)圓柱的網(wǎng)格變形示意圖Fig.5 Schematic of meshes deformation for two vibrating cylinders

理想情況下，網(wǎng)格運(yùn)動(dòng)不應(yīng)該污染流場(chǎng)。最簡(jiǎn)單的情況便是：動(dòng)網(wǎng)格下定常自由來(lái)流應(yīng)能始終保持定常，稱之為自由流保持(Free-stream Preservation)。為了滿足自由流保持，將自由來(lái)流參數(shù)代入式(2)，經(jīng)過(guò)變換，可以得到如下方程組：

(15)

式中:ξx、ξy、ξt、ηx、ηy和ηt均為變換雅克比矩陣J的元素。

式(15)只與網(wǎng)格的幾何變量相關(guān)，與流場(chǎng)參數(shù)無(wú)關(guān)，所以通常稱之為幾何守恒律(Geometric Conservation Law, GCL)[30]。對(duì)式(15)方程組進(jìn)行積分處理后會(huì)發(fā)現(xiàn)：前2個(gè)方程的物理意義是一個(gè)閉合的單元且必須一直保持閉合；第3個(gè)方程的物理意義是物理空間的體積變化率等于其邊界擴(kuò)展的速度。

對(duì)于GCL的處理，國(guó)內(nèi)外學(xué)者已經(jīng)做了很多研究[31-35]。由映射關(guān)系式(14)可知，和空間變換相關(guān)的量(|J|ξx等)可以表示成拉格朗日插值多項(xiàng)式形式，而SD格式中空間離散是對(duì)變量直接進(jìn)行微分，因此，若構(gòu)造曲線網(wǎng)格的多項(xiàng)式階數(shù)不高于構(gòu)造守恒變量的多項(xiàng)式階數(shù)，式(15)中的前兩個(gè)等式能夠自動(dòng)滿足[31,35]。但SD格式中時(shí)間推進(jìn)并不是解析的，導(dǎo)致式(15)中第3個(gè)方程并不能自動(dòng)滿足。一種可行的做法是將|J|視為未知量，采用和流場(chǎng)方程相同的時(shí)間及空間離散格式進(jìn)行更新，取代由網(wǎng)格坐標(biāo)直接計(jì)算得到|J|的做法[33]。

在SD方法中的具體計(jì)算過(guò)程為：將GCL第3個(gè)方程寫(xiě)成如下形式:

(16)

2.3 滑移網(wǎng)格分塊

并行計(jì)算首先要將全部網(wǎng)格分配給各個(gè)處理器，由于計(jì)算網(wǎng)格已經(jīng)被滑移界面分成了幾個(gè)子區(qū)域，因此并行時(shí)網(wǎng)格的分塊策略值得詳細(xì)研究。假設(shè)一共有NP個(gè)處理器，滑移界面將網(wǎng)格分為NM個(gè)子區(qū)域，將子區(qū)域合起來(lái)視為一個(gè)大的區(qū)域，在這個(gè)單獨(dú)但間斷的區(qū)域上進(jìn)行網(wǎng)格分塊。步驟如下：

步驟1根處理器讀入NM個(gè)子區(qū)域的網(wǎng)格信息并存儲(chǔ)。

步驟2通過(guò)加上相應(yīng)的偏移量，根處理器將各個(gè)子區(qū)域上的網(wǎng)格單元、節(jié)點(diǎn)以及邊界等的編號(hào)轉(zhuǎn)成全局編號(hào)。偏移量是之前子區(qū)域上相應(yīng)項(xiàng)目數(shù)量之和。例如，對(duì)于一個(gè)具有3個(gè)子區(qū)域的問(wèn)題而言，若每個(gè)子區(qū)域上的網(wǎng)格單元數(shù)量都是100，那么子區(qū)域1、2和3的網(wǎng)格單元偏移量分別是0、100和200。

步驟3根處理器調(diào)用Metis庫(kù)[36]里的函數(shù)，將整個(gè)大區(qū)域上的網(wǎng)格分成NP塊。

步驟4根處理器將分塊信息播散給其他子處理器。

步驟5所有處理器根據(jù)分塊文件讀入各自塊的網(wǎng)格信息。

這種方法的優(yōu)勢(shì)是其實(shí)際的分塊數(shù)量通常是最小的，減少了處理器之間的信息傳遞。另外，可以實(shí)現(xiàn)全局的負(fù)載平衡。比如圖6是6個(gè)處理器時(shí)圖5(c)對(duì)應(yīng)網(wǎng)格的分塊情況，其中粗黑線表示分塊邊界。從圖中可以看出，P0和P5的分塊都不連續(xù)，分別在兩個(gè)子區(qū)域上，但對(duì)于全局而言實(shí)現(xiàn)了負(fù)載平衡。實(shí)際上，對(duì)于網(wǎng)格量很大的情況(相對(duì)于處理器數(shù)量而言)，分區(qū)通常都是連續(xù)的，因此這種方法在達(dá)到負(fù)載平衡的同時(shí)也能使分塊數(shù)量最小。

圖6 6個(gè)處理器對(duì)繞兩個(gè)振動(dòng)圓柱網(wǎng)格的分塊示意圖Fig.6 Schematic of partitions of a mesh around two vibrating cylinders via 6 processors

2.4 滑移界面通信

滑移界面兩側(cè)的網(wǎng)格不均勻而且網(wǎng)格點(diǎn)有錯(cuò)位，導(dǎo)致了兩側(cè)網(wǎng)格實(shí)際上是間斷的，本文采用粘接元方法進(jìn)行滑移界面兩側(cè)的信息傳遞[37]。粘接元方法的主要思想是：將解/通量從網(wǎng)格面Ω投影到粘接元Ξ(如圖7(a))，計(jì)算得粘接元上的解/通量后，將數(shù)據(jù)投影回網(wǎng)格面(如圖7(b))，以保證流場(chǎng)的守恒。滑移界面處，每一個(gè)網(wǎng)格面對(duì)應(yīng)一個(gè)或者多個(gè)粘接元，而每一個(gè)粘接元左右各對(duì)應(yīng)一個(gè)網(wǎng)格面。在時(shí)間推進(jìn)的每一步都需要更新這種網(wǎng)格面-粘接元的對(duì)應(yīng)關(guān)系。

圖8(a)為串行求解器中兩塊滑移網(wǎng)格之間的粘接元分布示意圖，影線部分表示粘接元。對(duì)于并行求解器而言，情況有所差別。圖8(b)顯示了5個(gè)處理器時(shí)兩塊滑移網(wǎng)格之間的粘接元分布，其中黑色影線部分表示當(dāng)?shù)卣辰釉?，灰色陰影部分表示全局粘接元。每個(gè)處理器都對(duì)應(yīng)有當(dāng)?shù)卣辰釉颐總€(gè)粘接元只和一個(gè)網(wǎng)格面相對(duì)應(yīng)。如圖8(b)中處理器P1上的當(dāng)?shù)卣辰釉缓推渥髠?cè)的網(wǎng)格面對(duì)應(yīng)，因此左側(cè)的解/通量可以從當(dāng)?shù)刂苯荧@得，而右側(cè)的解/通量則需要從處理器P3或P4其上對(duì)應(yīng)粘接元發(fā)出。對(duì)于其他處理器上的粘接元也類似。通過(guò)當(dāng)?shù)?全局粘接元對(duì)應(yīng)關(guān)系，當(dāng)?shù)卣辰釉梢哉业狡湓谄渌幚砥魃蠈?duì)應(yīng)的粘接元。

圖7 網(wǎng)格面和粘接元之間的投影關(guān)系Fig.7 Projection between a cell face and mortars

圖8 兩塊滑移網(wǎng)格之間的粘接元分布示意圖Fig.8 Schematic of distribution of mortars between two sub-domain meshes with relative sliding motion

為了保證處理器之間信息傳遞的正確，網(wǎng)格面-粘接元對(duì)應(yīng)關(guān)系以及當(dāng)?shù)?全局粘接元對(duì)應(yīng)關(guān)系需要在每一迭代步更新。由于滑移界面處的網(wǎng)格節(jié)點(diǎn)數(shù)是固定的、在仿真過(guò)程中總粘接元數(shù)量也保持不變，因此，為了加快速度，可以用兩個(gè)處理器分別從同一個(gè)滑移界面的兩端開(kāi)始更新網(wǎng)格面-粘接元以及當(dāng)?shù)?全局粘接元對(duì)應(yīng)關(guān)系。

3 數(shù)值算例

本節(jié)驗(yàn)證非均勻滑移網(wǎng)格方法對(duì)無(wú)黏和有黏流動(dòng)的空間精度，并證明高階精度格式的低耗散特性，隨之，模擬4個(gè)圓柱的VIV算例來(lái)驗(yàn)證并行求解器的加速效果。算例采用5步4階顯式龍格-庫(kù)塔格式進(jìn)行時(shí)間推進(jìn)[38]。所有算例都是在無(wú)量綱化后進(jìn)行計(jì)算的。

3.1 Euler vortex流動(dòng)

在Euler vortex流動(dòng)[39]中，均勻來(lái)流輸運(yùn)一個(gè)等熵渦。Euler vortex流動(dòng)在無(wú)窮大區(qū)域內(nèi)的解析解為

(17)

式中：u、v、ρ和p分別為Euler vortex流動(dòng)中x和y方向的速度、密度和壓力；U∞、ρ∞、p∞和Ma∞分別為均勻流的速度、密度、壓力和馬赫數(shù)；γ為比熱比，通常取為1.4；θ為均勻流方向；ε和rc分別為渦的強(qiáng)度和尺寸；相對(duì)坐標(biāo)(xr,yr)的定義為

(18)

在本仿真中，均勻來(lái)流的流動(dòng)參數(shù)是：(U∞,ρ∞)=(1,1)，馬赫數(shù)Ma∞=0.3，渦隨均勻流沿θ=π/6方向運(yùn)動(dòng)。計(jì)算域?yàn)?≤x≤10,0≤y≤10。仿真初始時(shí)刻，將強(qiáng)度和尺寸分別為ε=1和rc=1的渦置于計(jì)算域中心，即(x0,y0)=(5,5)。在x=3和x=7兩處將計(jì)算域分為3部分。兩側(cè)的子區(qū)域用三角形單元?jiǎng)澐?，中間用四邊形單元?jiǎng)澐帧Ｓ镁W(wǎng)格量分別為44/15、156/50、568/200和2 216/800的4套三角形/四邊形混合網(wǎng)格進(jìn)行精度測(cè)試。令中間子區(qū)域的水平中線(y=5，3≤x≤7)以方程Δy=sin(0.5t)隨時(shí)間運(yùn)動(dòng)，變形區(qū)域限定為|y-5|≤4。兩側(cè)子區(qū)域網(wǎng)格則保持靜止。

圖9顯示了在568/200混合網(wǎng)格上用8個(gè)處理器并行求解四階精度格式得到的流場(chǎng)密度云圖，圖9(a)對(duì)應(yīng)初始時(shí)刻，圖9(b)為無(wú)量綱時(shí)間t=2.3時(shí)刻渦心剛好運(yùn)動(dòng)到滑移界面時(shí)的情況，圖9中淺灰色為網(wǎng)格線，深灰色為處理器邊界?？梢钥闯?，渦的層次清晰，滑移界面完全不會(huì)污染流場(chǎng)。

通過(guò)密度的L1和L2誤差，進(jìn)一步驗(yàn)證了計(jì)算精度。兩個(gè)誤差定義為

(19)

(20)

(21)

圖9 Euler vortex流動(dòng)對(duì)應(yīng)的密度云圖Fig.9 Density contours for Euler vortex flow

表1給出了t=2.3時(shí)刻上述4套網(wǎng)格在三階和四階格式下的誤差和精度。精度的計(jì)算依據(jù)為

(22)

式中：order為精度階；error為誤差值；下標(biāo)coarse和fine分別表示粗網(wǎng)格和密網(wǎng)格。

從表1中可見(jiàn)，誤差隨著網(wǎng)格的加密逐漸減小。當(dāng)使用三階精度格式時(shí)，數(shù)值計(jì)算得到的精度階數(shù)略小于3；四階精度格式得到的結(jié)果略大于4。很明顯，盡管存在不連續(xù)的滑移界面，滑移網(wǎng)格方法還是能夠在無(wú)黏流動(dòng)中保持SD方法的高精度。

表1Eulervortex流動(dòng)仿真的誤差和精度的階數(shù)(基于密度計(jì)算)

Table1Errorsandordersofaccuracy(basedondensity)forEulervortexflowsimulation

格式單元數(shù)L1誤差L1誤差精度階L2誤差L2誤差精度階三階44/151.58×10-32.97×10-3156/503.08×10-42.625.69×10-42.65568/2004.32×10-52.979.17×10-52.762216/8006.58×10-62.751.47×10-52.67四階44/156.95×10-41.37×10-3156/505.76×10-54.001.14×10-44.00568/2003.44×10-64.257.43×10-64.122216/8002.00×10-74.154.34×10-74.14

3.2 平板Couette流動(dòng)

對(duì)于在兩個(gè)相距H的無(wú)限長(zhǎng)平板之間的黏性平板Couette流，如果上板的溫度是Tt且以速度U運(yùn)動(dòng)，下板溫度為T(mén)b且靜止，那么穩(wěn)態(tài)流場(chǎng)為

(23)

式中：e=cVT為內(nèi)能,cV=R/(γ-1)為定容比熱容,T為溫度，R為氣體常數(shù)；下標(biāo)b和t分別表示下板和上板；Pr為普朗特?cái)?shù)，本文取為0.72。

在本算例中，設(shè)H=1，U=1.0，Tb=Tt=1.0。上板處的馬赫數(shù)Mat=0.1。雷諾數(shù)的值依賴于H、U以及流體黏性(假設(shè)為常量)，取為Re=100。計(jì)算域?yàn)?≤x≤2和0≤y≤H，并被x=0.6和x=0.4兩條直線分為3個(gè)子區(qū)域。令中間子區(qū)域的水平方向的中線以Δy=0.1sin(0.5t)運(yùn)動(dòng)，且其網(wǎng)格變形區(qū)域限制在0.1≤y≤0.9之間。剩余兩個(gè)子區(qū)域網(wǎng)格保持靜止。上下采用無(wú)滑移恒溫壁面邊界條件，左右視為周期性邊界條件。測(cè)試中用到了3套網(wǎng)格，分別含有4/4、16/16和64/64個(gè)三角形/四邊形混合網(wǎng)格單元。

圖10顯示了在64/64混合網(wǎng)格上用四階精度格式仿真得到的穩(wěn)態(tài)平板Couette流馬赫數(shù)(Ma)云圖,圖中淺灰色線為網(wǎng)格線，深灰色線為處理器界面。很明顯，網(wǎng)格的運(yùn)動(dòng)以及滑移界面網(wǎng)格的錯(cuò)位并沒(méi)有引入任何可見(jiàn)的擾動(dòng)。馬赫數(shù)云圖在縱向?qū)哟畏置?，且沿水平方向保持不變，與流場(chǎng)解析解一致。

誤差和精度的計(jì)算方法與Euler vortex流中一致，表2顯示了用速度u算得的結(jié)果。同樣地，隨著網(wǎng)格密度的增加，L1和L2誤差都相應(yīng)地逐漸減小，精度也在3和4左右。因此，對(duì)于黏性流動(dòng)而言，該滑移網(wǎng)格方法依然可以達(dá)到理想的精度。

為了便于說(shuō)明高精度求解器的低耗散特性，將圖10所示的流場(chǎng)區(qū)域全部用四邊形單元進(jìn)行劃分，中間子區(qū)域網(wǎng)格的運(yùn)動(dòng)規(guī)律保持不變。采用二階精度格式分別在單元數(shù)為200和800的兩套網(wǎng)格上進(jìn)行仿真，四階精度格式在單元數(shù)為50和200的兩套上進(jìn)行仿真，從而兩種精度格式兩次仿真的自由度分別都為800和3 200。圖11顯示了4次仿真密度的L1誤差隨迭代步數(shù)的變化曲線。由圖可知，在相同自由度的情況下，四階精度格式誤差始終小于二階精度格式；甚至自由度為800的四階精度格式仿真誤差也明顯小于自由度為3 200的二階精度格式仿真誤差，驗(yàn)證了高階精度格式的低耗散性。

圖10 平板Couette流馬赫數(shù)云圖Fig.10 Mach number contours for planar Couette flow

表2平板Couette流動(dòng)仿真的誤差和精度的階數(shù)(基于速度u計(jì)算)

Table2Errorsandordersofaccuracy(basedonuvelocity)forplanarCouetteflowsimulation

格式單元數(shù)L1誤差L1誤差的精度階L2誤差L2誤差的精度階三階4/42.03×10-52.36×10-516/162.54×10-63.002.89×10-63.0364/643.05×10-73.063.47×10-73.06四階4/42.46×10-73.91×10-716/161.77×10-83.792.71×10-83.8564/641.18×10-93.901.66×10-94.03

圖11 平板Couette流動(dòng)仿真密度的L1誤差隨迭代步數(shù)變化Fig.11 Variation of L1 errors of density with iteration steps for planar Couette flow simulation

3.3 定常均勻來(lái)流數(shù)值模擬

如圖12(a),計(jì)算區(qū)域?yàn)?0×10的矩形，劃分為三角形和四邊形混合網(wǎng)格，上下左右都設(shè)為周期性邊界條件。兩條黑色的滑移界面將流場(chǎng)分為3部分：中間子區(qū)域的水平中線(y=5，3≤x≤7)以方程Δy=sin(0.5t)隨時(shí)間運(yùn)動(dòng)，則周期T=4π,上下兩側(cè)區(qū)域相應(yīng)地拉伸或壓縮，變形區(qū)域限定為|y-5|≤4；左右兩側(cè)子區(qū)域內(nèi)部網(wǎng)格點(diǎn)按式(24)所示規(guī)律運(yùn)動(dòng)(其中，IV為網(wǎng)格點(diǎn)的序號(hào))。圖12顯示了仿真初始時(shí)刻和仿真過(guò)程中兩個(gè)不同時(shí)刻的網(wǎng)格情況，其中圖12(a)為初始時(shí)刻網(wǎng)格。

(24)

流場(chǎng)初始化為均勻來(lái)流，馬赫數(shù)取0.3。采用四階精度格式離散Euler方程，時(shí)間步長(zhǎng)取為1.0×10-3。圖13給出了滿足GCL(|J|作為未知量在時(shí)間推進(jìn)中迭代更新)和不滿足(|J|由網(wǎng)格坐標(biāo)信息直接算得)情況下，密度的L1和L2誤差對(duì)比。從圖中可以看出，在網(wǎng)格有壓縮和扭曲的情況下，滿足GCL的自由流仿真算例能夠維持在機(jī)器零附近沒(méi)有增加，證明了求解器的自由流保持特性。

圖12 3個(gè)不同時(shí)刻對(duì)應(yīng)的流場(chǎng)網(wǎng)格Fig.12 Meshes of flow field at three different time instants

圖13 四階精度格式自由來(lái)流仿真密度誤差變化Fig.13 Density error history of 4th order scheme for free stream flow simulation

3.4 單個(gè)圓柱的VIV數(shù)值模擬

此算例對(duì)彈性連接的單個(gè)圓柱進(jìn)行VIV仿真，并與其他文獻(xiàn)中的結(jié)果進(jìn)行比較，驗(yàn)證方法的有效性。計(jì)算狀態(tài)為：自由來(lái)流馬赫數(shù)Ma∞=0.2，雷諾數(shù)Re=100，圓柱的質(zhì)量比為m*=10.0，速度比U*=4.92，阻尼比ζ=0(變量定義見(jiàn)文獻(xiàn)[12])。

圓柱直徑D=1，計(jì)算區(qū)域長(zhǎng)與寬為60D×20D，圓柱圓心距入流邊界10D，兩個(gè)滑移界面將流場(chǎng)分為3個(gè)子區(qū)域。流場(chǎng)劃分成6 696個(gè)混合網(wǎng)格單元，圓柱周邊和尾跡區(qū)網(wǎng)格進(jìn)行了加密。邊界條件設(shè)置為：狄利克雷入口邊界條件、壓力出口邊界條件、上下為對(duì)稱邊界條件、圓柱表面為無(wú)滑移絕熱壁面。時(shí)間推進(jìn)步長(zhǎng)為Δt=5.0×10-4。圓柱沿來(lái)流方向(x方向)和垂直來(lái)流方向(y方向)都可以振動(dòng)，x方向振動(dòng)的控制方程也為式(4)，控制參數(shù)和y方向參數(shù)一致。因幅值很小，單純的解析變形函數(shù)便可滿足x方向的變形要求。

表3 單個(gè)圓柱VIV數(shù)值模擬響應(yīng)Table 3 Response of single cylinder VIV simulation

圖14 t從220到600過(guò)程中單個(gè)圓柱仿真CL和CD的利薩茹曲線Fig.14 Lissajous curve of CL and CD for single cylinder from t=220 to t=600

3.5 一列4個(gè)圓柱的VIV數(shù)值模擬

此算例對(duì)置于一列的4個(gè)圓柱進(jìn)行VIV仿真，目的是測(cè)試求解器的加速效果，并且驗(yàn)證求解器處理流場(chǎng)中存在多個(gè)物體，相互之間具有較大相對(duì)運(yùn)動(dòng)情況的能力。相鄰的兩個(gè)圓柱間隔ΔL/D=1.0，自由來(lái)流馬赫數(shù)Ma∞=0.1，雷諾數(shù)Re=200，每一個(gè)圓柱的質(zhì)量比都為m*=5.0，速度比U*=9，阻尼比ζ=0。

計(jì)算域的范圍為74D×50D，前面的圓柱距離入流邊界20D，用5個(gè)滑移界面將流場(chǎng)分為6個(gè)子區(qū)域。流場(chǎng)一共劃分為58 436個(gè)混合網(wǎng)格單元，其中中間4個(gè)動(dòng)網(wǎng)格區(qū)域分別有7 116個(gè)單元，前面的靜止網(wǎng)格區(qū)域有8 708個(gè)單元，最后的靜止網(wǎng)格區(qū)域有21 262個(gè)單元，對(duì)圓柱周邊和尾跡區(qū)網(wǎng)格進(jìn)行了加密。邊界條件設(shè)置為：狄利克雷入口邊界條件、壓力出口邊界條件、上下為對(duì)稱邊界條件、圓柱表面是無(wú)滑移絕熱壁面。時(shí)間推進(jìn)步長(zhǎng)為Δt=2.0×10-4。

圖15顯示了一定時(shí)間內(nèi)段內(nèi)4個(gè)圓柱的位移隨時(shí)間變化曲線。由于流場(chǎng)的渦結(jié)構(gòu)比較復(fù)雜(見(jiàn)圖16)，所以圓柱的位移響應(yīng)也比較豐富。平均振幅的相對(duì)大小在0.6～1.0之間，相對(duì)于圓柱間隔而言是比較大的。圖16是某一時(shí)刻流場(chǎng)的渦量云圖，灰線表示處理器的分區(qū)界面，水平的點(diǎn)畫(huà)線表示4個(gè)圓柱初始所在的位置。從圖中可以看出求解器能夠很清晰地刻畫(huà)流場(chǎng)，證明了本文動(dòng)網(wǎng)格方法的可靠性。

圖15 4個(gè)圓柱的位移隨時(shí)間變化曲線Fig.15 Displacement history of four cylinders

測(cè)試并行求解器加速效果時(shí)用到了多達(dá)240個(gè)處理器，圖17顯示了加速效率(Speedup)隨所用處理器數(shù)量(Number of processors)的變化關(guān)系，其中Linear表示理想情況下所能達(dá)到的線性加速效果。對(duì)于具有NM-1個(gè)滑移界面的網(wǎng)格，最多可以有2(NM-1)個(gè)處理器去更新滑移界面處的網(wǎng)格面-粘接元對(duì)應(yīng)關(guān)系以及當(dāng)?shù)?全局粘接元對(duì)應(yīng)關(guān)系，剩余的處理器需要等待，造成了一點(diǎn)資源浪費(fèi)。由圖17可見(jiàn)，在120個(gè)處理器以內(nèi)，加速基本都能保持線性。隨著處理器數(shù)量的繼續(xù)增加，處理器之間的信息傳遞時(shí)間以及更新粘接元對(duì)應(yīng)關(guān)系的等待時(shí)間占總時(shí)間的比重也在增加，影響了加速效果。

同時(shí)，對(duì)比3條加速曲線可以發(fā)現(xiàn)，精度越高，加速效果越好。因?yàn)殡S著精度增加，處理器內(nèi)部需要處理的自由度以O(shè)(N2)增加，而邊界交互自由度只以O(shè)(N)增加，使得處理器之間信息傳遞以及等待時(shí)間所占比重也相對(duì)減小。同理，對(duì)網(wǎng)格進(jìn)行加密時(shí)，也會(huì)使得線性加速的效果能夠保持到處理器數(shù)量更大的時(shí)候。因此該求解器適用于進(jìn)行大規(guī)模并行仿真。

圖16 某一時(shí)刻流過(guò)4個(gè)渦激振蕩圓柱流場(chǎng)的渦量云圖Fig.16 Vorticity contours at one instant for a row of four vibrating cylinders

圖17 并行求解器加速曲線Fig.17 Speedup curves of parallel solver

4 結(jié) 論

1) 提出了一種適用于渦激振蕩仿真的并行高精度方法，可用于研究大振幅、小間距等極端VIV問(wèn)題。

2) SD方法在三角形/四邊形非結(jié)構(gòu)混合網(wǎng)格的應(yīng)用，使得求解器對(duì)于復(fù)雜外形處理能力增強(qiáng)；利用了非均勻滑移網(wǎng)格方法，克服傳統(tǒng)協(xié)調(diào)網(wǎng)格方法所受到的變形幅度的限制，可以方便地實(shí)現(xiàn)網(wǎng)格變形，保證良好的網(wǎng)格質(zhì)量，并滿足GCL條件。

3) 數(shù)值算例表明，該方法對(duì)于無(wú)黏流動(dòng)和黏性流動(dòng)都能保持原SD方法的高階精度，證明了求解器具有低耗散性的特性，能夠精確捕捉流場(chǎng)細(xì)節(jié)，從而使得VIV仿真準(zhǔn)確。

4) 并行測(cè)試證明該求解器具有較好的加速效果，且所用數(shù)值格式精度越高、網(wǎng)格數(shù)量越大，加速效果越好，因此該求解器適合于大規(guī)模并行計(jì)算。