CN109190140B - 基于浸入式边界方法的连续元音生成方法 - Google Patents

Abstract

本发明涉及气体动力学和声学技术，为处理任意的移动复杂边界，避免贴体网格的使用，以及因网格重构和网格畸变带来的数值发散问题，本发明，基于浸入式边界方法的连续元音生成方法，步骤如下：步骤一，初始化包括计算域、网格间距、计算域网格即欧拉网格数、拉格朗日点数、计算步长、声速、空气密度、采样时间、声源点位置的相关变量和运行参数；步骤二，生成欧拉网格和拉格朗日点信息，设置观察点；步骤三，初始化完美匹配层；步骤四，在线性欧拉方程中加入一个附加力来代替复杂边界或者浸入边界，然后再利用时域有限差分方法来离散欧拉方程并迭代计算；步骤五，输出结果。本发明主要应用于连续元音生成场合。

Description

基于浸入式边界方法的连续元音生成方法

技术领域

本发明涉及气体动力学和声学技术领域，具体是涉及一种基于浸入式边界方法的连续元音生成技术。

背景技术

声音是由声源在声道中振动产生的。声波在声道中的传播模拟关键在于欧拉方程的求解。目前有关于声学问题的研究，大多使用传统网格方法，只关注边界固定的静态元音仿真，在对连续发音过程中的动态边界的处理上仍然缺乏有效的数值方法。而浸入式边界方法则为包含移动复杂边界的声传播处理提供了理论依据。时域有限差分方法凭借其高效、准确的特点，在人体发音过程的模拟中得到了广泛应用。在时域有限差分方法的基础上结合浸入式边界方法，在针对移动复杂边界时，避免了因网格重构和网格畸变产生的数值发散问题，实现了人体连续元音生成的数值仿真。

发明内容

为克服现有技术的不足，本发明旨在提出一种基于浸入式边界方法的连续元音生成技术。可以处理任意的移动复杂边界，避免贴体网格的使用，以及因网格重构和网格畸变带来的数值发散问题。为此，本发明采用的技术方案是，基于浸入式边界方法的连续元音生成方法，步骤如下：

步骤一，初始化包括计算域、网格间距、计算域网格即欧拉网格数、拉格朗日点数、计算步长、声速、空气密度、采样时间、声源点位置的相关变量和运行参数；

步骤二，生成欧拉网格和拉格朗日点信息，设置观察点；

步骤三，初始化完美匹配层；

步骤四，在线性欧拉方程中加入一个附加力来代替复杂边界或者浸入边界，然后再利用时域有限差分方法来离散欧拉方程并迭代计算；

步骤五，输出结果：

1)每一个采样周期结束时，保存观察点结果，并输出欧拉网格点计算结果；

2)结束时间循环，并输出最终观察点结果。

步骤一，初始化系统的相关变量和运行参数，具体包括：

计算域为x∈[-0.01,0.2]m,y∈[-0.03,0.03]m,网格间距Δx＝Δy＝0.001m,计算域网格即欧拉网格数大约为20000，拉格朗日点数为329，计算步长Δt＝2×10^-7s,总时间步为800000，模拟时长T＝0.16s。声速为c＝346.3m/s。空气密度为ρ＝1.17kg/m³，采样时间为4×10^-6s。声源点设在(0.0005,0)m处。

步骤三具体地：将麦克斯韦尔Maxwell方程引入到计算区域的人工层，并在人工层中加入阻尼，以此来吸收计算域中的外逃波，人工层为人为的添加在计算域四周的一小块计算空间。

步骤四具体地：

麦克斯韦尔Maxwell旋度方程如下：

在只考虑二维情况下，由Maxwell旋度公式得二维标量方程公式：

其中E,H,ε,σ,μ,σ^*分别表示电场，磁场，介电常数，电导率，磁导率，磁耗率；

将时域有限差分方法应用在声学问题时，声压使用标量P(x,y,t)表示,声速使用矢量u(x,y,t)表示，将上述公式按照以下方式替换：

整理并用矢量形式表示为：

其中P,u,ρ,c,κ,α,α^*分别表示介质压强，介质速度，空气密度，声速，与介质压缩相关的衰减系数，与介质密度相关的衰减系数；

浸入式边界方法是通过添加一个附加力到动量方程(4)来代替边界作用，使用拉格朗日点表示边界，在处理复杂边界时，只需在正交网格上计算控制方程，附加力由离散的几项组成，对于固定边界上的附加力表示为：

