[发明专利]一种高超声速飞行器滑翔段弹道解析求解方法

权利要求书

1.一种高超声速飞行器滑翔段弹道解析求解方法，包括以下几个步骤：

步骤1：定义当地地心坐标系；

当地地心坐标系S_e1的坐标原点位于地心，x_e1由地心指向飞行器；y_e1在由x_e1与V的平面内，且与V同一侧，z_e1由右手定则确定；

地心赤道旋转坐标系S_e的坐标原点也为地心，z_e由地心北极，x_e处在赤道平面内指向本初子午线，y_e由右手定则确定；

（1）获取转换矩阵；

当地地心坐标系S_e1与地心赤道旋转坐标系S_e之间的转换矩阵为：

Le1e=Lz(ω1)Lx(i1)Lz(Ω1)=cosω1cosi1cosΩ1-sinω1sinΩ1cosω1cosi1sinΩ1+sinω1cosΩ1-cosω1sini1-sinω1cosi1cosΩ1-cosω1sinΩ1-sinω1cosi1sinΩ1+cosω1cosΩ1sinω1sini1sini1cosΩ1sini1sinΩ1cosi1---(1)]]>

其中：表示当地地心坐标系与地心赤道旋转坐标系之间的转换矩阵；L_z(ω₁)表示绕z轴旋转ω₁对应的转换矩阵；L_x(i₁)表示绕x轴旋转i₁对应的转换矩阵；L_z(Ω₁)表示绕z轴旋转Ω₁对应的转换矩阵；ω₁为近地点幅角，Ω₁为升交点赤经，i₁为轨道倾角，具体为：

i₁＝arccos(cosφsinψ)               (2)

Ω1=arccos(cosθcosψ+sinφsinψsinθcos2ψ+sin2φsin2ψ)sinθcosψ-sinφsinψcosθcos2ψ+sin2φsin2ψ≥02π-arccos(cosθcosψ+sinφsinψsinθcos2ψ+sin2φsin2ψ)sinθcosψ-sinφsinψcosθcos2ψ+sin2φsin2ψ<0---(3)]]>

ω1=arccos(cosφcosψcos2ψ+sin2φsin2ψ)judge>02π-arccos(cosφcosψcos2ψ+sin2φsin2ψ)judge≤0---(4)]]>

θ为在地心赤道旋转坐标系下的当地经度；φ为地心赤道旋转坐标系下的纬度；ψ为地心赤道旋转坐标系下的行向角；judge为判断近地点幅角正负号的特征量；

（2）坐标转换

设在当地地心坐标系中存在一点[θ₁,φ₁]，并且该点的弹道偏角为ψ₁，假设该点在地心赤道旋转坐标系中的坐标为[θ,φ]，弹道偏角为ψ，则它们满足，

cosφcosθcosφsinθsinφ=(Le1e)Tcosφ1cosθ1cosφ1sinθ1sinφ1---(5)]]>

从而，获得在地心赤道旋转坐标系中的位置为，

φ=arcsin(l11cosφ1cosθ1+l21cosφ1sinθ1+l31sinφ1)θ=arcsin(l12cosφ1cosθ1+l22cosφ1sinθ1+l32sinφ1cosφ)kθ>0arcsin(l12cosφ1cosθ1+l22cosφ1sinθ1+l32sinφ1cosφ)+πkθ≤0---(6)]]>

其中：

l₁₁＝cosω₁cosi₁cosΩ₁-sinω₁sinΩ₁

l₂₁＝-sinω₁cosi₁cosΩ₁-cosω₁sinΩ₁

l₃₁＝sini₁cosΩ₁

l₁₂＝cosω₁cosi₁sinΩ₁+sinω₁cosΩ₁

l₂₂＝-sinω₁cosi₁sinΩ₁+cosω₁cosΩ₁

l₃₂＝sini₁sinΩ₁

k_θ＝l₁₁cosφ₁cosθ₁+l₂₁cosφ₁sinθ₁+l₃₁sinφ₁

单位速度矢量定义如下，

