基于虛擬彈性體的快速動網(wǎng)格方法

仲繼澤,徐自力,陶磊

(1.西安交通大學機械結構強度與振動國家重點實驗室,710049,西安;2.西安交通大學航天航空學院,710049, 西安; 3.國家電網(wǎng)浙江北侖第一發(fā)電有限責任公司,315800,浙江寧波)

?

基于虛擬彈性體的快速動網(wǎng)格方法

仲繼澤1,2,徐自力1,2,陶磊3

(1.西安交通大學機械結構強度與振動國家重點實驗室,710049,西安;2.西安交通大學航天航空學院,710049, 西安; 3.國家電網(wǎng)浙江北侖第一發(fā)電有限責任公司,315800,浙江寧波)

為減少流固耦合的計算時間,在現(xiàn)有彈性體方法的基礎上,發(fā)展了一種快速動網(wǎng)格方法。首先根據(jù)彈性體方法的基本假設,將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體,然后將結構與該虛擬彈性體視為一個整體系統(tǒng),并計算其固有振動的振型及頻率,最后以結構受到的流體作用力為激勵,通過振型疊加法計算結構網(wǎng)格及流場網(wǎng)格節(jié)點的位移。考慮到實際結構的流固耦合振動多為低階模態(tài)的振動,在流固耦合計算中可以通過低階模態(tài)的疊加計算流場網(wǎng)格節(jié)點的位移,從而達到快速更新流場網(wǎng)格的目的。采用該快速動網(wǎng)格算法,對某彈性梁顫振問題進行了流固耦合分析,計算結果與已有文獻的結果吻合很好,說明了該算法的正確性。與現(xiàn)有的彈性體方法相比,該算法使流固耦合計算時間減少了65.5%。對Wing 445.6模型的顫振問題進行了分析,得到顫振邊界與實驗值吻合良好,且與現(xiàn)有彈性體法相比,可以減少計算時間54.8%。

流固耦合;彈性體方法;快速動網(wǎng)格方法

在流固耦合計算中,為考慮結構振動對流場的影響,通常需要采用動網(wǎng)格方法更新流場網(wǎng)格[1]。目前所發(fā)展的網(wǎng)格變形方法主要有彈簧法[2]、彈性體方法[3]、溫度體方法[4]及徑向基函數(shù)方法[5]。與其他動網(wǎng)格方法相比,彈性體方法是一種基于物理模型的網(wǎng)格變形方法,把流場網(wǎng)格變形作為一個彈性體變形問題進行求解,更能貼近結構彈性變形的物理實際。文獻[6]率先提出了彈性體方法,將流場網(wǎng)格所圍繞的空間區(qū)域視為虛擬彈性體,流場邊界變形作為彈性體的位移載荷,通過求解靜力平衡方程計算網(wǎng)格節(jié)點的位移。

在早期的研究中,文獻[6-8]都假定虛擬彈性體的彈性模量是各向同性的,采用這種彈性體方法更新流場網(wǎng)格時,邊界層網(wǎng)格容易發(fā)生畸變,網(wǎng)格變形能力受到限制。文獻[9]引入了剛度系數(shù),將虛擬彈性體剛度與流場網(wǎng)格大小關聯(lián),通過增加小網(wǎng)格附近虛擬彈性體的剛度,使得網(wǎng)格變形能力和變形網(wǎng)格的質量得到顯著提高。文獻[10]提出了分層彈性體方法,采用分層的思想,逐層決定虛擬彈性體的剛度,提高了網(wǎng)格變形能力。但是,在流固耦合計算中,每一時間步,都需要更新流場網(wǎng)格,且流場網(wǎng)格的數(shù)量通常為100萬甚至是1 000萬,因此流場網(wǎng)格變形的靜力平衡方程規(guī)模大,反復求解會占用大量的計算時間。

本文在彈性體方法的基礎上,發(fā)展了一種快速動網(wǎng)格方法。將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體,推導了結構與虛擬彈性體同步振動的動力學方程。在流固耦合計算的每一時間步,通過低階模態(tài)的疊加,計算流場網(wǎng)格變形的節(jié)點位移,并更新流場網(wǎng)格。采用本文的動網(wǎng)格算法,對彈性梁顫振算例進行了流固耦合計算,計算結果與已有文獻的結果吻合很好,說明了本文算法的正確性。

