CN113128085B - 一种基于状态空间形式涡格法的阵风减缓分析方法 - Google Patents
一种基于状态空间形式涡格法的阵风减缓分析方法 Download PDFInfo
- Publication number
- CN113128085B CN113128085B CN202110297250.6A CN202110297250A CN113128085B CN 113128085 B CN113128085 B CN 113128085B CN 202110297250 A CN202110297250 A CN 202110297250A CN 113128085 B CN113128085 B CN 113128085B
- Authority
- CN
- China
- Prior art keywords
- aerodynamic
- vortex
- equation
- gust
- velocity
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 238000004458 analytical method Methods 0.000 title claims abstract description 35
- 230000008878 coupling Effects 0.000 claims abstract description 30
- 238000010168 coupling process Methods 0.000 claims abstract description 30
- 238000005859 coupling reaction Methods 0.000 claims abstract description 30
- 230000001133 acceleration Effects 0.000 claims abstract description 7
- 238000004088 simulation Methods 0.000 claims abstract description 7
- 239000013598 vector Substances 0.000 claims description 56
- 239000011159 matrix material Substances 0.000 claims description 47
- 238000004364 calculation method Methods 0.000 claims description 16
- 230000008859 change Effects 0.000 claims description 12
- 238000006073 displacement reaction Methods 0.000 claims description 10
- 230000004044 response Effects 0.000 claims description 7
- 230000006698 induction Effects 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims description 4
- 230000009471 action Effects 0.000 claims description 2
- 238000013016 damping Methods 0.000 claims description 2
- 238000005259 measurement Methods 0.000 claims description 2
- 230000008569 process Effects 0.000 claims description 2
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 230000005284 excitation Effects 0.000 abstract description 2
- 238000006243 chemical reaction Methods 0.000 abstract 1
- 238000013461 design Methods 0.000 description 6
- 238000012892 rational function Methods 0.000 description 5
- 238000011160 research Methods 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000000116 mitigating effect Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Fluid Mechanics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于状态空间形式涡格法的阵风减缓分析方法,依据给定模型建立气动网格,根据状态空间形式涡格法建立模型附着涡涡强、尾涡涡强与机翼表面诱导速度的关系,插值得到有限元模态及舵面模态广义坐标与机翼表面诱导速度的关系后形成气动弹性耦合方程,气动弹性耦合方程同样为状态空间形式,给定阵风模型作为外激励,可进行气动弹性阵风分析;在气动弹性耦合方程的基础上以机翼某点加速度作为系统输出,经过控制器后形成舵面偏转输入信号,搭建闭环控制系统,进行气动弹性阵风减缓仿真分析。本发明的方法能够提高阵风分析精度,避免频域时域转换,提高分析效率。
Description
技术领域
本发明属于结构动力学和气动弹性力学分析技术领域,尤其涉及一种基于状态空间形式涡格法的阵风减缓分析方法。
背景技术
气动弹性力学作为应用力学的分支,主要研究气动力、弹性力与惯性力之间耦合的问题。在气动力作用下弹性结构发生振动与变形,结构弹性运动反过来又会导致气动力的大小及分布发生变化,这种相互作用会带来多种多样的气动弹性现象,包括气动弹性变形、颤振、阵风响应等。飞行器在飞行中遭遇阵风扰动后机体结构会产生颠簸,干扰驾驶员操作,影响乘坐舒适性,增加气动弹性载荷,缩短结构疲劳寿命,严重威胁飞行安全。随着现代飞行器设计要求的逐渐提高,复合材料使用率及结构柔性的增大使得飞行器对阵风扰动的敏感性进一步提高,阵风扰动对于飞行器的影响更加复杂。阵风减缓分析研究对于提高现代飞行器的飞行性能,保障飞行安全具有重要意义。
利用势流理论面元法进行气动力计算是飞行器气动载荷分析的重要方法。面元法在求解精度上相较二维气动力建模方法提高很多且计算效率高,满足理论研究及工程应用需求。偶极子格网法和涡格法是两类常用的气动力计算方法,能够分别在频域和时域上给出气动力计算结果,在气动弹性分析特别是阵风问题分析上应用广泛。
基于偶极子格网法的阵风分析方法是飞行器气动弹性阵风分析的传统方法吗,由于偶极子格网法本身仅能提供频域气动力计算结果,为得到便于系统仿真及控制器设计的状态空间形式动力学方程,需要在偶极子格网法求得的广义非定常气动力影响系数矩阵的基础上进行有理函数拟合,得到气动力时域表达形式。即该方法需要先建立频域气动力模型,再进行频域到时域的拟合,计算较繁琐;同时,有理函数拟合的准确性取决于气动力滞后根的选取,选取依赖经验及反复试验,容易出现较大误差。
涡格法相较于偶极子格网法具有优势。涡格法是时域三维气动力分析方法,不需要简谐振动假设,可以直接计算任意运动下的非定常气动力。传统基于涡格法进行弹性飞行器阵风响应分析主要基于松耦合迭代算法进行结构与气动模型的耦合求解。在该分析框架下针对飞行器全机或部件进行阵风响应分析,计算结果精度令人满意。但随着分析模型复杂性的增加,基于涡格法的松耦合迭代求解方法计算量大,分析效率大幅度降低。同时,松耦合迭代算法也使得分析方法难以应用于控制器设计,无法进行阵风减缓研究。综上,目前气动弹性阵风分析中主流的基于(1)偶极子格网法- 有理函数拟合方法及(2)涡格法-松耦合迭代方法的两类方法均有不足。
状态空间形式的涡格法可以将涡格法气动力模型写为状态空间形式同时不损失精度,Breuker等人在其已发表论文中给出了状态空间形式涡格法的表达形式及基于简单梁模型与状态空间涡格法的气动弹性方程,并未将其应用于基于有限元方法的结构模型,使用范围受到很大限制。同时,已发表文献及专利中未有涉及状态空间形式涡格法能够将气动模型写为状态空间形式,在控制仿真及设计上具有很大优势。
发明内容
为了克服传统气动弹性阵风分析中(1)偶极子格网法-有理函数拟合方法需要进行有理函数拟合导致精度下降、步骤繁琐及(2)涡格法-松耦合迭代方法计算效率低,不便于控制系统设计的不足。已有的基于状态空间形式涡格法的阵风分析不能适用于有限元复杂模型同时未有扩展至阵风减缓控制问题,本发明基于状态空间形式的涡格法建立气动弹性阵风响应计算及气动弹性阵风减缓分析框架,提出一种高效准确的气动弹性阵风响应计算方法,建立气动弹性阵风减缓分析方法。本发明的具体技术方案如下:
一种基于状态空间形式涡格法的阵风减缓分析方法,包括以下步骤:
S1:计算初始化;
建立气动模型,在机翼中弧面沿弦向和展向划分四边形气动网格,气动网格包括在机翼表面的附着涡网格及沿来流方向拖出的尾涡网格;
S2:对步骤S1得到的气动网格建立状态空间形式的气动力方程;
S2-1:在气动网格每个1/4弦线处布置涡线段,四段等强度涡线首尾相连构成一个涡环;选择1/4弦线中点为力的作用点,选择3/4弦线中点为控制点;
S2-2:将气动模型中的涡分为三部分:机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡,设机翼表面附着涡强度列向量为Γb,机翼后缘第一排尾涡强度为Γw0,其他尾涡强度为Γwl,气动控制方程为:
KbΓb+Kw0Γw0+KwlΓwl=-w (1)
其中,w=(V∞+Vg-Vb)·n为机翼表面诱导速度,V∞为来流速度矢量,Vg为来流阵风扰动速度矢量及Vb为壁面运动速度矢量,n为控制点处的法向量列阵;Kb,Kw0,Kwl分别为附着涡、第一排尾涡及其他尾涡的诱导系数矩阵;
S2-3:机翼后缘的第一排尾涡在脱出的过程中保持涡强守恒,关系表达为:
S2-4:规定其他尾涡脱出后保持涡强不变,关系表达为:
S2-5:综合式(1)(2)(3)得到涡格法状态空间方程形式:
其中,Γw=[Γw0 Γwl]T,Aa,Ba为状态空间系数矩阵,仅与气动面几何形状及机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡三者的划分有关:
其中,I为单位阵,O为只包含0元素的零矩阵;
S2-6:给出网格上的气动力表达,气动网格上的压强差表达为:
其中,下标i,j表示沿展向第i个,沿弦向第j个网格中对应物理量矢量或矩阵的分量,Δpij为气动网格上压强差分量,Vl,ij为气动网格上当地来流速度分量,Γb,ij为气动网格上表面附着涡的涡强分量,ρ∞为来流大气密度,τ1,τ2分别为气动网格中控制点处沿当地速向及弦向方向的切向量;
假设气动面扰动速度远小于来流速度,即认为Vl=V∞,方程(7)中涡强的变化量由下述差分方程求解:
作用在气动网格上的气动力表示为:
Fij=ΔpijSijn (9)
其中,Δcij为气动网格沿弦向几何长度,Δbij为气动网格沿展向几何长度,n为控制点处的法向量列阵,Fij为气动网格上的气动力向量,Sij为气动网格面积;综合方程 (1)(7)(8)(9)得到气动力矢量F表达式:
其中,B1,B2,B3为系数矩阵;
S2-7:涡格法状态空间方程完整表达为:
状态空间形式涡格法气动力模型即方程(11)所示,给定输入量后即能够求解得到气动力数值,输入量中包含阵风扰动量,因此适用于计算阵风问题;
S3:建立结构有限元模型;
利用有限元方法建立结构模型,同时给出模态坐标下的结构动力学方程:
其中,Ms,Bs,Ks分别为广义质量阵、阻尼阵及刚度阵,qs为广义模态坐标,Fs为广义力;模态坐标与物理坐标的关系为:
us=Φqs (13)
其中,us为物理空间坐标,Φ为模态矩阵;
S4:利用三维样条插值理论得到结构-气动耦合关系;
步骤S2得到气动力方程,步骤S3得到结构动力学方程,气动弹性计算需要将气动与结构计算耦合在一起,利用三维样条插值理论耦合有限元模态广义自由度与气动面控制点法向量变化量,表征气动力方程与结构动力学方程之间的耦合关系;
气动面位移WA与结构位移WS关系为:
WA=EWS (14)
其中,E为位移插值矩阵,同样的,气动面速度VA插值能够由结构速度VS与映射关系求解:
VA=EVS (15)
由某一结构位移uS导致的气动面上任意点切向量τ的变化Δτ同样能够由映射关系求解:
Δτ=ET0uS (16)
其中,ET0为物理位移到切向量变化的插值矩阵,考虑物理位移与模态空间广义自由度关系(13),则有:
Δτ=ETqS (17)
其中,ET=ET0ΦuS为模态空间广义自由度到切向量变化的插值矩阵;
考虑结构变形后空间法向量与两个切向量的叉乘关系:
其中,n为控制点处的法向量列阵,n0,τ1,τ2分别为未变形状态下控制点处的法向量、气动网格中控制点沿当地速向的切向量及气动网格中控制点沿当地弦向方向的切向量,忽略二阶小量Δτ1×Δτ2,结合方程(17)得到变形后法向量与模态空间广义自由度关系:
n=n0+Δn=n0+ENqS (19)
其中,EN为模态空间广义自由度即结构模态坐标到法向量变化量的插值矩阵;
考虑结构变形后,诱导速度表达为:
忽略扰动速度与法向量增量相乘形成的高阶项,诱导速度表达为:
w=V∞n0+Vgn0-Vbn0+V∞Δn (21)
壁面运动速度Vb完全由弹性振动导致:
综合方程(19)(21)(22)得到耦合关系:
S5:给出气动弹性系统状态空间方程;
通过步骤S4的结构-气动耦合关系耦合步骤S2中的气动力方程及步骤S3中的结构动力学方程组装得到气动弹性系统状态空间方程;
改写方程(23)为:
其中,R1s=EΦ,R2s=V∞EN,R0为来流速度在初始构型下导致的诱导速度量,Rg为单位阵风速度导致的诱导速度系数矩阵,R1s为模态坐标变化单位量导致的诱导速度系数矩阵,R2s为模态坐标时间导数变化单位量导致的诱导速度系数矩阵,vg为空间流场中的阵风速度向量;
为便于引入控制面自由度,引入舵面偏转模态ΦC,舵面偏转模态与弹性振动模态处理方式一致,即将舵面偏转带来的气动力变化通过对诱导速度处法向量变化Δn 及运动速度Vb的影响体现,则考虑舵面偏转δ的诱导速度为:
结合方程(10),结构广义外力为:
其中,EF为力插值矩阵,综合方程(25)(26):
其中,K0,K1,K2,K3,K4,K5,K6,K7,K8,K9为对应的系数矩阵,结合涡格法状态空间方程(4),气动力控制方程为:
综合方程(12)(27)(28)得到气动弹性系统状态空间方程:
考虑结构信号输出,常用为传感器安置点加速度信号,则:
其中,Φr为传感器位置的模态矩阵;
S6:引入控制系统;
选择结构某位置处传感器测量的信号作为反馈信号,回传控制器,通过控制增益环节计算后转化为控制信号传给舵机,驱动控制面,控制系统方程为:
其中Kcon为控制增益矩阵;
方程(29)(30)(31)构成气动弹性阵风减缓控制系统,能够在Matlab软件Simulink模块中搭建实现。
本发明的有益效果在于:
1.传统气动弹性阵风减缓分析方法存在诸多不足。考虑到状态空间形式涡格法具有计算准确、状态空间形式便于控制器设计等优势,本发明对飞行器实际外形建立状态空间形式涡格法控制方程,采用有限元方法建立结构模型,使得本发明具有更好的模型适用性。
2.本发明给出了结构动力学方程与气动控制方程的耦合关系,将状态空间形式涡格法的尾涡涡强变量与有限元模型模态空间自由度耦合在一起,为气动弹性系统状态空间方程推导奠定了基础。
3.本发明给出了气动弹性系统状态空间方程,引入控制环节综合气动弹性系统状态空间方程得到气动弹性阵风减缓分析方程,将状态空间涡格法应用于气动弹性阵风减缓分析中,并适用于有限元方法建立的模型。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,通过参考附图会更加清楚的理解本发明的特征和优点,附图是示意性的而不应理解为对本发明进行任何限制,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,可以根据这些附图获得其他的附图。其中:
图1是本发明分析模型建立流程图;
图2是本发明控制思路;
图3是实施例的飞翼全机有限元模型示意图;
图4是实施例的拆分后机翼及机身有限元模型示意图;
图5是实施例的阵风减缓分析结果。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点,下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是,在不冲突的情况下,本发明的实施例及实施例中的特征可以相互组合。
在下面的描述中阐述了很多具体细节以便于充分理解本发明,但是,本发明还可以采用其他不同于在此描述的其他方式来实施,因此,本发明的保护范围并不受下面公开的具体实施例的限制。
如图1-2所示,本发明依据给定模型建立气动网格,根据状态空间形式涡格法建立模型附着涡涡强、尾涡涡强与机翼表面诱导速度的关系,插值得到有限元模态及舵面模态广义坐标与机翼表面诱导速度的关系后形成气动弹性耦合方程,气动弹性耦合方程同样为状态空间形式,给定阵风模型作为外激励,可进行气动弹性阵风分析;在气动弹性耦合方程的基础上以机翼某点加速度作为系统输出,经过控制器后形成舵面偏转输入信号,搭建闭环控制系统,进行气动弹性阵风减缓仿真分析。
为了方便理解本发明的上述技术方案,以下通过具体实施例对本发明的上述技术方案进行详细说明。
实施例1
采用一大展弦机翼模型,机翼模型主梁为翼尖到翼根线性增大的十字梁。主梁位置位于弦长的40%。主梁在MSC.Nastran中以CBEAM梁单元模拟,翼肋以大刚度梁单元模拟,质量特性以集中质量点模拟。模型部分参数如表1所示。
表1飞翼模型参数
参数名 | 参数值 |
机翼展长/m | 1.542 |
翼根弦长/m | 0.261 |
翼梢弦长/m | 0.069 |
扭转角/Deg | -2.0 |
展弦比 | 9.3 |
质量/kg | 2.8199 |
翼型 | 超临界翼型 |
S1:计算初始化;
建立气动模型,在机翼中弧面沿弦向划分5个网格,展向划分40个网格,拖出的尾涡网格4弦向划分50个网格,展向划分40个网格;机翼网格划分如图3所示。
S2:对步骤S1得到的气动网格建立状态空间形式的气动力方程;
涡格法状态空间方程完整表达为:
S3:建立结构有限元模型;
利用有限元方法建立结构模型,有限元模型如图4所示。有限元模型在MSC.Nastran软件中建立,计算中采用前6阶模态给出模态空间结构动力学方程,模态信息如表2所示。
表2模态信息表
S4:利用三维样条插值理论得到结构-气动耦合关系;
综合方程(19)(21)(22)得到耦合关系:
S5:给出气动弹性系统状态空间方程;
为便于引入控制面自由度,引入舵面偏转模态ΦC,面积为0.0000168m2,平均弦长为0.04m,据翼根距离0.8045m。
本实施例中采用“Sin”型离散阵风模型,数学表达式为:
其中,wg为阵风速度,wg0为阵风速度幅值,Lg为阵风尺度,t为仿真时间,xg为飞行器位置,x0为飞行器初始位置,一般可置零,V为来流速度。离散阵风指定为一维阵风,在展向保持均匀不变,阵风速度垂直于来流方向,考虑前述状态空间形式表达,阵风加载形式为:
vg=[0 0 wg(l)]T (36)
综合方程(12)(27)(28)得到气动弹性系统状态空间方程:
考虑结构信号输出,为加速度传感器信号,加速度传感器布置于距翼尖1/3翼展处,则:
其中,Φr为加速度传感器位置的模态矩阵。
S6:引入控制系统;
选择加速度传感器信号作为反馈信号,回传控制器,通过控制增益环节计算后转化为控制信号传给舵机,驱动控制面。本实施例采用经典PID控制律进行阵风减缓分析。
本实施例阵风减缓分析在Matlab软件Simulink模块中搭建实现。图5为阵风减缓分析结果,开环状态表示为引入控制系统的机翼阵风响应分析结果,闭环状态表示引入控制系统后的机翼阵风减缓分析结果,由图可知,本发明基于状态空间形式涡格法及结构有限元模型构建气动弹性系统阵风响应分析的方法,能够提高阵风分析精度,避免频域时域转换,提高分析效率。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (1)
1.一种基于状态空间形式涡格法的阵风减缓分析方法,其特征在于,包括以下步骤:
S1:计算初始化;
建立气动模型,在机翼中弧面沿弦向和展向划分四边形气动网格,气动网格包括在机翼表面的附着涡网格及沿来流方向拖出的尾涡网格;
S2:对步骤S1得到的气动网格建立状态空间形式的气动力方程;
S2-1:在气动网格每个1/4弦线处布置涡线段,四段等强度涡线首尾相连构成一个涡环;选择1/4弦线中点为力的作用点,选择3/4弦线中点为控制点;
S2-2:将气动模型中的涡分为三部分:机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡,设机翼表面附着涡强度列向量为Γb,机翼后缘第一排尾涡强度为Γw0,其他尾涡强度为Γwl,气动控制方程为:
KbΓb+Kw0Γw0+KwlΓwl=-w (1)
其中,w=(V∞+Vg-Vb)·n为机翼表面诱导速度,V∞为来流速度矢量,Vg为来流阵风扰动速度矢量及Vb为壁面运动速度矢量,n为控制点处的法向量列阵;Kb,Kw0,Kwl分别为附着涡、第一排尾涡及其他尾涡的诱导系数矩阵;
S2-3:机翼后缘的第一排尾涡在脱出的过程中保持涡强守恒,关系表达为:
S2-4:规定其他尾涡脱出后保持涡强不变,关系表达为:
S2-5:综合式(1)(2)(3)得到涡格法状态空间方程形式:
其中,Γw=[Γw0 Γwl]T,Aa,Ba为状态空间系数矩阵,仅与气动面几何形状及机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡三者的划分有关:
其中,I为单位阵,O为只包含0元素的零矩阵;
S2-6:给出网格上的气动力表达,气动网格上的压强差表达为:
其中,下标i,j表示沿展向第i个,沿弦向第j个网格中对应物理量矢量或矩阵的分量,Δpij为气动网格上压强差分量,Vl,ij为气动网格上当地来流速度分量,Γb,ij为气动网格上表面附着涡的涡强分量,ρ∞为来流大气密度,τ1,τ2分别为气动网格中控制点处沿当地速向及弦向方向的切向量;
假设气动面扰动速度远小于来流速度,即认为Vl=V∞,方程(7)中涡强的变化量由下述差分方程求解:
作用在气动网格上的气动力表示为:
Fij=ΔpijSijn (9)
其中,Δcij为气动网格沿弦向几何长度,Δbij为气动网格沿展向几何长度,n为控制点处的法向量列阵,Fij为气动网格上的气动力向量,Sij为气动网格面积;综合方程(1)(7)(8)(9)得到气动力矢量F表达式:
其中,B1,B2,B3为系数矩阵;
S2-7:涡格法状态空间方程完整表达为:
状态空间形式涡格法气动力模型即方程(11)所示,给定输入量后即能够求解得到气动力数值,输入量中包含阵风扰动量,因此适用于计算阵风问题;
S3:建立结构有限元模型;
利用有限元方法建立结构模型,同时给出模态坐标下的结构动力学方程:
其中,Ms,Bs,Ks分别为广义质量阵、阻尼阵及刚度阵,qs为广义模态坐标,Fs为广义力;模态坐标与物理坐标的关系为:
us=Φqs (13)
其中,us为物理空间坐标,Φ为模态矩阵;
S4:利用三维样条插值理论得到结构-气动耦合关系;
步骤S2得到气动力方程,步骤S3得到结构动力学方程,气动弹性计算需要将气动与结构计算耦合在一起,利用三维样条插值理论耦合有限元模态广义自由度与气动面控制点法向量变化量,表征气动力方程与结构动力学方程之间的耦合关系;
气动面位移WA与结构位移WS关系为:
WA=EWS (14)
其中,E为位移插值矩阵,同样的,气动面速度VA插值能够由结构速度VS与映射关系求解:
VA=EVS (15)
由某一结构位移uS导致的气动面上任意点切向量τ的变化Δτ同样能够由映射关系求解:
Δτ=ET0uS (16)
其中,ET0为物理位移到切向量变化的插值矩阵,考虑物理位移与模态空间广义自由度关系(13),则有:
Δτ=ETqS (17)
其中,ET=ET0ΦuS为模态空间广义自由度到切向量变化的插值矩阵;
考虑结构变形后空间法向量与两个切向量的叉乘关系:
其中,n为控制点处的法向量列阵,n0,τ1,τ2分别为未变形状态下控制点处的法向量、气动网格中控制点沿当地速向的切向量及气动网格中控制点沿当地弦向方向的切向量,忽略二阶小量Δτ1×Δτ2,结合方程(17)得到变形后法向量与模态空间广义自由度关系:
n=n0+Δn=n0+ENqS (19)
其中,EN为模态空间广义自由度即结构模态坐标到法向量变化量的插值矩阵;
考虑结构变形后,诱导速度表达为:
忽略扰动速度与法向量增量相乘形成的高阶项,诱导速度表达为:
w=V∞n0+Vgn0-Vbn0+V∞Δn (21)
壁面运动速度Vb完全由弹性振动导致:
综合方程(19)(21)(22)得到耦合关系:
S5:给出气动弹性系统状态空间方程;
通过步骤S4的结构-气动耦合关系耦合步骤S2中的气动力方程及步骤S3中的结构动力学方程组装得到气动弹性系统状态空间方程;
改写方程(23)为:
其中,R0=V∞n0,R1s=EΦ,R2s=V∞EN,R0为来流速度在初始构型下导致的诱导速度量,Rg为单位阵风速度导致的诱导速度系数矩阵,R1s为模态坐标变化单位量导致的诱导速度系数矩阵,R2s为模态坐标时间导数变化单位量导致的诱导速度系数矩阵,vg为空间流场中的阵风速度向量;
为便于引入控制面自由度,引入舵面偏转模态ΦC,舵面偏转模态与弹性振动模态处理方式一致,即将舵面偏转带来的气动力变化通过对诱导速度处法向量变化Δn及运动速度Vb的影响体现,则考虑舵面偏转δ的诱导速度为:
结合方程(10),结构广义外力为:
其中,EF为力插值矩阵,综合方程(25)(26):
其中,K0,K1,K2,K3,K4,K5,K6,K7,K8,K9为对应的系数矩阵,结合涡格法状态空间方程(4),气动力控制方程为:
综合方程(12)(27)(28)得到气动弹性系统状态空间方程:
考虑结构信号输出,常用为传感器安置点加速度信号,则:
其中,Φr为传感器位置的模态矩阵;
S6:引入控制系统;
选择结构某位置处传感器测量的信号作为反馈信号,回传控制器,通过控制增益环节计算后转化为控制信号传给舵机,驱动控制面,控制系统方程为:
其中Kcon为控制增益矩阵;
方程(29)(30)(31)构成气动弹性阵风减缓控制系统,能够在Matlab软件Simulink模块中搭建实现。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110297250.6A CN113128085B (zh) | 2021-03-19 | 2021-03-19 | 一种基于状态空间形式涡格法的阵风减缓分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110297250.6A CN113128085B (zh) | 2021-03-19 | 2021-03-19 | 一种基于状态空间形式涡格法的阵风减缓分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113128085A CN113128085A (zh) | 2021-07-16 |
CN113128085B true CN113128085B (zh) | 2022-10-14 |
Family
ID=76773642
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110297250.6A Active CN113128085B (zh) | 2021-03-19 | 2021-03-19 | 一种基于状态空间形式涡格法的阵风减缓分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113128085B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114818020B (zh) * | 2022-03-24 | 2024-08-13 | 北京航空航天大学 | 基于势流理论的阵风预测方法和装置 |
CN114564047B (zh) * | 2022-04-28 | 2022-08-16 | 北京航空航天大学 | 一种考虑气象条件的无人机等速飞行控制方法 |
CN114942595B (zh) * | 2022-07-25 | 2022-11-18 | 西安爱生技术集团有限公司 | 一种考虑降雨影响的无人机阵风响应建模和分析方法 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108052772A (zh) * | 2017-12-30 | 2018-05-18 | 北京航空航天大学 | 一种基于结构降阶模型的几何非线性静气动弹性分析方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20150377915A1 (en) * | 2013-02-07 | 2015-12-31 | King Abdullah University Of Science And Technology | Method and system for estimating and predicting airflow around air vehicles |
-
2021
- 2021-03-19 CN CN202110297250.6A patent/CN113128085B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108052772A (zh) * | 2017-12-30 | 2018-05-18 | 北京航空航天大学 | 一种基于结构降阶模型的几何非线性静气动弹性分析方法 |
Non-Patent Citations (5)
Title |
---|
基于曲面涡格法的柔性飞机静气动弹性分析;刘燚等;《工程力学》;20180228;第35卷(第02期);第249-255页 * |
大展弦比柔性机翼气动弹性风洞模型设计与试验验证;谢长川等;《工程力学》;20161130;第33卷(第11期);第249-256页 * |
弹性飞机阵风响应建模与减缓方案设计;吴志刚等;《中国科学:技术科学》;20110331;第41卷(第03期);第394-402页 * |
柔性飞机大变形曲面气动力计算及配平分析;刘燚等;《工程力学》;20151031;第32卷(第10期);第239-248页 * |
考虑几何非线性的复合材料机翼气动弹性分析;梁宇等;《力学季刊》;20190430(第04期);第60-68页 * |
Also Published As
Publication number | Publication date |
---|---|
CN113128085A (zh) | 2021-07-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113128085B (zh) | 一种基于状态空间形式涡格法的阵风减缓分析方法 | |
Murua et al. | Applications of the unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics | |
CN113111430B (zh) | 基于非线性气动力降阶的弹性飞机飞行动力学建模方法 | |
Murua et al. | Assessment of wake-tail interference effects on the dynamics of flexible aircraft | |
Potsdam et al. | Rotor airloads prediction using loose aerodynamic/structural coupling | |
CN108427322B (zh) | 一种大柔性飞行器基于在线辨识的建模方法 | |
Raveh | CFD-based models of aerodynamic gust response | |
CN105183996A (zh) | 面元修正与网格预先自适应计算方法 | |
Howcroft et al. | Aeroelastic modelling of highly flexible wings | |
CN112580241B (zh) | 一种基于结构降阶模型的非线性气动弹性动响应分析方法 | |
CN113868771B (zh) | 一种考虑结构和气动非线性的飞行动力学建模方法 | |
CN110287505B (zh) | 飞行器稳定性分析方法 | |
Murua | Flexible aircraft dynamics with a geometrically-nonlinear description of the unsteady aerodynamics | |
Preidikman et al. | Time-domain simulations of linear and nonlinear aeroelastic behavior | |
CN112558623B (zh) | 一种弹性高超声速飞行器多胞lpv系统建模方法 | |
Dussart et al. | Flight dynamic modelling and simulation of large flexible aircraft | |
Hang et al. | Analytical sensitivity analysis of flexible aircraft with the unsteady vortex-lattice aerodynamic theory | |
Ritter et al. | Comparison of nonlinear aeroelastic methods for maneuver simulation of very flexible aircraft | |
Fugate et al. | Aero-Structural Modeling of the Truss-Braced Wing Aircraft Using Potential Method with Correction Methods for Transonic Viscous Flow and Wing-Strut Interference Aerodynamics | |
Mauermann | Flexible aircraft modelling for flight loads analysis of wake vortex encounters | |
CN113536457B (zh) | 一种基于状态空间形式涡格法的气动力降阶方法 | |
van Schie et al. | Solver-independent aeroelastic coupling for large-scale multidisciplinary design optimization | |
Maza et al. | A cost-effective algorithm for the co-simulation of unsteady and non-linear aeroelastic phenomena | |
Ricci et al. | Active control of three-surface aeroelastic model | |
Nguyen et al. | Multi-point jig twist optimization of mach 0.745 transonic truss-braced wing aircraft and high-fidelity cfd validation |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |