CN116305523B - 基于离散伴随的层流翼型优化方法 - Google Patents
基于离散伴随的层流翼型优化方法 Download PDFInfo
- Publication number
- CN116305523B CN116305523B CN202310029337.4A CN202310029337A CN116305523B CN 116305523 B CN116305523 B CN 116305523B CN 202310029337 A CN202310029337 A CN 202310029337A CN 116305523 B CN116305523 B CN 116305523B
- Authority
- CN
- China
- Prior art keywords
- transition
- rans
- equation
- optimization
- disturbance
- 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 88
- 238000005457 optimization Methods 0.000 title claims abstract description 84
- 230000007704 transition Effects 0.000 claims abstract description 131
- 238000013461 design Methods 0.000 claims abstract description 42
- 230000008878 coupling Effects 0.000 claims abstract description 25
- 238000010168 coupling process Methods 0.000 claims abstract description 25
- 238000005859 coupling reaction Methods 0.000 claims abstract description 25
- 230000002829 reductive effect Effects 0.000 claims abstract description 24
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 10
- 230000008859 change Effects 0.000 claims abstract description 9
- 238000013519 translation Methods 0.000 claims abstract description 4
- 238000003199 nucleic acid amplification method Methods 0.000 claims description 15
- 230000036961 partial effect Effects 0.000 claims description 15
- 230000003321 amplification Effects 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000005516 engineering process Methods 0.000 claims description 9
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000013507 mapping Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 230000009467 reduction Effects 0.000 description 9
- 230000008569 process Effects 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 230000004069 differentiation Effects 0.000 description 4
- 239000000446 fuel Substances 0.000 description 4
- 230000035939 shock Effects 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 230000003111 delayed effect Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000002035 prolonged effect Effects 0.000 description 2
- 230000002441 reversible effect Effects 0.000 description 2
- 230000009194 climbing Effects 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000008014 freezing Effects 0.000 description 1
- 238000007710 freezing Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000010348 incorporation Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
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
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- 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
Abstract
本发明公开了一种基于离散伴随的层流翼型优化方法，方法中，通过自由变形参数化方法实现层流翼的几何参数化,基于逆距离权重的动网格通过给定的表面网格和对应的变形几何，计算每个网格变化前后的法向扭转角和对应的平移距离，从而获得变形后的空间网格；耦合基于简化eN方法的转捩预测方法和基于Spalart‑Allmaras一方程湍流模型的RANS方程求解流场状态变量Q和转捩位置状态变量Tr；通过基于RANS‑BLCode‑AFM转捩预测方法的耦合伴随方程求解各状态变量相对于设计变量的梯度，通过序列二次规划算法，根据所述梯度判断优化是否收敛，若未收敛，根据所述梯度确定下一步优化的设计变量值，然后重复直到收敛。
Description
技术领域
本发明属于航天航空技术领域，特别是一种基于离散伴随的层流翼型优化方法。
背景技术
燃油成本占当今商用飞机运营成本的很大一部分，减少燃油消耗能极大地提高民机的经济性，因此提高飞机燃油效率一直是一个重大挑战。研究表明对于常规布局，层流机翼技术能大约减阻10％，对应减少5％的燃油消耗，巨大的减阻潜力使其成为航空领域的研究热点。相比传统的机翼气动外形优化，层流翼气动外形优化需要同时在层流区域和激波控制之间权衡，因此更加复杂。为了能预测出准确的转捩位置，需要采用更高精度的CFD求解器以及合适的湍流模型。因此，高可信度转捩预测方法和高效、鲁棒、精确的考虑转捩的耦合离散伴随方法，是发展层流翼梯度优化设计方法的关键之一。
通常在中小后掠角(<20°)/中低雷诺数(≤20×106)条件下，层流翼的转捩主要是由TS波扰动和分离流失稳主导的，这与二维翼型的层流转捩物理现象相似。因此，发展一种基于离散伴随的层流翼型优化设计方法具有重要意义。
受制于层流到湍流转捩以及考虑转捩的耦合伴随方程构造与求解的复杂性，目前基于离散伴随的气动优化方法的应用主要集中在全湍流问题，针对层流翼型的优化设计，基于离散伴随的理论研究较少，发展还远不完善。目前考虑转捩的耦合伴随方程主要存在两大问题：冻结转捩模块，忽略转捩模块对梯度的影响；组装耦合伴随方程的偏导数矩阵时，采用了依赖于设计问题维度的有限差分或者复变量微分，降低了梯度求解效率和精度。
考虑到开发难度、可移植性以及拓展性要求，国际上针对基于离散伴随的梯度求解，开发了不同的自动微分工具进行偏导数求解。但是，自动微分工具有使用限制，如Tapenade在涉及MPI并行/串行交互计算等问题时失效，需要结合人工手动推导微分程序。同时，自动微分程序直接对应的微分代码，往往效率较低，需结合人工模块化处理，来提高偏导数的计算效率。
在背景技术部分中公开的上述信息仅仅用于增强对本发明背景的理解，因此可能包含不构成本领域普通技术人员公知的现有技术的信息。
发明内容
本发明的目的是通过以下技术方案予以实现，一种基于离散伴随的层流翼型优化方法，其特征在于，其包括以下步骤，
步骤S1：通过自由变形参数化方法实现层流翼的几何参数化,包括：建立自由变形参数化方法的控制框和物面的映射关系；针对层流翼的第一几何外形，利用自由变形控制点改变所述第一几何外形以获得变形后的层流翼的第二几何外形，其中，所述控制点在某个维度上对应的变化量为设计变量；
步骤S2：基于逆距离权重的动网格技术，分别对所述第一几何外形和第二几何外形进行处理，以获得对应的第一表面网格和变形后的第二表面网格，并计算第一表面网格因变形而导致的法向扭转角和对应的平移距离，从而获得变形后的空间网格；
步骤S3：基于简化eN方法的转捩预测方法和基于Spalart-Allmaras一方程湍流模型的RANS方程，以变形后的空间网格作为几何输入量，并基于给定的状态输入量，求解流场状态变量Q和转捩位置状态变量Tr；
步骤S4：通过基于RANS-BLCode-AFM转捩预测方法的耦合伴随方程,求解流场状态变量Q和转捩位置状态变量Tr相对于所述设计变量的梯度；
步骤S5：根据所述梯度和序列二次规划算法，判断优化是否收敛，若未收敛，则根据所述梯度和序列二次规划算法确定新的设计变量，然后以新的设计变量重复进行步骤S1-S4，直到收敛。
优选的，步骤3中，先利用RANS方程进行固定转捩计算流场以得到输出A、第一次的流场残差……第N次的流场残差，当第N次流场残差相比第一次流场残差降低10-8量级后，判定为内部收敛无需继续计算，将此时第N次计算所得的流场解和输入的空间网格作为简化eN方法的转捩预测方法中的转捩模块的输入，进一步获得转捩位置(x)，并采用间歇因子函数修正，之后将新的转捩位置返回RANS方程，循环迭代，直至流场解收敛。
优选的，步骤3中，RANS求解器的求解方程表示为：
其中为RANS方程的残差，Q是RANS方程的状态变量，转捩预测问题描述为：
其中，为转捩模块的残差，Tr为转捩位置状态变量，Tp为转捩位置预测值，转捩问题描述为以下形式：
其中X为所述变形后的空间网格对应的几何输入量，R为转捩残差。
优选的，步骤4中，针对全湍流状态，耦合伴随方程表示为：
其中为RANS方程的残差，Q为RANS方程对应的状态变量，其包含密度ρ、速度v、压力p及湍流变量ν，ψ为RANS方程对应的伴随向量，I为目标函数，将残差向量，状态变量以及对应的伴随向量分成关于RANS方程和转捩预测方法的两部分，
其中，残差变量状态变量Y及对应的伴随向量Ψ，φ为转捩模块对应的伴随向量，
RANS-BLCode-AFM的耦合伴随方程为：
代表RANS方程残差关于其状态变量Q的雅可比矩阵，/>表示转捩状态变量Tr的变化对RANS方程残差的影响，当转捩位置发生变化时，对应的间歇因子改变，从而改变了涡粘性系数，最终也导致了RANS方程其它状态变量的变化；/>为转捩模块的残差关于转捩位置的雅可比矩阵；/>表示转捩模块残差关于RANS状态变量Q的偏导数；/>代表目标函数对RANS方程状态的偏导数；/>为目标函数关于转捩位置的偏导数，当目标函数为升力系数Cl，阻力系数Cd以及力矩系数Cm时，其直接由RANS的状态变量决定，该项偏导数为0，而当目标函数为转捩位置时，该项值不为0。
优选的，基于简化eN方法的转捩预测方法具体包括如下步骤：
基于Navier-Stokes方程，在时均流动上，引入非定常的正弦扰动，非定常守恒变量表示为：
其中q代表三个方向的速度u，v，w，压力p或者密度ρ，x,y,z代表笛卡尔坐标系下的三个正交方向，x为流向站位，y在物面法向，z为展向，t为时间，公式右端第一项代表时均流，第二项表示扰动量，扰动量为
扰动相位取决于法向的流动，α和β分别代表x和z方向的扰动波数，ω为扰动频率，
流向扰动增长率表达为
对该公式两边分别取对数，并对其进行积分，得到
式中：A为波幅；为初始波幅；x0为在给定的扰动频率下，如果扰动放大率αi小于0，开始出现扰动的位置，N＝ln(A/A0)定义为扰动放大因子，由扰动放大率αi沿流向积分获得，eN方法通过对扰动放大率αi沿流向积分，得到不同流向位置的扰动放大因子，当某一位置的N值达到转捩阈值时，认为发生转捩。
优选的，所述扰动放大因子N与流向站位x的关系表示为：
其中，
式中：为动量厚度雷诺数；Hk为可压缩形状因子；Ue、μe分别为边界层边界处的合速度和运动学粘性系数；δ2为边界层动量厚度；ρ为密度；m(Hk)、l(Hk)均指与Hk相关的经验公式，对dNdx沿流向积分，得到不同流向站位的放大因子N值，对比不同流向站位处的N值与转捩阈值，得到转捩位置。
所述的基于离散伴随的层流翼型优化方法中，步骤3中，RANS方程中，先进行固定转捩计算，流场残差降低10-8量级后，将流场解和输入的空间网格作为简化eN方法的转捩预测方法中的转捩模块的输入，进一步获得转捩位置，并采用间歇因子函数进行修正，之后将新的转捩位置返回RANS方程，循环迭代，直至流场解收敛。
所述的基于离散伴随的层流翼型优化方法中，步骤3中，RANS求解器的求解方程表示为：
其中为RANS方程的残差，Q是RANS方程的状态变量，转捩预测问题描述为：
其中，为转捩模块的残差，Tr为转捩位置状态变量，Tp为转捩位置预测值，转捩问题描述为以下形式：
其中X是几何设计变量，R为转捩残差。
所述的基于离散伴随的层流翼型优化方法中，步骤4中，针对全湍流状态，耦合伴随方程表示为：
其中为RANS方程的残差，Q为RANS方程对应的状态变量，其包含密度ρ、速度v、压力p及湍流变量ν，ψ为RANS方程对应的伴随向量，I为目标函数，将残差向量，状态变量以及对应的伴随向量分成关于RANS方程和转捩预测方法的两部分，
其中，残差变量状态变量Y及对应的伴随向量Ψ，φ为转捩模块对应的伴随向量。
RANS-BLCode-AFM的耦合伴随方程为：
代表RANS方程残差关于其状态变量Q的雅可比矩阵，/>表示转捩状态变量Tr的变化对RANS方程残差的影响，当转捩位置发生变化时，对应的间歇因子改变，从而改变了涡粘性系数，最终也导致了RANS方程其它状态变量的变化；/>为转捩模块的残差关于转捩位置的雅可比矩阵；/>表示转捩模块残差关于RANS状态变量Q的偏导数；
代表目标函数对RANS方程状态的偏导数；/>为目标函数关于转捩位置的偏导数，当目标函数为Cl，Cd以及Cm时，Cl:升力系数，Cd：阻力系数，Cm：力矩系数，其直接由RANS的状态变量决定，该项偏导数为0，而当目标函数为转捩位置时，该项值不为0。
所述的基于离散伴随的层流翼型优化方法中，基于Navier-Stokes方程，在时均流动上，引入非定常的正弦扰动，非定常守恒变量表示为：
其中q代表三个方向的速度u，v，w，压力p或者密度ρ，x,y,z代表笛卡尔坐标系下的三个正交方向，x为流向，y在物面法向，z为展向，t为时间，公式右端第一项代表时均流，第二项表示扰动量，扰动量为
扰动相位取决于法向的流动，α和β分别代表x和z方向的扰动波数，ω为扰动频率，
流向扰动增长率表达为
对该公式两边分别取对数，并对其进行积分，得到
式中：A为波幅；为初始波幅；x0为在给定的扰动频率下，如果扰动放大率αi小于0，开始出现扰动的位置，N＝ln(A/A0)定义为扰动放大因子，由扰动放大率αi沿流向积分获得，eN方法通过对扰动放大率αi沿流向积分，得到不同流向位置的扰动放大因子，当某一位置的N值达到转捩阈值时，认为发生转捩。
所述的基于离散伴随的层流翼型优化方法中，基于AFM模型，对于特征解相似的边界层流动，放大因子N与流向站位x的关系表示为：
其中，
式中：为动量厚度雷诺数；Hk为可压缩形状因子；Ue、μe分别为边界层边界处的合速度和运动学粘性系数；δ2为边界层动量厚度；ρ为密度；m(Hk)、l(Hk)均指与Hk相关的经验公式。对dNdx沿流向积分，得到不同流向站位的放大因子N值，对比不同流向站位处的N值与转捩阈值，得到转捩位置。
阈值为：
其中，H12为形状因子，根据C1判据准则，当ReCF达到时，定义该点为转捩位置；
和现有技术相比，本发明具有以下优点：本发明构建了高效、可靠的层流翼型梯度优化设计方法。基于无矩阵存储技术，采用混合反向自动微分，构造了解析形式的、考虑转捩的耦合伴随方程，实现了梯度信息高效、精确求解，显著降低了存储要求，并通过CK数值算法，极大地提高了大规模耦合伴随线性方程组的求解效率。
附图说明
通过阅读下文优选的具体实施方式中的详细描述，本发明各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。说明书附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。显而易见地，下面描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。而且在整个附图中，用相同的附图标记表示相同的部件。
在附图中：
图1是本发明一个实施例提供的基于离散伴随的层流翼型优化方法的流程示意图；
图2是本发明一个实施例提供的基于离散伴随的层流翼型优化方法的基于RANS方程求解的转捩预测系统结构示意图；
图3是本发明实例一提供的基于离散伴随的层流翼型优化方法的RAE2822翼型初始网格示意图；
图4是本发明实例一提供的基于离散伴随的层流翼型优化方法的RAE2822翼型FFD控制点示意图；
图5是本发明实例一提供的基于离散伴随的层流翼型优化方法的单点优化翼型结果对比示意图；
图6是本发明实例一提供的基于离散伴随的层流翼型优化方法的单点优化压力分布结果对比图；
图7是本发明实例一提供的基于离散伴随的层流翼型优化方法的单点层流翼型优化的收敛历程示意图；
图8是本发明实例二提供的基于离散伴随的层流翼型优化方法的多点优化初始网格示意图；
图9是本发明实例二提供的基于离散伴随的层流翼型优化方法的多点优化FFD控制点示意图；
图10是本发明实例二提供的基于离散伴随的层流翼型优化方法的多点优化翼型结果对比示意图；
图11是本发明实例二提供的基于离散伴随的层流翼型优化方法的多点优化压力分布结果对比示意图；
图12是本发明实例二提供的基于离散伴随的层流翼型优化方法的多点层流翼型优化的收敛历程示意图；
以下结合附图和实施例对本发明作进一步的解释。
具体实施方式
下面将参照附图更详细地描述本发明的具体实施例。虽然附图中显示了本发明的具体实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。
需要说明的是，在说明书及权利要求当中使用了某些词汇来指称特定组件。本领域技术人员应可以理解，技术人员可能会用不同名词来称呼同一个组件。本说明书及权利要求并不以名词的差异来作为区分组件的方式，而是以组件在功能上的差异来作为区分的准则。如在通篇说明书及权利要求当中所提及的“包含”或“包括”为一开放式用语，故应解释成“包含但不限定于”。说明书后续描述为实施本发明的较佳实施方式，然所述描述乃以说明书的一般原则为目的，并非用以限定本发明的范围。本发明的保护范围当视所附权利要求所界定者为准。
为便于对本发明实施例的理解，下面将结合附图以具体实施例为例做进一步的解释说明，且各个附图并不构成对本发明实施例的限定。
如附图1所示的一种基于流场预测的气动外形优化方法，针对层流翼型的气动外形优化需要考虑层流区域和激波控制之间权衡的问题，将高可信度转捩预测方法和高效、鲁棒、精确的考虑转捩的耦合离散伴随方法相结合，建立了基于耦合离散伴随方程的层流翼优化设计框架，具体包括以下步骤：
S1：通过FFD(Free Form Deform)参数化方法实现几何参数化。通过建立FFD控制框和物面的映射关系，利用FFD控制点改变机翼几何外形
S2：利用基于IDW(Inverse Distance Weighting)的动网格技术，通过给定的表面网格和对应的变形几何，计算每个网格单元变化前后的法向扭转角和对应的平移距离，从而获得变形后的空间网格。
S3：利用变形后的空间网格及输入变量(如计算攻角、马赫数等)，耦合基于简化eN方法的转捩预测方法和基于Spalart-Allmaras一方程湍流模型的RANS方程，求解流场状态变量Q和Tr。
S4：通过基于RANS-BLCode-AFM转捩预测方法的耦合伴随方程求解各状态变量相对于设计变量的梯度。
S5：通过序列二次规划(Sequential Quadratic Programming：SQP)算法，根据S4提供的梯度判断优化是否收敛。若未收敛，根据梯度信息确定下一步优化的设计变量值，然后重复进行S1-S5，从而实现多模块的自动优化设计。
如附图2所示的基于RANS方程求解的转捩预测系统，将变形后的空间网格(即图2中的计算网格)及输入变量(如计算攻角、马赫数等)作为RANS求解器的输入，先进行固定转捩计算，流场残差降低10-8量级后，将流场解和输入的空间网格作为转捩模块的输入，进一步获得转捩位置，并采用间歇因子函数进行修正，之后将新的转捩位置返回RANS方程，循环迭代，直至流场解收敛。
基于建立的RANS-BLCode-AFM转捩预测方法，构造了考虑转捩的耦合伴随方程。并耦合FFD参数化方法、动网格技术以及梯度优化算法构建了层流翼型梯度优化设计框架。
实施例一
参考Cessna 172SP Skyhawk的设计状态，采用建立的层流翼/全湍流优化设计系统(附图1)进行优化设计研究。Cessna 172SP Skyhawk的设计状态为：马赫数M＝0.19，雷诺数Re＝5×106，升力系数Cl＝0.3。初始翼型为RAE2822，初始网格如附图3，流向一周为281个网格点，法向121个网格点。FFD控制点见附图4由于转捩对翼型头部比较敏感，所以FFD控制点在头部较密。针对该优化问题，湍流度约为0.07％，因而转捩阈值NLTS＝9.0。优化问题以及相关优化参数见表1，其中，Cd为阻力系数，AoA为攻角，优化过程中保证翼型面积S满足约束S≥Sinitial，厚度ty满足约束ty≥0.3tyinitial。
表1单点优化问题
优化结果见表2，其中第二行是模拟类型，即针对初始/优化构型，进行层流–湍流转捩数值模拟(“LT”)。例如“Optimized-LT”代表考虑转捩预测的优化构型的转捩预测结果。翼型的优化前后的对比见附图5，Initial为初始翼型，Optimized为层流翼型优化后的构型。附图6给出了Initial及Optimized的压力分布对比。
表2优化结果对比
表2给出了初始构型和优化构型的总阻力系数、压差阻力系数、粘性阻力系数、升力系数及对应该升力系数的攻角AoA。对于层流翼型优化，从表2中可以看到，相比初始翼型(Initial)，优化后的层流翼型(Optimized)总阻力系数减小32.468counts，即降低52.27％。其中压差阻力系数降低了13.385counts，相比初始构型降低63.87％，摩擦阻力系数降低了19.084counts，相比初始构型减小46.37％。摩擦阻力系数的减小量占了总阻力系数减小量的58.78％。在优化过程中，层流翼型的攻角从0.4499°变为0.6361°，变化量较小。附图5中，优化得到的层流翼型相对初始翼型头部半径减小，最大厚度显著后移。因而优化后的压力分布(见附图6)，有更长的顺压梯度区(顺压梯度有利于抑制TS波的扰动)，从而使得上翼面转捩位置由45.3％弦长推迟到75.7％弦长，下翼面转捩位置由47.9％弦长推迟到73.4％弦长。表明层流翼型优化系统可以有效的推迟转捩，降低摩擦阻力系数。
附图7给出了层流翼型优化的收敛过程，经过35次主迭代，满足了约束以及目标函数收敛条件。从主迭代的收敛过程可以看到，梯度优化在前面10步以后可以快速的降低阻力，在最后阶段，阻力降低较小，主要是压力分布形态发生变化，最终给出较光滑的压力分布。
实施例二
参考Honda Jet的设计状态，采用建立的层流翼/全湍流优化设计系统(附图1)进行优化设计研究。Honda Jet飞机的起飞重量为4082kg。其巡航状态(第一个设计状态)的升力系数为Cl＝0.26，雷诺数Re＝11.7×106，马赫数为M＝0.69。该飞机的爬升状态(第二个设计状态)的升力系数为Cl＝0.35，对应的Re＝13.6×106，M＝0.31。同时，考虑了另外了几个状态来降低优化结果对飞行条件变化的敏感性。第三个选择的巡航状态为Cl＝0.18，Re＝11.7×106，M＝0.69，目的是降低型阻。第四个选择的状态为：Cl＝0.38，Re＝7.93×106，M＝0.7，约束为俯仰力矩不小于-0.04，目的是在低雷诺数、高马赫数下最小化配平阻力。优化设计状态以及优化目标的权重见表3。优化过程对翼型面积的约束为S≥Sinit，对翼型厚度的约束为ty≥0.3tyinit。多点优化问题描述见表4。与单点的翼型优化相同，多点优化的TS波放大因子阈值定义为NLTS＝9.0。选择一个典型的层流翼型作为初始翼型进行优化，初始翼型具备典型的顺压梯度压力分布。翼型的空间网格见附图8。翼型法向网格点个数为121，周向点个数为281。初始翼型以及对应的FFD控制点见附图9。对多点优化，同样分别进行了层流翼型和全湍流翼型优化，并进行对比分析。
表3多点优化设计状态以及优化目标的权重
表4多点优化问题
初始翼型以及优化后翼型的气动力系数见表5，xtr/c(U)及xtr/c(L)分别代表翼型上表面和下表面的相对转捩位置，即转捩位置与翼型弦长比值。其中Point1、Point2、Point3以及Point4分别代表不同的设计状态。类似的，表中第3行和第14行为数值模拟类型，即针对初始/优化构型，进行转捩数值模拟(“LT”)。附图10给出了初始翼型(Initial)以及优化后的层流翼型(Optimzed)。附图11给出了初始翼型(Initial)以及优化翼型(Optimized)的压力分布对比。对于层流翼型优化，对比初始翼型(Initial)以及优化翼型(Optimized)的压力分布(附图11，黑色线和红色线)，优化翼型(Optimized)上翼面的吸力峰值在所有设计状态都提高，从而减小了顺压梯度。
表5多点优化结果对比
优化后，升力系数以及力矩系数都满足约束，同时，每一个设计状态下的阻力都减小(表5)。对于主设计点Point1，优化后的阻力(Optimized)相比初始构型(Initial)阻力降低了16.334counts，即33.2％。其中摩擦阻力占了10.18counts，压差阻力占了6.255counts。优化翼型(Optimized)上翼面的转捩位置由41.3％弦长延长到56.1％，而下翼面的转捩位置由58.3％弦长延长到70.9％。对比Point1，针对初始构型，可以看到对于初始构型(Initial)下翼面的转捩位置较对应的Point2靠下游。主要是Point2下翼面前缘存在逆压梯度，有利于TS波扰动的发展，从而提前了转捩位置。对于Point2，优化翼型(Optimized)相比初始构型(Initial)总阻力减小了28.64counts，即47.0％。其中摩擦阻力减小了19.115counts，压差阻力减小了10.637counts。上翼面的转捩位置由43.0％弦长推迟到56.3％弦长，下翼面的转捩位置由34.4％弦长推迟到70.6％弦长。对于爬升点的设计状态，即Point3，优化翼型(Optimized)相比初始构型(Initial)总阻力降低了4.525counts，即9.5％。摩擦阻力降低了2.845counts，压差阻力降低了1.681counts。在上翼面，转捩位置由初始的40.4％弦长提前到37.5％弦长。而下翼面的转捩位置由初始的63.8％推迟到73％。对于Point4，优化后的构型相比初始构型，总阻力降低了20.393counts，即33.85％。其中压差阻力降低了17.705counts，摩擦阻力降低了2.687counts。优化的层流翼型上翼面转捩位置由47.3％推迟到50.3％，而下翼面的转捩位置由61.8％推迟到74.8％。对于压力分布，优化之后消除了激波，因而激波阻力(压差阻力)大大降低。所以，Point4阻力的降低主要来源于压差阻力的降低。
附图12给出了优化的目标函数收敛图，经过30次主迭代之后，目标函数达到收敛标准，优化停止。目标函数由初始52.320counts的降低到35.629counts。
与现有技术相比，本发明基于RANS-BLCode-AFM转捩预测方法的耦合伴随方程，并耦合FFD参数化方法、动网格技术以及梯度优化算法构建了高效、可靠的层流翼梯度优化设计方法。基于无矩阵存储技术，采用混合反向自动微分，构造解析形式的、考虑转捩的耦合伴随方程，实现了梯度信息高效、精确求解，显著降低了存储要求，并通过CK数值算法，极大地提高了大规模耦合伴随线性方程组的求解效率。本发明提供的优化设计方法，可以针对低速、跨声速二维层流翼型有效地开展单点以及多点优化设计研究，对飞行器气动外形设计具有重要应用价值。
尽管以上结合附图对本发明的实施方案进行了描述，但本发明并不局限于上述的具体实施方案和应用领域，上述的具体实施方案仅仅是示意性的、指导性的，而不是限制性的。本领域的普通技术人员在本说明书的启示下和在不脱离本发明权利要求所保护的范围的情况下，还可以做出很多种的形式，这些均属于本发明保护之列。
Claims (3)
1.一种基于离散伴随的层流翼型优化方法，其特征在于，其包括以下步骤，
步骤S1：通过自由变形参数化方法实现层流翼的几何参数化,包括：建立自由变形参数化方法的控制框和物面的映射关系；针对层流翼的第一几何外形，利用自由变形控制点改变所述第一几何外形以获得变形后的层流翼的第二几何外形，其中，所述控制点在某个维度上对应的变化量为设计变量；
步骤S2：基于逆距离权重的动网格技术，分别对所述第一几何外形和第二几何外形进行处理，以获得对应的第一表面网格和变形后的第二表面网格，并计算第一表面网格因变形而导致的法向扭转角和对应的平移距离，从而获得变形后的空间网格；
步骤S3：基于简化方法的转捩预测方法和基于Spalart-Allmaras 一方程湍流模型的RANS 方程，以变形后的空间网格作为几何输入量，并基于给定的状态输入量，求解流场状态变量/>和转捩位置状态变量/>，RANS求解器的求解方程表示为：
，
其中为RANS方程的残差，/>是RANS方程的状态变量，转捩预测问题描述为：
，
其中，为转捩模块的残差，/>为转捩位置状态变量，/>为转捩位置预测值，转捩问题描述为以下形式：
，
其中为所述变形后的空间网格对应的几何输入量，/>为转捩残差；
步骤S4：通过基于RANS-BLCode-AFM转捩预测方法的耦合伴随方程,求解流场状态变量和转捩位置状态变量/>相对于所述设计变量的梯度，针对全湍流状态，耦合伴随方程表示为：
，
其中为RANS方程的残差，/>为RANS方程对应的状态变量，其包含密度/>、速度/>、压力/>及湍流变量/>，/>为RANS方程对应的伴随向量，/>为目标函数，将残差向量，状态变量以及对应的伴随向量分成关于RANS方程和转捩预测方法的两部分，
，
其中，残差变量、状态变量/>及对应的伴随向量/>，/>为转捩模块对应的伴随向量，
RANS-BLCode-AFM的耦合伴随方程为：
，/>代表 RANS 方程残差关于其状态变量/>的雅可比矩阵，/>表示转捩状态变量/>的变化对 RANS 方程残差的影响，当转捩位置发生变化时，对应的间歇因子改变，从而改变了涡粘性系数，最终也导致了 RANS 方程其它状态变量的变化；/>为转捩模块的残差关于转捩位置的雅可比矩阵；表示转捩模块残差关于 RANS 状态变量/>的偏导数；/>代表目标函数对RANS 方程状态的偏导数；/>为目标函数关于转捩位置的偏导数，当目标函数为升力系数/>，阻力系数/>以及力矩系数/>时，其直接由 RANS 的状态变量决定，该项偏导数为0，而当目标函数为转捩位置时，该项值不为 0；
步骤S5：根据所述梯度和序列二次规划算法，判断优化是否收敛，若未收敛，则根据所述梯度和序列二次规划算法确定新的设计变量，然后以新的设计变量重复进行步骤S1-S4，直到收敛，经过35次主迭代，满足了约束以及目标函数收敛条件，梯度优化在前面 10 步以后降低阻力，最终给出光滑的压力分布；
其中，步骤S3中，先利用RANS方程进行固定转捩计算流场以得到输出A、第一次的流场残差……第N次的流场残差，当第N次流场残差相比第一次流场残差降低量级后，判定为内部收敛无需继续计算，将此时第N次计算所得的流场解和输入的空间网格作为简化/>方法的转捩预测方法中的转捩模块的输入，进一步获得转捩位置(x)，并采用间歇因子函数修正，之后将新的转捩位置返回RANS方程，循环迭代，直至流场解收敛。
2.根据权利要求1所述的基于离散伴随的层流翼型优化方法，其中，基于简化方法的转捩预测方法具体包括如下步骤：
基于Navier-Stokes方程，在时均流动上，引入非定常的正弦扰动，非定常守恒变量表示为：
，
其中代表三个方向的速度/>，/>，/>，压力/>或者密度/>，/>代表笛卡尔坐标系下的三个正交方向，/>为流向站位，/>在物面法向，/>为展向，t为时间，公式右端第一项代表时均流，第二项表示扰动量，扰动量为
，
扰动相位取决于法向的流动，/>和/>分别代表/>和/>方向的扰动波数，/>为扰动频率，
流向扰动增长率表达为
，
对该公式两边分别取对数，并对其进行积分，得到
，
式中：为波幅；/>为初始波幅；/>为在给定的扰动频率下，如果扰动放大率/>小于0，开始出现扰动的位置，/>定义为扰动放大因子，由扰动放大率沿流向积分获得，/>方法通过对扰动放大率/>沿流向积分，得到不同流向位置的扰动放大因子，当某一位置的/>值达到转捩阈值时，认为发生转捩。
3.根据权利要求2所述的基于离散伴随的层流翼型优化方法，其中，所述扰动放大因子与流向站位/>的关系表示为：
，
其中，
，
，
式中：为动量厚度雷诺数；/>为可压缩形状因子; />、/>分别为边界层边界处的合速度和运动学粘性系数；/>为边界层动量厚度; />为密度；/>、/>均指与相关的经验公式，对/>沿流向积分，得到不同流向站位的放大因子/>值，对比不同流向站位处的/>值与转捩阈值，得到转捩位置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310029337.4A CN116305523B (zh) | 2023-01-09 | 2023-01-09 | 基于离散伴随的层流翼型优化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310029337.4A CN116305523B (zh) | 2023-01-09 | 2023-01-09 | 基于离散伴随的层流翼型优化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116305523A CN116305523A (zh) | 2023-06-23 |
| CN116305523B true CN116305523B (zh) | 2023-11-14 |
Family
ID=86798600
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310029337.4A Active CN116305523B (zh) | 2023-01-09 | 2023-01-09 | 基于离散伴随的层流翼型优化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116305523B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116956781B (zh) * | 2023-09-18 | 2024-01-05 | 西北工业大学 | 基于rans考虑转捩的静气动弹性分析方法 |
| CN117077298B (zh) * | 2023-10-17 | 2023-12-29 | 中国科学院工程热物理研究所 | 一种基于梯度增强随机Co-Kriging模型的飞行器稳健优化设计方法 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109033525A (zh) * | 2018-06-27 | 2018-12-18 | 浙江大学 | 一种基于简化三方程转捩模型的高超声速转捩预测方法 |
-
2023
- 2023-01-09 CN CN202310029337.4A patent/CN116305523B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109033525A (zh) * | 2018-06-27 | 2018-12-18 | 浙江大学 | 一种基于简化三方程转捩模型的高超声速转捩预测方法 |
Non-Patent Citations (4)
| Title |
|---|
| Natural Laminar-Flow Airfoil Optimization Design Using a Discrete Adjoint Approach;Yayun Shi 等;《AIAA Journal》;第58卷(第11期);全文 * |
| Transition prediction and sensitivity analysis for a natural laminar flow wing glove flight experiment;Tihao YANG 等;《Chinese Journal of Aeronautics》;第34卷(第8期);第34–47页 * |
| 一种基于动模态分解的翼型流动转捩预测新方法;韩忠华;《航空学报》;第38卷(第1期);第1-17页 * |
| 基于离散伴随的层流翼优化设计方法;杨体浩 等;《航空学报》;第43卷(第12期);第1-20页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116305523A (zh) | 2023-06-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116305523B (zh) | 基于离散伴随的层流翼型优化方法 | |
| EP2466288B1 (en) | Method of designing natural laminar flow wing for reynolds numbers equivalent to actual supersonic aircraft | |
| Chauhan et al. | Rans-based aerodynamic shape optimization of a wing considering propeller–wing interaction | |
| Keye | Fluid-structure coupled analysis of a transport aircraft and flight-test validation | |
| CN116227023B (zh) | 考虑横流的层流翼梯度优化方法 | |
| CN116502338B (zh) | 一种通用化的基于线性稳定性理论的工程转捩预测方法 | |
| KR102616901B1 (ko) | 넓은 속도 범위 극초음속기의 공기역학적 구성 설계 방법 및 시스템 | |
| Chao et al. | Effect of flexibility on unsteady aerodynamics forces of a purely plunging airfoil | |
| Chyu et al. | Calculation of unsteady transonic flow over an airfoil | |
| Manshadi et al. | Optimizing a two-element wing model with morphing flap by means of the response surface method | |
| Kaul et al. | Drag characterization study of variable camber continuous trailing edge flap | |
| Tseng et al. | Calculation of aerodynamic characteristics of airplane configurations at high angles of attack | |
| Smith | An assessment of the state-of-the-art from the 2019 ARO dynamic stall workshop | |
| Lei et al. | Aerodynamic performance of distributed electric propulsion with wing interaction | |
| Dhileep et al. | Numerical Study of Camber Morphing in NACA0012 Airfoil | |
| Nikkhoo et al. | Effect of different aero-structural optimization in the commercial airplane | |
| Agostinelli et al. | Propeller-Wing Interaction using Rapid Computational Methods | |
| Tu | Numerical study of steady and unsteady canard-wing-body aerodynamics | |
| Kraljic et al. | Stall of airfoils with blunt noses at low to moderately high Reynolds numbers | |
| Weston et al. | Exercise of Knowledge-Based Aerodynamic Design Process on the AFRL SUGAR Model | |
| Morgan | High-lift flaps for natural laminar flow airfoils | |
| Duan et al. | Aerodynamic Simulation of Aircraft Crusing Characteristics Based on FLUENT | |
| Guo et al. | Numerical Simulation Research on Static Aeroelastic Effect of the Transonic Aileron of a High Aspect Ratio Aircraft. | |
| Kefalas et al. | Aerodynamic analysis of a light aircraft using computational fluid dynamics | |
| Lei et al. | Numerical Optimization of Leading-Edge Deflection Angles for an SST Configuration at Low Speed |
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 |