基于分區交錯算法的高層建筑風振數值模擬平臺及實例

鄭德乾，顧 明，張愛社

(1.河南工業大學 土木建筑學院，鄭州 450001;2.同濟大學 土木工程防災國家重點實驗室，上海 200092;3.山東建筑大學 土木工程學院，濟南 250101)

隨著建筑高度的增加，以及輕質、高強材料在土木工程領域的應用，結構的風敏感性逐漸增強。在風荷載的作用下，高層建筑結構會發生彎曲或扭轉變形;在結構變形的同時，其周圍流場也會因為結構外形的改變而變化，風與結構的相互作用問題不是風的流動問題和結構運動問題的簡單迭加，二者之間存在著相互耦合的氣動彈性現象。以往對于風與結構相互作用問題的研究多采用風洞試驗方法［1－3］，隨著計算機技術的發展，數值模擬方法逐漸應用于風與結構的相互作用問題求解中［4－11］。

本文討論的風與結構的相互作用問題，其耦合作用僅發生結構表面上，由兩相耦合面上的動力學平衡及運動學協調條件來引入方程上的耦合［4］。直接求解的強耦合法和基于分區交錯算法的弱耦合法［4］是數值求解此類流固耦合 (Fluid Structure Interaction，FSI)問題的較常用方法。直接求解的強耦合法將流體域、結構域和耦合作用構造在同一個控制方程中，在同一時間步內同時求解所有變量［4］，這種方法概念清晰，但其分析求解不便于采用現有成熟程序及專業軟件?；诜謪^交錯算法的弱耦合法，在每一個時間步內分別依次對CFD(Computational Fluid Dynamics)方程和CSD(Computational Structural Dynamics)方程求解，通過中介交換固體域和流體域的計算結果數據，從而實現耦合求解［4］。在基于分區交錯算法的弱耦合法中，CFD和CSD求解過程相互獨立，便于采用各領域現有程序及軟件實現耦合系統的模塊化求解，并且數值模擬結果精度也相對較高，是目前使用相對較廣泛的流固耦合問題求解方法［5－11］。

FSI計算中CFD流體域及其邊界條件會隨時間發生變化，一般可采用達朗貝爾原理［10］或 ALE方法(Arbitrary Lagrangian-Eulerian Method)［11］等手段修改流體域Navier-Stokes(N-S)方程，實現氣動彈性問題的求解計算。ALE方法能夠較好地解決流體網格的歐拉描述與結構網格的拉格朗日描述的不一致問題，且在考慮網格移動速度的同時能夠對運動耦合界面準確描述，已被廣泛應用于求解工程領域 FSI問題［5－9］。

本文基于ALE方法和弱耦合分區交錯算法，通過對商業流體軟件Fluent和結構計算軟件Ansys的二次開發，搭建了一個較通用的氣動彈性問題數值模擬計算平臺。采用該平臺，對一寬高比為1∶6的方形截面高層建筑風振氣彈響應進行了數值模擬，通過與氣彈模型風洞試驗及文獻模擬結果的對比，以及基于有、無考慮氣動彈性時結構位移響應比較的氣動阻尼影響分析，驗證了本文方法的有效性。

1 數值模擬方法

1.1 控制方程

FSI問題求解中涉及的控制方程主要有CFD流體控制方程、CSD結構運動方程，以及CFDCSD耦合界面相容條件。基于ALE［11］描述的CFD流體域N-S控制方程為:

式中:ρ為流體密度，ν為運動粘性系數，u為流體速度矢量，ug為網格運動速度。在任一運動邊界控制體積V上，式(1)的積分形式可寫為:

式中:S表示控制體積V的邊界，(u－ug)為非定常條件下的對流速度，n為界面法向量。為保證守恒，每個控制體V上應滿足網格守恒律:

將式(3)對網格速度ug的積分項以通過控制體各面的質量通量(即由控制體各面的運動產生的網格通量)來代替，則式(2)包含網格速度的對流項可被離散為［6，8］:

式中:Φ為廣義標量，對于連續方程，Φ=1;對于動量方程，Φ ={u，v，w};c={w，e，s，n，b，t}。經過以上變換，式(2)的求解將基于相對通量(即通過控制體的由流體運動造成的通量減去網格通量)［6，8］，這樣即可在網格運動速度ug為未知情況下進行方程(2)的求解。

