仲繼澤,謝志強(qiáng),沈渡,徐自力

(1.中國艦船研究設(shè)計(jì)中心船舶振動噪聲重點(diǎn)實(shí)驗(yàn)室,430064,武漢;2.西安交通大學(xué)機(jī)械結(jié)構(gòu)強(qiáng)度與振動國家重點(diǎn)實(shí)驗(yàn)室,710049,西安)

基于動網(wǎng)格的流固耦合方法可以借助目前已成熟的CFD和CSD商業(yè)軟件進(jìn)行流場分析和結(jié)構(gòu)振動計(jì)算,從而省去繁重的代碼開發(fā)工作,目前的動網(wǎng)格技術(shù)已能夠處理二維及三維網(wǎng)格變形問題,已經(jīng)廣泛應(yīng)用在工程實(shí)際的具有復(fù)雜氣動外形的結(jié)構(gòu)流固耦合問題的研究中[1]。

在眾多流場動網(wǎng)格方法中,彈性體動網(wǎng)格法符合工程實(shí)際中結(jié)構(gòu)流固耦合振動多為彈性變形的物理實(shí)際,容易產(chǎn)生高質(zhì)量的動網(wǎng)格[2]。彈性體法最早由Tezduyar提出并發(fā)展,將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體,采用單一的彈性模量值,以流場網(wǎng)格作為虛擬彈性體的空間離散網(wǎng)格,然后通過有限元方法構(gòu)建虛擬彈性體變形的靜力平衡方程,通過求解該方程計(jì)算網(wǎng)格變形后網(wǎng)格節(jié)點(diǎn)新的位置坐標(biāo)[3]。但是,工程實(shí)際結(jié)構(gòu)的流場網(wǎng)格的規(guī)模一般都在百萬量級,所以彈性體法構(gòu)建的靜力平衡方程的規(guī)模一般比較龐大,網(wǎng)格變形的計(jì)算會占用大量的計(jì)算資源[4]。本課題組在Tezduyar的研究工作的基礎(chǔ)上,提出一種基于單一彈性模量的快速動網(wǎng)格方法[5-6],構(gòu)建了結(jié)構(gòu)-虛擬彈性體的整體動力學(xué)方程,然后通過振型截?cái)喾椒ㄓ?jì)算結(jié)構(gòu)和流場網(wǎng)格節(jié)點(diǎn)的位移,其計(jì)算量只有Tezduyar的彈性體法的萬分之一。本課題組在后續(xù)的研究中發(fā)現(xiàn),該快速動網(wǎng)格方法對微小變形問題的適用性很好,但是當(dāng)結(jié)構(gòu)變形較大時(shí),會造成流場邊界層網(wǎng)格的扭曲,使得網(wǎng)格質(zhì)量降低。

為了解決基于單一彈性體的快速動網(wǎng)格方法的網(wǎng)格扭曲問題,本文借鑒文獻(xiàn)[7-8]中提高彈性體動網(wǎng)格方法網(wǎng)格變形質(zhì)量的做法,引入空間分布彈性模量以改進(jìn)上述快速動網(wǎng)格方法,提出了一種基于空間分布彈性模量的快速動網(wǎng)格方法。采用改進(jìn)后的快速動網(wǎng)格方法計(jì)算了較大變形條件下的標(biāo)準(zhǔn)彈性梁和Agard Wing 445.6的流場網(wǎng)格變形,并對比了快速動網(wǎng)格方法改進(jìn)前后對流場網(wǎng)格變形質(zhì)量的影響。

1 基于單一彈性模量的快速動網(wǎng)格方法的網(wǎng)格扭曲問題

為了準(zhǔn)確模擬邊界層的流動,流場網(wǎng)格一般在結(jié)構(gòu)表面附近分布比較密集,網(wǎng)格尺寸較小,在遠(yuǎn)離結(jié)構(gòu)的地方網(wǎng)格尺寸較大?；趩我粡椥阅Ａ康目焖賱泳W(wǎng)格方法不考慮流場網(wǎng)格的大小隨空間分布發(fā)生變化的實(shí)際情況,對于整個虛擬彈性體設(shè)定單一的彈性模量,當(dāng)結(jié)構(gòu)變形較小時(shí)網(wǎng)格變形質(zhì)量是可以保證的,但是當(dāng)結(jié)構(gòu)變形較大時(shí)邊界層的尺寸較小的網(wǎng)格可能產(chǎn)生扭曲,使得網(wǎng)格變形質(zhì)量降低。圖1給出了采用基于單一彈性模量的快速動網(wǎng)格方法對標(biāo)準(zhǔn)彈性梁流固耦合問題[9]的流場網(wǎng)格變形情況,可以看出,隨著彈性梁最大位移的增加,變形后的流場網(wǎng)格最小角均呈現(xiàn)逐漸減小的趨勢,而且模態(tài)階數(shù)越高減小的速度越快。在發(fā)生第1階彎曲變形的情況下,彈性梁最大位移為0.07 m時(shí),變形后流場網(wǎng)格最小角為0°。

