基于空間分布彈性模量的快速動網格方法

仲繼澤,謝志強,沈渡,徐自力

(1.中國艦船研究設計中心船舶振動噪聲重點實驗室,430064,武漢;2.西安交通大學機械結構強度與振動國家重點實驗室,710049,西安)

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

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

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

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

為了準確模擬邊界層的流動,流場網格一般在結構表面附近分布比較密集,網格尺寸較小,在遠離結構的地方網格尺寸較大。基于單一彈性模量的快速動網格方法不考慮流場網格的大小隨空間分布發生變化的實際情況,對于整個虛擬彈性體設定單一的彈性模量,當結構變形較小時網格變形質量是可以保證的,但是當結構變形較大時邊界層的尺寸較小的網格可能產生扭曲,使得網格變形質量降低。圖1給出了采用基于單一彈性模量的快速動網格方法對標準彈性梁流固耦合問題[9]的流場網格變形情況,可以看出,隨著彈性梁最大位移的增加,變形后的流場網格最小角均呈現逐漸減小的趨勢,而且模態階數越高減小的速度越快。在發生第1階彎曲變形的情況下,彈性梁最大位移為0.07 m時,變形后流場網格最小角為0°。

圖1 變形后網格最小角與彈性梁最大位移的關系

彈性梁右端的流場動網格變形后會出現非法網格的情況,如圖2所示。

圖2 網格變形后出現非法網格的情況

2 基于空間分布彈性模量的快速動網格方法

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

(1)

對于第i個流場網格,其內部的虛擬彈性體的應力與應變之間滿足胡克定律

σ=2Gε+λδε

(2)

式中:σ為虛擬彈性體的應力張量;ε為虛擬彈性體的應變張量;δ為克羅內克爾符號;拉梅常數G、λ分別為

(3)

(4)

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

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

(5)

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

δΠ(Xg)=0

(6)

化簡后得到

KgXg=Fs

(7)

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

KgXg=0

(8)

參考基于單一彈性模量的快速動網格方法[5-6],將結構振動方程和虛擬彈性體的靜力平衡方程結合流固耦合面的位移插值方程,通過數學推導得到結構-虛擬彈性體系統的振動方程,然后通過模態疊加計算結構和流場網格節點的位移,以達到快速計算流場網格變形的節點位移的目的。

本文采用空間分布的彈性模量使得流場邊界層的體積較小的網格內部的虛擬彈性體的彈性模量較大,使得遠離邊界層的體積較大的網格內部的虛擬彈性體的彈性模量較小,而且虛擬彈性體的彈性模量隨網格體積的增加以指數規律減小。當流場網格隨結構振動發生變形時,邊界層網格會隨著結構振動作整體的運動,從而避免邊界層的流場網格在變形過程中由于擠壓、拉扯、扭轉引起的網格質量降低的問題。

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

3 算例驗證

3.1 彈性梁的流場動網格

采用本文改進后的快速動網格方法對標準彈性梁流固耦合問題[9]的流場網格變形進行計算得到彈性梁流場網格的變形情況,如圖3所示。

(a)網格變形情況

(b)網格變形質量對比圖3 彈性梁第1階彎曲變形時流場網格變形情況

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

3.2 Agard Wing 445.6的流場動網格

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

假定機翼以第1階模態振動,當最大位移為24 cm時,流場網格變形如圖4a所示。可以看出,流場網格變形整體比較平順。變形后流場網格最小角的分布情況如圖4b所示,采用本文方法變形后流場網格最小角維持在36°以上,而采用改進前的方法得到的網格變形后的網格最小角為0°,是非法網格。

(a)網格變形情況

(b)網格變形質量對比圖4 Agard Wing 445.6第1階彎曲變形時流場網格變形情況

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

(a)網格變形情況

(b)網格變形質量對比圖5 Agard Wing 445.6第1階扭轉變形時流場網格變形情況

4 結 論

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

[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] 張偉偉, 高傳強, 葉正寅. 氣動彈性計算中網格變形方法研究進展 [J]. 航空學報, 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] 仲繼澤, 徐自力, 陶磊. 基于虛擬彈性體的快速動網格方法 [J]. 西安交通大學學報, 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.