1 快速動網(wǎng)格方法

結構的網(wǎng)格節(jié)點可以分成流固耦合面節(jié)點和非流固耦合面節(jié)點兩個部分,據(jù)此將結構的振動控制方程寫成分塊矩陣的形式

(1)

式中:Ms、Cs、Ks分別是結構的質量矩陣、阻尼矩陣和剛度矩陣;Xs_in、Xs_n_in分別是結構流固耦合面節(jié)點位移和非流固耦合面節(jié)點位移;Fin是結構流固耦合面受到的流體作用力。

結構振動會改變流場的邊界,需要采用動網(wǎng)格方法對流場網(wǎng)格節(jié)點坐標進行更新。在彈性體方法中,流場網(wǎng)格所圍繞的空間區(qū)域被視為虛擬彈性體,首先設定虛擬彈性體的彈性模量Ep,然后采用有限元方法計算得到的整體剛度矩陣Kp。在虛擬彈性體變形時,流場網(wǎng)格節(jié)點位移Xp滿足靜力平衡方程[10]

(2)

式中:Xp_in、Xp_n_in分別是流場網(wǎng)格流固耦合面節(jié)點位移和非流固耦合面節(jié)點位移。對結構的網(wǎng)格流固耦合面節(jié)點的位移Xs_in進行插值,可以計算出流場網(wǎng)格流固耦合面節(jié)點的位移

Xp_in=NXs_in

(3)

式中:N是插值系數(shù)矩陣。流場網(wǎng)格變形的靜力平衡方程可以改寫為

(4)

將結構與流場網(wǎng)格域的虛擬彈性體視為一個整體系統(tǒng),通過疊加結構振動控制方程和虛擬彈性體的靜力平衡方程,得到該整體系統(tǒng)的振動方程

(5)

將該動力學系統(tǒng)的固有頻率和正則振型分別記為ω和Ψ。式(5)經(jīng)坐標變換后,得到正則坐標ξ下的動力學方程

(6)

式中:ζ是阻尼比;Q是正則坐標下的流體作用力。

(7)

(8)

(9)

將式(6)寫成分量的形式如下

(10)

模態(tài)位移可以采用Wilson-θ方法[11]計算得到。實際結構的流固耦合振動多與低階模態(tài)相關[12],因此通過低階模態(tài)的疊加可以快速計算出結構及流場網(wǎng)格節(jié)點的位移

(11)

式中:n是參與流固耦合振動的模態(tài)的數(shù)量。

將固有頻率矩陣寫成分塊矩陣的形式

(12)

式中:ωs是結構振動主導的模態(tài)的固有頻率;ωp是虛擬彈性體振動主導的模態(tài)的固有頻率。

固有模態(tài)是指振動系統(tǒng)自由振動時的固有頻率及振型。本文將結構與流場網(wǎng)格域的虛擬彈性體作為一個整體系統(tǒng)。在整體系統(tǒng)的固有模態(tài)中,如果結構部分基本保持了自身的固有模態(tài),將整體系統(tǒng)的固有模態(tài)稱為結構振動主導的模態(tài)。將振型矩陣寫成分塊矩陣的形式

(13)

式中:Ψss是整體系統(tǒng)振型對應的結構部分的振型;Ψpp是整體系統(tǒng)振型對應的虛擬彈性體部分的振型。

根據(jù)固有頻率與剛度矩陣之間的關系式(8),結構振動主導的模態(tài)的固有頻率可以采用下式表達

(14)

將結構部分的振動Ψss代入到結構的振動控制方程式(1)中,可以得到結構固有頻率的表達式為

(15)

與結構的剛度相比,虛擬彈性體的剛度必須是小量。即,虛擬彈性體的彈性模量必須滿足下式

(16)

式中:Es為結構的彈性模量。

2 算例驗證

采用文獻[13]的彈性梁顫振算例驗證本文提出的動網(wǎng)格高效算法。流場及彈性梁的具體尺寸見圖1。梁的密度為1.0×104kg/m3,彈性模量為1.4×106Pa,泊松比為0.4。流體的密度為1 000 kg/m3,動力學黏度為1 Pa·s。流場的左端面為流場的入口,采用速度進口邊界條件,流速方向與進口端面垂直,大小沿y方向呈拋物線分布

