祝志文 林君福 唐意 王欽華

摘要：采用任意拉格朗日一歐拉（Arbitrary Lagrangian-Eulerian，ALE）描述法，推導(dǎo)了二維流變區(qū)域不可壓流體流動的控制方程?；趧傂詳嗝孢吔缣攸c，簡化得到了薄平板強迫振動繞流場控制方程，采用有限差分法開展了薄平板大幅強迫扭轉(zhuǎn)振動的CFD模擬。研究發(fā)現(xiàn)，薄平板小幅扭轉(zhuǎn)振動氣動力呈線性特性，但大幅扭轉(zhuǎn)振動氣動力非線性顯著，且氣動力非線性隨來流折算風(fēng)速的增大而變得顯著。從一個周期內(nèi)氣流做功來看，雖然大幅扭轉(zhuǎn)振動時薄平板轉(zhuǎn)動軸到后緣是耗散氣流能量的氣動穩(wěn)定區(qū)，但當(dāng)折算風(fēng)速較高后，上游側(cè)的半個薄平板表面變?yōu)闅鈩硬环€(wěn)定區(qū)，此時氣流對整個薄平板所作總功由負變?yōu)檎∑桨鍖臍饬髦形漳芰?。研究認為，當(dāng)折算風(fēng)速較高時，大幅扭轉(zhuǎn)振動的薄平板將進入氣動不穩(wěn)定狀態(tài)。

關(guān)鍵詞：渦激振動;薄平板;大幅振動;氣動不穩(wěn)定;任意Lagrange-Eulerian描述

中圖分類號：TU311.3 文獻標(biāo)志碼：A 文章編號：1004-4523（2020）01-0012-12

DOI：10.16385/j.cnki.issn.1004-4523.2020.01.002

引言

處于大氣邊界層中的大跨度橋梁受自然風(fēng)的作用。在某些條件下，橋梁或構(gòu)件可能出現(xiàn)大幅度的振動，如拉索、吊桿和主梁的大幅渦激振動，風(fēng)雨振和尾流馳振、鄰近或進人顫振狀態(tài)的主梁大幅顫振等。這些大幅振動都是構(gòu)件表面作為流體邊界的流動域發(fā)生顯著變形的流體一固體耦合作用問題;對鈍體外形的橋梁構(gòu)件，或構(gòu)件因大幅振動產(chǎn)生較大的相對攻角效應(yīng)，其氣動力系統(tǒng)可能會因這種大幅振動呈現(xiàn)顯著的非線性。

在連續(xù)介質(zhì)力學(xué)中，流體微元是微觀上遠大于分子運動平均自由程而宏觀上遠小于所研究的流動域的一團流體，并在宏觀上將流體微元看作一個質(zhì)點。流體力學(xué)描述流體微元集合的運動狀態(tài)，并通過三種方法描述流體的運動，其中一種是拉格朗日法（Lagrangian），這種方法跟蹤流場每一質(zhì)點的運動，并獲得流場流動參數(shù)隨時問的變化，其特點是計算網(wǎng)格與流體質(zhì)點固定并同步運動，不存在質(zhì)點與網(wǎng)格問的相對運動。與歐拉描述法相比，這種方法不僅可跟蹤流體質(zhì)點的運動軌跡，而且能準確描述流體的移動界面。但當(dāng)流動出現(xiàn)大變形時，將直接導(dǎo)致計算網(wǎng)格的畸變，計算失敗。歐拉法（Euleri-an）通過將感興趣的流動區(qū)域離散為位置固定而數(shù)量有限的網(wǎng)格點，進而描述流動參數(shù)的空間分布隨時問的變化。因其計算網(wǎng)格不隨物體形位變化而改變，所以能容易地處理流動區(qū)域的大變形，但描述流變界面的運動需要引入非常復(fù)雜的映射關(guān)系，可能導(dǎo)致較大的計算誤差。任意拉格朗日一歐拉法就是將拉格朗日法和歐拉法組合起來，取長補短。其特點是描述流體運動的空問點位置既不固定，也不隨流體運動，而是以某種方式隨時問變化，這種質(zhì)點的變化方式能反映流動域的實際變化，因而可準確地描述流體的運動界面，并維持單元的合理形狀，因而可解決只用歐拉法或只用拉格朗日法解決不了的復(fù)雜運動邊界的大變形流動問題。

