CN112836416A - 一种用于抑制弹性波传播的声子晶体结构优化设计方法 - Google Patents

Abstract

本发明公开了一种用于抑制弹性波传播的声子晶体结构优化设计方法。该方法首先采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系，然后通过移动渐近线法(MMA)对声子晶体单胞构型进行基于变密度法的拓扑优化设计，使其能在保持一定力学性能的条件下带隙最大化，并在变密度法的基础上，通过特征驱动方法进行拓扑结构边界光顺，使其具有明显的边界特征，便于CAD重构，结合数值仿真手段和拓扑优化设计技术，最终得到符合工程实际需要的新型声子晶体，从而达到提高设计效率、减少工作量和提高实际应用价值的目的。

Description

一种用于抑制弹性波传播的声子晶体结构优化设计方法

技术领域

本发明属于噪声与振动控制技术领域，具体涉及一种声子晶体结构设计方法。

背景技术

传统周期性结构易于在外部激励下产生以弹性波形式传播的振动和噪声，而声子晶体是具有抑制弹性波传播的一种特殊周期性新材料。能带色散关系是声子晶体最基本的特性之一，在能带色散关系中对所有波矢均不出现色散曲线的频率段称为带隙。在带隙频率和波矢范围内，由于不存在相对应的振动模式，因此能对该频率和方向的弹性波起到隔振的效果。

声子晶体的传统设计思路有两种：一种是仿生设计，通过自然界中已有的结构进行仿生设计得到单胞构型，再通过经验进行晶体结构再设计，得到具有较好带隙特性的单胞；另一种是基于遗传算法的优化设计，通过模拟自然界生物群体中个体的行为模式来搜索最优解的方法，即给出目标期望，通过反演的逆向设计，给出最逼近带有目标物理特性的拓扑结构。其中，仿生设计存在单一无规律性、难以实现系统性设计的缺点；而遗传算法存在方法效率低，容易出现过早收敛等问题，不利于高效大规模设计。

中国专利CN103218529A《一种二维固一固声子晶体XY模态拓扑优化方法》通过快速平面波展开法计算二维固/固声子晶体XY模态的色散关系，获取相应的带隙值，然后应用遗传优化算法，根据带隙所要达到目标，搜索声子晶体最优材料拓扑布局，但改方法人为的施加了对称约束，导致优化结果可能不是最优的，并且该方法没有针对Z模态进行优化分析。并且遗传算法本身存在效率低下的缺点，优化得到的正方形像素型结构可制造性差，并且作为基体的环氧树脂刚度过小无法作为承力件。

发明内容

为了克服现有技术的不足，本发明提供了一种用于抑制弹性波传播的声子晶体结构优化设计方法。该方法首先采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系，然后通过移动渐近线法(MMA)对声子晶体单胞构型进行基于变密度法的拓扑优化设计，使其能在保持一定力学性能的条件下带隙最大化，并在变密度法的基础上，通过特征驱动方法进行拓扑结构边界光顺，使其具有明显的边界特征，便于CAD重构，结合数值仿真手段和拓扑优化设计技术，最终得到符合工程实际需要的新型声子晶体，从而达到提高设计效率、减少工作量和提高实际应用价值的目的。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤1：采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系；

步骤2：基于变密度法的声子晶体结构设计；

以声子晶体单胞体分比约束以及声子晶体单胞内部和单胞之间的连接性约束为约束条件，以声子晶体每个单元的相对密度作为拓扑优化过程中的设计变量，采用移动渐近线法对声子晶体结构设计的优化模型定义如下：

其中，N为声子晶体单元的个数，即设计域中有限单元的数目；ρ为声子晶体单元伪密度向量，ρ_e为第e个声子晶体单元伪密度；f(ρ)为目标函数；G^H为等效剪切模量，G^*为预设的单胞等效剪切模量下限值；

和

分别是声子晶体单胞施加周期性边界条件对应的不同边界上单元的伪密度值，B是边界连接性参数值，用以实现单胞间的连接；V表示单胞体分比，V_f代表预设的单胞体分比上限值；

目标函数f(ρ)被定义为相对带隙BandGap，同时表征带隙宽度和带隙位置，数学表达式为：

式中k为波矢，ω^* _n+1(k)与ω^* _n(k)分别为n+1和n阶标准化特征频率；

求解优化模型式(1)得到基于变密度法的声子晶体单胞优化结构；

步骤3：基于特征驱动的声子晶体结构设计；

以声子晶体结构特征的几何参数作为拓扑优化的特征设计变量，对步骤2得到的基于变密度法的声子晶体单胞优化结构进行改进，光滑声子晶体单元结构边界，得到声子晶体的基于特征驱动的优化模型为：

式中，m是特征数目；d表示特征设计变量所组成的向量，表征特征在拓扑优化中的位置布局和尺寸信息，d_i表示第i个特征设计变量；

表示特征设计变量d_i的上限，向量d_i表示特征设计变量d_i的下限；

采用四边形固定网格单元中心处的密度值来代替整个声子晶体单元的密度值，表达式为：

式中，x_e表示声子晶体单元e的中心点坐标，Φ(x_e)表示声子晶体单元结构边界，

表示对声子晶体单元结构边界光滑处理后的Heaviside函数，表达式为：

式中，Δ称为光顺因子，表示材料模型光顺边界的宽度的一半；α为非零常数；

求解优化模型式(3)，得到最终的声子晶体单胞优化结构；

步骤4：采用有限元方法进行能带计算以及传输分析计算，对最终的声子晶体单胞优化结构进行验证。

优选地，所述步骤1中采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系的具体方法如下：

定义离散形式的声子晶体动力学方程为如下形式的特征方程：

(K-ω²M)U＝0 (6)其中K为声子晶体刚度矩阵，M为声子晶体质量矩阵，U为声子晶体单胞各节点位移向量，U＝[U_左上 U_左下 U_右下 U_右上 U_左 U_右 U_上 U_下 U_内部]^T，U_左上、U_左下、U_右上、U_右下、U_左、U_右、U_上U_下和U_内部表示正方形单胞边界上四个顶点、四边中点以及内部中心点的位移向量；

根据Bloch定理，单胞边界满足如下关系：

式中，a_x和a_y为分别声子晶体单胞的基矢；k_x和k_y分别为波矢k在x和y方向的分量；

写成矩阵形式为：

其中，

式(8)定义了声子晶体单胞对应边界位移场的关系，被称为Floquet边界条件，是Bloch定理在单胞边界上的形式；

将式(8)代入特征值方程式(6)，并在两端同乘H^T，得到：

DU＝0 (11)其中

D＝H^T(K-ω²M)H (12)

矩阵D中包含了与波矢k有关Floquet边界条件，因此方程(12)被称为声子晶体单胞的频散方程，只需波矢k扫略不可约布里渊区边界并求解频散方程(12)即得到声子晶体单胞的能带色散关系。

优选地，所述步骤2中求解优化模型式(1)时，采用对低密度区的SIMP插值模型进行改进的MSIMP材料插值模型，表达式如下：

E＝ρ^penalE₀

其中，E为声子晶体单元弹性模量，E₀为声子晶体材料弹性模量，ρ₀为声子晶体材料密度，ρ为插值得到的声子晶体单元密度，ρ^penal是惩罚因子；

采用密度投影方法抑制和避免灰度单元的产生，利用基于双曲正切函数的光滑Heaviside公式：

式中，

表示通过密度过滤之后的声子晶体单元e的密度值，β为控制投影函数的陡峭程度参数，η∈[0,1]为阈值。

优选地，所述声子晶体结构特征的几何参数包括声子晶体单元中心位置坐标值、倾斜角、尺寸。

本发明的有益效果如下：

1、针对使用固定网格变密度法进行动力学连续体结构拓扑优化时，出现的棋盘格、灰度单元以及局部模态等数值不稳定现象，本发明能消除棋盘格现象；采用在SIMP材料插值模型中引入惩罚因子和密度投影法来抑制和避免灰度单元的产生；

