基于等效板的含離散源損傷機(jī)翼結(jié)構(gòu)分析*

（中國民航大學(xué)天津市民用航空器適航與維修重點(diǎn)實(shí)驗(yàn)室，天津 300300)

機(jī)翼和機(jī)身結(jié)構(gòu)的安全性直接關(guān)系到飛機(jī)的安全，因此相關(guān)的飛機(jī)設(shè)計(jì)規(guī)范或適航條例對(duì)機(jī)翼機(jī)身結(jié)構(gòu)的安全性均有明確的規(guī)定，其中一個(gè)重要的內(nèi)容就是要求機(jī)翼和機(jī)身結(jié)構(gòu)在帶有離散源損傷的情況下, 保證安全飛行[1]。因機(jī)翼和機(jī)身結(jié)構(gòu)復(fù)雜，設(shè)計(jì)難度大，而等效方法提出一種簡單的方案或設(shè)想，使某些特性和規(guī)律等方面的性能和指標(biāo)相同，即真實(shí)結(jié)構(gòu)與等效模型之間的某些特性可以彼此互換，而分析結(jié)論一致，從而通過簡化使問題化繁為簡，快速得到解決辦法。因此，簡化部分機(jī)翼結(jié)構(gòu)，使用等效設(shè)計(jì)方法建立精確的機(jī)翼結(jié)構(gòu)模型，對(duì)于快速有效的分析機(jī)翼特性具有重要作用。

本文把壓縮載荷作用下的含離散源損傷加筋板模型作為機(jī)翼結(jié)構(gòu)[2]，采用Nastran優(yōu)化設(shè)計(jì)方法，使優(yōu)化設(shè)計(jì)之后的含離散源損傷不加筋板模型的靜、動(dòng)力學(xué)特性和含離散源損傷加筋板模型的靜、動(dòng)力學(xué)特性達(dá)到等效，從而使優(yōu)化設(shè)計(jì)之后的不加筋板模型與機(jī)翼結(jié)構(gòu)是等效結(jié)構(gòu)，進(jìn)而在等效模型上展開進(jìn)一步分析和研究。這對(duì)于快速、有效的分析含離散源損傷機(jī)翼結(jié)構(gòu)的力學(xué)特性具有重要作用。

1 離散源損傷初始模型建立

為研究含離散源損傷機(jī)翼結(jié)構(gòu)的靜、動(dòng)力學(xué)特性，建立含離散源損傷不加筋板模型并作為初始模型；建立含離散源損傷加筋板模型即機(jī)翼結(jié)構(gòu)模型并作為目標(biāo)模型。使優(yōu)化設(shè)計(jì)之后的初始模型和目標(biāo)模型的靜、動(dòng)力學(xué)特性達(dá)到等效。含離散源損傷加筋板模型試件的幾何尺寸如圖1所示。

圖1 含離散源損傷加筋板幾何參數(shù)Fig. 1 Geometric parameters of stiffened plate with discrete source damage

模型均采用薄鋁板材料，具體材料參數(shù)如表1所示。

表1 材料參數(shù)

目標(biāo)模型和初始模型的離散源損傷位置和幾何尺寸相同，平面厚度相同且都不穿過兩側(cè)加強(qiáng)筋位置，具體如圖2所示。

兩個(gè)板模型的網(wǎng)格劃分方向和大小相同，為得到更精確的計(jì)算結(jié)果，在離散源損傷周圍對(duì)網(wǎng)格采取加密精細(xì)化處理。因?yàn)楸疚膬?yōu)化設(shè)計(jì)分析過程僅在線性范圍內(nèi)展開，因此對(duì)初始模型進(jìn)行預(yù)加載仿真分析，確保優(yōu)化設(shè)計(jì)過程不出現(xiàn)非線性變化等因素的影響。

圖2 含離散源損傷不加筋板模型Fig.2 Unstiffened plate model with discrete sources damage

通過PATRAN有限元分析可知，其自由端Z向變形位移和最大應(yīng)力水平如圖3、圖4所示：

在對(duì)初始模型的仿真分析中，由分析結(jié)果可知：節(jié)點(diǎn)最大變形位移的百分比函數(shù)是7.99%，安全系數(shù)0.0615＜ns，均符合相關(guān)的國家技術(shù)標(biāo)準(zhǔn)[3]要求。因此，初始模型符合小變形要求和強(qiáng)度理論要求，后續(xù)的優(yōu)化設(shè)計(jì)過程不會(huì)出現(xiàn)非線性等問題。

