随机声子晶体不确定性分析的直接概率积分法

2023-10-26陈国海杨迪雄

计算力学学报 2023年5期

李骜, 陈国海, 杨迪雄

(大连理工大学工程力学系,工业装备结构分析优化与CAE软件全国重点实验室,大连 116024)

1 引言

当外部激励作用于结构上,会引起结构的振动,并迫使周围的空气也随之振动,从而产生环境噪声。声子晶体PC(phononic crystal)作为新型人工周期性复合结构,因其具备的弹性波带隙特征,在带隙范围内能够很好地阻隔声波和弹性波的传播[1],从而起到减振降噪的作用,引起了广泛关注。

声子晶体是质量密度及弹性常数周期性分布的人工材料或结构[2,3],自提出以来已有近30年的发展,在船舶航海、传感器[4,5]和减振降噪设备[6,7]等领域有了广泛应用。同时研究者发展了多种分析声子晶体能带结构的方法,如平面波展开法、传递矩阵法[8]、有限元方法[9]和集中质量法[10]等。

对于声子晶体的研究,大多学者聚焦于声子晶体确定性分析,然而在实际的声子晶体制备中,由于加工误差和材料杂质等因素,声子晶体中不可避免存在随机性,从而导致其物理响应具有一定的变异性。因此,开展声子晶体的不确定性量化与传播分析,精确预测声子晶体随机响应与可靠度具有重要意义。对于含随机材料参数的声子晶体,使用蒙特卡洛模拟MCS(Monte Carlo simulation)方法虽然结果准确但计算代价昂贵。也有学者采用摄动有限元方法[11]对随机声学超材料进行分析,但受限于摄动法本身的局限性,只能对含小变异系数的随机参数的声子晶体模型进行物理响应的概率密度函数PDF(probability density function)的计算。针对这一问题,文献[12]提出的随机力学分析方法——直接概率积分法DPIM(direct probability integral method)拓展用于含随机材料参数的声子晶体不确定性量化,能够高效准确地预测随机声子晶体的物理响应,并就声子晶体带隙宽度的减振降噪可靠性这一重要问题进行了深入探讨。

本文分别考虑含随机材料参数的Bragg型和局域共振型[3,13]两种机理的声子晶体模型。利用直接概率积分法对两种模型的物理响应(色散曲线、带隙宽度和频率响应)进行不确定性量化,并与蒙特卡洛模拟方法对比,展示直接概率积分法计算随机声子晶体问题的准确性与高效性。而且提出一种用于判断声子晶体减振降噪能力的指标,即随机声子晶体可靠度,用来评估在考虑材料参数的随机不确定性下,弹性波带隙的大小变化对声子晶体减振降噪性能的影响。最后,考察了局域共振型声子晶体三种材料的随机性对其第一带隙可靠度的影响及有关机制。

2 声子晶体的弹性波波动方程

在各向同性的线弹性材料中,弹性波传播的控制方程可表示为

(1)

式中ω为弹性波的角频率,为微分算子,2为拉普拉斯算子,r为空间位置矢量,u(r)为位移矢量,λ和为与材料属性相关的拉梅常数,分别可用材料的弹性模量和泊松比表示为

(2)

式中E为材料的弹性模量,ν为材料的泊松比。

由于声子晶体为周期结构,依据Bloch定理[13],其位移场u(r)为

u(r)=ei(k·r)uk(r)

(3)

式中uk(r)为与声子晶体具有相同周期的函数,k为位于第一布里渊区FBZ(first Brillouin zone)内的波矢量。

基于有限元法求解式(1),可得到声子晶体一个晶胞内的广义特征值方程

[K-ω2M]U=0

(4)

式中K为单胞的刚度矩阵,M为单胞的质量矩阵,U为离散后各节点位移的矩阵。

将Bloch边界条件引入单胞,求解单胞边界上各节点的位移,其满足

U(r+a)=ei(k+a)U(r)

(5)

式中a为声子晶体的晶格常数。

