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

本发明属于超声速飞机飞行噪声领域，具体涉及一种用于音爆传播计算的多波系同步推进波形参数方法。

背景技术：

音爆是飞行器超声速飞行时，由于头波和结尾激波间的扰动传播到地面而形成的压力突跃信号，可引起地表结构、建筑物、人和动物耳膜等受迫振动，对地面环境、社区人群、生物群落产生非预期影响。在协和飞机上天后，美国禁止了陆地超声速民航活动，音爆成为超声速商用飞行的一大障碍。

音爆的形成包括两个关键环节，一是飞行器附近的近场流动扰动，在拓扑结构上表现为头部激波、机翼激波、肩部膨胀波、结尾激波、喷流激波等波系结构的相互作用；二是音爆的远场传播，流动拓扑结构在由飞行器附近向远场传播过程中发生合并、消融与畸变，到达地面形成压力突跃的短时信号——音爆。相应于这两个环节，在音爆预测与研究中，首先通过风洞试验或者流场计算，获得可用压力扰动表征的近场信号；然后，通过传播模型进行音爆传播计算，获得最终的地面音爆信号。目前，音爆传播计算途径主要有两类，一类是经典的、基于改进线性理论的传播计算，另外一类就是基于发展方程的非线性传播计算，二者均采用几何声学理论进行音爆传播的建模处理。其中，基于改进线性理论的音爆传播计算方法主要包括源自nasacr-1299报告中haye等发展的f函数方法和nasatnd-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处按过压比和距离描述的初始信号，需将初始信号转换为时域过压信号并根据需要把信号分段和线性近似，形成n个线性段δpi～τ,i＝1…n构成的折线型波形信号，置当前波形总段数n＝n；波形上每一段δpi(τ)采用该段斜率该段区间长度以及该段与前段之间的突跃来描述；

b.射线管单步跟踪；从当前位置开始，按照给定的射线管步进时间间隔δt，计算射线管推进的新位置，获得相应射线管参数、大气参数及音爆波形参数传播微分方程系数c1,c2，设置音爆信号传播推进的时间步长dt＝δt；

c.射线管单步步进中波形传播的预期步长参数计算；计算音爆信号传播推进的步长参数并将新步长参数置零，即tnew＝0；

d.步长参数t符合性的顺序递进检验及更新、区间长度因子f1,i和f2,i的计算,包括以下步骤：

d1.以步长参数t和当前波形参数，计算波形第一段区间长度表达式的第一因子若f1,1＜0，则按f1,1(tnew)＝0计算新步长参数tnew，更新步长参数t＝tnew；

d2.以步长参数t和当前波形参数，按从前至后顺序递进的方式对波形信号上的各段进行处理，首先计算波形信号上下一段的表达式第一区间长度因子若f1,i+1＜0，则按f1,i+1(tnew)＝0计算新步长参数tnew，更新步长参数t＝tnew；然后计算本段的区间长度表达式第二因子：

若f2,i(t)＜0则按f2,i(tnew)＝0计算[0,t]内的新步长参数tnew，更新步长参数t＝tnew；

d3.波形上所有段步长参数t符合性检验和更新完毕后，如果tnew存在，则新时间步长计算信号波形各段的区间长度因子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′推进后区间长度收缩为零，形成新生激波，其压力跃升为：

e2.如果f2,i＝0，则该段在δt′推进后区间长度收缩零，出现激波融合，其压力跃升为：

e3.如果f1,i＞0，f2,i＞0，则经δt′推进后，该段传播后区间长度大于零，波形正常传播，此时，增加新波形的段数，即inewseg＝inewseg+1，并记录新波形参数为然后将新间断压力跃升变量置为零，即shockdp＝0；

e4.波形所有的信号段传播处理完毕后，将新的波形作为当前波形供后续传播使用，即，,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；即，依据波形参数法中波形各段区间长度表达式的第一因子计算最高斜率段的区间长度第一因子如果f1,k＜0，则获得新的步长参数

