[发明专利]一种空天飞行器机翼振动响应的快速仿真方法

权利要求书

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_α为弹性坐标系的回转半径。