一种具有温度梯度的管道声学模拟方法与流程

1.本发明涉及一种具有温度梯度的管道声学模拟方法。


背景技术：

2.管道是动力机械中传输液体及气体的最基本结构单元，对管道中的声传播进行研究是探究排气消声装置声学计算及分析的核心内容。
3.现有的技术方案通常是研究定常温度场下管道声传播问题，以对定常均匀温度场下的声波动方程进行推导以作为本发明的参照对比。声振动属于宏观的物理问题，因此其传播过程一定符合牛顿第二定律、质量守恒定律及描述状态参数的物态方程。通过以上基本定律，可以得到：传播介质的连续方程，即密度ρ与质点振速u 的关系；介质的动量方程，即声压 p与质点振速u的关系；介质的物态方程，即声压 p 与密度ρ的关系。
4.现需做出一些假设以简化定常温度场声传播问题，这些假设为：1.声传播流体媒介为理想介质，即介质内无粘滞力，声传播理想流体介质传播过程中无能量损耗；2.声传播流体媒介为静态介质，即初始速度为0，且介质为均匀介质；3.声传播为绝热过程，即相邻介质不会因声传播而产生热交换；4.声传播过程中声波为小幅波，声压p的值远小于介质静态压力p0，质点振速u远小于声速c0，质点位移ξ远小于声波的波长λ。
5.一）定常温度场下声波方程的推导相对于环境状态，声扰动通常可以看作是小幅振动。对于流体介质，在没有声扰动时，环境状态可以用压力、速度、密度来表示，这些表示状态的变量满足流体动力学方程。在有声扰动时，状态变量如下式(1)所示：(1)式(1)中：、、分别是压力状态变量、质点振速状态变量、密度状态变量；p’(x,t)、u’(x,t)、ρ’(x,t)分别是声压变化量、质点振速变化量和密度变化量，分别代表对压力、速度和密度场的贡献；、、分别为无声扰动环境状态下的压力、流速与密度。
6.在很多情况下，把流体介质假设为理想化的各向同性静态介质，从而可以实现声学现象的定量分析。在各向同性静态介质中，状态变量、和满足如下式(2)所示的连续方程以及如下式(3)所示的运动方程： (2)
ꢀꢀꢀ
(3)
式(2)、(3)中：、、分别为状态变量、和；为全导数；代表对时间的偏导数；表示梯度。
7.对于各向同性静态介质，将式(1)代入式(2)与式(3)，忽略二阶以上声学小量得到如下式(4)、(5)所示的线性声学方程以及如下式(6)所示的物态方程：
ꢀꢀ
(4)
ꢀꢀ
(5)
ꢀꢀꢀꢀꢀꢀꢀ
(6)式(6)中：比热容c由公式计算，为气体定压比热容c
p
与定容比热容之比cv，即；表示气体压力；表示气体密度。
8.将式(6)代入式(4)消去ρ，然后对时间进行微分，再对式(5)取散度，二者相减得到下式(7)： (7)式(7)中，是 laplace 算子，即梯度的散度。式(7)即为声波方程或波动方程。
9.假设声压随时间的变化关系是简谐的，即声压表示成下式(8)： (8)式(8)中，w表示角频率。将式(8)代入声波方程(7)，得到只含有空间坐标的微分方程，如下式(9)所示： (9)即亥姆霍兹方程，也就是简谐声场的控制方程，k表示波数，。
10.二）定常温度场声波方程的求解在垂直于声波传播的平面上，如果所有位置处的声学量相等，则这类声波称为平面波。此时声学量只是沿着传播方向发生变化。沿声传播方向简谐声场的一维波动方程表示为下式(10)： (10)式(10)中，p(x)表示声压。
11.式(10)的一般解可以取正弦余弦的组合，也可以取复数组合，对于讨论声波在管道内传播的情况，式(10)的解的形式取为下式(11)： (11)式(11)中，a和b为常数。
12.将式(11)代入式(8)得到下式(12)：
ꢀ
(12)将式(12)代入运动方程(5)可得下式(13)：(13)式(13)中，为介质阻抗，。


技术实现要素：

