![大展弦比飞行器气动弹性分析的传递函数方法](https://wfqqreader-1252317822.image.myqcloud.com/cover/823/36511823/b_36511823.jpg)
2.3 传递函数方法求解
2.3.1 机翼沿展向离散
对于一般机翼,其物理参数和几何参数沿机翼展向是变化的,因而借鉴有限元方法的思想,将机翼沿展向划分为若干段,每段为一个机翼单元,即机翼离散,如图2-3所示。当划分为足够多的单元时,在机翼单元内假设机翼的物理参数和几何参数线性变化是合理。当然,进一步划分更多单元时,可以在单元内假设机翼的物理参数和几何参数是常数。首先,定义单元局部坐标系如下:
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-019-2.jpg?sign=1739452230-e0VKYyYnTRG5eFlJJwhDtlkZxP5v6MER-0-24ad706a5469b80761888149d03d6cba)
从而,在单元内机翼的物理参数变化规律为
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-019-3.jpg?sign=1739452230-NeINKpzfYUaOPaqIXXjwXZ8rJECkgywi-0-1ff1c0cc4702fc9dfc7e0c8166fa3aa3)
其中,ξ为单元局部坐标系的坐标;lj为第j个机翼单元的长度。
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-019-4.jpg?sign=1739452230-sdlrHo6ANjU1nfnEUuO0hMm6zBfRNfS2-0-4f5d1f112969333866106d3315f8d39a)
图2-3 机翼离散
除此之外,GJ、m、Iα、b、xα、等机翼的物理参数和几何参数在机翼单元内的变化也可写为式(2-5)的形式。
进一步,联立式(2-3)与式(2-5),机翼单元颤振的控制方程可写为
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-020-2.jpg?sign=1739452230-pKXSOzEU9SYCtjaqkdJLCg3xjZoLZ59V-0-93318f86fab42f4bab6c24108e1f122b)
对于整个机翼而言,可将其视为一根悬臂梁,两端的边界条件如下:
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-020-3.jpg?sign=1739452230-fDOE6lGlOocAlIC96WYA63CShQr6MKp5-0-5e23a4ec15bcbc31016ae650d082950e)
其中,M(L, t)、F(L, t)、Q(L, t)分别为机翼自由端的弯矩、剪力和扭矩。
假设初始条件为
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-020-4.jpg?sign=1739452230-vLHsEj4K5RorFxBwULhAnhw3pqBoQozB-0-7613aa323d8ba430f18c5b6c9a2778d3)
2.3.2 机翼颤振特性求解
本节利用传递函数方法求解机翼颤振速度。对式(2-6)进行Fourier变换,并利用初始条件式(2-8),整理可得
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-021-1.jpg?sign=1739452230-8hoT6PyyZFnjL2SmlFGLX57KOa8EzLVO-0-c611a3454d15f6cf9ce9a56dc888efbb)
其中,分别表示h、α的Fourier变换;ω为圆频率;
。系数An(ξ)、Bn(ξ, iω, V)、Cn(ξ, iω, V)(n=1, 2)的表达式如下:
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-021-4.jpg?sign=1739452230-TduD6O4u4WqrCCygogllcGQkpYMH2HH6-0-35c7dbe69ab52134f0b28e5a752a97ad)
根据传递函数方法[35],定义状态变量向量如下:
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-022-1.jpg?sign=1739452230-FWHsUkpDjImcvJwRNEf176y04XfGW5Fd-0-27cb9b4a27c0d40b82d348a47b2cc1a1)
其中,T表示向量转置。
从而,式(2-8)可写成如下状态空间方程的形式:
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-022-2.jpg?sign=1739452230-GBSPONwO0OYsZ6Sne2gY4VA8b3tA5li7-0-9c407e3a4061ed177cf072e9a6615230)
其中,
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-022-3.jpg?sign=1739452230-FLGVMThdfL6i4ULzhoXjPQeTfUKGzobE-0-6c9670d073863b68d474b14db9611020)
由于式(2-9)为齐次微分方程组,故式(2-12)中ge(ξ, ω)=0。
根据传递函数方法[35],机翼单元两端的变形协调条件可以写为如下选择矩阵形式:
Mbηe(0, ω)+Nbηe(1, ω)=γe(ω) (2-15)
其中,
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-022-4.jpg?sign=1739452230-8hT2OfkMsBPcZSwfksRf1epkApGcAdGn-0-ad5806302b34b8b516d52f00b6480edb)
其中,Mb为机翼根端边界条件选择矩阵;Nb为机翼尖端边界条件选择矩阵;γe(ω)为由机翼单元节点位移构成的向量。
根据传递函数理论[35],式(2-12)的解可写为
ηe(ξ, ω)=He(ξ, ω, V)γe(ω)+fe(ξ, ω) (2-17)
其中,
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-023-1.jpg?sign=1739452230-Kzh9juH9qi08GKeRWZt2yaEWfviiHEIm-0-3f22964df33c81973ffd5a55d35a8a20)
其中,G(ξ, ζ, ω, V)为状态空间方程的域内传递函数;He(ξ, ω, V)为状态空间方程的边界传递函数。
由于ge(ζ, ω)=0,故式(2-17)中fe(ξ, ω)=0,矩阵ΦF(ξ, ζ, ω, V)可采用高斯积分求值。
对于整个机翼,其总体求解方程可借鉴有限元方法的思想进行组集。具体推导如下。
机翼单元横截面上的内力σe(ξ, ω)为
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-023-2.jpg?sign=1739452230-5M1pHyU0ryZxCYTD3QH6lCMo8bUsybpW-0-b70d226944af1a0c8b85aa064419a1bc)
通过引入前先定义的状态向量ηe(ξ, ω), 可得
σe(ξ, ω)=Eeη(ξ)ηe(ξ, ω) (2-23)
其中,。
将式(2-17)代入式(2-23),并结合fe(ξ, ω)=0, 可得
σe(ξ, ω)=Eeη(ξ)He(ξ, ω, V)γe(ω) (2-24)
令ξ=0和ξ=1,由式(2-24)可得机翼单元两节点处的内力表达式,即
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-024-1.jpg?sign=1739452230-9Fs0a5BXj9HZiDsmN4l8vBJk2IkY80vE-0-245109e6254870d21b7f611371026986)
式(2-25)与有限元法中单元节点力的表达式十分相似。其中,可视为机翼单元广义刚度矩阵;γe(ω)可视为机翼单元广义节点位移向量。为了方便描述,令
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-024-3.jpg?sign=1739452230-jWxJP7qGWS0g2AFjVbJdyX74Ae3jLag9-0-e81023f2f2084f4387833ee83e947d5b)
因而,可按照有限元组集方法对各单元节点进行统一编号,并进行组集拼接,可得机翼整体的平衡方程
K(ω, V)γ(ω)=F(ω) (2-27)
其中,K(ω, V)可视为整体刚度矩阵;γ(ω)可视为整体节点位移向量;F(ω)为由各单元节点内力拼装成的向量。本书将机翼的气动力与机翼作为一个完整的系统来考虑,且忽略重力,除此之外,机翼没有受到其他外力作用,因此根据单元节点内力与外载荷平衡,可得出F(ω)=0。
根据机翼约束条件,按照有限元方法对整体刚度矩阵K(ω, V)进行边界条件处理。
当机翼颤振时,γ(ω)应有非零解,此时须满足条件
det[K(ω, V)]=0 (2-28)
由于K(ω, V)为复矩阵,所以其行列式值等于零的必要条件为矩阵行列式值的实部与虚部均为零,即
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-024-4.jpg?sign=1739452230-1HKaHMIGg5GvF3jZsDLBx7hoiU9mY8lV-0-cc494f1c1203e1cdde8bce8c95729a86)
矩阵K(ω, V)中有空速V和圆频率ω两个变量,而式(2-29)恰好有两个方程,可以有定解。求解上述方程组时,可能会得到多个解,即存在多组(V, ω)满足式(2-29)。根据机翼颤振时在某空速时由稳定转变为不稳定,空速V最小的一组解(V, ω)应为机翼的颤振速度和相应的颤振频率。整个求解过程如图2-4所示。
![](https://epubservercos.yuewen.com/3434AA/19391578301349906/epubprivate/OEBPS/Images/40359-00-025-1.jpg?sign=1739452230-08k3HeWFyoA2Jd3yubVZZWXOuMdVVXfN-0-f84ee4a67431f9ba8f1d66330b86051e)
图2-4 机翼颤振求解过程
2.3.3 机翼扭转发散求解
当求解机翼的发散特性时,仍然可利用上述方法,只需令式(2-29)中的圆频率ω→0,此时,复矩阵K(ω, V)将退化为实矩阵K(V),求解机翼扭转发散的特征方程将退化为
det[K(V)]=0 (2-30)
然而,在利用求解机翼颤振特性的程序退化后求解机翼发散时,当式(2-29)中令ω=0时,将会出现奇异从而导致计算无法进行。因而,在实际求解时可以令ω取接近于零的较小值,如0.001或更小,以避免求解奇异的现象。