CN104077479B - 一种基于守恒迎风格式获取参量阵声场空间分布的方法 - Google Patents
一种基于守恒迎风格式获取参量阵声场空间分布的方法 Download PDFInfo
- Publication number
- CN104077479B CN104077479B CN201410298902.8A CN201410298902A CN104077479B CN 104077479 B CN104077479 B CN 104077479B CN 201410298902 A CN201410298902 A CN 201410298902A CN 104077479 B CN104077479 B CN 104077479B
- Authority
- CN
- China
- Prior art keywords
- parametric array
- sound field
- source point
- conservation
- axial
- 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
Landscapes
- Stereophonic System (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明属于参量阵声场领域，具体涉及一种基于守恒迎风格式获取参量阵声场空间分布的方法。本发明包括：读取参量阵声源的几何尺度、轴对称信息、原波频率以及传播媒介的声速、密度和非线性系数，根据基础数据获取能匹配参量阵阵元形状的声场计算区域，进行离散网格化；读取参量阵发射系统的初始条件：轴向第一层网格各节点的源点频域信号；通过傅里叶反变换将轴向第一层网格各节点的源点频域信号变换为源点时域信号；利用描述参量阵非线性效应的无粘滞Burgers方程，通过守恒型迎风格式获取轴向第二层网格各节点的源点时域信号；获取参量阵声场空间分布。本发明利用守恒型迎风格式提高计算效率，对计算区域的网格划分更为规则、简单。
Description
技术领域
本发明属于参量阵声场领域，具体涉及一种基于守恒迎风格式获取参量阵声场空间分布的方法。
背景技术
1962年，Peter Westerwelt首先给出了声学参量阵(parametric array)的最初模型。1965年，H.O.Berkay给出了声学参量阵的更精确和完整的理论解释。不管是Westerwelt还是Berktay理论模型，在它们的推导过程中都做了很多近似，建立的物理模型不够精确。上个世纪70年代，Khokhlov、Zabolotskaya和Kuznetsov三位学者综合考虑原波传播过程中的非线性、吸收以及衍射效应，提供了一个更加精确的声学参量阵模型——KZK方程。该方程避免了利用体积阵模型对虚源求体积积分，是研究参量阵声场空间分布最有力的工具之一。
目前，众多求解KZK方程的数值模拟模型中，最为流行的是时域有限差分(timedomain method)、频域有限差分(frequency domain method)。采用频域方法对KZK方程进行数值模拟，非线性部分形象、直观的体现了各阶谐波的耦合，但是声场获取时间正比于N2(N为最大截断谐波阶次，衍射、吸收项的数值模拟时间正比于N)，故计算量非常大。程琴琴提出一种并行计算模型，对计算区域的进行径向空间分层，利用数个PC机、计算机网络和MPI并行编程分别计算子区域声场并最终求出整个声场。这种算法可以明显提高计算效率，但需要使用多台计算机才能实现并行运算，并且对计算机之间的协同有较高的要求。先永利对程琴琴等的算法加以改进，在一台PC机上设计了一套并行运算软件，采用多线程并行计算以提高计算效率，该算法需要多个处理器并发地执行计算步骤，要求声场计算规模和分割线程数能较好的匹配。本文借鉴时域求解的算子分离思想，将KZK理论模型分解为线性(包含衍射、吸收效应)和非线性(包含非线性效应)两部分，其中线性部分的求解采用后向隐式有限差分(IBFD)和Crank-Nicolsion有限差分(CNFD)相结合的方法，在频域进行；非线性部分则采用守恒型迎风格式积分求解，最后将两部分进行整合，即可得到参量阵的声场空间分布特性。
本发明基于KZK方程时频域结合的数值模拟方法，提出将守恒型迎风格式引入到参量阵非线性声场空间分布获取方法。该方法一方面利用守恒型迎风格式提高计算效率，使得计算量得到较大的减小。另一方面，形象直观地显示出参量阵声场分布情况，更加准确、全面地反映声场性能。
发明内容
本发明的目的在于提供一种提高轴向声场步进计算效率的提高轴向声场步进计算效率。
本发明的目的是这样实现的：
(1)读取参量阵声源的几何尺度、轴对称信息、原波频率以及传播媒介的声速、密度和非线性系数，根据基础数据获取能匹配参量阵阵元形状的声场计算区域，进行离散网格化；
(2)读取参量阵发射系统的初始条件：轴向第一层网格各节点的源点频域信号；
(3)通过傅里叶反变换将轴向第一层网格各节点的源点频域信号变换为源点时域信号；
(4)将源点时域信号限定在一个周期内，施加周期性边界条件，获取第一层网格各节点的源点时域信号时间间隔Δt、轴向步长Δz以及轴积分步数n；
(5)利用描述参量阵非线性效应的无粘滞Burgers方程，通过守恒型迎风格式获取轴向第二层网格各节点的源点时域信号；
(6)将第二层网格各节点的源点时域信号通过快速傅里叶变换成源点频域信号；
(7)将第二层网格各节点的源点频域信号作为第一层网格的初始条件，利用描述参量阵声场传播衍射、吸收效应的理论模型，获取第二层网格的源点频域信号；
(8)重复步骤(3)到(7)，以步长dz进行轴向声场逐步推演，获取参量阵声场空间分布。
本发明的有益效果是：一方面利用时域求解KZK方程的算子分离思想，将KZK方程分为线性部分和非线性部分，并利用守恒型迎风格式提高计算效率。另一方面，对计算区域的网格划分更为规则、简单、可对参量阵声场可视化研究中的能量累积过程、轴向传播声场性能，径向指向性等方面有更为直观、形象的了解，为参量阵的进一步应用提供相应的理论指导。
附图说明
图1利用守恒迎风格式获取参量阵声场的流程图；
图2参量阵有限计算区域的网格划分模型示意图；
图3a原波f1＝43kHz的声压幅值空间分布特性频域有限差分图；
图3b原波f1＝43kHz的声压幅值空间分布特性守恒型迎风格式图；
图3c原波f1＝43kHz的声压幅值空间分布特性轴向声压幅值对比图；
图3d原波f1＝43kHz的声压幅值空间分布特性径向声压幅值对比图；
图4a原波f2＝47kHz的声压幅值空间分布特性频域有限差分图；
图4b原波f2＝47kHz的声压幅值空间分布特性守恒型迎风格式图；
图4c原波f2＝47kHz的声压幅值空间分布特性轴向声压幅值对比图；
图4d原波f2＝47kHz的声压幅值空间分布特性径向声压幅值对比图；
图5a差频波fd＝4kHz的声压幅值空间分布特性频域有限差分图；
图5b差频波fd＝4kHz的声压幅值空间分布特性守恒型迎风格式；
图5c差频波fd＝4kHz的声压幅值空间分布特性轴向声压幅值对比图；
图5d差频波fd＝4kHz的声压幅值空间分布特性径向声压幅值对比图；
图6a两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用频域有限差分计算得到的轴向声压幅值特性图；
图6b两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用频域有限差分计算得到的径向声压幅值特性图；
图6c两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用守恒型迎风格式计算得到的轴向、径向声压幅值特性图；
图6d两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用守恒型迎风格式计算得到的径向声压幅值特性图。
具体实施方式
结合附图和实例对本发明进一步说明。
本发明对计算区域进行规则划分，将计算参量阵声场的KZK方程分为线性部分和非线性部分，并利用守恒型迎风格式获取参量阵非线性传播的声场空间分布，可在提高轴向声场步进计算效率的同时，更为准确的计算了参量阵辐射系统声场的分布特点。
本发明解决其技术问题所采用的技术方案包括以下步骤：
(a)读取参量阵声源的几何尺度、轴对称信息、原波频率以及传播媒介的声速、密度和非线性系数等基础数据，根据基础数据获取能匹配参量阵阵元形状的声场计算区域，并对其离散网格化；
(b)读取参量阵发射系统的初始条件，亦即轴向第一层网格各个节点的源点频域信号
(c)通过傅里叶反变换将轴向第一层网格各个节点的源点频域信号变换为源点时域信号p′(rj,z1,ti)；
(d)在获得第一层网格各个节点的源点时域信号时间间隔Δτ＝ti+1-ti、轴向步长Δz以及轴向积分步数n＝dz/Δz的基础上，施加周期性边界条件p′(rj,z1,ti)＝p′(rj，z1,ti+2π)，需要注意的是，为保证差分格式的收敛性，要对轴向步长Δz和时间步长Δτ之间的关系进行限制；
(e)利用描述参量阵非线性效应的无粘滞Burgers方程，通过守恒型迎风格式获取轴向第二层网格各个节点的源点时域信号p′non(rj,z2,ti)；
(f)将第二层网格各个节点的源点时域信号p′non(rj,z2,ti)通过快速傅里叶变换成源点频域信号pnon(rj,z2,k)；
(g)将第二层网格各个节点的源点频域信号pnon(rj,z2,k)作为第一层网格的初始条件，利用描述参量阵声场传播衍射、吸收效应的理论模型，获取第二层网格的源点频域信号pabs,dif,non(rj,z2,k)；
(h)重复步骤(c)到(g)，以pabs,dif,non(rj,z2,k)为第二层网格的虚拟激励声源声压，以步长dz进行轴向声场逐步推演，即可获取参量阵声场空间分布。
实施例：
(a)根据参量阵声源以及传播媒介的基础数据确定参量阵声场的有限计算区域，并将其离散化为网格模型；以圆形活塞声源辐射声场为例，化三维计算区域到二维roz平面，采用和边界拟合程度最佳的矩形网格对计算区域分割，建立参量阵声场的有限计算区域模型。沿z方向将0＜z＜zmax区间分成M段，每层对应一个轴向位置zm，把m＝1称为第1层，每层沿径向即r方向的径向坐标下标j从1变化到J，即每层都要计算J个声压值；
(b)根据参量阵的形成条件，读取z1层上各个径向节点对应的k阶谐波p(rj,z1,k)分量g(rj,z1,k)、h(rj,z1,k)：
其中：其中：1≤j≤J_，闭区间[1,J_]表示活塞声源所处的径向位置。
(c)将z1层上k阶谐波p(rj,z1,k)的三角函数表达分量g(rj,z1,k)、h(rj,z1,k)改写为复数形式P(rj,z1,k)：
(d)将声场频域信号的变形P(rj,z1,k)通过IFFT变换为声场时域信号p′(rj,z1,ti)：
p′(rj,z1,ti)＝K*Re{IFFT[P(rj,z1,k)]} (2)
其中：i＝1,2,...,2K
(e)将时域信号p′(rj,z1,ti)限制在一个周期2π内，得到时域波形的时间间隔为Δτ＝ti+1-ti＝2π/2K，轴向步长Δz：Δz＝Δτ·cfl/max(p′(rj,z1,ti))及在既定轴向间隔dz区间的时间积分步数满足：n＝dz/Δz；
(f)根据参量阵传播非线性效应理论模型，构造时域信号p′(rj,z1,ti)经轴向间隔轴向步长Δz传播的差分表达式：
其中：
(f.1)当p′(j,z1,t)≥0时：
当i＝1时，
(f.2)如果p′(rj,z1,ti)＜0：
如果i＝2K，有
(g)综合式(4)(5)，经过的n步计算，得到非线性效应作用下参量阵声场由z1→z2层的时域计算结果p′non(rj,z2,ti)；
(h)将p′non(rj,z2,ti)进行FFT变换，得到频域解pnon(rj,z2,k)：
其谐波分量g(rj,z2,k)、h(rj,z2,k)满足：
(i)将步骤(h)得到的声压幅值分量pnon(rj,z2,k)作为初始条件第一层网格的虚拟源点声压幅值，综合考虑参量阵传播过程中的衍射、吸收效应，获取第二层网格声压幅值pabs,dif,non(rj,z2,k)；
(j)设置合适轴向步长dz，重复(b)至(i)的步骤，利用第zm+1层上k阶谐波pabs,dif,non(rj,zm+1,k)即可计算出第zm+2层上的k阶谐波pabs,dif,non(rj,zm+2,k)，依次推演获取整个计算区域的声场分布。
实例一：参量阵声场特性分析
实例参数设置如下：以圆形活塞声源形成的参量阵辐射系统为例，基于活塞声源的轴对称特性，计算有限区域的参量阵声场空间分布(如图2所示)。设换能器的半径为a，径向计算区域为(0,rmax)，其中rmax＝31a，换能器的轴向计算区域为(0,zmax)，其中zmax＝2.5d，d＝πf0a2/c为换能器辐射原波中心频率对应的瑞利距离。活塞声源的位置[1,a]，径向计算范围为[1,31a]，为了减小边界反射对计算声场的干涉，设定区域[30a,31a]为PML区域，假定轴向每个单位瑞利距离划为120等分，径向每个单位半径长度划分为30等分，可知计算区域的网格数为930×300，两列原波通过一个换能器向轴向辐射声波，取原波频率为f1＝43kHz，f2＝47kHz，径向间隔dr＝3.337mm，轴向步进长度dz＝7.799mm，比较常规频域算法和守恒型迎风格式两种算法计算得到的参量阵声场空间分布特性，并将轴向的耗费时间。
图1利用守恒迎风格式获取参量阵声场的流程图
图2参量阵有限计算区域的网格划分模型示意图
图3原波f1＝43kHz的声压幅值空间分布特性(a)频域有限差分；(b)守恒型迎风格式；(c)轴向声压幅值对比；(d)径向声压幅值对比
图4原波f2＝47kHz的声压幅值空间分布特性(a)频域有限差分；(b)守恒型迎风格式；(c)轴向声压幅值对比；(d)径向声压幅值对比
图5差频波fd＝4kHz的声压幅值空间分布特性(a)频域有限差分；(b)守恒型迎风格式；(c)轴向声压幅值对比；(d)径向声压幅值对比
图6两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性(a)(c)分别为利用频域有限差分计算得到的轴向、径向声压幅值特性；(b)(d)分别为利用守恒型迎风格式计算得到的轴向、径向声压幅值特性
分析可知，有限频域差分算法和基于守恒迎风格式的时频域相结合算法求解KZK方程都可以有效获取参量阵原波、差频波的声场空间分布特性，并且二者的分布特性吻合度非常高。
实例二：计算所耗费时间分析
参数与实例一中相同，将轴向声场的延伸长度设定为2.5倍的瑞利距离，分为10等分，每等分约为0.25倍瑞利距离，我们将单位瑞利距离分别做60、120等分，根据声场轴向推演进程到达每个轴向位置节点所耗费时间绘制表格，比较不同声场获取方法的耗费时间。
表1、表2分布表示轴向单位瑞利距离划为60、120等分，两种不同方法获取参量阵声场空间分布在各个轴向位置耗费时间。
表1两种方法获取声场进程所耗费时间对比
表2两种方法获取声场进程所耗费时间对比
经对比分析可以看出：
本专利中给出的基于守恒迎风格式时频域结合获取参量阵声场空间分布的方法，不仅可以有效求解参量阵声场空间分布特性，同时还提高了计算效率。并且随着网格数越多、谐波阶数越大，计算效率越高。
Claims (1)
1.一种基于守恒迎风格式获取参量阵声场空间分布的方法，其特征在于：
(1)读取参量阵声源的几何尺度、轴对称信息、原波频率以及传播媒介的声速、密度和非线性系数，根据基础数据获取能匹配参量阵阵元形状的声场计算区域，进行离散网格化；
(2)读取参量阵发射系统的初始条件：轴向第一层网格各节点的源点频域信号；
(3)通过傅里叶反变换将轴向第一层网格各节点的源点频域信号变换为源点时域信号；
(4)将源点时域信号限定在一个周期内，施加周期性边界条件，获取第一层网格各节点的源点时域信号时间间隔Δt、轴向步长Δz以及轴积分步数n；
(5)利用描述参量阵非线性效应的无粘滞Burgers方程，通过守恒型迎风格式获取轴向第二层网格各节点的源点时域信号；
(6)将第二层网格各节点的源点时域信号通过快速傅里叶变换成源点频域信号；
(7)将第二层网格各节点的源点频域信号作为第一层网格的初始条件，利用描述参量阵声场传播衍射、吸收效应的理论模型，获取第二层网格的源点频域信号；
(8)重复步骤(3)到(7)，以步长Δz进行轴向声场逐步推演，获取参量阵声场空间分布。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410298902.8A CN104077479B (zh) | 2014-06-26 | 2014-06-26 | 一种基于守恒迎风格式获取参量阵声场空间分布的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410298902.8A CN104077479B (zh) | 2014-06-26 | 2014-06-26 | 一种基于守恒迎风格式获取参量阵声场空间分布的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104077479A CN104077479A (zh) | 2014-10-01 |
| CN104077479B true CN104077479B (zh) | 2017-02-22 |
Family
ID=51598731
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410298902.8A Active CN104077479B (zh) | 2014-06-26 | 2014-06-26 | 一种基于守恒迎风格式获取参量阵声场空间分布的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104077479B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109918821B (zh) * | 2019-03-15 | 2020-01-14 | 中国水利水电科学研究院 | 一种迎风守恒型河道漫溢出流数值模拟方法 |
| CN109933949B (zh) * | 2019-04-02 | 2022-08-02 | 哈尔滨工程大学 | 一种建立含气泡水介质中波动-振动非线性声场的方法 |
| CN112198515B (zh) | 2020-10-13 | 2021-06-29 | 湖南国天电子科技有限公司 | 一种参量阵浅剖差频转换性能优化方法 |
| CN115665633B (zh) * | 2022-12-26 | 2023-03-31 | 中国人民解放军海军工程大学 | 一种参量阵扬声器基波调制的方法、记录媒体及系统 |
| CN117932199B (zh) * | 2024-03-25 | 2024-05-28 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种时域声源定位方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101007230A (zh) * | 2006-12-22 | 2007-08-01 | 西安建筑科技大学 | 袋式除尘器脉冲喷吹清灰性能的数字实验方法 |
| CN101951355A (zh) * | 2010-09-16 | 2011-01-19 | 华南理工大学 | 一种降低ofdm前导码峰均比的方法 |
| CN101978627A (zh) * | 2008-01-25 | 2011-02-16 | 株式会社Ntt都科摩 | 移动通信系统中的发送装置及方法 |
| CN102289929A (zh) * | 2011-06-02 | 2011-12-21 | 西北工业大学 | 带耗散项lwr宏观交通流稳定性建模方法 |
| CN103143126A (zh) * | 2013-04-03 | 2013-06-12 | 南京大学 | 生物组织非线性hifu声场确定的方法 |
| CN103631992A (zh) * | 2013-11-07 | 2014-03-12 | 华南理工大学 | 一种自吸泵自吸过程流动模拟的计算方法 |
-
2014
- 2014-06-26 CN CN201410298902.8A patent/CN104077479B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101007230A (zh) * | 2006-12-22 | 2007-08-01 | 西安建筑科技大学 | 袋式除尘器脉冲喷吹清灰性能的数字实验方法 |
| CN101978627A (zh) * | 2008-01-25 | 2011-02-16 | 株式会社Ntt都科摩 | 移动通信系统中的发送装置及方法 |
| CN101951355A (zh) * | 2010-09-16 | 2011-01-19 | 华南理工大学 | 一种降低ofdm前导码峰均比的方法 |
| CN102289929A (zh) * | 2011-06-02 | 2011-12-21 | 西北工业大学 | 带耗散项lwr宏观交通流稳定性建模方法 |
| CN103143126A (zh) * | 2013-04-03 | 2013-06-12 | 南京大学 | 生物组织非线性hifu声场确定的方法 |
| CN103631992A (zh) * | 2013-11-07 | 2014-03-12 | 华南理工大学 | 一种自吸泵自吸过程流动模拟的计算方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104077479A (zh) | 2014-10-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104077479B (zh) | 一种基于守恒迎风格式获取参量阵声场空间分布的方法 | |
| Li et al. | Wavelet-based numerical analysis: A review and classification | |
| Lockard | A comparison of Ffowcs Williams-Hawkings solvers for airframe noise applications | |
| Song | The scaled boundary finite element method in structural dynamics | |
| Amiri et al. | XLME interpolants, a seamless bridge between XFEM and enriched meshless methods | |
| CN109883532B (zh) | 一种声源识别与声场预报方法 | |
| CN103970717B (zh) | 基于Associated Hermite 正交函数的无条件稳定FDTD方法 | |
| CN106354695A (zh) | 一种仅输出线性时变结构模态参数辨识方法 | |
| CN103617344A (zh) | 基于雷达后向散射实测数据对单层地表介电参数与粗糙度参数快速反演的联合优化算法 | |
| Hornikx et al. | openPSTD: The open source pseudospectral time-domain method for acoustic propagation | |
| CN106168942B (zh) | 一种基于奇异边界法的波动类型动态数据重构方法 | |
| CN104268322B (zh) | Weno差分方法的一种边界处理技术 | |
| Li et al. | Reduced-order thrust modeling for an efficiently flapping airfoil using system identification method | |
| CN107526105A (zh) | 一种波场模拟交错网格有限差分方法 | |
| CN104408295A (zh) | 一种大跨桥梁下部结构风-浪耦合作用荷载数值模拟方法 | |
| CN113158527A (zh) | 一种基于隐式fvfd计算频域电磁场的方法 | |
| CN104573376A (zh) | 一种时域有限差分法计算电磁散射的瞬态场远场外推方法 | |
| Bae et al. | On the use of Brinkman penalization method for computation of acoustic scattering from complex boundaries | |
| CN104597488A (zh) | 非等边长网格波动方程有限差分模板优化设计方法 | |
| CN104142908A (zh) | 电波传播抛物方程分步傅里叶变换解的上边界处理方法 | |
| CN105260615A (zh) | 一种粮食消费量预测方法 | |
| CN104949628A (zh) | 基于二维正交曲率的柔性板状结构复杂形态重建方法 | |
| CN109164488A (zh) | 一种梯形网格有限差分地震波场模拟方法 | |
| CN109657288A (zh) | 一种三维显隐时域电磁学数值方法 | |
| CN104215964A (zh) | 一种多列等差频率原波相互作用形成参量阵的声场获取方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant |