[发明专利]基于气动非线性与耦合的高速飞行器随控优化方法

权利要求书

1.一种基于气动非线性与耦合的高速飞行器随控优化方法，其特征在于包括以下步骤：

步骤一、将气动非线性评价指标定义为俯仰通道的气动非线性度，即

DONLmz(α)=|Cm(α)-Cm_L(α)||Cm(α)-Cm_L(α)|+|Cm_L(α)|---(1)]]>

=|mz(α)-mz(0)-mz′(0)α||mz(α)-mz(0)-mz′(0)α|+|mz(0)+mz′(0)α|]]>

式中，C_m(α)是飞行器的稳定俯仰力矩系数；C_{m_L}(α)为C_m(α)中的线性部分；

1）C_m(α)计算公式；

C_m(α)＝C_弹身+K_αα·C_翼面＝C_{m_zhi1}+C_{m_zhi2}+C_{m_zhu}+K_αα·C_{m_w}（2）

其中，

Cm_zhi1=4sinαcosαsin2θ1(L1tanθ1+L2tanθ2)2Lref(12xcL12-13L13)]]>

Cm_zhi2=4sinαcosαsinθ2cosθ2(L1tanθ1+L2tanθ2)2Lref[L1L2tanθ1(xc-L1-12L2)+L22tanθ2(12(xc-L1)-13L2)]]]>

Cm_zhu=8sin2αL3[xc-L3/2-(L1+L2)3π(L1tanθ1+L2tanθ2)Lref]]>

式中：L₁为尖锥长度，L₂为截锥长度，L₃为柱段长度，θ₁为尖锥锥角，θ₂为截锥锥角，α为攻角，x_c为质心位置，L_ref为参考长度；

当翼的平面形状为双后掠翼时，其稳定俯仰力矩系数C_{m_w}为：

Cm_w=(γ+1)α2SrefLref1+[4(γ+1)M∞α]2{cotλI[b0I2(xc-xs)-23b0I3]-]]>

23cotλII(b0II-b1)3+[cotλII(xc-xs)-b0I(cotλI+cotλII)](b0II-b1)2]]>

+2b0IcotλI(b0II-b1)(xc-xs-b0I)]]>

+[(b0II-b1)cotλII+b0IcotλI][2(xc-b0I-b0II-xs)b1+b12]}]]>

式中：γ为绝热系数，S_ref为参考面积，L_ref为参考长度，λ_I为内翼前缘后掠角，λ_II为外翼前缘后掠角，b_0I为内翼根弦长度，b_0II为外翼根弦长度，b₁为翼面梢弦，x_s为外露翼安装位置距离弹身顶点的距离；

当翼的平面形状为单后掠翼时，其稳定俯仰力矩系数C_{m_w}为：

Cm_w=Mzq∞SrefLref=(γ+1)α2SrefLrefcotλI1+[4(γ+1)M∞α]2[(b0I-b1)2(xc-xs)]]>

-23(b0I-b1)3+2b1(b0I-b1)(xc-xs-b0I+b1)-(b0I-b1)b12]]]>

式中：γ为绝热系数，S_ref为参考面积，L_ref为参考长度，λ_I为内翼前缘后掠角，b_0I为内翼根弦长度，b₁为翼面梢弦，x_s为外露翼安装位置距离弹身顶点的距离，K_αα是翼身干扰系数；

2）C_{m_L}(α)计算公式；

C_{m_L}(α)＝C_{m_L_zhi1}+C_{m_L_zhi2}+C_{m_wb_L}（3）

其中，

Cm_zhi1_L=4αsin2θ1(L1tanθ1+L2tanθ2)2Lref(12xcL12-13L13)]]>

Cm_zhi2_L=4αsinθ2cosθ2(L1tanθ1+L2tanθ2)2Lref[12L1L2tanθ1(2xc-2L1-L2)]]>

+L22tanθ2(12(xc-L1)-13L2)]]>

Cm_wb_L=-0.035αξ(1-ηt2)Lref(xc-Lsh-Rmaxηtcotθ3+23Rmaxcotθ3(1-ηt4)1-ηt3)]]>

式中：L₁为尖锥长度，L₂为截锥长度，θ₁为尖锥锥角，θ₂为截锥锥角，θ₃为尾部收缩角，α为攻角，x_c为质心位置，L_ref为参考长度，η_t为尾部收缩比，R_max为圆柱部分半径；