圖1 變形后網(wǎng)格最小角與彈性梁最大位移的關(guān)系

彈性梁右端的流場動網(wǎng)格變形后會出現(xiàn)非法網(wǎng)格的情況,如圖2所示。

圖2 網(wǎng)格變形后出現(xiàn)非法網(wǎng)格的情況

2 基于空間分布彈性模量的快速動網(wǎng)格方法

為了解決基于單一彈性體的快速動網(wǎng)格方法的網(wǎng)格扭曲問題,本文借鑒文獻(xiàn)[7-8]中提高彈性體動網(wǎng)格方法網(wǎng)格變形質(zhì)量的做法,引入空間分布彈性模量以改進(jìn)基于單一彈性模量的快速動網(wǎng)格方法,從而提出了一種基于空間分布彈性模量的快速動網(wǎng)格方法。本文通過流場網(wǎng)格的大小在空間的分布設(shè)定虛擬彈性體的彈性模量的空間分布

(1)

對于第i個流場網(wǎng)格,其內(nèi)部的虛擬彈性體的應(yīng)力與應(yīng)變之間滿足胡克定律

σ=2Gε+λδε

(2)

式中:σ為虛擬彈性體的應(yīng)力張量;ε為虛擬彈性體的應(yīng)變張量;δ為克羅內(nèi)克爾符號;拉梅常數(shù)G、λ分別為

(3)

(4)

其中ν為虛擬彈性體的泊松比。

將每一個網(wǎng)格單元視為一個有限單元,那么虛擬彈性體的整體的總位能等于所有有限單元位能的總和。虛擬彈性體總位能的有限元形式為

(5)

對于整個虛擬彈性體,其總位能滿足最小位能原理

δΠ(Xg)=0

(6)

化簡后得到

KgXg=Fs

(7)

對于流固耦合問題,通常將結(jié)構(gòu)對虛擬彈性體的作用視為位移條件,所以上述Fs=0,那么基于空間分布彈性模量的虛擬彈性體變形的靜力平衡方程可以寫成如下形式

KgXg=0

(8)

參考基于單一彈性模量的快速動網(wǎng)格方法[5-6],將結(jié)構(gòu)振動方程和虛擬彈性體的靜力平衡方程結(jié)合流固耦合面的位移插值方程,通過數(shù)學(xué)推導(dǎo)得到結(jié)構(gòu)-虛擬彈性體系統(tǒng)的振動方程,然后通過模態(tài)疊加計(jì)算結(jié)構(gòu)和流場網(wǎng)格節(jié)點(diǎn)的位移,以達(dá)到快速計(jì)算流場網(wǎng)格變形的節(jié)點(diǎn)位移的目的。

本文采用空間分布的彈性模量使得流場邊界層的體積較小的網(wǎng)格內(nèi)部的虛擬彈性體的彈性模量較大,使得遠(yuǎn)離邊界層的體積較大的網(wǎng)格內(nèi)部的虛擬彈性體的彈性模量較小,而且虛擬彈性體的彈性模量隨網(wǎng)格體積的增加以指數(shù)規(guī)律減小。當(dāng)流場網(wǎng)格隨結(jié)構(gòu)振動發(fā)生變形時(shí),邊界層網(wǎng)格會隨著結(jié)構(gòu)振動作整體的運(yùn)動,從而避免邊界層的流場網(wǎng)格在變形過程中由于擠壓、拉扯、扭轉(zhuǎn)引起的網(wǎng)格質(zhì)量降低的問題。

對于不會產(chǎn)生非法網(wǎng)格的工況,本文方法與單一彈性模量的快速動網(wǎng)格方法的計(jì)算速度和精度基本一致。對于產(chǎn)生了非法網(wǎng)格的工況,基于單一彈性模量的快速動網(wǎng)格方法產(chǎn)生的非法網(wǎng)格會使計(jì)算無法進(jìn)行,而本文方法可以避免產(chǎn)生非法網(wǎng)格,順利進(jìn)行計(jì)算,同時(shí)也能保證原有的計(jì)算精度和速度。

3 算例驗(yàn)證

3.1 彈性梁的流場動網(wǎng)格

采用本文改進(jìn)后的快速動網(wǎng)格方法對標(biāo)準(zhǔn)彈性梁流固耦合問題[9]的流場網(wǎng)格變形進(jìn)行計(jì)算得到彈性梁流場網(wǎng)格的變形情況,如圖3所示。

(a)網(wǎng)格變形情況

(b)網(wǎng)格變形質(zhì)量對比圖3 彈性梁第1階彎曲變形時(shí)流場網(wǎng)格變形情況