d2.对波形上所有的段，按从前至后的顺序，以区间长度表达式的第二个因子检验步长参数t的符合性及其更新，即，以当前步长参数t为输入，计算其区间长度表达式的第二个因子：

如果f2,i(t)＜0，采用数值方法计算[0,t]区间使得f2,i(tnew)＝0的新步长参数tnew，并更新当前步长参数t＝tnew，为后续各段中步长参数t关于区间长度第二因子的符合性检验提供输入；

d3.对确定的步长参数t，按照以及f2,i(t)的表达式计算信号波形上需要更新的各段区间长度因子；如果tnew存在，则计算推进到新激波生成或者激波融合所需时间步长否则按照预期的时间步长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中，曲线为具有三段等斜率的初始信号。

图3b为多波系同步推进物理过程的传统的波形参数法处理过程中各段压力跃升图；

图3b中，+曲线为传统的波形参数法在第8步推进形成的各段压力跃升；□曲线为传统的波形参数法在第9步推进形成的各段压力跃升；×曲线为传统的波形参数法在第10步推形成的各段压力跃升；曲线为传统的波形参数法在第11步推进形成的各段压力跃升。

图3c为多波系同步推进物理过程的本发明的用于音爆传播计算的多波系同步推进波形参数方法处理过程中各段压力跃升图；

图3c中，◇曲线为本发明的用于音爆传播计算的多波系同步推进波形参数方法在第8步推进形成的各段压力跃升；曲线为本发明的用于音爆传播计算的多波系同步推进波形参数方法在第9步推进形成的各段压力跃升。

具体实施方式

下面结合附图和实施例详细说明本发明。

本发明的用于音爆传播计算的多波系同步推进波形参数方法的流程图见图1。

本发明的用于音爆传播计算的多波系同步推进波形参数方法包括以下步骤：

a.初始化；飞行条件、初始信号、终止高度输入，从飞行轨迹开始进行音爆射线管计算，直至初始信号所在高度，置为当前位置，ts＝1，设定射线管步进时间间隔δt；对于给定距纵轴r1处按过压比和距离描述的初始信号，需将初始信号转换为时域过压信号并根据需要把信号分段和线性近似，形成n个线性段δpi～τ,i＝1…n构成的折线型波形信号，置当前波形总段数n＝n；波形上每一段δpi(τ)采用该段斜率该段区间长度以及该段与前段之间的突跃来描述；

b.射线管单步跟踪；从当前位置开始，按照给定的射线管步进时间间隔δt，计算射线管推进的新位置，获得相应射线管参数、大气参数及音爆波形参数传播微分方程系数c1,c2，设置音爆信号传播推进的时间步长dt＝δt；

c.射线管单步步进中波形传播的预期步长参数计算；计算音爆信号传播推进的步长参数并将新步长参数置零，即tnew＝0；

d.步长参数t符合性的顺序递进检验及更新、区间长度因子f1,i和f2,i的计算,包括以下步骤：

d1.以步长参数t和当前波形参数，计算波形第一段的区间长度表达式第一因子若f1,1＜0，则按f1,1(tnew)＝0计算新步长参数tnew，更新步长参数t＝tnew；

d2.以步长参数t和当前波形参数，按从前至后顺序递进的方式对波形信号上的各段进行处理，首先，计算波形信号下一段的区间长度表达式第一因子若f1,i+1＜0，则按f1,i+1(tnew)＝0计算新步长参数tnew，更新步长参数t＝tnew；然后，计算本段的区间长度表达式第二因子：

若f2,i(t)＜0则按f2,i(tnew)＝0计算[0,t]内的新步长参数tnew，更新步长参数t＝tnew；

d3.波形上所有段步长参数t符合性检验和更新完毕后，如果tnew存在，则新时间步长计算信号波形各段的区间长度因子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′推进后区间长度收缩为零，形成新生激波，其压力跃升为：

e2.如果f2,i＝0，则该段在δt′推进后区间长度收缩零，出现激波融合，其压力跃升为：

