王吉飛,萬德成
(1. 上海交通大學 船舶海洋與建筑工程學院 海洋工程國家重點實驗室,上海 200240; 2. 高新船舶與深海開發裝備協同創新中心,上海 200240)
三維頂板斜向驅動方腔流有限元并行計算
王吉飛1,2,萬德成1,2
(1. 上海交通大學 船舶海洋與建筑工程學院 海洋工程國家重點實驗室,上海 200240; 2. 高新船舶與深海開發裝備協同創新中心,上海 200240)
基于PETScFEM開源代碼,采用分步有限元算法和區域分解法,并行計算了不同雷諾數下的三維頂板斜向驅動方腔流問題。計算結果表明,當頂板沿其對角線方向運動時,流體與下游側壁發生斜向碰撞后在下游對角處匯聚并形成射流,該射流在與底面、上游側壁碰撞后形成較為復雜的渦流結構。雷諾數大小對三維頂板斜向驅動方腔流的渦流場結構形態具有重要影響。并行性能分析表明區域分解法能有效地提高三維粘性方腔流的計算速度。
三維方腔流;斜向角;有限元;區域分解;并行計算
頂板驅動方腔流可模擬由某一壁面運動或外部流場運動所引起的方腔內部環流。這類流動在航空、環境、工業、交通運輸、海洋工程等領域廣泛存在,研究其流動特性具有重要的工程應用價值。Aidun 等[1]指出頂板驅動方腔流在涂鍍和熔紡等工業裝置中起到重要的作用。在液壓管路運輸中,適當部位設置一個較大的方腔可減少或吸收液壓泵的壓力和流量脈動對系統的影響,防止或減少液壓沖擊現象的發生。在航空領域,前緣縫翼、后緣縫翼、起落架箱和釘狀洞穴等所形成的方腔流動是機身噪音的重要來源。在潛艇噪聲中,弦外孔和排水孔等所產生的內部環流也是其水動力噪聲的來源之一。研究方腔流對流動機理的探討也有非常重要的意義,因為其簡單的設置能展示幾乎所有重要的流體力學現象,如角渦、展向渦、分叉、轉捩、湍流等。方腔流模型的簡單性使得其結果的分析和比較更容易,同時極易擴展到全雷諾數范圍。
隨著實驗技術和計算機性能的不斷提高,人們對三維頂板驅動方腔流進行了越來越多的研究,其中對頂板沿平行于側壁方向運動引起的方腔流研究得最為廣泛。Prasad和Koseff[2]對頂板驅動方腔流進行了開拓性的實驗研究。Jiang等[3],Wong等[4]以及其他學者對中、低雷諾數下的方腔流進行了數值模擬,Leriche等[5],Hachem等[6],Zang等[7]對高雷諾數的方腔流進行了數值模擬等。此外,Shankar等[8]對方腔流的研究進展做了很好的總結。
但在現實情況中,頂板運動方向或外部流場方向并不一定平行于方腔側壁方向,而是與其有一定偏角。該偏角的存在使得方腔內部渦流場結構發生重大變化。研究頂板斜向驅動方腔流具有工程實用價值,能更好地防止或減少管路輸運中的液壓沖擊現象,更好地降低飛機機身或潛艇噪聲等,同時對流動機理的探討也有重要的意義。國內外學者對頂板斜向驅動方腔流的研究較少,主要集中在飛機機身降噪方面。Disimile等、Czech等分別對飛機著陸裝置(可簡化為頂板驅動方腔流)的噪聲進行了實驗研究,結果表明當外部流方向與方腔側壁成一定偏角時,不管流動是否振蕩,方腔中均存在較強的聲學壓力。Pocitsky[9]基于商業軟件Fluent研究了頂板沿其對角線方向運動所引起的方腔流,其雷諾數包括100,400,700,2 000和40 000,結果顯示頂板斜向驅動方腔流具有更復雜的渦流結構。為了便于與前人的數值計算結果進行比較,本文著重研究頂板沿其對角線方向運動所引起的方腔流。
根據控制方程的離散方式不同,計算流體力學數值離散方法大體上可分為三個分支[10]:有限差分法、有限體積法和有限元法。相比于其他方法,有限元方法具有如下特點和優勢[11,12]:1)有限元法具有理論完整可靠,形式單純、規范,程序標準化,精度和收斂性得到保證等優點;2)相比于有限差分法,有限元法對于求解區域的單元剖分沒有特別的限制,因此特別適合處理具有復雜邊界流場的區域;3)相比于有限體積法,有限元法可以通過提高單元插值多項式的次數來提高解的精度,且對單元正則性要求不高。基于此,本文采用有限元法對不可壓粘性流場進行數值模擬。
在對不可壓粘性流體進行有限元分析時,速度和壓力插值空間函數需滿足inf-sup相容性條件,如不滿足,則可能產生虛假的壓力振蕩。為解決壓力穩定性問題,學者們相繼提出了不同的方法,如速度和壓力的不同階插值,罰函數法,附加壓力穩定項等方法。由Chorin等和Temam等提出的基于POISSON投影的分步算法則可采用無需滿足相容性條件的速度壓力插值空間函數。Guermond等[13]進一步驗證,當速度和壓力采用同階插值函數時,只有在時間步長不太小時才能得到穩定的壓力場。Codina等[14]得到了類似的結論,并進一步提出了具有壓力穩定性的分步算法。由于有限元分步算法具有計算效率高,精度得到保證,編程較為簡單等優點[15],文中采用該方法對三維方腔流問題進行數值模擬。
計算能力不足是計算流體力學發展所面臨的一大挑戰。為解決大規模計算問題,并行化是大幅度提高計算效率的最有效手段。常用的并行化方法有線性代數方程組的并行計算和區域分解法兩種。并行計算工具箱PETSc(Portable, Extensible Toolkit for Scientific Computation)提供了并行計算線性代數方程組的庫函數。軟件包 PETScFEM 基于并行計算工具箱 PETSc,將非重疊區域分解法應用于計算流體力學的有限元分析中,實現了不可壓粘性流體力學問題的并行計算。
本文基于PETScFEM開源代碼[16-17],采用分步有限元算法對不可壓粘性流體進行求解,并結合區域分解法并行計算了不同雷諾數下的頂板斜向驅動方腔流問題。數值模擬時首先對計算區域進行有限元離散,不可壓 Navier-Stokes 方程中的速度和壓力采用同階的線性形函數和權函數,時間離散采用后退歐拉法;采用基于 Poisson 投影的分步算法對速度和壓力進行解耦,并形成三個子方程;計算網格通過非重疊的區域分解法進行剖分,并分配到不同的CPU中并行計算,MPI(Message Passing Interface)庫函數用于各子區域間并行計算數據通信,離散后的線性代數方程組采用 GMRES(Generalized Minimal RESidual)結合 Jacobi 預處理方法進行計算。在之前的工作中[18],已采用該方法對頂板正向驅動方腔流問題進行了數值模擬,計算結果表明該方法計算效率高,精度較好。在此基礎上,文中的主要工作是對雷諾數為100、400、1 000和2 000的頂板斜向驅動方腔流進行并行數值模擬。
本文結構安排如下:首先描述不可壓粘性流體控制方程及分步有限元算法,其次介紹非重疊區域分解并行算法,然后對不同雷諾數下的頂板斜向驅動方腔流進行數值模擬,給出速度矢量圖、壓力等值面、流線、渦量等值線以及中心垂線上的速度剖面等流場結構信息,并進行并行性能分析,最后給出結論。