假定彈性梁發(fā)生第1階彎曲振動,當(dāng)最大位移為0.07 m時(shí),其附近流場網(wǎng)格變形如圖3a所示。可以看出,流場網(wǎng)格變形整體比較平順,彈性梁右端未出現(xiàn)非法網(wǎng)格。變形后流場網(wǎng)格最小角的分布情況如圖3b所示,采用本文方法變形后流場網(wǎng)格最小角維持在67°以上,而采用改進(jìn)前的快速動網(wǎng)格方法計(jì)算得到的彈性梁最大位移為0.07 m時(shí)網(wǎng)格變形后的網(wǎng)格最小角為0°,此時(shí)的流場網(wǎng)格已經(jīng)成為非法網(wǎng)格。

3.2 Agard Wing 445.6的流場動網(wǎng)格

采用本文方法對三維的Agard Wing 445.6顫振問題[10]的流場網(wǎng)格變形進(jìn)行計(jì)算得到機(jī)翼流場網(wǎng)格的變形情況,如圖4、圖5所示。

假定機(jī)翼以第1階模態(tài)振動,當(dāng)最大位移為24 cm時(shí),流場網(wǎng)格變形如圖4a所示?？梢钥闯?流場網(wǎng)格變形整體比較平順。變形后流場網(wǎng)格最小角的分布情況如圖4b所示,采用本文方法變形后流場網(wǎng)格最小角維持在36°以上,而采用改進(jìn)前的方法得到的網(wǎng)格變形后的網(wǎng)格最小角為0°,是非法網(wǎng)格。

(a)網(wǎng)格變形情況

(b)網(wǎng)格變形質(zhì)量對比圖4 Agard Wing 445.6第1階彎曲變形時(shí)流場網(wǎng)格變形情況

假定機(jī)翼以第2階模態(tài)振動,當(dāng)最大位移為12 cm時(shí),其附近流場網(wǎng)格變形如圖5a所示,變形后流場網(wǎng)格最小角的分布情況如圖5b所示,采用本文方法變形后所有網(wǎng)格最小角依然維持在36°以上,而采用改進(jìn)前的快速動網(wǎng)格方法計(jì)算得到的機(jī)翼最大位移為12 cm時(shí)網(wǎng)格變形后的網(wǎng)格最小角為0°,此時(shí)的流場網(wǎng)格已經(jīng)成為非法網(wǎng)格。

(a)網(wǎng)格變形情況

(b)網(wǎng)格變形質(zhì)量對比圖5 Agard Wing 445.6第1階扭轉(zhuǎn)變形時(shí)流場網(wǎng)格變形情況

4 結(jié) 論

本文針對基于單一彈性模量的快速動網(wǎng)格方法的網(wǎng)格扭曲問題,參照彈性體方法中提高網(wǎng)格變形質(zhì)量的方法,提出了一種基于空間分布彈性模量的快速動網(wǎng)格方法。采用本文的方法對標(biāo)準(zhǔn)彈性梁流固耦合問題和Agard Wing 445.6顫振問題的流場網(wǎng)格變形進(jìn)行了計(jì)算,發(fā)現(xiàn)當(dāng)單一彈性模量的快速動網(wǎng)格方法的網(wǎng)格變形出現(xiàn)非法網(wǎng)格(網(wǎng)格最小角為0°)時(shí),本文方法得到的彈性梁流場網(wǎng)格變形后的最小角均維持在67°以上,得到的機(jī)翼流場網(wǎng)格變形后的最小角均維持在36°以上,解決了原快速動網(wǎng)格方法的網(wǎng)格扭曲問題。

[1] DOWELL E H. Some recent advances in nonlinear aeroelasticity: fluid-structure interaction in the 21st century [C]∥Proceedings of the AIAA Structures, Structural Dynamics, and Materials Conference. Reston, VA, USA: AIAA, 2010: 1-7.

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

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

[4] ZHOU X, LI S. A new mesh deformation method based on disk relaxation algorithm with pre-displacement and post-smoothing [J]. Journal of Computational Physics, 2013, 235(2): 199-215.

[5] 仲繼澤, 徐自力, 陶磊. 基于虛擬彈性體的快速動網(wǎng)格方法 [J]. 西安交通大學(xué)學(xué)報(bào), 2016, 50(10): 132-138. ZHONG Jize, XU Zili, TAO Lei. An efficient dynamic mesh method based on pseudo elastic solid [J]. Journal of Xi’an Jiaotong University, 2016, 50(10): 132-138.

[6] ZHONG Jize, XU Zili. A modal approach for coupled fluid structure computations of wing flutter [J]. Proceedings of the Institution of Mechanical Engineers: Part G Journal of Aerospace Engineering, 2017, 231(1): 72-81.

[7] 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.

[8] 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.

[9] 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.

[10]YATES E C J. AGARD standard aeroelastic configurations for dynamic response [J]. European Psychiatry, 1987, 25(10): 1596-1597.