CN112668252B - 一种系留气球动态响应特性计算方法 - Google Patents

Abstract

本发明公开了一种系留气球动态响应特性计算方法，包括：建立系留气球刚体运动方程，包括：把系留气球系统作为质点，其运动为刚体运动，确定系留气球对于稳定性坐标系中的质心的运动学和动力学方程；对系留气球进行受力分析，建立系留气球的动力学模型；确定外力和力矩的表达式求和并整理，得出绕气球质心的六自由度运动方程；对离散突风、大气紊流两种外界风场进行建模，为系留气球六自由度仿真提供风场环境输入；将建立的大气风场模型以函数的形式加载到系留气球的运动方程中得到状态空间函数，使其成为已知的并对系留气球运动状态产生影响的控制参量；求解此状态空间函数，得到系留气球在大气风场扰动情况下的动态响应过程。

Description

一种系留气球动态响应特性计算方法

技术领域

本发明涉及一种系留气球稳定性的计算方法，适用于所有系留气球的动态响应特性计算分析。

背景技术

系留气球系统通过携带光学、微波等遥感载荷和无线通信载荷，为战场侦察监视、实时监视、反恐维稳、防灾减灾、环境监测和高速通信等应用需求提供崭新的技术手段。在设计过程中系留气球动态响应特性若无法有效评估，气球驻空过程中受到风扰后的姿态及缆绳状态都处于未知状态，气球及任务载荷在抗风过程中将面临着损坏的风险。

目前，国内对系留气球稳定性的研究以定常静平衡特性仿真计算为主，针对系留气球动态响应特性计算方法的研究较少。气球动态响应特性的计算要求建立准确的数学模型，与固定翼飞行器不同，系留气球还受到浮力及缆绳牵引力的作用，尤其缆绳的动力学建模难度较大，模型不准将导致计算的精度下降。

发明内容

本发明的目的是提供一种系留气球动态响应特性计算方法，用以克服现有算法因模型不准确导致计算精度下降的问题，从而为工程设计及优化提供有效参考。

为了实现上述任务，本发明采用以下技术方案：

一种系留气球动态响应特性计算方法，包括：

建立系留气球刚体运动方程，包括：把系留气球系统作为质点，其运动为刚体运动，确定系留气球对于稳定性坐标系中的质心的运动学和动力学方程；

对系留气球进行受力分析，建立系留气球的动力学模型；确定外力和力矩的表达式求和并整理，得出绕气球质心的六自由度运动方程；

对离散突风、大气紊流两种外界风场进行建模，为系留气球六自由度仿真提供风场环境输入；

将建立的大气风场模型以函数的形式加载到系留气球的运动方程中得到状态空间函数，使其成为已知的并对系留气球运动状态产生影响的控制参量；求解此状态空间函数，得到系留气球在大气风场扰动情况下的动态响应过程。

进一步地，所述建立系留气球的动力学模型，包括：

附加质量，包括囊体附加质量以及尾翼的附加质量，通过理论计算方法求得；其中计算囊体附加质量时，将系留气球用一个三轴椭球体替代；计算尾翼的附加质量时，尾翼按当量平板来计算其附加质量；

气动力和力矩，首先通过由质心运动引起的参考点上的运动，然后计算由于作用在参考点上的力和力矩而在质心处产生的力和力矩，来确定作用在质心上的力和力矩系数；

系留缆绳力和力矩，首先确定系留缆绳力的基本表达式，然后确定在稳定性坐标系上绕气球质心的系留缆绳力矩，用气球质心的位移加缆绳铰链点绕质心的转动表示出系留绳顶端的位移；在小扰动情况下，求解并忽略高阶项，结合纵向、横向缆绳的导数表达式，得到系留缆绳力和力矩的表达式；

浮力和力矩，在假设有小扰动角的情况下，建立浮力和力矩的表达式；

重力和力矩，仅考虑因气球结构重量产生的分力，附加质量和气体的升力已分别包括在加速度项和浮力项系数中；在小扰动角的情况下，建立重力和力矩的表达式。

进一步地，所述绕气球质心的六自由度运动方程，表示为：

x向力方程：

z向力方程：

滚转力矩方程：

俯仰力矩方程：

偏航力矩：

其中，m_x、m_y、m_z分别表示系留气球在x向、y向、z向上的质量，x′,y′,z′分别表示地面坐标系上的扰动位移；ρ表示大气密度，V_∞表示定常风速，S表示系留气球特征面积，c_D表示阻力系数，

表示气动阻力对速度的导数，x,y,z表示系留气球在稳定性坐标系上的扰动位移，k_xx表示在缆绳铰链点每单位X方向位移的系留缆绳的力，

表示气动阻力对迎角的导数，c_L表示升力系数，k_xz表示系留气球每单位Z方向位移的系留缆绳的力，k_xθ表示系留气球每单位俯仰角位移的系留缆绳的力，B表示净浮力(总浮力-气体总重量)，W_S表示系留气球结构重量，T₁表示系留缆绳在上端的张力，γ₁表示上端的系留缆绳与水平线之间的夹角，θ表示俯仰角的扰动角；

表示系留气球长度，

表示侧力对侧滑角速率的导数，

表示侧力系数导数，k_yy表示在缆绳铰链点每单位Y方向位移的系留缆绳的力，

表示滚转角速度引起的侧力系数导数，φ表示滚转的扰动角，k_yφ表示每单位滚转角位移的系留缆绳的力，ψ表示偏航的欧拉角，k_yψ表示每单位偏航角位移的系留缆绳的力；

表示迎角变化率引起的升力系数导数，

表示气动升力对速度的导数，k_zx表示每单位X方向位移的系留缆绳的力，

表示升力系数导数，k_zz表示在缆绳铰链点每单位Z方向位移的系留缆绳的力，

表示俯仰角速度引起的升力系数导数,k_zθ表示每单位俯仰角位移的系留缆绳的力；

I_x、I_y、I_z分别表示在稳定性坐标系上分别绕气球质心的滚转，俯仰和偏航惯性矩，

表示滚转力矩对侧滑角速率的导数，I_xy、I_xz、I_yz分别表示在XOY平面，XOZ平面和YOZ平面绕气球质心的惯性积；

表示滚转力矩系数导数，k_φy表示每单位Y方向位移的缆绳滚转力矩，

表示滚转角速度引起的滚转力矩系数导数，h_k1、

表示方程定义的坐标，k_φφ表示绕质心滚转的体轴系上每单位滚转角位移的缆绳滚转力矩，

表示自定义力矩；

表示偏航角速度引起的滚转力矩系数导数，k_φψ表示每单位偏航角位移的缆绳滚转力矩；

表示迎角变化率引起的俯仰力矩系数导数，C_m表示俯仰力矩系数，

表示俯仰力矩对速度的导数，k_θx、k_θz每单位X方向、Z方向位移的缆绳俯仰力矩，

表示俯仰力矩系数导数，

表示俯仰角速度引起的俯仰力矩系数导数，k_θθ表示绕质心俯仰的体轴系上每单位俯仰角位移的缆绳俯仰力矩；

表示偏航力矩对侧滑角速率的导数，

表示偏航力矩系数导数，k_ψy表示每单位Y方向位移的缆绳偏航力矩，

表示滚转角速度引起的偏航力矩系数导数，k_ψφ表示每单位滚转角位移的缆绳偏航力矩，

表示偏航角速度引起的偏航力矩系数导数，k_ψψ表示绕质心偏航的体轴系上每单位偏航角位移的缆绳偏航力矩。

进一步地，所述系留气球对于稳定性坐标系中的质心的运动学和动力学方程表示为：

其中，m表示系留气球的质量，F_X,F_Y,F_Z表示作用在系留气球分别平行于X，Y和Z轴上的外力，M_X,M_Y,M_Z表示分别绕X，Y和Z轴的滚转，俯仰和偏航力矩，x′,y′,z′分别表示地面坐标系上的扰动位移，φ表示滚转的扰动角，ψ表示偏航的扰动角，θ表示俯仰的扰动角，参数中上标原点表示相对于时间的导数；I_x,I_y,I_z表示在稳定性坐标系上分别绕气球质心的滚转，俯仰和偏航惯性矩，I_xz表示在XOZ平面绕气球质心的惯性积。

进一步地，对所述外界风场中的离散突风进行建模的过程包括：

突风形状和强度定义如下：

其中：

U_m表示突风速度(m/s)；X表示进入突风区距离，0≤X≤2H；H表示突风梯度长度，L/4≤X≤243.8；L表示飞艇长度；

突风强度均采用《飞艇适航标准》对突风强度的规定，全波长型突风和半波长型突风采用《飞艇适航标准》的尺寸的规定。

进一步地，所述外界风场中的大气紊流模型采用Von Karman紊流模型。

进一步地，所述求解此状态空间函数，得到系留气球在大气风场扰动情况下的动态响应过程，包括纵向动态响应求解以及横航向动态响应求解，其中：

纵向动态响应求解包括纵向本征方程求解以及纵向加载风场方程求解：

所述纵向动态响应求解中，设运动方程的平衡部分等于零，可得出稳定性方程；对于所述六自由度运动方程，其中三个方程仅有纵向变量，另外三个方程仅有横向变量，因此，得出绕气球质心的稳定性方程的工作形式；借助MATLAB对其特征矩阵求解，即可解出其六个特征根；

将运动方程写成矩阵形式，将所建立的外界风场的模型以函数的形式加载到系留气球的运动方程中，求解可得到系留气球在大气风场扰动情况下的动态响应过程。

