三維鈍體結構風壓分布與繞流特性的大渦模擬研究

汪大洋，周 云，譚 平

(1.廣州大學土木工程學院，廣東 廣州 510006;2.廣州大學工程抗震研究中心，廣東 廣州 510405)

0 引言

鈍體繞流與風壓分布特性研究是結構風工程中非常重要的內容，是認識結構周圍風環境、掌握結構表面風壓分布的主要途徑之一，對結構抗風設計具有重要的參考、指導價值。近年來，國內外研究人員在模擬預測全尺寸或縮尺比例結構模型的繞流狀態、風壓分布、相互干擾效應等方面進行了較多的研究工作。Selvan［1］等采用LES(Large Eddy Simulation)模型對一足尺風洞實驗進行了研究，結果表明對平均值的預測與實測結果吻合較好，但對峰值壓力，只有根據實測數據生成脈動來流條件的模擬結果與實測吻合較好。Tutar［2］等人采用了大渦模擬法對兩棟平行建筑物進行了數值研究，給出了平行布置的建筑物的風荷載和流場的分布特點。Muralami［3］、Yu Da-hai［4］應用 LES 模型對鈍體繞流進行的研究表明，LES模型雖能準確反映繞流的復雜特征，但因計算量大而難以在工程中應用。Cheng［5］等采用LES和RANS(Reynolds-Averaged Navier-Stokes)對比研究了群體結構繞流特性，結果表明二者均能很好預測群體結構平均氣流的主要特性，但LES預測結果與試驗的吻合性優于 RANS。Shah和 Ferziger［6-7］通過研究得出LES模型能為風工程研究人員提供風壓分布特性、繞流狀態等重要信息。Hee［8］等采用LES研究了一位于湍流邊界層的矩形建筑，結果表明合理的來流及計算域邊界可大大提高LES與試驗結果的吻合性。李正良［9］等針對復雜高層進行了風洞和數值研究，表明RSM模型能較好地模擬結構周圍的流場，試驗與數值結果吻合性較好。王小華［10］等對方形鈍體受限繞流進行了模擬，結果表明均勻來流下湍流場沿槽道軸面對稱，LES能捕捉到豐富的渦系及渦動時變過程。吳文權［11］、陳斌［12］等采用離散渦法對多個鈍體繞流的研究表明，即使來流是穩定的，繞多個鈍體流動也會出現不穩定的流動。

雖然已有學者對鈍體結構的繞流問題進行了相關研究，但由于鈍體繞流的復雜性，流場與結構之間的很多問題還沒有清晰解決，且已有研究大多針對均勻來流進行的，較少考慮來流脈動特性，而實際工程顯然處于湍流脈動風場中，存在復雜的湍流脈動現象。因此，有必要對繞流問題進一步開展研究。本文針對三維矩形鈍體結構，以平均風速和脈動風速作為入口邊界條件，采用大渦模型研究該鈍體的分壓分布與繞流特性，再現鈍體繞流的沖撞、分離、重附著及尾流漩渦不對稱脫落現象，揭示鈍體繞流分離、重附著、漩渦生成、發展等的形成機理，并闡述鈍體繞流對其來流方向、上部及尾流速度場及表面風壓場分布的影響規律。

1 控制方程與離散方法

1.1 LES控制方程

大渦模擬的基本方法是對大尺度渦旋進行直接模擬，對小尺度渦旋作模式假定，因此大渦模擬必須引入濾波函數來修改Navier-Stokes方程(N-S方程)，使得高波數的波被截斷，但能量傳遞過程仍保留。經濾波函數濾波后，LES控制方程為:

Smagorinsky認為，對于局部各向同性湍流，如果假定湍流的產生項和耗散項達到平衡，則亞網格尺度應力可定義為:

式(3)即為代數壓網格尺度模型，又稱Smagorinsky模型［13］，式中應變率張量為:

1.2 數值離散與求解

采用有限體積法離散過濾后的不穩定N-S方程，采用時間積分算法進行時間積分，進而從速度項中分離壓力項。采用二階精確中心差分逼近法來離散空間LES控制方程，以 QUICK(Quadratic Upstream Interpo-lation Convective Kinematics)方式采用二階上游有限差分法離散空間對流項［14，15］。若 n△t時刻變量已知，則(n+1)△t時刻變量的計算如下:

1)將二階后差歐拉法與顯示二階Adams Bashforth算法相結合，考慮已知壓力場，采用二階時間法計算中心速度項U*:

2)采用Helmholtz分解定理將中心速度場投射到離散自由向量場的子空間上:

令:φ =pn+1-pn。

3)利用超松弛迭代法通過求解泊松方程計算壓力場，并得到增量壓力:

4)計算(n+1)△t時刻的變量:

2 計算模型與邊界條件

2.1 計算模型

基于ADINA軟件建立計算模型，采用三維流體單元模擬流場，空氣流體質量密度取1.29kg/m3，粘度取1.79×10-5kg/(m·s)，采用FCBI-C單元對流體域進行離散［16］。鈍體縮尺比例為1∶200，風速相似系數取1/5，參考高度風速vref=5.66m/s。計算域的選取應足夠大，以使邊界條件對繞流的影響降至最低［17-18］，本文取l(x)×l(y)×H=0.6m×0.4m×1.0m，計算域為24H×18H×10H，其中鈍體前方 L1(x)=7H、后方 L2(x)=16H。分析工況及單元劃分情況見表1。圖1為工況1流體域與鈍體尺寸布置圖。先對計算域進行穩態分析，然后進行瞬態分析，瞬態分析的時間步長為0.001s，總時間為3.6s。

圖1 流體域與結構模型示意圖Fig.1 Sketch of fluid field and bluff body

2.2 邊界條件

為提高數值模擬的精確性，在數值風洞入口邊界條件上同時施加平均風速和脈動風速。平均風速剖面采用指數律:

式中:zref、vref為參考高度和相應的平均風速;z、v(z)為任意高度和相應的平均風速;α為地面粗糙度指數，取α=0.22。

脈動風速采用 AR法(Auto-Regressive)生成［19-21］，模型階數為5階，時間步長為0.001s，脈動風速譜為Davenport譜。脈動風速的計算公式如下(模擬功率譜與目標功率譜的對比如圖2所示):

式中:X={x1，x2，…，xm}，Y={y1，y2，…，ym}，Z={z1，z2，…，zm}，(xi，yi，zi)為空間第 i點坐標，(i=1.2，…，m);△t為時間步長;p為模型階數;φk為AR模型的自回歸系數矩陣，為m×m階矩陣，(k=1.2，…，p);N(t)為獨立隨機過程向量。

圖2 功率譜函數對比Fig.2 Contrast of power spectrums

入口湍流強度采用日本規范建議［22］:

表2給出了計算域與鈍體的參數設置。

表1 分析工況Table1 Analysis cases

表2 計算模型參數設置Table2 Calculation parameters

3 數值分析驗證

為驗證上述數值計算方法、邊界條件選取等的可行性及計算結果的可信度，首先進行一正方體鈍體的繞流特性分析，并將分析結果與已有研究對比［23-25］。該正方形鈍體所處流場尺寸為15L×5L×10L(L為鈍體特征尺寸)，模型X、Y、Z三個方向單元網格數量依次為48、43和38(如圖3所示)。

圖4顯示了流場橫風向風速剖面與已有研究的對比，圖5顯示了鈍體表面壓力系數分布與已有研究的對比。由圖4、圖5可知，雖然由于模型邊界條件、湍流度、阻塞效應等原因引起本文研究結果與已有試驗和數值研究結果局部存在一定的差異(如鈍體頂面、側面風壓分布)，但總體流場特征和風壓系數分布與已有研究仍具有較好的吻合性，說明通過本文方法進行流場數值模擬是可信的，進而為下文分析提供可靠保證。

同時，針對工況1采用不同湍流模型進行對比研究。圖6顯示了模型中剖面速度分布圖。可見，三種湍流模型風速剖面的吻合性很好，進一步證實了本文方法的可行性。

4 結果分析與討論

4.1 繞流特性

圖7顯示了工況1的切平面A的豎向風速剖面隨時間的分布。圖8顯示了工況1的豎向風速剖面圖。圖9顯示工況5的豎向風速剖面圖。由圖7、圖8和圖9可知:

①來流區風速剖面隨時間變化較小，在x/H＜-0.25范圍內，風速剖面因受鈍體影響較小而隨時間基本不變，但在-0.25＜x/H＜0內，風速剖面受鈍體阻礙作用隨時間有一定的變化，但變化幅度不大;屋頂區和尾流區風速剖面在初始時刻隨時間變化較大，但隨時間增長而逐漸減小，1.5s后基本無變化，說明流場已發展到穩定狀態;

②來流區z/H＜-0.25內，因鈍體對該區域流場影響較小，風速剖面沿高度基本按照指數律變化。在-0.25＜x/H＜0內，鈍體對豎向風速剖面的影響分為四個階段:在0＜z/H＜0.3內，受地面粗糙度及回流渦旋的影響，風速上升較快;在0.3＜z/H＜0.9內，因鈍體繞流緩沖的影響風速增長緩慢;在0.9＜z/H＜1.4內，因流場在z/H=1處分離，風速因阻礙作用瞬間減小而急劇上升;在z/H＞1.4時，因鈍體對流場的影響較小并逐漸減弱，風速剖面將重新按指數律變化;

③因屋頂區x/H=0處尚未形成渦旋狀態，使得該處風速剖面與其它位置截然不同。在0＜x/H＜0.6(工況1)、0＜x/H＜0.4(工況5)內，屋頂繞流渦旋已經生成，風速剖面沿高度方向先負向增長后正向上升，在屋頂形成負風壓區。屋頂鈍體對風速剖面的影響主要集中在1＜z/H＜1.4區域，并在該區域發生繞流分離、重附著、漩渦生成、分離等復雜現象;

④尾流區風速剖面變化比較復雜，受尾流渦旋的影響，在0.6＜x/H＜2.6(工況1)、0.4＜x/H＜2.4(工況5)范圍內呈“S”型變化。在x/H＞2.6(工況1)、x/H＞2.4(工況5)內，風速剖面隨鈍體對流場影響的逐漸減弱而向指數律轉變。鈍體對尾流豎向風場的影響集中在0＜z/H＜1.4內，并隨與鈍體之間距離的增大而逐漸減小。如在x/H=4.6處(工況1)，其影響范圍0＜z/H ＜1.05。

⑤鈍體對切平面A流場的影響要大于切平面C，如工況5的來流區在x/H=-0.125處，切平面C的風速剖面沿高度方向的變化趨勢比切平面A、B平緩;同樣在尾流區x/H=4.6處，切平面C的風速剖面由于受尾流影響較小而更加趨于指數律，切平面B次之，切平面A最差，原因在于鈍體對中心區域流場的阻礙作用大于兩側流場。

圖3 計算模型Fig.3 Computational domain and grid system

圖4 橫向風速剖面對比(x=2L，z=0.5L)Fig.4 Lateral velocity profile(x=2L，z=0.5L)

圖5 鈍體表面壓力系數分布對比Fig.5 Surface pressure coefficients

圖6 風速剖面分布對比圖Fig.6 Comparison of wind velocity profiles in different turbulence models

圖7 豎向風速剖面隨時間變化圖Fig.7 Vertical velocity profiles with time

圖8 工況1的豎向風速剖面圖Fig.8 Vertical velocity profiles of the case 1

圖9 工況5的豎向風速剖面圖Fig.9 Vertical velocity profiles of the case 5

圖10 切平面示意圖Fig.10 Sketch map of cross section

圖11顯示了工況1和工況2的鈍體繞流跡線示意圖。表3給出了工況1的鈍體繞流隨時間變化的特征值。表4給出了不同工況的鈍體繞流特征值(表中“∞”表示渦旋尚未形成)。為更好地揭示鈍體繞流沖撞、分離、重附著及尾流渦旋的生成、發展情況，定義:將流場經過鈍體拐角處發生分離的點稱為分離點，將流場形成渦旋后再次與鈍體或地面接觸的點稱為重附著點，將流場與鈍體接觸后而向上、下或左、右反向流動的點稱為變向點，將流場跨越鈍體后而在尾流區重新匯合的點稱為交匯點，將流場經過鈍體時直接產生的渦旋區稱為主渦循環區，將由于主渦循環區對流場再次擾動而形成的渦旋區稱為次渦循環區。由圖11及表3、表4可知:

圖11 鈍體繞流跡線示意圖Fig.11 Flow trace sketch

①來流區鈍體繞流特性比較穩定，隨時間基本無變化。如切平面A在不同時刻Z2/H均為0.8(工況1)，切平面D在不同時刻均為0.3(工況5)。

②不同切平面來流區的變向點特征值由中心向兩側減小，其原因在于流場流經鈍體時，鈍體對流場的阻礙作用由中心向兩側降低。如3.0s時刻切平面A、B和C的特征值Z2/H依次為0.8、0.75、0.7(工況1)。

③不同工況下尾流區均形成明顯的非對稱漩渦脫落現象，這與高層結構橫風向振動往往大于順風向振動的現象是一致的。其原因在于尾流渦旋的非對稱脫落使結構左右兩側流場形成非對稱流動，并產生不平衡的作用力，從而使結構產生橫風向振動，當這種不平衡作用力上升到一定程度時，便會產生橫風向振動，甚至共振現象的發生。

④圖11清晰再現了鈍體繞流的分離、沖撞、重附著及渦旋形成等現象，以圖11(b)為例加以說明。流場經過鈍體形成主渦循環區IV、V，并產生流線K和J。流線K向左上方a點運動，受鈍體的影響在變向點c處分流，分別向b、d流動。在a-d區段內由于流線K與鈍體發生沖撞、分離現象，而使得部分流場背離鈍體向e流動，并在流線K、J之間形成次渦循環區III;另一部分在流線J和渦旋V的影響下向g運動，但由于f點為分離點，且流線K在d點處仍具有一定的速度，所以使得流線K不會平行鈍體運動，而是產生弧線流d-g-h，并形成重附著點i和次渦循環區I。由于流線K與鈍體在i點產生沖撞重附著，此時如果流線K的速度較小，其流動方向將會平行鈍體流動;如果流線K仍具有一定速度，則會導致流線K的流動方向仍然不會平行鈍體表面，而是再次產生弧線流j-k-l，并形成次渦循環區II。

⑤交匯點繞流特征值Xn/H隨風向角的變化而變化，工況3和工況7的特征值最小(如兩種工況下切平面G的特征值分別為1.45和1.42)，而工況5的特征值最大(Xn/H=2.25)。說明對于矩形鈍體而言，來流方向與其成45°或135°角時，鈍體對尾流的影響最小，并隨來流方向逐漸垂直于鈍體，其對尾流的影響也逐漸增大，如工況 7、8、9 的特征值依次為 1.42、1.7、2.05(G平面)。鈍體對尾流影響區域均隨高度的增大而減小，如工況6的切平面G、F、E、D的特征值Xn/H依次為1.85、1.7、1.6、1.2。

4.2 空間風壓分布特性

為探討不同工況下鈍體繞流的空間風壓分布特性，引入平均壓力分布系數:

式中:pi為第i點相對壓力，qref為參考風壓。

圖12位鈍體壓力參考線示意圖。圖13為工況1在不同時刻各參考線壓力系數的分布情況，圖14為不同工況下各參考線的平均壓力系數分布情況。由圖13、圖14可知:

圖12 參考線選取Fig.12 Selection of reference lines

①在1.5s時刻之前，各參考線的壓力系數波動較大，流場未達到穩定狀態;1.5s之后流場逐漸進入穩定狀態，各參考線上壓力系數分布漸趨穩定，波動幅度明顯減小。

②迎風面壓力系數分布比較均勻(如工況1、工況5、工況9)，沿高度方向均呈現由大變小再逐漸增大并在頂部再次回落的趨勢，最大值出現在變向點附近。其余工況迎風面壓力系數亦有類似結論，但壓力系數相對較小。

③屋頂區壓力系數在各工況下均為負值(參考線B、E)，平均風壓為吸力。當風場垂直鈍體時，屋頂區風場在迎風面前緣分離，分離處風速動能增加，使得屋頂上部一定范圍內的風場趨于較低的分離流線，并沿分離流線方向形成柱狀渦旋，進而造成了屋頂區沿分離流線方向的壓力系數基本相等，如工況1、工況9參考線E上的壓力系數;而當風場與鈍體存在一定夾角時，迎風面前緣分離流線變為兩條，沿每條分離流線方向的流場分離速度逐漸增大，進而沿著分離流線方向形成錐形渦，這就使得平行分離流線方向的其它參考線上的壓力系數不相等，如工況7的參考線B、E的壓力系數。

④隨著風向角的增加，鈍體迎風面逐漸變為側風面(工況5)、背風面(工況9)，壓力系數分布也隨之變化，由壓力(正壓區)逐漸變為吸力(負壓區)，在來流區(0°～60°)由于流場分布比較規則均勻，壓力系數隨風向角的增加而逐漸降低，在側向繞流區和尾流區(60°～180°)由于渦旋生成、脫落的不規則性，壓力系數分布也呈現不規則變化趨勢(如參考線A、G)。此外，對于工況4，雖然參考線A、G所在面仍處于來流區，但由于其與流場方向的夾角過大，使得該面上形成微小渦旋而產生負壓區，但壓力系數相對其它參考線要小得多。

表3 工況1的繞流特征值Table3 Flow characteristic dimensions of the case 1

表4 交匯點繞流特征值(2.0s)Table5 Flow characteristic dimensions of the intersection points at 2.0s

圖13 工況1的平均壓力系數分布Fig.13 Mean pressure coefficients of case 1

圖14 平均壓力系數分布Fig.14 Mean pressure coefficients

圖15顯示了工況1各面參考線的脈動風壓分布情況。可見，脈動風壓系數直接反應繞流渦旋運動的強烈程度，渦旋運動越強烈，脈動壓力系數也越大，如鈍體負壓面的脈動壓力系數均高于正壓面，說明脈動風壓反映鈍體表面壓力脈動能量的大小，對正確認識鈍體表面脈動特性具有參考價值。

圖15 脈動風壓系數分布Fig.15 Fluctuating pressure coefficients

5 結論

針對三維矩形鈍體采用大渦模型進行了數值模擬，再現了鈍體繞流的沖撞、分離、重附著等現象，揭示了鈍體繞流主、次渦循環區的形成、發展機理，得到以下結論:

(1)矩形鈍體對來流區0.25、屋頂區0.4、尾流區2.0內流場的影響最大，對中心區域流場的影響要大于兩側流場，影響程度由中心向兩側逐漸減小，且其強弱程度依次為尾流區、屋頂區、來流區。

(2)來流區風場與鈍體建筑物之間的繞流特性比較穩定，同一工況下任一切平面的繞流特征值隨時間基本無變化，但不同切平面的繞流特征值不同，由中心向兩側逐漸減小。鈍體與風場垂直時，其對尾流的影響最大;呈45°時對尾流的影響最小。

(3)尾流渦旋的非對稱生成、脫落使鈍體左右兩側流場與鈍體之間形成不對稱的沖撞、分離與重附著現象，導致鈍體兩側形成不平衡壓力分布，是結構產生橫風向振動的重要原因。

(4)當風場與鈍體垂直時，迎風面豎向壓力系數均呈現由大變小再逐漸增大并在頂部再次回落的變化趨勢，最大值出現在變向點附近，而其余工況下雖壓力系數亦有類似結論，但壓力系數值相對較小。脈動風壓系數直接反應繞流渦旋運動的強烈程度，渦旋運動越強烈，脈動壓力系數也越大。

［1］SELVAN N R.Computation of pressures on texas tech university building using large eddy simulation［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，67＆68:647-657.

［2］TUTAR M，OGUZ G.Large eddy simulation of wind flow around parallel buildings with varying configuration［J］.J.Fluid Dynamics Research，2002，9:289-315.

［3］MURALAMI S，MOCHIDA A.On turbulent vortex shedding flow past2D square cylinder predicted by CFD［J］.Journal of Wind Engineering and Industrial Aerodynamics，1995，54(3):191-211.

［4］YU D H，ASHAN K.Numerical simulation of flow around rectangular prism［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，67(2):195-208.

［5］CHENG Y，LIEN F S，YEE E，SINCLAIR R.A comparison of large eddy simulation with a standard k-εReynolds-averaged Navier-Stokes model for the prediction of a fully developed turbulent flow over a matrix of cubes［J］.Journal of Wind Engineering and Industrial Aerodynamics，2003，91:1301-1328.

［6］SHAH K B，FERZIGER J H.Large-eddy simulations offlow past a cubical obstacle［R］.Thermosciences Division Report TF-70.Department of Mechanical Engineering，Stanford U-niversity，1996.

［7］SHAH K B，FERZIGER J H.A fluid mechanician＇s view of wind engineering:large-eddy simulation of flow past a cubical obstacle［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，67＆68:211-224.

［8］HEE C L，THOMAS T G，CASTRO I P.Flow around a cube in a turbulent boundary layer:LES and experiment［J］.Journal of Wind Engineering and Industrial Aerodynamics，2009，97:96-109.

［9］李正良，石承啟，趙仕興，魏奇科，劉紅軍.復雜體型高層建筑風洞試驗及數值模［J］.土木建筑與環境工程，2009，31(5):69-73.

［10］王小華，朱文芳，何鐘怡.方形鈍體受限繞流的三維數值模擬［J］.計算力學學報，2008，25(5):671-675.

［11］吳文權.流體繞多個鈍體不穩定分離流動數值仿真［J］.工程熱物理學報，1997，19(3):1-8.

［12］陳斌，郭烈錦，楊曉剛.圓柱繞流的離散渦數值模擬［J］.自然科學進展，2002，12(6):964-969.

［13］GERMANO M，PIOMELLI U，MOIN P，CABOT W H.A dynamic subgrid-scale eddy viscosity model［J］.Phys.Fluids A，1991，3(7):1760-1765.

［14］LEONARD B P.A stable and accurate convective modeling procedure based on quadratic upstream interpolation［J］.Comput.Meth.Appl.Mech.Eng，1979，19:59-98.

［15］LAATARAD A H，BENAHMEDB M，BELGHITHC A，QUERE P L.2D large eddy simulation of pollutant dispersion around a covered roadway［J］.Journal of Wind Engineering and Industrial Aerodynamics，2002，90:617-637.

［16］ADINA R ＆ D，Inc.ADINA theory and modeling guide volume III:ADINA CFD ＆ FSI［M］.USA，2005.http://www.adina.Com.

［17］BAETKE F，WERNER H.Numerical simulation of turbulent flow over surface-mounted obstacles with sharp edges and corners［J］.Journalof Wind Engineering and Industrial Aerodynamics，1990(35):129-147.

［18］BEKELE S A，HANGAN H.A comparative investigation of turbulent of the TTU pressure-numerical versus laboratory and full scale results［J］.Wind and Structures，2002，5(2-4):337-346.

［19］AIANNUZAI P S.Artificial wind generation and structural response［J］.Journal of Structural Engineering，ASCE，1987，113(10):2382-2398.

［20］WANG D Y，ZHOU Y，DENG X S，CHEN L.Numerical simulation of fluctuating wind field for high-rise buildings considering spatial correlation［A］.ISSEYE-10［C］，Beijing:Science Press，2008:1618-1624.

［21］JEONG S H，BIENKIEWICZ B.Application of autoregressive modeling in proper orthogonal decomposition of building wind pressure［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，69-71:685-695.

［22］Architechtural Institute of Japan.AIJ recommendations for loads on building［S］.2004.

［23］CASTRO I P，ROBINS A G.The flow around a surfacemounted cube in uniform and turbulent streams［J］.J.Fluid Mech.，1977，79(2):307-335.

［24］GOMES G M，RODRIGUES A M，MENDES P.Experimental and numerical study of wind pressures on irregular-plan shapes［J］.Journal of Wind Engineering and Industrial Aerodynamics，2005，93:741-756.

［25］ZHANG Y Q，HUBER A H，ARYA S P S，SNYDER W H.Numerical simulation to determine the effects of incident wind shear and turbulence level on the flow around a building［J］.J.Wind Eng.Ind.Aerodyn.，1993，46-47:129-134.