13.本发明的目的是：探究具有温度梯度的管道声传播问题。
14.为了达到上述目的，本发明的技术方案是提供了一种具有温度梯度的管道声学模拟方法，其特征在于，通过建立并对温度梯度的管道声波方程求解以及计算具有温度梯度的管道声传播数值实现管道声学模拟，其中：建立温度梯度的管道声波方程包括以下步骤：步骤101、对于理想的、不考虑粘性、绝热的介质，推导出具有轴向温度梯度的恒定区域管道的一维波动方程，如下式所示：区域管道的一维波动方程，如下式所示：区域管道的一维波动方程，如下式所示：式中：p=p(x,t)、u=u(x,t)、ρ=ρ(x,t)分别表示压力状态变量、质点振速状态变量、密度状态变量，x为管道距离变量，t为时间变量，且有、、，p’(x,t)、u’(x,t)、ρ’(x,t)分别是声压变化量、质点振速变化量和密度变化量，分别代表对压力、速度和密度场的贡献；、、分别为无声扰动环境状态下的压力、流速与密度；r表示摩尔气体常数；t= t(x)表示沿管道距离的温度分布函数；p’=p’(x,t)、u’=u’(x,t)、ρ’=ρ’(x,t)、、、；为气体定压比热容c
p
与定容比热容之比cv，即；步骤102、假设周期解，将一维波动方程化为二阶变系数常微分方程，获得忽略二阶及以上小量的运动方程与能量方程，其中，运动方程如下式所示：能量方程如下式所示：步骤103、将导出的二阶变系数常微分方程从空间x的导数转换到温度t的导数，获得线性温度梯度分布下的声波方程，如下式所示：
式中，ω表示角频率，表示线性的温度分布函数；对温度梯度的管道声波方程进行求解，获得线性温度梯度分布下声波动方程的精确解析解，步骤102中运动方程中的声压用所所求得的精确解析解代入，求得声质点振速u(x)表达式如下式所示：计算具有温度梯度的管道声传播数值具体包括以下步骤：步骤201、建立如下式所示的声学有限元方程：式中，表示梯度，{p}为待求解节点压力的列向量，{n}为形函数，k0表示波数，v表示体积分量，s为直管道声学域分析模型边界表面且有s=sr sv sz，其中，sr、sv和sz分别代表刚性壁面、法向质点振速边界表面、已知法向声阻抗的边界表面；步骤202、将三类边界条件带入步骤201获得的声学有限元方程，获得最终的声学有限元方程：式中，kw表示具有温度梯度的声波波数；[k]为单元刚度矩阵,；[m]为单元质量矩阵,；为单元阻尼矩阵,；[f]为单元载荷向量，；步骤203、对步骤202得到的最终的声学有限元方程进行求解即可得到单元节点声压值，随后，将局部单元坐标转换为整体坐标，并对各单元有限元方程进行组装可得到计算域的声压分布。
[0015]
优选地，对温度梯度的管道声波方程进行求解包括以下步骤：步骤301、将步骤103获得的线性温度梯度分布下的声波方程进行转化并化简得到下式：式中，m表示常数；步骤302、引入独立变量s，有：
步骤303、利用引入的独立变量s，将步骤301获得的公式转换为如下式所示的典型零阶贝塞尔方程：式中，p’表示声压的变化量；步骤304、对步骤303获得的典型零阶贝塞尔方程进行求解，得到声压的精确解析解。
[0016]
优选地，步骤202中，三类边界条件包括刚性壁面边界条件、已知法向质点振速边界条件以及已知法向声阻抗边界条件。
[0017]
优选地，步骤203中，有限元各单元的单元刚度矩阵[k]、单元质量矩阵[m]、单元阻尼矩阵[c]及单元声压载荷向量[f]是在局部坐标系(ξ,η,ζ)下定义的，在数值积分前，建立局部坐标系(ξ,η,ζ)与整体坐标系(x, y, z)的变换关系，再对方程进行数值积分计算，其中：对于单元刚度矩阵[k]以及单元质量矩阵[m]有：元刚度矩阵[k]以及单元质量矩阵[m]有：式中，[j]为jacobi矩阵，且有，其中，对于等参单元，坐标变换和场变量模式用同一个形函数ni，有；对于单元阻尼矩阵[c]及单元声压载荷向量[f]有：对于单元阻尼矩阵[c]及单元声压载荷向量[f]有：式中，；由此得到了整体坐标系(x, y, z)及局部坐标系(ξ,η,ζ)之间的转换关系。
[0018]
本发明在定常温度波动方程的基础上，推导出了考虑线性温度梯度的声波动方程，并引入适用于排气管路的声学边界条件，得到了线性温度梯度下的声传播解析解，结合解析解，分析了温度梯度对声压幅值和波长等声学性能的影响，通过对比不同温度梯度下
的计算结果，得到了声场整体波形会随温降梯度的增加而向声学入口端偏移量的特性。
附图说明
[0019]
图1为直管道声学域分析模型。
具体实施方式
[0020]
本说明书中附图所显示的大小、比例等只是示意性的，用以配合说明书所描述的内容，并非用以限定本发明的实施条件，不影响本发明所产生的功效。本说明书中所述的“上”、“下”、“内”、“外”等位置关系仅是为了方便描述，而非用以限定本发明的可实施范围，其相对关系的改变，在无实质变更技术内容下亦视为本发明的可实施范畴。
[0021]
本实施例公开的一种具有温度梯度的管道声学模拟方法，通过建立并对温度梯度的管道声波方程求解以及计算具有温度梯度的管道声传播数值实现管道声学模拟，具体包括以下内容：一）考虑温度梯度的管道声波方程研究具有轴向温度分布的管道中的声传播问题，所建立的具有温度梯度波动方程及推导解析解过程，包括以下步骤：首先，对于理想的、不考虑粘性、绝热的介质，推导出具有轴向温度梯度的恒定区域管道的一维波动方程。随后假设周期解，将导出的一维波动方程化为二阶变系数常微分方程。由于推导出的二阶变系数常微分方程中存在变量系数，求解温度分布下的精确解析解存在局限。为了得到精确解，将导出的二阶变系数常微分方程从空间x的导数转换到到温度t的导数，最终得到线性温度梯度分布下声波动方程的精确解析解。具体包括以下步骤：1、考虑温度梯度的声波方程的推导假设管道中的介质为无粘性并且无热交换的理想气体，则通过运动方程、能量方程与物态方程对考虑温度梯度的声波动方程进行推导，有下式(14)至式(16)： (14) (15)
ꢀꢀ
(16)式(14)至式(16)中：p=p(x,t)、u=u(x,t)、ρ=ρ(x,t)分别表示压力状态变量、质点振速状态变量、密度状态变量，x为管道距离变量，t为时间变量，且有、、，p’(x,t)、u’(x,t)、ρ’(x,t)分别是声压变化量、质点振速变化量和密度变化量，分别代表对压力、速度和密度场的贡献；、、分别为无声扰动环境状态下的压力、流速与密度；r表示摩尔气体常数；t= t(x)表示沿管道距离的温度分布函数。
[0022]
则有：
ꢀ
(17) (18)式(17)、(18)中：p’=p’(x,t)、u’=u’(x,t)、ρ’=ρ’(x,t)、、、；为气体定压比热容c
p
与定容比热容之比cv，即。
[0023]
将流体介质理想化为各向同性静态介质，假设流体马赫数小于0.1，则。由式(17)与式(18)导出忽略二阶及以上小量的运动方程与能量方程，如下式(19)、(20)所示：
ꢀꢀꢀ
(19)
ꢀꢀꢀ
(20)为了得到声压p’的微分方程，对式(19)求其对空间x的偏导数，对于式(20)求其对时间t的偏导数，则有：
ꢀꢀ
(21)
ꢀꢀ
(22)通过式(21)、式(22)并结合式(19)得到式(23)： (23)由于流体介质为理想化的各向同性静态介质，在式(23)中，，故有： (24)对式(16)进行微分得到： (25)式(25)中，表示线性的温度分布函数。
[0024]
式(25)中，，因此有：
ꢀꢀ
(26)将式(24)与式(26)联立可以得到：
ꢀꢀ
(27)
由状态方程(16)得：
ꢀꢀꢀ
(28)将式(28)代入式(27)化简得：
ꢀꢀꢀ
(29)将式(29)中的声压p’表示为形式，ω表示角频率，则获得考虑温度梯度的声波方程，如下式(30)所示：
ꢀꢀꢀ
(30)2、考虑温度梯度的声波方程的求解由于式(30)所示的声波方程中为关于声压p’的二阶微分方程关于温度分布t的一阶微分方程的多变量表达式，因此无法得到对于确定的温度分布t(x) 的声波动方程的精确解，因此将式(30)中的微分变量进行转化，有下式(31)：
ꢀꢀꢀ
(31)对于线性温度梯度分布的管道，温度随管道距离变化的关系式为：
ꢀꢀꢀ
(32)式(32)中，t0表示初始温度，m表示常数。
[0025]
则将式(32)代入式(31)，并化简得到：
ꢀꢀꢀ
(33)为了进一步简化式(32)，引入独立变量s，s与其他变量间关系如下式(34)所示： (34)即： (35) (36)将式(36)代入式(33)并化简，可转换为下式(37)：(37)式(37)中，p’表示声压变化量。
[0026]
式(37)为典型的零阶贝塞尔方程，其解的形式如下式(38)所示：
ꢀꢀ
(38)式(38)中，c1和c2为不同的常数，j0(s)和y0(s)分别表示不同的关于s的函数。
[0027]
质点振速表示为简谐形式，将解式(37)求得的声压代入运动方程(19)中可求得声质点振速，如下式(39)所示： (39)式(39)中，，则有：
ꢀꢀꢀ
(40)进而求得质点振速u(x)表达式为：
ꢀꢀꢀ
(41)二）具有温度梯度的管道声传播数值计算1、声学有限元方程的建立对于声学有限元数值计算求解过程，待求解的声压被插值函数来表示，即：p={n}
t
{p}
ꢀꢀ
(42)式(42)中，{p}为待求解节点压力的列向量，{n}为形函数。
[0028]
对于静态介质，helmholtz 方程是声传播的控制方程，通过伽辽金加权余量法可以得到下式： (43)式(43)中，k0表示波数，v表示体积分量。
[0029]
通过格林公式可以得到： (44)式(44)中，s为图1所示的直管道声学域分析模型边界表面。如图 1所示边界表面s=sr sv sz，其中，sr、sv和sz分别代表刚性壁面、法向质点振速边界表面、已知法向声阻抗的边界表面。管道声学有限元频域数值求解中，声压p为复数，其实部和虚部与声压p的振幅及相位有关。对于声学有限元频率求解的边界条件可描述为以下三种：a．dirichelt边界条件也称刚性壁面边界条件，即法向声质点振速为0： (45)b．neumann 边界条件，也称已知法向质点振速边界条件，即： (46)
式(46)中，un表示质点速度分量；c．robin 边界条件，也称已知法向声阻抗边界条件，即： (47)式(47)中，z表示声阻抗。
[0030]
确定边界条件类型之后，将三类边界条件 a、b、c 带入式(44)可得到： (48)式(48)中：kw表示具有温度梯度的声波波数；[k]为单元刚度矩阵,
ꢀꢀ
(49)；[m]为单元质量矩阵, (50)；[c]为单元阻尼矩阵,
ꢀꢀꢀꢀ
(51)；[f]为单元载荷向量，
ꢀꢀꢀꢀꢀꢀ
(52)。
[0031]
对式(48)进行求解即可得到单元节点声压值，随后，将局部单元坐标转换为整体坐标，并对各单元有限元方程进行组装可得到计算域的声压分布。
[0032]
2、单元系数矩阵的计算由于有限元各单元的声学刚度矩阵[k]、声学质量矩阵[m]、声学阻尼矩阵[c]及声压载荷向量[f]是在局部坐标系下定义的，在数值积分前需建立局部坐标系与整体坐标系的变换关系。因为单元的形函数是定义在局部坐标系下的，因此首先形函数在整体坐标系下的导数就需通过数学关系转换至其在局部坐标系下的导数，同时对单元矩阵的面积分或体积分也要从整体坐标系转换至局部坐标系，随后再对方程进行数值积分计算。
[0033]
对于整体坐标系(x, y, z)及局部坐标系(ξ,η,ζ)，由复合积分求导法则可使得形函数ni在局部坐标系ξ下的导数为： (53)对于局部坐标η和ζ，采用相同的处理方式可得到形函数ni在局部坐标系η和ζ下的导数。
[0034]
将三个方程联立成为矩阵形式： (54)式(54)中，[j]为jacobi矩阵。
[0035]
通过式(54)可以得到整体坐标系(x, y, z)下形函数ni的导数：
ꢀ
(55)式(55)中，[j]-1
为 jacobi 矩阵的逆矩阵。
[0036]
对于等参单元，根据其定义坐标变换和场变量模式可用同一个形函数ni，可得： (56)将式(56)代入到jacobi 矩阵可以得到：
ꢀꢀ
(57)通过jacobi矩阵的行列式可以由局部坐标系得到整体坐标系的体积微元dv，如下式(58)所示： (58)因此，式(49)所示的单元刚度矩阵和式(50)所示的单元质量矩阵可表示为：
ꢀꢀ
(59)
ꢀꢀꢀ
(60)使用同样的方法，对于面积分的阻尼矩阵和声压载荷向量，需要把整体坐标系下的面积微元转换至局部坐标系下的面积微元。边界的外法向向量可以通过矢量乘积的形式表示，如下式(61)所示：
ꢀꢀ
(61)因此，面积微元ds可表示为：
ꢀꢀꢀ
(62)其中：
ꢀꢀ
(63)因此，式(51)和式(52)可表示为：
ꢀꢀ
(64)
ꢀꢀꢀꢀ
(65)经过上述的转化，得到了整体坐标系(x, y, z)及局部坐标系(ξ,η,ζ)之间的转换关系，从而声学刚度矩阵、声学质量矩阵、声学阻尼矩阵及声压载荷向量在整体坐标系下的积分值可以转换为局部坐标系下并可通过高斯积分求得。