Deep learning aided topology optimization of phononic crystals

In this work, a novel approach for the topology optimization of phononic crystals based on the replacement of the computationally demanding traditional solvers for the calculation of dispersion diagrams with a surrogate deep learning (DL) model is proposed. We show that our trained DL model is ultrafast in the prediction of the dispersion diagrams, and therefore can be efficiently used in the optimization framework. The main novelty of the proposed approach relies on the use of non-uniform rational basis spline (NURBS) curves instead of pixels and/or mesh elements to control the shape of the unit cells of phononic crystals. The surrogate DL model is combined with a genetic algorithm serving as a topology optimization tool. The validity of the approach is shown in the case of phononic crystals made of a continuous matrix with cavities. Several objective functions have been tested as an alternative to the most common gap to mid-gap ratio. This allowed us to obtain interesting phononic crystal geometries which can be easily additively manufactured. The proposed method applies to problems involving inverse design and can open new avenues in the design of computer-assisted periodic structures.


Introduction
Phononic crystals and elastic metamaterials are well-known in literature to be a powerful means for achieving unconventional vibrational properties determined by the elastomorphic and material parameters of their unit cells.Properly conceiving their design in terms of size, shape, and arrangement, as well as choosing the disposition of their density and elastic properties demonstrated great potential for attaining exceptional dynamic behaviours, such as frequency bandgaps (BGs), wave steering, splitting and control, as well as negative refraction, to cite a few [1].In this context, frequency BGs, i.e., frequency ranges where the propagation of waves is strongly attenuated, received significant attention due to their potential applications [2] including but not limited to ultra-sensitive devices [3], large-scale periodic structures/foundations for vibration shielding [4,5], sub-wavelength imaging [6], elastic wave cloaking/lensing [7], etc.
To date, three main mechanisms exist to open frequency BGs, namely (i) the Bragg scattering [8], (ii) local resonances [9] and/or (iii) inertial amplification [10,11].In all these cases, the central frequency and the size of the BG are related to the shape and material properties of the unit cell.In the case of phononic crystals, two or more contrasting materials can be used to create heterogeneous periodic modulation of the elastic impedance, as well as a single solid material with air cavities (or voids) to achieve effective BGs.Pennec et al. [12] report a comprehensive overview of design examples and applications of two-dimensional phononic crystals.In the latter cases, resonators are responsible for the opening of the BGs.
In addition to the numerous designs based on the simple selection of geometry and the tuning of its geometrical/mechanical parameters proposed so far, optimization procedures have gained increasing interest in recent years, mainly focusing on maximizing the size of the BG [13][14][15][16][17][18][19][20].Gazonas et al. [13] used a genetic algorithm for the optimal design of phononic bandgaps in periodic elastic two-phase media (for the matrix and the inclusions).They used finite element formulation in which the unit cell is divided into triangular elements, each representing one phase of the material.The material choice for each triangle in a fixed mesh is encoded in a chromosome.They converted the generalized eigenvalue problem to the standard eigenvalue problem to reduce the computation time required to obtain (ω, k) pairs on which the fitness function is evaluated for each individual in the population at each generation.The problem of maximization of BG size for bending waves in a Mindlin plate with voids was considered in Ref. [14].The method of moving asymptotes was used in the optimization process.The 2D case of a phononic crystal made of solid material with voids was also investigated by Bilal and Hussein [15].They optimized unit cells by assigning material (silicon) and no-material (void) to pixels of a square image.They employed a specialized genetic algorithm to obtain optimized designs for out-of-plane, in-plane and combined outof-plane and in-plate elastic wave propagation.It should be noted that they focused on high-symmetry plain-strain square lattices achieving exceptional BG size normalized with respect to the mid-gap frequency exceeding 60%.In contrast, 2D asymmetrical phononic crystals were investigated by Dong et al. [16].The topology optimization was performed by using multiple elitist genetic algorithms with adaptive fuzzy fitness granulation.Multi-objective optimization approach was proposed by Hedayatrasa et al. [17] in which apart from the maximization of BG also in-plane stiffness of the phononic crystal unit cell was evaluated.Thick and relatively thin Polysilicon phononic crystal unit cells with square symmetry were studied and modelled by using the 3D finite elements.Selected Pareto front topologies were later manufactured by using a water jet and further investigated confirming the predicted BGs [19].The bidirectional evolutionary structural optimization was applied for the first time directly for the 3D shape optimization of phononic crystals by D'Alessandro et al. [18].The authors used 3D finite element models and adopted the model reduction technique to decrease the computational time for the calculation of the wave dispersion diagrams.
Recently, thanks to the increased computation power, the interest in topological optimization has considerably shifted from investigating phononic crystals with 2D unit cells to 3D unit cells.Nevertheless, the issues of (i) the manufacturability of the optimized geometries of phononic crystals and (ii) the computational time required to solve the inverse problem are still left open.
This results in the drawback of supplementing the current optimization process with additional steps so that the resulting geometry is feasible.For example, checkerboard suppression [27] or other filtering algorithms [17] must be implemented to avoid a design which will be difficult or unfeasible for manufacturing.The problem arises when mesh-based domain division is applied [13,18,20] or when the unit cell is divided into N × N pixels forming a logic matrix in which the value 1 represents one material and 0 represents another material [28].Another approach to tackle this problem is multi-objective optimization aiming at a wide BG as well as in-plane stiffness maximization [17].
We propose an alternative approach in which the shape of the cavity of the unit cell is represented by curves constructed using NURBS curves.This has numerous benefits such as better control of the design (suitable for instance for 3D printing), the possibility to set a minimal wall thickness and no need for any checkerboard suppression.
Another critical issue when dealing with topological optimization is that of reducing as much as possible the computation time at each iteration of the optimization procedure.So far, this problem has been partially alleviated by (i) implementing parallel computing, (ii) using a simpler model (adopting 2D plain strain or 2D plain stress assumptions, plate theories, approximate solutions, reduced order modelling, etc.) or (iii) homogenization techniques [21].However, parallel computing requires costly hardware and a large amount of energy to run it.The homogenization-based techniques have shown some limitations in the description of the phononic crystal mechanics at higher frequencies, including compromised accuracy.Finally, when simplified models are applied, not all the dispersion branches are computed [18].
To overcome these limitations, deep learning (DL) approaches have been adopted in the field of photonic crystals proposing to achieve the fitness function by replacing the dispersion diagram solver (forward model) with a surrogate DL model.This technique has been applied for the topology optimization of photonic crystals [22], Q-factor optimization of photonic nanocavities [23], transmission spectra [24] and band prediction [25] of photonic nanostructures.Nevertheless, DL-aided topology optimizations in phononic crystals are still scarce in the literature [26].Li et al. [27] designed phononic crystals with anticipated BGs through a DL-based data-driven method.An autoencoder-type neural network was trained to extract the topological features from sample images whereas finite element analysis was employed to study BGs of respective unit cells.The authors propose a multi-layer perceptron to establish the inherent relation between BGs and topological features.However, in this research, the prediction of the full dispersion diagrams was not considered.
Han et al. [28] proposed a multistep DL-based inverse design of phononic crystals for wave attenuation.Their approach includes the application of a convolutional neural network and generative adversarial network to establish mapping relationships between the structural topology and the band structure of the phononic crystals.They considered 1D phononic crystals composed of 2D unit cells and applied the semi-analytical periodic spectral finite element method to compute their complex band diagrams.The topological structures of the unit cell were randomly generated by using the method of twice mirror symmetry with discrete grids of 20 × 20 square blocks.The authors themselves concluded that such an approach is insufficient to unlock the potential of DL.
It should be noted that if dispersion curves are evaluated fast, it is possible to employ a reinforcement learning algorithm as presented in [29].The authors used a Q-learning algorithm to effectively and inversely design the structural parameters to maximize the bandgap width of the phononic crystal beam structure.

