CN103310049A - 一种二维固相声子晶体混合模带隙的优化方法 - Google Patents

Abstract

本发明属于声带隙材料设计领域，涉及一种二维固相声子晶体带隙的优化方法，包括：步骤一，初始化；步骤二，应用有限元法计算二维固相声子晶体XY模和Z模能带；步骤三，用声子晶体混合模带隙构造目标函数并计算遗传个体的适应度；步骤四，进行遗传操作，生成下一代种群，使种群向前进化，不断更新种群；步骤五，如果种群稳定，转下一步；否则返回步骤二；步骤六，以前面获得的最优结果对应的种群为初始种群，将每个声子晶体原胞扩展离散为2M×2M正方形单元结构，重新进行寻优计算迭代。本发明克服了传统声子晶体经验设计的局限性，提升了根据带隙需要主动设计声子晶体结构的能力，增强了声子晶体的可设计性，提高了声子晶体结构的性能。

Description

一种二维固相声子晶体混合模带隙的优化方法

技术领域

本发明属声带隙材料设计领域，涉及一种二维固相声子晶体带隙的优化方法，通过对二维固相声子晶体原胞材料拓扑分布进行优化设计，获得所需要的二维固相声子晶体XY和Z混合模带隙特征。

背景技术

在声子晶体带隙内，弹性波是不允许通过的。声子晶体的很多应用都是基于它的带隙特性，因此，通过人工设计特定频率范围的带隙非常必要。相对于三维声子晶体，二维声子晶体也具有相似的带隙特征，而且更易于人工制造。当弹性波在二维声子晶体（x,y）平面内传播时，可以解耦成面内耦合波模（XY模）和面外剪切波模（Z模）。如果所设计的二维声子晶体完全禁止某频率范围弹性波的传播，则XY模和Z模应同时在该位置具有带隙，即存在混合模带隙。二维声子晶体混合模带隙的存在与否及带隙大小都非常重要；而且，带隙越大，其应用范围就越广，应用价值也越高。所以，寻求具有最大带隙的声子晶体结构，是声子晶体理论研究的重点之一。

目前，已经有很多关于声子晶体带隙的研究，结果表明声子晶体带隙的位置和宽度受如下因素影响：声子晶体晶格类型和晶格常数，组成材料的拓扑布局，以及材料常数对比。通过选择不同的晶格类型，如正方晶格、三角晶格等，采用不同的散射体图形，如二维问题中，采用圆柱、正方柱等，调整散射体的几何参数以及散射体和基体材料参数等方法，可以增加带隙尺寸或打开新的带隙。也就是说传统的声子晶体设计是通过试错的方法来改进声子晶体带隙特性，最后得到了一些简单的声子晶体结构。然而，对于所得到的简单声子晶体结构是否具有进一步改进的余地，以便获得更优的带隙特性，是无法确定的。可以看到，经验性的总结和直觉式设计方法所能得到的晶格结构十分有限；另外，对于某些具有特定功能需求的声子晶体带隙，例如固态声子晶体XY和Z混合模带隙，传统设计法也更无能为力。这在很大程度上束缚了人们对声子晶体的应用。

发明内容

针对声子晶体结构设计中，经验性总结和直觉式设计方法受到的种种限制，本发明提出一种二维固相声子晶体带隙的优化方法，根据对二维固相声子晶体XY和Z混合模带隙的要求，自动寻找对应的声子晶体材料最优拓扑布局，得到具有最优带隙特性的新颖的声子晶体结构。

一种二维固相声子晶体混合模带隙的优化方法，其特征在于包括以下步骤：

步骤一，进行初始化操作，方法如下：

（1）将正方晶格声子晶体原胞设定为关于x、y轴反射对称和绕z轴90°旋转对称模型。声子晶体原胞被离散为M×M正方形单元结构，将每个单元配置所选定的弹性材料，此时，声子晶体的结构设计问题等价于单元中弹性材料的选择问题。

进行上述设定后仅需1/8部分即可描述整个原胞，如图1所示。由于只考虑两种材料，可采用二进制形式表示，从而每个声子晶体原胞可以用一个二进制字符串（染色体）来表示。原则上M值越大越好，但搜索规模往往呈指数形式增长，而且单元尺寸过小也给工艺制备带来了难度，因此，M的选取应折中考虑。

（2）随机生成N_pop个二进制数表示的染色体(每个染色体对应于一种声子晶体原胞)构成初始种群。

N_pop为种群规模，其值过大会增加计算量，过小会使种群不具有代表性。所以，N_pop的值一般根据具体问题通过测试获得。

