Topological Optimization of Phononic Crystal Thin Plate by a Genetic Algorithm

Genetic algorithm (GA) is used for the topological optimization of phononic crystal thin plate composed of aluminum and epoxy resin. Plane wave expansion (PWE) method is used for calculations of band gaps. Fourier displacement property is used to calculate the structure function in PWE. The crossover rate and the mutation rate are calculated according to the adaptive GA method. Results indicate that filling rates, symmetry, polymerization degree and material parameters are key factors for design of topological configurations. The relations between the key factors and different topologies are studied in detail.

The manufacturing issues can affect the final properties of PC. At present, manufacturing methods of metamaterials mainly includes Fused Deposition Modeling (FDM) 35,36 , Stereo Lithography Apparatus (SLA) 37,38 , Selected Laser Sintering/Melting (SLS/SLM) 39,40 , Direct Laser Writing (DLW) 41,42 and Polymer Injection Technology (PIT) 43,44 . The manufacturing issues can affect the final properties. However, theoretical test verification and functional device applications are fundamental issues for PC design regardless of manufacturing types. So, the current research is still focusing on theoretical test verification and functional device applications. Significant achievements have been made by combination of topological optimization with PC design. However, the topologies of PCs can be affected by the given parameters for design, e.g. the filling rates, symmetry, polymerization degree and material parameters. So, it is necessary to investigate how the given material parameters affect the final topology of PC. In this paper, the optimization design of the PC thin plate is studied when the influence factors are considered to obtain the maximum relative band gap based on the plane wave expansion method and the genetic algorithm. The unit cell is divided into N × N elements and each element is filled with epoxy resin or aluminum randomly. The optimal topological configuration of the PC plate can be obtained. The key factors, including filling rates, symmetry, polymerization degree and material parameter for the generation of band gaps, are correlated with different optimal topologies. By this design method on key factors including the filling rates, symmetry, polymerization degree and material, the optimal PC with much wider band gap features on given conditions can be obtained.

Model and Method
Analysis of the thin pC plate. A 2D PC thin plate is discussed. The lattice constant is a = 0.02 m. As we can see in Fig. 1(a), the unit cell is a square lattice which is divided into 10 × 10 elements. Epoxy resin and aluminum are chosen as the material of the substrate and the scatterer respectively. The material parameters of the epoxy and aluminum are given as follows: the density ρ e = 1180 kg/m 3 , the elasticity modulus E e = 0.435 × 10 10 Pa and the Poisson's ratio v e = 0.159 for epoxy resin; the density ρ a = 2730 kg/m 3 , the elasticity modulus E a = 7.76 × 10 10 Pa and the Possion's ratio v a = 2.87 for aluminum.
The study of the thin plate can be regarded as the plane stress state because the thickness of the plate is much less than the lattice constant. When the elastic wave propagates in the thin plate, the stress-strain relation can be expressed as 45 : www.nature.com/scientificreports www.nature.com/scientificreports/ The longitudinal wave equation is 45 : Bring the stress-strain relation into the wave equation: By the plane wave expansion, the characteristic equations are obtained: For the center element P 0 in Fig. 1(b), the Fourier expansion coefficient of the parameters is 0 when it is filled with the epoxy resin. When P 0 is filled with the aluminum, the Fourier expansion coefficient of the parameters can be expressed as: where g represents the parameters ρ, α, β, μ According to the Fourier displacement property, the Fourier expansion coefficient of the parameters of the element P r (P r ∈ P a , P a : the set of the elements including the aluminum) in Fig. 1(b): At the same time, the method combined with PWE and GA can be extended to multi-material PC structure. The computational costs can be increased according to 46 for optimization in multi-material case. When this method is extended to three-material, the stress-strain relation and wave equation remain unchanged. However, the material expression Eq. (1) need to be the linear combination of the material parameters and x i in Eq. (1) should changes to be 0 or 1 or 2. At the same time, Eq. (8), the Fourier expansion coefficient of the parameters which can influence the calculation of the characteristic equation Eq. (5), should change to be the combination of two other materials: The algorithm can also be extended to 3D PC optimization. When this method is extended to 3D PC, the governing equation Eq. (4) should be changed to add the vibration along z-direction according to 47 . So, the calculation of the band gap includes two conditions: in-plane and out-plane models. The calculation of the Fourier expansion coefficient of the parameters, Eq. (6) is also changed according to 48 in the plane wave expansion method. For optimization, the computational costs are obviously increased when PC changes from 2D to 3D structure according to 49 . the genetic algorithm. The genetic algorithm starts with initializing a population of Np = 20 chromosomes randomly. The chromosomes are then evaluated by the plane wave expansion method. The objective function is the relative band gap (RBG) width and the optimization can be simplified as: where X is the chromosome of the topological configuration, k is the wave vector, j is the number of band, f is the filling rate of the unit cell. GA operation includes selection, crossover and mutation after population initiation. The roulette wheel selection algorithm, called proportional selection algorithm, is mainly used in the optimization algorithm 50 . The probability of being selected for each individual is proportional to its fitness. In this paper, we use this method to select individuals with larger fitness as the parents for next generation in the selection operation. The algorithm of the roulette is as follows: www.nature.com/scientificreports www.nature.com/scientificreports/ The multi-point crossover is adopted according to the characteristics of the thin PC plate chromosome to make the population keep high diversity and avoid trapping the local optimum. Generation with new individuals can be produced by altering the values of some genes of the chromosome in the mutation operation. The mutation and the crossover parameter can directly affect the convergence of the algorithm. If the crossover rate Pc and the mutation rate Pm are fixed values, there are a lot of disadvantages like the destruction of the diversity of the population and the individuals with high fitness and the stagnant process. So the adaptive genetic algorithm proposed by Srinvivas 51 is adopted. The crossover rate P c and the mutation rate P m can change with the fitness. The P c and P m increases when the population traps in local optimum. For the individual with large fitness, the P c  www.nature.com/scientificreports www.nature.com/scientificreports/ and P m becomes small to make the individual move to the next generation. So the current mutation scheme from Srinvivas ensures the convergence of genetic algorithm while maintaining the diversity of population. The P c and the P m are defined as 51 :   www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/ where f max is the maximum fitness of the population, f avg is the average fitness of the population, f ′ is the bigger fitness of the two individuals in crossover operation, f is the fitness of the individual in mutation operation, P c1 = 0.9, P c2 = 0.6, P m1 = 0.1, P m2 = 0.001.
The genetic algorithm evolution flow chart is shown in Fig. 2.

Results and Discussions
The influence of the filling rate. The optimization of the relative band gap is performed when the filling rate is not fixed. The convergence criterion is the continuous equality of the best fitness and the average fitness.
In fact, we have tested to calculate 500 generations. But after 150 generations, the best fitness converges a stable value of 0.56. So only 200 generations are given in Fig. 3. For each time for our calculation, the generated different random state leads to the same optimal configurations. The convergence speed is rapid in the earlier stage of the evolution and it becomes slower in the following stage. The average fitness of the first generation is 0.02 so that the BG is very small or even did not exist for the most initial individuals produced randomly in the earlier stage of the evolution. But the best fitness of the first generation is 0.1 so that the band gap of the thin PC plate opens and can be found easily. The unit cell obtained from topological optimization is shown in Fig. 4(a). It can be seen that the unit cell is a simple lattice structure with square scatterer. The corresponding filling rate of the scatterer is 0.64. The corresponding band structure is plotted in Fig. 4(b). The absolute BG is ranged from 63.9 kHz to 113.9 kHz between the third and fourth bands. The optimal RBG is 0.56. To verify the current algorithm, the band structures of PC thin plate is recalculated by the finite element software Comsol when the filling rate is 0.64 in Fig. 4(c). Because of the periodicity of the structure, only one unit cell is calculated. Stress-free boundary conditions are imposed on  www.nature.com/scientificreports www.nature.com/scientificreports/ the free surfaces and the periodic boundary conditions on the two opposite boundaries of the unit cell according to the Bloch theorem. By sweeping the wave vectors in the first irreducible Brillouin zone, the Bloch calculation gives the eigen-frequencies. Then the band structure is obtained. The RBG is 0.58. The error is smaller than 0.05%.
The optimization above is studied when N = 10. For different N, the optimization results are calculated. As we can see in Fig. 5(a), the optimal unit cell for N = 8 is also a simple lattice structure with square scatterer. The corresponding filling rate of the scatterer is 0.56. The optimal RBG is 0.55 in Fig. 5(b). When N = 12, the optimal unit cell changes and the corresponding filling rate is 0.64. The optimal RBG is 0.51. With the calculations of band gap based on the obtained configurations with different element numbers, the wider band gap can be obtained when the element number is 10 × 10. So, we selected 10 × 10 as fixed mesh for optimization in current work. Fixed mesh strategy is widely used in shape and topology optimization, as shown in 52,53 .
When the material of the substrate or the scatterer is limited, it is necessary to obtain wider BG by PC plate with adaptive configuration. The optimization including the optimal configurations, the corresponding 3 × 3 unit cells and the corresponding band structures of the thin PC plates for different filling rate are studied and given in Fig. 6. The filling rate f is selected at regular intervals. When f = 0.2, the optimal configuration of the scatterer is a square with corners along the diagonal in Fig. 6(a,b). The band gap is ranged from 35.9 kHz to 40.9 kHz and the RBG = 0.13 in Fig. 6(c). When f = 0.32, the band gap is ranged from 36.5 kHz to 50.2 kHz and the RBG = 0.32 in   www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 6(f). When f = 0.44, the band gap is ranged from 56.3 kHz to 87.5 kHz and the RBG = 0.43 in Fig. 6(i). When f = 0.56, the band gap is from 62.5 kHz to 104.1 kHz and the RBG = 0.49 in Fig. 6(l).
The variation of the widest RBG for different filling rate is discussed. As shown in Fig. 7, the RBG increases firstly when the filling rate increases to 0.64. In the case where filling rate is less than 0.64, the reflection and refraction becomes stronger because of the increasing contact area between the aluminium and the epoxy resign when the filling rate increase. However, when the filling rate increases, the interaction between the adjacent scatterers becomes stronger so that the pass band becomes wider and the band gap becomes narrower in the case where filling rate is greater than 0.64.
The influences of the symmetry and material parameters. The variation of the RBG is discussed when the symmetry is changed. In initial optimization, filling rate is a parameter being selected to be state variable in topological optimization. Now, it is selected to be 0.36 randomly to study the effect of symmetry parameters on RBG. The first Brillouin zone is shown in Fig. 8. Three kinds of PCs with different symmetries are shown in Fig. 9 when filling rate is 0.36. For the four-axis symmetry PC in Fig. 9(a), the RBG is 0.33. The RBG of the biaxial symmetry PC in Fig. 9(b) is 0.27. For the uniaxial symmetry PC in Fig. 9(c), the RBG is 0.24. When the topological configuration of the filling material is highly symmetrical, wider band gap can be obtained. The width of the band gap decreases when it is symmetrical only along one or two axes.
Define the coefficient of the dispersion as the quantity of the blank elements in the scatterer to represent the extent of polymerization of the PC. As shown in Fig. 10, when the coefficient increases, the relative band gap decrease. It has no band gap when coefficient of the dispersion increase to 40. The structure has good band gap only when the filling material is highly concentrated. When the coefficient of the dispersion increases, the distance between the adjust scatterers decreases. The interaction between the adjacent scatterers becomes stronger so that the pass band becomes wider and the band gap becomes narrower. At the same time, the scatterer effective density becomes smaller when the dispersion coefficient is increased in constant filling rate. BG width decreases when the density ratio between the scatterer and the substrate decreases.
The influences of the material parameters including Young's modulus and density on the optimization results are considered when the filling rate is ignored. In Fig. 11(a), the optimal configuration which includes a square scatterer remains unchanged when the relative Yong's modulus increases. The optimized relative band gap becomes wider because of the high contrast Young's modulus. As shown in Fig. 11(b), the optimized relative BG increases when the relative density ρ 2 /ρ 1 increases. The corresponding optimal configurations are different.
When ρ 2 /ρ 1 = 3.2, the scatterer is a square with missing corners as shown in Fig. 12(a). As we can see in Fig. 12(b), the band gap is located between the third and the fourth bands. As for the ρ 2 /ρ 1 = 4.0, the optimal configuration consists a complete graph in center and quarters of the complete graph in four corners of the unit cell as shown in Fig. 12(c). The corresponding band gap is located between the seventh and the eighth bands in Fig. 12(d).

Conclusions
In this paper, the PC thin plate with maximum RBG is designed by the genetic algorithm in combination with plane wave expansion method. Filling rates, symmetry, polymerization degree and material parameters are key factors for design of topological configurations and can lead to different topologies of phononic crystals. When the filling rate is ignored, the optimal configuration is a unit cell with a square scatterer that the filling rate is 0.64. For given filling rate, the optimization is also considered. The topological relative band gap increases first and then decreases when the filling rate increases. The optimal configurations are also different for different filling rate. The structure has good band gap only when the filling material is highly concentrated. The effects of the material parameters on the optimization are also considered. The discussions can provide an effective tool for the design of the PC thin plate with desired band gap.