P. Kudela et al.
The aforementioned significant knowledge gap and open research issues related to high computation cost motivated us to develop a DL model capable of ultrafast prediction of dispersion diagrams in a wide frequency range for a given shape of a unit cell of a phononic crystal.
Here we present a novel approach for the topology optimization of phononic crystals replacing the computationally demanding traditional solvers for the calculation of the dispersion diagram with a surrogate DL model.We show that the developed DL model is ultrafast in predicting the dispersion diagrams but at the cost of the time for the supervised training.However, once the DL model is trained, it can be efficiently used in the optimization framework.Until now this methodology has only been applied in the field of photonics with a very limited number of examples.Moreover, our DL model is predicting the whole dispersion diagram instead of outputting lower and upper frequencies of the band gap.Additionally, our approach controls the cavity shape by NURBS curves instead of pixels and/or mesh elements.This aspect is crucial for the manufacturability of the optimized phononic crystals.Several objective functions were tested as an alternative to the most common gap to mid-gap ratio leading to unique designs.
The targeted application of the proposed DL-aided topology optimization of phononic crystals is the design of sensors for ultrasonic wave-based non-destructive testing (NDT) and structural health monitoring (SHM).Certain concepts of phononic crystal waveguide transducers for nonlinear elastic wave sensing [29], filtering and localization of the source of nonlinear waves were already demonstrated [3].These studies indicate that it is desirable to develop an objective function which results in a BG around the selected frequency.