进一步地，所述方法以计算机程序的形式装载于计算机中；所述计算机包括处理器、存储器，处理器执行计算机程序时，实现所述方法的步骤。

进一步地，所述方法以计算机程序的形式装载于计算机可读存储介质中；计算机程序被处理器执行时，实现所述方法的步骤。

与现有技术相比，本发明具有以下技术特点：

1.本发明提出了系留气球非线性动态响应特性计算方法，建立完整的飞行动力学模型、几何模型、风场模型，能够准确输出实时的各项飞行参数变化情况，绘制本体特性曲线，计算简便快捷，在实际操作中能够迅速给出评估结果，具有较好的操纵性、通用性和可扩展性。

2.应用本方法进行动态响应特性仿真计算，能够评估出气球在不同风场环境下的动态响应特性，为气动优化、稳定性设计提供了有效反馈，使设计的气球能保持理想姿态完成驻空任务。同时，本方法适用于所有类型系留气球，计算简便，在设计中可快速评估动态响应特性，缩短了研制周期。

附图说明

图1为本发明方法的计算过程示意图；

图2为系留气球缆绳受力分析图；

图3为稳定性坐标示意图；

图4为稳定系坐标系方位示意图；

图5为气球系统几何图；

图6为缆绳受力示意图；

图7为椭球体附加质量系数K₁₁；

图8为椭球体附加质量系数k₂₂(k₃₃)；

图9为椭球体附加转动惯量系数k₅₅(k₆₆)；

图10为实施例中的纵向本体响应曲线；

图11为实施例中横航向本体响应曲线；

图12为垂直突风扰动下缆绳上端张力仿真曲线；

图13为垂直突风扰动下气球俯仰角仿真曲线；

图14为仿真计算模型示意图。

具体实施方式

参见附图，本发明提供了一种系留气球动态响应特性计算方法，包括以下步骤。

本方案中，所涉及公式中的参数解释如下：

不同变量下标含义：

A 空气动力项

B 浮力项

C 系留缆绳力项

G 重力项

R 参考点(见图5)

t 配平状态

0 系留缆绳的下端(见图6)

1 系留缆绳的上端(见图6)

气动导数所用的下标含义：

α 相对于迎角α的导数

相对于

的导数

β 相对于β的导数

相对于

的导数

p 相对于

的导数

q 相对于

的导数

r 相对于

的导数

u 相对于uV_∞的导数

注：在一个符号上的圆点用来表示相对于时间(t)的导数。

步骤1，建立系留气球刚体运动方程

把系留气球系统作为质点，其运动为刚体运动，运动学和动力学方程相对于稳定性坐标系中的质心来确定，方程如下：

三轴方向力的方程为：

三轴方向的力矩为：

由于xoz轴对称，所以式中I_xy＝I_yz＝0

姿态角速率与体轴系的三个角速度分量之间的关系式：

在地面坐标系中的线速度如下：

上述方程中所包括的轴系、运动模态、力和力矩见图2～图4。

利用本发明中给出的假定，前面的运动方程可适用于描述系留气球的运动。这些假定要求气球质心在平衡状态下不应当有位移或速度。因此，气球表现出的唯一运动是一种绕其起始平衡位置的很小扰动。所以，在很小位移情况下获得以下线性关系式：

U＝u；V＝v；W＝w

Θ＝θ；Ψ＝ψ；Φ＝φ

P＝p；Q＝q；R＝r

式中，小写字母表示质心的小扰动速度和位移。

与普通飞机分析相比较，速度U的x分量的表达式没有包括平衡速度项V_∞。区别的原因在于，在传统的分析中飞机是以平衡速度加扰动速度的速度通过静止空气，而在本分析中气球是由地面固定点系留，气球的唯一运动是一种很小的扰动运动。在本方案分析中，以不变的气流流过气球，但气流速度只会影响到气动力的大小并不是质心运动的明显部分。

通过上式方程可得：

通过带入并整理，可得：

步骤2，对系留气球进行受力分析，建立系留气球的动力学模型；确定外力和力矩的表达式求和并整理，得出绕气球质心的运动方程。

2.1附加质量

物体在流体中作加速运动时，总要促使一部分流体也作加速度运动，从而相当于物体所受的惯性力增加了，或相当于物体的质量增加了，这增加的假想质量，即为附加质量。由于系留气球的体积大，飞行速度慢，产生的附加质量与系留气球本身的质量量级相当，所以研究系留气球的非定常运动时，就必须考虑附加质量。附加质量可通过理论计算方法求得。

2.1.1囊体附加质量

系留气球囊体一般可以用一个三轴椭球体替代。椭球体的三轴半周长分别为a，b，c。由于椭球体λ_ij＝0(i≠j时)，因而所要计算的附加质量只有λ_ij(i＝j时)。椭球体的附加质量可以根据椭球体的半轴长a，b，c在图7、图8和图9所示的曲线上查得，该曲线是兰姆(Lamb)根据势流理论计算得到的，图中K₁₁,K₂₂,K₃₃分别是椭球体附加质量λ₁₁,λ₂₂,λ₃₃的相对值，采用椭球体排开的空气质量m作为基值。

