一种空天飞行器机翼振动响应的快速仿真方法.pdf

本发明公开了一种空天飞行器机翼振动响应的快速仿真算法，属于航空航天技术领域，涉及一种气动弹性结构的振动仿真技术。通过对空天飞行器在高速飞行过程中机翼在流场中流固耦合运动问题进行建模，并且提出新的微分方程求解方法进行求解。定量地研究空天飞行器二元机翼在高速流场中的颤振、极限环等动力学特征。为空天飞行器机翼的设计建立精确可靠的数学模型，并且将一种新型的快速有效的时间离散法应用于非线性二元机翼的动力学模型研究中，进而研究机翼的各种动力学特征。该方法能够适用于强非线性二元机翼振动仿真，具有高效、高精度的特点。

1.一种空天飞行器机翼振动响应的快速仿真方法，其特征在于，包括如下步骤：
1)建立空天飞行器在流固耦合下的二元机翼颤振动力学模型；
2)通过对机翼颤振动力学模型进行求解建立时间离散法代数方程组；
3)建立紧凑型的时间离散法代数方程组；
4)求解上述紧凑型时间离散法代数方程组，得到机翼系统的振动响应曲线。
2.根据权利要求1所述的一种空天飞行器机翼振动响应的快速仿真方法，其特征在于，
步骤1)具体为：
结合二元机翼含有线性弹簧时的气动弹性模型，并考虑机翼系统在俯仰和沉浮两个自
由度的结构非线性，建立机翼系统方程：
$\stackrel{\·\·}{\mathrm{\ξ}}+x_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+2{\mathrm{\ζ}}_{\mathrm{\ξ}}\frac{\stackrel{\‾}{\mathrm{\ω}}}{U*}\stackrel{\·}{\mathrm{\ξ}}+{\left(\frac{\stackrel{\‾}{\mathrm{\ω}}}{U*}\right)}^{2}G\left(\mathrm{\ξ}\right)=-\frac{1}{\mathrm{\π}\mathrm{\μ}}{C}_L\left(\mathrm{\τ}\right),---\left(1\right)$
$\frac{x_{\mathrm{\α}}}{r_{\mathrm{\α}}^2}\stackrel{\·\·}{\mathrm{\ξ}}+\stackrel{\·\·}{\mathrm{\α}}+2{\mathrm{\ζ}}_{\mathrm{\α}}\frac{1}{U*}\stackrel{\·}{\mathrm{\α}}+{\left(\frac{1}{U*}\right)}^{2}M\left(\mathrm{\α}\right)=\frac{2}{{\mathrm{\π\μr}}_{\mathrm{\α}}^2}{C}_M\left(\mathrm{\τ}\right),---\left(2\right)$
其中，x_α是机翼气弹坐标系原点到质心的无量纲距离，ξ＝h/b是无量纲沉浮量，(·)表
示对无量纲时间τ的导数，τ＝Ut/b，t为时间，U^*为无量纲速度，定义为U^*＝U/(bω_α)，U为空
气来流速度，其中ω_ξ和ω_α分别是不耦合方程沉浮和俯仰自由度的固有频率，
ζ_ξ和ζ_α是阻尼比，r_α为绕弹性轴的转矩，α是俯仰角，h是偏转角，μ＝m/πρb²，m是机翼质量，ρ
为空气密度，b为机翼的半弦长；
M(α)和G(ξ)分别是俯仰和沉浮自由度的非线性项，表达式为：
M(α)＝α+βα³，G(ξ)＝ξ+γξ³， (3)
其中β和γ为非线性项系数；
C_L(τ)和C_M(τ)是线性气动力和气动力矩，表达式为：
$\begin{array}{c}{C}_L\left(\mathrm{\τ}\right)=\mathrm{\π}(\stackrel{\·\·}{\mathrm{\ξ}}-a_h\stackrel{\·\·}{\mathrm{\α}}+\stackrel{\·}{\mathrm{\α}})+2\mathrm{\π}\[\mathrm{\α}\left(0\right)+\stackrel{\·}{\mathrm{\ξ}}\left(0\right)+(\frac{1}{2}-a_h)\stackrel{\·}{\mathrm{\α}}\left(0\right)\]\mathrm{\φ}\left(\mathrm{\τ}\right)\\ +2\mathrm{\π}{\∫}_0^{\mathrm{\τ}}\mathrm{\φ}\left(\mathrm{\τ}-\mathrm{\σ}\right)\[\stackrel{\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)+\stackrel{\·\·}{\mathrm{\ξ}}\left(\mathrm{\σ}\right)+(\frac{1}{2}-a_h)\stackrel{\·\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)\]d\mathrm{\σ},\end{array}$
$\begin{array}{c}{C}_M\left(\mathrm{\τ}\right)=\mathrm{\π}(\frac{1}{2}+a_h)\[\mathrm{\α}\left(0\right)+\stackrel{\·}{\mathrm{\ξ}}\left(0\right)+(\frac{1}{2}-a_h)\stackrel{\·}{\mathrm{\α}}\left(0\right)\]\mathrm{\φ}\left(\mathrm{\τ}\right)+\frac{\mathrm{\π}}{2}(\stackrel{\·\·}{\mathrm{\ξ}}-a_h\stackrel{\·\·}{\mathrm{\α}})-\frac{\mathrm{\π}}{16}\stackrel{\·\·}{\mathrm{\α}}\\ -\left(\frac{1}{2}-a_h\right)\frac{\mathrm{\π}}{2}\stackrel{\·}{\mathrm{\α}}+\mathrm{\π}\left(\frac{1}{2}+a_h\right){\∫}_0^{\mathrm{\τ}}\mathrm{\φ}\left(\mathrm{\τ}-\mathrm{\σ}\right)\[\stackrel{\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)+\stackrel{\·\·}{\mathrm{\ξ}}\left(\mathrm{\σ}\right)+(\frac{1}{2}-a_h)\stackrel{\·\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)\]d\mathrm{\σ},\end{array}$
其中Wagner函数φ(τ)为ψ₁，ψ₂，ε₁，ε₂是Wagner常
数，a_h是机翼中轴线到气弹坐标原点的无量纲距离；
然后，通过引入一组积分变换式，将上式中C_L和C_M包含的积分项消除，从而将积分微分
方程转化为微分方程组，记为如下形式：
$c_0\stackrel{\·\·}{\mathrm{\ξ}}+c_1\stackrel{\·\·}{\mathrm{\α}}+c_2\stackrel{\·}{\mathrm{\ξ}}+c_3\stackrel{\·}{\mathrm{\α}}+c_4\mathrm{\ξ}+c_5\mathrm{\α}+c_6{w}_1+c_7{w}_2+c_8{w}_3+c_9{w}_4+c_{10}G\left(\mathrm{\ξ}\right)=f\left(\mathrm{\τ}\right),$
$d_0\stackrel{\·\·}{\mathrm{\ξ}}+d_1\stackrel{\·\·}{\mathrm{\α}}+d_2\stackrel{\·}{\mathrm{\ξ}}+d_3\stackrel{\·}{\mathrm{\α}}+d_4\mathrm{\ξ}+d_5\mathrm{\α}+d_6{w}_1+d_7{w}_2+d_8{w}_3+d_9{w}_4+d_{10}M\left(\mathrm{\α}\right)=g\left(\mathrm{\τ}\right),$
${\stackrel{\·}{w}}_1=\mathrm{\α}-{\∈}_1{w}_1,$
${\stackrel{\·}{w}}_2=\mathrm{\α}-{\∈}_2{w}_2,$
${\stackrel{\·}{w}}_3=\mathrm{\ξ}-{\∈}_1{w}_3,$
${\stackrel{\·}{w}}_4=\mathrm{\ξ}-{\∈}_2{w}_4.---\left(4\right)$
其中，c₀＝1+1/μ，c₁＝x_α-a_h/μ，
c₃＝[1+(1-2a_h)(1-ψ₁-ψ₂)]/μ，c₄＝2(ε₁ψ₁+ε₂ψ₂)/μ，
c₅＝2[1-ψ₁-ψ₂+(1/2-a_h)(ε₁ψ₁+ε₂ψ₂)]/μ，
c₆＝2ε₁ψ₁[1-ε₁(1/2-a_h)]/μ，c₇＝2ε₂ψ₂[1-ε₂(1/2-a_h)]/μ，
$c_8=-2{\mathrm{\ϵ}}_1^2{\mathrm{\ψ}}_1,c_9=-2{\mathrm{\ϵ}}_2^2{\mathrm{\ψ}}_2/\mathrm{\μ},c_{10}={(\stackrel{\‾}{\mathrm{\ω}}/U^{*})}^2,$
$d_0=x_{\mathrm{\α}}/r_{\mathrm{\α}}^2-a_h/{\mathrm{\μr}}_{\mathrm{\α}}^2,d_1=1+(1+8a_h^2)/\left(8{\mathrm{\μr}}_{\mathrm{\α}}^2\right),$
$d_2=-(1+2a_h)(1-{\mathrm{\ψ}}_1-{\mathrm{\ψ}}_2)/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$

