CN113536457A - 一种基于状态空间形式涡格法的气动力降阶方法 - Google Patents

一种基于状态空间形式涡格法的气动力降阶方法 Download PDF

Info

Publication number
CN113536457A
CN113536457A CN202110746621.4A CN202110746621A CN113536457A CN 113536457 A CN113536457 A CN 113536457A CN 202110746621 A CN202110746621 A CN 202110746621A CN 113536457 A CN113536457 A CN 113536457A
Authority
CN
China
Prior art keywords
aerodynamic
vortex
equation
representing
state space
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.)
Granted
Application number
CN202110746621.4A
Other languages
English (en)
Other versions
CN113536457B (zh
Inventor
安朝
谢长川
杨澜
杨超
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202110746621.4A priority Critical patent/CN113536457B/zh
Publication of CN113536457A publication Critical patent/CN113536457A/zh
Application granted granted Critical
Publication of CN113536457B publication Critical patent/CN113536457B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • 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)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Automation & Control Theory (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于状态空间形式涡格法的气动力降阶方法，在机翼中弧面建立气动网格，根据状态空间形式的涡格法建立机翼表面附着涡强度、尾涡强度与机翼表面诱导速度的关系，形成气动力状态空间方程；在气动力状态空间方程的基础上，给出给定工况下一定运动周期内的尾涡强度时域变化过程数据，将尾涡强度时域变化过程数据作为POD降阶训练样本，以尾涡强度作为POD降阶系统状态量，求解POD模态，得到POD气动力降阶方程；再根据非定常伯努利方程得到机翼表面气动力分布结果。上述方法与基于CFD的气动力降阶方法相比，计算量低，能够更快速生成大量样本，快速计算气动力，并且在建模方面更加简单，适用于相关分析研究的快速应用。

Description

一种基于状态空间形式涡格法的气动力降阶方法
技术领域
本发明涉及空气动力学分析技术领域，尤其涉及一种基于状态空间形式涡格法的气动力降阶方法。
背景技术
在飞行器设计问题中，许多分析例如气动弹性阵风响应分析、气动弹性优化设计等所涉及的模型复杂度高，其计算耗费往往主要是由于气动力计算所带来的，针对单点状态的数值模拟方法很难直接应用于控制模型综合、多变量优化、稳定性预测和实时仿真等多学科设计领域。同时，尽管高精度数值模拟可以提供离散化流场变量的详尽时空信息，但如果缺乏其他辅助工具和分析方法，单靠数值模拟提供的高阶模型和海量数据并不足以让人们深入理解和描述系统的复杂动力学行为。建立气动力降阶模型成为降低问题求解规模、深入理解系统动力学行为的主要研究方向，对于飞行器设计问题分析具有重要意义。
构造降阶模型的目标主要有以下两点，一是以远少于原数值模型的阶数和计算耗费提供系统主要动力学特征较精确的数学描述；二是为研究者解释系统动力学特征提供工具。在保留全阶高精度模型的可信度和高保真度的同时，减小计算量，且能够方便地与其他学科模型进行耦合以用于多学科耦合分析与优化设计。通过降阶技术在单点仿真的高可信度数值模型和多学科耦合系统仿真设计之间架起桥梁。
目前，非定常气动力的降阶技术主要分为两类，一类是基于数据驱动或信号的方法，另一类是基于流场特征模态的方法。基于信号的降阶方法是利用流体系统的输入和输出之间的关系，建立系统的一个低阶传递函数或状态空间模型，以代替原始的全阶模型。从本质上讲是一种系统辨识方法，主要以Volterra级数模型及ARMA模型为代表。各类非线性气动载荷的代理模型也可归入此类。基于流场特征模态的降阶模型就是用一组低维流场变量的特征模态来描述总的流场运动，再将整个全阶模型通过Galerkin方法或Krylov方法投影到特征模态空间。本征正交分解(Proper Orthogonal Decomposition,POD)方法作为该类降阶模型的代表，近年来受到研究人员的关注。POD方法可以客观的、不带偏见的从实验或数值分析中获取在均方根意义下的POD正交模态，用来反映对象数据库中数据集合的特征。将原系统方程在POD基上进行Galerkin投影，就可以获得一组低维模型。由于可以事先获知各阶POD模态及其略去的高阶模态所占有的能量比重，因此应用POD方法可以比较客观的构造和分析原系统，适用于非定常气动力建模与分析。
降阶模型能够有效降低系统自由度，提高计算效率。现有的气动力降阶模型通常基于计算流体力学(Computational Fluid Dynamics,CFD)计算方法开展研究。近年来基于CFD数据的POD降阶模型得到广泛发展，而针对势流理论面元法的POD降阶模型研究相对较少。势流理论面元法是飞行器气动载荷分析的重要方法，其在生成降阶模型训练样本时所花费的时间少，气动力网格建模上更加简单，适用于快速应用，并且，在求解精度上相较二维气动力建模方法提高很多，虽然其计算效率比CFD方法高，但处理多变量优化、稳定性预测和实时仿真等多学科设计问题时同样会因为计算量过大而力不从心。涡格法是势流理论面元法中的一种重要方法，作为时域三维气动力分析方法，不需要简谐振动假设，可以直接计算任意运动下的非定常气动力。状态空间形式的涡格法可以将涡格法气动力模型写为状态空间形式同时不损失精度，但目前尚无基于状态空间形式涡格法的POD降阶模型的研究。
发明内容
有鉴于此，本发明提供了一种基于状态空间形式涡格法的气动力降阶方法，用以克服传统飞行器设计中气动分析模型采用CFD方法及涡格法在处理计算量较大的多学科分析问题时分析效率低的问题。
本发明提供的一种基于状态空间形式涡格法的气动力降阶方法，包括如下步骤：
S1：在机翼的中弧面沿弦向和展向划分出若干个四边形的气动网格，所述气动网格包括机翼表面的附着涡网格和沿来流方向拖出的尾涡网格；
S2：基于划分的气动网格，以尾涡强度为状态变量，机翼表面诱导速度为输入量，气动力为输出量，计算系数矩阵，得到涡格法状态空间方程；
S3：在所述涡格法状态空间方程的基础上，给定尾涡强度初值和预先选定的来流扰动形式，求解得到该来流扰动形式下尾涡强度在预设时间内的时域变化过程数据；
S4：利用POD降阶方法，基于所述涡格法状态空间方程和预先选定的来流扰动形式下尾涡强度在预设时间内的时域变化过程数据，建立尾涡强度与广义坐标的降阶关系及POD气动力降阶方程；
S5：基于所述POD气动力降阶方程，给定任意一种来流扰动形式，求解得到任意来流扰动形式下广义坐标的时域变化过程，结合尾涡强度与广义坐标的降阶关系，恢复出任意来流扰动形式下尾涡强度的时域变化过程；
S6：将所述广义坐标的时域变化过程代入所述涡格法状态空间方程的输出方程中，得到时域气动力。
在一种可能的实现方式中，在本发明提供的上述基于状态空间形式涡格法的气动力降阶方法中，步骤S2，具体包括：
在每个气动网格的四条边上布置涡线段，每个气动网格上布置的四条涡线段的强度相等且四条涡线段首尾相连构成一个涡环；选择每个气动网格的1/4弦线的中点为气动力的作用点，选择每个气动网格的3/4弦线的中点为气动力的控制点；
将气动网格中的涡分为机翼表面附着涡、机翼后缘第一排尾涡以及其他尾涡三部分，设机翼表面附着涡的强度列向量为Γb，机翼后缘第一排尾涡的强度列向量为Γw0，其他尾涡的强度列向量为Γwl，则气动控制方程为：
KbΓb+Kw0Γw0+KwlΓwl＝-w (1)
其中，w＝(V∞+Vg)·n，表示机翼表面的诱导速度；V∞表示来流的速度矢量，Vg表示来流扰动的速度矢量，n表示控制点处的法向量列阵；Kb表示机翼表面附着涡的诱导系数矩阵，Kw0表示机翼后缘第一排尾涡的诱导系数矩阵，Kwl表示其他尾涡的诱导系数矩阵；
机翼后缘第一排尾涡在脱出的过程中保持强度守恒，表达式为：
Figure BDA0003144536110000041
其中，
Figure BDA0003144536110000042
表示机翼后缘第一排尾涡的强度关于时间的导数；Δt表示时间步长；C1表示保证机翼后缘第一排尾涡与机翼表面附着涡对应关系正确的系数矩阵，包含0和1两种元素；
规定其他尾涡在脱出后保持强度不变，表达式为：
Figure BDA0003144536110000043
其中，
Figure BDA0003144536110000044
表示其他尾涡的强度关于时间的导数；C2和C3为表征其他尾涡与机翼后缘第一排尾涡位置对应关系的常数提取矩阵，包含0和1两种元素；I表示单位矩阵；
综合方程(1)、(2)和(3)得到涡格法状态空间方程中的状态方程为：
Figure BDA0003144536110000045
其中，Γw＝[Γw0 Γwl]T，表示尾涡强度；Aa和Ba表示状态空间系数矩阵，仅与气动面几何形状以及机翼表面附着涡、机翼后缘第一排尾涡和其他尾涡的划分有关，表达式为：
Figure BDA0003144536110000046
Figure BDA0003144536110000051
其中，O表示只包含0元素的零矩阵；
气动网格上的压强差表示为：
Figure BDA0003144536110000052
其中，下标ij表示沿展向第i个、沿弦向第j个气动网格中对应的物理量，Δpij表示沿展向第i个、沿弦向第j个气动网格中的压强差，Vl,ij表示沿展向第i个、沿弦向第j个气动网格中的当地来流速度，Γb,ij表示沿展向第i个、沿弦向第j个气动网格中机翼表面附着涡的涡强，ρ∞表示来流大气密度，τ1表示气动力作用点沿当地速度方向的切向量，τ2表示气动力作用点沿弦向的切向量，t表示仿真时间；假设Vl＝V∞，Vl表示当地速度方向矢量，V∞表示来流速度方向矢量；
方程(7)中沿展向第i个、沿弦向第j个气动网格中机翼表面附着涡的涡强的变化量由下述差分方程求解得到：
Figure BDA0003144536110000053
Figure BDA0003144536110000054
其中，Γb,i-1,j表示沿展向第i-1个、沿弦向第j个气动网格中机翼表面附着涡的涡强，Γb,i,j-1表示沿展向第i个、沿弦向第j-1个气动网格中机翼表面附着涡的涡强，Δcij表示气动网格沿弦向的几何长度，Δbij表示气动网格沿展向的几何长度；
作用在气动网格上的气动力表示为：
Fij＝ΔpijSijn (10)
其中，Sij表示气动网格的面积；
综合方程(1)、(7)、(8)和(9)得到涡格法状态空间方程中的输出方程为
Figure BDA0003144536110000061
其中，F表示气动力矢量，B1表示机翼表面诱导速度w与气动力矢量F之间关系的系数矩阵，B2表示机翼表面诱导速度w关于时间的导数
Figure BDA0003144536110000062
与气动力矢量F之间关系的系数矩阵，B3为表示尾涡强度Γw与气动力矢量F之间关系的系数矩阵；方程(4)和方程(11)构成涡格法状态空间方程。
在一种可能的实现方式中，在本发明提供的上述基于状态空间形式涡格法的气动力降阶方法中，步骤S3中，给定离散阵风作为来流扰动形式：
wt＝(V∞+Vgt)·n (12)
其中，
Figure BDA0003144536110000063
表示离散阵风的扰动速度；wgm表示离散阵风的幅值，ω表示离散阵风的频率；
当给定离散阵风作为扰动形式时，w＝wt，尾涡强度Γw＝Γwt。
在一种可能的实现方式中，在本发明提供的上述基于状态空间形式涡格法的气动力降阶方法中，步骤S4，具体包括：
将尾涡强度在预设时间内的时域变化过程数据Γwt在时间域内的离散点[t1 t2… ts … tm]下的数据作为快照数据，整理成如下快照矩阵：
Q＝[Γwt(t1) Γwt(t2) … Γwt(ts) … Γwt(tm)] (13)
其中，Γwt(ts)表示尾涡强度在ts时刻下的数值向量，s＝1,2,…,m，m表示时间离散点总数；
将快照矩阵中的数据列向量取均值：
Figure BDA0003144536110000064
将快照矩阵中的每一列元素减去均值
Figure BDA0003144536110000065
得到新的快照矩阵：
Figure BDA0003144536110000066
定义相关矩阵：
Figure BDA0003144536110000067
计算相关矩阵的非零特征值及特征向量：
Figure BDA0003144536110000071
其中，g＝1,2,…,l，l表示保留的特征向量总数，l≤m，λg表示第g阶特征值，λ1≥λ2≥…λl＞0，
Figure BDA0003144536110000072
表示第g阶特征向量；
构造最优POD基底φg如下：
Figure BDA0003144536110000073
其中，φg是一组标准正交基；
选择前k个POD基底构成k维POD基向量矩阵：
Φ＝[φ1 φ2 … φk] (19)
其中，k<l；基于涡格法状态空间方程，得到尾涡强度与广义坐标的降阶关系：
Γw＝Φp (20)
其中，p表示广义坐标列向量；将方程(20)代入涡格法状态空间方程中的状态方程(4)后左乘ΦT，得到POD气动力降阶方程如下：
Figure BDA0003144536110000074
其中，A＝ΦTAaΦ，B＝ΦTBa，A和B为POD气动力降阶方程的状态空间系数矩阵，Aa和Ba为涡格法状态空间方程中的状态方程的系数矩阵。
在一种可能的实现方式中，在本发明提供的上述基于状态空间形式涡格法的气动力降阶方法中，步骤S6中，时域气动力的表达式为：
Figure BDA0003144536110000075
本发明提供的上述基于状态空间形式涡格法的气动力降阶方法，在机翼中弧面建立气动网格，根据状态空间形式的涡格法建立机翼表面附着涡强度、尾涡强度与机翼表面诱导速度的关系，形成气动力状态空间方程；在气动力状态空间方程的基础上，给出给定工况下一定运动周期内的尾涡强度时域变化过程数据，将尾涡强度时域变化过程数据作为POD降阶训练样本，以尾涡强度作为POD降阶系统状态量，求解POD模态，得到POD气动力降阶方程；进一步根据非定常伯努利方程得到机翼表面气动力布结果。上述方法与基于CFD的气动力降阶方法相比，计算量低，能够更快速生成大量样本，快速计算气动力，并且，在建模方面更加简单，适用于相关分析研究的快速应用。
附图说明
图1为本发明提供的一种基于状态空间形式涡格法的气动力降阶方法的流程图；
图2为本发明实施例1中机翼气动面模型示意图；
图3为本发明实施例1中机翼POD降阶模型特征值分布图；
图4为本发明实施例1中机翼POD降阶模型气动力计算结果图。
具体实施方式
下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限制本发明。
传统基于CFD方法或涡格法进行气动力分析计算，虽然能够得到精确的气动力计算结果，但计算效率不高，尤其是在处理多学科分析等需要反复迭代计算的问题时，计算效率成为限制问题分析的主要因素。
基于此，本发明提供的一种基于状态空间形式涡格法的气动力降阶方法，如图1所示，包括如下步骤：
第一步：计算初始化。
在机翼的中弧面沿弦向和展向划分出若干个四边形的气动网格，气动网格包括机翼表面的附着涡网格和沿来流方向拖出的尾涡网格。
第二步：对第一步得到的气动网格建立状态空间形式的气动力方程。
基于划分的气动网格，以尾涡强度为状态变量，机翼表面诱导速度为输入量，气动力为输出量，计算系数矩阵，得到涡格法状态空间方程。
(1)在每个气动网格的四条边上布置涡线段，每个气动网格上布置的四条涡线段的强度相等且四条涡线段首尾相连构成一个涡环；选择每个气动网格的1/4弦线的中点为气动力的作用点，选择每个气动网格的3/4弦线的中点为气动力的控制点；
(2)将气动网格中的涡分为机翼表面附着涡、机翼后缘第一排尾涡以及其他尾涡三部分，设机翼表面附着涡的强度列向量为Γb，机翼后缘第一排尾涡的强度列向量为Γw0，其他尾涡的强度列向量为Γwl，则气动控制方程可以为：
KbΓb+Kw0Γw0+KwlΓwl＝-w (1)
其中，w＝(V∞+Vg)·n，表示机翼表面的诱导速度；V∞表示来流的速度矢量，Vg表示来流扰动的速度矢量，n表示控制点处的法向量列阵；Kb表示机翼表面附着涡的诱导系数矩阵，Kw0表示机翼后缘第一排尾涡的诱导系数矩阵，Kwl表示其他尾涡的诱导系数矩阵；
(3)机翼后缘第一排尾涡在脱出的过程中保持强度守恒，表达式可以为：
Figure BDA0003144536110000091
其中，
Figure BDA0003144536110000092
表示机翼后缘第一排尾涡的强度关于时间的导数；Δt表示时间步长；C1表示保证机翼后缘第一排尾涡与机翼表面附着涡对应关系正确的系数矩阵，包含0和1两种元素；
(4)规定其他尾涡在脱出后保持强度不变，表达式可以为：
Figure BDA0003144536110000093
其中，
Figure BDA0003144536110000094
表示其他尾涡的强度关于时间的导数；C2和C3为表征其他尾涡与机翼后缘第一排尾涡位置对应关系的常数提取矩阵，包含0和1两种元素；I表示单位矩阵；
(5)综合方程(1)、(2)和(3)得到涡格法状态空间方程中的状态方程可以为：
Figure BDA0003144536110000101
其中，Γw＝[Γw0 Γwl]T，表示尾涡强度；Aa和Ba表示状态空间系数矩阵，仅与气动面几何形状以及机翼表面附着涡、机翼后缘第一排尾涡和其他尾涡的划分有关，表达式可以为：
Figure BDA0003144536110000102
Figure BDA0003144536110000103
其中，O表示只包含0元素的零矩阵；
(6)给出气动网格上的气动力表达。气动网格上的压强差可以表示为：
Figure BDA0003144536110000104
其中，下标ij表示沿展向第i个、沿弦向第j个气动网格中对应的物理量，Δpij表示沿展向第i个、沿弦向第j个气动网格中的压强差，Vl,ij表示沿展向第i个、沿弦向第j个气动网格中的当地来流速度，Γb,ij表示沿展向第i个、沿弦向第j个气动网格中机翼表面附着涡的涡强，ρ∞表示来流大气密度，τ1表示气动力作用点沿当地速度方向的切向量，τ2表示气动力作用点沿弦向的切向量，t表示仿真时间；假设气动面扰动速度远小于来流速度，可认为Vl＝V∞，Vl表示当地速度方向矢量，V∞表示来流速度方向矢量；
方程(7)中沿展向第i个、沿弦向第j个气动网格中机翼表面附着涡的涡强的变化量可以由下述差分方程求解得到：
Figure BDA0003144536110000105
Figure BDA0003144536110000106
其中，Γb,i-1,j表示沿展向第i-1个、沿弦向第j个气动网格中机翼表面附着涡的涡强，Γb,i,j-1表示沿展向第i个、沿弦向第j-1个气动网格中机翼表面附着涡的涡强，Δcij表示气动网格沿弦向的几何长度，Δbij表示气动网格沿展向的几何长度；
作用在气动网格上的气动力可以表示为：
Fij＝ΔpijSijn (10)
其中，Sij表示气动网格的面积；
综合方程(1)、(7)、(8)和(9)可以得到涡格法状态空间方程中的输出方程为
Figure BDA0003144536110000111
其中，F表示气动力矢量，B1表示机翼表面诱导速度w与气动力矢量F之间关系的系数矩阵，B2表示机翼表面诱导速度w关于时间的导数
Figure BDA0003144536110000112
与气动力矢量F之间关系的系数矩阵，B3为表示尾涡强度Γw与气动力矢量F之间关系的系数矩阵；方程(4)和方程(11)构成涡格法状态空间方程。
第三步：计算一定运动周期内的尾涡强度变化过程。
在第二步得到的涡格法状态空间方程的基础上，给定尾涡强度初值Γw0和预先选定的来流扰动形式wt，可以求解得到该来流扰动形式下尾涡强度在预设时间0～tf内的时域变化过程数据Γwt。通常对于飞行器设计问题，可以给定指定离散阵风作为扰动形式，得到一定运动周期内的尾涡强度变化过程。利用阵风扰动形式计算降阶所需样本，更加贴合飞行器分析的实际工况要求。
给定离散阵风作为来流扰动形式：
wt＝(V∞+Vgt)·n (12)
其中，
Figure BDA0003144536110000113
表示离散阵风的扰动速度；wgm表示离散阵风的幅值，ω表示离散阵风的频率；当给定离散阵风作为扰动形式时，w＝wt，尾涡强度Γw＝Γwt。
第四步：建立POD降阶模型。
第二步得到涡格法状态空间方程，第三步得到一定运动周期内的尾涡强度时域变化过程数据，即可利用POD降阶方法，基于涡格法状态空间方程和预先选定的来流扰动形式下尾涡强度在预设时间内的时域变化过程数据，建立尾涡强度与广义坐标的降阶关系及POD气动力降阶方程。
(1)取快照(Snapshot)数据。
第三步得到的一定时间0～tf内尾涡强度的时域变化过程数据为Γwt，将其在时间域内的离散点[t1 t2 … ts … tm]下的数据作为快照数据，整理成如下快照矩阵：
Q＝[Γwt(t1) Γwt(t2) … Γwt(ts) … Γwt(tm)] (13)
其中，Γwt(ts)表示尾涡强度在ts时刻下的数值向量，s＝1,2,…,m，m表示时间离散点总数；
(2)求POD降阶基底。
将快照矩阵中的数据列向量取均值：
Figure BDA0003144536110000121
将快照矩阵中的每一列元素减去均值
Figure BDA0003144536110000122
得到新的快照矩阵：
Figure BDA0003144536110000123
定义相关矩阵：
Figure BDA0003144536110000124
计算相关矩阵的非零特征值及特征向量：
Figure BDA0003144536110000125
其中，g＝1,2,…,l，l表示保留的特征向量总数，l≤m，λg表示第g阶特征值，λ1≥λ2≥…λl＞0，
Figure BDA0003144536110000126
表示第g阶特征向量；
构造最优POD基底φg如下：
Figure BDA0003144536110000131
其中，φg是一组标准正交基；
选择前k个POD基底构成k维POD基向量矩阵：
Φ＝[φ1 φ2 … φk] (19)
其中，k<l；
(3)给出降阶方程。
根据涡格法状态空间方程和POD基向量矩阵Φ，可以得到尾涡强度与广义坐标的降阶关系：
Γw＝Φp (20)
其中，p表示广义坐标列向量；将方程(20)代入涡格法状态空间方程中的状态方程(4)后左乘ΦT，可以得到气动力POD降阶方程如下：
Figure BDA0003144536110000132
其中，A＝ΦTAaΦ，B＝ΦTBa，A和B为POD气动力降阶方程的状态空间系数矩阵，Aa和Ba为涡格法状态空间方程中的状态方程的系数矩阵。根据方程(21)可以在不同工况条件下给出气动力降阶计算结果。
第五步：计算气动力。
基于步骤四得到的POD气动力降阶方程(21)，给定扰动输入(即任意一种来流扰动形式)后，可以求解得到任意来流扰动形式下广义坐标p的时域变化过程，结合尾涡强度与广义坐标的降阶关系，可以恢复出任意来流扰动形式下尾涡强度Γw的时域变化过程。
将广义坐标的时域变化过程代入涡格法状态空间方程的输出方程(11)中，可以得到时域气动力，表达式为：
Figure BDA0003144536110000133
方程(21)和(22)构成一种基于状态空间形式涡格法的气动力降阶模型，能够高效进行气动力分析计算。
下面通过一个具体的实施例对本发明提供的一种基于状态空间形式涡格法的气动力降阶方法的具体实施进行详细说明。
实施例1：
采用一大展弦机翼模型，该机翼模型的主梁为从翼尖到翼根线性增大的十字梁。主梁位置位于弦长的40％处。主梁在MSC.Nastran中以CBEAM梁单元模拟，翼肋以大刚度梁单元模拟，质量特性以集中质量点模拟。机翼模型的部分参数如表1所示。
表1机翼模型参数
参数名称 参数值
机翼展长/m 1.542
翼根弦长/m 0.261
翼梢弦长/m 0.069
扭转角/Deg -2.0
展弦比 9.3
质量/kg 2.8199
翼型 超临界翼型
第一步：计算初始化。
在机翼的中弧面沿弦向划分5个网格，展向划分40个网格，拖出的尾涡网格4弦向划分50个网格，展向划分40个网格；机翼网格划分如图2所示，图2中坐标系的X轴沿来流方向，Y轴沿翼展方向，满足右手定则。
第二步：对第一步得到的气动网格建立状态空间形式的气动力方程。
涡格法状态空间方程可表达为：
Figure BDA0003144536110000141
其中，尾涡强度Γw＝[Γw0 Γwl]T，Aa和Ba表示状态空间系数矩阵，w表示包含扰动在内的机翼表面诱导速度。
气动力矢量F的表达式为：
Figure BDA0003144536110000151
其中，B1,B2,B3为系数矩阵。
第三步：计算一定运动周期内的尾涡强度变化过程。
在第二步得到的涡格法状态空间方程(23)的基础上，给定尾涡强度初值Γw0＝0，来流风速31m/s，阵风频率ω＝4Hz，阵风幅值wgm＝0.505m，仿真时间tf＝10s，时间步长0.0025s。计算得到一定运动周期(0s～10s)内的尾涡涡强Γwt。
第四步：建立POD降阶模型。
第三步得到了0s～10s内的尾涡强度时域变化过程数据为Γwt，时间步长为0.0025s，因此10s内共有4000个数据点，取全部4000个数据离散点下的数据作为快照数据，整理成快照矩阵：
Q＝[Γwt(t1) Γwt(t2) … Γwt(ti) … Γwt(t4000)] (25)
其中，Γwt(ti)表示尾涡强度在第i个时间离散点的数值向量。
将快照矩阵中的数据列向量取均值：
Figure BDA0003144536110000152
将快照矩阵中的每一列元素减去均值
Figure BDA0003144536110000153
得到新的快照矩阵：
Figure BDA0003144536110000154
定义相关矩阵并计算相关矩阵的非零特征值及特征向量，按照方程(18)构造最优POD基底。前4阶基底特征值之和已经占到所有特征值之和的99.99％，特征值分布如图3所示。前4阶基底能够满足气动力分析计算要求，选择前4个POD基构成4维POD基向量矩阵：
Φ＝[φ1 φ2 φ3 φ4] (28)
尾涡强度与广义坐标的降阶关系可表达为：
Γw＝Φp (29)
气动力POD降阶方程为：
Figure BDA0003144536110000155
第五步：计算气动力。
气动力计算表达式为：
Figure BDA0003144536110000161
给出第三步所使用的初值及扰动，气动力在0～4s内的计算结果如图4所示。建立的气动力降阶模型能够高效进行气动力分析计算。
本发明提供的上述基于状态空间形式涡格法的气动力降阶方法，在机翼中弧面建立气动网格，根据状态空间形式的涡格法建立机翼表面附着涡强度、尾涡强度与机翼表面诱导速度的关系，形成气动力状态空间方程；在气动力状态空间方程的基础上，给出给定工况下一定运动周期内的尾涡强度时域变化过程数据，将尾涡强度时域变化过程数据作为POD降阶训练样本，以尾涡强度作为POD降阶系统状态量，求解POD模态，得到POD气动力降阶方程；进一步根据非定常伯努利方程得到机翼表面气动力分布结果。上述方法与基于CFD的气动力降阶方法相比，计算量低，能够更快速生成大量样本，快速计算气动力，并且，在建模方面更加简单，适用于相关分析研究的快速应用。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.一种基于状态空间形式涡格法的气动力降阶方法，其特征在于，包括如下步骤：
S1：在机翼的中弧面沿弦向和展向划分出若干个四边形的气动网格，所述气动网格包括机翼表面的附着涡网格和沿来流方向拖出的尾涡网格；
S2：基于划分的气动网格，以尾涡强度为状态变量，机翼表面诱导速度为输入量，气动力为输出量，计算系数矩阵，得到涡格法状态空间方程；
S3：在所述涡格法状态空间方程的基础上，给定尾涡强度初值和预先选定的来流扰动形式，求解得到该来流扰动形式下尾涡强度在预设时间内的时域变化过程数据；
S4：利用POD降阶方法，基于所述涡格法状态空间方程和预先选定的来流扰动形式下尾涡强度在预设时间内的时域变化过程数据，建立尾涡强度与广义坐标的降阶关系及POD气动力降阶方程；
S5：基于所述POD气动力降阶方程，给定任意一种来流扰动形式，求解得到任意来流扰动形式下广义坐标的时域变化过程，结合尾涡强度与广义坐标的降阶关系，恢复出任意来流扰动形式下尾涡强度的时域变化过程；
S6：将所述广义坐标的时域变化过程代入所述涡格法状态空间方程的输出方程中，得到时域气动力。
2.如权利要求1所述的基于状态空间形式涡格法的气动力降阶方法，其特征在于，步骤S2，具体包括：
在每个气动网格的四条边上布置涡线段，每个气动网格上布置的四条涡线段的强度相等且四条涡线段首尾相连构成一个涡环；选择每个气动网格的1/4弦线的中点为气动力的作用点，选择每个气动网格的3/4弦线的中点为气动力的控制点；
将气动网格中的涡分为机翼表面附着涡、机翼后缘第一排尾涡以及其他尾涡三部分，设机翼表面附着涡的强度列向量为Γb，机翼后缘第一排尾涡的强度列向量为Γw0，其他尾涡的强度列向量为Γwl，则气动控制方程为：
KbΓb+Kw0Γw0+KwlΓwl＝-w (1)
其中，w＝(V∞+Vg)·n，表示机翼表面的诱导速度；V∞表示来流的速度矢量，Vg表示来流扰动的速度矢量，n表示控制点处的法向量列阵；Kb表示机翼表面附着涡的诱导系数矩阵，Kw0表示机翼后缘第一排尾涡的诱导系数矩阵，Kwl表示其他尾涡的诱导系数矩阵；
机翼后缘第一排尾涡在脱出的过程中保持强度守恒，表达式为：
Figure FDA0003144536100000021
其中，
Figure FDA0003144536100000022
表示机翼后缘第一排尾涡的强度关于时间的导数；Δt表示时间步长；C1表示保证机翼后缘第一排尾涡与机翼表面附着涡对应关系正确的系数矩阵，包含0和1两种元素；
规定其他尾涡在脱出后保持强度不变，表达式为：
Figure FDA0003144536100000023
其中，
Figure FDA0003144536100000024
表示其他尾涡的强度关于时间的导数；C2和C3为表征其他尾涡与机翼后缘第一排尾涡位置对应关系的常数提取矩阵，包含0和1两种元素；I表示单位矩阵；
综合方程(1)、(2)和(3)得到涡格法状态空间方程中的状态方程为：
Figure FDA0003144536100000025
其中，Γw＝[Γw0 Γwl]T，表示尾涡强度；Aa和Ba表示状态空间系数矩阵，仅与气动面几何形状以及机翼表面附着涡、机翼后缘第一排尾涡和其他尾涡的划分有关，表达式为：
Figure FDA0003144536100000031
Figure FDA0003144536100000032
其中，O表示只包含0元素的零矩阵；
气动网格上的压强差表示为：
Figure FDA0003144536100000033
其中，下标ij表示沿展向第i个、沿弦向第j个气动网格中对应的物理量，Δpij表示沿展向第i个、沿弦向第j个气动网格中的压强差，Vl,ij表示沿展向第i个、沿弦向第j个气动网格中的当地来流速度，Γb,ij表示沿展向第i个、沿弦向第j个气动网格中机翼表面附着涡的涡强，ρ∞表示来流大气密度，τ1表示气动力作用点沿当地速度方向的切向量，τ2表示气动力作用点沿弦向的切向量，t表示仿真时间；假设Vl＝V∞，Vl表示当地速度方向矢量，V∞表示来流速度方向矢量；
方程(7)中沿展向第i个、沿弦向第j个气动网格中机翼表面附着涡的涡强的变化量由下述差分方程求解得到：
Figure FDA0003144536100000034
Figure FDA0003144536100000035
其中，Γb,i-1,j表示沿展向第i-1个、沿弦向第j个气动网格中机翼表面附着涡的涡强，Γb,i,j-1表示沿展向第i个、沿弦向第j-1个气动网格中机翼表面附着涡的涡强，Δcij表示气动网格沿弦向的几何长度，Δbij表示气动网格沿展向的几何长度；
作用在气动网格上的气动力表示为：
Fij＝ΔpijSijn (10)
其中，Sij表示气动网格的面积；
综合方程(1)、(7)、(8)和(9)得到涡格法状态空间方程中的输出方程为
Figure FDA0003144536100000041
其中，F表示气动力矢量，B1表示机翼表面诱导速度w与气动力矢量F之间关系的系数矩阵，B2表示机翼表面诱导速度w关于时间的导数
Figure FDA0003144536100000042
与气动力矢量F之间关系的系数矩阵，B3为表示尾涡强度Γw与气动力矢量F之间关系的系数矩阵；方程(4)和方程(11)构成涡格法状态空间方程。
3.如权利要求2所述的基于状态空间形式涡格法的气动力降阶方法，其特征在于，步骤S3中，给定离散阵风作为来流扰动形式：
wt＝(V∞+Vgt)·n (12)
其中，
Figure FDA0003144536100000043
表示离散阵风的扰动速度；wgm表示离散阵风的幅值，ω表示离散阵风的频率；
当给定离散阵风作为扰动形式时，w＝wt，尾涡强度Γw＝Γwt。
4.如权利要求3所述的基于状态空间形式涡格法的气动力降阶方法，其特征在于，步骤S4，具体包括：
将尾涡强度在预设时间内的时域变化过程数据Γwt在时间域内的离散点[t1 t2…ts…tm]下的数据作为快照数据，整理成如下快照矩阵：
Q＝[Γwt(t1) Γwt(t2)…Γwt(ts)…Γwt(tm)] (13)
其中，Γwt(ts)表示尾涡强度在ts时刻下的数值向量，s＝1,2,…,m，m表示时间离散点总数；
将快照矩阵中的数据列向量取均值：
Figure FDA0003144536100000044
将快照矩阵中的每一列元素减去均值
Figure FDA0003144536100000045
得到新的快照矩阵：
Figure FDA0003144536100000046
定义相关矩阵：
Figure FDA0003144536100000051
计算相关矩阵的非零特征值及特征向量：
Figure FDA0003144536100000052
其中，g＝1,2,…,l，l表示保留的特征向量总数，l≤m，λg表示第g阶特征值，λ1≥λ2≥…λl＞0，
Figure FDA0003144536100000053
表示第g阶特征向量；
构造最优POD基底φg如下：
Figure FDA0003144536100000054
其中，φg是一组标准正交基；
选择前k个POD基底构成k维POD基向量矩阵：
Φ＝[φ1 φ2…φk] (19)
其中，k<l；基于涡格法状态空间方程，得到尾涡强度与广义坐标的降阶关系：
Γw＝Φp (20)
其中，p表示广义坐标列向量；将方程(20)代入涡格法状态空间方程中的状态方程(4)后左乘ΦT，得到POD气动力降阶方程如下：
Figure FDA0003144536100000055
其中，A＝ΦTAaΦ，B＝ΦTBa，A和B为POD气动力降阶方程的状态空间系数矩阵，Aa和Ba为涡格法状态空间方程中的状态方程的系数矩阵。
5.如权利要求4所述的基于状态空间形式涡格法的气动力降阶方法，其特征在于，步骤S6中，时域气动力的表达式为：
Figure FDA0003144536100000056
CN202110746621.4A 2021-07-02 2021-07-02 一种基于状态空间形式涡格法的气动力降阶方法 Active CN113536457B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110746621.4A CN113536457B (zh) 2021-07-02 2021-07-02 一种基于状态空间形式涡格法的气动力降阶方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110746621.4A CN113536457B (zh) 2021-07-02 2021-07-02 一种基于状态空间形式涡格法的气动力降阶方法

Publications (2)

Publication Number Publication Date
CN113536457A true CN113536457A (zh) 2021-10-22
CN113536457B CN113536457B (zh) 2023-05-16

Family

ID=78097566

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110746621.4A Active CN113536457B (zh) 2021-07-02 2021-07-02 一种基于状态空间形式涡格法的气动力降阶方法

Country Status (1)

Country Link
CN (1) CN113536457B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114638133A (zh) * 2022-02-28 2022-06-17 桂林理工大学 一种瞬变电磁小波伽辽金2.5维正演算法
CN117236104A (zh) * 2023-08-22 2023-12-15 西南科技大学 一种基于环境模拟时序仿真的特种阀高阶模型降阶方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100250205A1 (en) * 2009-03-31 2010-09-30 Airbus Espana, S.L. Method and system for a quick calculation of aerodynamic forces on an aircraft in transonic conditions
ES2356788A1 (es) * 2007-12-18 2011-04-13 Airbus España S.L. Método y sistema para un cálculo rápido de las fuerzas aerodinámicas en una aeronave.
CN108052772A (zh) * 2017-12-30 2018-05-18 北京航空航天大学 一种基于结构降阶模型的几何非线性静气动弹性分析方法
CN110162822A (zh) * 2019-03-19 2019-08-23 北京机电工程研究所 耦合结构模态的时域快速非定常气动力计算方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
ES2356788A1 (es) * 2007-12-18 2011-04-13 Airbus España S.L. Método y sistema para un cálculo rápido de las fuerzas aerodinámicas en una aeronave.
US20100250205A1 (en) * 2009-03-31 2010-09-30 Airbus Espana, S.L. Method and system for a quick calculation of aerodynamic forces on an aircraft in transonic conditions
CN108052772A (zh) * 2017-12-30 2018-05-18 北京航空航天大学 一种基于结构降阶模型的几何非线性静气动弹性分析方法
CN110162822A (zh) * 2019-03-19 2019-08-23 北京机电工程研究所 耦合结构模态的时域快速非定常气动力计算方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114638133A (zh) * 2022-02-28 2022-06-17 桂林理工大学 一种瞬变电磁小波伽辽金2.5维正演算法
CN114638133B (zh) * 2022-02-28 2024-04-30 桂林理工大学 一种瞬变电磁小波伽辽金2.5维正演算法
CN117236104A (zh) * 2023-08-22 2023-12-15 西南科技大学 一种基于环境模拟时序仿真的特种阀高阶模型降阶方法
CN117236104B (zh) * 2023-08-22 2024-04-02 西南科技大学 一种基于环境模拟时序仿真的特种阀高阶模型降阶方法

Also Published As

Publication number Publication date
CN113536457B (zh) 2023-05-16

Similar Documents

Publication Publication Date Title
Reuther et al. Aerodynamic shape optimization of wing and wing-body configurations using control theory
Raveh CFD-based models of aerodynamic gust response
CN113536457B (zh) 一种基于状态空间形式涡格法的气动力降阶方法
Ritter et al. An enhanced modal approach for large deformation modeling of wing-like structures
Jacobson et al. Flutter Analysis with Stabilized Finite Elements based on the Linearized Frequency-domain Approach
Sadeghi et al. Application of three-dimensional interfaces for data transfer in aeroelastic computations
Gumbert et al. Effect of random geometric uncertainty on the computational design of a 3D flexible wing
Kim Parametric model reduction for aeroelastic systems: Invariant aeroelastic modes
Martins et al. Aero-structural wing design optimization using high-fidelity sensitivity analysis
CN113128085B (zh) 一种基于状态空间形式涡格法的阵风减缓分析方法
Jacobson et al. Multiscale mesh adaptation for transonic aeroelastic flutter problems
Bisson et al. Adjoint-based aerodynamic optimization framework
Bekemeyer et al. Reduced order transonic aeroelastic gust response simulation of large aircraft
Patel et al. Aeroelastic Wing Design Sensitivity Analysis with SU2-Nastran Coupling in OpenMDAO
Raul et al. Multifidelity modeling similarity conditions for airfoil dynamic stall prediction with manifold mapping
Rozov et al. Antisymmetric boundary condition for small disturbance CFD
Mishra et al. Time-dependent adjoint-based optimization for coupled aeroelastic problems
Da Ronch et al. Benchmarking CEASIOM software to predict flight control and flying qualities of the B-747
Mishra et al. Helicopter rotor design using adjoint-based optimization in a coupled cfd-csd framework
Yang et al. Improved Fluid-Structure Interface for Aeroelastic Computations with Non-Matching Outer Mold Lines
Cross et al. Continuum shape sensitivity with spatial gradient reconstruction of nonlinear aeroelastic gust response
Lipták et al. LPV model reduction methods for aeroelastic structures
Widhalm et al. Comparison of Reduced Order Models for Evaluating Stability Derivatives for the DLR-F22 ONERA model
Agostinelli et al. Flexible Wing Twist Optimisation using Rapid Computational Methods
Bekemeyer Rapid Computational Aerodynamics for Aircraft Gust Response Analysis

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