K₅₅,K₆₆分别是椭球体附加转动惯量λ₅₅,λ₆₆的相对值，分别采用椭球体排开的空气的转动惯量I_yy、I_zz作为基值。

(λ_11,22,33)_B＝(k_11,22,33)_B·m_B

(λ_55,66)_B＝(k_55,66)_B·(I_yy,I_zz)_B

式中(m，I_yy，I_zz)_B分别是椭球体所排开空气体积的质量和对y，z轴的转动惯量，可以用下式计算：

式中，a为椭球体的长半轴，b，c为椭球体的短半轴。

当b＝c时，K₄₄＝0。

2.2.2尾翼的附加质量

尾翼按当量(面积S_i，展弦比λ_i)平板来计算它们的附加质量。由于平板的左右、上下对称，故除λ_ij(i＝j)，λ₃₅，λ₂₆外，其他附加质量(λ_ij)_w＝0(i≠j)，考虑到λ₃₅，λ₂₆值较小，所以忽略不计。

系留气球的尾翼为“倒Y”形尾翼，由于垂尾尾翼平板纵向绕流时的附加质量等于零，所以(λ₁₁)_w＝0。同理，沿z轴方向的m₃₃和沿y轴的(λ₂₂)_w，绕z轴的(λ₅₅)_w，绕y轴的(λ₆₆)_w，其余的(λ_ij)_w(i≠j)的项也均为零。

各尾翼的附加质量计算，以每片尾翼作为当量平板，公式如下：

式中：l——平板的展长，即尾翼的展长

b——平板的弦长，即尾翼的弦长

μ(λ)——有限翼展的修正值，按巴布斯特(ΠaδcTa)半经验公式计算

式中，λ＝l/b＝l²/s为尾翼的展弦比，s为单片外露尾翼面积。

由于系留气球的尾翼为“倒Y”型尾翼，所以

(λ₃₃)_w＝2·(λ₃₃)_swcos(30°)

(λ₂₂)_w＝(λ₂₂)_sw+(λ₂₂)_swsin(30°)

式中，x_p为尾翼压心对于坐标原点的坐标值，由于尾翼为当量矩形，所以尾翼压心即为矩形的形心。

2.2气动力和力矩

首先通过研究由质心运动引起的参考点上的运动，然后计算由于作用在参考点上的力和力矩而在质心处产生的力和力矩，来确定作用在质心上的力和力矩系数。

参考点上的力和力矩项转换到质心处为：

F_Z＝F_Z,R

F_X＝F_X,R

M_Y＝M_Y,R-z_tF_X,R+x_tF_Z,R

在z方程中所需使用的系数为：

(F_Z,A)_R＝-(L+Dα_R)

转换关系如下：

在x方程中使用所需的系数是：

(F_X,A)_R＝Lα_R-D

俯仰力矩系数如下：

式中，

项被忽略。

对有关质心的横向扰动，在参考点上的扰动为：

φ_R＝φ

ψ_R＝ψ

参考点上的力和力矩项转到质心上为：

F_Y＝F_Y,R

M_Z＝M_Z,R-x_tF_Y,R

M_X＝M_X,R+z_tF_Y,R

Y-向力的系数如下：

偏航力矩系数如下：

除

本身外，忽略涉及

的转换。

滚转力矩系数如下：

除

本身外，忽略涉及

的转换。

空气动力和力矩表示如下：

扰动迎角α和侧滑角β的一阶近似值表示如下：

把该近似值应用到上式并把α扩展成

和θ项，再把β扩展成

和ψ项，便获得以下结果：

对于很小的扰动迎角α，在空气动力升力和阻力项中可以把纵向力表示为：

F_X,A＝(Lα-D)

F_Z,A＝-(L+Dα)

通过以上方法，在平衡状态下评估这些方程中得出的空气动力系数并省略高阶扰动项，便可得到以下关系式：

2.3系留缆绳力和力矩

系留缆绳力的基本表达式在上文中已推导出。在稳定性坐标系上缆绳力为：

F_X,C＝F_X',C+F_Y',Cψ-F_Z',Cθ

F_Y,C＝-F_X',Cψ+F_Y',C+F_Z',Cφ

F_Z,C＝F_X',Cθ-F_Y',Cφ+F_Z',C

因此，在稳定性坐标系上绕气球质心的系留缆绳力矩为：

可用气球质心的位移加缆绳铰链点绕质心的转动表示出系留绳顶端的位移x_c，y_c和z_c。在小扰动情况下，其可表示为：

把上述方程求解并忽略高阶项，可得：