$d_4=-(1+2a_h)({\mathrm{\ϵ}}_1{\mathrm{\ψ}}_1+{\mathrm{\ϵ}}_2{\mathrm{\ψ}}_2)/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$
$d_5=-(1+2a_h)(1-{\mathrm{\ψ}}_1-{\mathrm{\ψ}}_2)/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right)-(1-4a_h^2)({\mathrm{\ϵ}}_1{\mathrm{\ψ}}_1+{\mathrm{\ϵ}}_2{\mathrm{\ψ}}_2)/\left(2{\mathrm{\μr}}_{\mathrm{\α}}^2\right),$
$d_6=-(1+2a_h){\mathrm{\ψ}}_1{\mathrm{\ϵ}}_1\[1-{\mathrm{\ϵ}}_1(1/2-a_h)\]/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$
$d_7=-(1+2a_h){\mathrm{\ψ}}_2{\mathrm{\ϵ}}_2\[1-{\mathrm{\ϵ}}_2(1/2-a_h)\]/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$
$d_8=(1+2a_h){\mathrm{\ψ}}_1{\mathrm{\ϵ}}_1^2/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),d_9=(1+2a_h){\mathrm{\ψ}}_2{\mathrm{\ϵ}}_2^2/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),d_{10}={(1/U^{*})}^2,$
${\mathrm{\ω}}_1={\∫}_0^{\mathrm{\τ}}{e}^{-c_1(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\α}\left(\mathrm{\σ}\right)d\mathrm{\σ},{\mathrm{\ω}}_2={\∫}_0^{\mathrm{\τ}}{e}^{-c_2(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\α}\left(\mathrm{\σ}\right)d\mathrm{\σ},$
$\begin{array}{cc}{\mathrm{\ω}}_3={\∫}_0^{\mathrm{\τ}}{e}^{-c_1(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\ξ}\left(\mathrm{\σ}\right)d\mathrm{\σ}& {\mathrm{\ω}}_4={\∫}_0^{\mathrm{\τ}}{e}^{-c_2(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\ξ}\left(\mathrm{\σ}\right)d\mathrm{\σ},\end{array}$
其中，r_α为弹性坐标系的回转半径。
3.根据权利要求2所述的一种空天飞行器机翼振动响应的快速仿真方法，其特征在于，
步骤2)具体为：
将步骤1)得到的微分方程组(4)中的6个待求函数α(τ)，ξ(τ)，ω_i(τ)，i＝1,2,3,4，首先
假设成Fourier级数形式：
$\mathrm{\α}\left(\mathrm{\τ}\right)={\mathrm{\α}}_0+\underset{n=1}{\overset{N}{\Σ}}\[{\mathrm{\α}}_{2n-1}cos\left(n\mathrm{\ω}\mathrm{\τ}\right)+{\mathrm{\α}}_{2n}sin\left(n\mathrm{\ω}\mathrm{\τ}\right)\],$
$\mathrm{\ξ}\left(\mathrm{\τ}\right)={\mathrm{\ξ}}_0+\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[{\mathrm{\ξ}}_{2n-1}\cos \left(n\mathrm{\ω}\mathrm{\τ}\right)+{\mathrm{\ξ}}_{2n}\sin \left(n\mathrm{\ω}\mathrm{\τ}\right)\],$
$w_i\left(\mathrm{\τ}\right)=w_0^i+\underset{n=1}{\overset{N}{\Σ}}\[w_{2n-1}^{i}cos\left(n\mathrm{\ω}\mathrm{\τ}\right)+w_{2n}^{i}sin\left(n\mathrm{\ω}\mathrm{\τ}\right)\],i=1,2,3,4.---\left(5\right)$
其中，α_j，ξ_j，ω_j，j＝0,...,2N为相应的谐波系数；
将假设的近似解(5)代入微分方程组(4)中，得到残差函数：
$R_1=c_0\stackrel{\·\·}{\mathrm{\ξ}}+c_1\stackrel{\·\·}{\mathrm{\α}}+c_2\stackrel{\·}{\mathrm{\ξ}}+c_3\stackrel{\·}{\mathrm{\α}}+(c_4+c_{10})\mathrm{\ξ}+c_5\mathrm{\α}+c_6{w}_1+c_7{w}_2+c_8{w}_3+c_9{w}_4\≠0,$
$R_2=d_0\stackrel{\·\·}{\mathrm{\ξ}}+d_1\stackrel{\·\·}{\mathrm{\α}}+d_2\stackrel{\·}{\mathrm{\ξ}}+d_3\stackrel{\·}{\mathrm{\α}}+d_4\mathrm{\ξ}+d_5\mathrm{\α}+d_6{w}_1+d_7{w}_2+d_8{w}_3+d_9{w}_4+d_{10}{\mathrm{\α}}^3\≠0,$
$R_3={\stackrel{\·}{w}}_1+{\∈}_1{w}_1-\mathrm{\α}\≠0,$
$R_4={\stackrel{\·}{w}}_2+{\∈}_2{w}_2-\mathrm{\α}\≠0,$
$R_5={\stackrel{\·}{w}}_3+{\∈}_1{w}_3-\mathrm{\ξ}\≠0,$
$R_6={\stackrel{\·}{w}}_4+{\∈}_2{w}_4-\mathrm{\ξ}\≠0.---\left(6\right)$
其中，R_j表示
最后，迫使残差函数R_j在一个周期上的2N+1个等距时间点τ_i上为零即得到时间离散法
代数方程组，该机翼系统含有6×(2N+1)个方程，6×(2N+1)+1个未知数。
4.根据权利要求3所述的一种空天飞行器机翼振动响应的快速仿真方法，其特征在于，
步骤3)具体为：
对时间离散法代数方程组的每项进行时间离散处理，对α(τ)在2N+1个时间点τ_i离散，可
得：
$\mathrm{\α}\left({\mathrm{\τ}}_i\right)={\mathrm{\α}}_0+\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[{\mathrm{\α}}_{2n-1}cos\left({\mathrm{n\ω\τ}}_i\right)+{\mathrm{\α}}_{2n}sin\left({\mathrm{n\ω\τ}}_i\right)\].$
其中θ_i＝ωτ_i，将上式写为矩阵形式：

简单表示为：
${\tilde{Q}}_{\mathrm{\α}}\≡\left(\begin{array}{c}\mathrm{\α}\left({\mathrm{\τ}}_0\right)\\ \mathrm{\α}\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \mathrm{\α}\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)={\mathrm{FQ}}_{\mathrm{\α}},$
其中，
$F=\left(\begin{array}{cccccc}1& {\mathrm{cos\θ}}_0& {\mathrm{sin\θ}}_0& ...& \cos \left({\mathrm{N\θ}}_0\right)& \sin \left({\mathrm{N\θ}}_0\right)\\ 1& {\mathrm{cos\θ}}_1& {\mathrm{sin\θ}}_1& ...& \cos \left({\mathrm{N\θ}}_1\right)& \sin \left({\mathrm{N\θ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 1& {\mathrm{cos\θ}}_{2N}& {\mathrm{sin\θ}}_{2N}& ...& \cos \left({\mathrm{N\θ}}_{2N}\right)& \sin \left({\mathrm{N\θ}}_{2N}\right)\end{array}\right)$
Q_α＝[α₀,α₁,...,α_2N]^T
类似的其中Q_ξ分别为与ξ有关的配点
和谐波系数，和分别为与ω_i有关的配点和谐波系数；
然后，对在2N+1个时间点进行离散可得：
$\stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_i\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[-{\mathrm{n\ω\α}}_{2n-1}sin\left({\mathrm{n\ω\τ}}_i\right)+{\mathrm{n\ω\α}}_{2n}cos\left({\mathrm{n\ω\τ}}_i\right)\]$
将上式写为矩阵形式：

令矩阵
因此有：
${\tilde{Q}}_{\stackrel{\·}{\mathrm{\α}}}={\mathrm{\ωFAQ}}_{\mathrm{\α}},---\left(8\right)$
其中
${\tilde{Q}}_{\stackrel{\·}{\mathrm{\α}}}\≡\left(\begin{array}{c}\stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_0\right)\\ \stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)$
同理有：其中和分别是
关于ξ和ω_j导数的配点；
然后，对在2N+1个时间点离散，有：
$\stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_i\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[-n^2{\mathrm{\ω}}^2{\mathrm{\α}}_{2n-1}cos\left({\mathrm{n\ω\τ}}_i\right)-n^2{\mathrm{\ω}}^2{\mathrm{\α}}_{2n}sin\left({\mathrm{n\ω\τ}}_i\right)\]---\left(9\right)$
将上式(9)写为矩阵形式：

上式中的方阵用FA²表示为
${\mathrm{FA}}^2=\left(\begin{array}{cccccc}0& -{\mathrm{cos\θ}}_0& -{\mathrm{sin\θ}}_0& ...& -N^2\cos \left({\mathrm{N\θ}}_0\right)& -N^2\sin \left({\mathrm{N\θ}}_0\right)\\ 0& -{\mathrm{cos\θ}}_1& -{\mathrm{sin\θ}}_1& ...& -N^2\cos \left({\mathrm{N\θ}}_1\right)& -N^2\sin \left({\mathrm{N\θ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 0& -{\mathrm{cos\θ}}_{2N}& -{\mathrm{sin\θ}}_{2N}& ...& -N^2\cos \left({\mathrm{N\θ}}_{2N}\right)& -N^2\sin \left({\mathrm{N\θ}}_{2N}\right)\end{array}\right)$
因此，
${\tilde{Q}}_{\stackrel{\·\·}{\mathrm{\α}}}\≡\left(\begin{array}{c}\stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_0\right)\\ \stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)={\mathrm{\ω}}^2{\mathrm{FA}}^2{Q}_{\mathrm{\α}}$
同理有其中和分别是关于ξ和ω_j二阶导数的配点；
由此时间离散法代数方程组经变换后得到如下形式：
$(c_0{\mathrm{\ω}}^2{D}^2+c_2\mathrm{\ω}D+c_{4}I+c_{10}I){\tilde{Q}}_{\mathrm{\ξ}}+(c_1{\mathrm{\ω}}^2{D}^2+c_3\mathrm{\ω}D+c_{5}I){\tilde{Q}}_{\mathrm{\α}}+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_{i+5}{\tilde{Q}}_{w_i}=0,$
$(d_0{\mathrm{\ω}}^2{D}^2+d_2\mathrm{\ω}D+d_{4}I){\tilde{Q}}_{\mathrm{\ξ}}+(d_1{\mathrm{\ω}}^2{D}^2+d_3\mathrm{\ω}D+d_{5}I){\tilde{Q}}_{\mathrm{\α}}+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{d}_{i+5}{\tilde{Q}}_{w_i}+d_{10}{\tilde{M}}_{\mathrm{\α}}=0,$
$(\mathrm{\ω}D+{\∈}_{1}I){\tilde{Q}}_{{\mathrm{\ω}}_1}-{\tilde{Q}}_{\mathrm{\α}}=0,$
$(\mathrm{\ω}D+{\∈}_{2}I){\tilde{Q}}_{{\mathrm{\ω}}_2}-{\tilde{Q}}_{\mathrm{\α}}=0,$
$(\mathrm{\ω}D+{\∈}_{1}I){\tilde{Q}}_{{\mathrm{\ω}}_3}-{\tilde{Q}}_{\mathrm{\ξ}}=0,$
$(\mathrm{\ω}D+{\∈}_{2}I){\tilde{Q}}_{{\mathrm{\ω}}_4}-{\tilde{Q}}_{\mathrm{\ξ}}=0,---\left(10\right)$
其中D＝FAF^-1，且
${\tilde{M}}_{\mathrm{\α}}=\left(\begin{array}{c}{\mathrm{\α}}^3\left({\mathrm{\τ}}_0\right)\\ {\mathrm{\α}}^3\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\α}}^3\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)$
将和用表示，然后将他们代入到非线性方程中得：
$(A_{2\mathrm{\α}}-B_{2\mathrm{\ξ}}{B}_{1\mathrm{\ξ}}^{-1}{A}_{1\mathrm{\α}}){\tilde{Q}}_{\mathrm{\α}}+d_{10}{\tilde{M}}_{\mathrm{\α}}=0,---\left(11\right)$
其中A_1α，B_1ξ，A_2α，B_2ξ分别为：
A_1α＝c₁ω²D²+c₃ωD+c₅I+c₆(ωD+∈₁I)^-1+c₇(ωD+∈₂I)^-1
B_1ξ＝c₀ω²D²+c₂ωD+(c₄+c₁₀)I+c₈(ωD+∈₁I)^-1+c₉(ωD+∈₂I)^-1
A_2α＝d₁ω²D²+d₃ωD+d₅I+d₆(ωD+∈₁I)^-1+d₇(ωD+∈₂I)^-1
B_2ξ＝d₀ω²D²+d₂ωD+d₄I+d₈(ωD+∈₁I)^-1+d₉(ωD+∈₂I)^-1.
式(11)为以时域变量为变量的紧凑型时间离散法代数方程组，使用时-频转换关系式
和将上式转换为以频率变量为变量的紧凑型时间离散法代
数方程组(12)：
$(A_2-B_2{B}_1^{-1}{A}_1)Q_{\mathrm{\α}}+d_{10}{F}^{-1}{\left({\mathrm{FQ}}_{\mathrm{\α}}\right)}^3=0,---\left(12\right)$
式(12)含有2N+2个未知数，其中，A₂＝A_2α，A₁＝A_1α，B₂＝B_2ξ，B₁＝B_1ξ。
5.根据权利要求1所述的一种空天飞行器机翼振动响应的快速仿真方法，其特征在于，
步骤4)中使用标量同伦算法求解上述紧凑型时间离散法代数方程组。

一种空天飞行器机翼振动响应的快速仿真方法

技术领域

本发明属于航空航天技术领域，涉及一种机翼翼型气动弹性振动仿真技术，具体
涉及一种空天飞行器机翼振动响应的快速仿真方法。

背景技术

空天飞行器是一种可重复使用的新一代天地往返系统，对于降低空间运输成本具
有重要的意义。其特点之一是能在外太空和大气层之间自由起落和飞行。要发展空天飞行
器,就要解决以下关键技术:即兼具超轻质量、高强韧、耐热和抗冲击性能的空天器结构技
术,高超声速技术,高机动飞行技术,长距离空天飞行技术,高隐形技术,精准打击及可靠性
技术。在这些关键技术中，气动弹性问题是设计、制造空天飞行器的基础问题，而与该问题
密切相关的飞行器机翼的设计则是当前的研究热点和技术难点。

空天飞行器的机翼在高速飞行过程中的振动形式直接影响机翼结构体的强度。一
般而言，机翼可能存在的振动形式包括：稳定态、极限环振动、周期振动、混沌运动和发散。
机翼的设计尽量使其在整个飞行过程中处于稳定状态。这样对机翼结构体的强度损伤是最
小的。极限环振动是指机翼做类简谐振动。而周期振动是指机翼的振动稍微复杂但是振动
幅值具有重复性。不同翼型参数和飞行速度下，机翼颤振边界和振动形式是不同的。因此有
必要对机翼颤振系统进行分析。原始的机翼模型富含积分项。因为积分项的存在，不管是对
于直接的数值求解算法还是摄动法或者加权残余法来说，求解起来都需要事先进行大量的
符号处理，消耗大量的计算资源。现有的方法如谐波平衡法，摄动法等往往不能快速、全面
的分析系统的响应。这很大程度上影响了机翼结构的参数设计。

发明内容

本发明的目的在于提供一种空天飞行器机翼振动响应的快速仿真方法，以克服上
述现有技术存在的缺陷，本发明基于时间离散法，通过在时间域的离散获得飞行器机翼气
弹系统的代数方程，避免了一般的非线性系统分析方法中大量的符号运算，同时简化了系
统方程，使得求解效率得到较大提高。

为达到上述目的，本发明采用如下技术方案：

一种空天飞行器机翼振动响应的快速仿真方法，包括如下步骤：

1)建立空天飞行器在流固耦合下的二元机翼颤振动力学模型；

2)通过对机翼颤振动力学模型进行求解建立时间离散法代数方程组；

3)建立紧凑型的时间离散法代数方程组；

4)求解上述紧凑型时间离散法代数方程组，得到机翼系统的振动响应曲线。

进一步地，步骤1)具体为：

结合二元机翼含有线性弹簧时的气动弹性模型，并考虑机翼系统在俯仰和沉浮两
个自由度的结构非线性，建立机翼系统方程：

其中，x_α是机翼气弹坐标系原点到质心的无量纲距离，ξ＝h/b是无量纲沉浮量，
(·)表示对无量纲时间τ的导数，τ＝Ut/b，t为时间，U^*为无量纲速度，定义为U^*＝U/(bω_α)，
U为空气来流速度，其中ω_ξ和ω_α分别是不耦合方程沉浮和俯仰自由度的固有
频率，ζ_ξ和ζ_α是阻尼比，r_α为绕弹性轴的转矩，α是俯仰角，h是偏转角，μ＝m/πρb²，m是机翼质
量，ρ为空气密度，b为机翼的半弦长；

M(α)和G(ξ)分别是俯仰和沉浮自由度的非线性项，表达式为：

M(α)＝α+βα³，G(ξ)＝ξ+γξ³， (3)

其中β和γ为非线性项系数；

C_L(τ)和C_M(τ)是线性气动力和气动力矩，表达式为：