步骤二，应用有限元法计算二维固相声子晶体XY模和Z模能带。

二维固相声子晶体XY模所满足的波动方程为：

$\begin{array}{c}p\frac{{\∂}^{2}u}{{\∂t}^2}=\frac{\∂}{\∂x}((2\mathrm{\μ}+\mathrm{\λ})\frac{\∂u}{\∂x}+\mathrm{\λ}\frac{\∂v}{\∂y})+\frac{\∂}{\∂y}\left(\mathrm{\μ}(\frac{\∂u}{\∂y}+\frac{\∂v}{\∂x})\right)\\ p\frac{{\∂}^{2}v}{{\∂t}^2}=\frac{\∂}{\∂x}\left(\mathrm{\μ}(\frac{\∂v}{\∂x}+\frac{\∂u}{\∂y})\right)+\frac{\∂}{\∂y}((2\mathrm{\μ}+\mathrm{\λ})\frac{\∂v}{\∂y}+\mathrm{\λ}\frac{\∂u}{\∂x})\end{array}---\left(1\right)$

二维固相声子晶体Z模所满足的波动方程为：

$p\frac{{\∂}^{2}w}{{\∂t}^2}=\frac{\∂}{\∂x}\left(\mathrm{\μ}\frac{\∂w}{\∂x}\right)+\frac{\∂}{\∂y}\left(\mathrm{\μ}\frac{\∂w}{\∂y}\right)---\left(2\right)$

其中，ρ为材料密度；λ和μ为材料弹性常数（拉梅常数）；u、v、w分别为x、y、z方向的位移。

应用有限元法求解方程组(1)和(2)。注意到方程(2)与流体波动方程具有相似性，可以分别借助多物理场耦合软件COMSOL Multiphysics3.5a中的力学模块和声学模块，根据问题的要求，对其脚本语言进行二次开发，得到matlab环境下的程序代码，启动COMSOL withMATLAB，从而可以分别求解方程组(1)和(2)，得到二维固相声子晶体XY模和Z模能带。

步骤三，用声子晶体混合模带隙构造目标函数并计算遗传个体的适应度。

研究发现：固相声子晶体XY模，两相邻能级间不一定都能形成完全带隙。但固相声子晶体Z模，每两相邻能级间的带隙都能打开。所以在对固相声子晶体混合模带隙优化时，先计算指定的XY模的带隙，然后将Z模的带隙和XY模的带隙作比较，从而确定混合模带隙。根据优化的目标，用声子晶体混合模带隙构造目标函数，然后以此目标函数值来度量遗传个体的适应度。

本拓扑优化的目标为二维固相声子晶体XY模带隙所决定的混合模带隙相对值最大，即：

$f=\max \left(\frac{{\min \mathrm{\ω}}_{i+1}\left(k\right)-{\max \mathrm{\ω}}_i\left(k\right)}{{\max \mathrm{\ω}}_i\left(k\right)}\right)---\left(3\right)$

式中，f为目标函数，k为波矢，i =1，2，…，maxω_i(k)和minω_i+1(k)分别为二维固相声子晶体XY模第i能级和第i+1能级所决定的混合模带隙下边缘本征频率的最大值和上边缘本征频率的最小值。

步骤四，进行遗传操作，生成下一代种群，使种群向前进化，不断更新种群，方法如下：

（1）进行选择操作：根据遗传个体适应度值的大小，采用精英选择与轮盘赌相结合机制；

（2）进行交叉操作：采用均匀交叉方式，随机选择染色体实施行交叉或列交叉，以增强算法的全局搜索能力；

（3）进行变异操作：个体变异采用位变异机制，合适地选择变异率是平衡进化趋于局部最优解和随机搜索的关键，故每隔一定进化次数将考察种群平均适应度值和最佳染色体适应度值，然后调整变异率。

根据“优胜劣汰,适者生存”自然选择和遗传进化机理，随着优化迭代的进行，最优个体的拓扑结构将越来越逼近目标要求的声子晶体原胞最优拓扑布局。

步骤五，如果种群稳定，转下一步；否则返回步骤二。

步骤六，以前面获得的最优结果对应的种群为初始种群，将每个声子晶体原胞扩展离散为2M×2M正方形单元结构，重新进行寻优计算迭代，直到再次找到最优结果，保存优化结果并退出计算。