F_X,C＝-k_xxx'-k_xzz'-(k_xθ+T₁sinγ₁)θ+T₁cosγ₁

F_Z,C＝-k_zxx'-k_zzz'+(T₁cosγ₁-k_zθ)θ+T₁sinγ

式中，k_xx、k_xz、k_zx和k_zz为纵向缆绳的导数(弹簧常数),表达式为：

其中：

δ＝x₁(sinγ₁-sinγ₀)+z₁(cosγ₀-cosγ₁)-lsin(γ₁-γ₀)

k_yy为横向缆绳导数(弹簧常数)，表达式为：

其中：

τ₁为τ(γ＝γ₁)的值。

x₁、z₁、γ₁、γ₀见系留气球静平衡方程。

k_φy＝k_yφ

k_φφ＝k_yφ

k_ψy＝k_yψ

k_ψφ＝k_φψ

2.4浮力和力矩

假设有小扰动角，其为：

F_X,B＝Bθ

F_Y,B＝-Bφ

F_Z,B＝-B

M_X,B＝-B[(h_cg-h_br)cosα_t+(l_br-l_cg)sinα_t]φ

M_Y,B＝B[(l_br-l_cg)cosα_t-(h_cg-h_br)sinα_t]-B[(h_cg-h_br)cosα_t+(l_br-l_cg)sinα_t]θ

M_Z,B＝-B[(l_br-l_cg)cosα_t-(h_cg-h_br)sinα_t]φ

2.5重力和力矩

在推导运动方程中的重力时，只需考虑因气球结构重量产生的分力。附加质量和气体的升力已分别包括在加速度项和浮力项系数中。

在小扰动角情况下，重力产生的力和力矩可表示为：

F_X,G＝-W_Sθ

F_Y,G＝W_Sφ

F_Z,G＝W_S

M_X,G＝-W_S[(h_sr-h_cg)cosα_t+(l_sr+l_cg)sinα_t]φ

M_Y,G＝W_S[(l_sr+l_cg)cosα_t-(h_sr-h_cg)sinα_t]-W_S[(h_sr-h_cg)cosα_t+(l_sr+l_cg)sinα_t]θ

M_Z,G＝-W_S[(l_sr+l_cg)cosα_t-(h_sr-h_cg)sinα_t]φ

2.6系留气球六自由度运动方程

将外力和力矩的表达式求和并整理，则可得出如下绕气球质心的六自由度运动方程。

x向力方程：

y向力方程：

z向力方程：

滚转力矩方程：

俯仰力矩方程：

偏航力矩：

步骤3，由于系留气球本身没有操纵系统，所以外界输入即为风场闭环输入。系统的运动时间过渡过程就是系统受到外界风场扰动后的运动响应过程，对以下两种外界风场进行建模，为系留气球六自由度仿真提供风场环境输入：

3.1离散突风

突风又称阵风，表现为确定性风速变化，在工程研究中，突风可表征离散的风切变、大气紊流中的峰值、地形诱导的气流等，突风可单独使用，也可以叠加到平均风和/或紊流上。

一般来说，在实际飞行中，半波长(1-cosin)型突风数学模型较为常见，且《飞艇适航标准》(AC-21-09)飞行载荷3.9节突风载荷对突风形状和强度的定义如下：

假设飞艇在平飞时承受遇到下列大气突风所产生的载荷：

1)当以V_H速度飞行时，7.62m/s的离散突风。

2)当以V_B速度飞行时，10.67m/s的离散突风。

3)突风形状和强度定义如下：

其中：

U_m＝上述规定的突风速度(m/s)；

X＝进入突风区距离，0≤X≤2H(m)

H＝突风梯度长度，L/4≤X≤243.8(m)；和

L＝飞艇长度(m)

V_H：设计最大平飞速度

V_B：对应最大突风强度的设计空速。

以上离散突风强度均采用《飞艇适航标准》对突风强度的规定，全波长(1-cosin)型突风和半波长(1-cosin)型突风采用《飞艇适航标准》的尺寸的规定，其他类型的突风，根据实际情况确定时间间隔。以上规定同样适用于系留气球对突风的要求。

2.2大气紊流

1)大气紊流的性质

大气的总速度场是既随空间而变化又随时间而变化的。一般，大气的速度可以看作由平均风速和紊流速度组成的。平均风速主要影响系留气球的航迹和性能以及导航等问题。本方案只关心大气紊流对系留气球动态特性的影响。

2)简化假设

实际的大气紊流是十分复杂的物理现象。为了使研究适当简化，通常采用以下假设：

i.认为大气紊流是高斯型的，即速度是正态分布的；

ii.认为大气紊流是平稳的。亦即其统计特性不随时间而变化。而且还认为大气紊流是均匀的，亦即其统计特性不随位置而变。把以上两个性质结合起来，就导致这样的结果：系留气球所经受的大气紊流作用是平稳随机过程；

