鄭德乾，顧 明，張愛社

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

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

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

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

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

1 數(shù)值模擬方法

1.1 控制方程

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

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

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

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

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

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

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

1.2 數(shù)值模擬方法

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

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

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

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

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

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

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

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

2 算例分析

2.1 計算對象及參數(shù)

文獻［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)方法求解流場，對該模型進行了氣彈數(shù)值模擬計算。

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

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

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

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

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

2.2 氣彈數(shù)值模擬結果比較分析

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

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

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

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

圖5 有、無考慮氣動彈性時模型頂部節(jié)點位移響應時程數(shù)值模擬結果比較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 有、無考慮氣動彈性時模型頂部節(jié)點位移響應根方差數(shù)值模擬結果比較Fig.6 Comparisons of displacement responses of points at the top of the building model between present simulation results with and without considering FSI

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

3 結論

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

［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］方平治.典型建筑結構氣動彈性問題的數(shù)值模擬研究［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.