(17)

流場的右端面為流場的出口,采用壓力出口邊界,設定出口的平均靜壓值為0。流場的上端面和下端面、彈性梁的流固耦合面及左端的圓形壁面均采用無滑移邊界。

圖1 彈性梁及附近流場區(qū)域

流場及彈性梁的網(wǎng)格見圖2。梁的網(wǎng)格共有4節(jié)點矩形單元153個,節(jié)點208個;流場共有4節(jié)點矩形單元19 627個,節(jié)點20 034個。

圖2 彈性梁及流場網(wǎng)格

2.1 流場網(wǎng)格域虛擬彈性體的彈性模量

本文將虛擬彈性體和結構作為一個整體系統(tǒng)進行研究,結構的固有振動特性勢必受到虛擬彈性體的影響。為了保證計算結果的正確性,必須使結構基本上保持原有的固有振動特性。如果設置的虛擬彈性體彈性模量不恰當,會使彈性梁振動主導的模態(tài)頻率及振型偏離其固有振動頻率和振型,從而產(chǎn)生錯誤的計算結果,因此必須對流場網(wǎng)格設定合適的彈性模量。本文以彈性模量比對數(shù)為參數(shù),分析虛擬彈性體對彈性梁振動主導的模態(tài)頻率及振型的影響,以確定虛擬彈性體的彈性模量。彈性模量比對數(shù)定義如下

R=lg(Es/Ep)

(18)

頻率偏差定義如下

(19)

式中:fs-p是彈性梁振動主導的模態(tài)的固有頻率;fs是彈性梁的固有頻率。

頻率偏差與彈性模量比對數(shù)的關系見圖3。從圖中可以看出:在虛擬彈性體的彈性模量一定的情況下,隨著模態(tài)階數(shù)的增加,頻率偏差會下降;隨著彈性模量比對數(shù)的增加,頻率偏差逐漸減小;當彈性模量比對數(shù)大于等于6,虛擬彈性體對彈性梁振動主導的模態(tài)頻率的影響基本上可以忽略。

圖3 頻率偏差與彈性模量比對數(shù)的關系

彈性模量比對數(shù)與彈性梁振動主導的模態(tài)振型之間的關系見圖4,其中位移和彈性梁長度均為無量綱量。隨著彈性模量比對數(shù)的增大,虛擬彈性體對彈性梁振型的影響逐漸減小。當彈性模量比對數(shù)大于等于5,虛擬彈性體對彈性梁振動主導的模態(tài)振型的影響基本上可以忽略。

(a)第1階振型

(b)第2階振型

(c)第3階振型

本文的研究結果表明,當虛擬彈性體的彈性模量小于等于結構彈性模量的1.0×106倍時,虛擬彈性體對彈性梁固有頻率及振型的影響均可忽略。如果虛擬彈性體彈性模量取值過小,會超出計算機的精度范圍,造成浮點溢出。因此,本文建議虛擬彈性體彈性模量的取值如下

(d)第4階振型圖4 模態(tài)振型與彈性模量比對數(shù)的關系

(20)

此時,彈性梁及虛擬彈性體整體系統(tǒng)的各階固有模態(tài)的振型見圖5。

(a)第1階模態(tài)

(b)第2階模態(tài)

(c)第3階模態(tài)

(d)第4階模態(tài)圖5 不同模態(tài)下的流場網(wǎng)格變形

2.2 流固耦合響應分析

本文以定常流場的結果作為瞬態(tài)流場的初值,采用該動網(wǎng)格算法對彈性梁顫振問題進行分析。設定時間步長為0.001 s,流場變量的相對殘差隨迭代步變化的曲線見圖6。流場變量的相對殘差減小,在進行25步迭代之后,殘差基本保持不變。

圖6 流場變量相對殘差與迭代次數(shù)的關系

(a)第1階

(b)第2階

(c)第3階

(d)第4階圖7 彈性梁的模態(tài)位移時間曲線