联立式(4,5),原本难以求解的波动方程转化为特征值问题,只需将波矢k依次遍历第一布里渊区内不可约区域的高对称点的连线,即可得到声子晶体的能带结构。图1(a)给出了声子晶体正方形晶格的示意图,正方形晶格的第一布里渊区如图1(b)所示,其中阴影部分为第一布里渊区内的不可约区域。

图1 正方形晶格和相应的第一布里渊区Fig. 1 Square lattice and corresponding FBZ

3 基于直接概率积分法的声子晶体不确定性量化

对于实际应用而言,需要考虑随机性的影响,不确定性参数会使声子晶体的物理响应发生变化,可能导致声子晶体的减振降噪能力减弱。Chen等[12]提出的直接概率积分法能够统一对时变和时不变随机系统的物理响应进行高效和准确的不确定性量化,适用于大变异性随机参数问题。本文采用直接概率积分法对含随机材料参数的声子晶体的物理特性展开研究。下面介绍直接概率积分法刻画随机性传播的基本方程、概率密度积分方程及其数值求解过程,并建立随机声子晶体的可靠度计算公式。

3.1 概率密度积分方程

对于时不变随机系统,可以用映射G来描述输入的随机向量Θ和输出的随机向量Y之间的关系

G∶Y=g(Θ)

(6)

基于概率守恒原理,Chen等[12]推导了式(6)所示的物理映射下,随机不确定性传播服从的概率密度积分方程PDIE(probability density integral equation)

(7)

式中pΘ(θ)和pY(y)分别为输入向量Θ和输出向量Y的联合概率密度函数,ΩΘ为输入随机向量的样本空间,δ(·)为狄拉克函数。概率密度积分方程从随机积分的新视角对随机系统输入和输出之间的不确定性传播过程进行了定量描述。

3.2 随机声子晶体的概率密度积分方程数值求解

在随机声子晶体问题中,狄拉克函数的奇异性对PDIE的解析求解带来不便。为此,Chen等[12]引入了两个关键技术来求解PDIE,(1) 输入概率空间剖分,借助基于GF偏差[15]的两步选点法来生成代表点和获得赋得概率; (2) 狄拉克函数的光滑化处理,采用高斯函数对PDIE的狄拉克函数进行光滑化处理以提高直接概率积分法的计算精度。采用这两项技术,以及式(6)中输入的随机向量Θ(本文为随机材料参数)和输出的随机向量Y(本文为随机声子晶体各物理响应)满足的映射关系,可以确定随机声子晶体各物理响应的概率密度函数

(8)

通过联立式(6)所示的物理映射与式(8)的概率密度积分方程的数值求解,获得随机响应的概率密度函数的方法即为直接概率积分法。借助直接概率积分法,随机声子晶体问题的求解解耦为概率密度积分方程和声子晶体物理响应的计算,进而进行不确定性分析。

3.3 随机声子晶体减振降噪可靠度计算

由于材料的随机性,声子晶体的带隙会发生随机变化。为更好地评价随机材料参数对声子晶体带隙宽度带来的影响,借鉴结构工程分析中可靠度的概念[16],对声子晶体可靠度做如下定义,当随机声子晶体的带隙宽度小于某一值时,相较于原确定性声子晶体而言,声子晶体的减振降噪能力达不到其原本的减振降噪功能,这种未达到基本功能的随机声子晶体是失效的,反之则是可靠的。基于以上论述,本文将随机声子晶体满足减振降噪功能的概率称作其减振降噪可靠性的度量,也即随机声子晶体可靠度。

对于考虑带隙宽度的随机声子晶体问题,带隙宽度的功能函数可以写作如下时不变映射

(i=1,2,…,m)(9)

针对第i个带隙下的减振降噪性能,随机声子晶体的可靠度可表示为

(10)

式中Pr[·]为中括号内随机事件发生的概率,Ωr为系统功能函数Zi的安全域。

根据直接概率积分法,声子晶体系统功能函数的概率密度积分方程可表示为

(11)

将式(11)代入式(10),利用Dirac函数与Heaviside函数的关系为

(12)

声子晶体减振降噪可靠度计算公式为

(13)

