基于小波有限元模型的二维声子晶体板结构带隙设计方法与流程

本发明属声学功能材料结构设计领域，具体涉及一种基于小波有限元模型的二维声子晶体板结构带隙设计方法。

背景技术：

近年来，以自然晶体中电子能带理论为基础，学者对周期结构中弹性波的传播产生了浓厚的兴趣，势必找到振动控制的良策。1993年，Kushwaha等人在研究材料周期变化的结构时第一次使用了声子晶体的概念。并指出声子晶体具有的带隙特性能够应用到高精密、无振动环境中。1995年，R.Martinez-Sala等人对雕塑“流动的旋律”做声学试验测试时，首次实验验证了存在弹性波带隙。从这之后，声子晶体得到了国内外众多学者的关注。

平面波展开法(PWE)是声子晶体研究中最常用的算法之一，适用于各维声子晶体板结构带隙计算。其基本理论是：因为声子晶体结构是材料周期性的，可以把材料参数等按傅里叶级数展开，并利用Bloch定理，在倒格矢空间中，可以将波的动力学方程以平面波叠加的形式展开，此时动力学方程就转化为特征值问题，经过计算特征值得到能带结构。PWE法在计算固/固、液(气)/液(气)等形成的多种声子晶体时十分有效，但是对于计算液(气)/固体组成的声子晶体时面临很大困难。当材料构成差异很大时，计算收敛很慢，耗费大量时间，而且其结果也不精确。

传递矩阵法(TM)是目前计算一维声子晶体带隙特性使用较多的方法。此方法首先推导出单个周期的传递矩阵。然后通过施加周期边界条件，从而得出频散曲线的解析解；在限定的周期结构振动传输特性计算中，只要利用传递矩阵乘积就能获得其结果。虽然TM法计算量很小，但目前较难分析二维和三维问题。

多重散射法(MST)可用于二维和三维声子晶体的带隙特性计算。该方法在求解声子晶体带结构和弹性波反射/透射系数时，是先通过分析散射弹性波与各个散射体的入射关系而得出的。但该方法存在一 定的局限，主要解决二维圆柱和三维球状散射体构成的声子晶体。

有限元法(FEM)已广泛用于复杂工程结构对象的定量分析。尤其在固体力学和结构分析领域，许多通用程序直接用于工程应用。求解的基本思想是将连续的求解区域进行离散处理，然后根据变分原理及弹性力学基本方程，将未知场函数用单元内假设的近似函数来分片的表示，再结合平衡条件得出有限元求解方程，引入边界条件，就可以通过插值函数计算出整个求解域的近似值。除了上述常用的算法外，国内外学者还提出了时域有限差分法、边界元法、谱有限元法、变分法、集中质量法等，在这里不作详述。

虽然上述有些算法已经得到了大量应用，但共性缺点在于精度不高、收敛慢等问题，这制约了声子晶体板结构带隙设计，应用于工程实践。

技术实现要素：

为了克服以上的技术不足，本发明提供一种基于小波有限元模型的二维声子晶体板结构带隙设计方法。

本发明提供一种基于小波有限元模型的二维声子晶体板结构带隙设计方法，其包括以下步骤：

一、通过将区间B样条小波与有限元法相结合，建立二维声子晶体板结构带隙特性计算模型；

二、采用一中所构造的计算模型，并在频域内，结合单胞结构和周期边界条件PBCs获得声子晶体的带隙特性；

三、为获得满足特定频带要求的带隙特性，通过不断计算调整二维声子晶体板结构尺寸，在固定晶格常数的前提下，通过获得最佳填充率，确定散射体几何尺寸关系，最终完成二维声子晶体板结构带隙设计，获得二维声子晶体板结构尺寸。

一中包括以下步骤：

1)获得板势能泛函

通过变分原理，令δΠ_p＝0，可得到单元求解方程

其中，单元载荷列阵为

获得单元刚度矩阵为

并得到一致质量矩阵

并获得声子晶体板单胞的特征频率方程：

(K-ω²M)v＝Dv＝0

其中D＝(K-ω²M)表示动力刚度矩阵，ω为角频率，K和M分别表示总体刚度矩阵与质量矩阵，v为总体自由度排列。

二中包括以下步骤：

首先将单胞节点分成9组，即4个边界节点、4个角节点和内部节点，则将带隙计算模型转换为

其中，v₄为内部节点自由度，v₃,v₇,v₈和v₉表示4个角点自由度，v₁,v₂,v₅和v₆表示4个边界节点自由度；

其次将波矢k作为横坐标x、特征频率为纵坐标y，当具体波矢k＝[k_x k_y]，k_x k_y为波数，当其在第一Brillouin边界取值时，就可得到声子晶体板结构带隙特性，可用简约波矢M、Γ、X做为横坐标，用频率作为纵坐标描述。

本发明的有益效果：本发明由于将BSWI与有限元法相结合，用BSWI尺度函数取代传统有限元的多项式插值，进而建立二维声子晶体板结构带隙特性计算模型。结合单胞结构和周期边界条件PBCs就能计算声子晶体的带隙特性。具有下列区别于传统有限元求解方法的显著优势：

1)BSWI结合传统有限元的多功能性与B样条函数的优良逼近性进行结构分析，在声子晶体带隙计算过程中，用精确的BSWI尺度函数代替传统的多项式插值来形成形状函数，进而构成单元，因此，可以方便计算出刚度与质量矩阵，能够利用较少的单元与自由度数获得较高的计算精度。这使得BSWI小波有限元较传统有限元法具有更好的计算效率和收敛性；

2)本发明结合ω(k)技术将复杂的波矢问题转化为纯粹的实频域内来解决，并通过动力凝聚将频散曲线简化为三次多项式特征值问题；

3)通过构建二维声子晶体板结构带隙设计的小波有限元模型，不断调整声子晶体板结构尺寸，可高精度、快速收敛地获得所需要的带隙特性，最终完成声子晶体板结构带隙设计。

附图说明

图1是本发明的板单元求解域Ω_e。

图2a是本发明的单胞的几何布局图

图2b是本发明的单胞的离散布局图。

图3是本发明的小波有限元模型计算的带隙特性(实线)和传统有限元模型计算的带隙特性(点与虚线)。

图4是本发明的长方晶格的带隙特性。

图5是本发明的填充率为f＝32.6％时长方晶格的带隙特性。

图6是本发明的各组元材料的弹性常数表。

具体实施方式

下面结合附图对本发明实施例作进一步说明：

如图所示，本发明提供如下方法：

1)将区BSWI与有限元法相结合，用BSWI尺度函数取代传统有限元的多项式插值，进而建立二维声子晶体板结构带隙特性计算模型。

板势能泛函为

式中，Ω_e为单元求解域，t为单元厚度，f＝{f_x f_y}^T为体力向量，面力向量p＝{p_xp_y}^T，位移场向量u＝{u v}^T，u_i＝{u_i v_i}^T为集中载荷作用点的位移，弹性阵D为

式中，E为弹性模量，μ为泊松比，应变阵ε＝{ε_x ε_y γ_xy}^T为

应力应变关系为

σ＝{σ_x σ_y τ_xy}^T＝Dε (4)

采用二维BSWI尺度函数插值时，标准区间Ω_s中的未知位移场函数表示为

式中，为某尺度下BSWI的尺度函数，u^e,v^e为物理空间自由度列向量，表示为

式中，i＝1,2,…,n+1。

单元矩形求解域Ω_e和单元节点及自由度排列如图1所示，单元 边长分别为l_ex和l_ey。

对势能泛函式(1)，首先将单元求解域Ω_e映射到单元标准求解域Ω_s，然后将式(2)、(3)和式(5)代入式(1)，并由变分原理，令δΠ_p＝0，可得到单元求解方程

其中，单元载荷列阵为

单元刚度矩阵为

K^e,3＝(K^e,2)^T；

以上各式中积分项为

将i,j＝0,1中l_ex，dξ和分别用l_ey，dη和替换，可得到 i,j＝0,1表达式。