e3.如果f1,i＞0，f2,i＞0，则经δt′推进后，该段传播后区间长度大于零，波形正常传播，此时，增加新波形的段数，即inewseg＝inewseg+1，并记录新波形参数为然后将新间断压力跃升变量置为零，即shockdp＝0；

e4.波形所有的信号段传播处理完毕后，将新的波形作为当前波形供后续传播使用，即，,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；即，依据波形参数法中波形各段区间长度表达式的第一因子计算最高斜率段的区间长度第一因子如果f1,k＜0，则获得新的步长参数t＝tnew；

d2.对波形上所有的段，按从前至后的顺序，以区间长度表达式的第二个因子检验步长参数t的符合性及其更新，即，以当前步长参数t为输入，计算其区间长度表达式的第二个因子：

如果f2,i(t)＜0，采用数值方法计算[0,t]区间使得f2,i(tnew)＝0的新步长参数tnew，并更新当前步长参数t＝tnew，为后续各段中步长参数t关于区间长度第二因子的符合性检验提供输入；

d3.对确定的步长参数t，按照以及f2,i(t)的表达式计算信号波形上需要更新的各段区间长度因子；如果tnew存在，则计算推进到新激波生成或者激波融合所需时间步长否则按照预期的时间步长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°、马赫数变化率航迹角变化率指向角变化率地面大气压器风洞试验模型长度lmodel＝10ft，风洞试验时音爆测量位置(距速度轴的距离与飞行器长度比值)测量的方位角φ＝47.0°，初始信号见图2中的△曲线。

按照本发明的用于音爆传播计算的多波系同步推进波形参数方法对信号进行参数化。然后，以0.2s的射线管跟踪步长(δt＝0.2)，循环处理射线管推进，循环计算确定步长参数tnew、各信号段的区间因子和多波系传播，最后得到如图2中的○曲线所示的音爆传播到地面的终止结果。

图2中的*曲线给出了传统的波形参数法的外插结果。

由图2可知，对于单步单波系推进情况，本发明的用于音爆传播计算的多波系同步推进波形参数方法的音爆外插结果与传统的波形参数法一致。

实施例2

本实施例的飞行条件与实施例1相同，但初始信号为含有三段等斜率的过压信号，其中，第一、三、十三段斜率相同，见图3a。分别按照thomas波形参数法和本发明的用于音爆传播计算的多波系同步推进波形参数方法进行了图3a初始信号的传播计算。

图3b给出了传统的波形参数法中激波异步形成过程的各段压力跃升结果。可见，采用传统的波形参数法，在传播的第8步推进后各段压力跃升仍然为零，没有激波形成；第9步推进使得第三段压力突跃大于零(□曲线)，表明初始信号第三段收缩为激波，第9步开始产生新激波；第10步推进使得第十二段出现了压力突跃(×曲线)，表明初始信号第十三段收缩为激波；第11步推进又使得第一段出现压力突跃(曲线)，表明初始信号第一段收缩为激波，至此，三道激波形成。从斜率数据上看，第9步使得第三段形成新激波，第一和十三段的斜率从10³增大到了10⁷量级；第10步，使得初始信号第十三段形成新激波，而第一段的斜率则达到了10¹¹量级，第11步使得初始信号第一段形成新激波，新信号中各段斜率均为负值。

图3c为本发明的用于音爆传播计算的多波系同步推进波形参数方法的各段压力跃升结果，◇曲线和曲线分别为本发明的用于音爆传播计算的多波系同步推进波形参数方法在第8、9步推进结果。可见，采用本发明的用于音爆传播计算的多波系同步推进波形参数方法，在第9步，信号的第一、三、十三段同时形成前、中、后三道激波。

由图3可知，本发明的用于音爆传播计算的多波系同步推进波形参数方法复现了多激波同步推进的传播机制，避免了传统的波形参数法异步处理激波同步形成过程中信号斜率急剧增大的问题。