數值求解FSI問題時，由于CFD和CSD對計算網格要求的不同，導致耦合界面上的網格不匹配［4］，為保證流體－結構耦合系統的一致性，耦合界面的數據傳遞需要滿足力學平衡條件和運動學連續的相容條件，即［6－7］

式(6)為耦合界面力的平衡條件，其中σS是結構應力張量，σF是流體粘性應力張量，p是流體壓力。式(7)為耦合界面的位移相容條件，其中uS和 uF分別指耦合界面上結構位移矢量和流體域的位移矢量。

1.2 數值模擬方法

基于弱耦合分區交錯算法，本文搭建的高層建筑風振數值模擬平臺由四部分組成，分別為:CFD計算模塊、CSD計算模塊、網格更新模塊和數據傳遞模塊。圖1a為計算平臺構成及各模塊的軟件實現方法。

(1)CFD計算:基于ALE法［11］的流體控制方程的求解采用有限體積法。基于商業流體軟件Fluent，采用Scheme語言編程實現流體計算的參數化、程序化執行。為獲得結構表面的平均和脈動風荷載，流體計算采用基于空間平均的大渦模擬 (Large Eddy Simulation，LES)方法。采用動態亞格子模型求解亞格子應力項。LES入流脈動的合成采用文獻［12］方法。由于LES方法計算量相對較大，為提高CFD計算效率，流體計算采用并行算法。針對Fluent流體并行計算中流體網格分塊和并行計算各進程的特點［13］，采用UDF編程實現流體并行計算時結構表面風荷載的準確提取。本文FSI計算時間步長與CSD計算一致，而CFD大渦模擬時間步長相對較小，為避免時間步長不匹配，CFD計算中采用子迭代技術。

圖1 氣動彈性數值模擬平臺構成、軟件實現以及求解流程示意圖Fig.1 Sketches of composition and solution procedure of present numerical simulation platform

(2)CSD計算:結構運動方程的求解采用Newmark法［14］。基于結構有限元計算軟件 Ansys，通過APDL(Ansys Parametric Design Language)語言實現結構計算的參數化編程。FSI計算的每個時間步，執行CSD計算的Ansys均需要退出，為此采用Ansys重啟計算技術來滿足結構瞬態響應計算需要。

(3)數據傳遞:在流固耦合界面上，因CFD和CSD計算對網格密度的要求不同，存在非匹配網格的搜索配對和相應的數據傳遞問題。本文采用分區搜索算法進行耦合界面網格配對，即對CFD計算中定義的結構各表面分別進行搜索配對，以提高網格匹配效率?；谑睾悴逯捣ǎ?］，采用有限元形函數插值實現流固耦合界面上數據傳遞，來確保CFD網格傳遞給CSD網格的荷載平衡，以及CSD網格傳遞給CFD網格的位移相容。通過UDF編程實現耦合界面非匹配網格的分區搜索配對及數據傳遞。

(4)網格更新:基于Fluent動網格技術［13］，通過UDF的DEFINE_GRID_MOTION宏編程實現耦合界面的運動變形，采用彈簧光滑法進行流體網格更新。

(5)主進程:基于分區交錯算法，采用VC++編制主進程實現上述四個模塊之間的調用(圖1(a))，本文FSI數值模擬平臺求解流程如圖1(b)所示。

需要指出的是，當采用本文FSI平臺計算時，若每個FSI時間步均不進行CFD模塊流場網格的更新，還可實現不考慮氣動彈性的結構風致振動數值模擬。

2 算例分析

2.1 計算對象及參數

文獻［1］采用一通用超高層建筑結構的多自由度結構模型［2］(如圖2所示)，進行了B類1/500縮尺比風場的氣彈模型風洞試驗，模型的寬度和高度分別為0.1 m 和0.6 m，模型總質量 1.64 kg，實測一階阻尼比0.0075，實測前二階基頻均為 20.5 Hz。文獻［8］基于RANS(Reynolds Averaged Navier-Stockes)方法求解流場，對該模型進行了氣彈數值模擬計算。

