鬼佬大哥大
  • / 21
  • 下載費用:30 金幣  

非定常氣動力最小狀態有理近似的非線性優化算法.pdf

摘要
申請專利號:

CN201510526199.6

申請日:

2015.08.25

公開號:

CN105046021A

公開日:

2015.11.11

當前法律狀態:

授權

有效性:

有權

法律詳情: 授權|||實質審查的生效IPC(主分類):G06F 17/50申請日:20150825|||公開
IPC分類號: G06F17/50 主分類號: G06F17/50
申請人: 西北工業大學
發明人: 劉祥; 孫秦; 李亮
地址: 710072陜西省西安市友誼西路127號
優先權:
專利代理機構: 西北工業大學專利中心61204 代理人: 慕安榮
PDF完整版下載: PDF下載
法律狀態
申請(專利)號:

CN201510526199.6

授權公告號:

||||||

法律狀態公告日:

2017.12.05|||2015.12.09|||2015.11.11

法律狀態類型:

授權|||實質審查的生效|||公開

摘要

一種非定常氣動力最小狀態有理近似的非線性優化算法。本發明對D矩陣的關鍵模態行初始化再順序求解E矩陣的各列和D矩陣的其余行,并充分利用當前迭代點處的函數值信息,有效提高了逆Hessian陣的近似精度,提高了外層非線性優化算法的效率。本發明有效提高了廣義氣動力矩陣中關鍵模態行元素的擬合精度,有利于提高時域顫振分析結果精度和突風響應分析結果精度。將頻域下離散的非定常氣動力系數矩陣以MS法的形式延拓至拉氏域,且能有效提高計算效率和關鍵模態項的精度。

權利要求書

