李裕龍，儲南洋，胡?湛，姚?鵬，鄧?銳

?

Common-Refinement方法在三維流固耦合數(shù)值計算中誤差分析與并行效率

李裕龍1，儲南洋1，胡?湛1，姚?鵬1，鄧?銳2

(1. 中山大學海洋科學學院，廣州 119077；2. 中山大學海洋工程與技術學院，廣州 119077)

本文成功實現(xiàn)了基于非匹配網(wǎng)格的三維Common-Refinement方法，可將其應用于離散后流體與彈性體的非重疊子區(qū)域之間的接觸面，用于數(shù)值求解三維非定常流固耦合問題．首先，為了滿足兩者之間液體和彈性耦合界面間牽引力的平衡條件，本文系統(tǒng)研究了Common-Refinement方法的空間插值的準確性和可靠性．本文根據(jù)一系列的不同拓撲結(jié)構(gòu)的網(wǎng)格劃分方案，通過改流體和結(jié)構(gòu)網(wǎng)格之間的網(wǎng)格匹配系數(shù)系統(tǒng)地評估Common-Refinement方法的準確性和精度，證明了Common-Refinement方法沿界面的流體和固體網(wǎng)格的數(shù)據(jù)傳輸具有穩(wěn)定的二階精度．進而將Common-Refinement方法應用于非定常流固耦合算例，研究了非定常計算中的誤差累積以及Common-Refinement方法的并行效率．最后本文將本方法應用于三維標準的圓柱體-彈性板問題，并與文獻中標準解進行了對比．求解結(jié)果證明了Common-Refinemen方法在三維非定常流固耦合數(shù)值計算中能夠保持全場二階精度．

流固耦合；誤差分析；Common-Refinement；非匹配網(wǎng)格

涉及多個物理場相互作用的科學和工程數(shù)值模擬中，通常需要精確和守恒的方法，通過非匹配離散網(wǎng)格來傳輸物理數(shù)據(jù)．特別在流固耦合的應用方面，需要流體和固體中各自網(wǎng)格離散要能夠精確地捕捉相互物理作用，其中涉及多尺度和復雜的多模態(tài)耦合動力學．這種非匹配網(wǎng)格的要求在機翼、深海立管、系泊管線的流固耦合數(shù)值模擬應用中非常常見．非守恒的數(shù)據(jù)傳輸以及局部不精確的插值和投影將導致數(shù)值不穩(wěn)定性以及較差的模擬結(jié)果．特別是當激勵頻率接近彈性結(jié)構(gòu)的固有頻率時尤為明顯．因此，準確的流固耦合模擬需要對存在非匹配網(wǎng)格的耦合界面的邊界條件進行精確而守恒的處理．

對于基于Eulerian-Lagrangian流固耦合系統(tǒng)，其數(shù)值模擬存在兩種主要方法，即整體模擬和分區(qū)域模擬[1-3]．在整體模擬方法中，流體和結(jié)構(gòu)方程使用完全耦合的方式一起聯(lián)立求解．整體求解方式穩(wěn)定度較高，即使存在非常強的附加質(zhì)量效應時也如此，但是這種方法缺少靈活性和可塑性．另一方面，分區(qū)域算法則是按照一定的次序分別求解流體和結(jié)構(gòu)兩個子區(qū)域的控制方程[4-5]，進而可以促進現(xiàn)有的流體和固體數(shù)值求解器之間進行合適的選擇搭配，進行空間和時間的離散，但也間接造成網(wǎng)格劃分的非匹配度的提高．

在一個典型的三維流固耦合問題分區(qū)域數(shù)值模擬中，耦合界面上的面網(wǎng)格拓撲一般并不相同．這意味著其非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)也并不相同，由于離散化要求，它們的幾何坐標也可能不一致．這種非匹配網(wǎng)格和相關的數(shù)據(jù)傳輸問題也存在于其他情形中，如自適應網(wǎng)格劃分和多重網(wǎng)格法[6-7]．