Deep learning approach to forward modelling
The general concept of the proposed method is presented in Fig. 1.It consists of three main blocks: dataset computation, supervised learning and topology optimization.A large dataset of phononic crystal examples is computed by using the finite element method (FEM).The dataset is used for the supervised training of a surrogate DL model.It learns the mapping between inputs (cavity shapes) and outputs (dispersion diagrams).The trained DL model is then used for ultrafast prediction of dispersion diagrams in the topology Fig. 1.Flowchart of the proposed DL-aided topology optimization approach.
P. Kudela et al. optimization framework based on a genetic algorithm.An initial population of cavity shapes represented by NURBS curves is evaluated in parallel leading to initial fitness scores.The standard procedure of selection, crossover and mutation is performed to produce offspring of the next generation (new cavity shapes).Again, the trained DL model is used for the prediction of dispersion diagrams.Based on it the objective function value is calculated.The termination criterion is usually set by the maximum number of generations and the number of stall generations.
The process of dataset preparation, the development of a surrogate DL model and its performance are described in this section whereas the topology optimization of the unit cell is described in section 3.

General concept
The harmonic wave equation in a homogeneous and isotropic solid is given by: where λ and μ are Lamè constants, ρ is the mass density, u is the displacement vector, ω is the angular frequency and ∇ is the gradient operator.The band structure for an infinite phononic crystal can be extracted by considering the unit cell infinitely replicated in space and applying the Bloch conditions to the solution of Eq. ( 1) in the form: where u k (r) is a periodic function for the spatial position vector r with the same periodicity as the structure and k = (k x , k y ) is the Bloch vector.Substituting Eq. ( 2) into Eq.( 1) leads to the eigenvalue problem of the form: which must be accompanied by appropriate Floquet boundary conditions [8].The mass and stiffness matrices can be derived by using the finite element method (FEM) formalism.A detailed description of the finite element procedure can be found in Ref. [30].
The efficient solution of the eigenvalue problem of Eq. ( 3) represents one of the most challenging aspects of the topology optimization problem of phononic crystals.COMSOL Multiphysics is a very popular commercial software with an embedded function calculating the dispersion diagrams k(ω) of periodic structures.Although an excellent tool for the calculation of dispersion curves, their computation for the case of 3D unit cells may take some time depending on the mesh density (for instance, in the case of a mesh consisting of 14,072 tetrahedral elements and 134,040 degrees of freedom, the computation of 24 eigenmodes on a computer equipped with AMD Ryzen 7 1800X Eight-Core 3.8 GHz processor and 64 GB RAM took about 47 min).Therefore, it makes it unfeasible for topology optimization problems.
The idea here is to replace the FEM forward solver for computing dispersion diagrams with the surrogate DL model as schematically depicted in Fig. 2. Such a surrogate DL model can predict dispersion curves in a time shorter than 100 ms.It takes as an input binary Fig. 2. Schematic representation of the concept of the FEM solver replacement for computing dispersion diagrams by the surrogate DL model in the case of a phononic crystal made by cross-like holes in a matrix.The idea is that the DL model takes an image of the cavity as an input in which white pixels (ones) represent solid and black pixels (zeros) represent air.
P. Kudela et al. image in which white pixels represent solid material (the aluminium alloy in our test case) whereas black pixels represent cavity (air).Of course, the approach can be extended to two any materials with high impedance mismatch.Since the input to the DL model is an image, it cannot represent the whole 3D geometry as it is in the case of the FEM solver.Therefore, some assumptions and simplifications must be made.