vxvyvz=(Le1e)T[-sinψ1sinθ1-cosψ1sinφ1cosθ1sinψ1cosθ1-cosψ1sinφ1sinθ1cosψ1cosφ1---(7)]]>

另外，在地心赤道旋转坐标系中，

-sinψsinθ-cosψsinφcosθsinψcosθ-cosψsinφsinθcosψcosφ=vxvyvz]]>

由式(7)得，

cosψ=vzcosφsinψ=vx+cosψsinφcosθ-sinθ0∈{0,π}vy+cosψsinφsinθcosθ0∈(0,π)∪(π,2π)---(8)]]>

从而可得，

ψ=arccos(vzcosφ)sinψ≥02π-arccos(vzcosφ)sinψ<0---(9)]]>

其中：v_y为单位速度矢量y轴分量；v_z为单位速度矢量z轴分量；

步骤2：获取在当地地心坐标系下的运动方程

无动力滑翔飞行器单位质量机械能表述为e＝μ/r-V²/2，忽略地球自转，在当地地心坐标系中三自由度运动模型可表述为，

h·=Vsinγθ·1=Vcosγsinψ1(h+R0)cosφ1φ·1=Vcosγcosψ1h+R0e·=DVmγ·=1V[Lcosσm+(V2h+R0-g)cosγ]ψ·1=1V[Lsinσmcosγ+V2h+R0cosγsinψ1tanφ1]s·=Vcosγh+R0---(10)]]>

其中：R₀＝6378km，θ₁和φ₁为当地地心坐标系下的经度和纬度；V为飞行器与地球的相对速度；ψ₁为当地地心坐标系下的弹道偏角；γ为地心赤道旋转坐标系下的弹道倾角；h为地心赤道旋转坐标系下的高度；s为滑翔段的飞行路程；和分别为高度、当地地心坐标系下的经度、当地地心坐标系下的纬度、负比能量、弹道倾角和当地地心坐标系下的弹道偏角及射程对时间的一阶导数；D和L分别为以过载形式表示的升力和阻力，即D＝0.5ρV²SC_D，L＝0.5ρV²SC_L，其中ρ为大气密度，S为气动参考面积，m为飞行器的质量；C_L和C_D分别为升力系数和阻力系数；g为当地重力加速度；σ为速度滚转角；

由于φ₁≈0、γ≈0、ψ₁≈π/2，并假设0.5ρV²SC_N0＝-m(V²/r-g)cosγ且C_Lcosσ＝C_N0+C_N1，其中C_N0为纵向平衡滑翔升力系数、C_N1为升力系数坡度率调整量，式(10)可进一步化为：

h·=Vsine·=DVmθ·1=Vh+R0γ·=ρVSCN12mφ·1=V(π/2-ψ1)h+R0ψ·1=ρVSCY2m+Vh+R0φ1s·=Vcosγh+R0---(11)]]>

其中：C_N1为升力系数坡度率调整量，C_Y为升力系数横向调整量，满足C_Y＝C_Lsinσ；

取ψ₂＝π/2-ψ₁，并对上式进行如下的变换，

dhsinγ=Vdt---(12)]]>

(h+R₀)dθ₁＝Vdt              (13)

(h+R0)dφ1ψ2=Vdt---(14)]]>

2mdγρSCN1=Vdt---(15)]]>

dψ2-ρSCY2m-1h+R0φ1=Vdt---(16)]]>

mdeD=Vdt---(17)]]>

(h+R₀)ds＝Vdt              (18)

其中：dh、dθ₁、dφ₁、de、dγ、dψ₁和ds分别为高度、当地地心坐标系下的经度、当地地心坐标系下的纬度、负比能量、弹道倾角和当地地心坐标系下的弹道偏角及射程的微分，dt为时间的微分；

步骤3：获取滑翔弹道的解析解

