劉祖軍，葛耀君，楊詠昕

(同濟大學 橋梁工程系，上海 200092)

目前橋梁風工程中的氣動導數(shù)一般通過風洞節(jié)段模型試驗來獲得，采用的方法主要有自由振動法和強迫振動法。前者通過給模型一個初始位移或脈沖激勵，使其在氣動力和結構阻尼下做自由衰減運動，后者采用機械裝置驅(qū)動模型作可控頻率、振幅的強迫運動，測得模型位移、加速度或表面壓力等時程，通過各種參數(shù)識別方法獲得氣動導數(shù)。自由振動法試驗設備簡單，識別精度低，折減風速范圍小，強迫振動法雖然可以克服這些缺陷，參數(shù)識別精度高，能得到大范圍折算風速下的氣動導數(shù)，但是由于其試驗設備要求高，在工程中應用并不很廣泛［1，2］。隨著計算機硬件能力的提高和計算流體動力學的發(fā)展，采用數(shù)值模擬方法進行氣動導數(shù)識別在斷面選型中應用越來越廣泛。

數(shù)值模擬［3］采用動網(wǎng)格方法識別氣動導數(shù)如果僅控制結構表面一層的網(wǎng)格節(jié)點運動就要求網(wǎng)格劃分得足夠細，結構振幅不易過大，否則容易出現(xiàn)網(wǎng)格節(jié)點間的間距過小，導致網(wǎng)格形成很小的銳角或網(wǎng)格空間位置發(fā)生翻轉而產(chǎn)生錯誤無法計算，而振幅太小不能模擬結構體的真實振動情況，為了兼顧振幅與網(wǎng)格密度，本文采用多層動網(wǎng)格技術來解決這一問題。為考慮湍流影響采用大渦模擬方式(LES)求解N-S方程，通過對流體計算軟件fluent的二次開發(fā)，利用其用戶自定義函數(shù)描述結構的運動狀態(tài)，采用最小二乘法擬合氣動力時程曲線來識別斷面的氣動導數(shù)。

1 基于多層動網(wǎng)格的氣動導數(shù)識別

1.1 流體動力學計算

流體域內(nèi)計算結構在一個振動周期內(nèi)作用在其表面的氣動力是通過將結構表面的流體網(wǎng)格節(jié)點的簡諧振動位移等分的離散到N個時間步上，將第i個時間步上變形后的節(jié)點坐標作為第i+1步的網(wǎng)格文件，每個時間步上求解N－S方程。

湍流［4］的存在使得求解上述N－S方程變得十分困難，通過大量學者的不懈努力已提出較多的湍流模型，較為常見的有:雷諾時均方法(RANS)，直接數(shù)值模擬(DNS)以及介于二者之間的大渦模擬(LES)。由于RANS方法平均的結果是將時空變化的細節(jié)一概抹平，丟失了包括在脈動運動中的大量信息，并且為了使N－S方程封閉，提出了Reynolds應力模型和渦粘模型，這些模型有一定的局限性，都存在著對經(jīng)驗數(shù)據(jù)過分依賴和預報精度較差的特點。DNS采用計算機直接數(shù)值求解三維非定常N－S方程，對湍流瞬時運動直接模擬，可認為是一種精確的方法。主要用做湍流基礎研究，目前僅限于科學研究領域。

大渦模擬［5－7］是放棄了對全部尺度范圍內(nèi)的渦運動進行直接模擬的夢想，改為只對尺度大的渦運動通過數(shù)值方法直接求解N－S方程，對尺度小的渦運動不直接求解而是通過建立模型模擬小渦運動對大渦的影響。其具體實現(xiàn)過程為:首先通過濾波方法將小尺度渦和大尺度渦分離開來，即將流動劃分為低頻可解和高頻不可解(亞格子部分)兩部分。對不可解部分建立亞格子尺度模型來模擬。目前常用的亞格子模型主要可以分為三類:渦粘模型，相似性模型和混合模型。本文CFD計算時采用Smagorinsky-Lilly渦粘形式的亞格子雷諾應力模型其表達式如下:

Smagorinsky模型相當于混合長度形式的渦粘模型，其中CsΔ相當于混合長度，Δ為過濾尺度，Cs為Smagorinsky常數(shù)，Lilly于1987年通過計算建議 Cs≈0.18。