非匹配網(wǎng)格的數(shù)據(jù)插值和投影，存在各種守恒與非守恒的方式．在流固耦合數(shù)值模擬中，數(shù)據(jù)在傳輸后必須保證足夠的準確和物理守恒，尤其是在非定常模擬中尤為重要．這是由于誤差可能隨著迭代和非定常模擬不斷累積，如果插值算法足夠準確且守恒，相比非守恒且不夠準確的格式，誤差便不會放大且會加快數(shù)值模擬的收斂．但是國內(nèi)和國外均缺少相關?研究．

Common-Refinement方法是一種能夠保證數(shù)據(jù)傳輸?shù)臏蚀_性和計算模擬的守恒性的，針對非匹配網(wǎng)格一種網(wǎng)格數(shù)據(jù)結(jié)構(gòu)的構(gòu)造方法．由于實現(xiàn)難度較大與計算技術問題，之前沒有得到足夠的重視．在當今計算資源與計算技術飛速發(fā)展的今天，Common-Refinement方法逐漸受到了流固耦合研究人員的關注[8-12]．

本文的主要目標是三維Common-Refinement方法的并行數(shù)值實現(xiàn)，以及對非匹配網(wǎng)格的誤差進行系統(tǒng)的量化和分析，并探討其在三維非定常流固耦合數(shù)值計算中的應用，目的是減少匹配網(wǎng)格的必要性，進而設計一個更具有一般性的針對三維流固耦合交界曲面上的非匹配網(wǎng)格節(jié)點的投影方案．本文結(jié)合不同的拓撲結(jié)構(gòu)非匹配網(wǎng)格處理耦合界面，對非匹配網(wǎng)格的空間誤差與時間誤差的累積進行系統(tǒng)的分析，并對三維Common-Refinement的并行效率進行了系統(tǒng)的統(tǒng)計與分析．基于此，本文進一步基于標準的圓柱體-彈性板問題探討了三維Common-Refinement方法在非定常流固耦合數(shù)值計算中的全局誤差，系統(tǒng)地劃分了多套非匹配網(wǎng)格并進行了網(wǎng)格收斂性分析與全局誤差分析，計算結(jié)果表明Common-Refinement方法在三維非定常流固耦合數(shù)值計算中并不會造成誤差的累積，同時能夠保證全局二階精度，具有非常好的穩(wěn)定性與準確性．

1?三維Common-Refinement方法與耦合界面邊界條件的數(shù)學描述

本文首先給出基于非匹配網(wǎng)格的三維流固耦合問題的邊界條件，以及Common-Refinement方法的數(shù)學描述．在三維流固耦合計算中，流體和結(jié)構(gòu)方程由速度和牽引力沿流體—結(jié)構(gòu)耦合界面的連續(xù)性來進行數(shù)據(jù)交換，進而互相滿足各自區(qū)域的邊界條件．由于本文采用了區(qū)域分解策略，本文首先各自生成了流體區(qū)域和結(jié)構(gòu)域的三維有限元非結(jié)構(gòu)網(wǎng)格，且兩者網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)互相獨立．流體作用在結(jié)構(gòu)表面上的力作為表面牽引力矢量來處理，同時結(jié)構(gòu)體的位移給流體域的變形提供了邊界條件．在本文中，耦合界面的插值多項式的次數(shù)在流體域與結(jié)構(gòu)域保持一致．

?(1)

?(2)

?(3)

?(4)

?(5)

?(6)

當質(zhì)量矩陣的構(gòu)造只需要結(jié)構(gòu)一側(cè)的形函數(shù)時，載荷矢量積分由流體和固體的形函數(shù)共同組成．對于匹配的網(wǎng)格，這不會造成任何問題．然而對于非匹配網(wǎng)格，形函數(shù)的不一致將導致兩側(cè)的不連續(xù)性．

Common-Refinement方法是一類重要且特殊的網(wǎng)格數(shù)據(jù)結(jié)構(gòu)，兩套網(wǎng)格間Common-Refinement曲面和網(wǎng)格是由對流體和固體子域的輸入邊界網(wǎng)格同時進行細分的多邊形來構(gòu)成．借助有限元形式，流體和固體界面網(wǎng)格的空間坐標可以寫為

