申躍奎 李 洋 劉晨迪 王佳佳

(西安建筑科技大學土木工程學院，西安 710055)

1 概 述

充氣膜結構是一種典型的曲面張力空間結構，廣泛應用于大型體育場館、工作車間、倉庫等建筑上，具有極大的發(fā)展空間和前景[1]。然而，充氣膜結構作為一種柔性結構，對風十分敏感[2]。同時，結構的流固耦合效應和脈動風效應對結構的風壓分布和風致響應具有較大影響[3]，這在大多數(shù)的實際設計中并未同時考慮，此外，國內(nèi)對于截橢球形大跨空間結構的抗風設計相關研究也相對較少。

在諸多充氣膜結構的風振特性的研究方法中，計算流體動力學技術(CFD)因適應性強、周期短、費用低等優(yōu)勢而得到廣泛應用[4]。本文基于CFD技術和計算結構動力學技術(CSD)研究截橢球形充氣膜結構的風荷載響應特性。首先，利用線性濾波器法通過MATLAB得到考慮脈動風的風速時程曲線，并編制UDF將其導入到Flunet。其次，利用Workbench 17.0平臺數(shù)值模擬風致響應穩(wěn)定后不同風向角、高度、內(nèi)壓下剛性結構和柔性結構的風壓系數(shù)及位移響應分布，并對得到的相關結果進行分析。

2 數(shù)值模擬原理

2.1 湍流的數(shù)值模擬

一般情況下，工程上大多將常態(tài)風假設為不可壓縮黏性低速流體，并需要遵循質量守恒定律和能量守恒定律。同時，為了考慮并簡化流體的湍流運動，人們給出了幾種不同的湍流數(shù)值模擬方法。其中，求解時均化Reynolds方程的Reynolds平均法(RANS)計算量較小，因而在工程上得到了廣泛的應用[5]。本文采用了RNGk-ε模型，因為RNGk-ε模型在修正后能更好地處理彎曲流線流動[6]。該湍流模型的控制方程如下：

C1ε=1.42，C2ε=1.68，η0=4.377

β=0.012，αk=αε=1.39

式中：ρ為流體密度；k為湍動能；μi為湍動黏度；μeff為湍動能有效黏性系數(shù)；ε為湍動能耗散率；Gk為由平均速度梯度而導致的湍動能產(chǎn)生項。

2.2 流固耦合數(shù)值模擬方法

對于充氣膜結構這類柔性大變形結構，需要考慮結構和流體之間的相互耦合作用。本文采用的是分區(qū)弱耦合算法，即在每一時間步內(nèi)分別獨立求解流體和結構的控制方程，然后將流體域中求得的風荷載通過耦合交界面?zhèn)鬟f給結構域，得到結構的變形，然后再將結構的變形通過耦合交界面?zhèn)鬟f給流體域，往復進行計算，直至結果趨于穩(wěn)定。Glück等[7]較早地采用結構和流體的商業(yè)計算軟件實現(xiàn)了薄膜結構的流固耦合數(shù)值模擬。Nayer等[8]在文獻[9]氣彈模型試驗的基礎上，基于大渦模擬的分區(qū)耦合算法進行了在不同雷諾數(shù)下充氣截球形膜結構流固耦合效應的數(shù)值模擬研究，并根據(jù)渦脫落頻率和斯特勞哈爾數(shù)對渦脫落進行了分類。文獻[10-12]基于分區(qū)弱耦合算法分別對二維張拉膜結構屋面、山脈形膜結構及典型細長索膜結構進行了流固耦合數(shù)值模擬，得到了結構的位移響應時程。

3 模型參數(shù)設置

3.1 流體域設置

建立2 500 m×1 000 m×400 m的長方體為流場域，結構縱剖面與風場剖面的面積阻塞率最大值僅為1%，能夠滿足阻塞率小于3%的要求[13]，并且結構距邊界入口5倍長度。根據(jù)我國GB 50009—2012《建筑結構荷載規(guī)范》和GB 55001—2021《工程結構通用規(guī)范》，該工程處于B類地貌，基本風速為20 m/s，相當于基本風壓0.4 kN/m2，梯度風高度350 m，標準高度10 m，地面粗糙度0.16。流體域入口采用指數(shù)率平均風剖面，相關算式如下：

(2a)

(2b)

ω=0.090.75k0.5/L