圖3 最大變形位移圖Fig.3 The biggest deformation displacement diagram

圖4 最大應(yīng)力圖Fig.4 The biggest stress diagram

2 Nastran優(yōu)化設(shè)計(jì)原理

Nastran優(yōu)化設(shè)計(jì)程序基本原理[4]如下：

尋找一組合適的設(shè)計(jì)變量：

{x}={x1,x2,...xn}，

使目標(biāo)函數(shù)的值達(dá)到最?。?/p>

minimize=F(x)，

滿足約束條件：

區(qū)間約束條件：XLi≤Xi≤XUii=1,2...n,n，

不等式約束條件：Gj(x)≤0j=1,2,...,L，

等式約束條件：

Hk(x)=0k=1,2,...,k，

Nastran求解器對(duì)于優(yōu)化設(shè)計(jì),設(shè)置了3種優(yōu)化算法。其中第3種方法是3種算法的核心，也是默認(rèn)的優(yōu)化算法[5]。

其中為第k+1次迭代的搜索方向，是滿足條件θj和可行條件為標(biāo)量的步長值。Nastran軟件對(duì)可行性方向閥的主要改進(jìn)在于：在約束邊界條件上，搜索方向的選擇是下面這個(gè)子優(yōu)化分析問題的解：

通過以上約束方程得到的搜索方向s緊貼著臨界約束的邊界移動(dòng)，使目標(biāo)函數(shù)逐漸下降。自動(dòng)迭代運(yùn)算直至出現(xiàn)最優(yōu)結(jié)果，則停止計(jì)算。

3 靜力學(xué)等效分析

基于Nastran優(yōu)化設(shè)計(jì)原理展開對(duì)初始模型的優(yōu)化設(shè)計(jì)分析，在靜力學(xué)特性的基礎(chǔ)上，通過優(yōu)化設(shè)計(jì)計(jì)算，建立與目標(biāo)模型靜力學(xué)特性一致的靜力學(xué)等效模型，若達(dá)到等效，則稱優(yōu)化設(shè)計(jì)之后的不加筋板模型為靜力學(xué)等效模型。

設(shè)計(jì)變量：設(shè)計(jì)變量取板的厚度，設(shè)置如圖5所示。通過優(yōu)化T1、T2、T3三部分板厚，使初始模型與目標(biāo)模型的靜力學(xué)特性相同。初始厚度設(shè)置均為0.006m。

圖5 優(yōu)化設(shè)計(jì)變量Fig.5 Optimization of the design variables

優(yōu)化設(shè)計(jì)的目標(biāo)函數(shù)：

式中，N是板模型的節(jié)點(diǎn)數(shù)量；wiE是優(yōu)化設(shè)計(jì)之后板模型第i節(jié)點(diǎn)位移；wiF是與優(yōu)化設(shè)計(jì)之后板模型第i節(jié)點(diǎn)相對(duì)應(yīng)的目標(biāo)模型第i節(jié)點(diǎn)的變形位移[6]。

將包含初始模型幾何尺寸、設(shè)計(jì)變量、目標(biāo)函數(shù)等參數(shù)組成的優(yōu)化設(shè)計(jì)程序提交Nastran求解器SOL200迭代計(jì)算[7]，設(shè)計(jì)變量T的變化如表2所示。

如表2中所示，經(jīng)過5次迭代運(yùn)算，找到了使目標(biāo)函數(shù)Φ趨于最小的T值。把厚度最優(yōu)的設(shè)計(jì)變量數(shù)值分別賦予T1、T2、T3，在Patran中重構(gòu)初始模型并進(jìn)行計(jì)算分析。驗(yàn)證優(yōu)化設(shè)計(jì)后的不加筋板模型與目標(biāo)模型之間各對(duì)應(yīng)節(jié)點(diǎn)變形位移和模態(tài)頻率是否相同。

本文選取圖5的靜力優(yōu)化后不加筋板模型和目標(biāo)模型互相對(duì)應(yīng)的15個(gè)節(jié)點(diǎn)的變形位移作對(duì)比分析，此15個(gè)節(jié)點(diǎn)均為同一平面直線上最大的位移變形節(jié)點(diǎn)。驗(yàn)證對(duì)比如表3所示。頻率結(jié)果對(duì)比如表4所示。

