趙張峰,鄧洪洲

(同濟大學 建筑工程系,上海 200092)

空氣動力分析中動網(wǎng)格技術的數(shù)值阻尼

趙張峰,鄧洪洲*

(同濟大學 建筑工程系,上海 200092)

利用FLUENT進行空氣動力分析，并采用常加速度Newmark法計算物體運動參數(shù)時，動網(wǎng)格宏模塊由于數(shù)據(jù)傳遞方式的限制，修改了常加速度Newmark法的原有算法，動網(wǎng)格宏模塊的軟件算法其有效性存在質疑。針對以上問題，首先給出數(shù)值算例以顯示軟件算法的缺陷特征，提出軟件算法會引入數(shù)值阻尼的假定，而后通過數(shù)學手段證明數(shù)值阻尼的存在，并給出數(shù)值阻尼的理論計算公式。之后，給出算例，得出算例的數(shù)值阻尼，并利用數(shù)值阻尼修正軟件算法，同時驗證理論公式的有效性。最后給出文章理論的工程應用。

FLUENT；動網(wǎng)格；Newmark法；數(shù)值阻尼

0 引 言

FLUENT是常用的空氣動力分析軟件，軟件進行流固耦合分析的計算步驟如下：在每一個時間步內，先求解流體控制方程，得到速度場、壓力場以及作用于固體的升力和阻力，通過UDF提取升力和阻力，并結合動力方程計算固體的振動響應，將振動響應中的速度項通過FLUENT動網(wǎng)格宏模塊DEFINE_CG_MOTION傳遞到FLUENT計算數(shù)據(jù)中，速度與時間步之積即為固體近似的位移響應變化，根據(jù)位移響應變化確定固體下個時間步的空間坐標，同時更新固體附近的網(wǎng)格劃分，進行下一個時間步的計算(以上參考FLUENT用戶自定義函數(shù)手冊[1])。

由文獻[2-4]可知，在FLUENT中流固耦合的計算方法應用廣泛，同時可以采用多種數(shù)學方法求解運動微分方程(新型顯示積分法[5]、四階Runge-Kutta、Newmark-β算法等)，本文著重對采用常加速度Newmark法求解運動微分方程的流固耦合計算方法進行分析。

常加速度Newmark法獲取固體的振動響應，振動響應包括位移、速度、加速度，其計算結果是穩(wěn)定的，同時具有二階精度[6]，因此在時間步長較小的條件下，可以認為這三個振動響應是精確的，然而在FLUENT中采用動網(wǎng)格宏模塊DEFINE_CG_MOTION將速度項傳遞到計算數(shù)據(jù)中，通過速度與時間步的乘積確定位移響應變化，從而確定新時間步的位移響應，這個位移響應與常加速度Newmark法直接計算得出的位移響應不存在等價性。因此，在FLUENT中，用速度響應結合時間步長計算得出的新位移響應來替代用常加速度Newmark法直接計算得出的位移響應，這樣的做法將改變常加速度Newmark法原有的計算過程，是一種近似于常加速度Newmark法的算法(下文將這種算法統(tǒng)一命名為軟件算法，同時常加速度Newmark法簡稱為Newmark法)，需要對軟件采用的算法流程進行分析，以確定算法流程的可靠性。

本文首先給出數(shù)值算例，提出軟件采用的算法流程會引入數(shù)值阻尼的假定。而后通過數(shù)學手段證明數(shù)值阻尼的存在，并給出數(shù)值阻尼的理論計算公式。之后，給出算例，得出算例的數(shù)值阻尼，并利用數(shù)值阻尼修正軟件采用的算法，同時驗證理論公式的有效性。最后，給出文章理論的工程應用。

1 Newmark算法與軟件算法算例比較

1.1 算例信息

文獻[7]闡述了輸電塔典型節(jié)點鋼管構件渦激振動的研究。其中鋼管構件直徑70 mm，壁厚3.5 mm，鋼管桿件長度在2354 mm～3766 mm之間，自振頻率在7.06 Hz～13.36 Hz之間。在風洞中改變來流風速的大小，對鋼管在空氣中渦激振動的現(xiàn)象進行研究。

由于輸電塔鋼管構件渦激振動現(xiàn)象是在微風下產(chǎn)生的，來流風速較小，其振動由一階自振頻率控制，因此采用單自由度體系進行分析。其力學模型的數(shù)學表達式為：

式中：k為物體剛度；c為阻尼系數(shù)；p(t)為外荷載；x為物體的位移。