其中Wagner函数φ(τ)为ψ₁，ψ₂，ε₁，ε₂是
Wagner常数，a_h是机翼中轴线到气弹坐标原点的无量纲距离；

然后，通过引入一组积分变换式，将上式中C_L和C_M包含的积分项消除，从而将积分
微分方程转化为微分方程组，记为如下形式：

其中，c₀＝1+1/μ，c₁＝x_α-a_h/μ，

c₃＝[1+(1-2a_h)(1-ψ₁-ψ₂)]/μ，c₄＝2(ε₁ψ₁+ε₂ψ₂)/μ，

c₅＝2[1-ψ₁-ψ₂+(1/2-a_h)(ε₁ψ₁+ε₂ψ₂)]/μ，

c₆＝2ε₁ψ₁[1-ε₁(1/2-a_h)]/μ，c₇＝2ε₂ψ₂[1-ε₂(1/2-a_h)]/μ，

其中，r_α为弹性坐标系的回转半径。

进一步地，步骤2)具体为：

将步骤1)得到的微分方程组(4)中的6个待求函数α(τ)，ξ(τ)，ω_i(τ)，i＝1,2,3,
4，首先假设成Fourier级数形式：

其中，α_j，ξ_j，ω_j，j＝0,...,2N为相应的谐波系数；