圖1 三維頂板斜向驅動方腔流計算模型Fig. 1 Computational model of 3D lid-driven cubic cavity flow at yaw
粘性流體流動可用 Navier-Stokes 方程進行求解。對于不可壓縮流體流動,其控制方程包含動量守恒方程和質量守恒方程:
?tu+u·u+

為了方便描述有限元空間離散式,引入如下記號:
a(u,v):=ν(u,v),b(q,v):=(q,v),c(u,v,w):=(u·v,w)
(?tu,v)+c(u,u,v)+a(u,v)-b(p,v)=0 ?v∈V,
b(q,u)=0 ?q∈Q
采用梯形法則對方程進行整體時間離散,即同時求解速度和壓力。簡單起見,將時間 [0,T] 劃分為N個均勻的時間步,時間步長為δt。令f為關于時間的通用函數,fn為函數f在tn=nδt時刻的值,fn+θ:=θfn+1+(1-θ)fn,δtfn:=(fn+1-fn)/δt,其中θ∈[0,1],則弱解方程的時間離散式為:


當θ=1/2時,該方程對應于二階的 Crank-Nicolson 時間離散,當θ=1時,該方程對應于向后 Euler 時間離散。
令Vh為V的近似有限元空間,Qh為Q的近似有限元空間,Vh空間中的函數要求為分段連續多項式,Qh空間中函數的連續性不要求,則方程(3)~(4)的有限元空間離散式為:


