本發(fā)明涉及一種基于傳遞函數(shù)的細(xì)長(zhǎng)體顫振速度確定方法，屬于飛行器氣動(dòng)彈性
技術(shù)領(lǐng)域：
。
背景技術(shù)：
：對(duì)于飛機(jī)機(jī)身、導(dǎo)彈彈體、機(jī)翼外掛火箭等，當(dāng)長(zhǎng)細(xì)比較大時(shí)(一般大于10左右)，稱為細(xì)長(zhǎng)體。細(xì)長(zhǎng)體顫振是指發(fā)生在飛行器飛行中的動(dòng)不穩(wěn)定性，此時(shí)的飛行速度稱為顫振速度。在飛行達(dá)到顫振速度時(shí)所發(fā)生的自激振動(dòng)，大多數(shù)都會(huì)造成災(zāi)難性的后果。顫振分析就是求出顫振發(fā)生的條件，亦即求出顫振速度，并尋求在飛行器飛行速度范圍內(nèi)避免顫振發(fā)生的措施，進(jìn)而尋找提高顫振臨界速度的方法。進(jìn)行細(xì)長(zhǎng)體顫振分析和計(jì)算時(shí)，一方面需要知道細(xì)長(zhǎng)體結(jié)構(gòu)的動(dòng)力學(xué)特性，另一方面需要知道細(xì)長(zhǎng)體結(jié)構(gòu)表面非定常流動(dòng)的空氣動(dòng)力特性。目前，工程中為了避免分析計(jì)算過(guò)于復(fù)雜，在細(xì)長(zhǎng)體結(jié)構(gòu)動(dòng)力學(xué)計(jì)算時(shí)盡量減小細(xì)長(zhǎng)體振動(dòng)自由度，通常忽略高階模態(tài)，僅利用細(xì)長(zhǎng)體振動(dòng)的低階模態(tài)來(lái)近似表示細(xì)長(zhǎng)體的彎曲振動(dòng)位移。然后，再結(jié)合細(xì)長(zhǎng)體準(zhǔn)定常氣動(dòng)理論，求解細(xì)長(zhǎng)體顫振速度。實(shí)際上，細(xì)長(zhǎng)體是無(wú)限多自由度的彈性體，在顫振分析中也沒(méi)有對(duì)自由度的限制。因此，上述計(jì)算方法的缺點(diǎn)是：為了避免過(guò)于復(fù)雜的計(jì)算而減少了細(xì)長(zhǎng)體振動(dòng)的自由度，對(duì)細(xì)長(zhǎng)體實(shí)際振動(dòng)作出了近似，使計(jì)算結(jié)果的準(zhǔn)確性降低了。技術(shù)實(shí)現(xiàn)要素：本發(fā)明所要解決的技術(shù)問(wèn)題是提供一種計(jì)算方法簡(jiǎn)捷、計(jì)算結(jié)果較準(zhǔn)確的基于傳遞函數(shù)的細(xì)長(zhǎng)體顫振速度確定方法。本發(fā)明解決其技術(shù)問(wèn)題所需采用的技術(shù)方案：一種基于傳遞函數(shù)法的細(xì)長(zhǎng)體顫振速度確定方法，所述方法利用梁的橫向彎曲振動(dòng)微分方程和細(xì)長(zhǎng)體準(zhǔn)定常氣動(dòng)力模型，得出細(xì)長(zhǎng)體顫振微分方程，然后對(duì)所述細(xì)長(zhǎng)體離散為細(xì)長(zhǎng)體單元，在單元內(nèi)對(duì)細(xì)長(zhǎng)體單元顫振微分方程進(jìn)行Fourier變換，再運(yùn)用傳遞函數(shù)法求解細(xì)長(zhǎng)體顫振速度；一、所述方法需要依據(jù)的計(jì)算公式：a.利用梁的橫向彎曲振動(dòng)微分方程和細(xì)長(zhǎng)體準(zhǔn)定常氣動(dòng)力模型建立細(xì)長(zhǎng)體顫振微分方程：本發(fā)明利用梁的橫向彎曲振動(dòng)微分方程式(1)來(lái)描述細(xì)長(zhǎng)體的彎曲振動(dòng)：∂2(EI∂2h∂y2)∂y2+m∂2h∂t2+Δp=0---(1)]]>式中，h為細(xì)長(zhǎng)體彎曲振動(dòng)位移，單位為米，EI為細(xì)長(zhǎng)體抗彎剛度，單位牛頓·米2，m為細(xì)長(zhǎng)體單位長(zhǎng)度質(zhì)量，單位為千克，△p為細(xì)長(zhǎng)體單位長(zhǎng)度的氣動(dòng)力，單位為牛頓/米，y為細(xì)體長(zhǎng)軸向坐標(biāo)，單位為米，t為時(shí)間，單位為秒，為細(xì)長(zhǎng)體彎曲振動(dòng)位移h對(duì)機(jī)軸向坐標(biāo)y的二階偏導(dǎo)數(shù)，為細(xì)長(zhǎng)體彎曲振動(dòng)位移h對(duì)時(shí)間t的二階偏導(dǎo)數(shù)；根據(jù)細(xì)長(zhǎng)體準(zhǔn)定常氣動(dòng)力模型，得到細(xì)長(zhǎng)體單位長(zhǎng)度的氣動(dòng)力△p為下式(2)：Δp=-{ρA(∂2h∂t2+V∂2h∂t∂y)+ρV∂A∂y(∂h∂t+V∂h∂y)+ρAV(∂2h∂t∂y+V∂2h∂y2)}---(2)]]>式中，V為空速，單位為米/秒，A為細(xì)長(zhǎng)體橫截面積，單位為平方米，ρ為空氣密度，單位為千克/立方米；將(2)式代入(1)式得到細(xì)長(zhǎng)體顫振微分方程式(3)：∂2(EI∂2hay2)∂y2+m∂2h∂t2-ρA(∂2h∂t2+V∂2h∂t∂y)-ρV∂A∂y(∂h∂t+V∂h∂y)-ρAV(∂2h∂t∂y+V∂2h∂y2)=0---(3)]]>b.沿細(xì)長(zhǎng)體軸向劃分n個(gè)細(xì)長(zhǎng)體單元，建立細(xì)長(zhǎng)體單元的顫振微分方程：一般地，細(xì)長(zhǎng)體的物理參數(shù)EI、m、A等沿其軸向是變化的，為了方便求解，按照有限元法，沿細(xì)長(zhǎng)體軸向劃分n個(gè)單元，為了使細(xì)長(zhǎng)體單元的形式統(tǒng)一顫振微分方程，令：ξ=y-yjlj---(4)]]>式中，yj為第j個(gè)細(xì)長(zhǎng)體單元前節(jié)點(diǎn)距離細(xì)長(zhǎng)體前端點(diǎn)的距離，lj為第j個(gè)細(xì)長(zhǎng)體單元的長(zhǎng)度，ξ為細(xì)長(zhǎng)體單元內(nèi)軸向無(wú)量綱坐標(biāo)，由(4)式可知ξ∈(0,1)；當(dāng)劃分單元足夠多時(shí)，假設(shè)EI、m，A在每個(gè)細(xì)長(zhǎng)體單元內(nèi)近似線性變化是可理的，則在第j個(gè)細(xì)長(zhǎng)體單元內(nèi)有如下線性表達(dá)式：EI(ξ)=EI(0)+EI(1)-EI(0)ljξm(ξ)=m(0)+m(1)-m(0)ljξA(ξ)=A(0)+A(1)-A(0)ljξ---(5)]]>式中，EI(0)、m(0)、A(0)為EI、m、A在第j個(gè)細(xì)長(zhǎng)體單元前節(jié)點(diǎn)處的值，EI(1)、m(1)、A(1)為EI、m、A在第j個(gè)細(xì)長(zhǎng)體單元后節(jié)點(diǎn)處的值，EI(ξ)、m(ξ)、A(ξ)為EI、m、A在第j個(gè)細(xì)長(zhǎng)體單元內(nèi)ξ處的值；將(4)式與(5)式代入(3)式，可得到細(xì)長(zhǎng)體單元的顫振微分方程：EI(ξ)lj4∂4h∂ξ4+2(EI(1)-EI(0))lj4∂3h∂ξ3-ρA(ξ)V2lj2∂2h∂ξ2-2ρA(ξ)Vlj∂2h∂t∂ξ-ρV2A(1)-A(0)lj2∂h∂ξ+[m(ξ)-ρA(ξ)]∂2h∂t2-ρVA(1)-A(0)lj∂h∂t=0---(6)]]>c.利用傳遞函數(shù)法求解細(xì)長(zhǎng)體顫振速度①對(duì)(6)式中有關(guān)時(shí)間t的變量作Fourier變換，可得到下式(7)：EI(ξ)lj4∂4h∂ξ4+2(EI(1)-EI(0))lj4∂2h∂ξ3-ρA(ξ)V2lj2∂2h∂ξ2-[iω2ρA(ξ)Vlj+V2ρA(1)-ρA(0)lj2]∂h∂ξ+[(iω)2[m(ξ)-ρA(ξ)]-(iω)VρA(1)-ρA(0)lj]h=0---(7)]]>式中，虛數(shù)ω為角頻率，單位為弧度/秒；②將細(xì)長(zhǎng)體單元的顫振微分方程改寫(xiě)為狀態(tài)空間方程形式；為了便于應(yīng)用傳遞函數(shù)理論，定義一個(gè)狀態(tài)變量的向量ηe(ξ,ω)如下式(8)：ηe(ξ,ω)=h∂h∂y∂2h∂y2∂3h∂y3T---(8)]]>式中，T表示向量轉(zhuǎn)置；基于狀態(tài)變量向量ηe(ξ,ω)，可將(7)式寫(xiě)成狀態(tài)方程的形式如下式(9)：∂ηe(ξ,ω)∂ξ=Fe(ξ,ω,V)ηe(ξ,ω)+ge(ξ,ω)---(9)]]>式中，F(xiàn)e(ξ,ω,V)、ge(ξ,ω)的表達(dá)式可根據(jù)(7)式得出，具體如下Fe(ξ,ω,V)=010000100001A4A3A2A1ge(y,ω)=0---(10)]]>式中，A1=-2EI(1)-EI(0)EI(ξ)A2=ρA(ξ)V2EI(ξ)lj2A3=iω2ρA(ξ)VEI(ξ)lj2+V2ρA(1)-ρA(0)EI(ξ)lj2A4=-(iω)2[m(ξ)-ρA(ξ)]EI(ξ)lj4+(iω)VρA(1)-ρA(0)EI(ξ)lj3---(11)]]>根據(jù)傳遞函數(shù)理論，方程(9)式的標(biāo)準(zhǔn)解為：ηe(ξ,ω)＝He(ξ,ω,V)γe(ω)+de(ξ,ω)(12)式中，γe(ω)為由邊界條件給定的位移或力組成的向量，可根據(jù)傳遞函數(shù)方法理論得到表達(dá)式如下：γe(ω)=h(0,ω)∂h(0,ω)∂ξh(1,ω)∂h(1,ω)∂ξ---(13)]]>式中，He(ξ,ω,V)、fe(ξ,ω)的表達(dá)式也可根據(jù)傳遞函數(shù)方法理論得到，如下：He(ξ,ω,V)=ΦF(ξ,0,ω,V)[MbΦF(0,0,ω,V)+NbΦF(1,0,ω,V)]-1de(ξ,ω)=∫01G(ξ,ζ,ω,V)ge(ζ,ω)dζG(ξ,ζ,ω,V)=He(ξ,ω,V)MbΦF(0,ζ,ω,V)ζ<ξ-He(ξ,ω,V)NbΦF(1,ζ,ω,V)ζ>ξΦF(ξ,ζ,ω,V)=eFe(ξ,ω,V)ζ---(14)]]>式中，變量ζ∈(0,1)，Mb為細(xì)長(zhǎng)體單元前端邊界條件選擇矩陣，Nb為細(xì)長(zhǎng)體單元后端邊界條件選擇矩陣，可根據(jù)傳遞函數(shù)方法理論得到，如下：Mb=1000010000000000,Nb=0000000010000100---(15)]]>由于(10)式中給出ge(y,ω)＝0，由(14)式可知de(ξ,ω)＝0，從而(12)式可簡(jiǎn)化為：ηe(ξ,ω)＝He(ξ,ω,V)γe(ω)(16)③根據(jù)細(xì)長(zhǎng)體內(nèi)力關(guān)系建立單元平衡方程；由材料力學(xué)梁的彎矩和剪力公式可知，細(xì)長(zhǎng)體截面上彎矩和剪力的表達(dá)式可寫(xiě)為：M(y,ω)=EI∂2h∂y2F(y,ω)=∂(EI∂2h∂y2)∂y---(17)]]>在細(xì)長(zhǎng)體單元內(nèi)，并利用EI在單元內(nèi)的線性變化假設(shè)，(17)式可寫(xiě)為：σe(ξ,ω)=M(ξ,ω)F(ξ,ω)=EI(ξ)li2∂2h∂ξ2EI(ξ)li3∂3h∂ξ3+EI(1)-EI(0)li3∂2h∂ξ2---(18)]]>引入狀態(tài)變量向量(8)式，從而(18)可寫(xiě)成如下矩陣形式：σe(ξ,ω)＝Eeη(ξ)ηe(ξ,ω)(19)式中，Eeη(ξ)可根據(jù)(18)式得到其表達(dá)式如下：Eeη(ξ)=00EI(ξ)li2000EI(1)-EI(0)li3EI(ξ)li3---(20)]]>將(16)式代入(19)式，可得到σe(ξ,ω)＝Eeη(ξ)He(ξ,ω,V)γe(ω)(21)式中，γe(ω)可視為細(xì)長(zhǎng)體單元兩節(jié)點(diǎn)處的廣義位移，如果ξ取0和1，則可得到細(xì)長(zhǎng)體單元兩節(jié)點(diǎn)的廣義內(nèi)力與廣義位移之間的關(guān)系，從而建立單元平衡方程，即：令ξ＝0，ξ＝1，則有σe(ω)=σe(0,ω)σe(1,ω)=Eeη(0)He(0,ω,V)Eeη(1)He(1,ω,V)γe(ω)---(22)]]>上式即為表征細(xì)長(zhǎng)體單元節(jié)點(diǎn)廣義內(nèi)力與廣義位移關(guān)系的單元平衡方程，從而可視為單元的廣義剛度矩陣，可記為：Ke(ω,V)=Eeη(0)He(0,ω,V)Eeη(1)He(1,ω,V)---(23)]]>④將細(xì)長(zhǎng)體單元平衡方程組裝形成細(xì)長(zhǎng)體整體平衡方程；按照有限元法進(jìn)行單元組裝，可得到細(xì)長(zhǎng)體整體平衡方程為K(ω,V)γ(ω)＝F(ω)(24)式中，K(ω,V)可視為整體剛度矩陣，γ(ω)可視為整體節(jié)點(diǎn)位移向量，F(xiàn)(ω)為各單元節(jié)點(diǎn)內(nèi)力拼裝成的向量。⑤確定顫振速度；對(duì)于細(xì)長(zhǎng)體，顫振時(shí)不需慮重力，本專利中將細(xì)長(zhǎng)體所受氣動(dòng)力的影響囊括在剛度矩陣K(ω,V)，除此以外細(xì)長(zhǎng)體不受其它力的作用，從而外載荷為零，根據(jù)單元節(jié)點(diǎn)內(nèi)力與外載荷平衡，可得出F(ω)＝0(25)當(dāng)細(xì)長(zhǎng)體顫振時(shí)，γ(ω)應(yīng)有非零解，此時(shí)整體平衡方程必須滿足條件為det[K(ω,V)]＝0(26)由于K(ω,V)為復(fù)矩陣，其行列式值等于零的必要條件為矩陣行列式值的實(shí)部與虛部均為零，即Re{det[K(ω,V)]}=0Im{det[K(ω,V)]}=0---(27)]]>矩陣K(ω,V)中有空速V和圓頻率ω兩個(gè)變量，而(27)式恰好有兩個(gè)方程，可以定解。求解上述方程組時(shí)，可能會(huì)得到滿足方程組多個(gè)解，即存在多組(V,ω)能滿足方程組(27)式。根據(jù)細(xì)長(zhǎng)體顫振時(shí)在某一空速時(shí)由穩(wěn)定轉(zhuǎn)變?yōu)椴环€(wěn)定，因而，空速V最小的一組解(V,ω)應(yīng)為細(xì)長(zhǎng)體的顫振速度和相應(yīng)的顫振頻率。將得到滿足方程組的空速V和圓頻率ω，分別記為Vcz和ωcz。其中，Vcz即為細(xì)長(zhǎng)體的顫振速度，ωcz即為細(xì)長(zhǎng)體的顫振圓頻率。二、所述方法的具體步驟如下：步驟(一)：將細(xì)長(zhǎng)體劃分為n個(gè)單元，測(cè)量細(xì)長(zhǎng)體各單元節(jié)點(diǎn)的下列物理參數(shù)：細(xì)長(zhǎng)體各單元軸向長(zhǎng)度[l1l2…lj…ln]，單位為米；細(xì)長(zhǎng)體單元節(jié)點(diǎn)處的橫截面積[A1A2…Aj…An+1]，單位為米2細(xì)長(zhǎng)體單位節(jié)點(diǎn)處單位長(zhǎng)度質(zhì)量[m1m2…mj…mn+1]，單位為千克/米；細(xì)長(zhǎng)體單位節(jié)點(diǎn)處抗彎剛度[EI1EI2…EIj…EIn+1]，單位為牛頓·米2；空氣密度ρ，單位為千克/米3；步驟(二)：確定細(xì)長(zhǎng)體飛行的空速V和圓頻率ω的大致范圍，假設(shè)為步驟(三)：在范圍內(nèi)劃分合適步長(zhǎng)△V和△ω，并進(jìn)行離散，空速V和圓頻率ω的取值為：步驟(四)：取空速V＝V0，圓頻率ω依次取ω0+j△ω(j＝0,1,2,3…)，將空速V和圓頻率ω的取值與步驟(一)中的細(xì)長(zhǎng)體物理參數(shù)代入(11)式，計(jì)算各單元系數(shù)A1、A2、A3、A4的值；步驟(五)：將系數(shù)A1、A2、A3、A4的值代入(10)式，得到細(xì)長(zhǎng)體各單元矩陣Fe(ξ,ω,V)的值；步驟(六)：將矩陣Fe(ξ,ω,V)的值代入(14)式，結(jié)合(15)式，計(jì)算He(ξ,ω,V)的值；步驟(七)：將He(ξ,ω,V)的值代入(23)式，結(jié)合(20)式，計(jì)算單元?jiǎng)偠染仃嘖e(ω,V)的值；步驟(八)：利用Ke(ω,V)組裝整體剛度矩陣K(ω,V)；步驟(九)：計(jì)算Re{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；步驟(十)：再依次取空速V＝V0+j△V(j＝1,2,3…)，重復(fù)步驟(四)至步驟(十)，計(jì)算Re{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；步驟(十一)：確定滿足(27)式的空速V的值，即可確定細(xì)長(zhǎng)體的顫振速度Vcz。本發(fā)明的有益效果如下：(1)由于本發(fā)明利用梁的橫向彎曲振動(dòng)微分方程來(lái)準(zhǔn)確描述細(xì)長(zhǎng)體的橫向振動(dòng)，而沒(méi)有采用梁橫向振動(dòng)的低階模態(tài)來(lái)近似描述細(xì)長(zhǎng)體振動(dòng)，所以計(jì)算的細(xì)長(zhǎng)體顫振速度更加準(zhǔn)確。(2)本發(fā)明求解細(xì)長(zhǎng)體顫振速度的方法與現(xiàn)有技術(shù)相比更加簡(jiǎn)捷。附圖說(shuō)明圖1為細(xì)長(zhǎng)體正視圖；圖2為細(xì)長(zhǎng)體截面圖；圖3為細(xì)長(zhǎng)體結(jié)構(gòu)示意圖；圖4為細(xì)長(zhǎng)體單元剖分示意圖；圖5為實(shí)施例1的Re[detA]和Im[detA]的等值線圖(V∈(0,50),ω∈(0,200π))、在附圖中：1細(xì)長(zhǎng)體、2第i個(gè)細(xì)長(zhǎng)體、3前節(jié)點(diǎn)、4后節(jié)點(diǎn)。具體實(shí)施方式如附圖1-5所示，本發(fā)明的實(shí)施例1如下：某細(xì)長(zhǎng)體長(zhǎng)度為1.8米，重量22千克，最大截面直徑為0.025平方米。本實(shí)施例1的具體計(jì)算步驟如下：步驟(一)：將細(xì)長(zhǎng)體劃分為10個(gè)單元，測(cè)量細(xì)長(zhǎng)體各單元節(jié)點(diǎn)的下列物理參數(shù)：各單元軸向長(zhǎng)度：l1＝l2＝…＝l10＝0.18，單位為米；各單元節(jié)點(diǎn)處的橫截面積：A1＝0,A2＝A3＝…＝A10＝0.025,A11＝0.02，單位為米2各單元節(jié)點(diǎn)處的單位長(zhǎng)度質(zhì)量m1＝10,m2＝m3…＝m10＝20，單位為千克/米；各單元節(jié)點(diǎn)處的抗彎剛度EI1＝0.5×103,EI2＝…＝EI11＝1.1×103，單位為牛頓·米2；空氣密度ρ＝1.225，單位為千克/米3；步驟(二)：確定飛機(jī)機(jī)翼的空速V和圓頻率ω的大致范圍為劃分步長(zhǎng)依據(jù)
發(fā)明內(nèi)容部分中步驟(三)至步驟(十一)，計(jì)算Re[detA]和Im[detA]的值。圖3給出了在范圍內(nèi)Re[detA]和Im[detA]的等值線圖。從圖3中可以發(fā)現(xiàn)，在給定范圍內(nèi)存在滿足(26)式的點(diǎn)A，其坐標(biāo)為(14.7,3020)，它表示細(xì)長(zhǎng)體顫振速度Vcz＝3020m/s，顫振頻率為14.7Hz(顫振圓頻率ωcz＝92.32rad/s)。當(dāng)前第1頁(yè)1 2 3