1.2 動網(wǎng)格空間守恒條件

采用動網(wǎng)格方法必須滿足幾何守恒條件，Thomas和Lombard明確提出了用來確定網(wǎng)格速度的法則也稱為空間法則，該法則直接從連續(xù)方程中推得:

式中流體速度為常數(shù)，因此幾何守恒法則表達為:

因此只要是動網(wǎng)格問題都必須保證網(wǎng)格速度滿足幾何守恒法則，否則可能在動網(wǎng)格計算過程中導致人為的質(zhì)量源項。

本文提出的多層動網(wǎng)格方法是通過設定每層動網(wǎng)格做剛體的平動或扭轉運動，相鄰動網(wǎng)格層之間通過彈簧近似方法連接形成變形網(wǎng)格層，變形網(wǎng)格區(qū)，網(wǎng)格只發(fā)生變形沒有網(wǎng)格速度。因此只需證明發(fā)生剛體平動的動網(wǎng)格層滿足幾何守恒條件即可。

如圖1所示的控制體中，w，s，e，n分別為運動前網(wǎng)格的四個邊界，當網(wǎng)格發(fā)生豎向平動時，邊界s平移到s'的位置，n平移到n'位置，由于是剛體平移，故Δy1=Δy2，邊界w和e沒有變化。

圖1 網(wǎng)格運動圖示Fig.1 The motion of mesh

在圖示控制體上，對連續(xù)方程進行離散;

式中:u，v代表笛卡爾坐標系下的水平速度分量和豎向速度分量，n和n+1對應兩個速度時刻，ΔV為控制體體積，假定流體運動速度不變，因此有關流體速度項可以抵消，上式簡化為:

按照前面定義網(wǎng)格做豎向剛體運動時:

代入上式則有:

故每層網(wǎng)格的運動都能保證幾何守恒法則。

1.3 多層動網(wǎng)格的實現(xiàn)

在工程實際問題中經(jīng)常遇到流體與固體相互耦合作用的情況，通常計算域的邊界是運動的，從而使計算域的形狀隨時間變化。動網(wǎng)格技術隨時間變化對流場模型進行更新。在某一時刻，邊界發(fā)生一定位移后，邊界內(nèi)部的網(wǎng)格應該如何確定，已經(jīng)有很多研究成果。較常用的有彈簧近似法，彈性體法，采用Laplace方程控制的方法，其中以彈簧近似方法應用最多。

彈簧近似方法原理簡單［8］(圖2)，使用方便。在實際應用中當邊界位移較大時，容易出現(xiàn)網(wǎng)格節(jié)點間的間距過小，從而導致網(wǎng)格形成很小的銳角(圖3)或網(wǎng)格點的空間位置發(fā)生翻轉(圖4)，導致計算失敗。因此在CFD計算時僅控制結構體表面一層網(wǎng)格節(jié)點運動具有一定的局限性。如果要模擬出流體域中的附面層不同大小的渦時就要將流體的網(wǎng)格劃分得足夠細，一般結構表面的網(wǎng)格在0.00001 m的數(shù)量級，所以結構表面節(jié)點空間位置的變化對其周圍流體單元的形態(tài)影響較大。結構振幅較大，容易導致網(wǎng)格變形過大，而產(chǎn)生錯誤無法計算;振幅太小，不能模擬結構的真實振動情況，振幅與網(wǎng)格密度不能兼顧，因此需要通過多層動網(wǎng)格技術來實現(xiàn)。

多層動網(wǎng)格方法將結構的振動位移按照一定的比例加到結構周圍的N層網(wǎng)格節(jié)點上，一定程度上加大了結構的可移動范圍，使流體網(wǎng)格在N層范圍內(nèi)實現(xiàn)運動。通過設定每層動網(wǎng)格層做剛體的平動或扭轉運動，相鄰動網(wǎng)格層之間通過彈簧近似方法連接形成變形網(wǎng)格層，由于采用多層動網(wǎng)格后相鄰網(wǎng)格的平動速度差值變小，因此變形區(qū)的網(wǎng)格不容易出現(xiàn)網(wǎng)格翻轉或過小尖角，易保證網(wǎng)格質(zhì)量。多層動網(wǎng)格的具體做法是:按照各層間的比例關系給出各層節(jié)點的運動位移(如圖5)，結構表面距離為Ri的第i層網(wǎng)格節(jié)點的位移最大等于結構的振幅即Si=Am，結構表面距離為Rn的第n層網(wǎng)格位移為0即Sn=0，位于這兩層之間的第j層網(wǎng)格節(jié)點位移進行線性插值，這樣就實現(xiàn)了運動結構周圍流體網(wǎng)格的多層運動。這樣既提高了網(wǎng)格可動范圍，又減小了相鄰兩層動網(wǎng)格之間的相對位移，避免了網(wǎng)格位移過大造成的計算失敗。