本文所采用的流動條件正是文獻[13]中彈性梁顫振的流動條件。在此流動條件下,本文計算得到的彈性梁振動的模態(tài)位移與時間曲線見圖7,模態(tài)位移為質量歸一化量。隨著時間的推進,第1階,第3階和第4階振動的模態(tài)位移幅值逐漸減小,即不會發(fā)生顫振。第2階模態(tài)振動的位移幅值不隨時間變化,即第2階模態(tài)的振動處于顫振臨界點。該彈性梁的顫振為第2階彎曲顫振。

(d)T圖8 振動彈性梁附近的流場

采用流固耦合的方法,計算得到彈性梁右端中點的位移和時間曲線見圖9,從中可以看出,本文的計算結果與文獻[13]的結果吻合很好,驗證了本文動網(wǎng)格算法的正確性。本文計算得到圖9所示的彈性梁顫振的位移時間曲線,總耗時為19.8

h

。在相同計算條件下,采用傳統(tǒng)彈性體網(wǎng)格變形方法

[10]

的計算時間為57.4

h

,而本文方法可以減少計算時間65.5%。

為進一步驗證本文算法,我們對Wing445.6模型[14]的顫振問題進行了流固耦合分析,計算得到的顫振邊界與實驗值吻合良好,計算結果見圖10,與已有彈性體方法相比,計算時間減少了54.8%。

(a)x方向位移

(b)y方向位移圖9 彈性梁右端中點的位移和時間的關系

圖10 Wing 445.6模型的顫振邊界

3 結 論

本文提出了一種基于虛擬彈性體的快速動網(wǎng)格方法。將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體,并與結構一起,視為一個整體系統(tǒng),以結構受到的流體作用力為激勵,通過振型疊加法計算結構網(wǎng)格及流場網(wǎng)格節(jié)點的位移。考慮到實際結構的流固耦合振動多為低階模態(tài)的振動,在流固耦合計算中可以通過低階模態(tài)的疊加計算流場網(wǎng)格節(jié)點的位移,從而達到快速更新流場網(wǎng)格的目的。采用該方法對彈性梁顫振問題進行了流固耦合分析,計算結果與文獻的結果吻合很好。與已有彈性體方法相比,本文算法可以減少合計算時間65.5%。研究發(fā)現(xiàn),流場渦旋的周期性形成、發(fā)展及脫落是彈性梁顫振的主要原因。對Wing445.6模型顫振問題進行了分析,得到顫振邊界與實驗值吻合良好,且與已有彈性體法相比,可以減少計算時間54.8%。