在基于壓力 Poisson 方程的分步算法中,當時間步長較大時,速度和壓力插值函數不需要滿足 LBB 條件,因此在本文的計算中速度和壓力插值函數均采用一階線性多項式。
方程(5)~(6)的矩陣形式可表示為:
式中:U和P分別表示結點速度和壓力數組,M為質量矩陣,K為包含擴散項和對流項的矩陣,G為梯度矩陣,D為散度矩陣。
方程(7)~(8)等價于如下的方程組:




圖2 非重疊區域分解法計算模型(2個子區域)Fig. 2 Computational model of non-overlapping domain decomposition
令Aii,i=1,2,…,n為各子區域Ωi內部結點所形成的矩陣,ALL=diag[A11,A22,…,Ann]為分塊對角矩陣,AII為各子區域間交界面結點所形成的矩陣。相應地,令x=(xL,xI)T,y=(yL,yI)T,則代數方程組可分裂為:
式中:ALI和AIL為子區域內部結點和交界面結點的連接關系矩陣。該方程組的求解可以分為兩步完成,首先求解方程(17)獲得交界面結點的變量值xI,各子區域間的數據通信通過 MPI 庫函數來完成,然后再并行求解方程(18),得到各子區域內部結點的變量值xL:
首先模擬雷諾數為1 000時的頂板斜向驅動方腔流,進行了網格收斂性驗證,并給出速度矢量圖、壓力等值面、流線、渦量等值線等流場信息。其次模擬雷諾數為100,400,2 000等不同頂板斜向驅動方腔流,比較不同雷諾數對頂板斜向驅動方腔流的影響。
首先模擬雷諾數為1 000時的頂板斜向驅動方腔流。計算采用三套不同精度的均勻六面體網格進行網格收斂性驗證,粗、中、細網格的單元數分別為 36×36×36、48×48×48、64×64×64。計算時間步長取為 0.05,迭代1 000步均能得到收斂的結果。計算結果表明,中等精度網格已能得到較好的結果,而計算時間卻比細網格少很多。圖3所示結果為中等精度網格計算結果。
雷諾數為1 000時頂板斜向驅動方腔流的流線分布如圖3 所示。從正視圖可以看出,流體由頂板驅動,以45°角撞擊x=1的下游側壁和y=1的下游側壁,撞擊后的流體在由兩側壁構成的角落處匯聚,并形成一股強大的射流射向底面。從后視圖可以看出,射流撞向底面后在底面呈90°扇形散開。流體在遇到x=0的上游側壁和y=0的上游側壁后轉而向上流動,最終匯入頂板附近流體,完成一次循環。從側視圖可以看出,由于頂板驅動流體的匯聚和射流的形成,使得主渦變得傾斜而狹長,且有流體從下游側壁中心附近匯入主渦中。

圖3 雷諾數為1 000時流線分布的不同視角(左:正視圖,中:后視圖,右:側視圖)Fig. 3 Streamlines at Re=1 000: front view (left); back view (middle); oblique view (right)

從圖5可看到方腔流不同剖面的速度矢量投影。在CP剖面中,由于流體的堆積和射流的形成,可以看到主渦在平面的左下方,而右下角可看到次級渦的存在。在CP2剖面中,渦心位于剖面中央,平均速度較CP平面小。在CP3剖面中,渦心位于平面右上角,且大部分區域的速度矢量向上。從PP剖面和MP剖面的速度矢量投影中,可以看到二次渦的存在。本文結果與Povitsky的計算結果[9]基本一致。
圖6 所示為雷諾數為1 000時頂板斜向驅動方腔流速度大小等值面。從圖6(a)中可以看到,由頂板帶動的流體在由x=1的下游側壁和y=1的下游側壁形成的垂直角落里匯聚形成一股速度較大的射流。圖6(b)則顯示了由頂板帶動的流體在撞擊下游側壁后匯聚的過程。圖6(c)顯示了射流在撞擊底面之后朝兩個方向散開的低速射流。圖6(d)、6(e)和6(f)分別顯示了頂板驅動流體在到達底面之后呈扇形面散開,在遇到x=0的上游側壁和y=0的上游側壁后轉而向上流動,由于流體運動空間的增大和壁面粘性的阻滯作用,這一過程中流體的速度逐漸降低。