步骤二、定义气动耦合度评价指标之方向舵偏转对滚转通道的耦合度

mxδymxδx=p2kTScwyr(Kδ0)cwrScwr(Kδ0)WSL(Kδ0)WrSWrσ(ηk+1)(1-D‾)2ηk+1-2D‾[D‾+(1-D‾)f]---(4)]]>

式中，飞行器为单垂尾布局时，p＝1；飞行器为双垂尾布局时，p是关于(α,Ma,h)的函数；k_T为气流阻滞系数，S_cw为垂尾面积，S为飞行器特征面积，L为飞行器特征长度，y_r为方向舵面心到弹身纵轴的距离，S_cwr为方向舵面积，S_Wr为垂尾面积，(K_δ0)_cwr为方向舵与翼身之间的干扰系数，(K_δ0)_W为翼身干扰系数，(K_δ0)_Wr为垂尾与翼身之间的干扰系数，η_k为根梢比，为径展比，f为外露翼根弦剖面到压心的距离与半展长之彼岸，σ为弹身相对直径变化对滚动力矩产生影响所引进的修正系数；

定义气动耦合度评价指标之方向舵偏转对俯仰通道的耦合度

mzδymzδz=qa0(1+a1h)(1+a2α)(1+a3Ma)(1+a4ηkew)(1+a5λkew)(1+a6cosχcw)---(5)]]>

(1+a7λkcwr)(1+a8ycwr)(1+a9Lcw)(1+a10D‾cw)(1+a11ScwrS)(1+a12ScwS)]]>

式中，飞行器为单垂尾布局时，q＝1；飞行器为双垂尾布局时，q是关于(α,Ma,h)的函数；S_cw、η_kcw、λ_kcw、χ_cw分别为垂直尾翼的面积，外露垂尾根梢比，外露垂尾展弦比和后掠角；S_cwr、λ_kcwr、L_cw分别表示方向舵的面积，方向舵展弦比，飞行器垂尾径展比和出位外露根弦中点到飞行器头部的距离；

定义气动耦合度评价指标之差动舵偏转对偏航通道的耦合度

myδxmyδy=mb0(1+b1h)(1+b2α)(1+b3Ma)(1+b4ηkw)(1+b5λkw)(1+b6cosχw)---(6)]]>

(1+b7λkwr)(1+b8ywr)(1+b9Lw)(1+b10D‾cw)(1+b11SwrS)(1+b12SwS)]]>

式中，飞行器为单后掠布局时，m＝1；飞行器为双后掠布局时，m是关于(α,Ma,h)的函数；S_w翼面面积，η_kw外露垂尾根梢比，λ_kw为翼的展弦比，χ_w为垂尾后掠角，S_wr为差动舵面积，λ_kwr为差动舵展弦比，为飞行器弹翼径展比，L_w为外露翼根弦中点到飞行器头部的距离，y_wr为差动舵面心到弹身纵轴的距离；

定义气动耦合度评价指标之差动舵偏转对俯仰通道的耦合度

mzδxmzδz=nc0(1+c1h)(1+c2α)(1+c3Ma)(1+c4ηkw)(1+c5λkw)(1+c6cosχw)---(7)]]>

(1+c7λkwr)(1+c8ywr)(1+c9Lw)(1+c10D‾w)(1+c11SwrS)(1+c12SwS)]]>

式中，飞行器为单后掠布局时，n＝1；飞行器为双后掠布局时，n是关于(α,Ma,h)的函数；S_w翼面面积，η_kw外露垂尾根梢比，λ_kw为翼的展弦比，χ_w为垂尾后掠角，S_wr为差动舵面积，λ_kwr为差动舵展弦比，为飞行器弹翼径展比，L_w为外露翼根弦中点到飞行器头部的距离，y_wr为差动舵面心到弹身纵轴的距离；

式（5）、式（6）和式（7）中的系数a₀～a₁₂、b₀～b₁₂、c₀～c₁₂是基于大量CFD计算数据并利用多元非线性最小二乘拟合方法计算得到的拟合值；

a0=0.08263,a1=6.138e-05,a2=-0.08484,a3=-0.03714,a4=0.02908a5=0.04677,a6=0.05188,a7=0.01160,a8=0.05631,a9=0.00876a10=0.06555,a11=0.59681,a12=0.15504]]>