f(x)＝f_a(x)+f_i(x)+f_v(x)+f_p(x) 5)

公式的右边分别代表了加速度力项、惯性力项、粘性力项和压力梯度项，具体表现形式为：

其中Re为雷诺数；

浸入式边界方法的实现是在动量方程中添加附加力，该力在拉格朗日点上计算得到，在声学问题中，惯性力项和粘性力项忽略不计，公式(4)改写为：

附加力表示为：

f(x)＝f_a(x_k)+f_p(x_k) 8)

其中,f_a(x_k)和f_p(x_k)分别表示为：

如果正交网格点与最近的拉格朗日点的距离小于先计算该拉格朗日点的附加力，再通过以下公式将附加力插值到相应的正交网格点上，公式如下：

F(i,j)＝(1-ds/h)f(x) 10)

其中i和j为正交网格坐标，ds为正交网格点到拉格朗日点的距离，h为正交网格的对角距离。

将附加力F(i,j)添加到动量方程中，则公式(4)改写为：

时域有限差分和浸入式边界相结合的全离散格式写成如下形式：

由上式，未来某一时刻的压力值取决于相同位置前一时间点的压力值以及与该位置相邻的前后两个速度值，未来某一时间的速度值取决于相同位置前一时间点的速度值以及与该位置相邻的前后两个压力值,依次迭代求解。

本发明的特点及有益效果是：

本发明解决的技术问题是提供一种基于浸入式边界方法的连续元音生成技术，该技术采用时域有限差分方法求解正交坐标系统下的欧拉方程，通过在动量方程中增加一项附加力项来代替流场中的障碍物，可以处理任何的复杂边界和移动边界，避免了使用贴体网格，以及由于网格重构和网格畸变而带来的数值发散问题。

附图说明：

图1基于浸入式边界方法的连续元音生成技术的计算流程图。

图2声道横切面侧视图。

图3元音/a/和/i/的二维声道轮廓。

图4元音/a/到/i/的三个时间点的声波传播模拟结果图。

图5元音/a/到/i/的声波传播的时域图和频率-振幅图。

图6元音/a/和/i/的传递函数(频谱图)。

具体实施方式

本发明解决的技术问题是提供一种基于浸入式边界方法的连续元音生成技术，该技术采用时域有限差分方法求解正交坐标系统下的欧拉方程，通过在动量方程中增加一项附加力来代替流场中障碍物作用，可以处理任意的移动复杂边界，避免贴体网格的使用，以及因网格重构和网格畸变带来的数值发散问题。

为了解决上述问题，本发明的具体技术方案主要包括以下步骤：

步骤一，初始化系统的相关变量和运行参数，具体包括：

实验模拟的计算域为x∈[-0.01,0.2]m,y∈[-0.03,0.03]m,网格间距Δx＝Δy＝0.001m,计算域网格即欧拉网格数大约为20000，拉格朗日点数为329。计算步长Δt＝2×10^-7s,总时间步为800000，模拟时长T＝0.16s。声速为c＝346.3m/s。空气密度为ρ＝1.17kg/m³。采样时间为4×10^-6s。声源点设在(0.0005,0)m处。

步骤二，生成欧拉网格和拉格朗日点信息，设置观察点，具体包括：

欧拉网格信息根据声源公式，初始化声源处压强，计算域中速度信息初始化为0。读入声道面积函数，根据采样点初始化拉格朗日点位置信息，压强和速度信息均初始化为0。在(0.185,0)m处设置观察点。

步骤三，初始化完美匹配层，具体包括：

将Maxwell方程引入到计算区域的人工层，并在人工层中加入阻尼，以此来吸收计算域中的外逃波。人工层为人为的添加在计算域四周的一小块计算空间。具体为在整个计算域边界布置8层人工层作为计算域的完美匹配层，并且设置边界吸收率α＝100。

步骤四，列出求解方程并迭代计算：

时域有限差分方法的基本思想是使用中心差分将安培定理和法拉第定理，即Maxwell旋度方程，在时间和空间上进行离散，分别得到相应的差分方程。中心差分保证了时域有限差分方法的数值解达到二阶精度。Maxwell旋度方程如下：

在只考虑二维情况下，由Maxwell旋度公式可得二维标量方程公式：

将时域有限差分方法应用在声学问题时，声压使用标量P(x,y,t)表示,声速使用矢量u(x,y,t)表示，将上述公式按照以下方式替换：

整理并用矢量形式表示为：