(2c)

將MATLAB中按照達文波風速譜模擬生成的脈動風速時程通過UDF輸入到Fluent環(huán)境中，參考的線性濾波器法原理詳見文獻[14]。風場高度10 m處的風速時程如圖1所示。同時，將該風速譜進行FFT變換，與Davenport進行對比，如圖2所示?？梢园l(fā)現(xiàn)，兩者功率譜較為吻合，說明了該方法具有可行性。

圖1 10 m高度處的風速時程Fig.1 Time history of wind speed at 10 m height

流體域兩側及頂部采用滑移對稱邊界條件，底部采用無滑移壁面邊界條件，出口采用自由出流邊界條件[15]。模型采用四面體網(wǎng)格劃分，同時要注意對膜面及附近流場局部進行加密。通過試算和網(wǎng)格無關性檢驗，確定該處最小網(wǎng)格尺寸0.5 m，增長率1.2，網(wǎng)格總數(shù)達約385萬個，底部網(wǎng)格劃分情況詳見圖3。

圖3 流體域底部網(wǎng)格劃分Fig.3 Mesh generation at the bottom of fluid domain

求解器設置方面，本文對流體控制方程采用SIMPLEC算法，選用二階迎風格式對控制方程進行離散，同時模擬過程中迭代計算的收斂標準為無量綱連續(xù)性殘差達到 10-5并且結構監(jiān)測點風壓系數(shù)趨于穩(wěn)定，前后迭代誤差在5%內(nèi)。

3.2 結構域設置

實際工程中充氣膜結構膜材采用1.0 mm厚的PVDF涂層膜材和0.8 mm厚PVC保溫內(nèi)膜，膜材表面附加有外包PE的φ10鋼索。為了簡化模型，根據(jù)文獻[2]給出的方法對模型進行等效處理。本文假定索與膜面無相對滑移，將膜材和索等效為均勻材料，等效后模型在均布荷載作用下所產(chǎn)生的位移相同。最后得到模型的等效密度為2 500 kg/m3，膜厚1 mm，等效彈性模量為1×109Pa，泊松比為0.34。結構底面長軸150 m，短軸100 m，底邊采用固定約束，網(wǎng)格采用Shell 181劃分，結構具體尺寸如圖4所示。

a—俯視圖；b—正視圖。圖4 結構尺寸示意 mFig.4 Sizes of the air-supported structure

3.3 流固耦合面

將結構與流場的交界面定義為流固耦合交界面，該面作為流體域和結構域數(shù)據(jù)傳遞的交互面。

3.4 變量設置

充氣膜結構的剛度主要由內(nèi)壓決定[16]，同時，相關風洞試驗[17-18]也證實了結構的高度和風向角會對結構表面風壓分布產(chǎn)生較大影響。因此，本文主要的研究變量為風向角、高度和內(nèi)壓。其中，風向角為變量時，高度取30 m，內(nèi)壓取400 Pa，風向角分別取0°、15°、30°、45°、60°和90°；高度為變量時，風向角取0°，內(nèi)壓取400 Pa，高度分別取26,30,34 m；內(nèi)壓為變量時，風向角取45°，高度取30 m，內(nèi)壓分別取400，500，600，700 Pa。

4 準確性驗證

4.1 實例模型

以某鞍形張拉膜為例，基于上述方法數(shù)值模擬得到了不同風向角下的風壓系數(shù)分布，并與同濟大學TJ2風洞的試驗結果[19]進行對比，來說明這一方法的可行性。

該鞍形張拉膜具體參數(shù)和模型見表1和圖5，流體域尺寸大小為1 600 m×800 m×80 m，模型中心距離流場入口處 200 m，阻塞率僅為 0.5%。流場域如圖6所示。網(wǎng)格選用四面體單元，膜面及附近流場局部需要加密。流場邊界條件、湍流模型和求解器設置與前文一致。

表1 算例模型參數(shù)Table 1 Model parameters of the example

a—算例模型；b—模型平面圖[19]；c—模型立面。圖5 算例模型尺寸Fig.5 Example model sizes

圖6 流場域Fig.6 Flow field

4.2 結果對比

B類地貌、15 m/s基本風速下，風向角為0°、45°和90°時數(shù)值模擬和風洞試驗方法下結構的風壓系數(shù)分布情況詳見圖7～9。