Assumptions
It is assumed that the unit cell is made of a solid isotropic material (aluminium alloy) in which a cavity of a certain shape is formed.It is assumed that the cavity extends through the whole thickness of the unit cell.Moreover, the thickness-to-unit cell size ratio is fixed.In our case, the dimensions of the unit cell were as follows: • thickness 3 mm • unit cell width 8 mm.
Because we aim to make a proof of the concept, we have made further simplifications which however do not influence the generality of the proposed framework for topology optimization.The simplifications are related to the geometry of the cavity which reduces the computation time needed to generate a relatively large dataset of examples for supervised learning.Since the input to the DL model is encoded into a pixel-wise image, the DL model can take more complex geometry shapes than the simplified one presented in the next section.
First, for the generation of a dataset, we use only unit cells with a single cavity.Furthermore, the square lattice symmetry is taken into consideration which means that the design domain is 1/8 of the unit cell.The calculations of the dataset were performed by using COMSOL Multiphysics and covered the first irreducible Brillouin zone Г-X-M (see Fig. 3).

Dataset description
Cavity geometry is designed by using NURBS.A NURBS curve is defined by its order, a set of weighted control points, and a knot vector.NURBS curves are generalizations of both B-splines and Bézier curves, the primary difference being the weighting of the control points, which makes NURBS curves rational.In our implementation, we used a third-order curve which is represented by a quadratic polynomial.
The scheme of construction of cavity shapes within a unit cell by using NURBS and mirroring is presented in Fig. 3.The NURBS curve is spanned over four segments, and its shape is determined by the position of 9 control points P, 9 weights w and knot vector t.But some constraints regarding positions of control points are enforced so that the NURBS curve is within the design space.More details about the implementation of NURBS and cavity geometry are given in Appendix A. The shape of the NURBS curve is given by the cyan solid line in Fig. 3.However, we selected only 4 evaluation points per segment which resulted in 16 NURBS points denoted by diamond Fig. 3. Construction of cavity shapes within a unit cell by using mirrored NURBS curves; Г-X-M is the irreducible Brillouin zone for a unit cell with square lattice symmetry, superimposed here on the unit cell for the sake of synthetic representation; d is the minimum wall thickness measured from the edge of the unit cell.
P. Kudela et al. markers.The final geometry of the cavity is produced by appropriate mirroring of NURBS points.In the final step, a random scaling factor is applied to NURBS points.These points are then used for mesh generation, so the final geometry of the cavity is simplified by linear interpolation.In this way, the computation time by using COMSOL Multiphysics can be significantly reduced.On the other hand, the flexibility of geometry refinement is kept within NURBS representation and can be used during DL model evaluation or for the final manufacturing.
A 2D mesh of the unit cell generated by using GMSH software [31] is shown in Fig. 4. The GMSH software was used for mesh generation instead of COMSOL Multiphysics due to the ease of automation and slightly better quality of resulting meshes.The mesh was generated by using NURBS points defined in Fig. 3.In the next step, the mesh is imported to the COMSOL Multiphysics and extruded to the thickness direction so that the full 3D model is created.Then periodic boundary conditions are applied, and dispersion diagrams are computed.The COMSOL Multiphysics runs through the MATLAB Live-Link so that the consecutive meshes can be imported and the computation process can be automated.
The dataset consisting of 9000 unit cell geometries was generated in which the inputs are binary images and the outputs are dispersion diagrams stored as (k*, f) vector pairs in which k* is the reduced wave vector with wavenumber components k * = [k x a/π, k y a/π] with a being the lattice parameter (unit cell size) and f is the vector of frequencies in Hz.The dispersion diagrams are represented by 1464 points corresponding to frequency-wavenumber pairs.It covers about 24 modes.
Three types of cavity shapes examined are shown in Fig. 5 and were constructed by random NURBS manipulation (see Appendix A for more details).The cavities in the first row of Fig. 5 are spanned over horizontal and vertical directions, the cavities in the second row of Fig. 5 are spanned over diagonals and the last row has blot-like cavity shapes.The same number of examples (3000) was generated for each type of cavity shape so that the dataset is well-balanced.
The computation of the dataset in COMSOL Multiphysics took about one month on a standard desktop computer.Because the computation time needed to create the large dataset is very long, it would be worth exploring alternative methods which require fewer data for supervised training.Good candidates are physics-informed neural networks [32] or other similar approaches [33] in which partial differential equations along with boundary conditions describing the underlying physics are incorporated in the loss function used for the neural network training by using a stochastic gradient descent optimization algorithm.
The dataset is available in an open repository [34].