圖4 不同剖面的速度矢量和流線投影及其法向渦量分量(Re=1 000)Fig. 4 Perspective 3D solution summary at Re=1 000

圖5 不同剖面的速度矢量投影分布(Re=1 000)Fig. 5 2D planar projections of velocity vector at Re=1 000 on different planes
圖7 所示為雷諾數為1 000時頂板斜向驅動方腔流的壓力等值面。圖8所示為雷諾數為1 000時不同精度網格計算結果(沿 CP 剖面中垂線的水平速度分量)比較。從圖8中可以看到,從粗網格到細網格,計算結果逐漸收斂,中等精度網格與細網格計算結果基本保持一致。由于中等精度網格已能得到較好的計算精度,而網格單元數卻比細網格少很多,可以減少大量的計算時間,因此在本文其它的工況計算中均采用中等精度網格進行計算。從速度剖面看,水平速度分量在底面附近存在一個峰值,該峰值由底面射流產生。在距底面0.1~0.3的區域速度絕對值均較大,由流體的匯聚擠壓引起。

圖6 頂板斜向驅動方腔流速度大小等值面(Re=1 000)Fig. 6 Iso-surfaces for different velocity magnitudes at Re=1 000

圖7 壓力等值面(Re=1 000)Fig. 7 Pressure iso-surfaces at Re=1 000

圖8 沿CP 剖面中垂線的水平速度分量(Re=1 000)Fig. 8 Horizontal component velocity profile at Re=1 000
本文還模擬了雷諾數為100、400和2 000時的頂板斜向驅動方腔流,并考察不同雷諾數對頂板斜向驅動方腔流的影響。本節計算采用48×48×48的均勻六面體網格,單元特征尺度為0.020 83。計算時間步長取為0.05,迭代1 000步均能得到收斂的結果。
圖9所示為雷諾數為100、400和2 000時頂板斜向驅動方腔流的流線分布。從圖中可以看出,隨著雷諾數的增加,流線由簡單變為復雜,特別是主渦流線,當雷諾數為2 000時變得雜亂。
圖10所示為不同雷諾數下CP剖面內速度矢量投影。從圖中可以看出,隨著雷諾數的增加,主渦渦心位置不斷變化,當雷諾數為100時,主渦位于中部偏上的位置,當雷諾數為400時,主渦中心在右上角,當雷諾數為1 000時,主渦變為傾斜狹長的卵形,渦心在平面的左下方,當雷諾數達到2 000時,渦心位置更靠近底面。此外,當雷諾數為400時,可以看到較為明顯的下游二次渦。當雷諾數達到2 000時,可以明顯地觀察到上游二次渦的形成。

圖9 不同雷諾數(從上至下:Re=100,400,2 000)流線分布的不同視角(從左至右:正視圖,后視圖,側視圖)Fig. 9 Streamlines for different Reynolds numbers

圖10 不同雷諾數下CP剖面內速度矢量投影Fig. 10 2D planar projections of velocity vector on CP plane for different Reynolds numbers
圖11所示為不同雷諾數下 MP 剖面內速度矢量投影(本文結果)。從圖中可以看出,隨著雷諾數的增加,二次流變得越來越復雜,當雷諾數為1 000時,可以觀察到3對二次渦,而當雷諾數增加到2 000時,可以觀察到6對二次渦。圖12所示為Povitsky的計算結果[9]。將本文結果與參考結果進行比較,可以看出對應的漩渦個數及漩渦形態均保持一致。

圖11 不同雷諾數下MP剖面內速度矢量投影(本文結果)Fig. 11 2D planar projections of velocity vector on MP plane for different Reynolds numbers (our results)