1.一種非定常氣動力最小狀態有理近似的非線性優化算法,其特征在于,具體過程是:
步驟1,建立機翼有限元模型和氣動力模型,計算機翼在給定減縮頻率下的廣義氣
動力系數矩陣,具體是:
Ⅰ建立機翼有限元模型,通過Nastran軟件對建立的機翼有限元模型進行模態分
析,模態分析中取機翼的前9階彈性模態;
Ⅱ建立氣動力模型,并將得到的模態分析結果導入到Zaero軟件中,得到機翼在
給定減縮頻率下的廣義氣動力系數矩陣:
Q(ik)=F(ik)+iG(ik)(1)
其中Q(ik)為給定減縮頻率下的廣義氣動力系數矩陣,k=ωb/V為減縮頻率,
ω為機翼振蕩的圓頻率,b為機翼參考弦長,V為來流速度;F(ik)和G(ik)分別表
示廣義氣動力系數矩陣的實部和虛部;
對Q(ik)擬合時采用MS法的擬合公式,即
Q a p ( s ) = F a p ( i k ) + iG a p ( i k ) = A 0 + A 1 s + A 2 s 2 + D ( s I - R ) - 1 E · s - - - ( 2 ) ]]>
其中s為拉氏變量;表示擬合得到的廣義氣動力系數矩陣,Fap(ik)
和Gap(ik)分別表示的實部和虛部;A0∈Rn×n、A1∈Rn×n、A2∈Rn×n、D∈Rn×m、
R∈Rm×m、E∈Rm×n均為待定系數矩陣,n為結構模態數,m為氣動滯后根數量;I
為m階單位矩陣;R為由氣動滯后根組成的對角矩陣,表示為
R=diag(x)=diag([x1x2…xm])(3)
其中x表示氣動滯后根組成的向量;
步驟2,給定氣動滯后根向量的初值x0,對MS法的擬合公式進行總體擬合并計算
該總體擬合的誤差f(x0);
給定氣動滯后根向量的初值x0;
根據給定的氣動滯后根向量的初值x0得到由氣動滯后根組成的對角矩陣R;對MS
法的擬合公式進行總體擬合并計算總體擬合誤差f(x0),具體過程是:
Ⅰ取結構的第r階模態為關鍵模態,并令D矩陣的第r行元素為任意非零常數,
其中D矩陣的第r行元素代表結構的第r階模態對應的廣義氣動力;此時方程(2)
變為
Q a p , r ( s ) = A 0 , r + A 1 , r s + A 2 , r s 2 + D r ( s I - R ) - 1 E · s - - - ( 4 ) ]]>
其中下標r均表示矩陣的第r行;
Ⅱ約束方程在s=0處的有理逼近值等于氣動力系數矩陣列表值,在s=ikf處的實部
有理逼近值等于矩陣的列表值,在s=ikg處的虛部有理逼近值等于矩陣的列表值,
其中kf和kg均為指定的減縮頻率值;之后從D矩陣第r行出發擬合E矩陣各列元
素,所述E矩陣第j列的加權最小二乘求解式如下
Σ l = 1 L ( P l T W r j 2 P l + Q l T W r j 2 Q l ) E j = Σ l = 1 L { P l T W r j 2 F j ( ik l ) + Q l T W r j 2 G j ( ik l ) } - - - ( 5 ) ]]>
其中Ej表示E矩陣第j列,
P l = k l 2 D r [ ( k l 2 I + R 2 ) - 1 - ( k f 2 I + R 2 ) - 1 ] Q l = k l D r [ ( k l 2 I + R 2 ) - 1 - ( k g 2 I + R 2 ) - 1 ] R - - - ( 6 ) ]]>
W r j = 1 max l { | Q r j ( ik l ) | , 1 } - - - ( 7 ) ]]>
其中L為減縮頻率個數;Qrj(ikl)表示Q(ik)第r行第j列元素在減縮頻率kl處的值,
Wrj表示Q(ik)的加權矩陣W的第r行第j列元素;和分別表示F(ikl)和
G(ikl)的第j列;
因當前機翼在顫振點處的減縮頻率為0.07,故取kf=kg=0.05以使顫振點附近的擬
合精度更高;
Ⅲ求解出E矩陣后,再逐行用加權最小二乘法求解D矩陣的各元素,第r行除外;
D矩陣第i行的加權最小二乘求解式如下
Σ l = 1 L ( P l * T W i 2 P l * + Q l * T W i 2 Q l T ) D i T = Σ l = 1 L ( P l * T W i 2 F ^ i ( ik l ) + Q l * T W i 2 G ^ i ( ik l ) ) - - - ( 8 ) ]]>
其中Di表示D矩陣第i行,
P l * = k l 2 E T [ ( k l 2 I + R 2 ) - 1 - ( k f 2 I + R 2 ) - 1 ] Q l * = k l E T R [ ( k l 2 I + R 2 ) - 1 - ( k g 2 I + R 2 ) - 1 ] - - - ( 9 ) ]]>
Wi2表示加權矩陣W第i行各元素的平方值構成的對角矩陣,和分別表
示F(ikl)和G(ikl)的第i行;
Ⅳ計算對MS法的擬合公式進行總體擬合結果的總誤差f(x0)
f ( x 0 ) = Σ l = 1 L Σ j = 1 n Σ i = 1 n W i j 2 { [ F a p , i j ( k l ) - F i j ( k l ) ] 2 + [ G a p , i j ( k l ) - G i j ( k l ) ] 2 } - - - ( 10 ) ]]>
其中i和j分別表示各矩陣的行和列;
步驟3,開始對滯后根向量進行非線性優化;
給定初始Hessian陣的逆陣H0,為m階單位陣;令迭代下標kk=0;
步驟4,計算第一個迭代點x0處的梯度值其第i項的算式如下
( f 0 ) i = f ( x 0 + α · e i ) - f ( x 0 - α · e i ) 2 α , ( i = 1 , 2 , ... m ) - - - ( 11 ) ]]>
其中ei為第i個元素為1的m階單位向量,α為一個充分小的實數,此處取為0.001;
f(x0+α·ei)和f(x0-α·ei)的計算方法同步驟2中f(x0)的計算方法;
步驟5,確定搜索方向;
通過公式(12)確定搜索方向
d k k = - H k k f k k - - - ( 12 ) ]]>
步驟6,沿dkk線性搜索步長因子αkk;具體過程是:
Ⅰ給定參數0<ξ<0.5和0<β<1;取線性搜索時計算目標函數值的最大許可次數
為N1;令mm=0代表本次循環目標函數值的計算次數,令nn=0代表本次循環自
標量的越界次數;
Ⅱ令xkk,st=xkk+βmm+nndkk,xkk,st代表從xkk出發得到的試探點;
Ⅲ判斷 s kk , st i > x i , ( i = 1 , . . . , m ) ]]> s kk , st i < x i , ( i = 1 , . . . , m ) ]]>是否滿足,其中表示xkk,st的第
i項,xi和分別表示自變量第i項的下界和上界;若都滿足則轉至步驟Ⅳ;否則
令nn+1→nn并轉至步驟Ⅱ,其中“→”表示賦值運算;
Ⅳ計算試探點xkk,st處的目標函數值f(xkk,st);
Ⅴ判斷mm是否超過N1;若不超過,則轉到步驟Ⅵ;否則取mm為使f(xkk,st)最小
的數,令步長因子αkk=βmm+nn,fkk+1=f(xkk,st),停止線性搜索;
Ⅵ判斷試探點處Armijo不等式的滿足情況
f ( x k k , s t ) - f ( x k k ) ξβ m m + n n f k k T d k k - - - ( 13 ) ]]>
f ( x k k , s t ) - f ( x k k ) ( 1 - ξ ) β m m + n n f k k T d k k - - - ( 14 ) ]]>
若滿足,則令步長因子αkk=βmm+nn,fkk+1=f(xkk,st),停止線性搜索;否則令
mm+1→mm,轉至步驟Ⅴ繼續搜索計算,直至mm超過N1或者在試探點處滿足
Armijo不等式(13)和(14);
步驟7,求下一個迭代點xkk+1:
xkk+1=xkk+αkkdkk(15)
步驟8,計算點xkk+1處的梯度值計算方法同步驟4;
步驟9,判斷非線性優化過程是否收斂,具體是:
Ⅰ判斷||xkk||>ε6是否成立,若不成立則轉步驟Ⅱ,否則再判斷||xkk+1-xkk||/||xkk||≤ε1
是否成立,若成立則停止優化迭代,否則轉步驟Ⅲ;
Ⅱ判斷||xkk+1-xkk||≤ε2是否成立,若成立則停止優化迭代,否則轉步驟Ⅲ;
Ⅲ判斷|fkk|>ε7是否成立,若不成立則轉步驟Ⅳ,否則再判斷|fkk-fkk+1|/|fkk|≤ε3是
否成立,若成立則停止優化迭代,否則轉步驟Ⅴ;
Ⅳ判斷|fkk-fkk+1|≤ε4是否成立,若成立則停止優化迭代,否則轉步驟Ⅴ;
Ⅴ判斷是否成立,若成立則停止優化迭代,否則轉步驟10,繼續進行
優化迭代;
在上述的收斂準則中,ε1、ε2、ε3、ε4、ε5為誤差限,取為10-6;ε6、ε7分別用于
判斷||xk||和|fk|的數量級,取為10-4;
步驟10,計算Hkk+1;具體是:
Ⅰ令自變量變化為skk=xkk+1-xkk,梯度差為
Ⅱ通過改進的BFGS校正公式(16)和(17)計算逆Hessian陣Hkk+1:
H k k + 1 = ( I - s k k y k k T y k k T s k k ) H k k ( I - y k k s k k T y k k T s k k ) + s k k s k k T y k k T s k k , y k k T s k k > 0 H k k , y k k T s k k 0 - - - ( 16 ) ]]>
其中
y k k = y k k + 12 ( f k k - f k k + 1 ) + 5 f k k + 1 T s k k + ( 7 - α k k ) f k k T s k k s k k T θ k k θ k k - - - ( 17 ) ]]>
式中,表示修正的梯度差,θkk∈Rm是滿足的任意向量,此處取為skk;
步驟11,令迭代下標kk→kk+1,轉到步驟5繼續迭代,直至滿足步驟9中的收斂
條件;至此,完成了對NACA0012對稱翼型的M6機翼的非定常氣動力最小狀態
有理近似的非線性優化算法。
2.如權利要求1所述一種非定常氣動力最小狀態有理近似的非線性優化算法,其特征
在于,所述建立的機翼有限元模型中,機翼下翼面的弦向剖面共有8個網格節點分
別位于距前緣0.0%、5.0%、10.0%、27.0%、45.5%、63.5%、82.0%和100%處;沿
機翼展向各網格節點數量與翼肋的數量相同,并與各翼肋的展向位置一致。
3.如權利要求1所述一種非定常氣動力最小狀態有理近似的非線性優化算法,其特征
在于,所述建立的氣動力模型的弦向被均分為15等份,展向被均分為24等份。
4.如權利要求1所述一種非定常氣動力最小狀態有理近似的非線性優化算法,其特征
在于,在對得到的廣義氣動力系數矩陣進行擬合前,需先給定氣動滯后根向量的初
值,初始化方法為:當m=1時,x0=[-0.3];當m=2時,x0=[-0.3-0.5];當m=3時,
x0=[-0.3-0.5-0.7];當m=4時,x0=[-0.3-0.5-0.7-0.9]。