橋梁構(gòu)件大幅振動中，構(gòu)件周圍空氣所占據(jù)的實際空問一般是隨時問變化的，并具有運動的界面?？梢?，離散點空問分布固定的歐拉方法顯然不適合描述此類問題，因流動區(qū)域扭曲而引起計算網(wǎng)格的畸形也會導(dǎo)致計算失敗，因而合理描述應(yīng)該采用任意拉格朗日一歐拉法。文獻[11-12]通過CFD模擬獲得的振動剛性模型上的氣動力和強迫振動位移時程，建立氣動力離散時問氣動模型，開展了顫振導(dǎo)數(shù)識別，但模型豎彎和扭轉(zhuǎn)振動的幅值很小。文獻開展了薄平板豎彎和扭轉(zhuǎn)大幅振動CFD模擬和氣動力小波分析，發(fā)現(xiàn)在氣動力時程中出現(xiàn)了比薄平板振動頻率高的二階高頻成分，也即大幅振動的薄平板氣動力為非線性，但因CFD基于Fluent開展，而Fluent的動網(wǎng)格即使在均勻網(wǎng)格上也僅有一階精度，因此研究結(jié)果的準確性需要驗證。文獻[14]基于Volterra理論開展了薄平板非線性氣動力系統(tǒng)識別和CFD模擬，通過豎彎振動幅值為1%和5%薄平板寬度的薄平板強迫振動CFD分析，發(fā)現(xiàn)氣動力的非線性效應(yīng)并不顯著，原因可能是薄平板的振動幅度偏小。下面就橋梁氣彈特性分析常采用的薄平板為研究對象，針對二維不可壓流動特征，首先推導(dǎo)能適應(yīng)任意流變區(qū)域的流體運動控制方程，再根據(jù)剛性斷面振動的特點，給出合適的空問點運動方式和相應(yīng)的流體運動控制方程。并以薄平板為例研究其大幅振動的氣動特征。

1二維流體運動的任意拉格朗日-歐拉描述方程

基于式（18）-（22）可在二維計算域上實現(xiàn)對任意流變區(qū)域內(nèi)流體的運動描述，也即為描述任意流變區(qū)域不可壓流動的基本方程，其中映射關(guān)系式（1）需根據(jù)所研究的流動區(qū)域的具體變形形式來確定。當(dāng)物理域變形特征簡單時，對應(yīng)的映射關(guān)系可簡化，此時，式（18）-（22）可能會得到明顯簡化，比如橋梁斷面的氣動彈性問題，因橋梁斷面的絕對剛性假設(shè)，使得映射關(guān)系明顯簡化，如下文所述。

2薄平板強迫振動繞流場控制方程

式（1）和（18）-（22）能描述任意變形物理域的流體流動，比如帶移動邊界的流體動力學(xué)問題。移動邊界可為內(nèi)流場或外流場流體外部邊界的變形或移動，也可為外流場繞流物體在流體力作用下其形位的改變，這些均能導(dǎo)致與繞流物體接觸的流體邊界的移動和變形。在進行橋梁斷面一類繞流場模擬時，由于作用在橋梁斷面上的氣動壓力和摩阻力很小，因此橋梁斷面的變形很小，可近似將橋梁斷面看作剛體。這樣，橋梁斷面繞流場內(nèi)邊界的運動可用剛性橋梁斷面的整體運動來描述，也即采用剛體運動的描述方法。在橋梁氣動彈性分析中，通常用橋梁斷面的豎彎和繞某一參考軸的轉(zhuǎn)動來描述橋梁斷面上任意點的運動。

對振動的剛性斷面，如果其計算域的遠場邊界不變，將能恒定地給出遠場邊界條件，但由于計算域內(nèi)、外邊界之問的物理區(qū)域在不斷發(fā)生變化，如果進行數(shù)值模擬，則需在每一時問步形成一次網(wǎng)格，如果計算的時問步數(shù)量很大，由網(wǎng)格生成所帶來的附加計算量將很大，且難以保證每一步上重新生成的網(wǎng)格的高質(zhì)量。相比較，如果在離剛性斷面足夠遠處劃定外邊界，并對計算域劃分高質(zhì)量網(wǎng)格，數(shù)值模擬時將計算域網(wǎng)格和剛性斷面一起同步振動，就可只在計算開始時形成一次網(wǎng)格，而在每一時問步上給出相應(yīng)的邊界條件，如此處理，不僅使得數(shù)值模擬的網(wǎng)格質(zhì)量明顯提高，也能省去每一時問步上網(wǎng)格的重新生成。這樣，需求解的物理區(qū)域只隨時問發(fā)生位置改變，形狀并不改變。由于計算域事先給定，計算域和物理域的對應(yīng)關(guān)系也能明確。