?(7)

?(8)

在本文中，Common-Refinement方法生成的網(wǎng)格拓撲由兩側(cè)輸入的界面網(wǎng)格的交線來定義，或者稱作子單元．三維情況下子單元的構(gòu)成參見圖1．

圖1?三維Common-Refinement方法示意

?(9)

2?三維Common-Refinement方法的靜態(tài)數(shù)值傳輸誤差分析

圖2?誤差分析圓柱體算例示意

Fig.2 Schematic of error andysis based cylinder case

這里定義單元數(shù)量為

?(10)

?(11)

?(12)

表1?非匹配網(wǎng)格相對誤差計算結(jié)果

Tab.1?Relative error results of non-matching mesh

圖4?四邊形非匹配網(wǎng)格誤差分析結(jié)果

為了更進一步研究與量化Common-Refinement方法的數(shù)據(jù)傳輸?shù)目臻g誤差與提高Common-Refinement方法的應用范圍，本文研究了兩側(cè)不同幾何拓撲的非匹配網(wǎng)格數(shù)據(jù)靜態(tài)傳輸?shù)目臻g誤差．數(shù)據(jù)或載荷傳輸?shù)男问胶椭跋嗤?，唯一不同的是傳輸一?cè)，也就是流體一側(cè)的網(wǎng)格設置為三角形非結(jié)構(gòu)網(wǎng)格，而數(shù)據(jù)接收一側(cè)，也就是結(jié)構(gòu)體一側(cè)的網(wǎng)格為四邊形非匹配網(wǎng)格．此時非匹配網(wǎng)格見圖5．

圖5?典型的非匹配網(wǎng)格Ⅱ

Fig.5 Typical non-matching meshⅡ of solid(blue line)and fluid(red line)

流體一側(cè)傳輸至結(jié)構(gòu)一側(cè)的載荷形式同樣定義為式(11)，空間誤差定義為式(12)，此時網(wǎng)格尺度范圍為0.04,m～0.20,m，4套不同非匹配網(wǎng)格的空間誤差結(jié)果見圖6．

圖6?三角形-四邊形非匹配網(wǎng)格誤差分析結(jié)果

Fig.6 Spatial error analysis results for tri-quad non-matching meshes

從圖6能夠看出，從三角形-四邊形非匹配非結(jié)構(gòu)網(wǎng)格的設置中，Common-Refinement方法同樣能夠保持數(shù)據(jù)傳輸?shù)目臻g二階精度，在證明了Common-Refinement方法精度的同時也提高了此方法的應用范圍．

3?Common-Refinement方法在非定常計算中的誤差累積與并行效率分析

在典型的非定常流固耦合模擬中，流體和固體的表面總是互相完整重合的．流體的載荷作用在到結(jié)構(gòu)體表面上，而結(jié)構(gòu)的位移作用在流體表面，進而在每一個時間步重復這樣的數(shù)據(jù)傳輸．本文分析了Common-Refinement方法在非定常流固耦合計算中重復傳遞數(shù)據(jù)時引入的累積誤差．在非定常流固耦合計算中，本文采用連續(xù)Galerkin(CG)有限元離散化方法對線彈性體結(jié)構(gòu)模型進行離散，采用Petrov-Galerkin有限元離散求解不可壓縮流動．

圖7?Common-Refinement方法的瞬時誤差分析問題示意

圖8?瞬時誤差分析用匹配網(wǎng)格剖面

圖9?瞬時誤差分析用非匹配網(wǎng)格剖面

圖10?Common-Refinement網(wǎng)格的子單元示意

?(13)

?(14)

?(15)

表2?特征響應的瞬時誤差計算結(jié)果

Tab.2?Transient error results of characteristic response

從瞬時誤差的結(jié)果可以看出，每組響應之間只存在非常小的差別．結(jié)果表明，即使在網(wǎng)格匹配度不高的情況下，本文提出的三維Common-Refinement方法能夠在非定常計算中保證誤差不會隨時間的增加而增大，在數(shù)值上能夠保證數(shù)值計算過程的準確性與穩(wěn)定性．