為方便分析，對構件的參數(shù)進行簡化。物體的質量取為m=5.74 kg，自振頻率取為10 Hz，則k取為22660.6 N/m，空氣密度ρ=1.225 kg/m3，來流風速為3.2 m/s，迎風面積為0.07 m2，阻尼比取ξ為0.015。

由于FLUENT計算量較大，用C語言模擬軟件算法的過程，程序借鑒于文獻[8]。C語言計算結果經(jīng)過驗證與FLUNET穩(wěn)定后的計算結果一致，因此下文采用C語言對FLUENT軟件的算法流程以及Newmark法的算法流程進行編譯計算，并對結果進行比較分析。

根據(jù)結構動力學理論[9]，當阻尼比不為0時，物體振動振幅公式為：

式中：p0是荷載最大值，取為0.307 N；β為荷載頻率與固有頻率比，取為0.97。

不同阻尼比條件下，計算結果如表1所示。

表1 不同阻尼比對應振幅Table 1 Amplitude of different damping ratios

1.2 有阻尼時Newmark算法與軟件算法比較

時間步長為0.001 s，計算9000個時間步，阻尼為0.015。圖1為阻尼比為0.015時，Newmark算法和軟件算法的位移時程圖。采用Newmark算法時，穩(wěn)定的位移最大值為2.07×10-4，與理論值的誤差在0.5%以內，可以認為Newmark算法是精確的。采用軟件算法時，穩(wěn)定的位移最大值為1.60×10-4，與理論值的誤差在22.3%，軟件算法計算誤差很大。

1.3 無阻尼時Newmark算法與軟件算法比較

在有阻尼條件下，軟件算法比Newmark算法的位移幅值要小的多。結合文獻[10]中對數(shù)值阻尼的定義(這種阻尼并不是材料本身的內阻尼，而是由于數(shù)值計算方法產(chǎn)生的“數(shù)值阻尼”，或稱為“算法阻尼”)，筆者提出在軟件算法中引入數(shù)值阻尼的假定。為驗證以上假定，筆者將分析無阻尼條件下的位移時程圖(見圖2)。

為了顯示阻尼經(jīng)典的衰減特征，在無阻尼算例中，外力為0，初始位移為0.0005 m。

在圖2中，采用Newmark算法時，位移沒有隨著時間出現(xiàn)衰減的現(xiàn)象，證明Newmark算法是不會引入數(shù)值阻尼的。然而在采用軟件算法時，阻尼為0的條件下，位移隨著時間出現(xiàn)了衰減的現(xiàn)象，同時根據(jù)阻尼比與位移關系式[9]：

式中：xn、xn+m是物體位移按時間等間距排列后序號為n和n+m的數(shù)值，ξ為阻尼比，Δt為時間間隔，T為自振周期。

根據(jù)式(3)計算得阻尼比為1.57%。計算得到數(shù)值阻尼后，將初始阻尼與數(shù)值阻尼疊加，計算疊加后阻尼比的理論位移值(根據(jù)式(2))，并與軟件算法的數(shù)值算例值進行比較，如表2所示。

表2 考慮數(shù)值阻尼后理論位移與數(shù)值算例位移比較Table 2 Comparison of theoretical displacement and numerical displacement considering numerical damping

在表2中總阻尼為初始阻尼與數(shù)值阻尼之和?？傋枘岜雀鶕?jù)式(2)計算得理論位移，并與算例位移進行比較。根據(jù)表2可知，理論位移與數(shù)值算例位移擬合性很好，誤差在0.9%左右。

因此數(shù)值算例驗證了軟件算法引入數(shù)值阻尼的假定，而數(shù)值阻尼的最終確定需要用數(shù)學方法進行理論證明。

2 算法流程說明

2.1 常加速度Newmark法計算流程說明

根據(jù)文獻[11]給出的常加速度Newmark法采用的基本假定：

為求解t+Δt時刻位移、速度和加速度，引入t+Δt時刻的平衡方程：

式中：m、c、k分別為結構的質量、阻尼系數(shù)、剛度，pt+Δt為t+Δt時刻的外荷載。

5) 此時t+Δt時刻的位移、速度和加速度均為已知，可以進行下一個時間步未知參數(shù)的求解。

2.2 軟件算法計算流程說明

圖3兩種算法的流程比較，從中可以看出，軟件算法在Newmark算法的基礎上，根據(jù)式(7)用速度值對位移值進行了修正。

3 軟件算法數(shù)值阻尼理論證明

3.1 軟件算法位移遞推關系式求解