圖2 高層建筑結構多自由度氣動彈性模型風洞試驗結構模型［2］Fig.2 Multi-degree-of-freedom aeroelastic tall building model［2］

采用文獻［1］相關參數，基于搭建的FSI數值模擬平臺，本文進行了圖2所示結構模型0°風向角 (來流垂直于結構迎風面)情況的氣彈響應數值模擬。CFD計算采用與風洞試驗相同的縮尺比，網格劃分如圖3(a)所示，圖中 B=D=0.1 m，H=0.6 m，網格總數896000，近壁面最小網格尺度為1/2000D，對應壁面y+≤5。CSD建模中分別采用梁單元和殼單元模擬結構的剛度分布和外形，在Ansys中建立了結構的有限元模型，如圖3(c)所示。CSD模型兩個主軸方向的基頻均為20 Hz，非常接近試驗模型相應值。耦合界面上的CFDCSD網格分別如圖3(b)和圖3(c)所示，共計6800個CFD面網格和2500個CSD面網格。

本文FSI計算中，先進行結構靜止不動時的CFD大渦模擬計算，待流場充分發展后，釋放結構運動，接著進行FSI數值模擬計算。表1為本文FSI計算工況，表中折減風速U*=UH/(f1D)，其中UH為模型高度H處來流平均風速，f1為結構一階頻率。FSI計算中不改變來流風速，通過改變結構固有頻率實現不同折減風速變化。表1中所有FSI計算工況的結構質量均為1.64 kg，與文獻［1］風洞試驗和文獻［8］氣彈數值模擬保持一致。為分析氣動阻尼對結構風振響應的影響，本文還進行了不考慮氣動彈性時結構風振響應數值模擬。

圖3 CFD、CSD模型及網格示意圖Fig.3 CFD meshes ＆ boundary conditions，and sketches of CSD model

表1 高層建筑風振響應數值模擬工況Tab.1 Case details for study on wind-structure interaction of a tall building

2.2 氣彈數值模擬結果比較分析

FSI計算所得高層建筑模型頂部節點位移響應根方差與文獻［1］氣彈模型風洞試驗和文獻［8］數值模擬結果的比較如圖4所示。由圖可見:① 本文數值模擬所得結構的順風向響應與模型風洞試驗結果總體上吻合較好，只是數值模擬結果數值偏小;② 對于結構的橫風向響應，本文氣彈數值模擬結果與風洞試驗結果也具有較好的一致性，并且精度明顯高于基于RANS的氣彈數值模擬結果。

2.3 有、無考慮氣動彈性數值模擬結果比較分析

數值模擬所得有、無考慮氣動彈性時，不同折減風速下結構頂部節點位移響應時程的比較如圖5所示，圖6為位移響應根方差隨折減風速的變化比較。由圖可見:① 折減風速較小時，有、無考慮氣動彈性時結構順風向位移響應振幅及根方差比較接近，不考慮氣動彈性結果稍大;隨著折減風速U*的增加，結構順風向位移響應振幅均有所增大，不考慮氣動彈性情況的振幅增大相對顯著，由圖6(a)可明顯看出兩種情況下位移響應根方差數值上的差異隨U*的增加逐漸增大。② 折減風速U*＜12.0時，有、無考慮氣動彈性時結構橫風向位移響應的比較結果與順風向類似;但折減風速增大至U*=12.0，橫風向渦激共振發生時，結構橫風向氣彈位移響應振幅及其根方差值明顯大于不考慮氣動彈性結果，此時氣動阻尼對結構風振響應的影響不可忽略。

圖4 模型頂部節點氣彈位移響應根方差的比較Fig.4 Comparison results of r.m.s values of aeroelastic displacement responses of points at the top of the building model

圖5 有、無考慮氣動彈性時模型頂部節點位移響應時程數值模擬結果比較Fig.5 Comparisons of time histories of displacement responses of points at the top of the building model between present simulation results with and without considering FSI

圖6 有、無考慮氣動彈性時模型頂部節點位移響應根方差數值模擬結果比較Fig.6 Comparisons of displacement responses of points at the top of the building model between present simulation results with and without considering FSI

由以上分析結果可推斷:本文計算折減風速范圍內，順風向氣動阻尼均為正值;而橫風向氣動阻尼則隨著折減風速的增加，由正值逐漸變為負值。這與文獻［15－16］的試驗研究結果基本一致。