为了加快优化设计收敛速度，同时使优化的声子晶体材料边界更加光滑，本发明采用两阶段优化方法：第一阶段的优化是指，经过步骤二～五循环迭代直到达到收敛条件的优化过程；第二阶段的优化是指，当第一阶段优化完成后步骤六进行的优化过程。

与现有技术相比，本发明的有益效果是：

本发明基于遗传算法，借助于两阶段优化设计方法，克服了传统声子晶体经验设计的局限性；根据对二维固相声子晶体XY和Z混合模带隙的要求，自动快速寻找到具有全局最优的声子晶体拓扑形状，极大地提升了根据带隙需要主动设计声子晶体结构的能力，增强了声子晶体的可设计性；通过优化设计，极大地提高了声子晶体结构的性能，并将进一步推动声子晶体的研究及其工程应用。

附图说明

图1是本发明中关于x、y轴反射对称和绕z轴90°旋转对称的声子晶体原胞10×10单元结构示意图和基因表达示意图；

图2是本发明所涉及的方法流程图；

图3是本发明实施例20×20单元离散的最优声子晶体3×3原胞图；

图4是本发明实施例40×40单元离散的最优声子晶体3×3原胞图；

图5是本发明实施例20×20和40×40单元离散的最优声子晶体能带图；

图6是本发明实施例遗传算法进化过程中的种群平均适应度及最优个体适应度变化图；

图7是本发明实例遗传算法进化过程中的最优声子晶体XY模第三能级和第四能级所决定的混合模带隙上下限及带隙宽度变化图。

具体实施方式

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

本实施例研究由两种材料构成的二维正方晶格固相声子晶体。材料A为钢，其密度ρ_A=7780kg/m3，材料拉梅常数分别为：λ_A=122GPa，μ_A=81GPa；材料B为环氧树脂，其密度ρ_B=1180kg/m³，材料拉梅常数分别为：λ_B=4.43GPa，μ_B=1.59GPa。

本实施例的方法流程图如图2所示，包括以下步骤：

步骤一，进行初始化操作：声子晶体原胞设为关于x、y轴反射对称和绕z轴90°旋转对称，此时仅需原胞的1/8即可描述整个原胞。声子晶体原胞被离散为20×20正方形单元结构，遗传个体二进制字符串长度对应为55。将每个单元随机配置所选定的钢或环氧树脂，采用二进制字形式表示，1和0分别表示材料钢和材料环氧树脂。根据以上设置，随机生成二进制字符串表示的染色体构成初始种群，初始种群规模N_pop=60。

步骤二，应用有限元法计算二维固相声子晶体XY模和Z模能带。

启动COMSOL with MATLAB，采用有限元法解方程组（1）和（2），计算每个遗传个体对应的二维固相声子晶体原胞XY模和Z模的能带，得到对应的带隙数值。

步骤三，用声子晶体混合模带隙构造目标函数并计算遗传个体的适应度。

研究发现：钢-环氧树脂声子晶体XY模前八个能级内，第三和第四能级以及第六和第七能级间的带隙比较容易打开，而且带隙很宽，而Z模的每相邻能级间的带隙都容易打开。所以本实施例在对钢-环氧树脂声子晶体混合模带隙优化时，先固定XY模的第三和第四能级间带隙，然后将Z模的带隙和XY模的第三和第四能级间带隙作比较，从而确定混合模带隙。此实例拓扑优化的目标设为钢-环氧树脂声子晶体XY模第三和第四能级间的带隙所决定的混合模带隙相对值最大，即：

$f=\max \left(\frac{{\min \mathrm{\ω}}_4\left(k\right)-{\max \mathrm{\ω}}_3\left(k\right)}{{\max \mathrm{\ω}}_3\left(k\right)}\right)$

式中，maxω₃(k)和minω₄(k)分别为二维固相声子晶体XY模第三能级和第四能级所决定的混合模带隙下边缘本征频率的最大值和上边缘本征频率的最小值。

然后以，此目标函数值来表示遗传个体的适应度。

步骤四，产生下一代种群：依次执行选择、交叉和变异遗传操作生成下一代种群，使种群向前进化，不断更新种群。交叉率为0.5、变异率根据收敛情况取值为0.05～0.5。

步骤五，如果种群稳定，转到步骤六；否则返回步骤二。

步骤六，以前面获得的最优结果对应的种群为初始种群，将每个声子晶体原胞扩展离散为40×40正方形单元结构，遗传个体二进制字符串长度对应为210。重新进行寻优计算迭代，直到再次找到最优结果，保存优化结果并退出计算。