借鑒文獻[12]中Newmark法解的穩(wěn)定性中相關的推導方式，進行軟件算法的遞推關系式推導。為了避免出現(xiàn)參數(shù)表示的歧義，將臨時xt+Δt用at+Δt表示。由前面可知數(shù)值阻尼與初始阻尼比沒有關系，同時外荷載僅僅影響振幅的大小，不會改變數(shù)值阻尼，因此在求解通項公式的過程中忽略初始阻尼系數(shù)以及外荷載，簡化平衡方程式如式(8)。將簡化平衡方程式與Newmark算法基本假定式(9)(10)以及修正公式(11)聯(lián)立，如下：

將式(8)代入式(9)(10)可得式(12)(13)：

將式(13)代入式(12)可得式(14)：

將式(14)代入式(11)可得式(15):

式(14)為t+Δt時刻的速度表達式，易知t時刻的速度表達式為式(16)：

將式(16)代入式(13)可得式(17):

at+Δt=xt+at-xt-Δt-0.25w2(at-Δt+at)Δt2-

由式(15)易得式(18)(19):

將式(18)(19)代入式(17)可得式(20):

由式(20)易得式(21):

將式(20)(21)代入式(15)可得式(22):

式(22)即為軟件算法的遞推關系式。

3.2 軟件算法位移通項式求解

根據(jù)軟件算法位移的遞推關系式，利用特征值[13]的方法求解位移的通項式。

由遞推關系式(22)可知：

式中：a=[1/2+1/8(w2Δt2)]-1[1-3/8(w2Δt2)]，b=[1/2+1/8(w2Δt2)]-1[-1/2-1/8(w2Δt2)]=-1，c=[1/2+1/8(w2Δt2)]-1[1/8(w2Δt2)]

求解矩陣A，得矩陣A的三個特征值，記為λ1、λ2、λ3，對應的特征向量分別為V1、V2、V3，則可以推得式(24):

式中，ci為待定參數(shù)，可以根據(jù)初始值確定。

由式(24)可以得出軟件算法的位移通項式(25):

式中，v1i為特征向量Vi第一行的數(shù)值。

由數(shù)值算例的求解可知，結構振動形式具有周期性，且振幅隨著時間的推移出現(xiàn)衰減性，因此特征值λ1和λ2必定為一對共軛復數(shù)，同時c1v11=c2v12(保證位移不為虛數(shù)，消除虛數(shù)項)，記λ1,2=α±βi，并假定物體的初始位移為0(忽略通項式中的第三項)，則可以進一步簡化式(25)，得式(26):

3.3 軟件算法數(shù)值阻尼計算

根據(jù)結構動力學理論[9]，位移按e-ξwt衰減，即：在一個周期內，位移將衰減為原來的e-ξw2π/w=e-2πξ。一個周期內，在復平面上，復數(shù)需要旋轉的次數(shù)為2π/θ，投影長度將變?yōu)樵瓉淼膃-2πξ倍，計算式如式(27):

因此，由式(27)可以推出阻尼比公式：

因此，軟件算法引入數(shù)值阻尼的命題得到證明。

3.4 數(shù)值阻尼與采樣步長、自振周期的關系

由于軟件算法將引入數(shù)值阻尼，數(shù)值阻尼與矩陣A的特征值有關，同時矩陣A中的特征值與矩陣A中行列的數(shù)值有關，矩陣A中行列的數(shù)值與wΔt有關，因此軟件算法引入的數(shù)值阻尼僅僅與wΔt相關。

由于wΔt在分析中不直觀，因此筆者對wΔt進行如下變換：wΔt=(2π/T)Δt=2π(Δt/T)。Δt為采樣時間間隔，T為結構自振周期，該比值直觀地反映了單位周期內采樣密度，在之后的分析中將采用Δt/T作為主要參數(shù)。

取一系列Δt/T值，計算對應的wΔt值，將wΔt值代入矩陣A，并計算矩陣A的特征值，根據(jù)復數(shù)特征值實部和虛部大小，結合式(28)計算阻尼比。計算結果如圖5。

由圖5可知，當Δt/T小于0.01時，數(shù)值阻尼呈現(xiàn)增長趨勢，數(shù)值計算的結果是穩(wěn)定的；當Δt/T大于0.3時，數(shù)值阻尼為負，數(shù)值計算的結果是不穩(wěn)定的。表3給出Δt/T介于0.0001～0.01時數(shù)值阻尼的大小，從表3中可以看出數(shù)值阻尼隨Δt/T的變化是接近線性的，同時若結構初始的阻尼值較小，則數(shù)值阻尼將引起較大的誤差。

表3 Δt/T介于0.0001～0.01時數(shù)值阻尼的大小Table 3 Numerical damping while the value of Δt/T is between 0.0001～0.01

4 算例驗證

4.1 算例的數(shù)值阻尼計算