iii.认为大气紊流是遍历性的随机过程，因而可以用时间平均值代替集合平均值；

iv.认为大气紊流是各向同性的，即所有的统计特性是与坐标轴的方向无关。因而对于任何坐标系，紊流强度(即紊流速度均方根值)在三个方向上是相等的；

v.认为大气紊流速度场是“冻结的”。

3)大气紊流数学模型

工程上对连续紊流模型，主要有两种：Dryden模型和Von Karman模型。它们的建立的理论体系不同，Dryden模型是先建立大气紊流的相关函数，再推导出频谱表达式，而VonKarman模型则是根据大量的测量统计结果建立大气紊流的频谱函数，再推导出相关函数，显然Von Karman模型更能真实反映大气紊流实际情况；对于工程上而言，Dryden模型与VonKarman模型差别很小，如果再考虑测量的误差，可以认为，这两种模型的差别对于工程来说是微不足道的，本方案使用Von Karman紊流模型作为风场扰动模型。

Von Karman模型的速度频谱函数为：

Von Karman模型的角速度频谱函数为：

上式中σ_u、σ_v、σ_w代表坐标轴(x,y,z)方向上的紊流强度，L_u、L_v、L_w代表坐标轴(x,y,z)方向上的紊流尺度，a为常量值1.339、b代表横航向的参考长度，Ω代表速度频谱。

步骤4，将步骤3所建立的以速度表示的大气风场模型以函数u的形式加载到系留气球的运动方程中得到状态空间函数，使其成为已知的并对系留气球运动状态产生影响的控制参量；求解此状态空间函数，即可得到系留气球在大气风场扰动情况下的动态响应过程。

4.1纵向动态响应求解

1)纵向本征方程求解

设运动方程的平衡部分等于零，可得出稳定性方程。其余表达式给出六个线性稳定性扰动方程，其中三个方程仅有纵向变量，另外三个方程仅有横向变量。因此，得出绕气球质心的稳定性方程的工作形式。

纵向运动方程可按以下形式写出：

将上述方程降阶：

则，

且V_K＝V_A+V_W，将上式整理为

的形式，则变为：

其中，

x＝[u_k w_k q_k x′ z′ θ]^-1

则

借助MATLAB对其特征矩阵求解，即可解出其六个特征根，其表现为振荡运动模态的共轭复根(即λ＝η±iw)，或是表现为非周期运动模态的实根(即λ＝η)。因此，在纵向情况下，根据方程根是实根或共轭复根，气球可能呈现三个到六个运动模态，对于任何一种运动模态，当η<0时气球系统是稳定的，当η＝0时是中立稳定的，而η>0时是不稳定的。

2)纵向加载风场方程求解

将系统加入控制阵，写成

的形式，其中A、B阵和x如上所述，C为单位矩阵，D为零阵，E为风场系数阵，如下所示：

u为外界信号输入，在本方案中，u即为外界扰动风场，可分为两种情况，一种情况为离散突风，一种情况为连续大气紊流，即可将前文所建立的以速度表示的大气风场模型以函数u的形式加载到系留气球的运动方程中，使其成为已知的并对系留气球运动状态产生影响的控制参量。

求解此状态空间函数，即可得到系留气球在大气风场扰动情况下的动态响应过程。

考虑缆绳在平衡位置，如果缆绳顶端在xoz平面内从其初始位置(x₁，z₁)位置移到新位置(x₁+Δx,z₁+Δz)，那么x向和z向的合力增量为：

ΔF_x＝k_xxΔx+k_xzΔz＝(T₁+ΔT₁)cos(γ₁+Δγ₁)-T₁cosγ₁

≈ΔT₁cosγ₁-T₁sinγ₁·Δγ₁

ΔF_z＝k_zxΔx+k_zzΔz＝(T₁+ΔT₁)sin(γ₁+Δγ₁)-T₁sinγ₁

≈ΔT₁sinγ₁+T₁cosγ₁·Δγ₁

由上式可得：

ΔT₁＝(k_xxcosγ₁+k_zxsinγ₁)Δx+(k_zzsinγ₁+k_xzcosγ₁)Δz

至此，已解出了在大气风场扰动下的纵向运动参数Δx、Δz、Δθ、ΔT₁和Δγ₁的时间响应过程。

4.2横航向动态响应求解

1)横航向本征方程求解

设运动方程的平衡部分等于零，可得出横向稳定性方程。因此，得出绕气球质心的横向稳定性方程的工作形式。

横向运动方程可按以下形式写出：

将方程降阶：

则

经整理，将上式整理为

的形式，则变为：

其中，

则