取絕對坐標(biāo)系OXY平面內(nèi)任意點為物理域之參考點，剛性斷面的振動方式為相對于參考軸的平動和繞參考點的轉(zhuǎn)動。如取計算域的坐標(biāo)原點在物理域的映射點為參考點，并在此參考點上建立計算域的相對坐標(biāo)系oxy，當(dāng)相對坐標(biāo)系隨剛性斷面作同步運動時，計算域Dcomp（x，y）和物理域Dphys（x，y，t）存在下述映射關(guān)系

通過選擇合適的時問和空問步離散格式，對式（36）-（38）和式（41）-（42）進行離散，在合理計算域網(wǎng)格基礎(chǔ)上，基于給定的初始和邊界條件下，就可開展剛性斷面強迫振動的數(shù)值模擬。

3薄平板扭轉(zhuǎn)振動繞流場控制方程的數(shù)值求解

3.1CFD數(shù)值實現(xiàn)

對圖6所示的長度為26的無厚度扭轉(zhuǎn)振動剛性薄平板（為便于查看，圖示有厚度），其繞流控制方程式本文采用非定常不可壓流計算的二階Projec-tion算法，即將上述控制方程分裂為依次求解動量方程得到中問速度場，由中問速度場求解壓力場的Poisson型方程，再由獲得的壓力場更新中問速度場，得到修正的速度場方程，如此反復(fù)求解，反復(fù)修正，并依據(jù)修正速度場在連續(xù)方程式（37）中產(chǎn)生的殘余值的大小來判斷當(dāng)前時問步流場變量的收斂情況。當(dāng)修正速度對連續(xù)方程的殘余值滿足收斂條件時，求解進入下一時問步。

由二階Projection算法得到的偏微分方程，空問離散采用基于交錯網(wǎng)格的有限差分法。在關(guān)于中問速度場的動量方程求解上，左端動量項采用二階外插的Adams-Bashforth格式，右端黏性項采用二階中心差分格式;壓力場方程左端項采用二階中心差分格式;時問導(dǎo)數(shù)采用二階中心格式，得到的差分方程均采用三對角矩陣算法求解，并通過Fortran程序編程實現(xiàn)。

因剛性薄平板厚度為零，計算域網(wǎng)格采用直角坐標(biāo)系下的結(jié)構(gòu)化網(wǎng)格，物面上下第一個網(wǎng)格高度為0.001b，并沿垂直物面方向網(wǎng)格尺度逐漸增大，但增長率不大于1.05，以捕捉邊界層的復(fù)雜流動，并適度控制網(wǎng)格規(guī)模和非定常計算量。薄平板沿流向均勻劃分32個網(wǎng)格，但從薄平板前緣往計算域入口，以及從薄平板后緣往計算域出口，采用不大于1.05的網(wǎng)格增長率逐漸增大網(wǎng)格。計算域總網(wǎng)格數(shù)為168×168，圖7是計算域整體和薄平板周圍的網(wǎng)格劃分，其中紅色方框內(nèi)有模擬的薄平板。

設(shè)圖6所示的薄平板繞垂直O(jiān)XY平面的薄平板轉(zhuǎn)動軸轉(zhuǎn)動，振動方式為正弦振動，即

式中fa為薄平板扭轉(zhuǎn)振動的頻率;a（t）和ao分別為扭轉(zhuǎn)振動位移角和位移角幅值（°），以薄平板前緣抬頭為正向，如圖6所示。

從控制方程式（34）-（36）和式（39），（40）的推導(dǎo)可知，網(wǎng)格系統(tǒng)和薄平板剛性固定，同步振動。因此在每一時問步，基于式（41）可得到網(wǎng)格和薄平板的同步振動速度;薄平板表面采用無滑移條件，因此薄平板物面流體的運動速度等于所在薄平板表面扭轉(zhuǎn)速度;對計算域的人口和上、下邊界，因來流平行于絕對參考系OXY的OX軸，需在每一時問步上分別投影到運動參考系oxy的兩個坐標(biāo)軸方向，成為計算域?qū)?yīng)邊界的速度邊界條件;計算域出口設(shè)置為紐曼邊界條件。對壓力邊界條件，在薄平板物面和計算域人口、上下邊界也采用紐曼邊界條件。式中B=26?？紤]作用在振動薄平板上的氣動力對雷諾數(shù)不敏感，從降低不可壓流非定常迭代求解的計算量考慮，本文所有模擬的雷諾數(shù)均為1200。