Dataset pre-processing
The synthetically generated dataset contains 9000 samples of unit cells that are classified into three equal subsets of 3000 samples each.The shape of a unit cell is represented by a binary image of the size of 512 × 512 pixels.As mentioned earlier, the shape of the unit cell is diagonally symmetrical.To reduce the complexity of the computation, we used only the upper left quarter of the unit cell shape with a size of 256 × 256 pixels for developing inputs to the DL model.Consequently, the total shape of the dataset is (3, 3000, 256, 256, 1).For training purposes, 2900 samples of unit cells were randomly selected from each subset, resulting in a training set with a shape of (8700, 256, 256, 1).The labels (ground truths) of the dataset are the dispersion diagrams (frequency versus wavenumber) values with a range exceeding 500 kHz (the maximum frequency in our dataset was equal to 668.6 kHz).It is important to mention that we have found that normalizing the frequency values to a range between (0, 1) gives better-predicted outputs.Fig. 4. Mesh for case no 1004 constructed by using NURBS points (diamonds) as shown in Fig. 3.

Deep learning model
In general, DL techniques such as artificial neural networks (ANNs) and convolution neural networks (CNNs) are considered universal approximation functions.In essence, the developed model is an approximation function, which can map the input (the shape of the unit cell) to an output of 1464 frequency values.All dispersion diagrams in the synthetically generated dataset have fixed wavenumber values, therefore only the frequency vector was used as an output (prediction).Accordingly, our developed DL model is a multi-output regressor.Fig. 6 presents the general architecture of the developed DL model, which is composed of two main parts: (1) the encoder (CNN) and ( 2) the dense layers (ANN).
The encoder is a CNN model that performs convolution operations followed by subsampling (pooling).Consequentially, the encoder is responsible for extracting the features from the shape of the unit cell.Such features contain all the information required regarding the dispersion diagram (the predicted output).The ANN part consists of several layers with different numbers of neurons and is responsible for mapping the extracted features by the encoder into 1464 predicted frequency values.
The task of hyperparameter optimization is crucial, especially, when developing a DL model to solve a complex task.Hence, the quality of a predictive model depends largely on its hyperparameter configuration.For this purpose, we employed the hyperband optimisation technique [35], which is an extension of the successive halving algorithm [36].The hyperband technique initially runs a random configuration of hyperparameters on a specific schedule of iterations (a small number of iterations) per configuration.Then it takes the best performers and again runs them through a larger number of iterations.This process is repeated until we get the best configuration of hyperparameters.
Table 1 presents the hyperparameters that were tuned with the hyperband technique, the initial range of values of the hyperparameters, and the optimised value.Regarding the pooling operation hyperparameter, convolution blocks {1, 2, 4, 5, 6, 7} have average pooling, and convolution block {3} has max pooling.The output dense layer has 1464 units representing the predicted  frequencies corresponding to the uniform sweep over wavenumbers used in COMSOL Multiphysics.Furthermore, the mean square error loss function was used as our objective loss function with the Adam optimizer.Additionally, the early stopping technique was used to reduce the model overfitting.

Results of the forward modelling
The results obtained by using the DL model were compared with the results obtained by using the FEM (COMSOL Multiphysics).The following errors were calculated for quantitative assessment of DL model predictions: mean squared error (MSE) and maximum frequency error ε max : where f is the frequency obtained by using the FEM and f i is the frequency vector obtained by using the DL model.The selected results for random test cases (cases no 1029, 4672 and 10425) are presented in Fig. 7.The DL model predictions are represented by blue circular markers whereas the FEM results are given by red triangular markers.It can be seen qualitatively that the results agree well, especially for the first eight propagating modes.For each example, the BG is properly predicted.The lowest frequency errors were obtained for the cross-like cavity shape (first column of Fig. 7) where ε max is 9.1 kHz and the highest errors are for blot-like cavity shapes (third column of Fig. 7) where ε max is 15.4 kHz.
The results of error analysis for all training and testing examples related to three types of cavity shapes are given in Table 2.The training dataset consisted of 2900 samples whereas the testing dataset consisted of 100 samples per cavity shape (cross-like, diagonallike and blot-like).It can be noticed that the errors are on a similar level irrespective of cavity shapes.Also, errors on the testing dataset are higher than on the training dataset, as expected.
In conclusion, the results are satisfactory, and the trained surrogate DL model can be used for BG size prediction.It is expected that further improvement in the accuracy of the DL model can be achieved by using a larger dataset for supervised learning.

