CN104077479B - 一种基于守恒迎风格式获取参量阵声场空间分布的方法 - Google Patents

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方程进行数值模拟，非线性部分形象、直观的体现了各阶谐波的耦合，但是声场获取时间正比于N²(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原波f₁＝43kHz的声压幅值空间分布特性频域有限差分图；

图3b原波f₁＝43kHz的声压幅值空间分布特性守恒型迎风格式图；

图3c原波f₁＝43kHz的声压幅值空间分布特性轴向声压幅值对比图；

图3d原波f₁＝43kHz的声压幅值空间分布特性径向声压幅值对比图；

图4a原波f₂＝47kHz的声压幅值空间分布特性频域有限差分图；

图4b原波f₂＝47kHz的声压幅值空间分布特性守恒型迎风格式图；

图4c原波f₂＝47kHz的声压幅值空间分布特性轴向声压幅值对比图；

图4d原波f₂＝47kHz的声压幅值空间分布特性径向声压幅值对比图；

图5a差频波f_d＝4kHz的声压幅值空间分布特性频域有限差分图；

图5b差频波f_d＝4kHz的声压幅值空间分布特性守恒型迎风格式；

图5c差频波f_d＝4kHz的声压幅值空间分布特性轴向声压幅值对比图；

图5d差频波f_d＝4kHz的声压幅值空间分布特性径向声压幅值对比图；

图6a两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用频域有限差分计算得到的轴向声压幅值特性图；

图6b两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用频域有限差分计算得到的径向声压幅值特性图；

图6c两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用守恒型迎风格式计算得到的轴向、径向声压幅值特性图；

图6d两列原波和差频波的轴向、一倍瑞利距离处径向声压幅值空间分布特性利用守恒型迎风格式计算得到的径向声压幅值特性图。

具体实施方式

结合附图和实例对本发明进一步说明。

本发明对计算区域进行规则划分，将计算参量阵声场的KZK方程分为线性部分和非线性部分，并利用守恒型迎风格式获取参量阵非线性传播的声场空间分布，可在提高轴向声场步进计算效率的同时，更为准确的计算了参量阵辐射系统声场的分布特点。

本发明解决其技术问题所采用的技术方案包括以下步骤：

(a)读取参量阵声源的几何尺度、轴对称信息、原波频率以及传播媒介的声速、密度和非线性系数等基础数据，根据基础数据获取能匹配参量阵阵元形状的声场计算区域，并对其离散网格化；

(b)读取参量阵发射系统的初始条件，亦即轴向第一层网格各个节点的源点频域信号

(c)通过傅里叶反变换将轴向第一层网格各个节点的源点频域信号变换为源点时域信号p′(r_j,z₁,t_i)；

(d)在获得第一层网格各个节点的源点时域信号时间间隔Δτ＝t_i+1-t_i、轴向步长Δz以及轴向积分步数n＝dz/Δz的基础上，施加周期性边界条件p′(r_j,z₁,t_i)＝p′(r_j，z₁,t_i+2π)，需要注意的是，为保证差分格式的收敛性，要对轴向步长Δz和时间步长Δτ之间的关系进行限制；

(e)利用描述参量阵非线性效应的无粘滞Burgers方程，通过守恒型迎风格式获取轴向第二层网格各个节点的源点时域信号p′_non(r_j,z₂,t_i)；

(f)将第二层网格各个节点的源点时域信号p′_non(r_j,z₂,t_i)通过快速傅里叶变换成源点频域信号p_non(r_j,z₂,k)；

(g)将第二层网格各个节点的源点频域信号p_non(r_j,z₂,k)作为第一层网格的初始条件，利用描述参量阵声场传播衍射、吸收效应的理论模型，获取第二层网格的源点频域信号p_abs,dif,non(r_j,z₂,k)；

(h)重复步骤(c)到(g)，以p_abs,dif,non(r_j,z₂,k)为第二层网格的虚拟激励声源声压，以步长dz进行轴向声场逐步推演，即可获取参量阵声场空间分布。

实施例：

(a)根据参量阵声源以及传播媒介的基础数据确定参量阵声场的有限计算区域，并将其离散化为网格模型；以圆形活塞声源辐射声场为例，化三维计算区域到二维roz平面，采用和边界拟合程度最佳的矩形网格对计算区域分割，建立参量阵声场的有限计算区域模型。沿z方向将0＜z＜z_max区间分成M段，每层对应一个轴向位置z_m，把m＝1称为第1层，每层沿径向即r方向的径向坐标下标j从1变化到J，即每层都要计算J个声压值；

(b)根据参量阵的形成条件，读取z₁层上各个径向节点对应的k阶谐波p(r_j,z₁,k)分量g(r_j,z₁,k)、h(r_j,z₁,k)：

其中：其中：1≤j≤J_，闭区间[1,J_]表示活塞声源所处的径向位置。

(c)将z₁层上k阶谐波p(r_j,z₁,k)的三角函数表达分量g(r_j,z₁,k)、h(r_j,z₁,k)改写为复数形式P(r_j,z₁,k)：

(d)将声场频域信号的变形P(r_j,z₁,k)通过IFFT变换为声场时域信号p′(r_j,z₁,t_i)：

p′(r_j,z₁,t_i)＝K*Re{IFFT[P(r_j,z₁,k)]} (2)

其中：i＝1,2,...,2K

(e)将时域信号p′(r_j,z₁,t_i)限制在一个周期2π内，得到时域波形的时间间隔为Δτ＝t_i+1-t_i＝2π/2K，轴向步长Δz：Δz＝Δτ·cfl/max(p′(r_j,z₁,t_i))及在既定轴向间隔dz区间的时间积分步数满足：n＝dz/Δz；

(f)根据参量阵传播非线性效应理论模型，构造时域信号p′(r_j,z₁,t_i)经轴向间隔轴向步长Δz传播的差分表达式：

其中：

(f.1)当p′(j,z₁,t)≥0时：

当i＝1时，

(f.2)如果p′(r_j,z₁,t_i)＜0：

如果i＝2K，有

(g)综合式(4)(5)，经过的n步计算，得到非线性效应作用下参量阵声场由z₁→z₂层的时域计算结果p′_non(r_j,z₂,t_i)；

(h)将p′_non(r_j,z₂,t_i)进行FFT变换，得到频域解p_non(r_j,z₂,k)：

其谐波分量g(r_j,z₂,k)、h(r_j,z₂,k)满足：

(i)将步骤(h)得到的声压幅值分量p_non(r_j,z₂,k)作为初始条件第一层网格的虚拟源点声压幅值，综合考虑参量阵传播过程中的衍射、吸收效应，获取第二层网格声压幅值p_abs,dif,non(r_j,z₂,k)；

(j)设置合适轴向步长dz，重复(b)至(i)的步骤，利用第z_m+1层上k阶谐波p_abs,dif,non(r_j,z_m+1,k)即可计算出第z_m+2层上的k阶谐波p_abs,dif,non(r_j,z_m+2,k)，依次推演获取整个计算区域的声场分布。

实例一：参量阵声场特性分析

实例参数设置如下：以圆形活塞声源形成的参量阵辐射系统为例，基于活塞声源的轴对称特性，计算有限区域的参量阵声场空间分布(如图2所示)。设换能器的半径为a，径向计算区域为(0,r_max)，其中r_max＝31a，换能器的轴向计算区域为(0,z_max)，其中z_max＝2.5d，d＝πf₀a²/c为换能器辐射原波中心频率对应的瑞利距离。活塞声源的位置[1,a]，径向计算范围为[1,31a]，为了减小边界反射对计算声场的干涉，设定区域[30a,31a]为PML区域，假定轴向每个单位瑞利距离划为120等分，径向每个单位半径长度划分为30等分，可知计算区域的网格数为930×300，两列原波通过一个换能器向轴向辐射声波，取原波频率为f₁＝43kHz，f₂＝47kHz，径向间隔dr＝3.337mm，轴向步进长度dz＝7.799mm，比较常规频域算法和守恒型迎风格式两种算法计算得到的参量阵声场空间分布特性，并将轴向的耗费时间。

图1利用守恒迎风格式获取参量阵声场的流程图

图2参量阵有限计算区域的网格划分模型示意图

图3原波f₁＝43kHz的声压幅值空间分布特性(a)频域有限差分；(b)守恒型迎风格式；(c)轴向声压幅值对比；(d)径向声压幅值对比

图4原波f₂＝47kHz的声压幅值空间分布特性(a)频域有限差分；(b)守恒型迎风格式；(c)轴向声压幅值对比；(d)径向声压幅值对比

图5差频波f_d＝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进行轴向声场逐步推演，获取参量阵声场空间分布。