2、本发明采用对低密度区的SIMP插值模型进行适当改进的MSIMP材料插值模型能有效地处理局部模态问题；

3、本发明采用特征驱动进行边界光顺，提高CAD可重构性，使其具有易于制造的拓扑构型。

4、本发明所用的优化方法与遗传算法相比，具有计算效率高，系统性强，计算结果可靠的优点，并且使用单一材料优化得到的单胞具有全局最优解，即带隙性能最好的特点。

附图说明

图1为本发明实施例中声子晶体单胞拓扑优化设计流程图。

图2为本发明实施例中有限元法正方形单元示意图。

图3为本发明实施例中SIMP和MSIMP模型插值函数比曲线。

图4为本发明实施例中基于变密度法的XY模态声子晶体构型设计第三阶带隙优化结果：(a)优化结果单胞(b)优化结果4×4结构(c)优化结果相对带隙。

图5为本发明实施例中基于变密度法的Z模态声子晶体构型设计第三阶带隙优化结果：(a)优化结果单胞(b)优化结果4×4结构(c)优化结果相对带隙。

图6为本发明实施例中第三阶最优带隙与刚度约束的关系：(a)XY模态(b)Z模态。

图7为本发明实施例中基于特征驱动的XY模态声子晶体构型设计第三阶带隙优化结果：(a)优化结果单胞(b)优化结果4×4结构。

图8为本发明实施例中基于特征驱动的Z模态声子晶体构型设计第三阶带隙优化结果：(a)优化结果单胞(b)优化结果4×4结构。

图9为本发明实施例中弹性波ΓX方向传输谱计算的示意图。

图10为本发明实施例中横波输入和纵波输入对ΓX和ΓM两个方向的传输谱：(a)横波输入下ΓM方向的传输谱(b)纵波输入下ΓX方向的传输谱(c)横波输入下ΓX方向的传输谱(d)纵波输入下ΓX方向的传输谱。

图11为正方形单胞示意图。

具体实施方式

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

为克服现有声子晶体微结构设计的盲目性、效率低下、可制造性差等缺点；同时，为了避免使用固定网格变密度法进行动力学连续体结构拓扑优化时，出现棋盘格图案、灰度单元以及局部模态等数值不稳定现象。本发明利用有限元方法分别计算二维声子晶体XY模态和Z模态的色散关系并获得相应的相对带隙，然后通过变密度法和特征驱动方法进行优化设计，获得具有最大相对带隙的声子晶体单胞结构。

一种用于抑制弹性波传播的声子晶体结构优化设计方法，包括以下步骤：

步骤1：采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系；采用施加Floquet边界条件来模拟周期性结构中弹性波场的情况，并由有限元法计算各个波矢下的特征频率，从而获得周期性结构的XY模态和Z模态的能带色散关系。

步骤2：基于变密度法的声子晶体结构设计；

变密度法适用于各向同性材料，以每个单元的相对密度作为拓扑优化过程中的设计变量，从而将拓扑优化问题转化为设计域内材料最优分布设计问题。

以声子晶体单胞体分比约束以及声子晶体单胞内部和单胞之间的连接性约束为约束条件，以声子晶体每个单元的相对密度作为拓扑优化过程中的设计变量，采用移动渐近线法对声子晶体结构设计的优化模型定义如下：

其中，N为声子晶体单元的个数，即设计域中有限单元的数目；ρ为声子晶体单元伪密度向量，ρ_e为第e个声子晶体单元伪密度；f(ρ)为目标函数；G^H为等效剪切模量，G^*为预设的单胞等效剪切模量下限值；

和

分别是声子晶体单胞施加周期性边界条件对应的不同边界上单元的伪密度值，B是边界连接性参数值，用以实现单胞间的连接；V表示单胞体分比，V_f代表预设的单胞体分比上限值；

目标函数f(ρ)被定义为相对带隙BandGap，同时表征带隙宽度和带隙位置，数学表达式为：

式中k为波矢，ω^* _n+1(k)与ω^* _n(k)分别为n+1和n阶标准化特征频率；

求解优化模型式(1)得到基于变密度法的声子晶体单胞优化结构；

步骤3：基于特征驱动的声子晶体结构设计；

以声子晶体结构特征的几何参数作为拓扑优化的特征设计变量，对步骤2得到的基于变密度法的声子晶体单胞优化结构进行改进，光滑声子晶体单元结构边界，得到声子晶体的基于特征驱动的优化模型为：

式中，m是特征数目；d表示特征设计变量所组成的向量，表征特征在拓扑优化中的位置布局和尺寸信息，d_i表示第i个特征设计变量；

表示特征设计变量d_i的上限，向量d_i表示特征设计变量d_i的下限；

采用四边形固定网格单元中心处的密度值来代替整个声子晶体单元的密度值，表达式为：

式中，x_e表示声子晶体单元e的中心点坐标，Φ(x_e)表示声子晶体单元结构边界，

表示对声子晶体单元结构边界光滑处理后的Heaviside函数，表达式为：

式中，Δ称为光顺因子，表示材料模型光顺边界的宽度的一半；α为非零常数；

求解优化模型式(3)，得到最终边界清晰的声子晶体单胞优化结构；

步骤4：采用有限元方法进行能带计算以及传输分析计算，对最终的声子晶体单胞优化结构进行验证。

优选地，所述步骤1中采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系的具体方法如下：

定义离散形式的声子晶体动力学方程为如下形式的特征方程：

(K-ω²M)U＝0 (6)其中K为声子晶体刚度矩阵，M为声子晶体质量矩阵，U为声子晶体单胞各节点位移向量，以正方形单胞为例，单胞上的点按位置分共有九类，如图11所示。U＝[U_左上 U_左下 U_右下 U_右上 U_左 U_右U_上 U_下 U_内部]^T，U_左上、U_左下、U_右上、U_右下、U_左、U_右、U_上U_下和U_内部表示正方形单胞边界上四个顶点、四边中点以及内部中心点的位移向量；

根据Bloch定理，单胞边界满足如下关系：

式中，a_x和a_y为分别声子晶体单胞的基矢；k_x和k_y分别为波矢k在x和y方向的分量；

写成矩阵形式为：

其中，

式(8)定义了声子晶体单胞对应边界位移场的关系，被称为Floquet边界条件，是Bloch定理在单胞边界上的形式；

将式(8)代入特征值方程式(6)，并在两端同乘H^T，得到：

DU＝0 (11)其中

D＝H^T(K-ω²M)H (12)

矩阵D中包含了与波矢k有关Floquet边界条件，因此方程(12)被称为声子晶体单胞的频散方程，只需波矢k扫略不可约布里渊区边界并求解频散方程(12)即得到声子晶体单胞的能带色散关系。

优选地，所述步骤2中求解优化模型式(1)时，采用对低密度区的SIMP插值模型进行改进的MSIMP材料插值模型，表达式如下：

E＝ρ^penalE₀

其中，E为声子晶体单元弹性模量，E₀为声子晶体材料弹性模量，ρ₀为声子晶体材料密度，ρ为插值得到的声子晶体单元密度，ρ^penal是惩罚因子；

采用密度投影方法抑制和避免灰度单元的产生，利用基于双曲正切函数的光滑Heaviside公式：

式中，

表示通过密度过滤之后的声子晶体单元e的密度值，β为控制投影函数的陡峭程度参数，η∈[0,1]为阈值。

优选地，所述声子晶体结构特征的几何参数包括声子晶体单元中心位置坐标值、倾斜角、尺寸。

具体实施例：

本发明力图从微观尺度上分析结构对弹性波的控制原理及特性，通过对结构单元的创新设计，制备出具有天然材料所不具备的超常物理性能的新材料，提高对弹性波传播的控制和操纵能力。

如图1所示，本实施例优化算法采用MMA算法。拓扑优化设计域取a＝1m的正方形区域，有限元分析基于平面应变假设。材料选取硅，材料属性设置为弹性模量E＝18.5GPa，泊松比v＝0.27，密度ρ＝2330kg/m3。设计域以50×50的正方形网格进行离散划分。

包括以下步骤：

1、如图2所示，按照步骤1用有限元方法在MATLAB软件中编写二维能带色散关系计算程序。

2、基于变密度法的声子晶体构型设计。

在特定体分比约束下，以单胞相对带隙最大化为目标函数，施加单胞内部和单胞之间的连接性约束，对声子晶体的单胞进行优化设计。对于要满足体分比约束和单胞内单胞间连接性约束的蜂窝型声子晶体拓扑优化问题，其优化模型如式(1)

下面对目标函数、约束以及数值不稳定现象的处理方法进行详细的说明。

首先，为了消除单胞尺寸的影响，对能带图中的特征频率进行标准化处理，记ω*＝ωa/2πc，c为材料横波波速。带隙被定义为：

BandGap＝minω^* _n+1(k)-maxω^* _n(k) (16)

式(16)中仅仅描述宽窄的带隙称为绝对带隙，实际应用对带隙的要求是尽可能地低频和宽。因此，用带隙中心的频率值来表征带隙的位置，用绝对带隙和带隙中心频率来定义相对带隙，从而得到目标函数：

其次，本实施例包含三个约束条件，分别为单胞内部连接约束、单胞胞间连接约束和体积约束。

单胞内部连接约束是为了解决在优化过程中会出现单胞内固体材料分布不连续的问题，其表达为：G^H≥G^*。

单胞胞间连接约束要求单胞对应边界上有限单元伪密度值相等，在实际优化过程中，考虑数值计算中浮点数无法比较相等，放宽该约束，并考虑后续灵敏度计算，将其改为：

式中B为控制参数，其取值影响优化的结果：如果B取值太小，优化过程中该约束将很难被满足，容易陷入局部最优解；若B取值太大，那么该约束将起不到保证单胞胞间连接的作用。体积约束即体分比约束，是单胞中固体材料体积与单胞总体积的比值，对于蜂窝型声子晶体等效于质量约束能在优化过程中有效的控制结构总体质量，使得结果的质量小于预定值。体积为所有单元伪密度之和，因此体积对某一单元伪密度的导数应恒等于1，即体积约束灵敏度为：

如图3所示，本实施例采用对低密度区的SIMP插值模型进行适当改进的MSIMP(Modified Solid Isotropic Microstructures with Penalization)材料插值模型，其表达式如式(14)。

本实施例采用密度投影方法来抑制和避免灰度单元的产生。除了在SIMP材料插值模型中引入惩罚因子外，同时采用，利用基于双曲正切函数的光滑Heaviside公式。随着β值的不断增大，投影函数曲线上数值在η附近的中间密度区域就会变得更加狭窄。当β→∞时，单元e的密度值

就会趋近于离散的0-1解。

其中，MSIMP插值模型中取ρ^penal＝3来抑制中间密度单元的产生。同时，对于密度投影法，取密度投影阈值η＝0.5。声子晶体带隙性能拓扑优化问题具有高度非凸性，设计过程中容易落入局部解中，为了保证优化算法的稳定性，避免最终优化结果过于依赖于初始值，可对密度投影法中的参数β采用渐进式优化策略(Continuation approach)。β初值取为β₀＝4，当步长小于一定值时β＝β+0.5，最终β值取为β_max＝64。渐进式优化策略能有效地减少拓扑优化结果中灰度单元的出现和保证拓扑优化的稳定性。

如图4-图6所示，通过面内模式和面外模式的拓扑优化结果可以发现，具有宽频带隙的蜂窝型弹性超材料呈现出集中质量-韧带的结构形式。为了验证这种结构形式在宽频带隙的蜂窝型弹性超材料中是否具有普遍性，通过改变拓扑优化过程中的刚度约束，对面外模式和面内模式第三阶带隙进行优化，得到刚度约束-最优刚度的关系，两种模式下最优的带隙值和刚度大致呈线性负相关，即随着刚度的增大，最优带隙值减小。

3、基于特征驱动的声子晶体构型设计