b0=1.01895,b1=4.422e-05,b2=6.5758,b3=-0.03879,b4=-0.03635b5=-0.24995,b6=-1.19032,b7=-0.05226,b8=-0.25926,b9=-0.10399b10=-0.74664,b11=-2.13691,b12=-0.06801]]>

c0=-0.0269,c1=-1.2473e-05,c2=-0.6602,c3=-0.0121,c4=0.0124c5=0.0737,c6=0.4357,c7=0.0178,c8=0.0772,c9=0.0343c10=0.2224,c11=0.6530,c12=0.0233]]>

步骤三、飞行器在全弹道上任意特征点i处都存在一组气动非线性与耦合度评价指标：

{(DONLmz(α))i,(mxδy/mzδz)i,(mzδy/mzδz)i,(myδx/myδy)i,(mzδx/mzδz)i}]]>

其中，i＝1,2...N，N为全弹道上特征点的总数；

由于在全弹道多特征点处进行随控优化的结果需要满足全弹道上任意特征点处的非线性与耦合度的评价指标要求，因此需要构造以下函数：

f1(X)=max{|DONLmz(α)|i},]]>i＝1,2...N

f2(X)=max{|mxδy/mxδx|i},]]>i＝1,2...N

f3(X)=max{|mzδy/mzδz|i},]]>i＝1,2...N

f4(X)=max{|myδx/myδy|i},]]>i＝1,2...N

f5(X)=max{|mzδx/mzδz|i},]]>i＝1,2...N

式中，X是由飞行器的外形结构参数与飞行状态参数构成的向量(θ₁,θ₂,...h,Ma,α)^T；

根据高速飞行器随控优化的目的以及飞行器的气动非线性度、耦合度评价指标的定义，高速飞行器全弹道多特征点随控优化的目标函数为：

minf1(X)=max{|DONLmz(α)|i}]]>

minf2(X)=max{|mxδy/mxδx|i}]]>

minf3(X)=max{|mzδy/mzδz|i}]]>

minf4(X)=max{|myδx/myδy|i}]]>

minf5(X)=max{|mzδx/mzδz|i}]]>

其中，i＝1,2...N，N为全弹道上特征点的总数；

步骤四、建立高速飞行器全弹道多特征点随控优化模型：

min f₁(X)

min f₂(X)

；（8）

min f₅(X)

s.t.g_i(X)≤0,i＝1,2,...m

式中，g_i(X)≤0,i＝1,2,...m是随控优化模型的约束条件；

通过理想点法将多目标优化问题转换为单目标优化问题，即将多个目标函数合成为单一目标函数；根据理想点法的计算公式可得转换后单目标函数h(X)的表达式：

h(X)=[Σj=15(fj(x)-fj0)2fj0]1/2---(9)]]>

其中，为函数f_j(X)的理想值；

利用理想点法得到单目标函数的随控优化模型如下：

min h(X)（10）

s.t.g_i(X)≤0,i＝1,2,...m

利用模拟退火算法在单目标函数随控优化模型的可行域内搜索满足气动非线性度与耦合度评价指标要求的优化解；模拟退火算法利用模拟热力学中经典粒子系统的降温过程来求解规划问题的极值，其基本步骤为：

①给定冷却进度表参数及其迭代初始解X₀和h(X₀)；

②参数t＝t_k时作L_k次试探搜索；

根据当前解X_k产生一个随机向量Z_k得到X_k邻域的新的试探点X′_k：

X′_k＝X_k+Z_k

③产生一个在（0,1）上均匀分布的随机数η，计算出在给定当前迭代点X_k和温度t_k下与Metropolis接受准则相对应的转移概率P，有：

P=1,h(Xk′)≤h(Xk)eh(Xk)-h(Xk′)tk,h(Xk′)>h(Xk)]]>

若η≤P，则接受新解X_k＝X′_k，h(X_k)＝h(X′_k)，否则X_k不变；

④试探点搜索小于L_k次，转步骤②，否则转步骤⑤；

⑤满足迭代终止条件？

是，则停止，当前解为近似全局最优解；否，则转步骤⑥；

⑥根据给定的温度衰减函数，产生新的温度控制参数t_k+1及Markov链的长度L_k+1；⑦重复步骤②-⑥，直到找到最优解。