将假设的近似解(5)代入微分方程组(4)中，得到残差函数：

其中，R_j表示

最后，迫使残差函数R_j在一个周期上的2N+1个等距时间点τ_i上为零即得到时间离
散法代数方程组，该机翼系统含有6×(2N+1)个方程，6×(2N+1)+1个未知数。

进一步地，步骤3)具体为：

对时间离散法代数方程组的每项进行时间离散处理，对α(τ)在2N+1个时间点τ_i离
散，可得：

其中θ_i＝ωτ_i，将上式写为矩阵形式：

简单表示为：

其中，

Q_α＝[α₀,α₁,...,α_2N]^T

类似的其中Q_ξ分别为与ξ有关
的配点和谐波系数，和分别为与ω_i有关的配点和谐波系数；

然后，对在2N+1个时间点进行离散可得：

将上式写为矩阵形式：

令矩阵

因此有：

其中

同理有：其中和分
别是关于ξ和ω_j导数的配点；

然后，对在2N+1个时间点离散，有：

将上式(9)写为矩阵形式：

上式中的方阵用FA²表示为

因此，

同理有其中和分别是关于ξ和ω_j二阶导数的配点；

由此时间离散法代数方程组经变换后得到如下形式：

其中D＝FAF^-1，且

将和用表示，然后将他们代入到非线性方程中得：