[1] 張偉偉, 高傳強, 葉正寅. 氣動彈性計算中網(wǎng)格變形方法研究進展 [J]. 航空學報, 2014, 35(2): 303-319. ZHANG Weiwei, GAO Chuanqiang, YE Zhengyin. Progress on mesh deformation in computational aeroelasticity [J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(2): 303-319.

[2] 史愛明, 楊青, 楊永年. 非結構運動網(wǎng)格下的三維機翼顫振數(shù)值分析 [J]. 振動與沖擊, 2006, 24(6): 27-28. SHI Aiming, YANG Qing, YANG Yongnian. Numerical flutter analysis of a 3-D wing using unstructured dynamic mesh Euler method [J]. Journal of Vibration and Shock, 2006, 24(6): 27-28.

[3] TEZDUYAR T E. Stabilized finite element formulations for incompressible flow computations [J]. Advances in Applied Mechanics, 1991, 28: 1-44.

[4] 陳炎, 曹樹良, 梁開洪, 等. 基于溫度體模型的動網(wǎng)格生成方法及在流固耦合振動中的應用 [J]. 振動與沖擊, 2010, 29(4): 1-5. CHEN Yan, CAO Shuliang, LIANG Kaihong, et al. A new dynamic grids based on temperature analogy and its application in vibration engineering with fluid-solid interaction [J]. Journal of Vibration and Shock, 2010, 29(4): 1-5.

[5] 謝亮, 徐敏, 張斌, 等. 基于徑向基函數(shù)的高效網(wǎng)格變形算法研究 [J]. 振動與沖擊, 2013, 32(10): 141-145. XIE Liang, XU Min, ZHANG Bin, et al. Space points reduction in grid deforming method based on radial basis functions [J]. Journal of Vibration and Shock, 2013, 32(10): 141-145.

[6] TEZDUYAR T E, BEHR M, MITTAL S, et al. A new strategy for finite element computations involving moving boundaries and interfaces-the deforming-spatial-domain/space-time procedure: I The concept and the preliminary tests [J]. Computer Methods in Applied Mechanics and Engineering, 1992, 94(3): 339-351.

[7] TEZDUYAR T E, BEHR M, MITTAL S, et al. A new strategy for finite element computations involving moving boundaries and interfaces-the deforming-spatial-domain/space-time procedure: II Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders [J]. Computer Methods in Applied Mechanics and Engineering, 1992, 94(3): 353-371.

[8] BAR-YOSEPH P Z, MEREU S, CHIPPADA S, et al. Automatic monitoring of element shape quality in 2-D and 3-D computational mesh dynamics [J]. Computational Mechanics, 2001, 27(5): 378-395.

[9] STEIN K, TEZDUYAR T, BENNEY R. Mesh moving techniques for fluid-structure interactions with large displacements [J]. Journal of Applied Mechanics, 2003, 70(1): 58-63.

[10]HUO S H, WANG F S, YAN W Z, et al. Layered elastic solid method for the generation of unstructured dynamic mesh [J]. Finite Elements in Analysis and Design, 2010, 46(10): 949-955.

[11]ZIENKIEWICZ O C, PAUL D K, HINTON E. Cavitation in fluid-structure response with particular reference to dams under earthquake loading [J]. Earthquake Engineering & Structural Dynamics, 1983, 11(4): 463-481.

[12]MARSHALL J G, IMREGUN M. A review of aeroelasticity methods with emphasis on turbomachinery applications [J]. Journal of Fluids and Structures, 1996, 10(3): 237-267.

[13]TUREK S, HRON J. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow [M]. Berlin, Germany: Springer, 2006: 371-385.

[14]YATES E C, Jr. AGARD standard aeroelastic configurations for dynamic response I-wing 445.6 [R]. Neuilly Sur Seine, France: Advisory Group for Aerospace Research and Development Neuilly-Sur-Seine, 1988: 1-74.

(編輯 杜秀杰 趙煒)

An Efficient Dynamic Mesh Method Based on Pseudo Elastic Solid

ZHONG Jize1,2,XU Zili1,2,TAO Lei3

(1. State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University,Xi’an 710049, China; 2. School of Aerospace, Xi’an Jiaotong University, Xi’an 710049, China;3. State Grid Zhejiang Beilun No.1 Power Generation Co. Ltd., Ningbo, Zhejiang 315800, China)

To reduce the computing time of fluid structure coupling, an efficient dynamic mesh method based on the pre-existing elastic solid method is developed. The flow mesh domain is assumed to be a pseudo elastic solid according to the basic hypotheses of the elastic solid method, then the structure and the pseudo elastic solid are considered together as one holistic system. Subsequently the natural frequencies and vibration modes are calculated for the system and the flow force acted on the structure is considered as the excitation of the holistic system. The nodal displacements for the structure and the flow mesh are computed by mode superposition. In fact, the actual fluid structure coupled vibration for structures often appears associated with low order modes, the nodal displacements of the flow mesh can be calculated by modal superposition of the first few of low order modes, thus the flow mesh can be updated efficiently. A beam flutter problem is discussed with the present dynamic mesh method. The results coincide well with the data reported in the existing reference verifying the validation of the present method. The computing time is reduced by 65.5% compared with the pre-existing elastic solid method. The flutter of wing 445.6 is also analyzed and the calculated flutter boundary agrees with the experimental data, and the computing time is reduced by 54.8%.

fluid structure interaction; elastic solid method; efficient dynamic mesh method

2016-03-22。

仲繼澤(1988—),男,博士生;徐自力(通信作者),男,教授,博士生導師。

國家自然科學基金資助項目(51275385)。

時間:2016-07-14

http:∥www.cnki.net/kcms/detail/61.1069.T.20160714.1719.002.html

10.7652/xjtuxb201610020

O323;O351.2

A

0253-987X(2016)10-0132-07