表3中，分別賦予優(yōu)化迭代計(jì)算后的3個(gè)板厚值T之后，靜力學(xué)優(yōu)化設(shè)計(jì)后的不加筋板模型各節(jié)點(diǎn)變形位移與目標(biāo)模型對(duì)應(yīng)各節(jié)點(diǎn)變形位移數(shù)值基本相同，靜力學(xué)特性基本達(dá)到等效，可稱此優(yōu)化設(shè)計(jì)后的模型為靜力學(xué)等效模型，即建立了靜力學(xué)等效。表4中的動(dòng)力學(xué)特性模態(tài)頻率方面，靜力學(xué)等效模型和目標(biāo)模型之間的模態(tài)頻率仍相差較大，靜力學(xué)等效模型在動(dòng)力學(xué)方面的等效性則較差。因此，需要在靜力學(xué)等效的基礎(chǔ)上，繼續(xù)展開對(duì)靜力學(xué)等效模型的模態(tài)頻率特性作進(jìn)一步優(yōu)化分析。

表2 靜力學(xué)優(yōu)化設(shè)計(jì)變量變化m

表3 靜力學(xué)優(yōu)化設(shè)計(jì)變形位移對(duì)比

表4 靜力學(xué)優(yōu)化設(shè)計(jì)模態(tài)頻率對(duì)比Hz

4 動(dòng)力學(xué)等效分析

在靜力學(xué)等效模型和目標(biāo)模型的各節(jié)點(diǎn)變形位移相同的基礎(chǔ)上，同時(shí)對(duì)模態(tài)頻率展開優(yōu)化設(shè)計(jì)分析，從而使動(dòng)力學(xué)優(yōu)化設(shè)計(jì)后的不加筋板模型與目標(biāo)模型的變形位移和模態(tài)頻率都相同，達(dá)到靜、動(dòng)力學(xué)特性同時(shí)等效的目的。根據(jù)Nastran優(yōu)化設(shè)計(jì)原理和動(dòng)力學(xué)特性，改善優(yōu)化設(shè)計(jì)程序中目標(biāo)函數(shù)、響應(yīng)參數(shù)和限制條件的設(shè)置。

優(yōu)化設(shè)計(jì)變量：板厚變化T1、T2、T3

MPC約束：所有節(jié)點(diǎn)的變形位移的最小二乘最小

優(yōu)化設(shè)計(jì)目標(biāo)函數(shù)：

其中，N是板模型的節(jié)點(diǎn)數(shù)量；是目標(biāo)模型的頻率；動(dòng)力學(xué)優(yōu)化設(shè)計(jì)后板模型的頻率[8]。

把修改之后的Nastran優(yōu)化設(shè)計(jì)程序重新提交SOL200求解器求解計(jì)算，在靜力學(xué)等效的基礎(chǔ)上進(jìn)行動(dòng)力學(xué)優(yōu)化設(shè)計(jì)計(jì)算。

動(dòng)力學(xué)優(yōu)化計(jì)算結(jié)果如表5所示。

表5 動(dòng)力學(xué)優(yōu)化設(shè)計(jì)變量變化m

表6 動(dòng)力學(xué)優(yōu)化設(shè)計(jì)變形位移對(duì)比

表5中的板厚度變量T在經(jīng)過5次迭代計(jì)算后，設(shè)計(jì)變量趨于穩(wěn)定，目標(biāo)函數(shù)Φ的最小二乘最小，目標(biāo)函數(shù)取得最小值。把厚度最優(yōu)的T值分別賦予板厚，在Patran仿真分析程序中對(duì)動(dòng)力學(xué)優(yōu)化模型進(jìn)行靜、動(dòng)力學(xué)分析計(jì)算，提取數(shù)據(jù)文件。檢驗(yàn)動(dòng)力學(xué)優(yōu)化設(shè)計(jì)后不加筋板模型節(jié)點(diǎn)變形位移和模態(tài)頻率與目標(biāo)模型各節(jié)點(diǎn)變形位移和頻率是否相等。檢驗(yàn)結(jié)果如表6所示。模態(tài)頻率對(duì)比如表7所示。

表7 動(dòng)力學(xué)優(yōu)化設(shè)計(jì)模態(tài)頻率對(duì)比

其中，動(dòng)力學(xué)優(yōu)化設(shè)計(jì)板和目標(biāo)模型的前3階的彎曲和扭轉(zhuǎn)振型如圖6所示。

圖6 有限元模型振型對(duì)比圖Fig.6 Finite element model of vibration mode comparison chart