3.2CFD方法驗證

薄平板顫振導(dǎo)數(shù)有Theodorsen解析解，通常用于檢驗氣彈分析CFD方法和程序的正確性。對單自由度扭轉(zhuǎn)強迫振動的薄平板，其非定常氣動力與相應(yīng)顫振導(dǎo)數(shù)的表達式為：

分別開展不同折算風(fēng)速的CFD模擬，得到作用在扭轉(zhuǎn)振動薄平板上的氣動力，利用式（46）和（47）基于最小二乘法識別薄平板的4個顫振導(dǎo)數(shù)，并將結(jié)果和Theodorsen解析解作了比較，如圖8所示。可見二者具有一致的趨勢線，高折算風(fēng)速二者出現(xiàn)較小偏差，可能是Theodorsen解析解是基于無黏、無流動分離的勢流理論，但本文流動是真實的有黏空氣，且薄平板振動會導(dǎo)致流動分離，而黏性效應(yīng)和流動分離可能隨折算風(fēng)速的增大而增大。研究表明，本文CFD方法能較準確得到振動薄平板上的氣動力，數(shù)值方法是正確的。

4薄平板大幅扭轉(zhuǎn)強迫振動的能量特征

4.1振動能量的定義

在來流中繞扭轉(zhuǎn)軸作扭轉(zhuǎn)強迫振動的薄平板，定義薄平板表面的壓力系數(shù)為

為對比薄平板大幅振動的氣動力特征，下面給出了薄平板小幅扭轉(zhuǎn)振動（振幅0.2°）的氣動力時程。圖9（a）和（b）分別是折算風(fēng)速Vr=4和24時，一個周期內(nèi)的氣動力系數(shù)時程曲線。可見無論是升力還是扭矩系數(shù)時程，其曲線具有明顯的諧波特性，其波形頻率等于扭轉(zhuǎn)強迫振動頻率，可知當(dāng)薄平板小幅扭轉(zhuǎn)振動時，其氣動力系統(tǒng)具有良好的線性特性;另外，從不同折算風(fēng)速下的曲線來看，升力和扭矩問的相位差隨折算風(fēng)速不同而有變化。

圖10為小幅扭轉(zhuǎn)振動薄平板在一個完整周期內(nèi)4個關(guān)鍵時刻的壓力等高線，對應(yīng)從平衡位置出發(fā)后達到最大正攻角位置（圖10（a）），然后回到平衡位置（圖10（b）），再達到最大負攻角位置（圖10（c）），最后回到零攻角的平衡位置（圖10（d））?？梢姀囊粋€狀態(tài)到另外一個狀態(tài)，壓力等高線呈現(xiàn)有規(guī)律地變化，薄平板前緣為高風(fēng)壓區(qū)，后緣為低風(fēng)壓區(qū)。等高線分布規(guī)律的改變反映了薄平板的扭轉(zhuǎn)振動方向。

圖11（a）和（b）分別為在一個扭轉(zhuǎn)振動周期內(nèi)，氣流在振動薄平板上表面各點、以及整個薄平板所作的無量綱總功?？梢姳∑桨迩熬墸↙.E.）局部為穩(wěn)定區(qū)外，其后到轉(zhuǎn)動軸（A.R.）為不穩(wěn)定區(qū)，也即氣流對薄平板作正功，薄平板從氣流中吸收能量，因此為激勵區(qū);扭轉(zhuǎn)軸到薄平板后緣（T.E.）氣流做負功，因而是耗散氣流能量的穩(wěn)定區(qū)，因此薄平板上激勵區(qū)與穩(wěn)定區(qū)共存。圖11（b）給出了Vr=4-40，一個扭轉(zhuǎn)周期內(nèi)氣流對薄平板所作的總功?？梢娫谒嬎愕恼鬯泔L(fēng)速范圍內(nèi)，一個周期內(nèi)氣流對扭轉(zhuǎn)振動的薄平板均做負功，也即振動薄平板消耗氣流的能量，因此小幅扭轉(zhuǎn)振動的薄平板是穩(wěn)定的。