第一阶段优化设计，经过110次循环迭代，得到离散为20×20正方形单元的二维钢-环氧树脂声子晶体原胞最优拓扑结构，其3×3原胞见图3。以第一阶段最优结果对应的种群为初始种群，将每个声子晶体原胞扩展离散为40×40正方形单元结构，进行第二阶段优化设计，再经过120次循环迭代，得到离散为40×40正方形单元的二维钢-环氧树脂声子晶体原胞最优拓扑结构，其3×3原胞见图4；可以发现两阶段最优声子晶体结构均为简单格子的正方晶格结构，钢(刚性材料)周期排列但并不相连，为散射体(图3和图4中黑色部分)，而环氧树脂(柔性材料)为基体(图3和图4中白色部分)；散射体均近似为正方形图案，但第二阶段优化结果对应的散射体边界更加光滑；两阶段最优声子晶体的散射体填充率分别为58%和57.25%，即材料构成比例差别不大。两阶段优化的最优声子晶体能带图见图5，图中横轴坐标为波矢，纵轴坐标为频率（单位：ωa/2πc_s,其中c_s为环氧树脂的剪切波波速，c_s=1160.8m/s）；由于采用1/8模型，用有限元法求二维钢-环氧树脂声子晶体能带时，波矢沿不可约布里渊区ΓXM所有边(即Γ→X→M→Γ方向)扫描；离散为20×20正方形单元的原胞最优拓扑结构的XY模及Z模能带分别为空心点图和实点图，而离散为40×40正方形单元的原胞最优拓扑结构的XY模及Z模能带分别为粗实线图和细实线图；两阶段优化结果对应的由二维钢-环氧树脂声子晶体XY模第三能级和四能级所决定的混合模带隙完全打开，且带隙相对宽度都为最大，而其他位置没有出现混合模带隙；两阶段优化结果的能带图形状相似，混合模带隙以下能带差别很小，而混合模带隙以上能带随着能级的增加，差别也在提高；混合模带隙下边为XY模第三能级，而上边为Z模第二能级，即此混合模带隙是由XY模和Z模共同构成；两阶段最优声子晶体混合模带隙下边缘最大值分别为0.6747和0.6734，上边缘最小值分别为1.532和1.5872，带隙宽度分别为0.8573和0.9138，可见经过第二阶段的优化，混合模带隙上边能级有所提高，从而增加了混合模带隙宽度。种群的平均适应度和最优个体的适应度随优化过程而变化，如图6所示，图中横轴坐标为遗传优化进化次数（单位：次），纵轴坐标为相对带隙值表示的适应度值；种群的平均适应度为种群每个个体适应度的平均值，其中适应度最大的个体为该种群的最优个体；由图6可以发现，遗传算法优化进化前期，进化速度很快，而在后期，特别是接近最优解时，速度很慢，两阶段最优原胞所对应的个体适应度分别为2.4207和2.5534。进化过程中声子晶体XY模第三能级和四能级所决定的混合模带隙上下限（分别指带隙上边缘最小值和带隙下边缘最大值）及带隙宽度（带隙上限与带隙下限的差值）变化见图7；由图7可以发现，随着优化的进行，混合模带隙的上边缘值趋于变大，下边缘值趋于变小，带隙的宽度在变大。综上所述，两阶段优化设计可以快速找到声子晶体的最佳拓扑结构，同时经过第二阶段的优化，又使得散射体材料边界更加光滑，优化目标（即相对带隙值）也进一步提高。

上述优化结果表明，随着迭代的进行，种群向前进化，最后得到了最优的二维钢-环氧树脂声子晶体结构，由钢-环氧树脂声子晶体XY模第三和第四能级间的带隙所决定的混合模带隙完全打开，且相对带隙值达到最大，实现了优化的目标；这种全局拓扑优化方法超越了经验设计的局限，实现了声子晶体的主动设计。

Claims (2)

1.一种二维固相声子晶体混合模带隙的优化方法，其特征在于包括以下步骤：

步骤一，进行初始化操作，方法如下：

（1）将正方晶格声子晶体原胞设定为关于x、y轴反射对称和绕z轴90°旋转对称模型；声子晶体原胞被离散为M×M正方形单元结构，将每个单元配置所选定的弹性材料；

经上述设定后，声子晶体的结构设计问题等价于单元中弹性材料的选择问题；仅需1/8部分即可描述整个原胞；只考虑两种材料，采用二进制形式表示，每个声子晶体原胞用一个二进制字符串来表示；M的值原则上越大越好，但搜索规模往往呈指数形式增长，而且单元尺寸过小也给工艺制备带来了难度，因此，M的选取应折中考虑；