從表6、表7中可以看出，優(yōu)化設(shè)計(jì)迭代計(jì)算之后，動(dòng)力學(xué)優(yōu)化設(shè)計(jì)后的不加筋板模型的各節(jié)點(diǎn)變形位移與目標(biāo)模型各節(jié)點(diǎn)變形位移基本一致，模態(tài)頻率也基本一致。靜、動(dòng)力學(xué)特性基本達(dá)到等效，可稱此優(yōu)化設(shè)計(jì)后的板模型為靜、動(dòng)力學(xué)等效模型。

5 算例

美國NASA蘭利研究中心研發(fā)的AGARD 445.6機(jī)翼是一個(gè)國際上公認(rèn)的用于在風(fēng)洞中進(jìn)行顫振試驗(yàn)的跨音速標(biāo)準(zhǔn)顫振計(jì)算模型。其試驗(yàn)結(jié)果可以和利用仿真分析軟件計(jì)算的跨音速顫振結(jié)果進(jìn)行對(duì)比驗(yàn)證[9]。

本文選取這一標(biāo)模試驗(yàn)，驗(yàn)證Nastran動(dòng)力學(xué)優(yōu)化設(shè)計(jì)程序在模型等效處理上的有效性，進(jìn)而分析本文編制的Nastran優(yōu)化程序的有效性。

AGARD 445.6機(jī)翼翼型為NACA65A004,是具有明顯跨音速氣動(dòng)特性的變厚度薄形機(jī)翼,展長為762mm,1/4弦線的后掠角為45°,機(jī)翼展弦比(展長與平均弦長的比值)為1.62, 機(jī)翼根稍比(翼稍與翼根部的比值)為0.66。

在Nastran優(yōu)化設(shè)計(jì)程序中，圖7建立的機(jī)翼模型是殼單元。AGARD 445.6機(jī)翼用勻質(zhì)大的桃花芯木層合薄板制成，其二維模型如圖7所示。

在Nastran優(yōu)化設(shè)計(jì)程序中，圖7建立的機(jī)翼模型是殼單元，AGARD 445.6機(jī)翼用勻質(zhì)大的桃花芯木層合薄板制成，其具體材料參數(shù)如表8所示。

用桃花心木層合板制成的AGARD445.6機(jī)翼試驗(yàn)?zāi)Ｐ妥匀荒B(tài)分布范圍較寬,取前4階計(jì)算結(jié)果。

圖7 AGARD 445.6 機(jī)翼二維模型Fig.7 AGARD 445.6 two-dimensional wing model

表8 AGARD 445.6機(jī)翼的材料參數(shù)

表9 AGARD 445.6機(jī)翼的固有頻率表

表9為本文的計(jì)算結(jié)果與Goura[10]、Kolnay[11]、Ryan[12]計(jì)算結(jié)果和試驗(yàn)結(jié)果比較的數(shù)據(jù)。

AGARD 445.6機(jī)翼的固有頻率比較，如表9所示。

取機(jī)翼前4階模態(tài)的固有頻率數(shù)值與Goura、Kolonay和Ryan的計(jì)算結(jié)果進(jìn)行比較分析,本文的計(jì)算結(jié)果與試驗(yàn)值較為接近。固有模態(tài)云圖如圖8（a）～圖10（a）所示,與文獻(xiàn)[13]給出的模態(tài)圖8（b）～圖10（b）比較。

需要指出的是:圖8～圖11中的固有頻率值是文獻(xiàn)給出的計(jì)算參考結(jié)果,不是試驗(yàn)結(jié)果。所以這個(gè)值與表9中列舉的值并不相同。圖8中本文計(jì)算的模態(tài)圖8（a）是向下,文獻(xiàn)[13]給出圖8（b）則是向上,這只是相位相差180°,振型仍然是一致的。圖8～圖11中振型沒有相位差,本文計(jì)算的固有模態(tài)結(jié)果與文獻(xiàn)[13]給出的模態(tài)是吻合的[14]。

圖8 AGARD 445.6機(jī)翼第1階模態(tài)對(duì)比圖Fig.8 AGARD 445.6 wing for the first modal contrast

圖9 AGARD 445.6機(jī)翼第2階模態(tài)對(duì)比圖Fig. 9 AGARD 445.6 wing for the second modal contrast

圖10 AGARD 445.6機(jī)翼第3階模態(tài)對(duì)比圖Fig.10 AGARD 445.6 wing for the third modal contrast

圖11 AGARD 445.6機(jī)翼第4階模態(tài)對(duì)比圖Fig.11 AGARD 445.6 wing for the fourth modal contrast