浸入式边界方法主要是通过添加一个附加力到动量方程(4)来代替边界作用。使用拉格朗日点表示边界，在处理复杂边界时，不需要贴体网格，只需在正交网格上计算控制方程。附加力由离散的几项组成，每一项都包含了特定的物理意义。对于固定边界上的附加力可以表示为：

f(x)＝f_a(x)+f_i(x)+f_v(x)+f_p(x) 5)

公式的右边分别代表了加速度力项、惯性力项、粘性力项和压力梯度项。具体表现形式为：

浸入式边界方法的实现是在动量方程中添加附加力，该力在拉格朗日点上计算得到。在声学问题中，惯性力项和粘性力项可以忽略不计，公式(4)改写为：

附加力(拉格朗日力)表示为：

f(x)＝f_a(x_k)+f_p(x_k) 8)

其中,f_a(x_k)和f_p(x_k)分别表示为：

如果正交网格点与最近的拉格朗日点的距离小于先计算该拉格朗日点的附加力，再通过以下公式将附加力插值到相应的正交网格点上，公式如下：

F(i,j)＝(1-ds/h)f(x) 10)

其中i和j为正交网格坐标，ds为正交网格点到拉格朗日点的距离，h为正交网格的对角距离。

将附加力F(i,j)添加到动量方程中，则公式(4)可改写为：

在基于浸入式边界方法的连续元音生成技术中，浸入式边界方法是在线性欧拉方程中加入一个附加力来代替复杂边界或者浸入边界，然后再利用时域有限差分方法来离散欧拉方程。时域有限差分和浸入式边界相结合的全离散格式可写成如下形式：

由上式可知，未来某一时刻的压力值取决于相同位置前一时间点的压力值以及与该位置相邻的前后两个速度值，未来某一时间的速度值取决于相同位置前一时间点的速度值以及与该位置相邻的前后两个压力值,可以依次迭代求解。

在本发明所需解决的问题中使用到的声源公式如下所示：

其中a＝π/t₁，

根据方程(12)计算不同时刻每个网格点的速度和压强信息。

具体计算过程为：

1)对时间变量循环。

2)在每个时间步中，遍历整个计算域的欧拉网格点，如果其与边界的距离小于则计算距其最近拉格朗日点的附加力，并将附加力插值到欧拉网格点上，反之则作为非边界点，忽略声道边界对其的影响。

3)计算并更新网格的速度和压强信息。

4)每间隔一定时间步，在声道中声源位置添加声源波信息。

步骤五，输出结果：

1)每一个采样周期(一个采样周期一般为100个时间步)结束时，保存观察点结果，并输出欧拉网格点上计算出的速度，压强等信息。

2)当结束时间循环后，输出最终观察点保存的整个发音过程中的该点的速度，压强等信息。

进一步地，在上述方案中，步骤一所述初始化系统的相关变量，步骤二所述的生成网格和拉格朗日点信息以及步骤三初始化完美匹配层的具体操作包括：

1)初始化与问题相关的变量信息和运行参数。

2)生成正交网格信息，并添加初始的声源信息。根据声道边界信息，生成拉格朗日坐标系，并初始化拉格朗日坐标信息，设置观察点。

3)在计算域的边界上添加完美匹配层信息，生成完美匹配层。

进一步地，在上述方案中，所述初始化变量信息和运行参数具体设置如下：实验模拟的计算域为x∈[-0.01,0.24]m,y∈[-0.04,0.04]m,网格间距Δx＝Δy＝0.001m,计算域网格即欧拉网格数约为20000，拉格朗日点数为329。计算步长Δt＝2×10^-7s,总时间步为800000，模拟时长T＝0.16s。声速为c＝346.3m/s。空气密度为ρ＝1.17kg/m³。采样时间为4×10^-6s。声源点设在(0.0005,0)m处。

进一步地，在上述方案中，所述初始化欧拉网格、拉格朗日点和观察点信息具体设置如下：根据声源公式初始化声源处压强，计算域中速度信息初始化为0。读入声道面积函数，根据采样点初始化拉格朗日点位置信息，压强和速度信息均初始化为0。在(0.185,0)m处设置观察点。

进一步地，在上述方案中，所述完美匹配层的初始化信息参数如下：在计算域边界布置8层完美匹配层，同时边界吸收率α＝100。

下面结合附图对本发明做进一步详细地描述。