设起滑点坐标为(θ₀,φ₀)，起滑高度为h₀，起滑弹道倾角为γ₀，终端高度为h_f，纵向升力系数余量为C_N1，横向升力系数为C_Y；利用式(1)可获得该点处当地地心坐标系与地心赤道旋转坐标系之间的转换关系，并且该点在当地地心坐标系中满足θ₁₀＝0、φ₁₀＝0以及ψ₁₀＝π/2；则在当地地心坐标系下滑翔段弹道的解析解如下；

（1）滑翔弹道倾角解析解

由式(12)和式(15)可得，

dhsinγ=2mdγρSCN1---(19)]]>

对式(19)进行积分可得，

γf-cosγ0=SCN12mβh(ρf-ρ0)---(20)]]>

其中：ρ_f为预测终点大气密度；γ_f为预测终点弹道倾角；β_h为指数大气模型常数；对式(20)进一步处理可得，

SCN12mβhρf-cosγf=SCN12mβhρ0-cosγ0---(21)]]>

从而得到常数，

K*=SCN1ρ02mβh-cosγ0---(22)]]>

上述的K^*即为决定滑翔段纵向弹道弯曲程度的常数；

（2）飞行路程解析解

由式(15)和式(18)可得：

(h+R0)ds=2mdγρSCN1---(23)]]>

将式(21)和式(22)带入式(23)可得：

βh(h+R0)ds=1cosγ+K*dγ---(24)]]>

设为地心平均距离；对式(24)进行积分可得：

βhR‾sf=11+K*1+K*1-K*ln|tanγf2+1+K*1-K*tanγf2-1+K*1-K*|-ln|tanγ02+1+K*1-K*tanγ02-1+K*1-K*|K*|<12(1+K*)1+K*K*-1atan(K*-1K*+1tanγf2)-atan(K*-1K*+1tanγ02)|K*|>12K*-1(cotγf2-cotγ02)|K*|=1---(25)]]>

其中：s_f为预测终点航程；

（3）滑翔经度解析解

由式(13)和式(15)可得，

(h+R0)dθ1=2mdγρSCN1---(26)]]>

式(26)与式(23)比较可得，

dθ₁＝ds             (27)

从而可得，

θ_1f＝s_f             (28)

其中：θ_1f为当地地心坐标系下的预测终点经度；

（4）滑翔纬度与弹道偏角解析解

由式(14)、式(15)和式(16)可得，

dφ1dγ=2m(h+R0)ρSCN1ψ2dψ2dγ=-CYCN1-2m(h+R0)ρSCN1φ1---(29)]]>

将式(21)和式(22)带入式(29)，可得，

dφ1dγ=ψ2R‾βh(K*+cosγ)dψ2dγ=-CYCN1-φ1R‾βh(K*+cosγ)---(30)]]>

取，

K1=1γf-γ0∫γ0γf1R‾βh(K*+cosγ)dγ=θ1fγf-γ0---(31)]]>

则，式（30）可化为，

dφ1dγ=K1ψ2dψ2dγ=-CYCN1-K1φ1---(32)]]>

由于φ₁₀＝0，ψ₂₀＝0，对上式进行拉式变换可得，

sφ1(s)=K1ψ2(s)sψ2(s)=-CYCN11s-K1φ1(s)---(33)]]>

由上式可得，

ψ2(s)=-CYCN1(s2+K12)φ1(s)=-CYCN1K1s(s2+K12)---(34)]]>

对式(34)进行拉式反变换可得，

ψ2f=-CYCN11K1sinθ1fφ1f=-CYCN11K1(1-cosθ1f)---(35)]]>

其中：φ_1f为当地地心坐标系下的预测终点纬度；由上式还可得当地地心坐标系下的预测终点弹道偏角ψ_1f为ψ_1f＝π/2-ψ_2f；利用式(35)和式(1)给出的坐标转换关系，就可以获得在地心赤道旋转坐标系下预测终点经度θ_f预测终点纬度φ_f和预测终点弹道偏角ψ_f；

（5）滑翔速度解析解

假定在滑翔过程中，保持定常纵向升阻比K_N，K_N＝C_Lcosσ/C_D，则滑翔阻力为，