圖12 不同雷諾數下MP剖面內速度矢量投影(參考結果)Fig. 12 2D planar projections of velocity vector on MP plane for different Reynolds numbers (reference results)
圖13所示為不同雷諾數下沿CP剖面中垂線的水平速度分量。從圖中可以看到,隨著雷諾數的增加,邊界層的厚度逐漸減小。當雷諾數為100時,粘性作用較強,壁面的阻滯效應較大,速度剖面變化較為平緩。當雷諾數為400時,底面附近已形成射流,因此出現了較大的峰值。當雷諾數為1 000時,底面附近的射流峰值有所減小,但流體的堆積擠壓現象較為明顯,導致距底面0.1~0.3附近的速度絕對值較大。當雷諾數達到2 000時,流體堆積擠壓現象更為嚴重,使得流體在距底面0.2附近出現峰值。

圖13 不同雷諾數下沿CP剖面中垂線的水平速度分量Fig. 13 Horizontal velocity component along the vertical midline of the CP plane at different Reynolds numbers
不同雷諾數下沿 CP 剖面中垂線的水平速度分量數值如表1所示。
本文采用網格非重疊區域分解法結合消息傳遞的模式進行并行計算,非重疊網格區域劃分的方式如圖14所示。以雷諾數為1 000的頂板斜向驅動方腔流計算為例,取時間步長為0.05,迭代1 000步,不同進程數(從1到8)的計算時間及加速比如表2所示,相對應的曲線如圖15所示。從圖中可以看出,隨著進程數的增加,計算時間不斷減少,并行加速比不斷增大,但由于進程間數據傳遞量的增加,加速比的增長速度比線性增長速度為小。

表1 不同雷諾數下沿 CP 剖面中垂線的水平速度分量Tab. 1 Horizontal velocity component along the vertical midline of the CP plane at different Reynolds numbers


圖14 不同塊數的區域分解網格Fig. 14 Domain decomposition meshes for different domains

進程數12345678時間/s4323020520156601286010540951791288171加速比12.10672.76053.36164.10154.54244.73605.2907

