張其林 閆雁軍 李 晗
(同濟大學建筑工程系, 上海 200092)
張拉結構具有質量輕、剛度柔、阻尼小等特點,風荷載作用下,維護結構易發(fā)生振動,可能引發(fā)局部撕裂或破壞,甚至導致主體結構失效.對風致振動的研究已有很多,大部分集中在大跨橋梁和高層建筑[1-6].對于大跨張拉結構的研究工作相對滯后,研究成果較少,目前尚沒有成熟的抗風設計方法.本文采用風洞試驗、FSI數(shù)值模擬方法和基于隨機振動理論的頻域方法,對結構風致振動響應結果進行了分析比較,驗證了理論方法的可靠性.
流固耦合(FSI)是研究變形固體在流場作用下的各種行為以及固體位形與流場分布之間相互作用的一門交叉學科.流固耦合的一個重要特征是流體與固體兩相介質的相互作用,固體在流體作用下發(fā)生變形、運動;固體變形或運動反過來又影響流體流動,從而改變流體載荷的分布和大?。?/p>
耦合問題的整體方程可表示為
(1)
式中,Ms,Mf分別為結構質量矩陣和荷載質量矩陣;Mfs,Msf分別為流固耦合界面上固體和流體的質量;ρ,p,uf,us,ufs分別為流體的密度、壓力、流體速度、結構速度和流固耦合面速度;fρ,fp,fu分別為流體的密度向量、壓力向量和速度矢量;fufs和fusf分別為流固耦合面上分別作用在流體和固體表面的荷載矢量;fs為作用在結構上的荷載矢量;L為荷載矩陣;σfs為界面上的流體荷載.
流固耦合方程求解方法主要分為2大類[7]:① 直接耦合算法(direct coupling method).該算法將結構、流場、耦合界面的物理量統(tǒng)一在一個方程組里進行直接求解,適用性寬泛,但由于計算需求龐大,發(fā)展比較緩慢.② 迭代耦合算法(iterative coupling method),又稱為分離式算法(partitioned method).其基本思想是流場、固體在各自的CFD和CSD程序中,通過流固耦合界面,完成雙向數(shù)據(jù)傳遞、交替更新,直至達到收斂.本文采用迭代耦合算法.
由于流體方程是非線性的,流固耦合方程也具有非線性特性.求解方程的過程實質上是一個反復迭代逼近真實解的過程.在迭代過程中,需要設立應力準則或位移準則來判斷是否收斂,即
(2)
(3)
依據(jù)隨機振動理論,結構在平穩(wěn)隨機脈動風荷載作用下的動力方程為
(4)
(5)
采用振型分解法對方程(4)進行解耦,可得廣義坐標下的運動方程為
(6)
求解結構的脈動響應可轉換為求解脈動位移根方差,具體步驟如下:
① 計算節(jié)點廣義荷載譜.依據(jù)隨機振動理論,可得脈動風壓荷載譜Spp和廣義荷載譜Sfjfk間的關系式為
(7)
空間中任意2點w,q間的脈動風壓荷載互譜為
(8)
式中,coh(r,n)為空間相干函數(shù);r為2點間距離;Spw(n)和Spq(n)分別為節(jié)點w和q的脈動風壓譜.
② 依據(jù)廣義荷載,計算廣義坐標功率譜.第k階和第j階模態(tài)響應的互功率譜為
(9)
式中,Hj(n)為第j階模態(tài)的頻響函數(shù).
③ 計算結構響應功率譜.節(jié)點i處的動力位移響應根方差σyi為
(10)
式中,φji為節(jié)點i的第j階振型位移值.
式(10)中包含了所有參振模態(tài)交叉項,計算量大,計算效率低.隨著結構節(jié)點數(shù)目的增加,與振型分解組合法(SRSS)計算方法相比,CQC的計算耗時呈指數(shù)型增長.因此,目前大多采用SRSS法,具體表達式為
(11)
作用在結構上的風荷載可以分為平均風荷載和脈動風荷載2個部分.在平均風荷載作用下,結構變形到達新平衡位置,結構剛度和頻率隨之改變.風荷載作用下,結構在新平衡位置隨機振動,剛度近似不變,根據(jù)式(11),可以計算得到結構脈動響應.
式(11)可轉化為
(12)
式(12)中等式右邊第1部分可記為背景分量σB,第2部分記為共振分量σR,即
(13)
(14)
定義共振分量占總風致響應的比值為
(15)
評價各階模態(tài)對結構響應貢獻的大小,可根據(jù)各振型響應的應變能在結構總應變能中所占比例來衡量[8].實際計算中,通??紤]有限階計算振型,通過下式計算各振型應變能比例,衡量計算振型對結構響應貢獻的大小:
(16)
風洞試驗不僅可以測試剛性建筑表面的風壓和體型系數(shù),對橋梁、高層等結構風致振動的測試也很有效.張拉結構動力相似準則在風洞試驗中難以滿足,本文中風洞試驗僅用于模型本身的風致振動特性研究,同時為FSI數(shù)值模擬分析和頻域法分析結果的正確性提供必要的佐證.
在風洞中,對邊長為1m的方形索膜進行氣彈模型試驗.膜材選擇702 Fluotop T2,厚度為0.52mm,密度為750g/m2,經(jīng)向、緯向抗拉強度分別為3000N/5cm,2800N/5cm,經(jīng)向彈性模量為558MPa,緯向彈性模量為521MPa.試驗流場選擇均勻紊流場,紊流強度為14%.風速選取為10,16,22m/s.膜表面預張力選取為0.75,1.50,2.25,3.00kN/m.分別對膜面測點的位移和加速度進行測量.圖1為試驗測點布置及試驗現(xiàn)場照片.
圖1 單片索膜氣彈風洞試驗
采用移動錘擊跑點法,根據(jù)測點加速度時程信號識別膜片振動特性.傾角α=45°的張拉膜在4個預應力水平下,不同測點識別出的前5階頻率結果見表1.
表1 不同預張力膜結構頻率(α=45°) Hz
當風速為22m/s、膜面傾角α=45°時,不同預張力水平下膜面測點的平均位移和根方差見圖2.由圖可知,隨著預張力水平的下降,張拉膜面各測點的位移均值呈非線性增加,且增加幅度逐漸變大.位移根方差則隨預張力水平降低而逐漸增大.
圖2 測點位移均值及根方差隨預張力的變化曲線
膜結構振動與周圍氣體形成的耦合效應會產(chǎn)生附加質量,影響結構動力特性參數(shù)的求解.考慮空氣,建立勢流體流固耦合的計算模型(見圖3).
圖3 單片索膜勢流體模型
當膜面傾角α=45°、預張力為2.24kN/m時,索膜結構數(shù)值模型頻率的計算結果見表2.由于結構振動模態(tài)較為密集,扣除重根,計算值與試驗結果較為吻合.
表2 索膜結構頻率試驗的識別及數(shù)值計算結果 Hz
當風速為16m/s、膜面傾角α=45°、膜面預張力為2.24kN/m時,測點D3的試驗平均值為11.0mm,位移根方差為0.86mm.提取結構模型測點D3的位移時程曲線,結果見圖4.數(shù)值計算測點D3的位移均值為11.2mm,位移根方差為0.81mm,這與試驗結果基本吻合.
圖4 測點D3數(shù)值模擬位移時程
根據(jù)風速16m/s時的測點位移、加速度時程曲線,通過FFT變化得位移功率譜(見圖5)和加速度功率譜(見圖6).位移譜可以體現(xiàn)出結構振動的低頻部分,而加速度譜則體現(xiàn)出結構振動的高頻部分.由位移時程曲線及自功率譜曲線可以看出,結構的振動頻率包含2個部分,一部分為風荷載的頻率成分,另一部分為結構自身固有的頻率成分.結構以比其自振頻率低得多的頻率做隨機振動,振動形式主要為脈動風引起的強迫振動.而加速度反應譜表明結構的振動是一個寬帶過程,脈動風場會激起結構的某些振型.
圖5 測點D3的位移時程及功率譜(v=16m/s)
采用單向流固耦合方法,計算結構的風致動力響應.先將流場中膜結構表面設置為固壁,進行鈍體繞流的瞬態(tài)分析,獲取結構表面的風壓,存儲到文件中.然后,在結構表面對應的節(jié)點上,完成膜結構在風壓作用下的隱式動力分析.
在風速為16m/s、膜面預張力為2.24kN/m、膜面傾角α=45°的工況下,利用FSI計算單向耦合和雙向耦合并對比其區(qū)別.由圖7可知,2種方法得到的結構響應差異很小,這主要是因為結構體型較為簡單,膜結構預張力較大,風作用下膜結構表面變形較小,結構表面的風壓分布主要受到結構宏觀尺寸的影響,局部形狀的變化對大尺度渦的影響較?。?/p>
圖6 測點D3的加速度時程及功率譜(v=16m/s)
圖7 單向耦合與雙向耦合比較
采用頻域法計算單片索膜脈動位移.由公式y(tǒng)=μσy計算脈動位移峰值,其中μ=2.2為保證系數(shù),結果見圖8.脈動位移峰值最大為1.94mm,對應位移根方差為0.88mm,計算結果與試驗結果(見圖4)接近.第1階模態(tài)所占振型能量比較大,達到98%.共振分量占脈動位移百分數(shù)約為16%.結合圖5可知,脈動響應以背景分量為主.
圖8 頻域法計算結果
圖9為一單層平面索網(wǎng),高24.64m,寬26m,作為玻璃幕墻的支撐體系,位于宜興東氿大廈入口.玻璃采用8mm+8mm的雙層夾膠玻璃,分格列數(shù)為17,行數(shù)為16.第1列和最后1列的分格尺寸為1750mm×1540mm,中間部分的分格尺寸為1500mm×1540mm.
圖9 東氿大廈及平面索網(wǎng)測點布置
建立了包含玻璃面板、索網(wǎng)、爪件、密封膠在內(nèi)的玻璃-索網(wǎng)結構整體計算模型.在基本風壓w0=450Pa作用下,索網(wǎng)玻璃幕墻的最大靜力位移為0.1621m.對結構進行模態(tài)分析,頻率計算結果見表3.
表3 索網(wǎng)幕墻前8階頻率和周期
考慮玻璃幕墻及S形裙房,建立流場,流體模型網(wǎng)格劃分結果見圖10,單元數(shù)為8.4×105,節(jié)點數(shù)為4.1×105.
圖10 流場網(wǎng)格劃分示意圖
結合CSD和CFD模型,進行00風向角(即來流方向與玻璃面板垂直)流固耦合計算.采用AR模型生成脈動風速,模型階數(shù)P=4,目標譜采用湍流尺度沿高度不變的Davenport風速譜,按照基本風壓換算成平均風速為27m/s的脈動風速時程,入口風速剖面沿高度不變.截取60s風速時程施加在風速入口(見圖11),玻璃面板設置為流固耦合面.
圖11 入口風速時程及功率譜
在玻璃面板上布置9個測點(見圖9),提取測點位移和壓力時程,進行統(tǒng)計分析,相關結果見表4.平均位移最大值為0.156m,位于玻璃面板中部;在基本風壓為w0=450N/m2的荷載作用下,結構最大靜力位移為0.162m.這二者基本相等.
表4 測點位移及風壓流固耦合結果
圖12為測點S5的位移、加速度和節(jié)點壓力時程,圖13為頻域分析得到的位移幅值譜和加速度功率譜.由圖12可知,在脈動風荷載作用下,索網(wǎng)幕墻主要做受迫振動,以背景響應為主,共振響應比重較大.由圖13可知,在脈動風作用下,結構的前幾階模態(tài)被激發(fā)出來,振動為窄帶過程.
圖12 測點S5的位移、加速度和節(jié)點壓力時程曲線
圖13 測點S5的位移、加速度功率譜
忽略玻璃玻璃面板,對單層索網(wǎng)進行頻域法分析,索網(wǎng)位移根方差最大值為54mm,位于索網(wǎng)中部.脈動位移峰值見圖14,最大值為118.8mm.第1階模態(tài)所占的振型能量比較大,達到87.3%;索網(wǎng)中部共振分量占脈動位移的百分數(shù)約為30%,共振響應比重較大,這與圖13的結論相同.
圖14 頻域法計算結果
如圖15所示,選擇世博軸1號陽光谷和2號陽光谷及其之間的膜結構,建立1∶40氣彈模型試驗,平面尺寸為5.812 m×3.044 m.
圖15 膜結果平面圖和試驗模型
風洞試驗時,來流分為均勻來流和模擬C類地貌的紊流.通過柵欄、尖劈和粗糙元模擬了C類地貌的邊界層風速剖面和湍流度.
試驗位移和加速度測點布置如圖 16所示.不同的風向、預張力水平、風速和來流的試驗工況列于表5.
圖16 試驗測點布置圖
表5 試驗工況
根據(jù)風洞測點加速度時程數(shù)據(jù),采用隨機子空間法識別低預應力膜結構頻率,結果見表6.根據(jù)風洞模型,建立包括膜面及其支撐系統(tǒng)的數(shù)值模型,將計算頻率與識別頻率對比可知,數(shù)值模型基本反映了結構的自振特性.
在均勻流、v=10m/s、α=0°、低膜面應力的工況下,膜面測點1、2、7、8的位移和加速度時程見圖17和圖18.
表6 數(shù)值模擬及風洞識別頻率 Hz
圖17 位移時程(均勻流,v=10m/s,α=0o,低膜面應力)
由于膜結構的空間造型比較復雜,均勻來流受建筑物干擾形成特征湍流,引起膜面振動.對測點位移時程和加速度時程進行功率譜分析,結果見圖19和圖20.由圖19可知,結構振動主要以低頻部分為主,做受迫振動.由圖20可知,結構振動是一個寬帶過程,某些振型會被脈動風激發(fā)出來.
采用混合網(wǎng)格建立CFD模型,單元數(shù)為660690,節(jié)點數(shù)為240761.圖21為結構及流域網(wǎng)格劃分示意圖.
圖18 加速度時程(均勻流,v=10m/s,α=0o,低膜面應力)
表7為測點位移平均值與位移根方差結果.由表可知,數(shù)值模擬的測點位移平均值與試驗值比較接近,測點1~測點6的位移根方差相差較大,測點7和測點8的根方差與試驗值較為接近.分析試驗現(xiàn)象,特征湍流可能引起部分膜片發(fā)生顫振的現(xiàn)象[10],測點振動幅度增大,FSI模擬無法再現(xiàn)試驗中類似顫振的現(xiàn)象,有待在網(wǎng)格精細度、模型參數(shù)、算法等方面進一步改進.
表7 0o風向角的試驗與數(shù)值模擬測點數(shù)據(jù)對比 mm
圖19 位移功率譜(均勻流,v=10m/s,α=0o,低膜面應力)
選取結構前200階的頻率和振型,采用頻域法計算結構脈動位移峰值分布,結果見圖22.由圖可知,最大值為0.58mm,位于試驗測點2,3之間,對應的位移根方差為0.26mm.表7中,該位置附近的脈動位移根方差為0.29~0.84mm.與數(shù)值模擬結果相比,頻域法能夠考慮結構高階模態(tài)對脈動位移的貢獻,由頻域法得到的位移根方差分布與試驗結果更加接近.
1) 在紊流風作用下,單片索膜做隨機受迫振動,振動以背景響應為主,振動為寬帶過程,脈動風會激起結構的某些振型.對于單片索膜結構,氣彈模型試驗、FSI數(shù)值模擬、頻域法結果較為接近.膜面局部振動對風場分布影響較小.
圖20 加速度功率譜(均勻流,v=10m/s,α=0o,低膜面應力)
圖21 網(wǎng)格劃分示意圖
2) 單層索網(wǎng)玻璃幕墻在風荷載作用下振動為窄帶過程,共振分量所占比重較大.不考慮玻璃面板,對單層索網(wǎng)進行頻域法分析,由于忽略了玻璃的剛度,與索網(wǎng)幕墻數(shù)值模型相比,單層索網(wǎng)測點的平均位移、脈動位移都有不同程度的增大.工程中將玻璃幕墻簡化為索網(wǎng),通過頻域法計算結構風振響應偏于保守.
圖22 脈動位移峰值分布
3) 均勻來流在膜結構屋蓋表面形成特征湍流,使得某些膜片振動幅度較大,發(fā)生了類似顫振的現(xiàn)象.FSI數(shù)值模擬無法再現(xiàn)試驗中類似的現(xiàn)象.基于隨機振動理論的頻域法能夠考慮高階模態(tài)對振動的貢獻,脈動位移結果更接近試驗值.
4) 對于張拉結構的風致振動響應研究,仍需在現(xiàn)場實測、風洞試驗、數(shù)值模擬和理論方法等各方面開展長期、深入的研究.
)
[1] Vickery B J, Majowiecki M. Wind induced response of a cable supported stadium roof [J].JournalofWindEngineeringandIndustrialAerodynamics, 1992,42(1/2/3): 1447-1458.
[2] Nakamura O, Tamura Y, Miyashita K, et al. A case study of wind pressure and wind-induced vibration of a large span open-type roof [J].JournalofWindEngineeringandIndustrialAerodynamics,1994,52(1/2/3):237-248.
[3] Yasushi Uematsu, Theodore Stathopoulos, Eri Iizumi. Wind loads on free-standing canopy roofs. Part 1: local wind pressures [J].JournalofWindEngineeringandIndustrialAerodynamics, 2008,96(6/7): 1015-1028.
[4] Kinya Ando, Atush Ishi, Toshio Suzuki, et al. Design and construction of a double membrane air-supported structure [J].EngineeringStructures, 1999,21(8): 786-794.
[5] Michalski A, Kermel P D, Haug E, et al. Validation of the computational fluid—structure interaction simulation at real-scale tests of a flexible 29 m umbrella in natural wind flow [J].JournalofWindEngineeringandIndustrialAerodynamics, 2011,99(4): 400-413.
[6] 孫曉穎, 武岳, 沈世釗. 薄膜結構流固耦合效應的簡化數(shù)值模擬方法[J]. 土木工程學報, 2010, 43(10): 30-35.
Sun Xiaoying, Wu Yue, Shen Shizhao. A combined numerical approach for wind-structure interaction of membrane structures[J].ChinaCivilEngineeringJournal, 2010,43(10): 30-35.(in Chinese)
[7] Bathe K J. Adina theory and modeling guide [EB/OL]. (2008-02-02)[2013-03-01]. http://www.adina.com.
[8] Masanao Nakayama, Yasuhito Sasaki, Keiji Masuda, et al. An efficient method for selection of vibration modes contributory to wind response on dome-like roofs [J].JournalofWindEngineeringandIndustrialAerodynamics, 1998,73(1): 31-43.
[9] 中國建筑科學研究院. JGJ 102—2003玻璃幕墻工程技術規(guī)范[S]. 北京:中國建筑工業(yè)出版社, 2003.
[10] Vitale A M, Letchford C M. Experimental study of wind effects on flat porous fabric roofs [C]//ProceedingsofWindEngineeringintothe21stCentury. Rotterdam, the Netherlands, 1999:1545-1555.