本文進而分析了Common-Refinement方法在流固耦合非定常計算過程中的計算效率．在每一個時間步，流固耦合求解器每次迭代所花費的時間由3個主要部分組成，即流體求解器、結(jié)構(gòu)求解器和Common-Refinement方法的數(shù)據(jù)傳輸．本文通過統(tǒng)計每次迭代所花費的時間來評估每個組成部分的性能．在不失一般性的情況下，本文繼續(xù)使用本節(jié)的中彈性圓柱體算例進行性能評估．此時流體域被離散為具有35,984節(jié)點的33,450個六面體單元，而結(jié)構(gòu)域被離散成具有2,314節(jié)點的1,800個六面體單元．流固耦合界面由具有832節(jié)點的800個四邊形單元組成．數(shù)值計算使用雙核Intel(R)Xeon(R)E5-2630L v2@2.40,GHz的CPU．每個組成部分的時間占用是在經(jīng)過500個時間步長后統(tǒng)計的．表3總結(jié)了數(shù)值計算中每個部分的平均時間占用．平均而言，大約在每一個時間步需要用4次迭代滿足10-5量級的收斂容差．

表3?數(shù)值模擬中時間占用結(jié)果

Tab.3 Elapsed time for each component with numerical simulation

從表3中能夠看出，大部分的時間都花在流體求解和結(jié)構(gòu)求解計算，而Common-Refinement方法的數(shù)據(jù)傳輸并沒有對總計算時間有著明顯貢獻．與一般計算成本相比，Common-Refinement方法數(shù)據(jù)的計算成本是微不足道的．

基于此，本文進一步研究了Common-Refinement方法的并行實現(xiàn)，并評估了其并行性能．并行效率通過使用不同數(shù)量的處理器在流體表面與結(jié)構(gòu)表面?zhèn)鬏數(shù)臄?shù)據(jù)所占用的時間來統(tǒng)計．此時流體和結(jié)構(gòu)表面設置為單位長度100×100，表面離散為5040×1000個四邊形單元，共有5041×1001個節(jié)點．在不失一般性的情況下，將任意數(shù)據(jù)分配給流體表面上的每個節(jié)點，并通過Common-Refinement方法將其傳輸?shù)浇Y(jié)構(gòu)表面．在并行設置中，網(wǎng)格將根據(jù)數(shù)據(jù)傳輸使用的處理器的數(shù)量進行劃分，由MPI實現(xiàn)并行數(shù)據(jù)傳輸過程．

圖11?Common-Refinement方法的并行效率

圖11顯示了在不同數(shù)量的處理器的并行數(shù)據(jù)傳輸?shù)乃璧臅r間，處理器數(shù)目從2～16不等．圖中最佳擬合線的梯度為-1，表示時間與處理器數(shù)成反比．這證明了本文的Common-Refinement方法良好的并行性能．表明此方法可以擴展到更多處理器，提高流固耦合非定常數(shù)值模擬的并行性能．

4?Common-Refinement方法在流固耦合數(shù)值計算中全局誤差分析

為了進一步評估Common-Refinement方法在非定常流固耦合計算中的全局誤差，本文基于標準算例進行誤差測試．本文采用連續(xù)Galerkin(CG)有限元離散化方法對非線性超彈性體結(jié)構(gòu)模型進行離散，采用Petrov-Galerkin有限元離散求解不可壓縮流動．流固耦合的計算流程見圖12．

圖12?3D流固耦合計算流程

在數(shù)值實現(xiàn)過程中，本文采用OpenMP結(jié)合MPI實現(xiàn)并行計算流程，使用C語言編寫了完整的三維Common-Refinement方法與并行三維流固耦合求解代碼．流體域的載荷通過Common-Refinement網(wǎng)格傳輸給結(jié)構(gòu)域，作為超彈性體的邊界條件，超彈性體在當前載荷下產(chǎn)生的位移同樣通過Common-Refinement網(wǎng)格傳輸給流體域，作為ALE網(wǎng)格求解和N-S方程求解的邊界條件，期間使用非線性力的迭代方式對當前流動載荷逐次修正．

