《1 引言》

1 引言

在声学工程领域中, 计算机模拟预报已成为越来越重要的手段和工具。无论是家用电器, 如空调器、冰箱等的开发设计, 还是国防工业中的武器装备的研制, 都离不开计算机数值分析与模拟。在这些产品的研制过程中, 通过计算机模拟, 可以大大提高产品的性能, 降低成本, 缩短产品研制周期, 产生巨大的经济效益。

子波分析在数值分析计算中的应用是较新的研究领域, 在偏微分方程和边界积分方程的数值解中令人瞩目。这主要是子波有这样几个性质:一是子波具有紧支撑, 由它组成的基函数既有全域基的特点, 又能体现函数的局部特性, 能通过尺度的变化自适应地刨分函数, 或展开系数能自适应地适应函数的局部变化;二是子波的消失矩特性, 它决定了子波变换的数据压缩率和求解精度。理论上, 消失矩越大, 数据压缩率越高。另外, 子波作为基函数有较快的收敛性。即对相同的函数, 用子波展开比用其他基函数展开所须的展开系数少。

数学家最早用子波求解本质光滑无振荡核的积分算子[1], 证实了用子波求解边界积分方程可获得稀疏的系数矩阵。其算法由N维满阵时的O (N3) 或O (N2) 降为O (Nlb N) 。更重要的是:由于系数矩阵的大量稀疏化, 将大大节省计算机内存, 提高边界积分方程求解问题的规模。在应用方面, 子波分析主要在求解电动力学积分方程中受到了极大的关注, 并吸引了大量的研究者[2,3]。这些研究表明, 子波在求解积分方程中能使方程系数矩阵大量稀疏化, 收敛快, 精度良好。这些研究也显示了子波在积分方程的数值求解中的巨大潜力和广阔的应用前景。

在计算声学领域中, 边界积分方程谱方法的研究已令人瞩目。但无论是传统的边界元方法, 还是目前的边界谱方法, 他们的系数矩阵都是非对称满阵。这是边界积分方程方法应用于大型工程问题的瓶颈。因此, 根据边界积分方程的特点, 寻找更合适的插值方法或谱级数, 既能发挥边界积分方程方法的优势, 又能提高他们的收敛速度, 并使他们的系数矩阵成为稀疏矩阵, 已成为众多科学家和工程师们的追求目标。迄今, 就笔者所掌握的资料, 国内外极少见到子波分析用于求解声学中积分方程的报道。现笔者介绍自1996年以来所做的有关子波分析在声辐射和声散射中的应用研究工作。

《2 子波及多分辨分析》

2 子波及多分辨分析

设:尺度空间Vm=Span{φm, n (x) , nZ}, 子波空间Wm=Span{Ψm, n (x) , nZ}。这里φm, n (x) , Ψm, n (x) 分别表示尺度函数和子波函数。对于二进离散子波可表示为:

Ψm,n(x)=2-m2Ψ(2-mx-n)(1)φm,n(x)=2-m2φ(2-mx-n),(2)

式中:m, nZ, Z为整数空间。对于f (x) ∈L2 (R) 的函数, 在尺度空间Vm和子波空间Wm同时展开:

f(x)=Ρm0f(x)+m=-m0Qmf(x),(3)

式中:

Ρmf(x)=nΖ<f,φm,n>φm,n(x),(4)Qmf(x)=nΖ<f,Ψm,n>ψm,n(x)(5)

根据多分辨分析有:

Ρm+1f(x)=Ρmf(x)+Qmf(x)(6)

式 (6) 为正交子波分解的基本结构。Pmf (x) 包含了尺度大于2m的有关f (x) 的全部特性, 是一个模糊的像; Qmf (x) 给出了f (x) 的某些细节。他刻画了分解过程在两个不同尺度2m与2m+1之间的差别, 也即是两个不同分辨率之间的差别。他把函数刻画从一个分辨率提高到下一个更高的分辨率。

《3 用子波形成的边界谱方法》

3 用子波形成的边界谱方法

《3.1积分方程》

3.1积分方程

声学中一般的积分方程形式为:

CΦ(Ρ)+ΩΦ(Q)G(Ρ,Q)dΩ=ΩΦ(Q)G(Ρ,Q)dQ+Φin(7)