其中A_1α，B_1ξ，A_2α，B_2ξ分别为：

A_1α＝c₁ω²D²+c₃ωD+c₅I+c₆(ωD+∈₁I)^-1+c₇(ωD+∈₂I)^-1

B_1ξ＝c₀ω²D²+c₂ωD+(c₄+c₁₀)I+c₈(ωD+∈₁I)^-1+c₉(ωD+∈₂I)^-1

A_2α＝d₁ω²D²+d₃ωD+d₅I+d₆(ωD+∈₁I)^-1+d₇(ωD+∈₂I)^-1

B_2ξ＝d₀ω²D²+d₂ωD+d₄I+d₈(ωD+∈₁I)^-1+d₉(ωD+∈₂I)^-1.

式(11)为以时域变量为变量的紧凑型时间离散法代数方程组，使用时-频转换关
系式和将上式转换为以频率变量为变量的紧凑型时间离散
法代数方程组(12)：

式(12)含有2N+2个未知数，其中，A₂＝A_2α，A₁＝A_1α，B₂＝B_2ξ，B₁＝B_1ξ。

进一步地，步骤4)中使用标量同伦算法求解上述紧凑型时间离散法代数方程组。

与现有技术相比，本发明具有以下有益的技术效果：

本发明方法基于时间离散法，通过在时间域的离散获得飞行器机翼气弹系统的代
数方程，避免了一般的非线性系统分析方法中大量的符号运算，同时简化了系统方程，使得
求解效率得到较大提高，通过将时间离散方案与参数扫描法配合，成功避免了一般分析方
法中可能出现的非物理解，提高了机翼振动仿真的可靠性，能够对机翼系统的稳态周期响
应进行直接估计，而无需对瞬态振动进行模拟，针对性更强，节省了计算资源，同时也便于
分析系统在变参数下动力学响应的变化情况。

附图说明

图1为本发明方法的流程图；

图2为二元机翼模型示意图；

图3为使用低阶时间离散法计算的机翼振动频率和飞行速度响应曲线，其中a)为
使用9阶时域配点法得到的结果，b)为使用11阶时域配点法得到的结果；