解此方程，得6个特征根，这些特征根表现为振荡运动模态的共轭复根(即λ＝η±iw)，或是表现为非周期运动模态的实根(即λ＝η)。因此，在横向情况下，根据方程根是实根或共轭复根，气球可能呈现三个到六个运动模态，对于任何一种运动模态，当η<0时气球系统是稳定的，当η＝0时是中立稳定的，而η>0时是不稳定的。

2)横航向加载风场求解

将系统加入控制阵，写成

的形式，其中A、B阵和x如上所述，C为单位矩阵，D为零阵，E为风场系数阵，如下所示：

u为外界信号输入，本方案中即为外界风场，其受外界风场扰动后，其运动求解过程同纵向运动方程的求解过程。

步骤5，系留气球动态响应特性仿真

以体积为2000m³的系留气球为例，计算了气球在15m/s的来流风速扰动下的动态响应特性及受到7.62m/s的垂直突风时缆绳上端张力及气球俯仰角的变化情况，如图10至图13所示。

基于步骤1～5中的公式编制系留气球动力学模型计算程序，加载风场模型后，导入步骤4中4.1-2)与4.2-2)的特征方程(见图14状态空间模块)，从而计算出在风场扰动下各状态量的动态响应特性(见图14动力学解算模块)，得出步骤5的动态响应特性曲线。

以上实施例仅用于说明本申请的技术方案，而非对其限制；尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行同等替换；而这些修改或替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的精神和范围，均应包含在本申请的保护范围之内。

Claims (8)

1.一种系留气球动态响应特性计算方法，其特征在于，包括：

建立系留气球刚体运动方程，包括：把系留气球系统作为质点，其运动为刚体运动，确定系留气球对于稳定性坐标系中的质心的运动学和动力学方程；

对系留气球进行受力分析，建立系留气球的动力学模型；确定外力和力矩的表达式求和并整理，得出绕气球质心的六自由度运动方程；

对离散突风、大气紊流两种外界风场进行建模，为系留气球六自由度仿真提供风场环境输入；

将建立的大气风场模型以函数的形式加载到系留气球的运动方程中得到状态空间函数，使其成为已知的并对系留气球运动状态产生影响的控制参量；求解此状态空间函数，得到系留气球在大气风场扰动情况下的动态响应过程；

所述建立系留气球的动力学模型，包括：

附加质量，包括囊体附加质量以及尾翼的附加质量，通过理论计算方法求得；其中计算囊体附加质量时，将系留气球用一个三轴椭球体替代；计算尾翼的附加质量时，尾翼按当量平板来计算其附加质量；

气动力和力矩，首先通过由质心运动引起的参考点上的运动，然后计算由于作用在参考点上的力和力矩而在质心处产生的力和力矩，来确定作用在质心上的力和力矩系数；

系留缆绳力和力矩，首先确定系留缆绳力的基本表达式，然后确定在稳定性坐标系上绕气球质心的系留缆绳力矩，用气球质心的位移加缆绳铰链点绕质心的转动表示出系留绳顶端的位移；在小扰动情况下，求解并忽略高阶项，结合纵向、横向缆绳的导数表达式，得到系留缆绳力和力矩的表达式；

浮力和力矩，在假设有小扰动角的情况下，建立浮力和力矩的表达式；

重力和力矩，仅考虑因气球结构重量产生的分力，附加质量和气体的升力已分别包括在加速度项和浮力项系数中；在小扰动角的情况下，建立重力和力矩的表达式。

2.根据权利要求1所述的系留气球动态响应特性计算方法，其特征在于，所述绕气球质心的六自由度运动方程，表示为：

x向力方程：

y向力方程：

z向力方程：

滚转力矩方程：

俯仰力矩方程：

偏航力矩：

其中，m_x、m_y、m_z分别表示系留气球在x向、y向、z向上的质量，x′,y′,z′分别表示地面坐标系上的扰动位移；ρ表示大气密度，V_∞表示定常风速，S表示系留气球特征面积，c_D表示阻力系数，

表示气动阻力对速度的导数，x,y,z表示系留气球在稳定性坐标系上的扰动位移，k_xx表示在缆绳铰链点每单位X方向位移的系留缆绳的力，

表示气动阻力对迎角的导数，c_L表示升力系数，k_xz表示系留气球每单位Z方向位移的系留缆绳的力，k_xθ表示系留气球每单位俯仰角位移的系留缆绳的力，B表示净浮力，W_S表示系留气球结构重量，T₁表示系留缆绳在上端的张力，γ₁表示上端的系留缆绳与水平线之间的夹角，θ表示俯仰角的扰动角；

表示系留气球长度，

表示侧力对侧滑角速率的导数，

表示侧力系数导数，k_yy表示在缆绳铰链点每单位Y方向位移的系留缆绳的力，

表示滚转角速度引起的侧力系数导数，φ表示滚转的扰动角，k_yφ表示每单位滚转角位移的系留缆绳的力，ψ表示偏航的欧拉角，k_yψ表示每单位偏航角位移的系留缆绳的力；