D=LcosσKN---(36)]]>

由于Lcosσ＝0.5ρV²S(C_N0+C_N1)，忽略C_N1后可得，

D=m[g-V2/(R0+h)]cosγKN+0.5ρV2SCN1KN---(37)]]>

其中：D为阻力；m为飞行器质量；V为飞行器相对地球的速度；C_N1为升力系数坡度率调整量；g为当地重力加速度，写成如下形式，

g=μ(R0+h)2]]>

其中：μ为由地球引力常数；

式(37)中，阻力可分为两个部分，即保持平衡滑翔的升力产生的阻力以及用于纵向转弯的升力产生的阻力，忽略第二部分的影响，第一部分阻力满足，

dedθ1=-μR0+h+2eKN---(38)]]>

由于h相对R₀为小量，采用它的算术平均值代替，则上式积分可得，

ln(2e-μ/R‾2e0-μ/R‾)=2θ1fKN---(39)]]>

其中：e₀为初始负比能量。从而可得，

2ef=μ/R‾+(2e0-μ/R‾)exp(2θ1fK)---(40)]]>

其中：e_f为预测终点负比能量；对于另一部分用于转弯的升力，其消耗的机械能为，

Δe=∫s0sVγ·KNrds=∫ttfV2γ·KNdt---(41)]]>

其中：为弹道倾角的导数；

取V‾2=(V02+Vf2+V0Vf)/3,]]>则可得，

Δe=V‾2KN(γf-γ0)---(42)]]>

其中：Δe为负比能量修正量，综上可得，终端速度预测为：

Vf=μ/R‾-(2e0-μ/R‾)exp(2θ1f/KN)-2V‾2KN(γf-γ0)---(43)]]>

其中：V_f为预测终点速度；为距离地心的平均距离；K_N为升阻比；

在式(43)的计算中利用了平均速度，需要通过迭代来不断修正同时，利用迭代修正对K_N的估计；

（6）滑翔攻角和倾侧角解析解

式(21)、式(28)、式(35)、式(43)分别给出了预测终点弹道倾角γ_f、经度θ_f、纬度φ_f、弹道偏角ψ_f、速度V_f与高度h_f之间的关系；利用平衡滑翔关系可知，

CN0=-2m(1/(R0+h)-g/V2)cosγρS---(44)]]>

从而可得攻角和倾侧角分别满足，

CL(α,Ma)=[(CN0+CN1)2+CY2]1/2σ=arcsinCY[(CN0+CN1)2+CY2]1/2---(45)]]>

其中：α为攻角；σ为速度滚转角；Ma为马赫数；

（7）最优初始下滑角度解析解

具体推导如下；

由拟平衡滑翔条件知，

12ρV2SCN0+mV2R0+hcosγ-mμ(R0+h)2cosγ=0---(46)]]>

其中：C_N0为拟平衡滑翔对应的纵向升力系数；假设当γ＝γ^*时，飞行器能够保持平衡滑翔状态，对上式求全微分可得，

(12V2SCN0∂ρ∂h)dhdt+(ρVSCN0)dVdt+(2mVR0+hcosγ)dVdt+(-mV2(R0+h)2cosγ)dhdt+(-mV2R0+hsinγ)dγdt+(2mμ(R0+h)3cosγ)dhdt+(mμ(R0+h)2sinγ)dγdt=0---(47)]]>

将式(11)带入上式可得，

γ*=DmLcosσV2βh2mg-Lcosσm-V42μ---(48)]]>

其中：γ^*为拟平衡滑翔最优下滑弹道倾角；Lcosσ为升力纵向分量；由于(Lcosσ/m)V²β_h(2g)>>(Lcosσm)[1+mV⁴/(2μLcosσ)]，从而上式可进一步简化为：

γ*=g(V2/2)βrcosσ(1KN)---(49)]]>

其中：K_N为纵向升阻比；式(20)、式(25)、式(28)、式(35)、式(43)、式(45)和式(49)共同构成了滑翔段弹道的解析解。