变密度法得到的拓扑优化结果具有重构困难的缺点，特征驱动的拓扑优化方法从CAD特征设计理念出发，可以很好的解决重构困难的问题。特征驱动的拓扑优化方法以结构特征的几何参数作为拓扑优化的设计变量，通过特征的形状优化和位置布局实现结构的拓扑变化，从而实现优化结构的力学性能和工程特征一体化设计。

声子晶体带隙性能的特征驱动拓扑优化问题的数学模型表示为式(8)。

由于固定网格边界和特征结构边界不重合，本发明采用四边形固定网格单元中心处的密度值来代替整个单元的密度值，α通常取为非零的、足够小的常数，以避免空材料所导致的刚度矩阵奇异。实质是将结构的水平集函数模型转化为可以进行有限元求解的连续材料模型。

如图7和图8所示，由于基于特征驱动的拓扑优化方法得到的结构边界轮廓清晰，不会产生棋盘格现象和大量灰度单元，无需用密度过滤法等方法进行处理，但动力学拓扑优化中的局部模态现象仍存在，同样地采用MSIMP材料插值模型进行规避。本实施例采用超椭圆作为优化设计中的基元，以基元的几何参数作为拓扑优化设计变量，从而通过改变基元的布局和形状优化实现结构的拓扑优化设计。超椭圆设计变量取值范围见表1。每个基元几何参数包括中心位置坐标x、y，长半轴长度a，短半轴长度b以及倾斜角θ。Heaviside函数中参数α＝1×10^-5，光顺因子Δ＝0.1。

表1超椭圆设计变量取值范围

4、优化结果能带计算与传输谱分析。

能带色散关系中的带隙表征的是无限大周期性结构的理想情况下的隔振情况，而实际应用中的结构都是有限尺寸的，这时需要用传输谱来表征隔振情况，也能从传输谱中各频段内振动的具体衰减量来验证能带带隙理论的正确性。为了验证优化结果的能带色散关系与带隙性能，利用弹性波不能在在声子晶体中带隙范围内传播的特性，采用有限元法对有限结构的传输谱进行计算。本发明采用有限元软件COMSOL Multiphysics，借助有限元软件COMSOL Multiphysics的Solid Mechanics模块来计算传输谱，分别计算横波输入和纵波输入对ΓX和ΓM两个方向的传输谱，对比两个方向上的传输频谱，共同点为在理论计算的带隙内横波和纵波均出现传输谷，即输出端的振幅急剧下降，说明弹性波在此频段传播受到抑制，这也从另一面说明了带隙理论计算的正确性。传输谱数值上看，在ΓX方向上的隔振效果要明显优于ΓM方向。

如图9和图10所示，通过对优化结果的有限元分析可以看出，以上优化方案得到的声子晶体具有低频宽带隙的特点，XY模态的归一化带隙宽度为0.44左右，相比优化前提高了121.5％，Z模态的归一化带隙宽度为0.7左右，相比优化前提高了112％，这种结构设计方案为设计特定带隙频段的声子晶体提供了一种高效、系统性的方法，且优化构型具有可制造性好，承载能力强的特点，可以应用于一般承力构件的减振降噪。

Claims (4)

1.一种用于抑制弹性波传播的声子晶体结构优化设计方法，其特征在于，包括以下步骤：

步骤1：采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系；

步骤2：基于变密度法的声子晶体结构设计；

以声子晶体单胞体分比约束以及声子晶体单胞内部和单胞之间的连接性约束为约束条件，以声子晶体每个单元的相对密度作为拓扑优化过程中的设计变量，采用移动渐近线法对声子晶体结构设计的优化模型定义如下：

其中，N为声子晶体单元的个数，即设计域中有限单元的数目；ρ为声子晶体单元伪密度向量，ρ_e为第e个声子晶体单元伪密度；f(ρ)为目标函数；G^H为等效剪切模量，G^*为预设的单胞等效剪切模量下限值；

和

分别是声子晶体单胞施加周期性边界条件对应的不同边界上单元的伪密度值，B是边界连接性参数值，用以实现单胞间的连接；V表示单胞体分比，V_f代表预设的单胞体分比上限值；

目标函数f(ρ)被定义为相对带隙BandGap，同时表征带隙宽度和带隙位置，数学表达式为：

式中k为波矢，ω^* _n+1(k)与ω^* _n(k)分别为n+1和n阶标准化特征频率；

求解优化模型式(1)得到基于变密度法的声子晶体单胞优化结构；

步骤3：基于特征驱动的声子晶体结构设计；

以声子晶体结构特征的几何参数作为拓扑优化的特征设计变量，对步骤2得到的基于变密度法的声子晶体单胞优化结构进行改进，光滑声子晶体单元结构边界，得到声子晶体的基于特征驱动的优化模型为：

式中，m是特征数目；d表示特征设计变量所组成的向量，表征特征在拓扑优化中的位置布局和尺寸信息，d_i表示第i个特征设计变量；

表示特征设计变量d_i的上限，向量d _i表示特征设计变量d_i的下限；

采用四边形固定网格单元中心处的密度值来代替整个声子晶体单元的密度值，表达式为：

式中，x_e表示声子晶体单元e的中心点坐标，Φ(x_e)表示声子晶体单元结构边界，

表示对声子晶体单元结构边界光滑处理后的Heaviside函数，表达式为：

式中，Δ称为光顺因子，表示材料模型光顺边界的宽度的一半；α为非零常数；

求解优化模型式(3)，得到最终的声子晶体单胞优化结构；

步骤4：采用有限元方法进行能带计算以及传输分析计算，对最终的声子晶体单胞优化结构进行验证。

2.根据权利要求1所述的一种用于抑制弹性波传播的声子晶体结构优化设计方法，其特征在于，所述步骤1中采用有限元方法计算声子晶体XY模态和Z模态下的能带色散关系的具体方法如下：

定义离散形式的声子晶体动力学方程为如下形式的特征方程：

(K-ω²M)U＝0 (6)其中K为声子晶体刚度矩阵，M为声子晶体质量矩阵，U为声子晶体单胞各节点位移向量，U＝[U_左上 U_左下 U_右下 U_右上 U_左 U_右 U_上 U_下 U_内部]^T，U_左上、U_左下、U_右上、U_右下、U_左、U_右、U_上U_下和U_内部表示正方形单胞边界上四个顶点、四边中点以及内部中心点的位移向量；

根据Bloch定理，单胞边界满足如下关系：

式中，a_x和a_y为分别声子晶体单胞的基矢；k_x和k_y分别为波矢k在x和y方向的分量；

写成矩阵形式为：

其中，

式(8)定义了声子晶体单胞对应边界位移场的关系，被称为Floquet边界条件，是Bloch定理在单胞边界上的形式；

将式(8)代入特征值方程式(6)，并在两端同乘H^T，得到：

DU＝0 (11)

其中

D＝H^T(K-ω²M)H (12)

矩阵D中包含了与波矢k有关Floquet边界条件，因此方程(12)被称为声子晶体单胞的频散方程，只需波矢k扫略不可约布里渊区边界并求解频散方程(12)即得到声子晶体单胞的能带色散关系。

3.根据权利要求1所述的一种用于抑制弹性波传播的声子晶体结构优化设计方法，其特征在于，所述步骤2中求解优化模型式(1)时，采用对低密度区的SIMP插值模型进行改进的MSIMP材料插值模型，表达式如下：

E＝ρ^penalE₀

其中，E为声子晶体单元弹性模量，E₀为声子晶体材料弹性模量，ρ₀为声子晶体材料密度，ρ为插值得到的声子晶体单元密度，ρ^penal是惩罚因子；

采用密度投影方法抑制和避免灰度单元的产生，利用基于双曲正切函数的光滑Heaviside公式：

式中，

表示通过密度过滤之后的声子晶体单元e的密度值，β为控制投影函数的陡峭程度参数，η∈[0,1]为阈值。

4.根据权利要求1所述的一种用于抑制弹性波传播的声子晶体结构优化设计方法，其特征在于，所述声子晶体结构特征的几何参数包括声子晶体单元中心位置坐标值、倾斜角、尺寸。

Images