3 結論

基于弱耦合分區交錯算法，通過對商業流體軟件Fluent和結構計算軟件Ansys進行二次開發，搭建了風與結構相互作用氣動彈性問題數值模擬計算平臺。采用本文方法，對三維方形截面高層建筑風振響應進行了數值模擬分析。結果表明，基于本文方法的氣彈數值模擬結果與試驗結果基本一致，并且精度上高于已有文獻數值模擬結果;有、無考慮氣動彈性時結構風振響應數值模擬結果的比較表明，折減風速較小時，氣動阻尼的存在有利于減弱結構的順、橫風向風振響應;但當折減風速增大至出現橫風向渦激共振現象時，氣動阻尼的存在明顯增強了結構的橫風向振動。本文方法可以較精確地數值求解高層建筑結構的風振氣彈響應問題。

［1］黃鵬.高層建筑風致干擾效應研究［D］.上海:同濟大學，2001.

［2］全涌，顧明.超高層建筑通用氣動彈性模型設計［J］.同濟大學學報，2001，29(1):122－126.QUAN Yong，GU Ming.Design of super-high rise buildings’global aeroelastic model［J］.Journal of Tongji University，2001，29(1):122－126.

［3］章李剛，樓文娟.復雜外形超高層建筑結構三維風致響應分析［J］.振動與沖擊，2013，32(24):38 －43.ZHANG Li-gang，LOU Wen-juan.Three dimensional windinduced responses of super-tall buildings with irregular geometric shapes［J］.Journal of Vibration and Shock，2013，32(24):38－43.

［4］錢若軍，董石麟，袁行飛.流固耦合理論研究進展［J］.空間結構，2008，14(1):3－15.QIAN Ruo-jun，DONG Shi-lin，YUAN Xing-fei.Advances in research on fluid-structure interaction theory［J］.Spatial Structures，2008，14(1):3－15.

［5］Braun A L，Awruch A M.Aerodynamic and aeroelastic analyses on the CAARC standard tall building model using numerical simulation［J］.Computers ＆ Structures，2009，87(9－10):564－581.

［6］Glück M.Computation of wind-induced vibrations of flexible shells and membranous structures［J］.Journal of Fluids and Structures，2003，17(5):739 －765.

［7］Farhat C，LeTallec M L P.Load and motion transfer algorithms for fluid/structure interaction problems with nonmatching discrete interfaces:Momentum and energy conservation，optimal discretization and application to aeroelasticity［J］.Comput.Methods Appl.Mech.Engrg.，1998，157:95 －114.

［8］方平治.典型建筑結構氣動彈性問題的數值模擬研究［D］.上海:同濟大學，2007.

［9］Michalski A，Kermel P D，Haug E，et al.Validation of the computational fluid-structure interaction simulation at realscale tests of a flexible 29m umbrella in natural wind flow［J］. Journal of Wind Engineering and Industrial Aerodynamics，2011，99(4):400 －413.

［10］Tutar M，Holdo A E.Large eddy simulation of a smooth circular cylinder oscillating normal to a uniform flow ［J］.Journal of Fluids Engineering，2000，122:694－702.

［11］Hirt C W，Amsden A A，Cook J L.An arbitrary Lagrangian-Eulerian computing method for all flow speeds［J］.Journal of Computational Physics，1997，135:203－216.

［12］Zheng D Q，Zhang A S，Gu M.Improvement of inflow boundary condition in Large Eddy Simulation of flow around a tall building［J］.Engineering Applications of Computational Fluid Mechanics，2012，6(4):646－660.

［13］Fluent.Inc.Fluent 6 Documentation，2005.

［14］Bathe K J.Finite element procedures［D］.Prentice-Hall，Englewood Cliffs，1996.

［15］Gu M，Quan Y.Across-wind loads of typical tall buildings［J］.Journal of Wind Engineering and Industrial Aerodynamics，2004，92(13):1147 －1165.

［16］Quan Y，Gu M，Tamura Y.Experimental evaluation of aerodynamic damping of square super high-rise buildings［J］.Wind and Structures，2005，8(5):309 －324.