取前文的數(shù)值算例對理論進行驗證。易知：Δt=0.001 s，T=0.1 s，Δt/T=0.01。根據(jù)表3可知，數(shù)值阻尼為0.0157，與圖2結合阻尼比與位移關系式計算所得的數(shù)值阻尼0.0157一致，說明理論公式的正確性以及精確性。

4.2 利用數(shù)值阻尼修正算法

軟件算法的數(shù)值阻尼為0.0157，初始阻尼0.015，因此在實際計算中將初始阻尼扣除數(shù)值阻尼，即為-0.0007，將扣除數(shù)值阻尼的軟件算法記為修正算法。

圖6為阻尼比為0.015時，Newmark算法和修正算法的位移時程圖。采用Newmark算法時，穩(wěn)定的位移最大值為2.07×10-4，與理論值的誤差在0.5%以內；采用修正算法時，穩(wěn)定的位移最大值為2.02×10-4，與理論值的誤差在2.0%。

修正算法計算精度可以滿足要求。數(shù)值阻尼的理論公式，可以修正軟件算法存在數(shù)值阻尼的缺陷。

4.3 Δt/T取值范圍拓展討論

文章公式主要應用于輸電塔鋼管構件渦激振動，因此采樣頻率范圍取為每周期20～2000次，可以涵蓋實際工程應用，對應Δt/T介于0.0005～0.05。因此，本節(jié)中Δt/T討論的范圍為0.0005～0.05。

當Δt/T=0.01時，算例的數(shù)值阻尼與理論公式計算所得數(shù)值阻尼一致。下面將給出Δt/T取值介于0.0005～0.05時，算例的數(shù)值阻尼與理論公式計算的數(shù)值阻尼的比較，以進一步驗證理論公式的正確性。

在表4中，算例數(shù)值阻尼的獲取方式和1.3節(jié)無阻尼時Newmark算法與軟件算法比較中所用方式一致，均采用假定外荷載為0，初始位移為0.0005 m，獲取位移衰減時程圖，根據(jù)公式(3)計算數(shù)值阻尼。理論數(shù)值阻尼根據(jù)表3獲得。

表4 算例數(shù)值阻尼與理論數(shù)值阻尼比較Table 4 Numerical damping of cases and theoretical formula

由表4可知，當Δt/T小于0.01時，理論公式的計算誤差在1%以內；當Δt/T為0.03以及0.05時，理論公式計算精度下降。計算精度下降的原因如下：采樣間隔過大，算法自身計算精度下降，計算精度導致的誤差比例上升，導致理論數(shù)值阻尼與算例數(shù)值阻尼誤差加大。

因此，實際應用中，采用理論公式求解數(shù)值阻尼并用數(shù)值阻尼修正軟件算法的必要條件是Δt/T小于0.01。

5 工程應用

FLUENT動網(wǎng)格宏模塊DEFINE_CG_MOTION在實際工程中被用于模擬空氣對彈性體的作用，文獻[2-3,14-15]均涉及動網(wǎng)格宏模塊在空氣動力學方面的應用。采用本文的數(shù)值阻尼理論公式修正軟件算法的數(shù)值阻尼，可以提高軟件的計算精度，對相關的研究工作及工程應用有實用價值。

本文的算例是輸電塔鋼管構件的渦振，鋼管構件渦振由一階振型控制，因此力學模型采用單自由度體系，之后的理論均是建立在單自由度體系的理論基礎上。鑒于以上原因，多自由度體系不在本文考慮之內?？紤]其它多自由度結構體系時模型時，因多自由度體系的復雜性，數(shù)值阻尼的求解較為困難，筆者建議可以從如下三種方法著手，嘗試處理數(shù)值阻尼造成的誤差：

1) 計算機運行能力足夠的條件下，可以采用減小Δt/T比值的方法，減小數(shù)值阻尼的影響。

2) 當物體的位移為主要參考變量時，可以返回修正速度，使得修正速度與時間步長計算得的位移為真值，確保位移的正確性。

3) 采用其他的數(shù)值積分方法，使得數(shù)值積分方法與軟件的數(shù)據(jù)傳遞模式相匹配。

6 結 論

1) 利用FLUENT進行空氣動力分析，并采用常加速度Newmark法計算物體運動參數(shù)時，常用的計算模塊是動網(wǎng)格宏模塊。動網(wǎng)格宏模塊由于自身數(shù)據(jù)傳遞方式的限制，修改了常加速度Newmark法的原有算法，引入了數(shù)值阻尼。