圖13?圓柱體-彈性板問題示意

表4?圓柱體-彈性板問題相關參數(shù)

Tab.4?Parameters for cylinder-foil problem

圓柱直徑/m平均流速/(m·s-1)流體密度/(kg·m-3)彈性板厚度/m彈性板長度/m彈性板密度/(kg·m-3)楊氏模量/Pa雷諾數(shù)密度比泊松比 0.11.01,0000.020.351,0005.6×1052001.00.4

???(16)

(17)

本文使用有限元網(wǎng)格對流體和固體區(qū)域進行分解．有一層邊界層網(wǎng)格劃分在結(jié)構(gòu)體四周．3套不同密度的網(wǎng)格用來進行網(wǎng)格收斂性的研究，同時進行全局誤差分析研究．圖14中給出了圓柱體-彈性板問題網(wǎng)格劃分示意．

圖14?圓柱體-彈性板問題網(wǎng)格劃分示意

Fig.14 Schematic of mesh configuration of cylinder-foil problem

在網(wǎng)格收斂性分析的同時，為了檢測全局誤差，使用彈性板的共振頻率定義全局誤差為

?(18)

式中：為彈性板方向的共振頻率；對應于網(wǎng)格M3的計算結(jié)果．表5給出了網(wǎng)格收斂性的計算結(jié)果，從表中可以看出M3網(wǎng)格的收斂性達到最好，其前端的最大位移和響應頻率和標準算例的結(jié)果吻合良好．

表5?網(wǎng)格收斂性計算結(jié)果

Tab.5?Calculation results of mesh convergence

網(wǎng)格流體單元結(jié)構(gòu)單元 M17,1711362.6832.945.15 M211,7705442.5531.895.34 M325,8821,2802.4532.895.42 FSI-III9,2161,0802.6835.345.30

圖16給出了全局誤差隨著網(wǎng)格數(shù)量變化的結(jié)果．可見在網(wǎng)格數(shù)量變動的情況下，全局誤差較好的保持了二階精度，數(shù)值計算的結(jié)果證明了Common-Refinement方法能夠在流固耦合非定常計算中保證全局二階精度．

圖15?彈性板前端位移時歷曲線

圖16?網(wǎng)格收斂性結(jié)果

5?結(jié)?論

本文實現(xiàn)了基于Common-Refinement方法和有限元法的非匹配網(wǎng)格之間的空間插值方式，自主開發(fā)了完整的數(shù)值求解器，并系統(tǒng)地進行了誤差分析和并行性能分析，且將其用于不可壓縮流體和非線性超彈性體的非定常流固耦合問題．通過本文的研究能夠得出以下結(jié)論：

(1) Common-Refinement方法在靜態(tài)數(shù)據(jù)傳輸中能夠非常好的保證空間插值的二階精度，同時可以任意劃分和使用不同拓撲形式的非匹配網(wǎng)格，具有足夠的靈活性和穩(wěn)健性．

(2) Common-Refinement方法在非定常流固耦合計算中表現(xiàn)穩(wěn)健，即使網(wǎng)格匹配度不高的情況下同樣能夠保證數(shù)值計算的精度，且空間插值誤差不會隨著時間的推移不斷積累導致計算發(fā)散．

(3) Common-Refinement方法具有非常好的并行效率，能夠應用于大規(guī)模非定常流固耦合計算中的耦合界面數(shù)據(jù)傳輸問題，且不會占用過多的計算資源.

(4) 本文通過三維低密度比的標準算例的一系列詳細而系統(tǒng)的數(shù)值實驗，證明了Common-Refinement方法能夠保證流固耦合數(shù)值計算中的全局二階精度，具有非常好的準確性與可靠,．

［1］ Hron J，Turek S.[M]. Germany：Springer，2006.

［2］ Felippa C，Park K. Staggered transient analysis procedures for coupled mechanical systems：Formulation[J].，1980，24：61-111.