Topology optimization of the unit cell
A genetic algorithm (GA) was employed for the problem of the topology optimization of the unit cell of phononic crystals.In particular, the implementation of GA in the Global Optimization Toolbox for MATLAB was used with default parameters such as crossover rate, mutation rate, elitism, etc. with the only changes applied to stopping criteria: • maximum no of stall generations: 40, • maximum no of generations: 50.
We also used 400 individuals in the initial population instead of the default 200 because the evaluation of the DL model can be realized more efficiently for the batch of inputs at once and leads to faster convergence.The DL model implemented in Keras was imported to MATLAB by using the Deep Learning Toolbox so that the optimization process was performed in one environment.
The total number of design variables describing the shape of the cavity was 21.The design variables were real numbers in a range between 0 and 1.The detailed implementation of the geometry representation by the design variables is given in Appendix A.

Objective functions
The Global Optimization Toolbox for MATLAB is implemented in such a way that the minimization problem is solved in which a custom objective function is delivered.We tested several versions of objective functions which are given in Appendix A. For further analysis, we selected three of them which were the most interesting and robust.
The first objective function is built based on the standard gap to mid-gap ratio: where f top is the minimum frequency at the top of the BG and f bot is the maximum frequency at the bottom of the BG.The next objective function was proposed for designing a unit cell with a middle of the BG which is around the selected frequency f sel : It is composed of two weighted parts: (i) gap to mid-gap ratio and (ii) relative difference between selected frequency and frequency corresponding to the middle of the BG.The weights were selected so that the preference is set towards the BG centered at the selected frequency.Such an objective function can be useful if we want to design a filtering device around a certain frequency for already fixed material properties and thickness to unit cell size ratio.Objective functions F 1 and F 4 take values in a range of 0-1.
The last objective function represents the cumulative BG of two successive BG separated by a narrow branch: where f max is the maximum frequency in our dataset.Often two successive BGs are separated by a non-propagating mode so that the effective BG obtained by using the objective function given by Eq. ( 8) can be larger than that with the use of the objective function given by Eq. ( 6).Additional conditions were added to all objective functions in the form of a penalty function for the invalid shape of cavity geometry.Specifically, if any coordinate of evaluated points is exceeding the design space (1/8th of the unit cell), the objective function takes value 1.