图4为使用高阶时间离散法计算的机翼振动频率和飞行速度响应曲线，其中a)为
使用15阶时域配点法得到的结果，b)为使用21阶时域配点法得到的结果。

具体实施方式

下面对本发明作进一步详细描述：

一种基于时间离散方法的空天飞行器机翼振动响应的快速仿真方法，该方法的特
点在于,采用积分变换方法,将空天飞行器机翼翼型颤振积分-微分方程转化为纯微分方
程。进一步，选取含有未知系数的傅里叶级数作为颤振响应近似解，将该近似解带入颤振微
分方程，然后使用时间离散方法将微分方程离散为以傅里叶系数和振动频率作为待求变量
的代数方程组。最后，使用标量同伦法求解待求系数，从而确定空天飞行器机翼颤振的周期
解。该方法给出的解释颤振系统的半解析周期解，该方法具有很高的计算效率和精度。

具体包括以下步骤：

1)空天飞行器机翼颤振动力学模型建立：根据拉格朗日方程，考虑结构非线性和
活塞理论，建立空天飞行器在流固耦合下的二元机翼振动原始积分微分方程，并进行无量
纲处理。引入一组积分变换式，将积分微分方程中的积分项消除，从而得纯微分方程形式的
动力学模型。

2)使用时间离散法对机翼颤振动力学方程进行离散化处理：将机翼模型中的待求
函数，即机翼俯仰角和沉浮位移，以及由积分变换式引入的变量假设为傅里叶级数的形式，
并将傅里叶级数的待求系数组装成向量。将假设的近似解代入上述微分方程组中，得到残
差函数。然后，迫使残差函数在一个周期上有限个点上为零得到时间离散法的代数方程组。

3)建立紧凑型的时间离散法代数方程组：由于步骤2中建立的方程组维数太大，不
利于计算，因此对其进行等价转换，以得到紧凑形式的代数方程。利用假设函数为傅里叶级
数这一特点，建立离散点与谐波系数之间的矩阵关系，然后通过矩阵变换，建立时间一阶导
数和二阶导数的离散点与谐波系数之间的关系。基于矩阵变换关系，对残差函数构成的非
线性方程组进行逐项处理，并将所得方程组中的线性部分带入到非线性部分中，就能够得
到形式简单，维数较低的紧凑型时间离散法代数方程组。

4)使用标量同伦方法求解上述紧凑型的时间离散法代数方程组：时间离散法对初
值非常敏感，随意给定初值一般无法得到真实解，因此，用初值不敏感的标量同伦法为提供
合理初值获得半解析周期解。

为了更好地说明本发明的目的和优点，下面结合附图和实例对本发明内容做进一
步说明：

1.建立二元机翼振动模型的数学模型：

图2是含有俯仰和沉浮两自由度的二元机翼振动模型。沉浮用h表示，向下为正方
向。关于弹性轴的俯仰用α表示，往上仰为正。弹性轴距翼型中心的距离为a_hb，质心离弹性
轴的距离为x_ab；两个距离的正方向指向机翼后缘。

考虑俯仰和沉浮的结构非线性的无量纲形式系统方程为:

其中x_α是机翼气弹坐标系原点到质心的无量纲距离，ξ＝h/b是无量纲沉浮量，
(·)表示对无量纲时间τ的导数，其中τ＝Ut/b，t为时间，U^*为无量纲速度，定义为U^*＝U/(b
ω_α)，U为空气来流速度。其中ω_ξ和ω_α分别是不耦合方程沉浮和俯仰自由度的
固有频率。ζ_ξ和ζ_α是阻尼比，r_α为绕弹性轴的转矩，α是俯仰角，h是偏转角，μ＝m/πρb²，m是机
翼质量，ρ为空气密度，b为机翼的半弦长。M(α)和G(ξ)分别是俯仰和沉浮自由度的非线性
项，他们的具体表达式为：

M(α)＝α+βα³，G(ξ)＝ξ+γξ³， (3)

其中β和γ为非线性项系数。C_L(τ)和C_M(τ)是线性气动力和气动力矩：

其中Wagner函数φ(τ)为ψ₁，ψ₂，ε₁，ε₂是
Wagner常数，a_h是机翼中轴线到气弹坐标原点的无量纲距离。为了便于模型求解，我们引入
下面四个积分变换关系式，原方程组可等价变换为如下形式：

其中：

c₀＝1+1/μ，c₁＝x_α-a_h/μ，

c₃＝[1+(1-2a_h)(1-ψ₁-ψ₂)]/μ，c₄＝2(ε₁ψ₁+ε₂ψ₂)/μ，

c₅＝2[1-ψ₁-ψ₂+(1/2-a_h)(ε₁ψ₁+ε₂ψ₂)]/μ，

c₆＝2ε₁ψ₁[1-ε₁(1/2-a_h)]/μ，c₇＝2ε₂ψ₂[1-ε₂(1/2-a_h)]/μ，

其中，r_α为弹性坐标系的回转半径。

2.二元机翼模型的时间离散法求解方案：

将方程组(4)中的6个待求函数α(τ)，ξ(τ)，ω_i(τ)，i＝1,2,3,4首先假设成
Fourier级数形式。

将假设的近似解(5)代入数学方程(4)中，得到残差函数：