如图3所示，为本发明实验模拟计算域为x∈[-0.01,0.24]m,y∈[-0.04,0.04]m,网格间距Δx＝Δy＝0.001m,计算域网格即欧拉网格数约为20000，拉格朗日点数为329。计算步长Δt＝2×10^-7s,总时间步为800000，模拟时长T＝0.16s。声速为c＝346.3m/s。空气密度为ρ＝1.17kg/m³。采样时间为4×10^-6s。声源点设在(0.0005,0)m处。

基于浸入式边界方法的连续元音生成技术为：

步骤一，初始化系统的相关变量和运行参数，具体包括：

实验模拟的计算域为x∈[-0.01,0.2]m,y∈[-0.03,0.03]m,网格间距Δx＝Δy＝0.001m,计算域网格即欧拉网格数大约为20000，拉格朗日点数为329。计算步长Δt＝2×10^-7s,总时间步为800000，模拟时长T＝0.16s。声速为c＝346.3m/s。空气密度为ρ＝1.17kg/m³。采样时间为4×10^-6s。声源点设在(0.0005,0)m处。

步骤二，生成欧拉网格和拉格朗日点信息，设置观察点，具体包括：

欧拉网格信息根据声源公式，初始化声源处压强，计算域中速度信息初始化为0。读入声道面积函数，根据采样点初始化拉格朗日点位置信息，压强和速度信息均初始化为0。在(0.185,0)m处设置观察点。

步骤三，初始化完美匹配层，具体包括：

将Maxwell方程引入到计算区域的人工层，并在人工层中加入阻尼，以此来吸收计算域中的外逃波。人工层为人为的添加在计算域四周的一小块计算空间。具体为在整个计算域边界布置8层人工层作为计算域的完美匹配层，并且设置边界吸收率α＝100。

步骤四，列出求解方程并迭代计算，具体包括：

时域有限差分方法的基本思想是使用中心差分将安培定理和法拉第定理，即Maxwell旋度方程，在时间和空间上进行离散，分别得到相应的差分方程。中心差分保证了时域有限差分方法的数值解达到二阶精度。Maxwell旋度方程如下：

在只考虑二维情况下，由Maxwell旋度公式可得二维标量方程公式：

将时域有限差分方法应用在声学问题时，声压使用标量P(x,y,t)表示,声速使用矢量u(x,y,t)表示，将上述公式按照以下方式替换：

整理并用矢量形式表示为：

浸入式边界方法主要是通过添加一个附加力到动量方程(4)来代替边界作用。使用拉格朗日点表示边界，在处理复杂边界时，不需要贴体网格，只需在正交网格上计算控制方程。附加力由离散的几项组成，每一项都包含了特定的物理意义。对于固定边界上的附加力可以表示为：

f(x)＝f_a(x)+f_i(x)+f_v(x)+f_p(x) 5)

公式的右边分别代表了加速度力项、惯性力项、粘性力项和压力梯度项。具体表现形式为：

浸入式边界方法的实现是在动量方程中添加附加力，该力在拉格朗日点上计算得到。在声学问题中，惯性力项和粘性力项可以忽略不计，公式(4)改写为：

附加力(拉格朗日力)表示为：

f(x)＝f_a(x_k)+f_p(x_k) 8)

其中,f_a(x_k)和f_p(x_k)分别表示为：

如果正交网格点与最近的拉格朗日点的距离小于先计算该拉格朗日点的附加力，再通过以下公式将附加力插值到相应的正交网格点上，公式如下：

F(i,j)＝(1-ds/h)f(x) 10)

其中i和j为正交网格坐标，ds为正交网格点到拉格朗日点的距离，h为正交网格的对角距离。

将附加力F(i,j)添加到动量方程中，则公式(4)可改写为：

在基于浸入式边界方法的连续元音生成技术中，浸入式边界方法是在线性欧拉方程中加入一个附加力来代替复杂边界或者浸入边界，然后再利用时域有限差分方法来离散欧拉方程。时域有限差分和浸入式边界相结合的全离散格式可写成如下形式：

由上式可知，未来某一时刻的压力值取决于相同位置前一时间点的压力值以及与该位置相邻的前后两个速度值，未来某一时间的速度值取决于相同位置前一时间点的速度值以及与该位置相邻的前后两个压力值,可以依次迭代求解。

在本发明所需解决的问题中使用到的声源公式如下所示：

其中a＝π/t₁，

根据方程(12)计算不同时刻每个网格点的速度和压强信息。

具体计算过程为：

1)对时间变量循环。

2)在每个时间步中，遍历整个计算域的欧拉网格点，如果其与边界的距离小于则计算距其最近拉格朗日点的附加力，并将附加力插值到欧拉网格点上，反之则作为非边界点，忽略声道边界对其的影响。