4.3大幅扭轉(zhuǎn)振動

強迫薄平板發(fā)生大幅度的扭轉(zhuǎn)振動，如扭轉(zhuǎn)振幅角為20°，圖12（a）-（d）分別為折算風(fēng)速Vr=4，24，44和68時，一個周期內(nèi)的氣動力系數(shù)時程曲線。與圖9相比，可見隨著折算風(fēng)速的提高，升力和扭矩系數(shù)時程逐漸失去諧波特性，特別是高折算風(fēng)速。由于扭轉(zhuǎn)振動為正弦強迫位移，可見當(dāng)薄平板大幅扭轉(zhuǎn)振動時，其氣動力系統(tǒng)不再為線性，且隨著折算風(fēng)速的提高，其非線性特性越來越顯著。

圖13為不同折算風(fēng)速下大幅扭轉(zhuǎn)振動的薄平板，在一個振動周期內(nèi)氣流在振動薄平板上表面各點所作的無量綱總功。與圖11對比可見，在計算的所有折算風(fēng)速范圍上，轉(zhuǎn)動軸到薄平板后緣氣流做負功，因而是薄平板上耗散氣流能量的穩(wěn)定區(qū)，且在低折算風(fēng)速（Vr=4-16）下，薄平板表面的能量消耗特征與小幅度扭轉(zhuǎn)振動類似。但當(dāng)折算風(fēng)速高于20后，迎風(fēng)側(cè)的半個薄平板上每個點在一個周期內(nèi)氣流做的功均為正，即為不穩(wěn)定區(qū)，氣流對薄平板作正功，薄平板從氣流中吸收能量，因此薄平板上的激勵區(qū)與穩(wěn)定區(qū)共存，而一個周期內(nèi)氣流對整個薄平板表面所有點的總功之和決定了一個周期內(nèi)氣流對薄平板所作的總功，該總功的正負反映了扭轉(zhuǎn)振動薄平板的氣動穩(wěn)定性。

圖14給出了一個扭轉(zhuǎn)振動周期內(nèi)氣流對薄平板所作總功隨折算風(fēng)速的變化，可見當(dāng)折算風(fēng)速小于16時，一個周期內(nèi)氣流對薄平板所作總功為負，也即薄平板消耗氣流的能量;但當(dāng)折算風(fēng)速大于20后，氣流對薄平板所作的總功由負變?yōu)檎?，此時一個周期內(nèi)扭轉(zhuǎn)振動薄平板從氣流中吸收能量，這表明大幅度扭轉(zhuǎn)振動的薄平板已進入氣動不穩(wěn)定狀態(tài)。

5結(jié)論

本文推導(dǎo)了二維流變區(qū)域的任意拉格朗日一歐拉描述方程，基于剛性薄平板簡化得到了薄平板強迫振動繞流場控制方程，采用數(shù)值方法開展了薄平板大幅強迫扭轉(zhuǎn)振動的氣動特征研究，可得到下述結(jié)論：

1）基于任意拉格朗日一歐拉描述的二維流體運動方程，通過建立計算域和物理域的隨時問變化的映射關(guān)系，可在不隨時問變化的二維計算域上，采用歐拉描述法實現(xiàn)對二維物理域上任意時變區(qū)域的流體流動的描述。對薄平板扭轉(zhuǎn)強迫振動，相應(yīng)的映射關(guān)系和控制方程明顯簡化。

2）采用有限差分法開展了薄平板扭轉(zhuǎn)振動的CFD模擬。從小幅扭轉(zhuǎn)與大幅扭轉(zhuǎn)振動氣動力時程的對比發(fā)現(xiàn)，大幅扭轉(zhuǎn)振動氣動力呈現(xiàn)非線性特征，折算風(fēng)速越高，氣動力的非線性特征越顯著。

3）大幅扭轉(zhuǎn)振動時，薄平板轉(zhuǎn)動軸到后緣是能量耗散區(qū)。但當(dāng)折算風(fēng)速高于20后，迎風(fēng)側(cè)的半個薄平板從氣流中吸收能量，為氣動不穩(wěn)定區(qū);扭轉(zhuǎn)振動一個周期內(nèi)氣流對整個薄平板作的總功由負變?yōu)檎?，此時薄平板從氣流中吸收能量，表明大幅扭轉(zhuǎn)強迫振動的薄平板已進入氣動不穩(wěn)的狀態(tài)。