Results and discussion
The GA was run several times for each objective function type, and we selected two exemplary results of optimal geometries of unit cells and corresponding dispersion diagrams.Furthermore, the DL model predictions for optimized geometries were verified by using the FEM and both results are shown in consecutive figures for qualitative comparison.
The results for the objective function F 1 given by Eq. ( 6) are presented in Fig. 8.The DL model predictions of dispersion curves agree well with the FEM results.The relative BG size predicted by the DL model is about 48%.The BG is highlighted in grey.It should be noted that due to the heuristic nature of the GA, the optimal cavity shapes differ for each run of the GA.Nevertheless, the obtained BGs for each run are of similar size.
The results for the objective function F 4 given by Eq. ( 7) are presented in Figs.9-10 for selected frequencies f sel = 150kHz and f sel = 200kHz, respectively.The objective function is formulated to impose the centre of the BG at these selected frequencies which was successfully achieved.The DL model predictions of dispersion curves agree mostly well with the FEM results.The exception is in the second run of the GA shown in Fig. 9 where the DL model failed to predict properly multiple BGs but the bottom of the largest BG is still properly indicated.The relative BG sizes predicted by the DL model are not as big as it is in the case of the objective function F 1, but the largest BG is still quite wide reaching 41.1% in the best case.It should be noted, however, that the largest BG in all these cases has a neighbouring BG which is separated by one dispersion branch which could be related to standing waves.Hence, the effective BG can be wider.The effect of successive BGs is explored next.
The results for the objective function F 9 given by Eq. ( 8) are presented in Fig. 11.Two successive BGs predicted by the DL model are highlighted in grey.The dispersion curves are shown only for the FEM but with colour-coded polarizations which helps to distinguish between in-plane and out-of-plane modes.The polarization was calculated according to the formula [37]: where V is the volume of the unit cell u x , u y , and u z are the displacement components along x, y, and z axes, respectively.The colours of points of the dispersion curves in Fig. 11 vary from p = 0 (blue) to p = 1 (dark red).Thus, colours close to dark red indicate vibration modes that are predominantly polarized out-of-plane, while colours close to blue indicate modes that are predominantly polarized inplane.
It is clear from Fig. 11 that the mode in between the two successive BGs has dominant out-of-plane displacements and is the nonpropagating mode.The cumulative BG of the two successive BGs separated by non-propagating mode can be as wide as 134 kHz (see Table 3).These results are accounted for the DL model prediction.However, the FEM results shown in Fig. 11 indicate that the cumulative BG is even larger.The results for all cases discussed above are summarized in Table 3.The largest relative gap to mid-gap ratio equal to 48.1% was achieved for the case of the objective function F 1 for the second run of the GA.However, in terms of the cumulative BG of two successive BGs, the best case was obtained for the objective function F 9 and the second run of GA in which BG cum = 134 kHz.These values are shown in black bold font in Table 3.It should also be underlined that the frequencies of the middle of the BGs align well with   3).
It should be noted that the optimal geometries obtained by the GA tend to lead to very thin wall thickness measured from the edge of the unit cell which is indicated by d in Fig. 3.Such a situation is present e.g., in Fig. 11 where the cavity extends almost to the edge of the unit cell.Such geometry can lead to manufacturing issues or is simply unfeasible.Therefore, we conducted another study in which we run GA with imposed fixed selected wall thicknesses in a range from 0.2 mm to 1.4 mm.For simplicity, we performed calculations for only one objective function, namely F 1 given by Eq. ( 6).The results are gathered in Fig. 12 which shows that the increase in the wall thickness leads to the closing of the BG.Interestingly, the shape of the cavity is gradually changing with the increase in the wall thickness.In conclusion, there is a trade-off between the wall thickness and the relative BG size.

Conclusions
In conclusion, in this work, a new topology optimization method based on deep DL algorithms for the calculation of the dispersion band diagrams in phononic crystals was presented.The conventional formula-driven approach to determine the dispersion curves was

Table 3
Summary of results obtained by using various objective functions and DL model predictions where F i (x) is the objective function value, x is the vector of design variables, BG% is the relative gap to mid-gap ratio, BG cum is the cumulative BG width of two successive BGs excluding the frequency range of non-propagating branch between them.replaced by a surrogate DL model capable of predicting dispersion curves in less than 100 ms, whereas ordinary FEM calculations would require approximately 0.5 -1 h.Several geometries have been tested and inputted as binary images in which white pixels represent solid material and black pixels cavities (or air).This implied the need for specific assumptions on the in-plane size and the thickness of the unit cell for reconstructing the whole 3D geometry of the FEM model.The validity of the approach is shown in the case of phononic crystals made of a continuous matrix with cavities, the shape of which is controlled by NURBS curves instead of pixels and/or mesh elements.Several objective functions have been tested as an alternative to the most common gap to mid-gap ratio.
Our results may be of interest for all the problems involving inverse design and open new avenues in the design of computerassisted periodic structures.

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: co-author, Prof. Wieslaw Ostachowicz is an Editorial Board Member of the MSSP journal.

Fig. 5 .Fig. 6 .
Fig. 5. Nine examples of the examined cavity shapes.The cavities span over horizontal and vertical directions in the first row, over diagonals in the second row and blot-like in the third row.

P
.Kudela et al.

Fig. 8 .
Fig. 8. Phononic crystals and corresponding dispersion diagrams obtained during 2 runs of GA by using objective function F 1 .

Fig. 9 .
Fig. 9. Phononic crystals and corresponding dispersion diagrams obtained during 2 runs of GA by using objective function F 4 for the selected frequency f sel = 150kHz.

Fig. 10 .
Fig. 10.Phononic crystals and corresponding dispersion diagrams obtained during 2 runs of GA by using objective function F 4 for the selected frequency f sel = 200kHz.

Fig. 11 .
Fig. 11.Phononic crystals and corresponding dispersion diagrams obtained during 2 runs of GA by using objective function F 9 .
Fig. 7. Comparison of DL model predictions with FEM dispersion diagrams for random test cases 1029, 4672 and 10425.P.Kudela et al.

Table 2
MSE averaged over all training and testing examples; Units: [kHz 2 ].