3)计算并更新网格的速度和压强信息。

4)每间隔一定时间步，在声道中声源位置添加声源波信息。

步骤五，输出结果：

1)每一个采样周期(一个采样周期一般为100个时间步)结束时，保存观察点结果，并输出欧拉网格点上计算出的速度，压强等信息。

2)当结束时间循环后，输出最终观察点保存的整个发音过程中的该点的速度，压强等信息。

尽管上面结合图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (3)

1.一种基于浸入式边界方法的连续元音生成方法，其特征是，步骤如下：

步骤一，初始化包括计算域、网格间距、计算域网格即欧拉网格数、拉格朗日点数、计算步长、声速、空气密度、采样时间、声源点位置的相关变量和运行参数；

步骤二，生成欧拉网格和拉格朗日点信息，设置观察点；

步骤三，初始化匹配层；

步骤四，在线性欧拉方程中加入一个附加力来代替复杂边界或者浸入边界，然后再利用时域有限差分方法来离散欧拉方程并迭代计算；

步骤五，输出结果：

1)每一个采样周期结束时，保存观察点结果，并输出欧拉网格点计算结果；

2)结束时间循环，并输出最终观察点结果，其中，步骤四具体如下：

麦克斯韦尔Maxwell旋度方程如下：

在只考虑二维情况下，由Maxwell旋度公式得二维标量方程公式：

其中E,H,ε,σ,μ,σ^*分别表示电场，磁场，介电常数，电导率，磁导率，磁耗率；

将时域有限差分方法应用在声学问题时，声压使用标量P(x,y,t)表示,声速使用矢量u(x,y,t)表示，将上述公式按照以下方式替换：

整理并用矢量形式表示为：

其中P,u,ρ,c,κ,α,α^*分别表示介质压强，介质速度，空气密度，声速，与介质压缩相关的衰减系数，与介质密度相关的衰减系数；

浸入式边界方法是通过添加一个附加力到动量方程(4)来代替边界作用，使用拉格朗日点表示边界，在处理复杂边界时，只需在正交网格上计算控制方程，附加力由离散的几项组成，对于固定边界上的附加力表示为：

f(x)＝f_a(x)+f_i(x)+f_v(x)+f_p(x) 5)

公式的右边分别代表了加速度力项、惯性力项、粘性力项和压力梯度项，具体表现形式为：

其中Re为雷诺数；

浸入式边界方法的实现是在动量方程中添加附加力，该力在拉格朗日点上计算得到，在声学问题中，惯性力项和粘性力项忽略不计，公式(4)改写为：

附加力表示为：

f(x)＝f_a(x_k)+f_p(x_k)8)

其中,f_a(x_k)和f_p(x_k)分别表示为：

如果正交网格点与最近的拉格朗日点的距离小于先计算该拉格朗日点的附加力，再通过以下公式将附加力插值到相应的正交网格点上，公式如下：

F(i,j)＝(1-ds/h)f(x)10)其中i和j为正交网格坐标，ds为正交网格点到拉格朗日点的距离，h为正交网格的对角距离；

将附加力F(i,j)添加到动量方程中，则公式(4)改写为：

时域有限差分和浸入式边界相结合的全离散格式写成如下形式：

由上式，未来某一时刻的压力值取决于相同位置前一时间点的压力值以及与该位置相邻的前后两个速度值，未来某一时间的速度值取决于相同位置前一时间点的速度值以及与该位置相邻的前后两个压力值,依次迭代求解。

2.如权利要求1所述的基于浸入式边界方法的连续元音生成方法，其特征是，步骤一，初始化系统的相关变量和运行参数，具体包括：计算域为x∈[-0.01,0.2]m,y∈[-0.03,0.03]m，网格间距Δx＝Δy＝0.001m,计算域网格即欧拉网格数大约为20000，拉格朗日点数为329，计算步长Δt＝2×10^-7s,总时间步为800000，模拟时长T＝0.16s，声速为c＝346.3m/s，空气密度为ρ＝1.17kg/m³，采样时间为4×10^-6s，声源点设在(0.0005,0)m处。

3.如权利要求1所述的基于浸入式边界方法的连续元音生成方法，其特征是，步骤三具体地：将麦克斯韦尔Maxwell方程引入到计算区域的人工层，并在人工层中加入阻尼，以此来吸收计算域中的外逃波，人工层为人为的添加在计算域四周的一小块计算空间。