由于u和v各自独立插值，可以将单元求解方程式(7)按照单元节点自由度{u_1,1v_1,1…u_1,n+1 v_1,n+1|…|u_n+1,1 v_n+1,1…u_n+1,n+1 v_n+1,n+1}^T进行排列，从而得到BSWI板单元有限元求解方程

式中

式中，各元素分别为式(7)中各子阵K^e,1、K^e,2、K^e,3和K^e,4中对应元素一致质量矩阵为

如图所示，板的单胞结构图，一个声子晶体板单胞由矩形B部分镶嵌于矩形A中，B部分的几何尺寸为长L_a宽L_b，矩形A部分几何尺寸为长L_x宽L_y。那么，声子晶体板结构就由单胞结构向XY平面延伸得到的。

如图2(b)，现将单胞结构离散成9个BSWI4_j(4表示B样条小波阶数,j为B样条小波尺度)有限元单元，分别编号为1,2,…,9。单胞内相应的自由度v也被分成9组，即v_i(i＝1,…,9)，其中，v₄为内部节点自由度，v₃,v₇,v₈和v₉表示4个角点自由度，v₁,v₂,v₅和v₆表示4个边界节点自由度。

二维声子晶体板单胞的特征频率方程为

(K-ω²M)v＝Dv＝0 (13)

这里，D＝(K-ω²M)表示动力刚度矩阵，ω为角频率，K和M分别表示总体刚度矩阵与质量矩阵，v为总体自由度排列。

总体刚度矩阵K和质量矩阵M是由9个BSWI4_j单元的刚度矩阵与质量矩阵叠加而得到的。和表示A部分的单元刚度矩阵与单元质量矩阵，和分别表示B部分的单元刚度矩阵与质量矩阵。K^A,j和M^A,j分别为A部分的局部总体刚度矩阵与质量矩阵；K^B,j和M^B,j分别表示B部分的局部总体刚度矩阵与质量矩阵。

为了简化计算，将单胞的节点分成9个子块，则式(13)变为

在边界节点x＝L_x、y＝L_y、x＝0和y＝0处施加周期边界条件，则矩阵将被分成以下4个部分，它们相应的表达式为

这里，k_x和k_y表示波数。边界节点自由度自然就与波数建立起了联系。内部节点自由度v₄将通过以下动力凝聚矩阵得到缩减

其中，

式中，n₁＝2(2^jm+2m+2^jn+2n-1)，n₂＝2(10m+10n+1)，n_i为单胞的内部节点自由度，m和n分别表示单胞在X和Y方向上的单元数，I和0分别表示相应维度上的单位矩阵与零矩阵。最终，二维声子晶体板的带隙计算方程为

其中，

其中，m₁＝2(2^jm+3m-1)，m₂＝2(2^jn+3n-1)，m₃＝2。方程呈现了一个关于二维声子晶体板离散结构的实对称特征值问题，可以用一个实矩阵表达式来替换复杂的矩阵方程。通过求解特征值方程，每一个波矢k可以得到与之对应的一组特征频率，每个特征频率有其对应的特征向量，特征向量表示在该频率下声子晶体结构的运动形式。将波矢k作为横坐标x、特征频率为纵坐标y，当具体波矢k＝[k_x k_y]，k_x k_y为波数，当其在第一Brillouin边界取值时，就可得到声子晶体板结构带隙特性，可用简约波矢M、Γ、X做为横坐标，用频率作为纵坐标描述。

实施例1：本实施例主要验证二维声子晶体板结构带隙计算的小波有限元数值求解模型计算精度。为了验证所构造的BSWI板单元在声子晶体带隙计算方面的正确性与有效性，本节将给出算例验证。如图2所示，散射体采用矩形铅板(B部分)周期排列镶嵌于环氧树脂(A部分)基体中构成了二维声子晶体板，其填充率为f＝(L_a×L_b)/(L_x×L_y)＝11％，限定晶格常数L_x与L_y均等于0.03m，铅与环氧树脂的材料参数见图6。

首先采用9(m＝3,n＝3)个BSWI4₃板单元去计算二维声子晶体带隙特性，这里使用的BSWI尺度函数为4阶3尺度，9个BSWI单元总体自由度为1992个，其计算结果见图3实线所示。同样，也采用了50×50，80×80个传统板单元去计算二维声子晶体带隙，其总体自由度分别为5202个与13122个，计算结果见图3虚线与点线表示。从图3通过对比得到，在低频带隙区域，9个BSWI4₃板单元与传统板单元吻合的很好，但是在高频带隙区域，小波有限元法较传统有限元 法具有更好地收敛性与稳定性。有限元法已经成功计算了二维声子晶体带隙且计算精度较好。图3中显示9个BSWI4₃板单元已经比900个传统板单元的精度还高，这表明了本节所构造的BSWI板单元在二维声子晶体带隙计算方面具有了很高的精度，能够以很少的计算自由度获取很高的计算精度，这势必减少了计算时间和提高了计算效率。

实施例2：本实施例主要验证二维长方板声子晶体带隙计算的小波有限元数值求解模型计算效率。长方板作为结构部件也广泛应用于工程实际，本文将利用BSWI小波有限元法研究不同晶格形式的声子晶体带隙。铅与环氧树脂的材料参数见图6。

我们采用长方晶格去研究声子晶体带隙的影响，其中L_x不等于L_y。相应的结构参数为：L_x＝0.03m、L_y＝0.02m、L_a＝L_b＝0.01m、填充率为f＝16.7％。分别采用9个BSWI4₃板单元(m＝3,n＝3)、20×20、40×40个传统长方板单元去计算声子晶体带隙，其结果见图4，分别用点、实线、正方框表示其结果。与正方晶格声子晶体具有相似的计算结果，在低频域内，BSWI小波有限元法的带隙计算结果与传统有限元法基本吻合。但是，在高频域内，9个BSWI4₃板单元1992个自由度明显小于传统板单元的7442个，且计算精度高。这说明了BSWI小波有限元在计算长方晶格时也具有较高的精度和较快的收敛速度。这大大降低了计算时间，提高了计算效率。

实施例3：本实施例主要给出利用基于小波有限元模型计算得到的具有较宽带隙的二维长方板声子晶体结构尺度范围。铅与环氧树脂的材料参数见图6。不失一般性，固定结构参数为：L_x＝0.03m、L_y＝0.02m。

采用9个BSWI4₃板单元(m＝3,n＝3)计算声子晶体带隙。为获得满足特定频带要求的带隙特性，通过不断计算调整二维声子晶体板结构尺寸，在固定晶格常数的前提下，通过获得最佳填充率，确定散射体几何尺寸关系，最终完成二维声子晶体板结构带隙设计，获得二维声子晶体板结构尺寸。当填充率为f＝32.6％时能获得最低的第一完全带隙起始频率，由f＝(L_a×L_b)/(0.03×0.02)＝32.6％，可得L_a×L_b＝1.956×10^-4m²。满足L_a×L_b＝1.956×10^-4m²要求的散射体(B部分)尺寸的带隙求解结果见图5，第一完全带隙如灰色区域所示，其起始 频率(灰色区域的下边界)低至13.1KHz，而截止频率(灰色区域的上边界)高至21.6KHz。禁带(该频率范围内声波不能通过)宽度为8.5KHz，而工程应用中高频噪声范围是：10000Hz～20000Hz。因此，该二维长方板声子晶体结构具有良好的抑制工程中高频噪声的能力。

以上数值算例表明，与传统有限元相比，BSWI平面板单元在计算二维声子晶体带隙方面具有计算精度高和收敛性快的特点。最后，可通过构建二维声子晶体板结构带隙设计的小波有限元模型，不断调整声子晶体板结构尺寸(填充率为f，晶格常数L_x与L_y，L_a与L_b)，可高性能(快速、稳定性、收敛)地获得所需要的带隙特性，最终完成声子晶体板结构带隙设计，获得符合特定要求带隙特性的声子晶体板结构。

实施例不应视为对本发明的限制，任何基于本发明的精神所作的改进，都应在本发明的保护范围之内。