式中:C为常数;Ф表示速度势;Φ′表示速度势对边界表面Q点外法线的导数;Фin 为入射声波声压速度势;G (P, Q) 及G′ (P, Q) 分别表示自由空间的格林函数及它对边界表面Q点外法线的导数。其中:

G(Ρ,Q)j4Η0(1)(kr)(8)G(ΡQ)=exp(-jkr)4πr(9)G(Ρ,Q)=02πexp(jkr)4πrcoslθρQdθl=0,1,2(10)

式中, H0(1)为零阶第一类Hankel函数;r为场点P和源点Q之间的距离;k为波数;ρQ为轴对称体Q点的回转半径;l为展开的傅立叶级数的阶。

将式 (7) 中的边界用有限元方法离散, 采用低阶插值函数就可得到传统的边界元方法。

《3.2边界积分方程的子波展开》

3.2边界积分方程的子波展开

在二维及轴对称问题中, 将式 (7) 中的边界变量沿闭曲线Г用子波gk展开[4,5,6,7,8]:

Φ(Γ)=k=02m-1ckgk(Γ)(11)Φ(Γ)=k=02m-1dkgk(Γ)(12)

式 (11) 和式 (12) 中ck, dk为展开系数。如果采用正交子波基, 则上式中的展开系数ck, dk也可以写为:

c=gikΤu(13)d=gikΤq(14)

式 (13) 和式 (14) 中:c, d分别是由ck, dk 组成的列向量;u, q分别为Φ (P) , Φ´ (P) 组成的列向量。将式 (11) 和式 (12) 代入式 (7) 积分方程:

k=02m-1[Cgk(Ρ)-Γgk(Q)G(Ρ,Q)dΓ]ck=k=02m-1[Γgk(Q)G(Ρ,Q)dΓ]dk+Φin(15)

在边界上取2m个不同点, 对每一点应用式 (15) , 可得2m个式 (15) 方程, 并用子波基矩阵gik前乘式 (15) , 则可以得到以子波展开系数为未知量的矩阵形式的代数方程

Ac=Bd+Φi(16)

其中: Φi为2m个Фin点组成的列向量;矩阵AB的元素具体表达式见参考文献[4,5]。如果将式 (13) 和式 (14) 代入式 (16) 则可以得到以边界物理量u, q为未知数的代数方程:

Ηu=Gq+ui(17)

其中, ui=gPkΦi。式 (17) 与通常的边界元有相同的形式, 将边界条件代入即可求解边界未知量。

方程系数中含有奇异积分、非奇异积分, 它们的计算可以分别采用解析积分和高斯数值积分。快速子波变换可用于方程的系数计算, 具体计算参阅文献[4,5,6,9], 这些文献还对方法的计算效率和矩阵的稀疏性做了详细和深入的讨论。

图1是无限长圆柱体的声散射时子波展开方法的系数矩阵实部和虚部的稀疏结构, ka=15。其压缩率分别为14.4 %和9.5 % (矩阵中的非零元素和矩阵元素之比) , 其结果与理论解的范数误差为2.98 %。图2是球面的声散射时, 子波展开方法的系数矩阵实部和虚部的稀疏结构, ka=20。其压缩率分别为10.1 %和6.3 %, 其结果与理论解的范数误差为5.7 %。

《图1》

图1圆柱体声散射系数矩阵的实部、虚部稀疏结构图, 压缩率分别为14.4%和9.5%

图1圆柱体声散射系数矩阵的实部、虚部稀疏结构图, 压缩率分别为14.4%和9.5%  

Fig.1 Sparseness structrues of the real and imaginary part of system matrices for scattering from an infinite cylinder at ka=15, the compression rates are 14.4% and 9.5% respectively in the real and imaginary part

《图2》

图2球面声散射系数矩阵的实部、虚部稀疏结构图, 压缩率分别为10.1%和6.3%

图2球面声散射系数矩阵的实部、虚部稀疏结构图, 压缩率分别为10.1%和6.3%  

Fig.2 Sparseness structures of the real and imaginary part of system matrices for scattering from a sphere at ka=20, the compression rates are 10.1% and 6.3% respectively in the real and imaginary part

《3.3有限元与子波谱方法的耦合》

3.3有限元与子波谱方法的耦合

为了将子波谱方法用于结构声耦合中, 笔者建立了轴对称体声振耦合的子波谱与有限元耦合方法。考虑任意边界条件的轴对称结构有限元形式, 把式 (17) 写为声压形式[4,7]