式中H[·]为Heaviside函数,pΘ(θ)为输入随机向量Θ的联合概率密度函数,Pq为概率空间剖分后,第q个代表点对应的赋得概率。

4 数值算例

基于直接概率积分法分别对二维Bragg型和局域共振型声子晶体模型的随机响应进行不确定性量化,并通过MCS验证直接概率积分法求解随机声子晶体问题的准确性,本文MCS样本点数目均为105。摄动法是声子晶体随机不确定性分析的一种方法[11],由于摄动法适用于较小的随机参数变异性问题,在大变异性情形下得到的随机响应会有较大的误差。为更贴近于随机变量大变异性的实际需要,本文算例中材料密度的变异系数取0.05(文献[11] 取约0.015),材料弹性模量的变异系数取0.1(文献[11] 取约0.02)。由于摄动法在大变异性问题的精度有限,文献[11]的变异系数都未超过0.02。

4.1 Bragg型声子晶体

考虑如图1(a)所示的单胞为正方形的Bragg型声子晶体,其中蓝色部分的散射体材料为铅,绿色部分为基体材料环氧树脂,单胞的晶格常数a= 0.085 m,散射体的直径d= 0.06 m。其材料参数概率分布列入表1。当材料参数取表1中均值进行确定性分析时,可以得到声子晶体的能带结构,如图2所示,其中,第一带隙区间为[6117.61,12529.98] Hz。对于声子晶体,本文关心的是带隙宽度的大小,但对于随机声子晶体问题,更关心带隙上下边界和带隙宽度的变化。为方便起见,将确定第一带隙下边界的色散曲线的点,即波矢沿Γ方向的第三阶自振频率称作点P1;将确定第一带隙上边界的色散曲线的点,即波矢沿X方向的第四阶自振频率称作点P2。

表1 声子晶体材料参数Tab.1 Material parameters of phononic crystals

图2 Bragg型声子晶体能带结构Fig.2 Energy band structure of Bragg-type PC

图3 点P1和点P2和第一带隙宽度的PDF曲线Fig.3 PDF curves of point P1,P2 and the first bandgap

表2 Bragg型随机声子晶体MCS和DPIM的结果比较Tab.2 Comparison of results for Bragg-type stochastic PCs by MCS and DPIM

基于直接概率积分法,表3和表4分别给出了材料弹性模量和密度的变异系数对Bragg型声子晶体第一带隙宽度的影响。由表3和表4可知,随着材料弹性模量和密度的随机性增大,第一带隙宽度的均值都呈现出减小,而标准差都呈现出增大的趋势,这是由于Bragg型声子晶体的带隙存在与否及其带宽主要依赖于散射体与基体间的物理特性差异(密度和弹性模量)[13]和组分材料的特性等。相较而言,材料弹性模量的随机性对Bragg型声子晶体能带结构的影响更大。

表3 Bragg型声子晶体弹性模量的变异系数对带隙宽度的影响Tab.3 Influence of coefficient of variation for ela-stic modulus of Bragg-type PCs on bandgap width

表4 Bragg型声子晶体密度的变异系数对带隙宽度的影响Tab.4 Influence of coefficient of variation for mass density of Bragg-type PCs on bandgap width

除了能带结构外,频率响应函数FRF(frequency response function)可用于表示波的传播行为[17],也是评价声子晶体减振降噪能力的重要参考。

为获得声子晶体的频率响应函数,考虑如图4所示的二维Bragg型声子晶体有限维模型。沿其底部法线方向向结构内侧施加位移激励,计算结构顶部中心点的法向位移,即可得到有限维声子晶体的频率响应函数。如图5所示,其带隙区间与能带结构的计算结果相近,说明了该有限维模型的合理性。考虑材料的各参数为表1的随机变量时,图6给出了频率为100 Hz时采用直接概率积分法和MCS方法计算得到频率响应的概率密度函数,两种方法的结果高度吻合,CPU耗时列入表2。

图4 有限维Bragg型声子晶体模型Fig.4 Finite-dimensional Bragg-type PCs model

图5 Bragg型声子晶体频率响应函数Fig.5 Frequency response function of Bragg-type PCs

图6 Bragg型随机声子晶体在100 Hz频率响应的概率密度函数Fig.6 PDF of the frequency response of Bragg-type stochastic PCs at 100 Hz

值得指出的是,能带结构的求解是基于声子晶体单胞无限周期性排列的假设,因此即便在考虑材料的随机性时,用于计算色散曲线的随机声子晶体仍为周期性排列。对于有限维声子晶体而言,本文考虑的是整个有限维模型材料参数的随机性,而并非各单胞间材料参数的无序,对这种非严格周期性的无序声子晶体的研究将作为未来的工作。

对Bragg型随机声子晶体进行减振降噪可靠性分析,仅考虑第一带隙,式(9)可写为

FBragg∶ZBragg=Δf1,Bragg(Θ)-5771.13

(14)

式中Δf1,Bragg(Θ)为Bragg型声子晶体不确定性分析时的第一带隙宽度。表4给出了分别基于MCS方法和直接概率积分法计算随机声子晶体的可靠度和时间。由表2可知,对于含不确定材料参数的Bragg型声子晶体模型,直接概率积分法是高效且准确的,表2中直接概率积分法样本点数的选取为自适应选取的结果[18]。

4.2 局域共振型声子晶体

考虑如图7(a)所示三元二维局域共振型声子晶体的单胞,其中散射体由黄色表示的橡胶包围蓝色表示的铝构成,基底由外围绿色代表的环氧树脂表示。单胞晶格常数a=0.02 m,散射体内径din=0.005 m,外径dout=0.008 m。基于表1各参数的均值求解,该局域共振型声子晶体的能带结构如图7(b)所示,带隙区间为[325.38,593.92] Hz。将决定第一带隙宽度的上边界的点称为点P3,第一带隙宽度的下边界由第三阶色散曲线近乎平坦的部分决定,所以将波矢沿X方向的第三阶固有频率称作点P4。

图7 局域共振型声子晶体和能带结构Fig.7 Locally resonant PC and energy band structure

图8 局域共振型声子晶体点P3、点P4和第一带隙宽度的PDF曲线Fig.8 PDF curves of point P3,P4 and the first bandgap of locally resonant PC

基于直接概率积分法,表5和表6分别给出了三种材料弹性模量和密度的变异系数对局域共振型声子晶体的带隙宽度的影响。和4.1节算例中Bragg型声子晶体不同,随着三种材料弹性模量和密度变异系数的增大,第一带隙宽度的均值几乎不发生变化,标准差呈现增大的趋势。这是因为局域共振型声子晶体带隙产生的机理和Bragg型声子晶体不同,其基本单元在特定频率处具有局域共振模式,共振模式将使得声子晶体在相应的本征频率处诱发平直能带或共振带隙,这些局域共振模式主要依赖于单个散射体单胞的自身特性[13],而受材料物理属性的影响较小。

表5 弹性模量变异系数对带隙宽度的影响Tab.5 Influence of coefficient of variation for elas-tic modulus of locally resonant PCs on bandgap width

表6 密度的变异系数对带隙宽度的影响Tab.6 Influence of coefficient of variation for mass density of locally resonant PCs on bandgap width

对于局域共振型声子晶体的频率响应函数,考虑如图9所示的有限维局域共振型声子晶体模型,沿底部法线方向向结构内侧施加位移激励。通过计算结构顶部中心点处的法向位移,可以得到如图10所示的频率响应函数,其带隙区间与能带结构的计算结果相近,说明了该有限维模型的合理性。考虑材料的各参数为表1的随机变量,图11给出了在频率为100 Hz时基于直接概率积分法和MCS方法计算得到的频率响应概率密度函数,两种方法的结果高度吻合,计算的CPU耗时列入表7。

图9 有限维局域共振型声子晶体模型Fig.9 Model of finite-dimensional locally resonant PC

图10 局域共振型声子晶体频率响应函数Fig.10 Frequency response function of locally resonant PC

图11 局域共振型随机声子晶体在100 Hz频率响应的概率密度函数Fig.11 PDF of the frequency response of locally resonant sto-chastic PC at 100 Hz