圖15 并行性能分析Fig. 15 Analysis of the parallel performance
采用分步有限元算法結合非重疊區域分解法并行計算了不同雷諾數(100、400、1 000 和2 000)下的頂板斜向(斜向角為 45°)驅動方腔流問題。從計算結果可以看出,頂板斜向驅動方腔流具有較為復雜的渦流結構。由頂板驅動的流體以45°角撞擊x=1的下游側壁和y=1的下游側壁后匯聚成一股速度較大的射流,該射流在底面呈扇形散射,在遇到x=0的上游側壁和y=0的上游側壁后轉向上流動,最終匯入頂板附近流體并進入循環。在這一循環過程中,由于主流運動方向與側壁呈45°偏角,使得二次流現象較為明顯和復雜。且隨著雷諾數的增加,流動更為劇烈,整個渦流場結構也更為復雜。
計算表明,采用基于Poisson投影的分步算法后,應用傳統的有限元方法即能很好地求解不可壓粘性方腔流問題。同時,結合網格區域分解和消息傳遞的并行計算模式,可以得到很好的并行計算性能。開源軟件PETScFEM為計算不可壓粘性流體問題提供了很好的工具,且為進一步開發提供了很好的平臺。
[1] AIDUN C K, TRIANTAFILLOPOULOS N G, BENSON J D. Global stability of a lid-driven cavity with through flow: flow visualization studies [J]. Physics of Fluids, 1991, 3: 2081-2091.
[2] PRASSAD A K, KOSEFF J R. Reynolds number and end-wall effects on a lid-driven cavity flow [J]. Physics of Fluids, 1989,1(2): 208.
[3] JIANG B N, LIN T L, POVINELLI L A. Large scale computation of incompressible viscous flows by least-squares finite element method [J]. Computer Methods in Applied Mechanics and Engineering, 1994, 114(3-4): 213-231.
[4] WONG K L, BAKER A J. A 3D incompressible Navier-Stokes velocity-vorticity weak form finite element algorithm [J]. International Journal for Numerical Methods in Fluids, 2002, 38(9): 9-23.
[5] LERICHE E. Direct numerical simulation in a lid-driven cubical cavity at high reynolds number by a chebyshev spectral method [J]. Journal of Scientific Computing, 2006: 335-345.
[6] HACHEM E, RIVAUX B, KLOCZKO T, et al. Stabilized finite element method for incompressible flows with high Reynolds number [J]. Journal of Computational Physics, 2010, 229: 643-665.
[7] ZANG Y, STREET RL, KOSEFF JR. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows [J]. Physics of Fluids, 1993, 5: 3186.
[8] SHANKAR PN, DESHPANDE MD. Fluid mechanics in the driven cavity [J]. Annual Review of Fluid Mechanics, 2000, 32: 93-136.
[9] POVITSKY A. Three-dimensional flow in cavity at yaw [J]. Nonlinear Analysis, 2005, 63(5-7): 1573-1584.
[10] 任安祿. 不可壓縮粘性流場計算方法 [M]. 北京: 國防工業出版社, 2003.(REN Anlu. Numerical methods in incompressible viscous flow [M]. Beijing: National Defence Industry Press, 2003.(in Chinese))
[11] 章本照,印建安,張宏基. 流體力學數值方法 [M]. 北京: 機械工業出版社, 2003. (ZHANG Benzhao, YIN Jianan, ZHANG Hongji. Numerical methods in fluid dynamics [M]. Beijing: China Machine Press, 2003.(in Chinese))
[12] 王獻孚,周樹信,陳澤梁,等. 計算船舶流體力學 [M]. 上海: 上海交通大學出版社, 1991. (WANG Xianfu, ZHOU Shuxin, CHEN Zeliang, et al. Computational ship hydrodynamics [M]. Shanghai: Shanghai Jiao Tong University Press, 1991.(in Chinese))
[13] GUERMOND J L, QUARTAPELLE L. On stability and convergence of projection methods based on pressure poisson equation [J]. International Journal for Numerical Methods in Fluids, 1998, 26: 1039-1053.
[14] CONDINA R. Pressure stability in fractional step finite element methods for incompressible flows [J]. Journal of Computational Physics, 2001, 170(1): 12-40.
[15] 劉淼兒,任玉新,張涵信. 數值求解不可壓縮流動的投影方法研究進展 [J]. 力學進展, 2006, 36(4): 591-598. (LIU Miaoer, REN Yuxin, ZHANG Hanxin. Review on the projection methods in the numerical solution of the incompressible flow [J]. Advances in Mechanics, 2006, 36(4): 591-598.(in Chinese))
[16] PAZ R R, NIGRO N M, STORTI M A. On the efficiency and quality of numerical solutions in CFD problems using the interface strip preconditioner for domain decomposition methods [J]. International Journal for Numerical Methods in Fluids, 2006, 52(1): 89-118.
[17] SONZOGNI V E, YOMMI A M, NIGRO N M, et al. A parallel finite element program on a Beowulf cluster [J]. Advances In Engineering Software, 2002, 33(7-10): 427-443.
[18] WANG J F, WAN D C. Parallel simulation of 3D lid-driven cubic cavity flows by finite element method [C]// Proceedings of the 21st International Offshore and Polar Engineering Conference. Hawaii:[s.n.], 2011: 644-651.
Parallel simulation of 3D lid-driven cubic cavity flows at yaw by finite element method
WANG Jifei1,2, WAN Decheng1,2
(1. State Key Laboratory of Ocean Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China; 2. Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, Shanghai 200240, China)
The fractional step finite element method and domain decomposition method are applied to simulate in parallel the 3D lid-driven cubic cavity flows at yaw based on the open source codes PETScFEM. When the lid moves along its diagonal, the driven fluid impinges in the spatial angle formed by the downstream side walls, which leads to formation of the jet flow. A system of vortices is caused when the jet flow impinges the bottom wall and upstream side walls. Different Reynolds numbers are investigated, which show the significant influence for this flow. Parallel performance analysis reveals that the domain decomposition method can efficiently speed up the simulation of 3D lid-driven cavity flows.
3D cavity flow, yaw angle, finite element method, domain decomposition, parallel simulation
O351
A
10.16483/j.issn.1005-9865.2015.02.001
1005-9865(2015)02-0001-12
2014-01-06
國家自然科學基金資助項目(51379125,11432009,51490675);上海高校特聘教授(東方學者)崗位跟蹤計劃項目(2013022);國家重點基礎研究發展計劃(973計劃)項目(2013CB036103);工信部高技術船舶科研項目;上海交通大學高性能計算中心(HPC)資助項目
王吉飛(1984-),男,重慶人,博士研究生,從事計算流體力學研究。E-mail:wangjifei2000@126.com
萬德成。E-mail:dcwan@sjtu.edu.cn