說明書

非定常氣動力最小狀態有理近似的非線性優化算法

技術領域

本發明涉及氣動彈性力學領域,具體是一種用于頻域非定常氣動力有理近似的非
線性優化方法。

背景技術

為便于氣動伺服彈性系統的多學科優化設計,需要將彈性飛行器的運動方程轉換
為時域的一階時不變狀態空間方程,而結構在任意運動下的非定常氣動力模型則是其
中一個關鍵組成部分。在亞音速或超音速情況下,以偶極子格網法為基礎計算得到的
廣義非定常氣動力系數矩陣(GAF)是一系列給定減縮頻率的函數,代表結構在簡諧
振蕩時所受的廣義氣動力。為了將廣義氣動力變換到時域空間,需要將其近似延拓為
拉氏域的有理函數形式,再經過整理,結合彈性飛行器的運動方程即可得到氣動彈性
系統的狀態空間模型。

現在工程中最常用的有理函數近似(RFA)方法均以最小二乘法為基礎,包括Roger
法、修正矩陣(MMP)法和最小狀態(MS)法等。從這些方法出發得到的狀態空間
方程中均包含由氣動滯后根產生的氣動狀態擴充項。Roger法產生的氣動擴充項數量
為結構模態數與氣動滯后根數的乘積,MMP法對應的氣動力擴充項數量為廣義氣動
力系數矩陣各列對應的氣動滯后根數量之和,而MS法對應的氣動力擴充項數等于氣
動滯后根的數量。大量的工程實踐表明,在氣動力擴充項的數量相同時,MS法的擬
合精度最高,但其中的非線性最小二乘迭代過程計算量很大。文獻1“陳青.一種建立
非定常氣動力頻域模型的簡單方法[J].空氣動力學學報,1988年,第6卷第4期”吸取
了Roger近似式中同一矩陣各元素獨立確定時精度較高的特點,對選定模態受到的非
定常氣動力進行了精確擬合,并保留了MS法的形式,改善了擬合精度,但對擬合項
無合理加權且對氣動滯后根的一維優化方法效果不佳。針對MS法計算量大、初值的
選取依賴于工程經驗等缺點,文獻2“何程.一種改進的非定常氣動力擬合方法[J].航
空學報,1993年,第14卷第7期”通過調整擬合式最后一項的形式并將其中氣動滯后
根取為大的互異負實數的方式將問題轉化為線性擬合問題,但要求滯后根的數量必須
與結構模態數一致且擬合精度不高。文獻3“E.Nissim.OntheFormulationof
Minimum-StateApproximationasaNonlinearOptimizationProblem1.JournalofAircraft,
Vol.43,No.4,2006,pp.1007–1013.”在保留MS法擬合形式的基礎上釋放了等式約束,
并通過解析法得到了自變量的梯度算法,再結合對廣義氣動力系數矩陣的比例縮放使
擬合效率和精度都得到了提高,但非線性優化項的數量遠大于氣動滯后根的個數,應
用難度較大。目前非線性優化算法已經被應用到各類擬合方法中來優化氣動滯后根以
減小擬合誤差,這就進一步對擬合方法和非線性優化算法的效率提出了要求。

發明內容

為克服已有方法精度有限或效率不高的缺點,本發明提出了一種非定常氣動力最
小狀態有理近似的非線性優化算法。

本發明的具體過程是:

步驟1,建立機翼有限元模型和氣動力模型,計算機翼在給定減縮頻率下的廣義
氣動力系數矩陣,具體是:

Ⅰ建立機翼有限元模型,通過Nastran軟件對建立的機翼有限元模型進行模態分
析,模態分析中取機翼的前9階彈性模態。所述建立的機翼有限元模型中,機翼下翼
面的弦向剖面共有8個網格節點,分別位于距前緣0.0%、5.0%、10.0%、27.0%、45.5%、
63.5%、82.0%和100%處。沿機翼展向各網格節點數量與翼肋的數量相同,并與各翼
肋的展向位置一致。

Ⅱ建立氣動力模型,并將得到的模態分析結果導入到Zaero軟件中,得到機翼在
給定減縮頻率下的廣義氣動力系數矩陣:

Q(ik)=F(ik)+iG(ik)(1)

其中Q(ik)為給定減縮頻率下的廣義氣動力系數矩陣,k=ωb/V為減縮頻率,ω
為機翼振蕩的圓頻率,b為機翼參考弦長,V為來流速度。F(ik)和G(ik)分別表示廣義
氣動力系數矩陣的實部和虛部。

對Q(ik)擬合時采用MS法的擬合公式,即

Q a p ( s ) = F a p ( i k ) + iG a p ( i k ) = A 0 + A 1 s + A 2 s 2 + D ( s I - R ) - 1 E · s - - - ( 2 ) ]]>

其中s為拉氏變量。表示擬合得到的廣義氣動力系數矩陣,Fap(ik)和
Gap(ik)分別表示的實部和虛部。A0∈Rn×n、A1∈Rn×n、A2∈Rn×n、D∈Rn×m、R∈Rm×m、
E∈Rm×n均為待定系數矩陣,n為結構模態數,m為氣動滯后根數量。I為m階單位矩
陣。R為由氣動滯后根組成的對角矩陣,表示為

R=diag(x)=diag([x1x2…xm])(3)

其中x表示氣動滯后根組成的向量。