表示迎角变化率引起的升力系数导数，

表示气动升力对速度的导数，k_zx表示每单位X方向位移的系留缆绳的力，

表示升力系数导数，k_zz表示在缆绳铰链点每单位Z方向位移的系留缆绳的力，

表示俯仰角速度引起的升力系数导数,k_zθ表示每单位俯仰角位移的系留缆绳的力；

I_x、I_y、I_z分别表示在稳定性坐标系上分别绕气球质心的滚转，俯仰和偏航惯性矩，

表示滚转力矩对侧滑角速率的导数，I_xy、I_xz、I_yz分别表示在XOY平面，XOZ平面和YOZ平面绕气球质心的惯性积；

表示滚转力矩系数导数，k_φy表示每单位Y方向位移的缆绳滚转力矩，

表示滚转角速度引起的滚转力矩系数导数，h_k1、

表示方程定义的坐标，k_φφ表示绕质心滚转的体轴系上每单位滚转角位移的缆绳滚转力矩，

表示自定义力矩；

表示偏航角速度引起的滚转力矩系数导数，k_φψ表示每单位偏航角位移的缆绳滚转力矩；

表示迎角变化率引起的俯仰力矩系数导数，C_m表示俯仰力矩系数，

表示俯仰力矩对速度的导数，k_θx、k_θz每单位X方向、Z方向位移的缆绳俯仰力矩，

表示俯仰力矩系数导数，

表示俯仰角速度引起的俯仰力矩系数导数，k_θθ表示绕质心俯仰的体轴系上每单位俯仰角位移的缆绳俯仰力矩；

表示偏航力矩对侧滑角速率的导数，C_nβ表示偏航力矩系数导数，k_ψy表示每单位Y方向位移的缆绳偏航力矩，

表示滚转角速度引起的偏航力矩系数导数，k_ψφ表示每单位滚转角位移的缆绳偏航力矩，

表示偏航角速度引起的偏航力矩系数导数，k_ψψ表示绕质心偏航的体轴系上每单位偏航角位移的缆绳偏航力矩。

3.根据权利要求1所述的系留气球动态响应特性计算方法，其特征在于，所述系留气球对于稳定性坐标系中的质心的运动学和动力学方程表示为：

其中，m表示系留气球的质量，F_X,F_Y,F_Z表示作用在系留气球分别平行于X，Y和Z轴上的外力，M_X,M_Y,M_Z表示分别绕X，Y和Z轴的滚转，俯仰和偏航力矩，x′,y′,z′分别表示地面坐标系上的扰动位移，φ表示滚转的扰动角，ψ表示偏航的扰动角，θ表示俯仰的扰动角，参数中上标原点表示相对于时间的导数；I_x,I_y,I_z表示在稳定性坐标系上分别绕气球质心的滚转，俯仰和偏航惯性矩，I_xz表示在XOZ平面绕气球质心的惯性积。

4.根据权利要求1所述的系留气球动态响应特性计算方法，其特征在于，对所述外界风场中的离散突风进行建模的过程包括：

突风形状和强度定义如下：

其中：U_m表示突风速度，m/s；X表示进入突风区距离，0≤X≤2H；H表示突风梯度长度，L/4≤X≤243.8；L表示飞艇长度；

突风强度均采用《飞艇适航标准》对突风强度的规定，全波长型突风和半波长型突风采用《飞艇适航标准》的尺寸的规定。

5.根据权利要求1所述的系留气球动态响应特性计算方法，其特征在于，所述外界风场中的大气紊流模型采用Von Karman紊流模型。

6.根据权利要求1所述的系留气球动态响应特性计算方法，其特征在于，所述求解此状态空间函数，得到系留气球在大气风场扰动情况下的动态响应过程，包括纵向动态响应求解以及横航向动态响应求解，其中：

纵向动态响应求解包括纵向本征方程求解以及纵向加载风场方程求解：

所述纵向动态响应求解中，设运动方程的平衡部分等于零，可得出稳定性方程；对于所述六自由度运动方程，其中三个方程仅有纵向变量，另外三个方程仅有横向变量，因此，得出绕气球质心的稳定性方程的工作形式；借助MATLAB对其特征矩阵求解，即可解出其六个特征根；

将运动方程写成矩阵形式，将所建立的外界风场的模型以函数的形式加载到系留气球的运动方程中，求解可得到系留气球在大气风场扰动情况下的动态响应过程。

7.一种终端设备，包括处理器、存储器以及存储在所述存储器中的计算机程序；其特征在于，处理器执行计算机程序时，实现所述根据权利要求1至6中任一权利要求所述方法的步骤。

8.一种计算机可读存储介质，其特征在于，所述存储介质中存储有计算机程序；计算机程序被处理器执行时，实现所述根据权利要求1至6中任一权利要求所述方法的步骤。

Images