6 結(jié)論

本文提出基于優(yōu)化設(shè)計(jì)的等效分析方法，在針對(duì)含離散源損傷的機(jī)翼結(jié)構(gòu)進(jìn)行靜、動(dòng)力學(xué)分析的過程中，本方法只需要知道含離散源損傷的機(jī)翼結(jié)構(gòu)的基本幾何參數(shù)即可分析機(jī)翼結(jié)構(gòu)的靜、動(dòng)力學(xué)特性。利用等效方法，使用參數(shù)化建模，內(nèi)部自動(dòng)迭代運(yùn)算，計(jì)算效率高，適合于機(jī)翼結(jié)構(gòu)初始設(shè)計(jì)階段的快速建模分析。

采用本文開發(fā)的Nanstran優(yōu)化設(shè)計(jì)程序?qū)x散源損傷機(jī)翼結(jié)構(gòu)特性進(jìn)行分析時(shí)，為準(zhǔn)確對(duì)含離散源損傷機(jī)翼結(jié)構(gòu)的仿真分析，要注意對(duì)約束、變量控制以及目標(biāo)函數(shù)的選擇和控制。

[1]杜凱，矯桂瓊，王翔. 含離散源損傷復(fù)合材料加筋板的拉伸特性[J]. 復(fù)合材料學(xué)報(bào), 2008,8:181-186.

DU Kai，JIAO Guiqiong， WANG Xiang. Tensile properties of stiffened composite panels with discretesource damage [J].Acta Materiae Compositae Sinica,2008,8:181-186.

[2]KRISHNAMURTHY T, BRIAN H. Mason.Equivalent plate analysis of aircraft wing with discrete source damage[C]// 7th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Confere,2006.AIAA 2006-2218.

[3]GB/T228.1-2010金屬材料拉伸試驗(yàn)第1部分：室溫試驗(yàn)方法[S].

GB/T228.1-2010 Test of metallic materials - Part 1 Tensile .RT Test Method [S].

[4]馬愛軍，周傳月. Patran和Nastran有限元分析專業(yè)教程[M].北京：清華大學(xué)出版社, 2005.

MA AiJun，ZHOU ChuanYue. Patran and Nastran finite element analysis professional tutorial [M]. Beijing:Tsinghua University Press, 2005.

[5]程鵬. MSC/NASTRAN優(yōu)化設(shè)計(jì)方法的討論[J]. 航天器工程,1996,5: 221-225.

CHENG Peng. The MSC/NASTRAN discuss of optimization design method [J].Spacecraft Engineering,1996,5: 221-225.

[6]KRISHNAMURTHY T, Frequency response of an aircraft wing with discrete source damage using equivalent plate analysis [C]// 8th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Confere,2007, AIAA, 2007-2144.

[7]李增剛. Nastran快速入門與實(shí)例[M].北京： 國防工業(yè)出版社，2007,6：128-130.

LI ZengGang. Nastran quick start and example[M].Beijing: National Defence Industry Press，2007.

[8]KRISHNAMURTHY T.Frequencies and flutter speed estimation for damaged aircraft wing using scaled equivalent plate analysis[C]// NASA Langley Research Center, Hampton, VA 23681, U.S.A,2010.

[9]史愛明,楊永年,葉正寅. 兩種跨聲速氣動(dòng)彈性問題分析研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2005,23: 414-418.

SHI AiMing, YANG YongNian, YE ZhengYin. Investigation of two aeroelasticity problems in transonic flow[J]. Acta Aerodynamica Sinica,2005,23: 414-418.

[10]GOURA, LAURE G S. Time marching analysis of flutter using computational fluid dynamics[D].Glasgow: University of Glasgow, 2001.

[11]Kolonay, R. M. Unsteady aeroelastic optimization in the transonic regime [D].U S A: Purdue University,1996.

[12]RYAN J,BEAUBIEN, FRED N, et al. Time and frequency domain fluttersolutions for the AGARD 445.6 wing[R]. Ottawa:Carleton University. 2005.

[13]YATES C E,AGARD standard aeroelastie configu-rations for dynamic response - AGARD 445.6 Wing[R].1985.

[14]袁鵬程，非定?？缫羲贆C(jī)翼的顫振分析[D]. 浙江:浙江大學(xué), 2011.

YUAN PengCheng，F(xiàn)lutter analysisof unsteady transonic wing [D].Zhejiang :Zhejiang University,2011.