所述建立的氣動力模型的弦向被均分為15等份,展向被均分為24等份。

步驟2,給定氣動滯后根向量的初值x0,對MS法的擬合公式進行總體擬合并計
算該總體擬合的誤差f(x0)。

在對得到的廣義氣動力系數矩陣進行擬合前,需先給定氣動滯后根向量的初值,
初始化方法為:當m=1時,x0=[-0.3];當m=2時,x0=[-0.3-0.5];當m=3時,
x0=[-0.3-0.5-0.7];當m=4時,x0=[-0.3-0.5-0.7-0.9]。

根據給定的氣動滯后根向量的初值x0得到由氣動滯后根組成的對角矩陣R;對
MS法的擬合公式進行總體擬合并計算總體擬合誤差f(x0),具體過程是:

Ⅰ取結構的第r階模態為關鍵模態,并令D矩陣的第r行元素為任意非零常數,
其中D矩陣的第r行元素代表結構的第r階模態對應的廣義氣動力。此時方程(2)變

Q a p , r ( s ) = A 0 , r + A 1 , r s + A 2 , r s 2 + D r ( s I - R ) - 1 E · s - - - ( 4 ) ]]>

其中下標r均表示矩陣的第r行。

Ⅱ約束方程在s=0處的有理逼近值等于氣動力系數矩陣列表值,在s=ikf處的實
部有理逼近值等于矩陣的列表值,在s=ikg處的虛部有理逼近值等于矩陣的列表值,
其中kf和kg均為指定的減縮頻率值。之后從D矩陣第r行出發擬合E矩陣各列元素,
所述E矩陣第j列的加權最小二乘求解式如下

Σ l = 1 L ( P l T W r j 2 P l + Q l T W r j 2 Q l ) E j = Σ l = 1 L { P l T W r j 2 F j ( ik l ) + Q l T W r j 2 G j ( ik l ) } - - - ( 5 ) ]]>

其中Ej表示E矩陣第j列,

P l = k l 2 D r [ ( k l 2 I + R 2 ) - 1 - ( k f 2 I + R 2 ) - 1 ] Q l = k l D r [ ( k l 2 I + R 2 ) - 1 - ( k g 2 I + R 2 ) - 1 ] R - - - ( 6 ) ]]>

W r j = 1 m a x l { | Q r j ( ik l ) | , 1 } - - - ( 7 ) ]]>

其中L為減縮頻率個數。Qrj(ikl)表示Q(ik)第r行第j列元素在減縮頻率kl處的值,Wrj
表示Q(ik)的加權矩陣W的第r行第j列元素。和分別表示F(ikl)和G(ikl)的
第j列。

因當前機翼在顫振點處的減縮頻率為0.07,故取kf=kg=0.05以使顫振點附近的擬
合精度更高。

Ⅲ求解出E矩陣后,再逐行用加權最小二乘法求解D矩陣的各元素,第r行除
外。D矩陣第i行的加權最小二乘求解式如下

Σ l = 1 L ( P l * T W i 2 P l * + Q l * T W i 2 Q l * ) D i T = Σ l = 1 L ( P l * T W i 2 F ^ i ( ik l ) + Q l * T G ^ i ( ik l ) ) - - - ( 8 ) ]]>

其中Di表示D矩陣第i行,

P l * = k l 2 E T [ ( k l 2 I + R 2 ) - 1 - ( k f 2 I + R 2 ) - 1 ] Q l * = k l E T R [ ( k l 2 I + R 2 ) - 1 - ( k g 2 I + R 2 ) - 1 ] - - - ( 9 ) ]]>

Wi2表示加權矩陣W第i行各元素的平方值構成的對角矩陣,和分別表示
F(ikl)和G(ikl)的第i行;

Ⅳ計算對MS法的擬合公式進行總體擬合結果的總誤差f(x0)

f ( x 0 ) = Σ l = 1 L Σ j = 1 n Σ i = 1 n W i j 2 { [ F a p , i j ( k l ) - F i j ( k l ) ] 2 + [ G a p , i j ( k l ) - G i j ( k l ) ] 2 } - - - ( 10 ) ]]>

其中i和j分別表示各矩陣的行和列。

步驟3,開始對滯后根向量進行非線性優化。

給定初始Hessian陣的逆陣H0,為m階單位陣;令迭代下標kk=0;

步驟4,計算第一個迭代點x0處的梯度值其第i項的算式如下

( f 0 ) i = f ( x 0 + α · e i ) - f ( x 0 - α · e i ) 2 α , ( i = 1 , 2 , ... m ) - - - ( 11 ) ]]>

其中ei為第i個元素為1的m階單位向量,α為一個充分小的實數,此處取為0.001。
f(x0+α·ei)和f(x0-α·ei)的計算方法同步驟2中f(x0)的計算方法。

步驟5,確定搜索方向。

通過公式(12)確定搜索方向

d k k = - H k k f k k - - - ( 12 ) ]]>

步驟6,沿dkk線性搜索步長因子αkk。具體過程是:

Ⅰ給定參數0<ξ<0.5和0<β<1。取線性搜索時計算目標函數值的最大許可次數
為N1。令mm=0代表本次循環目標函數值的計算次數,令nn=0代表本次循環自標量
的越界次數;

Ⅱ令xkk,st=xkk+βmm+nndkk,xkk,st代表從xkk出發得到的試探點;

Ⅲ判斷或是否滿足,其中表示xkk,st的
第i項,xi和分別表示自變量第i項的下界和上界。若都滿足則轉至步驟Ⅳ;否則令
nn+1→nn并轉至步驟Ⅱ,其中“→”表示賦值運算;

Ⅳ計算試探點xkk,st處的目標函數值f(xkk,st);

Ⅴ判斷mm是否超過N1。若不超過,則轉到步驟Ⅵ;否則取mm為使f(xkk,st)最
小的數,令步長因子αkk=βmm+nn,fkk+1=f(xkk,st),停止線性搜索;

Ⅵ判斷試探點處Armijo不等式的滿足情況

f ( x k k , s t ) - f ( x k k ) ξβ m m + n n f k k T d k k - - - ( 13 ) ]]>