a—數(shù)值模擬；b—風洞試驗[19]。圖7 0°風向角風壓系數(shù)Fig.7 Wind pressure coefficients under wind direction angle 0°

a—數(shù)值模擬；b—風洞試驗[19]。圖8 45°風向角風壓系數(shù)Fig.8 Wind pressure coefficients under wind direction angle 45°

a—數(shù)值模擬；b—風洞試驗[19]。圖9 90°風向角風壓系數(shù)Fig.9 Wind pressure coefficients under wind direction angle 45°

可以看出，在3個不同風向角下數(shù)值模擬和風洞試驗得到的風壓系數(shù)分布大致相同，僅在數(shù)值上和極少部分區(qū)域分布上有差異，這主要可能是風洞試驗相似性原理所產(chǎn)生的誤差。

5 結果分析

5.1 結果表示

本文主要關注的是風向角、高度和內(nèi)壓變化對結構表面風壓分布的影響，由于篇幅限制，所以僅給出了結構平均化的風壓系數(shù)分布及匯總情況的部分結果。

5.2 不同風向角下結構的風壓分布

表2給出了不同風向角下充氣膜結構柔性模型的風壓系數(shù)極值匯總。圖10和圖11給出了0°風向角30 m高剛性模型以及0°、45°、60°和90°風向角30 m高400 Pa內(nèi)壓柔性模型的風壓系數(shù)分布。圖12給出了各風向角下柔性模型風壓系數(shù)沿長軸與短軸的分布情況。

表2 不同風向角下剛性與柔性模型風壓系數(shù)極值Table 2 Extreme values of wind pressure coefficients of rigid model and flexible model under different wind directions

a—剛性；b—柔性。圖10 0°風向角下剛性模型和柔性模型的風壓系數(shù)分布Fig.10 Wind pressure coefficient distribution of rigid model and flexible model under wind direction angle 0°

從圖10可以發(fā)現(xiàn)：相同風向角下剛性和柔性模型的表面風壓分布大致相同，僅在數(shù)值上有差異，剛性模型中正壓和負壓系數(shù)極值分別為0.853和-0.738，柔性模型中正壓和負壓系數(shù)極值分別為0.784和-0.763；同時結合表2發(fā)現(xiàn)結構總體上受負壓極值控制，并且考慮流固耦合效應的柔性模型負壓系數(shù)極值更大，風壓分布更均勻[2]。

a—45°；b—60°；c—90°。圖11 45°、60°和90°風向角下柔性模型的風壓系數(shù)分布Fig.11 Wind pressure coefficient distribution of flexible model under wind direction angle 45°,60°and 90°

a—短軸；b—長軸。圖12 各風向角下柔性模型的風壓系數(shù)在短軸與長軸上的分布Fig.12 Wind pressure coefficient distribution of flexible model on short axis and long axis under different wind directions

從圖10b和圖11可以發(fā)現(xiàn)，隨著風向角的變化，結構風壓的大小和分布都發(fā)生了較大變化。不同風向角下等值線呈橢圓狀分布，風壓帶大致沿風向垂直方向分布，并且風向角0°和90°時沿結構的中軸線對稱分布；大部分區(qū)域以負壓為主，最大負壓出現(xiàn)在結構頂部中心區(qū)，此處容易產(chǎn)生強烈的流動分離；最大正壓出現(xiàn)在迎風面底部，同時由于流動再附[13]，在背風面也會存在正壓區(qū)域。隨著風向角的增大，風壓分布帶逐漸變得狹長。風向角90°時結構正壓和負壓系數(shù)極值達到最大，分別為1.35和-2.01，說明該工況下結構最容易因受到較大上拔力而發(fā)生破壞。

從圖12可以發(fā)現(xiàn)，各風向角下風壓系數(shù)值大致沿短軸與長軸呈對稱拋物線形分布。結構頂部中心處負壓最大，風壓梯度變化最大；風向角0°和45°時結構沿短軸與長軸兩端30 m和25 m范圍內(nèi)存在正壓；隨著風向角增大，風壓系數(shù)沿短軸與長軸的變化幅度增大，90°時變化最大。同時，風向角0°和45°時風壓系數(shù)沿長軸的變化幅度更大，60°和90°時風壓系數(shù)沿短軸變化幅度更大。

5.3 不同高度下結構的風壓分布