ΗlΡl=-ρω2GlYl+Ρin(18)

其中ρ表示流体密度, 则声振耦合的矩阵形式为:

[Κl+jωCl-ω2ΜlΤρω2GlΗl][YlΡl]=[FlΡin](19)

声振耦合也可以写成另外两种形式:一种以声学变量为未知量

(Ηl+DΤ)Ρl=DFl+Ρin(20)D=ρω2G1(-ω2Μl+jωCl+Κl)-1,(21)

另一种以结构位移为未知量

(-ω2Μl+jωCl+Κl-jωΤΖl)Yl=Fl-ΤΗl-1Ρin(22)

式 (22) 中Zl是声阻抗矩阵元素。阻抗矩阵又可以分解为Zl=Rl+jXl, 则式 (22) 又可以写成如下形式:

{[Κl+jω(Cl+ΤRl)]-ω2(Μl+ΤXlω)}Yl=Fl-ΤΗl-1Ρin(23)

在式 (23) 中忽略阻尼项和外力项, 即是声学中常用的声与结构相互作用的耦合模态计算式

(Κl-ω2Μl)φil=0(24)

式中:φil为模态向量;Ml为耦合质量矩阵元素。

有限元与子波谱方法的耦合形式与有限元与边界元的耦合形式相似, 可参阅文献[4,7]

《3.4系数矩阵的频率迭代计算》

3.4系数矩阵的频率迭代计算

在声响应计算中, 往往要在一个较宽的频带内求解频率响应。由于在声的系统方程中, 系数矩阵是频率的函数, 因而每计算一个频率点的响应, 都必须重新计算系数矩阵, 所以, 计算效率特别低。文献[10]发表了一个算例:一个球壳声振耦合的散射问题, 用有限元与边界元耦合方法, 在CONVEX CS-1 小型超级计算机上, 以步长0.01, 计算无量纲频率ka从0到100的91个远场点的频率响应, 将需要几十年的运算时间。如采用频率插值技术, 也需100多天的运算时间, 计算效率仍然不是很高。在边界元方法中, 目前除了频率插值技术外, 还没有其他的方法可用。如何提高求解声频率响应的效率仍是一个难题。频率插值技术仅适用于边界元的常元方法, 不能用于其他方法。为了提高频率响应的计算效率, 在边界谱方法中, 若把积分核函数用级数展开, 通过频率迭代方法求解声系统中的系数矩阵, 可大大提高计算效率。

在三维边界积分方程中, 积分核函数的形式为:

E=ejkr/r(25)

式中k为波速。当选取较小的频率间隔Δk, 这时k=k0k。式 (25) 中的指数函数可表示为:

e-jkr=e-j(k0+Δk)r=e-jk0re-jΔkr,(26)

将e-jΔkr展开成泰勒级数, 如取前三项, 我们可以得到一组以下方程:

e-jkrr=e-jk0r(1r-jΔk-Δk2r2)(27)

从上面的级数可以看出, 如果Δk取得较小, 级数的收敛速度很快。分别以e-jkr/r, e-jkr, re-jkr, r2e-j2kr为积分核函数形成的矩阵记为G, G1, G2, G3。以频率间隔Δk, 可导得这些系数矩阵的迭代计算形式。

其计算步骤如下:

a.首先以k0计算第一个频率点的矩阵G0, G01, G02, G03;

G=G0-jΔkG01-Δk22G02,(28)G1=G01-jΔkG02-Δk22G03,(29)G2=2Δk2(G-G0)+j2ΔkG1,(30)G3=2Δk2(G1-G01)+j2ΔkG2,(31)

b.选取Δk, 用式 (28) 和式 (29) 计算G, G1;

c.用式 (30) 、 (31) 计算G2, G3;

d.用求得的G, G1, G2, G3作为新的迭代初始值, 回到b迭代下一个Δk的系统矩阵。

在迭代过程中, 每隔一定的频率间隔做一次修正。例如, 取Δk = 0.01, 每隔100步回到a重新计算初始矩阵。式 (28) 至式 (31) 的级数展开, 可以根据精度增加级数的项数。例如, 4项、6项等, 其迭代方法同上。