其中R_j表示现在，迫使残差函
数R_j在一个周期上的2N+1个等距时间点τ_i上为零就可以得到时间离散法代数方程组，该系
统含有6×(2N+1)个方程6×(2N+1)+1个未知数。

3.建立紧凑型的时间离散法代数方程组：

对时间离散法代数方程组的每项进行时间离散处理。对α(τ)在2N+1个时间点τ_i离
散，可得：

其中θ_i＝ωτ_i，上式可写为矩阵形式：

也可简单表示为，

其中

Q_α＝[α₀,α₁,...,α_2N]^T

类似的有其中Q_ξ分别为与ξ有
关的配点和谐波系数，和分别为与ω_i有关的配点和谐波系数。

然后，对在2N+1个时间点进行离散可得：

上式易写为矩阵形式：

令矩阵

因此有：

其中

相同的道理，其中和
分别是关于ξ和ω_j导数的配点。

然后，对在2N+1个时间点离散，有：

方程(9)的矩阵形式为：

上式中的方阵可用FA2表示为

因此，

类似的其中和分别是关于ξ和ω_j二阶导数的配点。

上面，我们对方程进行了逐项处理，得到了时间离散后的矩阵形式。时间离散法代
数方程经变换后得到如下形式：

其中D＝FAF^-1，且

由于只有第二个方程是非线性的，其他五个方程均为线性方程，因此方程组(10)
可进一步简化。将和用表示，然后将他们代入到非线性方程中得：

其中A_1α，B_1ξ，A_2α，B_2ξ为：

A_1α＝c₁ω²D²+c₃ωD+c₅I+c₆(ωD+∈₁I)^-1+c₇(ωD+∈₂I)^-1

B_1ξ＝c₀ω²D²+c₂ωD+(c₄+c₁₀)I+c₈(ωD+∈₁I)^-1+c₉(ωD+∈₂I)^-1

A_2α＝d₁ω²D²+d₃ωD+d₅I+d₆(ωD+∈₁I)^-1+d₇(ωD+∈₂I)^-1

B_2ξ＝d₀ω²D²+d₂ωD+d₄I+d₈(ωD+∈₁I)^-1+d₉(ωD+∈₂I)^-1.

系统(11)为以时域变量为变量的紧凑型时间离散法代数方程组。

使用时-频转换关系式和系统(11)变换为：

方程(12)是以频率变量为变量的紧凑型时间离散法代数方程组，含有2N+2个未知
数，其中，A₂＝A_2α，A₁＝A_1α，B₂＝B_2ξ，B₁＝B_1ξ。

4.使用标量同伦方法求解上述紧凑型的时间离散法代数方程组，并结合参数扫描
法求解上述紧凑型的时间离散法代数方程组：

时间离散法对初值非常敏感，随意给定初值一般无法得到真实解，因此，用初值不
敏感的标量同伦法为提供合理初值获得半解析周期解。为了进行系统参数设计，需要求取
响应的幅频曲线。这里使用参数扫描法为时间离散法提供合理初值以获得幅频响应曲线。
图3和图4是时间离散法和RK4计算得出的关于振动基频比飞行速度的响应曲线。图3(a)显
示，对于正向扫描，在之前，TDC9和RK4的结果吻合，超过以后TDC9开始与RK4
分开。并且，TDC9的响应曲线是连续变化的，也就是说TDC9不能预测出二元机翼振动系统的
次级分岔。对于逆向扫描，可以看出TDC9与RK4在整个下支曲线一致。使用RK4预测的分岔点
为1.84，超过该分岔点后RK4的响应曲线从下支跳跃到上支。对于TDC9来说，在分岔点之前
TDC9能够得到整个下支响应曲线。但是超越分岔点后，TDC9响应曲线没有跳跃到上支响应
曲线，而是蜿蜒曲折的在下支附近得到其他一些曲线。这是因为到达分岔点后参数扫描法
无法提供合理的初值了，扫描法失效。由于TDC法对初值的极其敏感性，因此TDC9在分岔点
后的结果对应于非物理解。

图3(b)是TDC11计算出的频率--速度响应曲线图。正向扫描时，TDC11的结果和标
准解直到2.27都是很接近的，这要好于TDC9的结果。继续增大速度TDC11的结果出现偏离，
这意味着非物理解出现了。TDC11跟TDC9一样，都不能检测出次级分岔点。逆向扫描时，
TDC11可以得到整个响应曲线的下支，并且得到了次级Hopf分岔值比标准值1.84略小。由于
分岔点使扫描法失效，因此超过分岔点后TDC11给出的是非物理解。

图4(a)表明逆向扫描时TDC15可以精确的描绘出整个下支响应曲线，但是正向扫
描时不能给出整个上支。TDC21的结果见图4(b)。如图所示，正向扫描得到了整个上支曲线，
并且预测的分岔点为2.38接近标准值2.34。超越分岔点后，不可避免的出现了非物理解。逆
向扫描的结果与标准解吻合甚好，整个下支曲线高度一致，分岔点也精确的获得。超越分岔
点后，非物理解分支出现，直到才从非物理解分支跳跃到有物理意义的上
支曲线。

总之，使用时间离散法与参数扫描法相结合的办法可以得到非线性动力学系统的
各种响应曲线。并且，随着谐波个数的增加，计算精度越高。