f ( x k k , s t ) - f ( x k k ) ( 1 - ξ ) β m m + n n f k k T d k k - - - ( 14 ) ]]>

若滿足,則令步長因子αkk=βmm+nn,fkk+1=f(xkk,st),停止線性搜索;否則令mm+1→mm,
轉至步驟Ⅴ繼續搜索計算,直至mm超過N1或者在試探點處滿足Armijo不等式(13)
和(14)。

步驟7,求下一個迭代點xkk+1:

xkk+1=xkk+αkkdkk(15)

步驟8,計算點xkk+1處的梯度值計算方法同步驟4。

步驟9,判斷非線性優化過程是否收斂,具體是:

Ⅰ判斷||xkk||>ε6是否成立,若不成立則轉步驟Ⅱ,否則再判斷||xkk+1-xkk||/||xkk||≤ε1
是否成立,若成立則停止優化迭代,否則轉步驟Ⅲ;

Ⅱ判斷||xkk+1-xkk||≤ε2是否成立,若成立則停止優化迭代,否則轉步驟Ⅲ;

Ⅲ判斷|fkk|>ε7是否成立,若不成立則轉步驟Ⅳ,否則再判斷|fkk-fkk+1|/|fkk|≤ε3是
否成立,若成立則停止優化迭代,否則轉步驟Ⅴ;

Ⅳ判斷|fkk-fkk+1|≤ε4是否成立,若成立則停止優化迭代,否則轉步驟Ⅴ;

Ⅴ判斷是否成立,若成立則停止優化迭代,否則轉步驟10,繼續進行
優化迭代。

在上述的收斂準則中,ε1、ε2、ε3、ε4、ε5為誤差限,取為10-6;ε6、ε7分別用
于判斷||xk||和|fk|的數量級,取為10-4。

步驟10,計算Hkk+1。具體是:

Ⅰ令自變量變化為skk=xkk+1-xkk,梯度差為

Ⅱ通過改進的BFGS校正公式(16)和(17)計算逆Hessian陣Hkk+1:

H k k + 1 = ( I - s k k y k k T y k k T s k k ) H k k ( I - y k k s k k T y k k T s k k ) + s k k s k k T y k k T s k k , y k k T s k k > 0 H k k , y k k T s k k 0 - - - ( 16 ) ]]>

其中

y k k = y k k + 12 ( f k k - f k k + 1 ) + 5 f k k + 1 T s k k + ( 7 - α k k ) f k k T s k k s k k T θ k k θ k k - - - ( 17 ) ]]>

式中,表示修正的梯度差,θkk∈Rm是滿足的任意向量,此處取為skk。

步驟11,令迭代下標kk→kk+1,轉到步驟5繼續迭代,直至滿足步驟9中的收斂
條件。至此,完成了對NACA0012對稱翼型的M6機翼的非定常氣動力最小狀態有理
近似的非線性優化算法。

為求解MS擬合式中的各未知矩陣,傳統方法是先對整個D矩陣以規定的方式初
始化后再對D矩陣和E矩陣迭代求解直至收斂,這種方法計算效率較低。若在傳統
MS法的基礎上對滯后根進行非線性優化,則整個算法就包含了內外兩層迭代過程,
其計算時間將大大超過其他方法。而本發明在步驟2中則解決了這個問題,先對D矩
陣的關鍵模態行初始化再順序求解E矩陣的各列和D矩陣的其余行,其優點有兩個方
面:首先是避免了傳統方法中對D矩陣和E矩陣的迭代過程,節省了計算時間,有利
于外層非線性優化過程的實施;其次是有效提高了廣義氣動力矩陣中關鍵模態行元素
的擬合精度,這有利于提高時域顫振分析結果精度和突風響應分析結果精度。另外,
本發明在步驟10中采用了改進的BFGS校正公式進行改進,充分利用了當前迭代點處
的函數值信息,有效提高了逆Hessian陣的近似精度,這有利于提高外層非線性優化
算法的效率。

圖6給出了當m=1時廣義氣動力矩陣中關鍵模態行兩個元素的擬合結果,圖中共
比較了五種方法,分別為MS-Dr法、MS法、MS-opt法、Roger法和Roger-opt法。其
中“MS-Dr法”為本發明所使用的方法,表明本發明采用MS法的擬合形式,且起始于
對D矩陣第r行的初始化。“MS-opt法”表示采用現有技術的非線性算法對MS法的氣
動滯后根進行優化的方法。“Roger-opt法”表示采用現有技術的非線性算法對Roger法
的氣動滯后根進行優化的方法。圖中“Real”表示被擬合項的離散列表值。

為便于對比,“MS-opt法”和“Roger-opt法”中的非線性優化算法均與本發明中的優
化算法框架一致。從圖中可以看出,在未進行非線性優化時,Roger方法和MS法都
出現了誤差較大的情況,而在進行非線性優化后兩者的擬合精度都有明顯的改善,從
而驗證了本發明中非線性優化算法的效果。另外,本發明的精度接近優化后的Roger
方法,但因其保留了MS法的形式,因此由其出發得到的狀態空間模型能夠保持MS
法模型階數低的優點。

本發明能夠將頻域下離散的非定常氣動力系數矩陣以MS法的形式延拓至拉氏
域,且能有效提高計算效率和關鍵模態項的精度。

附圖說明

圖1是本發明的具體實施過程示意圖;

圖2是機翼的整體有限元模型;

圖3是機翼前梁腹板、后梁腹板和翼肋的有限元模型;

圖4是機翼有限元模型的弦向剖面;

圖5是機翼的氣動力模型;

圖6是部分元素的擬合結果,其中,6a是Q(ik)矩陣第一行第一列元素的擬合效果,
6b是第一行第九列元素的擬合效果;

圖7是本發明的流程圖。圖中:

1.實部;2.虛部;3.MS-Dr法;4.MS法;5.MS-opt法;6.Roger法;7.Roger-opt法;
8.Real;9.翼肋;10.前梁腹板;11.蒙皮;12.后梁腹板。

具體實施方式