［3］ Gurugubelli P，Jaiman R. Self-induced flapping dynamics of a flexible inverted foil in a uniform flow[J].，2015，78：657-694.

［4］ Jaiman R，Sen S，Gurugubelli P. A fully implicit combined field scheme for freely vibrating square cylinders with sharp and rounded corners[J].，2015，112：1-18.

［5］ Jaiman R，Geubelle P，Loth E，et al. Combined interface condition method for unsteady fluid-structure interaction[J].，2011，200：27-39.

［6］ Boer A，Zuijlen A，Bijl H. Review of coupling methods for non-matching meshes[J].，2007，8：1515-1525.

［7］ Boer A，Zuijlen A，Bijl H. Comparison of conservative and consistent approaches for the coupling of non-matching meshes[J].

，2009，11：415-455.

［8］ Jaiman R，Jiao X，Geubelle X，et al. Conservative load transfer along curved fluid-solid interface with nonmatching meshes[J].，2006，218：372-397.

［9］ Jaiman R，Jiao X，Geubelle X，et al. Assessment of conservative load transfer on fluid-solid interface with nonmatching meshes[J].，2005，64：2014-2038.

［10］ Jiao X，Heath M. Common-Refinement based data transfer between nonmatching meshes in multi-physics simulations[J].，2004，61：2402-2427.

［11］ Jiao X，Heath M. Overlaying surface meshes，partⅠ：Algorithms[J].，2004，14(6)：379-402.

［12］ Jiao X，Heath M. Overlaying surface meshes，part Ⅱ：Topology preservation and feature matching[J].，2004，14(6)：403-419.

［13］ Turek S，Hron J.，：[M]. Germany：Springer，2006.

(責任編輯：孫立華)

Error Analysis and Parallel Performance of 3D Common-Refinement Method for Unsteady Numerical Calculation of Fluid-Structure Interaction

Li Yulong1，Chu Nanyang1，Hu Zhan1，Yao Peng1，Deng Rui2

(1. College of Marine Science，Sun Yat-sen University，Guangzhou 119077，China；2. College of Ocean Engineering and Technology，Sun Yat-sen University，Guangzhou 119077，China)

A three-dimensional(3D)Common-Refinement method for non-matching meshes between discrete non-overlapping sub-domains of incompressible fluid and nonlinear elastic structure was presented in this paper. To begin with，the accuracy of Common-Refinement method(CRM)to satisfy traction equilibrium condition along the fluid-elastic interface with non-matching meshes was investigated. Then，the accuracy of CRM against the matching grid solution by varying grid mismatch between the fluid and solid meshes with different topological structures was systematically assessed. The second-order accuracy of CRM was acquired through uniform refinements of fluid and solid meshes along the interface. Meanwhile，error accumulation during unsteady fluid-structure interaction simulation as well as parallel performance of Common-Refinement method was investigated. In the end，a 3D benchmark problem of a cantilevered hyperelastic plate behind a circular bluff body was solved to verify the accuracy of coupled solutions with respect to the solution available in the literature. Numerical results show that CRM can maintain second-order accuracy for unsteady 3D fluid structure interaction numerical simulations.

fluid-structure interaction；error analysis；Common-Refinement；non-matching meshes

the National Key Research and Development Program of China(No.2016YFC0402601)，the National Natural Science Foundation of China(No.51609269，No.51679053，No.51209051，and No.51209048)，and the Foundation of Pre-Research of Military Equipment of the Chinese People's Liberation Army(No.6140241010116CB01005).

1234

A

0493-2137(2018)11-1201-09

2018-03-12；

2018-05-03.

李裕龍（1984—??），男，副研究員，liylong7@mail.sysu.edu.cn.

姚?鵬，yaop3@mail.sysu.edu.cn.

國家重點研發(fā)計劃資助項目(2016YFC0402601)；國家自然科學基金資助項目(51609269，51679053，51209051，51209048)；海軍裝備預研基金資助項目(6140241010116CB01005).

10.11784/tdxbz201803039