（2）随机生成N_pop个二进制数表示的染色体构成初始种群；

N_pop为种群规模，其值过大会增加计算量，过小会使种群不具有代表性；所以，N_pop的值一般根据具体问题通过测试获得；

步骤二，应用有限元法计算二维固相声子晶体XY模和Z模能带；

二维固相声子晶体XY模所满足的波动方程为：

$\begin{array}{c}p\frac{{\∂}^{2}u}{{\∂t}^2}=\frac{\∂}{\∂x}((2\mathrm{\μ}+\mathrm{\λ})\frac{\∂u}{\∂x}+\mathrm{\λ}\frac{\∂v}{\∂y})+\frac{\∂}{\∂y}\left(\mathrm{\μ}(\frac{\∂u}{\∂y}+\frac{\∂v}{\∂x})\right)\\ p\frac{{\∂}^{2}v}{{\∂t}^2}=\frac{\∂}{\∂x}\left(\mathrm{\μ}(\frac{\∂v}{\∂x}+\frac{\∂u}{\∂y})\right)+\frac{\∂}{\∂y}((2\mathrm{\μ}+\mathrm{\λ})\frac{\∂v}{\∂y}+\mathrm{\λ}\frac{\∂u}{\∂x})\end{array}---\left(1\right)$

二维固相声子晶体Z模所满足的波动方程为：

$p\frac{{\∂}^{2}w}{{\∂t}^2}=\frac{\∂}{\∂x}\left(\mathrm{\μ}\frac{\∂w}{\∂x}\right)+\frac{\∂}{\∂y}\left(\mathrm{\μ}\frac{\∂w}{\∂y}\right)---\left(2\right)$

其中，ρ为材料密度；λ和μ为材料弹性常数（拉梅常数）；u、v、w分别为x、y、z方向的位移；

应用有限元法求解方程组（1）和（2），得到二维固相声子晶体XY模和Z模能带；

步骤三，用声子晶体混合模带隙构造目标函数并计算遗传个体的适应度，方法如下：

首先计算指定的XY模的带隙，然后将Z模的带隙和XY模的带隙作比较，从而确定最优的混合模带隙；根据优化的目标，用声子晶体混合模带隙构造目标函数，然后以此目标函数值来度量遗传个体的适应度；

拓扑优化的目标为二维固相声子晶体XY模带隙所决定的混合模带隙相对值最大，即：

$f=\max \left(\frac{{\min \mathrm{\ω}}_{i+1}\left(k\right)-{\max \mathrm{\ω}}_i\left(k\right)}{{\max \mathrm{\ω}}_i\left(k\right)}\right)---\left(3\right)$

式中，f为目标函数，k为波矢，i=1，2，...，maxω_i(k)和minω_i+1(k)分别为二维固相声子晶体XY模第i能级和第i+1能级所决定的混合模带隙下边缘本征频率的最大值和上边缘本征频率的最小值；

步骤四，进行遗传操作，生成下一代种群，使种群向前进化，不断更新种群，方法如下：

（1）进行选择操作：根据遗传个体适应度值的大小，采用精英选择与轮盘赌相结合机制；

（2）进行交叉操作：采用均匀交叉方式，随机选择染色体实施行交叉或列交叉，以增强算法的全局搜索能力；

（3）进行变异操作：个体变异采用位变异机制，合适地选择变异率是平衡进化趋于局部最优解和随机搜索的关键，故每隔一定进化次数将考察种群平均适应度值和最佳染色体适应度值，然后调整变异率；；

步骤五，如果种群稳定，转下一步；否则返回步骤二；

步骤六，以前面获得的最优结果对应的种群为初始种群，将每个声子晶体原胞扩展离散为2M×2M正方形单元结构，重新进行寻优计算迭代，直到再次找到最优结果，保存优化结果并退出计算。

2.根据权利要求1步骤二所述的一种二维固相声子晶体混合模带隙的优化方法，其特征在于：二维固相声子晶体Z模所满足的波动方程（2）与流体波动方程具有相似性，用有限元法求解二维固相声子晶体混合模所满足的波动方程（1）和（2）时，借助多物理场耦合软件COMSOL Multiphysics3.5a中的力学模块和声学模块，根据问题的要求，对其脚本语言进行二次开发，得到matlab环境下的程序代码，启动COMSOL with MATLAB，从而可以求解方程（1）和（2），得到二维固相声子晶体XY模和Z模能带。

Images