本實施例是NACA0012對稱翼型的M6機翼的非定常氣動力最小狀態有理近似的
非線性優化算法。所述機翼根部完全固支。機翼的參數如下:根弦長為0.8139m,翼
尖弦長為0.4573m,展長為1.1963m,前緣后掠角為30°,后緣后掠角為15.8°,無扭
轉角,蒙皮厚度為0.002m,梁腹板厚度為0.0015m,翼肋厚度為0.0015m,材料為鋁
合金,彈性模量E=70GPa,泊松比μ=0.30,密度為ρ=2.7×103kg/m3。

本實施例的具體過程是:

步驟1,建立機翼有限元模型和氣動力模型,計算機翼在給定減縮頻率下的廣義
氣動力系數矩陣,具體是:

Ⅰ建立機翼有限元模型:在Patran軟件中建立機翼有限元模型并用Nastran軟件
進行模態分析。采用三角形和四邊形殼單元進行有限元建模,圖2為機翼的整體有限
元模型,圖3為機翼前梁腹板、后梁腹板和翼肋的有限元模型。機翼由翼肋9、梁腹
板和蒙皮11三部分組成,所述的梁腹板包括前梁腹板10和后梁腹板12。所述11個
翼肋沿展向均勻分布。所述前梁腹板10位于距離機翼前緣27.0%處,后梁腹板12位
于距離機翼前緣63.5%處。機翼下翼面的弦向剖面共有8個網格節點Zi,i=1~8。如圖
4所示,8個網格節點分別記為Z1、Z2、Z3、Z4、Z5、Z6、Z7、Z8,分別位于距前緣0.0%、
5.0%、10.0%、27.0%、45.5%、63.5%、82.0%和100%處。沿機翼展向各網格節點數
量與翼肋的數量相同,并與各翼肋的展向位置一致。通過Nastran軟件對建立的機翼
有限元模型進行模態分析,模態分析中取機翼的前9階彈性模態。

Ⅱ建立氣動力模型:在Zaero軟件中建立機翼氣動力模型并計算廣義氣動力系數
矩陣。機翼的氣動力模型如圖5所示,氣動力模型的弦向被均分為15等份,展向被均
分為24等份。將通過Nastran軟件得到的模態分析結果導入到Zaero中,得到機翼在
給定減縮頻率下的廣義氣動力系數矩陣:

Q(ik)=F(ik)+iG(ik)(1)

其中Q(ik)為給定減縮頻率下的廣義氣動力系數矩陣,k=ωb/V為減縮頻率,ω
為機翼振蕩的圓頻率,b為機翼參考弦長,V為來流速度。F(ik)和G(ik)分別表示廣義
氣動力系數矩陣的實部和虛部。

本實施例中,在獲得所述的廣義氣動力系數矩陣時,設定馬赫數為0.7,減縮頻率
為k=0.0、0.05、0.1、0.15、0.2、0.25、0.3、0.35、0.4、0.5、0.6、0.7、0.8、0.9、1.0。

對Q(ik)擬合時采用MS法的擬合公式,即

Q a p ( s ) = F a p ( i k ) + iG a p ( i k ) = A 0 + A 1 s + A 2 s 2 + D ( s I - R ) - 1 E · s - - - ( 2 ) ]]>

其中s為拉氏變量。表示擬合得到的廣義氣動力系數矩陣,Fap(ik)和
Gap(ik)分別表示的實部和虛部。A0∈Rn×n、A1∈Rn×n、A2∈Rn×n、D∈Rn×m、R∈Rm×m、
E∈Rm×n均為待定系數矩陣,n為結構模態數,m為氣動滯后根數量。I為m階單位矩
陣。R為由氣動滯后根組成的對角矩陣,可表示為

R=diag(x)=diag([x1x2…xm])(3)

其中x表示氣動滯后根組成的向量。

步驟2,給定氣動滯后根向量的初值x0,對MS法的擬合公式進行總體擬合并計
算該總體擬合的誤差f(x0)。

在對得到的廣義氣動力系數矩陣進行擬合前,需先給定氣動滯后根向量的初值,
初始化方法為:當m=1時,x0=[-0.3];當m=2時,x0=[-0.3-0.5];當m=3時,
x0=[-0.3-0.5-0.7];當m=4時,x0=[-0.3-0.5-0.7-0.9]。

根據給定的氣動滯后根向量的初值x0得到由氣動滯后根組成的對角矩陣R;對
MS法的擬合公式進行總體擬合并計算總體擬合誤差f(x0),具體過程是:

Ⅰ取結構的第r階模態為關鍵模態,并令D矩陣的第r行元素為任意非零常數,
其中D矩陣的第r行元素代表結構的第r階模態對應的廣義氣動力。此時方程(2)變

Q a p , r ( s ) = A 0 , r + A 1 , r s + A 2 , r s 2 + D r ( s I - R ) - 1 E · s - - - ( 4 ) ]]>

其中下標r均表示矩陣的第r行。

因在Zaero軟件中用P-K法進行頻域顫振分析時發現,顫振模態主要由一階彎曲
模態產生,故取一階模態為關鍵模態,即r=1,因為此模態不僅對結構的穩定性影響
最大,其對結構的響應特性也有顯著影響。

Ⅱ約束方程在s=0處的有理逼近值等于氣動力系數矩陣列表值,在s=ikf處的實部
有理逼近值等于矩陣的列表值,在s=ikg處的虛部有理逼近值等于矩陣的列表值,其
中kf和kg均為指定的減縮頻率值。之后從D矩陣第r行出發擬合E矩陣各列元素,其
中E矩陣第j列的加權最小二乘求解式如下

Σ l = 1 L ( P l T W r j 2 P l + Q l T W r j 2 Q l ) E j = Σ l = 1 L { P l T W r j 2 F j ( ik l ) + Q l T W r j 2 G j ( ik l ) } - - - ( 5 ) ]]>

其中Ej表示E矩陣第j列,

P l = k l 2 D r [ ( k l 2 I + R 2 ) - 1 - ( k f 2 I + R 2 ) - 1 ] Q l = k l D r [ ( k l 2 I + R 2 ) - 1 - ( k g 2 I + R 2 ) - 1 ] R - - - ( 6 ) ]]>

W r j = 1 m a x l { | Q r j ( ik l ) | , 1 } - - - ( 7 ) ]]>