现以一个水下钢球壳为算例, 用32个轴对称结构有限元离散结构和32项子波展开声场, 计算了它的耦合频率;以及在球壳上φ=0处作用一点激励时, φ=90°处的水下钢球壳的频率响应, 并与解析解进行了比较, 表明该方法有良好的精度。

频率响应计算采用笔者现提出的频率迭代算法。以Δk=0.01为步长, 一次迭代150次而成, 即从k=1迭代到k=2.5 (238.7~596.8 Hz) 。在奔腾133机上, 完成以上频率响应计算所需总的运算时间不到10 min。如果不采用频率迭代算法, 直接采用耦合方程计算150点的频率响应, 所需的运算时间约为750 min, 频率迭代算法大大提高了计算效率。

《4 工程应用》

4 工程应用

机械噪声是航行器的主要噪声之一, 且产生机理十分复杂。现主要模拟航行器动力系统轴向激励、侧向激励和组合激励, 研究航行器的声辐射特征。航行器主要由4段组成, 用来装载自导设备、战斗部、控制系统、动力燃料、发动机等。各段由轴对称壳体、环形肋骨、圆板等结构组成, 它们之间由楔环连接为一整体。航行器的实际结构十分复杂。笔者用有限元离散航行器结构, 声场用子波展开的边界积分方程建模, 用耦合的任意边界条件的轴对称有限元和子波谱方法, 对航行器进行了水下耦合的模态分析和在多种激励下的声辐射研究, 取得了良好的结果。图3和图4分别是几种组合激励下的声辐射指向性和辐射声功率图;图5是耦合的第一阶弯曲模态。

《图3》

图3k=3时, 沿航行器轴向的声辐射指向性图

图3k=3时, 沿航行器轴向的声辐射指向性图  

Fig.3 Magnitude of the radiated sound pressure on the surface of the submerged object versus polar angles at k=3

《图4》

图4辐射声功率级

图4辐射声功率级  

Fig.4 the radiated sound power level

《图5》

图5耦合的第一阶弯曲模态

图5耦合的第一阶弯曲模态  

Fig.5 The first coupled bending modal shape of the submerged object

《5 结论与展望》

5 结论与展望

提出求解二维及任意边界条件下, 轴对称声辐射和声散射的子波边界谱方法是有效的, 能处理较为复杂的边界几何形状。子波谱方法的系数矩阵是稀疏矩阵, 有较高的压缩率 (正交子波的非零元素的比例为10 %左右) 。在较高的系数矩阵压缩率下, 有较高的精度。由于子波谱方法的系数矩阵的稀疏性, 因而能大大提高边界积分方程方法求解工程问题的能力, 大大提高计算效率, 节约计算成本。因此, 子波谱方法将大大优于传统的边界元方法。

建立的轴对称弹性结构声振耦合的边界子波谱与有限元耦合方法, 能有效地处理任意边界条件下的声振耦合问题。频率响应函数的计算效率一直是人们关心的问题, 目前没有特别高效的算法。通过用级数展开积分方程中的积分核函数, 建立了频率响应函数计算的频率迭代算法, 大大提高了计算效率。这些结果也显示了子波分析在声学工程数值分析中的巨大应用价值。

子波分析及应用是当今世界异常活跃的研究领域, 子波分析的研究已深入到科学与工程的各个领域。随着子波分析研究的深入, 子波对曲面的几何描述与压缩取得了突破性的进展, 因而子波分析在有限元和边界元中的应用会取得重要的突破, 为大规模的科学与工程计算提供高效率的手段, 因此子波分析在这方面的研究前景将是十分光明。

作为一个新的研究方向, 笔者的工作仅仅是开始。尽管在某几个方面做了最基本的工作, 但每一部分均还不完善, 无论是基础理论还是应用方面的研究, 都还有许多艰巨而细致的工作要做。如从数学角度研究方法的数值稳定性, 误差分析, 子波空间的逼近精度等;研究各种子波基在声学积分方程中展开的特点, 如样条子波, 区间子波等, 寻找一种最佳子波基, 使方法具有最佳计算效率和最大压缩率;研究子波自适应数值算法, 以及开展子波谱方法在声学领域的各科学与工程中的应用研究, 并开发相应的计算机软件, 等等。当然须要众多学者在此方面的共同努力, 推动这一研究方向的发展。笔者认为这是一个充满光明, 有巨大应用前景的研究方向。