圖2 動網(wǎng)格彈簧節(jié)點連接Fig.2 Spring node of dynamic mesh

圖3 網(wǎng)格變形后出現(xiàn)過小的尖角Fig.3 Produce small cusp after deformation mesh

圖4 網(wǎng)格變形后出現(xiàn)翻轉Fig.4 Produce flip after deformation mesh

圖5 多層動網(wǎng)格示意圖Fig.5 The multi-layer dynamic mesh

1.4 氣動參數(shù)的識別原理

氣動力可以用氣動導數(shù)來表示，而氣動導數(shù)是系統(tǒng)運動折減頻率的函數(shù)，其函數(shù)關系決定于斷面的氣動外形，一般很難有解析表達式，大多通過節(jié)段模型風洞試驗來識別。氣動力用試驗測到的8個氣動導數(shù)表示為:

理想平板是完全理想化的流線型結構，可以通過Theodorson建立的理論公式得到振動平板所受到的氣動力，從而導出氣動導數(shù)理論值。平板振動時所受到的非定常氣動升力和升力矩可表示為:

式(8)～式(11)中:U為來流速度，B=2b為平板寬度;K=ωB/U為折算頻率;ω為振動頻率;h，a分別為平板豎向位移和扭轉角;

C(k)=F(k)+iG(k)為Theodorson函數(shù);

根據(jù)Theodorson建立的上述自激力公式可以導出平板斷面的8個氣動導數(shù)，并表達為折減頻率函數(shù):

強迫平板單獨作豎向正弦振動或扭轉正弦振動，將結構表面壓力值和摩阻力沿表面積分，可得到在絕對坐標系里的氣動升力和彎矩。然后由得到的氣動力和式(8)－(9)中對應的氣動導數(shù)之系數(shù)，根據(jù)最小二乘法，即可得到相應的氣動導數(shù)。豎向振動時采用上述方法可以提取H，4個氣動導數(shù)，做扭轉運動可以確定另外的4個氣動導數(shù)

2 平板的顫振導數(shù)識別

通過流體分析軟件Fluen［9］的用戶自定義函數(shù)來描述模型的運動狀態(tài)，并利用其大渦模擬的湍流模式進行數(shù)值計算(圖6)。

圖6 數(shù)值計算流程Fig.6 The process of Numerical

平板模型如圖7，寬高比22.5，模型前后緣尖銳，小攻角下流線形特征明顯，流場計算域外邊界為矩形，參考同濟大學TJ－1風洞 尺寸設為2×6 m。

圖7 平板模型圖Fig.7 The model of flat plate

平板豎向振動的振幅0.08m，扭轉幅角3°。入口為速度邊界條件，入口風速為U=5 m/s，出口為壓力邊界條件，表壓設為0，參考壓力點在出口邊界中心處;上下邊界設為對稱邊界，以加快計算速度，時間步長設為0.001 s。建立二維計算模型，網(wǎng)格線與物面平行，周向網(wǎng)格線在靠近物面附近加密，離開物面間距逐漸加大。交于物面的網(wǎng)格線在模型前后緣位置均加密.網(wǎng)格數(shù)為3.02萬。

模型從零攻角的平衡位置開始運動.圖8、圖10是低折算風速(Vr=5)下的運動，圖9、圖11是高折算風速(Vr=20)下的運動.平板分別作純豎向振動和純扭轉振動，氣動力數(shù)據(jù)在2個周期后開始采集。從圖8－圖11可知，在這兩種振動情況下，升力系數(shù)和升力矩系數(shù)均表現(xiàn)為對無量綱時間零均值的一次穩(wěn)態(tài)諧波曲線特性，升力、升力矩與位移之間的良好的線性關系表明流場與平板一起進入諧振動狀態(tài)，流場沒有出現(xiàn)不穩(wěn)定情況。

