CN111695264A - 一种用于音爆传播计算的多波系同步推进波形参数方法 - Google Patents

一种用于音爆传播计算的多波系同步推进波形参数方法 Download PDF

Info

Publication number
CN111695264A
CN111695264A CN202010545227.XA CN202010545227A CN111695264A CN 111695264 A CN111695264 A CN 111695264A CN 202010545227 A CN202010545227 A CN 202010545227A CN 111695264 A CN111695264 A CN 111695264A
Authority
CN
China
Prior art keywords
waveform
parameter
propagation
new
length
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
CN202010545227.XA
Other languages
English (en)
Other versions
CN111695264B (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.)
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Original Assignee
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
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 Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center filed Critical Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Priority to CN202010545227.XA priority Critical patent/CN111695264B/zh
Publication of CN111695264A publication Critical patent/CN111695264A/zh
Application granted granted Critical
Publication of CN111695264B publication Critical patent/CN111695264B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M9/00Aerodynamic testing; Arrangements in or on wind tunnels
    • G01M9/02Wind tunnels
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M9/00Aerodynamic testing; Arrangements in or on wind tunnels
    • G01M9/08Aerodynamic models
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/10Geometric CAD
    • G06F30/15Vehicle, aircraft or watercraft design
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10Noise analysis or noise optimisation
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • Fluid Mechanics (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Automation & Control Theory (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种用于音爆传播计算的多波系同步推进波形参数方法。该方法基于音爆射线管追踪和Thomas的波形参数模型,采用内外双循环流程进行传播计算,具体步骤为:a.初始化;b.射线管单步跟踪(外循环起始);c.计算射线管单步推进中波形传播预期步长参数T(内循环起始);d.步长参数T符合性检验及其更新、各区间长度因子F1,i和F2,i的递进计算;e.多波系同步推进的波形参数传播计算;f.判断波形传播是否达到射线管末端(内循环终止判断);g.判断射线管是否达到终止位置(外循环终止判断)。该方法可以在一个时间步内完成多道激波生成或者激波融合,提高了计算效率和计算稳定性。

Description

一种用于音爆传播计算的多波系同步推进波形参数方法
技术领域
本发明属于超声速飞机飞行噪声领域,具体涉及一种用于音爆传播计算的多波系同步推进波形参数方法。
背景技术
音爆是飞行器超声速飞行时,由于头波和结尾激波间的扰动传播到地面而形成的压力突跃信号,可引起地表结构、建筑物、人和动物耳膜等受迫振动,对地面环境、社区人群、生物群落产生非预期影响。在协和飞机上天后,美国禁止了陆地超声速民航活动,音爆成为超声速商用飞行的一大障碍。
音爆的形成包括两个关键环节,一是飞行器附近的近场流动扰动,在拓扑结构上表现为头部激波、机翼激波、肩部膨胀波、结尾激波、喷流激波等波系结构的相互作用;二是音爆的远场传播,流动拓扑结构在由飞行器附近向远场传播过程中发生合并、消融与畸变,到达地面形成压力突跃的短时信号——音爆。相应于这两个环节,在音爆预测与研究中,首先通过风洞试验或者流场计算,获得可用压力扰动表征的近场信号;然后,通过传播模型进行音爆传播计算,获得最终的地面音爆信号。目前,音爆传播计算途径主要有两类,一类是经典的、基于改进线性理论的传播计算,另外一类就是基于发展方程的非线性传播计算,二者均采用几何声学理论进行音爆传播的建模处理。其中,基于改进线性理论的音爆传播计算方法主要包括源自NASA CR-1299报告中Haye等发展的F函数方法和NASA TN D-6832报告中Thomas的波形参数方法,由于传统的波形参数法物理模型简单、计算效率高,在音爆的快速分析与低音爆设计中获得广泛应用。
传统的波形参数法采用折线来逼近音爆的波形信号,并以折线上各线段的区间长度、斜率以及压力跃升三个参数来实现音爆信号波形的完整描述。波形信号参数化后,根据几何声学理论和音爆信号在分层大气中的演化机制建立波形参数的时间变化率方程,以此复现音爆的传播现象。传统的波形参数法的特色是:(1)理论模型比较直观,由一组具有解析解的常微分方程描述;(2)该方法直接外插音爆的过压信号,可以直接从风洞试验或者数值计算结果进行外插,使用更为方便;(3)无需显式的面积平衡来确定激波位置,激波演化更为简洁。
传统的波形参数法主要工作包括音爆射线管追踪计算和波形参数演化计算,其再现音爆传播过程是:从航迹处开始,将音爆射线和信号波形按照离散时间或者空间步向地面推进,直至给定的高度为止。如果当前时间推进步中波形上某段区间长度表达式的因子出现负值,则意味着在该推进步的中间某个时刻出现了新激波或者激波融合,此时,以当前时间步长进行传播计算则出错。这种情况下,需要重新确定更短的新时间步长,使其刚好在推进末尾时刻信号波形上某一段形成新激波或者激波融合,并在必要时人工调节波形上其它特定段的斜率或者区间长度参数以确保整个波形上有且仅有一段形成新激波或者激波融合。新时间步长确定后,以新时间步长进行推进计算,处理新激波或者激波融合。因此,传统的波形参数法在每一推进步中只处理波形上的一处新激波或者激波融合。
传统的波形参数法的典型计算流程,除了初始化(音爆射线管初始位置计算、初始信号参数化)和音爆射线管追踪计算等步骤外,核心工作包括确定步长参数T和计算信号各段区间长度表达式的因子F1,i和F2,i,以及音爆传播两大工作。
传统的波形参数法采用两步双循环方式来确定步长参数T和计算信号各段区间长度表达式的因子F1,i和F2,i。第一步,计算信号波形各段的区间长度表达式第一因子F1,i和确定步长参数T。以当前步长参数T为输入,采用顺次推进方式,计算信号波形上每一段的区间长度表达式的第一因子F1,i,如果在某一段上该因子小于零(F1,k<0),则表明在该推进步内此段将出现新激波或者激波融合,由该段的F1,k(Tnew)=0来确定新的步长参数Tnew并赋给当前步长参数T,将与该段斜率相同的其他段的线段斜率减小0.01%,确保折线上只有一段的区间长度表达式第一因子为零(F1,k=0),然后回到本步开头重新处理,直到所有段的区间长度表达式因子F1,i中最多只有一个为零,其余大于零为止;第二步,计算各段区间长度表达式第二因子F2,i和确定步长参数T。在第一步的基础上,以当前步长参数T为输入,采用顺次推进方式,计算信号波形上每一段的区间长度表达式第二因子F2,i。如果某一段该因子小于零(F2,q<0),则由该段的F2,q(Tnew)=0来确定新的步长参数Tnew并赋给当前步长参数T,同时将其他段上第二因子为零的区间长度增大0.01%,确保折线上只有一段的第二因子为零(F2,q=0),回到第一步重新计算各段区间长度第一因子F1,i和确定步长参数T。采用这样的两步双循环方式确定步长参数T和波形信号各段区间长度因子F1,i、F2,i,可以确保在每一推进步中波形上最多只有一段收缩形成新激波或者激波融合。
在采用上面的两步双循环方式确定步长参数T和波形信号各段区间长度因子F1,i、F2,i后,传统的波形参数法就以步长参数T,进行音爆信号参数的传播处理。如果所有段的区间长度因子均大于零,则传播推进步后不会出现激波,波形所有段的参数按正常传播解计算,得到新的波形;如果出现某一段区间长度因子为零的情况,即F1,k=0或者F2,q=0,则该段收缩为新激波或者激波融合,产生的压力突跃计入其后面邻接段的压力跃升,其前后段的其余参数都采用正常传播方式计算,此时,传播推进步则得到只有一段信号收缩后的新波形。
传统的波形参数法采用四重循环进行波形传播计算,内部双循环用于步长参数T的检查、更新和各段区间长度因子F1,i、F2,i的计算,外部双循环分别用于射线管跟踪单个推进步内的信号传播和整个射线管在大气中的传播。传统的波形参数法在处理多波系同步推进的物理过程时,通过人工的斜率调整或者区间长度调整将多波系进行了异步推进处理。
传统的波形参数法存在以下不足:一是采用两步双循环来确定步长参数T,只要出现新的步长参数Tnew就要回到第一步循环迭代进行步长参数T和各段区间长度因子F1,i、F2,i的计算,一方面浪费了计算资源,另一方面从逻辑上存在死循环的潜在风险;二是在处理多波系同步推进的物理过程时,每一推进步只允许处理一道激波而人为给予其他同斜率段小扰动(0.01%)或者区间长度小扰动(0.01%),虽然避免了原所有斜率相同段区间长度因子等于0的问题,但也致使相关段的区间长度表达式第一因子接近于0,从而使得相应段的斜率急剧增大。这既没有模拟真实的物理过程,也存在数值发散的风险。
发明内容
本发明所要解决的技术问题是提供一种用于音爆传播计算的多波系同步推进波形参数方法。
本发明的用于音爆传播计算的多波系同步推进波形参数方法包括以下步骤:
a.初始化;飞行条件、初始信号、终止高度输入,从飞行轨迹开始采用几何声学方法进行音爆射线管计算,直至初始信号所在高度,置为当前位置,TS=1,设定射线管步进时间间隔Δt;对于给定距纵轴r1处按过压比和距离描述的初始信号,需将初始信号
Figure BDA0002540477120000051
转换为时域过压信号
Figure BDA0002540477120000052
并根据需要把信号分段和线性近似,形成n个线性段δpi~τ,i=1…n构成的折线型波形信号,置当前波形总段数N=n;波形上每一段δpi(τ)采用该段斜率
Figure BDA0002540477120000053
该段区间长度
Figure BDA0002540477120000054
以及该段与前段之间的突跃
Figure BDA0002540477120000055
来描述;
b.射线管单步跟踪;从当前位置开始,按照给定的射线管步进时间间隔Δt,计算射线管推进的新位置,获得相应射线管参数、大气参数及音爆波形参数传播微分方程系数C1,C2,设置音爆信号传播推进的时间步长dt=Δt;
c.射线管单步步进中波形传播的预期步长参数计算;计算音爆信号传播推进的步长参数
Figure BDA0002540477120000056
并将新步长参数置零,即Tnew=0;
d.步长参数T符合性的顺序递进检验及更新、区间长度因子F1,i和F2,i的计算,包括以下步骤:
d1.以步长参数T和当前波形参数,计算波形第一段区间长度表达式的第一因子
Figure BDA0002540477120000057
若F1,1<0,则按F1,1(Tnew)=0计算新步长参数Tnew,更新步长参数T=Tnew;
d2.以步长参数T和当前波形参数,按从前至后顺序递进的方式对波形信号上的各段进行处理,首先计算波形信号上下一段的表达式第一区间长度因子
Figure BDA0002540477120000058
若F1,i+1<0,则按F1,i+1(Tnew)=0计算新步长参数Tnew,更新步长参数T=Tnew;然后计算本段的区间长度表达式第二因子:
Figure BDA0002540477120000059
若F2,i(T)<0则按F2,i(Tnew)=0计算[0,T]内的新步长参数Tnew,更新步长参数T=Tnew;
d3.波形上所有段步长参数T符合性检验和更新完毕后,如果Tnew存在,则新时间步长
Figure BDA0002540477120000061
计算信号波形各段的区间长度因子F1,i和F2,i;否则按照预期的时间步长dt传播,即Tnew=T,Δt′=T-1(Tnew)=dt;
e.多波系同步推进波形参数传播计算;开辟线性存储结构以存储新的波形参数mi、Δpi、λi;按照确定的步长参数T,对波形的每一段进行音爆传播处理,包含以下过程,首先,置新波形的段数inewseg为零,即inewseg=0;然后,对于波形的每一段i,置新间断压力跃升变量shockdp为零,即shockdp=0,并按以下步骤进行传播处理:
e1.如果F1,i=0,则该段在Δt′推进后区间长度收缩为零,形成新生激波,其压力跃升为:
Figure BDA0002540477120000062
e2.如果F2,i=0,则该段在Δt′推进后区间长度收缩零,出现激波融合,其压力跃升为:
Figure BDA0002540477120000063
e3.如果F1,i>0,F2,i>0,则经Δt′推进后,该段传播后区间长度大于零,波形正常传播,此时,增加新波形的段数,即inewseg=inewseg+1,并记录新波形参数为
Figure BDA0002540477120000064
然后将新间断压力跃升变量置为零,即shockdp=0;
e4.波形所有的信号段传播处理完毕后,将新的波形作为当前波形供后续传播使用,即,
Figure BDA0002540477120000065
,i=1…,N=inewseg;
f.判断波形传播是否达到射线管单步推进的末端位置;计算射线管TS推进步内信号的剩余传播时间dt=dt-Δt′,如果dt>0,则将dt作为音爆下一步推进的预期时间步长,进入步骤d,继续完成射线管TS推进步内的信号传播;如果dt=0,表明射线管TS推进步完成,进入步骤g;
g.判断射线管跟踪循环是否终结;判断是否已经达到终止高度,终止高度默认为地面;若未达到终止高度则将计算的射线管新位置作为当前射线管位置,TS=TS+1,并转到步骤b进行下一射线管单步传播计算;否则,将信号参数还原为时域信号,得到最终结果;若终止位置为地面,则需乘以地面反射因子。
所述的步骤d可替换为以下步骤:
d1.采用波形上斜率最大的段进行步长参数T关于区间长度第一因子的符合性检验,必要时更新步长参数T;即,依据波形参数法中波形各段区间长度表达式的第一因子
Figure BDA0002540477120000071
计算最高斜率段的区间长度第一因子
Figure BDA0002540477120000072
如果F1,k<0,则获得新的步长参数
Figure BDA0002540477120000073
d2.对波形上所有的段,按从前至后的顺序,以区间长度表达式的第二个因子检验步长参数T的符合性及其更新,即,以当前步长参数T为输入,计算其区间长度表达式的第二个因子:
Figure BDA0002540477120000074
如果F2,i(T)<0,采用数值方法计算[0,T]区间使得F2,i(Tnew)=0的新步长参数Tnew,并更新当前步长参数T=Tnew,为后续各段中步长参数T关于区间长度第二因子的符合性检验提供输入;
d3.对确定的步长参数T,按照
Figure BDA0002540477120000075
以及F2,i(T)的表达式计算信号波形上需要更新的各段区间长度因子;如果Tnew存在,则计算推进到新激波生成或者激波融合所需时间步长
Figure BDA0002540477120000081
否则按照预期的时间步长dt传播,即Tnew=T,Δt′=T-1(Tnew)=dt。
本发明的用于音爆传播计算的多波系同步推进波形参数方法基于音爆射线管追踪和Thomas的波形参数模型,采用内外双循环流程进行传播计算,运用顺序递进方式确定内循环的波形传播步长参数T,确保波形上所有段的区间长度因子F1,i和F2,i具有非负性,保留了同步形成多道激波或者激波融合的可能性(即多个区间长度因子为零)。在波形参数传播计算中,采用区间长度因子F1,i和F2,i依次判别波形信号各段在传播步后的模态,分别按照新激波生成、激波融合和正常传播三种模态进行激波跃升、波形参数传播计算处理。采用线性存储结构记录新波形参数,实现了新时间步长确定过程的顺序递进计算和多波系的同步推进。
本发明的用于音爆传播计算的多波系同步推进波形参数方法的步骤要点为:a.初始化;b.射线管单步跟踪(外循环起始);c.计算射线管单步推进中波形传播预期步长参数T(内循环起始);d.步长参数T符合性检验及其更新、各区间长度因子F1,i和F2,i的递进计算;e.多波系同步推进的波形参数传播计算;f.判断波形传播是否达到射线管末端(内循环终止判断);g.判断射线管是否达到终止位置(外循环终止判断)。
本发明的用于音爆传播计算的多波系同步推进波形参数方法中的步骤d,采用顺序递进方式检验和更新步长参数T,确保了有限步内实现步长参数T的检验和更新,避免了传统的波形参数法采用双循环逻辑进行检验和更新步长参数T方式可能带来逻辑死循环的风险;步骤d允许音爆传播中波形上多段出现区间长度第一因子或者第二因子等于零,允许多激波同步生成;步骤d无需人工调节任何段的斜率、任何段的区间长度,避免了传统的波形参数法对于多道激波同步生成或者融合段的斜率或者区间长度人工调节需求。
本发明的用于音爆传播计算的多波系同步推进波形参数方法中的步骤e,对波形上每一段的传播分为新生激波、激波融合和正常传播三种模态进行处理。每当遇到新生激波、激波融合时,只需计算新间断的压力跃升参数,等待插入其后面的正常传播段。允许处理多激波的同步推进,降低了异步处理中的斜率剧增。
本发明的用于音爆传播计算的多波系同步推进波形参数方法处理具有多波系同步生成特征的音爆传播物理进程时,无需人工调整波形信号任何段的斜率或者区间长度,可以在一个时间步内完成多道激波生成或者激波融合,克服了传统的波形参数法将同步推进转化为异步推进而带来的激波段波形斜率迅速增大和区间长度急剧减小的问题,既可提高计算效率,又可提高计算的稳定性。
附图说明
图1为本发明的用于音爆传播计算的多波系同步推进波形参数方法的计算流程图。
图2为航天飞机助推器再入大气音爆外插算例验证结果对比图;
图2中,△曲线为初始信号,○曲线为本发明的用于音爆传播计算的多波系同步推进波形参数方法的外插结果,*曲线为传统的波形参数法的外插结果。
图3a为多波系同步推进传播过程的初始信号;
图3a中,
Figure BDA0002540477120000091
曲线为具有三段等斜率的初始信号。
图3b为多波系同步推进物理过程的传统的波形参数法处理过程中各段压力跃升图;
图3b中,+曲线为传统的波形参数法在第8步推进形成的各段压力跃升;□曲线为传统的波形参数法在第9步推进形成的各段压力跃升;×曲线为传统的波形参数法在第10步推形成的各段压力跃升;
Figure BDA0002540477120000101
曲线为传统的波形参数法在第11步推进形成的各段压力跃升。
图3c为多波系同步推进物理过程的本发明的用于音爆传播计算的多波系同步推进波形参数方法处理过程中各段压力跃升图;
图3c中,◇曲线为本发明的用于音爆传播计算的多波系同步推进波形参数方法在第8步推进形成的各段压力跃升;
Figure BDA0002540477120000102
曲线为本发明的用于音爆传播计算的多波系同步推进波形参数方法在第9步推进形成的各段压力跃升。
具体实施方式
下面结合附图和实施例详细说明本发明。
本发明的用于音爆传播计算的多波系同步推进波形参数方法的流程图见图1。
本发明的用于音爆传播计算的多波系同步推进波形参数方法包括以下步骤:
a.初始化;飞行条件、初始信号、终止高度输入,从飞行轨迹开始进行音爆射线管计算,直至初始信号所在高度,置为当前位置,TS=1,设定射线管步进时间间隔Δt;对于给定距纵轴r1处按过压比和距离描述的初始信号,需将初始信号
Figure BDA0002540477120000103
转换为时域过压信号
Figure BDA0002540477120000104
并根据需要把信号分段和线性近似,形成n个线性段δpi~τ,i=1…n构成的折线型波形信号,置当前波形总段数N=n;波形上每一段δpi(τ)采用该段斜率
Figure BDA0002540477120000105
该段区间长度
Figure BDA0002540477120000106
以及该段与前段之间的突跃
Figure BDA0002540477120000107
来描述;
b.射线管单步跟踪;从当前位置开始,按照给定的射线管步进时间间隔Δt,计算射线管推进的新位置,获得相应射线管参数、大气参数及音爆波形参数传播微分方程系数C1,C2,设置音爆信号传播推进的时间步长dt=Δt;
c.射线管单步步进中波形传播的预期步长参数计算;计算音爆信号传播推进的步长参数
Figure BDA0002540477120000111
并将新步长参数置零,即Tnew=0;
d.步长参数T符合性的顺序递进检验及更新、区间长度因子F1,i和F2,i的计算,包括以下步骤:
d1.以步长参数T和当前波形参数,计算波形第一段的区间长度表达式第一因子
Figure BDA0002540477120000112
若F1,1<0,则按F1,1(Tnew)=0计算新步长参数Tnew,更新步长参数T=Tnew;
d2.以步长参数T和当前波形参数,按从前至后顺序递进的方式对波形信号上的各段进行处理,首先,计算波形信号下一段的区间长度表达式第一因子
Figure BDA0002540477120000113
若F1,i+1<0,则按F1,i+1(Tnew)=0计算新步长参数Tnew,更新步长参数T=Tnew;然后,计算本段的区间长度表达式第二因子:
Figure BDA0002540477120000114
若F2,i(T)<0则按F2,i(Tnew)=0计算[0,T]内的新步长参数Tnew,更新步长参数T=Tnew;
d3.波形上所有段步长参数T符合性检验和更新完毕后,如果Tnew存在,则新时间步长
Figure BDA0002540477120000115
计算信号波形各段的区间长度因子F1,i和F2,i;否则按照预期的时间步长dt传播,即Tnew=T,Δt′=T-1(Tnew)=dt;
e.多波系同步推进波形参数传播计算;开辟线性存储结构以存储新的波形参数mi、Δpi、λi;按照确定的步长参数T,对波形的每一段进行音爆传播处理,包含以下过程,首先,置新波形的段数inewseg为零,即inewseg=0;然后,对于波形的每一段i,置新间断压力跃升变量shockdp为零,即shockdp=0,并按以下步骤进行传播处理:
e1.如果F1,i=0,则该段在Δt′推进后区间长度收缩为零,形成新生激波,其压力跃升为:
Figure BDA0002540477120000121
e2.如果F2,i=0,则该段在Δt′推进后区间长度收缩零,出现激波融合,其压力跃升为:
Figure BDA0002540477120000122
e3.如果F1,i>0,F2,i>0,则经Δt′推进后,该段传播后区间长度大于零,波形正常传播,此时,增加新波形的段数,即inewseg=inewseg+1,并记录新波形参数为
Figure BDA0002540477120000123
然后将新间断压力跃升变量置为零,即shockdp=0;
e4.波形所有的信号段传播处理完毕后,将新的波形作为当前波形供后续传播使用,即,
Figure BDA0002540477120000124
,i=1…,N=inewseg;
f.判断波形传播是否达到射线管单步推进的末端位置;计算射线管TS推进步内信号的剩余传播时间dt=dt-Δt′,如果dt>0,则将dt作为音爆下一步推进的预期时间步长,进入步骤d,继续完成射线管TS推进步内的信号传播;如果dt=0,表明射线管TS推进步完成,进入步骤g;
g.判断射线管跟踪循环是否终结;判断是否已经达到终止高度,终止高度默认为地面;若未达到终止高度则将计算的射线管新位置作为当前射线管位置,TS=TS+1,并转到步骤b进行下一射线管单步传播计算;否则,将信号参数还原为时域信号,得到最终结果;若终止位置为地面,则需乘以地面反射因子。
所述的步骤d可替换为以下步骤:
d1.采用波形上斜率最大的段进行步长参数T关于区间长度第一因子的符合性检验,必要时更新步长参数T;即,依据波形参数法中波形各段区间长度表达式的第一因子
Figure BDA0002540477120000131
计算最高斜率段的区间长度第一因子
Figure BDA0002540477120000132
如果F1,k<0,则获得新的步长参数
Figure BDA0002540477120000133
T=Tnew;
d2.对波形上所有的段,按从前至后的顺序,以区间长度表达式的第二个因子检验步长参数T的符合性及其更新,即,以当前步长参数T为输入,计算其区间长度表达式的第二个因子:
Figure BDA0002540477120000134
如果F2,i(T)<0,采用数值方法计算[0,T]区间使得F2,i(Tnew)=0的新步长参数Tnew,并更新当前步长参数T=Tnew,为后续各段中步长参数T关于区间长度第二因子的符合性检验提供输入;
d3.对确定的步长参数T,按照
Figure BDA0002540477120000135
以及F2,i(T)的表达式计算信号波形上需要更新的各段区间长度因子;如果Tnew存在,则计算推进到新激波生成或者激波融合所需时间步长
Figure BDA0002540477120000136
否则按照预期的时间步长dt传播,即Tnew=T,Δt′=T-1(Tnew)=dt。
实施例1
力p0g=2116.25psf,地面反射因子取1.9。助推器长度L=253ft,以Thomas报告给出的航天飞机助推器再入大气的音爆数据进行了本发明的用于音爆传播计算的多波系同步推进波形参数方法的验证实施。
助推器的飞行条件为:飞行马赫数M=1.2、飞行高度H=50400ft、航迹角γ=12.75°、指向角ψ=356.5°、马赫数变化率
Figure BDA0002540477120000141
航迹角变化率
Figure BDA0002540477120000142
指向角变化率
Figure BDA0002540477120000143
地面大气压器风洞试验模型长度Lmodel=10ft,风洞试验时音爆测量位置(距速度轴的距离与飞行器长度比值)
Figure BDA0002540477120000144
测量的方位角φ=47.0°,初始信号见图2中的△曲线。
按照本发明的用于音爆传播计算的多波系同步推进波形参数方法对信号进行参数化。然后,以0.2s的射线管跟踪步长(Δt=0.2),循环处理射线管推进,循环计算确定步长参数Tnew、各信号段的区间因子和多波系传播,最后得到如图2中的○曲线所示的音爆传播到地面的终止结果。
图2中的*曲线给出了传统的波形参数法的外插结果。
由图2可知,对于单步单波系推进情况,本发明的用于音爆传播计算的多波系同步推进波形参数方法的音爆外插结果与传统的波形参数法一致。
实施例2
本实施例的飞行条件与实施例1相同,但初始信号为含有三段等斜率的过压信号,其中,第一、三、十三段斜率相同,见图3a。分别按照Thomas波形参数法和本发明的用于音爆传播计算的多波系同步推进波形参数方法进行了图3a初始信号的传播计算。
图3b给出了传统的波形参数法中激波异步形成过程的各段压力跃升结果。可见,采用传统的波形参数法,在传播的第8步推进后各段压力跃升仍然为零,没有激波形成;第9步推进使得第三段压力突跃大于零(□曲线),表明初始信号第三段收缩为激波,第9步开始产生新激波;第10步推进使得第十二段出现了压力突跃(×曲线),表明初始信号第十三段收缩为激波;第11步推进又使得第一段出现压力突跃(
Figure BDA0002540477120000151
曲线),表明初始信号第一段收缩为激波,至此,三道激波形成。从斜率数据上看,第9步使得第三段形成新激波,第一和十三段的斜率从103增大到了107量级;第10步,使得初始信号第十三段形成新激波,而第一段的斜率则达到了1011量级,第11步使得初始信号第一段形成新激波,新信号中各段斜率均为负值。
图3c为本发明的用于音爆传播计算的多波系同步推进波形参数方法的各段压力跃升结果,◇曲线和
Figure BDA0002540477120000152
曲线分别为本发明的用于音爆传播计算的多波系同步推进波形参数方法在第8、9步推进结果。可见,采用本发明的用于音爆传播计算的多波系同步推进波形参数方法,在第9步,信号的第一、三、十三段同时形成前、中、后三道激波。
由图3可知,本发明的用于音爆传播计算的多波系同步推进波形参数方法复现了多激波同步推进的传播机制,避免了传统的波形参数法异步处理激波同步形成过程中信号斜率急剧增大的问题。

Claims (2)

1.一种用于音爆传播计算的多波系同步推进波形参数方法,其特征在于,包括以下步骤:
a.初始化;飞行条件、初始信号、终止高度输入,从飞行轨迹开始采用几何声学方法进行音爆射线管计算,直至初始信号所在高度,置为当前位置,TS=1,设定射线管步进时间间隔Δt;对于距纵轴r1处按过压比和距离描述的初始信号,需将初始信号
Figure FDA0002540477110000011
转换为时域过压信号
Figure FDA0002540477110000012
并根据需要把信号分段和线性近似,形成n个线性段δpi~τ,i=1…n构成的折线型波形信号,置当前波形总段数N=n;波形上每一段δpi(τ)采用该段斜率
Figure FDA0002540477110000013
该段区间长度
Figure FDA0002540477110000014
以及该段与前段之间的突跃
Figure FDA0002540477110000015
来描述;
b.射线管单步跟踪;从当前位置开始,按照给定的射线管步进时间间隔Δt,计算射线管推进的新位置,获得相应射线管参数、大气参数及音爆波形参数传播微分方程系数C1,C2,设置音爆信号传播推进的时间步长dt=Δt;
c.射线管单步步进中波形传播的预期步长参数计算;计算音爆信号传播推进的步长参数
Figure FDA0002540477110000016
并将新步长参数置零,即Tnew=0;
d.步长参数T符合性的顺序递进检验及更新、区间长度因子F1,i和F2,i的计算,包括以下步骤:
d1.以步长参数T和当前波形参数,计算波形第一段的区间长度表达式第一因子
Figure FDA0002540477110000017
若F1,1<0,则按F1,1(Tnew)=0计算新步长参数Tnew,更新步长参数T=Tnew;
d2.以步长参数T和当前波形参数,按从前至后顺序递进的方式对波形信号上的各段进行处理,首先计算波形信号下一段的区间长度表达式第一因子
Figure FDA0002540477110000018
若F1,i+1<0,则按F1,i+1(Tnew)=0计算新步长参数Tnew,更新步长参数T=Tnew;然后计算本段的区间长度表达式第二因子:
Figure FDA0002540477110000021
若F2,i(T)<0则按F2,i(Tnew)=0计算[0,T]内的新步长参数Tnew,更新步长参数T=Tnew;
d3.波形上所有段步长参数T符合性检验和更新完毕后,如果Tnew存在,则新时间步长
Figure FDA0002540477110000022
计算信号波形各段的区间长度因子F1,i和F2,i;否则按照预期的时间步长dt传播,即Tnew=T,Δt′=T-1(Tnew)=dt;
e.多波系同步推进波形参数传播计算;开辟线性存储结构以存储新的波形参数mi、Δpi、λi;按照确定的步长参数T,对波形的每一段进行音爆传播处理,包含以下过程,首先,置新波形的段数inewseg为零,即inewseg=0;然后,对于波形的每一段i,置新间断压力跃升变量shockdp为零,即shockdp=0,并按以下步骤进行传播处理:
e1.如果F1,i=0,则该段在Δt′推进后区间长度收缩为零,形成新生激波,其压力跃升为:
Figure FDA0002540477110000023
e2.如果F2,i=0,则该段在Δt′推进后区间长度收缩零,出现激波融合,其压力跃升为:
Figure FDA0002540477110000024
e3.如果F1,i>0,F2,i>0,则经Δt′推进后,该段传播后区间长度大于零,波形正常传播,此时,增加新波形的段数,即inewseg=inewseg+1,并记录新波形参数为
Figure FDA0002540477110000031
然后将新间断压力跃升变量置为零,即shockdp=0;
e4.波形所有的信号段传播处理完毕后,将新的波形作为当前波形供后续传播使用,即,
Figure FDA0002540477110000032
N=inewseg;
f.判断波形传播是否达到射线管单步推进的末端位置;计算射线管TS推进步内信号的剩余传播时间dt=dt-Δt′,如果dt>0,则将dt作为音爆下一步推进的预期时间步长,进入步骤d,继续完成射线管TS推进步内的信号传播;如果dt=0,表明射线管TS推进步完成,进入步骤g;
g.判断射线管跟踪循环是否终结;判断是否已经达到终止高度,终止高度默认为地面;若未达到终止高度则将计算的射线管新位置作为当前射线管位置,TS=TS+1,并转到步骤b进行下一射线管单步传播计算;否则,将信号参数还原为时域信号,得到最终结果;若终止位置为地面,则需乘以地面反射因子。
2.根据权利要求1所述的用于音爆传播计算的多波系同步推进波形参数方法,其特征在于,所述的步骤d替换为以下步骤:
d1.采用波形上斜率最大的段进行步长参数T关于区间长度第一因子的符合性检验,必要时更新步长参数T;即,依据波形参数法中波形各段区间长度表达式的第一因子
Figure FDA0002540477110000033
计算最高斜率段的区间长度第一因子
Figure FDA0002540477110000034
如果F1,k<0,则获得新的步长参数
Figure FDA0002540477110000041
T=Tnew;
d2.对波形上所有的段,按从前至后的顺序,以区间长度表达式的第二个因子检验步长参数T的符合性及其更新,即,以当前步长参数T为输入,计算其区间长度表达式的第二个因子:
Figure FDA0002540477110000042
如果F2,i(T)<0,采用数值方法计算[0,T]区间使得F2,i(Tnew)=0的新步长参数Tnew,并更新当前步长参数T=Tnew,为后续各段中步长参数T关于区间长度第二因子的符合性检验提供输入;
d3.对确定的步长参数T,按照
Figure FDA0002540477110000043
以及F2,i(T)的表达式计算信号波形上需要更新的各段区间长度因子;如果Tnew存在,则计算推进到新激波生成或者激波融合所需时间步长
Figure FDA0002540477110000044
否则按照预期的时间步长dt传播,即Tnew=T,Δt′=T-1(Tnew)=dt。
CN202010545227.XA 2020-06-16 2020-06-16 一种用于音爆传播计算的多波系同步推进波形参数方法 Active CN111695264B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010545227.XA CN111695264B (zh) 2020-06-16 2020-06-16 一种用于音爆传播计算的多波系同步推进波形参数方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010545227.XA CN111695264B (zh) 2020-06-16 2020-06-16 一种用于音爆传播计算的多波系同步推进波形参数方法

Publications (2)

Publication Number Publication Date
CN111695264A true CN111695264A (zh) 2020-09-22
CN111695264B CN111695264B (zh) 2023-03-03

Family

ID=72481036

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010545227.XA Active CN111695264B (zh) 2020-06-16 2020-06-16 一种用于音爆传播计算的多波系同步推进波形参数方法

Country Status (1)

Country Link
CN (1) CN111695264B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112924131A (zh) * 2021-01-28 2021-06-08 西北工业大学 考虑大气边界层湍流效应的远场声爆预测方法

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB1205865A (en) * 1968-06-24 1970-09-23 Heinz Preuss Apparatus for reducing sonic boom shockwave propagation in the direction towards the earth caused by an aircraft in supersonic flight
US4634970A (en) * 1983-12-30 1987-01-06 Norland Corporation Digital waveform processing oscilloscope with distributed data multiple plane display system
CA2245769A1 (en) * 1996-02-13 1997-08-21 James T. Sears Tactiley-guided, voice-output reading apparatus
US20040164984A1 (en) * 2003-02-25 2004-08-26 Pickerd John J. Method of constraints control for oscilloscope vertical subsection and display parameters
JP2005178491A (ja) * 2003-12-18 2005-07-07 Japan Aerospace Exploration Agency 超音速航空機の胴体形状の決定方法および胴体前胴部形状
US8145366B1 (en) * 2008-06-13 2012-03-27 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Real-time, interactive sonic boom display
TW201506915A (zh) * 2013-08-15 2015-02-16 Jack Kuei-Chung Shih 針對空間中多音源進行萃取出單一音源的方法及裝置
CN204925207U (zh) * 2015-09-21 2015-12-30 成都华兴智造科技有限公司 一种应用于电流电位测量的数字无线示波器
CN108170878A (zh) * 2016-12-08 2018-06-15 中国航空工业集团公司沈阳空气动力研究所 一种超音速飞行器音爆预测方法
US10209122B1 (en) * 2017-10-31 2019-02-19 Honeywell International Inc. Systems and methods for generating avionic displays including forecast overpressure event symbology
CN110132528A (zh) * 2019-06-27 2019-08-16 中国空气动力研究与发展中心高速空气动力研究所 一种暂冲式超声速风洞音爆测量试验装置及测定方法
US20200180780A1 (en) * 2018-12-07 2020-06-11 U.S.A. As Represented By The Administrator Of The National Aeronautics And Space Administration System and Method for Performing Multi-Point, Full-Mission Sonic Boom Prediction
CN112924131A (zh) * 2021-01-28 2021-06-08 西北工业大学 考虑大气边界层湍流效应的远场声爆预测方法

Patent Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB1205865A (en) * 1968-06-24 1970-09-23 Heinz Preuss Apparatus for reducing sonic boom shockwave propagation in the direction towards the earth caused by an aircraft in supersonic flight
US4634970A (en) * 1983-12-30 1987-01-06 Norland Corporation Digital waveform processing oscilloscope with distributed data multiple plane display system
CA2245769A1 (en) * 1996-02-13 1997-08-21 James T. Sears Tactiley-guided, voice-output reading apparatus
US20040164984A1 (en) * 2003-02-25 2004-08-26 Pickerd John J. Method of constraints control for oscilloscope vertical subsection and display parameters
JP2005178491A (ja) * 2003-12-18 2005-07-07 Japan Aerospace Exploration Agency 超音速航空機の胴体形状の決定方法および胴体前胴部形状
US8145366B1 (en) * 2008-06-13 2012-03-27 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Real-time, interactive sonic boom display
TW201506915A (zh) * 2013-08-15 2015-02-16 Jack Kuei-Chung Shih 針對空間中多音源進行萃取出單一音源的方法及裝置
CN204925207U (zh) * 2015-09-21 2015-12-30 成都华兴智造科技有限公司 一种应用于电流电位测量的数字无线示波器
CN108170878A (zh) * 2016-12-08 2018-06-15 中国航空工业集团公司沈阳空气动力研究所 一种超音速飞行器音爆预测方法
US10209122B1 (en) * 2017-10-31 2019-02-19 Honeywell International Inc. Systems and methods for generating avionic displays including forecast overpressure event symbology
US20200180780A1 (en) * 2018-12-07 2020-06-11 U.S.A. As Represented By The Administrator Of The National Aeronautics And Space Administration System and Method for Performing Multi-Point, Full-Mission Sonic Boom Prediction
CN110132528A (zh) * 2019-06-27 2019-08-16 中国空气动力研究与发展中心高速空气动力研究所 一种暂冲式超声速风洞音爆测量试验装置及测定方法
CN112924131A (zh) * 2021-01-28 2021-06-08 西北工业大学 考虑大气边界层湍流效应的远场声爆预测方法

Non-Patent Citations (8)

* Cited by examiner, † Cited by third party
Title
F.ALAUZET ET AL.: "High-order sonic boom modeling based on adaptive methods", 《JOURNAL OF COMPUTATIONAL PHYSICS》 *
KENNETH J. PLOTKIN: "State of the art of sonic boom modeling", 《THE JOURNAL OF THE ACOUSTICAL SOCIETY OF AMERICA》 *
MAGLIERI, DOMENIC J. ET AL.: "Focused and Steady-State Characteristics of Shaped Sonic Boom Signatures: Prediction and Analysis", 《NTRS-NASA TECHNICAL REPORTS SERVER》 *
乔建领 等: "基于广义Burgers方程的超声速客机远场", 《空气动力学学报》 *
冷岩 等: "均匀各向同性大气湍流对声爆传播特性的影响", 《航空学报》 *
刘刚 等: "超声速飞行器声爆/气动力综合设计技术研究", 《空气动力学学报》 *
刘少伟等: "考虑声爆特性的超声速客机气动优化设计", 《西北工业大学学报》 *
王刚 等: "典型标模音爆的数值预测与分析", 《航空学报》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112924131A (zh) * 2021-01-28 2021-06-08 西北工业大学 考虑大气边界层湍流效应的远场声爆预测方法

Also Published As

Publication number Publication date
CN111695264B (zh) 2023-03-03

Similar Documents

Publication Publication Date Title
CN104538028B (zh) 一种基于深度长短期记忆循环神经网络的连续语音识别方法
Wintzer et al. Under-track CFD-based shape optimization for a low-boom demonstrator concept
Marchiano et al. Numerical simulation of shock wave focusing at fold caustics, with application to sonic boom
Han et al. A variable-fidelity modeling method for aero-loads prediction
Ordaz et al. Full-carpet design of a low-boom demonstrator concept
CN111695264B (zh) 一种用于音爆传播计算的多波系同步推进波形参数方法
KR102349898B1 (ko) 적응적 학습률로 뉴럴 네트워크를 학습하는 방법 및 장치, 이를 이용한 테스트 방법 및 장치
Rogers et al. CFD simulations of the space launch system ascent aerodynamics and booster separation
CN108170878B (zh) 一种超音速飞行器音爆预测方法
Ueno et al. Multi-fidelity low-boom design based on near-field pressure signature
Campbell et al. Building a practical natural laminar flow design capability
Nakahashi et al. Numerical Simulations on Separation of Scaled Supersonic Experimental Airplane from Rocket Booster at Supersonic Speed
CN111695195B (zh) 一种基于长短时间记忆网络的空间物理运动体建模方法
Choi et al. Numerical and mesh resolution requirements for accurate sonic boom prediction
Kirz Surrogate based shape optimization of a low boom fuselage wing configuration
Reddy et al. Achieving Quieter Supersonic Flight Through Outer-Mold Line Modifications: An Optimization Study
CN105487110B (zh) 一种基于透射方程的各向异性参数反演方法
CN107869993A (zh) 基于自适应迭代粒子滤波的小卫星姿态估计方法
Haas et al. A multi-shock inverse design method for low-boom supersonic aircraft
Epstein et al. A new efficient technology of aerodynamic design based on CFD driven optimization
CN106601234A (zh) 一种面向货物分拣的地名语音建模系统的实现方法
Zlenko et al. Method of optimal aerodynamic design of the nacelle for the main propulsion system with a high bypass ratio
Housman et al. LAVA Simulations for the First AIAA Sonic Boom Prediction Workshop
KR101061188B1 (ko) 항공기 면적곡선 분석 방법
Minelli et al. Advanced optimization approach for supersonic low-boom design

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