其中L為減縮頻率個數。Qrj(ikl)表示Q(ik)第r行第j列元素在減縮頻率kl處的值,Wrj
表示Q(ik)的加權矩陣W的第r行第j列元素。和分別表示F(ikl)和G(ikl)的
第j列。

因當前機翼在顫振點處的減縮頻率為0.07,故取kf=kg=0.05以使顫振點附近的擬
合精度更高。

Ⅲ求解出E矩陣后,再逐行用加權最小二乘法求解D矩陣的各元素,第r行除外。
D矩陣第i行的加權最小二乘求解式如下

Σ l = 1 L ( P l * T W i 2 P l * + Q l * T W i 2 Q l * ) D i T = Σ l = 1 L ( P l * T W i 2 F ^ i ( ik l ) + Q l * T W i 2 G ^ i ( ik l ) ) - - - ( 8 ) ]]>

其中Di表示D矩陣第i行,

P l * = k l 2 E T [ ( k l 2 I + R 2 ) - 1 - ( k f 2 I + R 2 ) - 1 ] Q l * = k l E T R [ ( k l 2 I + R 2 ) - 1 - ( k g 2 I + R 2 ) - 1 ] - - - ( 9 ) ]]>

Wi2表示加權矩陣W第i行各元素的平方值構成的對角矩陣,和分別表示
F(ikl)和G(ikl)的第i行;

Ⅳ計算對MS法的擬合公式進行總體擬合結果的總誤差f(x0)

f ( x 0 ) = Σ l = 1 L Σ j = 1 n Σ i = 1 n W i j 2 { [ F a p , i j ( k l ) - F i j ( k l ) ] 2 + [ G a p , i j ( k l ) - G i j ( k l ) ] 2 } - - - ( 10 ) ]]>

其中i和j分別表示各矩陣的行和列。

步驟3,開始對滯后根向量進行非線性優化。

給定初始Hessian陣的逆陣H0,為m階單位陣;令迭代下標kk=0;

步驟4,計算第一個迭代點x0處的梯度值其第i項的算式如下

( f 0 ) i = f ( x 0 + α · e i ) - f ( x 0 - α · e i ) 2 α , ( i = 1 , 2 , ... m ) - - - ( 11 ) ]]>

其中ei為第i個元素為1的m階單位向量,α為一個充分小的實數,此處取為0.001。
f(x0+α·ei)和f(x0-α·ei)的計算方法同步驟2中f(x0)的計算方法。

步驟5,確定搜索方向

通過公式(12)確定搜索方向

d k k = - H k k f k k - - - ( 12 ) ]]>

步驟6,沿dkk線性搜索步長因子αkk,具體過程是:

Ⅰ給定參數0<ξ<0.5和0<β<1,所述ξ和β分別為0.3和0.7。取線性搜索時計
算目標函數值的最大許可次數為N1,本實施例中,所述N1為20。令mm=0代表本次循
環目標函數值的計算次數,令nn=0代表本次循環自標量的越界次數;

Ⅱ令xkk,st=xkk+βmm+nndkk,xkk,st代表從xkk出發得到的試探點;

Ⅲ判斷或是否滿足,其中表示xkk,st的第
i項,xi和分別表示自變量第i項的下界和上界,在此分別取為-3.0和-0.1。若都滿
足則轉至步驟Ⅳ;否則令nn+1→nn并轉至步驟Ⅱ,其中“→”表示賦值運算;

Ⅳ計算試探點xkk,st處的目標函數值f(xkk,st);

Ⅴ判斷mm是否超過N1。若不超過,則轉到步驟Ⅵ;否則取mm為使f(xkk,st)最小
的數,令步長因子αkk=βmm+nn,fkk+1=f(xkk,st),停止線性搜索;

f.判斷試探點處Armijo不等式的滿足情況

f ( x k k , s t ) - f ( x k k ) ξβ m m + n n f k k T d k k - - - ( 13 ) ]]>

f ( x k k , s t ) - f ( x k k ) ( 1 - ξ ) β m m + n n f k k T d k k - - - ( 14 ) ]]>

若滿足,則令步長因子αkk=βmm+nn,fkk+1=f(xkk,st),停止線性搜索;否則令mm+1→mm,
轉至步驟Ⅴ繼續搜索計算,直至mm超過N1或者在試探點處滿足Armijo不等式(13)
和(14)。

步驟7,求下一個迭代點xkk+1:

xkk+1=xkk+αkkdkk(15)

步驟8,計算點xkk+1處的梯度值計算方法同步驟4。

步驟9,判斷非線性優化過程是否收斂,具體是:

Ⅰ判斷||xkk||>ε6是否成立,若不成立則轉步驟Ⅱ,否則再判斷||xkk+1-xkk||/||xkk||≤ε1
是否成立,若成立則停止優化迭代,否則轉步驟Ⅲ;

Ⅱ判斷||xkk+1-xkk||≤ε2是否成立,若成立則停止優化迭代,否則轉步驟Ⅲ;

Ⅲ判斷|fkk|>ε7是否成立,若不成立則轉步驟Ⅳ,否則再判斷|fkk-fkk+1|/|fkk|≤ε3是
否成立,若成立則停止優化迭代,否則轉步驟Ⅴ;

Ⅳ判斷|fkk-fkk+1|≤ε4是否成立,若成立則停止優化迭代,否則轉步驟Ⅴ;

Ⅴ判斷是否成立,若成立則停止優化迭代,否則轉步驟10,繼續進行
優化迭代。

在上述的收斂準則中,ε1、ε2、ε3、ε4、ε5為誤差限,取為10-6;ε6、ε7分別用
于判斷||xk||和|fk|的數量級,取為10-4。

步驟10,計算Hkk+1。具體是:

Ⅰ令自變量變化為skk=xkk+1-xkk,梯度差為

Ⅱ通過改進的BFGS校正公式(16)和(17)計算逆Hessian陣Hkk+1:

H k k + 1 = ( I - s k k y k k T y k k T s k k ) H k k ( I - y k k s k k T y k k T s k k ) + s k k s k k T y k k T s k k , y k k T s k k > 0 H k k , y k k T s k k 0 - - - ( 16 ) ]]>

其中

y k k = y k k + 12 ( f k k - f k k + 1 ) + 5 f k k + 1 T s k k + ( 7 - α k k ) f k k T s k k s k k T θ k k θ k k - - - ( 17 ) ]]>