表7 局域共振型随机声子晶体MCS和DPIM的结果比较Tab.7 Comparison of results of locally resonant stochastic PCs by MCS and DPIM

对于局域共振型随机声子晶体的可靠度分析,仅考虑第一带隙时,式(9)可写为

FLR∶ZLR=Δf1,LR(Θ)-241.68

(15)

式中Δf1,LR(Θ)为局域共振型声子晶体不确定性分析的第一带隙宽度。表7给出了分别基于MCS和直接概率积分法计算随机局域共振型声子晶体各响应可靠度的详细比较结果,由表7可知,对于含不确定性材料参数的局域共振型声子晶体模型,直接概率积分法也能够高效且准确地对其随机响应进行不确定性量化,表7中直接概率积分法样本点数的选取为自适应选取的结果[18]。

由随机Bragg型和随机局域共振型声子晶体的计算结果分析可知,对于随机声子晶体问题,其材料参数的随机性会影响声子晶体的物理响应,并且该影响不可忽略。值得注意的是由于随机性的影响,声子晶体的带隙宽度可能增大也可能减小,换而言之,声子晶体的减振降噪能力可能变强也可能变弱。因此,在声子晶体的设计和制造中适当考虑随机性,或许能够带来声子晶体减振降噪能力的提升。其次,利用直接概率积分法能够精确高效地预测随机声子晶体物理响应的概率密度函数和随机声子晶体的可靠度,这对声子晶体的增材制造以及实际的工程应用具有一定的指导意义。

4.3 材料变异性对局域共振型声子晶体减振降噪可靠性的影响

声子晶体的材料中不可避免地会出现随机不确定性,而哪一种材料的随机性对声子晶体带隙宽度的影响最为显著成为本节关注的重点。本节从随机声子晶体可靠度的观点出发,探讨材料随机性对局域共振型声子晶体可靠度的影响。分别计算局域共振型声子晶体三种常用材料变异系数依次为零(此时为确定性材料)而剩余材料的变异性如表1情况下的声子晶体的可靠度,并将结果与4.2节中含9个随机变量的局域共振型声子晶体可靠度进行比较。

基于直接概率积分法,图12给出了局域共振型声子晶体在橡胶、环氧树脂和铅的材料变异系数分别为零时的可靠度,并与参考组(三种材料的变异性均存在)的可靠度比较。

图12 局域共振型声子晶体不同材料变异系数为零时可靠度的比较Fig.12 Comparison of reliabilities of locally resonant PCs when the coefficient of variance of different materials is zero

可以看出,当环氧树脂和铅的材料参数从随机性变为确定性后,局域共振型声子晶体的可靠度有一定程度的提高,而当橡胶的材料参数从随机性变为确定性后,局域共振型声子晶体的可靠度大幅提高,近似达到1。三种材料对局域共振型声子晶体可靠度的影响程度由大到小依次为橡胶、环氧树脂和铅,其中橡胶材料的随机性对可靠度的影响最为显著。

上述影响机制可以从橡胶的材料特性方面和局域共振型声子晶体的材料布置方式得到解释。橡胶作为一种近乎不可压缩的材料,相较于其余两种材料(表1),其杨氏模量小得多,泊松比却更大,这意味着橡胶的刚度小且质软,其作为散射体的包覆层包覆铅时,在整个声子晶体中最易发生变形,进而引起能量耗散,而局域共振型声子晶体能够具备低频减振的功能也正是由于声子晶体的散射体中具有包覆层橡胶[13],其能够使散射体在局域共振时发生能量耗散,从而达到低频减振降噪的目的。当橡胶材料参数存在随机性时,会对声子晶体的局域共振频率区间造成显著影响,因此当橡胶材料的随机性变小时,局域共振型声子晶体第一带隙对应的可靠度会大幅提高。

表8给出了局域共振型声子晶体橡胶弹性模量和密度的变异系数对第一带隙对应的可靠度的影响。由表8可知,随着材料参数变异性的增高,声子晶体可靠度降低,其中相较于密度而言,弹性模量的变异性对声子晶体可靠度的影响更大。