圖13給出了26 m和34 m柔性模型的風壓系數(shù)分布，圖14給出了26 m和34 m給出了各風向角下柔性模型風壓系數(shù)沿短軸與長軸的分布情況。

a—26 m；b—34 m。圖13 26 m和34 m柔性模型的風壓系數(shù)分布Fig.13 Wind pressure coefficient distribution for 26 m and 34 m flexible models

從圖10b和圖13可以看出，隨著結構高度的增加，結構頂部負壓系數(shù)極值分布帶逐漸變得狹長，負壓系數(shù)極值逐漸增大，26，30，34 m時分別為-0.691、-0.763和-0.865。從圖14可以看出，風壓系數(shù)值大致沿短軸與長軸呈對稱拋物線形分布。其中結構頂部中心處負壓最大，風壓梯度變化也最大；結構沿長軸兩端25 m范圍內(nèi)存在正壓，沿短軸方向都為負壓；隨著高度增加，風壓系數(shù)沿短軸與長軸的變化幅度逐漸增大，34 m時變化最大，并且結構沿長軸的風壓系數(shù)變化幅度更明顯。

5.4 不同內(nèi)壓下結構的風壓分布

圖15給出了內(nèi)壓為500，600，700 Pa時柔性模型的風壓系數(shù)分布，圖16給出了各內(nèi)壓下柔性模型在短軸與長軸的分布情況。

a—500 Pa；b—600 Pa；c—700 Pa。圖15 內(nèi)壓500，600，700 Pa的柔性模型的風壓系數(shù)分布Fig.15 Wind pressure coefficient distribution of flexible models with internal pressure of 500,600，700 Pa

a—短軸；b—長軸。圖16 各內(nèi)壓下柔性模型的風壓系數(shù)在短軸與長軸上的分布Fig.16 Wind pressure coefficient distribution of flexible model on short axis and long axis with different internal pressures

從圖11a和圖15可以看出，隨著內(nèi)壓的增加，負壓系數(shù)極值帶逐漸變得狹長，正壓區(qū)位于結構的迎風面和背風面。而從不同內(nèi)壓下結構正負壓系數(shù)極值變化趨勢發(fā)現(xiàn)，隨著內(nèi)壓的增加，結構正壓和負壓系數(shù)極值先增加后分別保持在1.41和-1.45左右不變。

從圖16可以看出，風壓系數(shù)值大致沿短軸與長軸呈對稱拋物線形分布。其中結構頂部中心負壓最大，風壓梯度變化也較大；結構沿長軸兩端25 m范圍內(nèi)存在正壓，沿短軸方向都為負壓；隨著內(nèi)壓增加(500～700 Pa)，風壓系數(shù)沿短軸與長軸的變化幅度較小，說明內(nèi)壓的增加不一定能有效地減少風壓系數(shù)極值，該工況下內(nèi)壓宜采用400 Pa。

6 結束語

本文研究了考慮流固耦合與脈動風效應下不同風向角和結構參數(shù)對截橢球形充氣膜結構的風壓系數(shù)的影響，得到了以下結論。

1)不同風向角下，結構風壓系數(shù)分布變化較大，分布帶基本都沿垂直風向呈橢圓狀分布；結構正壓區(qū)分布范圍較小，主要集中在迎風面和背風面；結構負壓區(qū)分布較大，最大負壓區(qū)在結構頂部中心位置，在風向角為90°時有最大負壓系數(shù)極值；各風向角下風壓系數(shù)值大致沿短軸與長軸呈對稱拋物線形分布，并且隨著風向角的增大變化幅度增大，0°和45°時風壓系數(shù)沿長軸的變化幅度更大，60°和90°時風壓系數(shù)沿短軸變化幅度更大。

2)不同高度下，結構表面風壓分布大致相同，頂部負壓系數(shù)極值分布帶逐漸變得狹長，并隨著高度的增加負壓極值增大；隨著高度增加，風壓系數(shù)沿短軸與長軸的變化幅度逐漸增大，并且結構沿長軸的風壓系數(shù)變化幅度更明顯。

3) 隨著內(nèi)壓的增加，結構負壓系數(shù)極值先增加后保持在-1.45左右不變，說明內(nèi)壓的增加不一定能有效地減少風壓系數(shù)極值，在其他條件相同的情況下內(nèi)壓宜采用400 Pa。