2) 數(shù)學方法證明了數(shù)值阻尼的存在，并得出數(shù)值阻尼的理論計算公式。軟件算法的數(shù)值阻尼僅與Δt/T的值有關(Δt為采樣步長，T為自振周期)，當Δt/T小于0.01時，數(shù)值阻尼呈現(xiàn)增長趨勢，數(shù)值計算的結果是穩(wěn)定的；當Δt/T大于0.3時，數(shù)值阻尼為負，數(shù)值計算的結果是不穩(wěn)定的。

3) 數(shù)值算例的驗證結果表明：利用數(shù)值阻尼理論公式可以修正軟件算法存在數(shù)值阻尼的缺陷。

4) 實際應用中，采用理論公式求解數(shù)值阻尼并用數(shù)值阻尼修正軟件算法的必要條件是Δt/T小于0.01。

[1]ANSYS Inc.ANSYS FLUENT 14.0 UDF Manual[M].PA,USA,2011:181-182.

[2]李田,張繼業(yè),張衛(wèi)華,等.二維彈性圓柱渦致振動的尾渦模態(tài)[J].空氣動力學學報,2010,28(6):689-695.

[3]徐楓,歐進萍.正三角形排列三圓柱繞流與渦致振動數(shù)值模擬[J].空氣動力學學報,2010,28(5):582-590.

[4]艾尚茂,孫麗萍,沙勇,等.時間步長對低質量比圓柱渦激振動數(shù)值結果的影響[J].船海工程,2011,40(5):168-171.

[5]翟婉明.車輛-軌道耦合動力學[M].北京:中國鐵道出版社,2001:398.

[6]邢譽峰 郭靜.與結構動特性協(xié)同的自適應Newmark方法[J].力學學報,2012,44(5):904-911.

[7]李峰.輸電塔典型節(jié)點鋼管桿件渦激振動研究[D].同濟大學,2008:57.

[8]韓占忠.FLUENT:流體工程仿真計算實例與分析[M].北京:北京理工大學出版社,2009:160.

[9]彭津J,克拉夫R.結構動力學第二版 (修訂版)[M].盛宏玉,譯.北京:高等教育出版社.2006:27-30.

[10]方秦,陳志龍.顯式 Newmark法求解波動問題精度的探討[J].巖土工程學報,1993,15(1):10-15.

[11]李鴻晶,王通,廖旭.關于Newmark-β法機理的一種解釋[J].地震工程與工程振動,2011,31(2):55-62.

[12]王勖成.有限單元法[M].北京:清華大學出版社,2003:491-495.

[13]居于馬.線性代數(shù)[M].北京:清華大學出版社,2002:223-230.

[14]車競,唐碩,謝長強.空射導彈發(fā)射初始彈道數(shù)值仿真[J].空氣動力學學報,2006,23(2):205-208.

[15]朱雄峰,郭正,侯中喜.基于動網(wǎng)格高空長航時機翼優(yōu)化[J].空氣動力學學報,2014,32(4):468-474.

Numericaldampingofdynamicmeshforaerodynamicanalysis

ZHAO Zhangfeng,DENG Hongzhou*

(DepartmentofStructuralEngineering,TongjiUniversity,Shanghai200092,China)

In order to carry out aerodynamic analysis with the FLUENT software and constant averaged acceleration Newmark method on dynamic mesh,the algorithm in the Newmark method needs to be changed due to the limitation of the data transfer mode.For the validation of this modification,numerical tests were simulated to show the deficiency in the algorithm.Based on these tests,we propose a hypothesis that numerical damping is introduced by the algorithm.This assumption was proved by mathematics method,and theoretical formula was provided to calculate the numerical damping.The effectiveness of this formula was validated by calculating numerical dumping in a test and calibrating the algorithm with the calculated numerical dumping.Moreover,this formula was applied to some engineering applications.

FLUENT; dynamic mesh; Newmark method; numerical damping

0258-1825(2017)06-0860-06

V211.3

A

10.7638/kqdlxxb-2015.0080

2015-06-25;

2015-08-23

國家自然科學基金(51578421)

趙張峰(1990-)，男，浙江嵊州人，碩士研究生，研究方向：輸電塔鋼管構件渦激振動.E-mail:zzfmmf@163.com

鄧洪州*(1960-)，男，教授，研究方向：輸電塔優(yōu)化設計等.E-mail:denghz@#edu.cn

趙張峰,鄧洪洲.空氣動力分析中動網(wǎng)格技術的數(shù)值阻尼[J].空氣動力學學報,2017,35(6):860-865.

10.7638/kqdlxxb-2015.0080 ZHAO Z F,DENG H Z.Numerical damping of dynamic mesh for aerodynamic analysis[J].Acta Aerodynamica Sinica,2017,35(6):860-865.