式中,表示修正的梯度差,θkk∈Rm是滿足的任意向量,此處取為skk。

所述改進的BFGS校正公式,即本實施例中的公式(16)和公式(17)2014年在
西北工業大學博士論文《大規模結構并行優化方法及其工程應用研究》中公開。

步驟11,令迭代下標kk→kk+1,轉到步驟5繼續迭代,直至滿足步驟9中的收斂
條件。至此,完成了對NACA0012對稱翼型的M6機翼的非定常氣動力最小狀態有理
近似的非線性優化算法。

圖6給出了當m=1時廣義氣動力矩陣中關鍵模態行兩個元素的擬合結果,圖中共
比較了五種方法,分別為MS-Dr法、MS法、MS-opt法、Roger法和Roger-opt法。
其中“MS-Dr法”為本發明所使用的方法,表明本發明采用MS法的擬合形式,且起始
于對D矩陣第r行的初始化。“MS-opt法”表示采用現有技術的非線性算法對MS法
的氣動滯后根進行優化的方法。“Roger-opt法”表示采用現有技術的非線性算法對
Roger法的氣動滯后根進行優化的方法。圖中“Real”表示被擬合項的離散列表值。

為便于對比,“MS-opt法”和“Roger-opt法”中的非線性優化算法均與本發明中的優
化算法框架一致。從圖中可以看出:

1.在采用相同數量滯后根的情況下,Roger方法的擬合精度高于MS法,但在未
進行非線性優化時兩種方法都會出現誤差較大的情況。而在進行非線性優化后,兩者
的擬合精度都有明顯的改善。這驗證了本發明中非線性優化算法的效果。

2.本實施例的精度接近優化后的Roger方法。但因其保留了MS法的形式,因此
由其出發得到的狀態空間模型能夠保持MS法模型階數低的優點。

表1從擬合誤差和計算效率兩個方面對比了MS-Dr法與MS法和MS-opt法。其
中fr表示廣義氣動力矩陣Q(ik)的關鍵模態行的總擬合誤差,如式(18)所示,f表示
所有元素的總擬合誤差,如式(10)所示。t表示算法的總耗時。所有結果均在一個
CPU為2.4GHz,內存4G的計算機上運行得到。

f r = Σ l = 1 L Σ j = 1 n W r j 2 { [ F a p , r j ( k l ) - F r j ( k l ) ] 2 + [ G a p , r j ( k l ) - G r j ( k l ) ] 2 } - - - ( 18 ) ]]>

從表1中可以看出MS-Dr法對關鍵模態行的擬合精度和擬合效率明顯優于MS法
和MS-opt法,且隨著氣動滯后根數量的增加,MS-Dr法的優勢不斷增大。當m=4時,
MS-Dr法的fr僅為MS-opt法的9%,且其計算時間僅為MS-opt法的0.5%。從f的對
比可以看出,MS-Dr法的總體誤差略大于MS法和MS-opt法,這是提高關鍵模態行
的擬合精度所不能避免的,但可以通過適當增加滯后根數量來彌補。

表1MS-Dr法與MS法、MS-opt法的擬合效果對比



為深入分析MS-Dr法的擬合效果,以下由其擬合結果出發,構造氣動彈性系統的
狀態空間方程,并用根軌跡法進行顫振分析。具體做法參考文獻4“Karpel,M.,Time
DomainAeroservoelasticModelingUsingWeightedUnsteadyAerodynamicForces,
JournalofGuidance,Control,andDynamics,Vol.13,No.1,1990,pp.30–37.”。同時,由
MS-opt法和Roger-opt法擬合結果得到的狀態空間方程也被用來進行顫振分析以進行
比較,前者狀態空間模型的構建方法同參考文獻4,后者參考文獻5“Tiffany,S.H.,and
Adams,W.M.,Non-LinearProgrammingExtensionstoRationalFunctionApproximation
MethodsforUnsteadyAerodynamicForces,NASATP-2776,July1988.”。

對當前模型P-K法計算得到的顫振速度為445.8m/s,顫振頻率為13.69Hz。下文將
以此為參考結果來檢驗擬合精度。表2給出了這3種方法時域狀態空間模型的顫振分
析結果,其中Vf表示顫振速度誤差,ωf表示顫振頻率誤差。從表2中可以看出,雖
然氣動滯后根數量的增加會提高模型的整體擬合精度,但由Roger-opt法和MS-opt
法得到顫振結果卻并沒有明顯改善,反而可能出現隨滯后根增加誤差明顯增大的情
況。這是因為整體模型精度的提高并不能保證關鍵模態相關項精度的改善,而這對模
型的穩定性分析和動態響應分析是至關重要的。但MS-Dr法則能一直保持較高的精
度,且隨著氣動滯后根數量的增加誤差有明顯減小的趨勢。

表2MS-Dr法與MS法、MS-opt法的顫振分析誤差對比



為進行氣動伺服彈性系統的控制律設計、時域仿真和迭代優化設計,系統的運動
方程需要從頻域模型轉換為一階時域狀態空間模型,而其中的關鍵技術就是將頻域的
非定常氣動力復延拓到拉氏域,本實施例對當前廣泛使用的MS法提出了進一步的改
進,在明顯提高計算效率和關鍵模態相關項的擬合精度的同時保證了模型的整體精度,
能夠準確地計算出模型的顫振特征,在氣動伺服彈性系統的分析設計中具有重要意義。

關 鍵 詞:
非定常 氣動力 最小 狀態 有理 似的 非線性 優化 算法
  專利查詢網所有資源均是用戶自行上傳分享,僅供網友學習交流,未經上傳用戶書面授權,請勿作他用。
關于本文
本文標題:非定常氣動力最小狀態有理近似的非線性優化算法.pdf
鏈接地址:http://www.wwszu.club/p-6401572.html
關于我們 - 網站聲明 - 網站地圖 - 資源地圖 - 友情鏈接 - 網站客服客服 - 聯系我們

[email protected] 2017-2018 zhuanlichaxun.net網站版權所有
經營許可證編號:粵ICP備17046363號-1 
 


收起
展開
鬼佬大哥大