在不同的折算風速下，分別計算平板做純豎向振動、純扭轉振動的氣動力，當氣動力的變化趨勢變得平穩(wěn)時，開始采集氣動力數(shù)據(jù)，并用最最小二乘法提取氣動導數(shù)。本文比較了計算得到的氣動導數(shù)與通過Theodorsen理論［10］推導出的氣動導數(shù)之間的差別，見圖12，從圖中各顫振導數(shù)值的比較來看，當折減風速較低時，各導數(shù)值吻合較好，隨著折減風速增大，各導數(shù)值偏差增大，其中和的誤差相對大些，但均與通過Theodorsen理論導出的氣動導數(shù)值有合理的一致性，其中幾個關鍵導數(shù)及吻合的很好。

圖8 低折算風速下純豎向運動產(chǎn)生的氣動升力系數(shù)和升力矩系數(shù)時程Fig.8 Lift and moment coefficients Vs time on heaving vibration at lower reduce wind speed

圖9 高折算風速純豎向運動產(chǎn)生的氣動升力系數(shù)和升力矩系數(shù)時程Fig.9 Lift and moment coefficients Vs time on heaving vibration at higher reduce wind speed

圖10 低折算風速純扭轉運動產(chǎn)生的氣動升力系數(shù)和升力矩系數(shù)時程Fig.10 Lift and moment coefficients Vs time on pitching vibration at lower reduce wind speed

圖11 高折算風速純扭轉運動產(chǎn)生的氣動升力系數(shù)和升力矩系數(shù)時程Fig.11 Lift and moment coefficients Vs time on pitching vibration at higher reduce wind speed

圖12 平板氣動導數(shù)Fig.12 Flutter derivatives of flat plate

3 結論

本文通過對流體計算軟件fluent的二次開發(fā)，利用其提供的用戶自定義函數(shù)模塊，結合動網(wǎng)格技術實現(xiàn)了平板的氣動導數(shù)識別。由于固體模型在流場中運動受到流體網(wǎng)格尺寸限制，當模型運動幅度較大時，常會造成網(wǎng)格的變形過大，導致計算無法進行。針對這一問題，本文采用了多層動網(wǎng)格方法較好地解決了這一問題。能較好的兼顧模型振幅和網(wǎng)格密度。采用本文方法所獲得的計算結果與通過Theodorsen理論導出的平板氣動導數(shù)具有良好的一致性。

［1］Scanlan R H，Tomko J J.Airfoil and bridge flutter derivatives［J］.Journal of Engineering Mechanics Division，1971，97:1717－1737.

［2］Yang Yongxin，Ge Yaojun，Xiang Haifan.Investigation on fluttermechanism oflong-span bridgeswith 2d-3DOF method. WINDANDSTRUCTURES，2007，10(5):421－435.

［3］Chen Z，Hong L，Xiaofeng W.Full function numerical methodfor flow over a self-excited vibrating Body［J］.Tsinghua Science and Technology，1999，4(4):1688 －1691.

［4］Fasel H F，Seidel J，Wernz S.A Methodology for simulations ofcomplex turbulentflows ［J］. Journalof Fluids Engineering，2002，124(5):933 －942.

［5］Kalttenbach J H，Schumann U，Gerz T.Large eddy simulation of turbulent diffusion in stably-stratified flow［J］.Fluids Mech，1994，280:1 －40.

［6］Kondo K，Murakami S，Mochida A.Generation of velocity fluctuations for inflow boundary condition of LES，［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，67:51－64.

［7］Tetsuro Tamura，Yoshiyuki Ono.LES analysis on aeroelastic instability of prisms in turbulent flow［J］.Journal of Wind Engineering and IndustrialAerodynamics，2003，91:51 －64.

［8］郭 正，李曉斌.用非結構動網(wǎng)格方法模擬有相對運動的多體繞流［J］.空氣動報，2001，19(3):28－35.

［9］王福軍.計算流體動力學分析——CFD軟件原理與應用［M］.北京:清華大學出版社，2004.

［10］Theodorsen T.General theory of aerodynamic instability and the mechanism of flutter［R］.NACA Report，1935.