Topology optimization of microstructures with perturbation analysis and penalty methods

Topology optimization at the continuum nano/microscale is of wide interest in designing and developing more efficient micro/nano electromechanical systems. This paper presents a new methodology for topology optimization of microstructures that is based on perturbation analysis and the penalty methods. The homogenized material coefficients are numerically computed based on perturbation analysis, and periodic boundary conditions are imposed by the penalty methods. The sensitivity analysis is implemented directly without the adjoint method. The extension of the proposed method to the design of components for multi-field analysis is straightforward. The capability and performance of the presented methodology are demonstrated through several numerical examples.


Introduction
Recent advances in additive manufacturing have allowed for complex architected materials. Geometrically sophisticated devices or structures at the micro/nanoscale, such as micro-electro-mechanical systems, are designed with thinfilm fabrication techniques (Villanueva and Maute 2014). Even direct manipulation of individual atoms to design materials in the field of nanotechnology has been possible (Chen et al. 2020). Topology optimization can be used to find the optimal material distribution or topology within a given domain so as to improve the performance, durability, and efficiency of materials or structures. On the other hand, topology optimization provides unrestricted and heuristic designs of micro/nanostructures. To fully utilize the new manufacturing approaches, it is critical to design microstructural materials with desired properties, which are determined by the topologies and compositions of the microstructure.
Homogenization is commonly used to compute the effective parameters for materials with complex microstructures. Topology optimization of material microstructures is the inverse problem of homogenization, aiming to find the interior topology of a base cell with specific properties (Sigmund 1994), such as negative Poisson's ratio (Clausen et al. 2015;Vogiatzis et al. 2017), maximum shear modulus (Du et al. 2017;Huang et al. 2011), maximum bulk modulus (Huang et al. 2011;Gibiansky and Sigmund 2000), maximum buckling strength (Thomsen et al. 2018), etc. Popular topology optimization methods included the density-based approach (Andreassen et al. 2011;Xia and Breitkopf 2015;Gao et al. 2019), level-set method (Wang et al. 2003;Vogiatzis et al. 2017;Gao et al. 2018;Nguyen et al. 2020), and BESO method (Huang and Xie 2009;Huang et al. 2011;Radman et al. 2013), to name a few. When the material consists of a periodically repeated microstructure, the asymptotic homogenization method is a good choice to extract the homogenized material properties. It corresponds to the energy-based approaches that employ average stress and strain theorems (Sigmund 1994). The homogenization method is often combined with an adjoint method (Bendsoe and Sigmund 2003) to conduct the sensitivity analysis, which can be difficult with increasing complexity for instance for multi-field problems.
In this study, we propose a new methodology for topology optimization of microstructures using perturbation analysis and penalty methods. The homogenized material coefficients are efficiently computed with multiple load cases based on perturbation analysis (Kaczmarczyk et al. 2008), which avoids complex mathematical derivations. Periodic boundary conditions are imposed with the penalty methods. The sensitivity analysis is implemented directly without the adjoint method. The main advantage of the method is that it can be easily extended to the multi-field coupling material design problems of first-order and even second-order homogenization. The derivation of the proposed method is more straightforward than the traditional method using the asymptotic expansion. In terms of optimization problems related to the stiffness tensor, this method may have no advantage. However, the adjoint method is used in the traditional direct manner, which is not intuitive for the multi-field coupling problem, especially the second-order homogenization problem. In addition, the calculation of the homogenized material parameter (piezoelectric tensor, dielectric tensor) based on the asymptotic expansion method is quite complex.
The remainder of the paper is organized as follows: Sect. 2 introduces the proposed method including a comparison with the energy-based method; Sect. 3 presents the formulations for the design of piezoelectric materials while Sect. 4 describes the numerical implementation. Three numerical examples are used to demonstrate the capability of the methodology in section 5 before Sect. 6 concludes our manuscript.

Homogenization with perturbation analysis
In the FE 2 (Feyel and Chaboche 2000) approach, the microscopic displacement field u m can be written as the sum of a macroscopic field and a periodic fluctuation field, where M is the macroscopic strain, y is the local base cell spatial coordinate, and r u is the displacement microfluctuation field. The volume average of the microstrain field is equal to the macroscopic strain, where the symmetric (1) u m = M ⋅ y + r u operator is dropped since it has been assumed (Kaczmarczyk et al. 2008).
The last integral of the RHS of Eq. (2) can be transformed into a boundary integral by the Gauss divergence theorem yielding where R , L , T , and B represent the right, left, top, and bottom boundaries of the base cell, respectively. Periodic boundary conditions finally lead to the equation where [[u]] denotes the jump in the displacement field ( Fig. 1), where the matrix W is determined by the node coordinate y: y 1 0 0 1 2 y 2 0 1 2 y 3 0 y 2 0 1 2 y 1 1 2 y 3 0 0 0 y 3 0 1 2 y 2

Fig. 1 Periodic boundary conditions
The Hill-Mandel macro-homogeneity condition (Hill 1963) requires the increment energy equivalence between the microscale and the macroscale: The macroscopic stress is obtained by The stress-strain relationship at the macroscale is not known a priori, and the linearized relation between stain increments and stress increments are The above homogenized stiffness tensor C H klmn ( C H ) in this study is computed by perturbation analysis. The so-called two-index notation is introduced for the material tangent matrix replacing kl or mn by i or j, respectively. After discretization of the base cell, the homogenized stiffness tensor can be computed as where u e ( j M ) is the solution of nodal unknowns corresponding to the unit test strain; I i n is the i-th row of the identity matrix I n , where n is equal to three (2D) or six (3D). B , C e are the shape function gradient and the material stiffness tensor, respectively.

Imposition of the periodic boundary conditions using the penalty methods
The variational formulation of the penalty methods and periodic boundary conditions is given by (Henyš et al. 2019;Nguyen et al. 2014;Sanders et al. 2012): (7) where is the penalty parameter, and u , v are the displacement variables of the solution and trial functions, respectively. With the interpolation on the boundary, the solution and trial displacement gap functions are approximated as: where N + and N − are the shape function matrices. To illustrate the imposition of periodic boundary constraints, we partition the displacement vector U into the interior displacements and displacements on the boundary Γ + and Γ − . The global system of equations can be obtained as where K b is the bulk stiffness matrix, and K p and F are the penalty coupling matrix and external force vector generated by penalty methods, respectively. The penalty parameter is chosen to make the pivot element of K p be orders of magnitude greater than the pivot element of K b . The above matrices are given by

Sensitivity analysis
The modified SIMP approach (Solid Isotropic Material with Penalization) (Andreassen et al. 2011) is applied here, where the design domain is discretized by square finite elements, and each element is assigned a pseudo-density e that determines its stiffness tensor C e : where C 0 is the material stiffness tensor, f 0 = 1 , and f min is a very small factor assigned to void regions to prevent the singularity of the stiffness matrix, and p is a penalization factor introduced to reduce the intermediate densities. The mathematical formulation of the optimization problem reads as follows: where c is the objective function, K is the global stiffness matrix defined in (14), U s and F s are the global displacement and force vectors corresponding to the unit test strain, respectively; v e , V, and are the element volume, the base cell volume, and the prescribed volume fraction, respectively. The sensitivity of the objection function c∕ e is computed based on the perturbation analysis and penalty methods. By differentiating (19), we obtain: where The penalty coupling matrix K p and external force vector F s are independent of the element density variable. G is the transformation matrix forming the global stiffness matrix, K k b and U k s are the element stiffness matrix and the vector containing the degrees of freedom (DOFs). The sensitivity of the homogenized stiffness tensor C H ij ∕ e consists of two terms: The first term on the right side of (24) can be written as (21) in the second term on the right side of (24) yields with The computation of J can be prepared in advance, and l is identical for each element and only needs to be calculated once in each optimization loop reducing the computational cost. Hence, the sensitivity of the homogenized stiffness tensor is given as:

Comparison
In the classic asymptotic homogenization method, the homogenized stiffness tensor is given by the volume integration over the base cell (Sigmund 1994) where * (kl) pq is the periodic solution to the variational problem where 0(kl) pq corresponds to the unit test strain fields on the base cell, and v is the Y-periodic admissible displacement field. There is an equivalent method to predict the effective material properties named the energy-based approach (Sigmund 1994;Xia and Breitkopf 2015;Gao et al. 2018), where the unit test strain fields are imposed on the boundary of the base cell (Sigmund 1994), and the homogenized stiffness tensor is calculated as After discretization of the base cell, (31) is approximated by is the element displacement solution corresponding to the unit test strain fields, and k e is the element stiffness matrix. The sensitivity of the homogenized stiffness tensor C H ijkl ∕ e is computed using the adjoint method (Xia and Breitkopf 2015), k 0 being the element stiffness matrix for elements with e = 1 . Periodic boundary conditions are often imposed in a direct manner for the above energy-based method in 2D (Xia and Breitkopf 2015) and 3D (Gao et al. 2018), where the boundary constraints are separated into the corner, edge, and face (3D) constraints to prevent overconstraints, then the redundant unknown freedoms are eliminated. It is simpler to impose the periodic boundary conditions using the penalty methods as constraints are applied only on the boundary elements, especially for 3D problems, while the direct approach requires dealing with additional constraints on the corner or the edge as shown in Fig. 2. However, the penalty method is variationally inconsistent, and discrete results are sensitive to the penalty parameter (Sanders et al. 2012).
A comparison of the proposed method and the energybased approach is found in Table 1. For the energy-based approach, the imposition of periodic boundary conditions is not limited to the direct manner, and the penalty methods and Lagrange multipliers can also be adopted. With the adjoint method, the formulations of homogenization and sensitivity analysis are more concise. However, the proposed method is more straightforward in derivation and has hardly been applied before. More importantly, it can be easily extended to other multi-field coupling microstructure optimization issues, which will be illustrated in the next section.
Extension to the design of piezoelectric materials

Homogenization
The piezoelectric constitutive relations are given by: with E = − ∇ . C , e , and indicating the elastic stiffness tensor, the piezoelectric tensor, and the dielectric tensor, respectively; D , E , and are the electric displacement vector, the electric field vector, and the electric potential, respectively. The microscopic displacement u m and m can be written as the sum of a macroscopic field and a periodic fluctuation field, where E M is the macroscopic electric field while r is the electric potential microfluctuation field. After the volume average of m and E m , the following boundary conditions can be obtained,

Periodic boundary conditions lead to
where [[u]] and [[ ]] are the jump in the displacement field and the electric potential (Fig. 3), and the matrix W (6) and L are determined by the node coordinate y: The homogenized stress and electric displacement can be obtained by averaging: (41) L = y 1 y 2 (2D problems), The updates of the macrolevel stress and electric displacement are performed by the following incremental relations: The homogenized stiffness tensor C H , the homogenized piezoelectric tensor C H E ( C H D ), and the homogenized dielectric tensor C H DE are computed by perturbation analysis associated with independent unit test macroscopic strain and electric field. After discretization of the base cell, the above homogenized tensors for 2D problems can be computed as: where u e and e are the solutions of unknown freedoms corresponding to the unit test strain and electric field. B u and B are the displacement and electric potential shape function gradient, respectively. C e , e e , and e are the material coefficient tensors of the elements.
The use of perturbation analysis makes the homogenization method easy to extend to multi-field coupling problems. In

Fig. 3 Periodic boundary conditions
contrast, the formulations of the asymptotic homogenization method are more complex and need to be re-derived to deal with such problems (Silva et al. 1997(Silva et al. , 1999.

Periodic boundary conditions
Periodic boundary conditions are also imposed with the penalty methods. The formulation can be derived in a similar way to the previous section. By dividing the unknown vector U into the interior unknowns, displacements, and electric potentials on the boundary Γ + and Γ − , the global system of equations can be obtained as: where K b is the bulk coupling stiffness matrix, K pu , K p , F U , and F Φ are, respectively, the penalty coupling matrices and external force vectors generated by penalty methods. The above matrices are given by where and are the penalty parameters.

Sensitivity analysis
The material properties based on the SIMP model are determined by the design variable e : The optimization objective c is a function of the homogenized tensors C H , e H , and H . Taking the homogenized piezoelectric tensor e H (50) as an example, its sensitivity analysis is as follows.
The first term on the right side of (60) can be written as (21) in the second term on the right side of (60) yields where (57) The matrix K in (62) can be rewritten as Hence, the sensitivity of the homogenized piezoelectric tensor is computed as

Implementation
The base cell is discretized into nelx × nely rectangular elements in 2D, and the mesh grids of 3 × 3 elements are given in Fig. 4 and Fig. 5 to illustrate the design for the 2D material stiffness and piezoelectric cases. The degrees of freedom of one corner are fixed. The stopping criterion is the change in terms of design variables between two consecutive designs becomes less than 0.01. Optimality criteria (OC) method (Andreassen et al. 2011) is used to solve the optimization problem for the 2D and 3D cases, while the method of moving asymptotes (MMA) (Svanberg 1987) is adopted for the optimization of piezoelectric material. To ensure the existence of solutions to the topology optimization problem and to avoid the formation of checkerboard patterns, the density filter (Andreassen et al. 2011) is applied. The proposed method is implemented in Matlab. The procedures of the topology optimization for the homogenized stiffness tensor and piezoelectric materials are given in Table 2. The penalty coupling matrix, external force vector, and the integration of shape function gradient are all independent of the element density variable and can be calculated in advance. The transformation matrices are included in the matrices K b , M s , J , and s , which can be computed using the sparse matrices. Their rows and columns are prepared ahead, and the values are updated in the optimization loop. The processes barely change when extending the proposed method to the design of the piezoelectric materials, as no other treatments are applied except perturbation analysis and the penalty methods. For the first-order multi-field coupling and even second-order homogenization problems, the corresponding modifications are straightforward. On the other hand, the extension to the above cases requires complex derivations for the asymptotic expansion method and adjoint method, especially the adjoint method is not intuitive to conduct.

Numerical examples
This section presents three examples of designs of microstructured materials with the proposed method, namely materials with negative Poisson's ratio (2D), maximum bulk modulus (3D), and maximum hydrostatic coupling coefficient (piezocomposite). There are educational MATLAB codes targeting topology optimization of microstructures to achieve extreme material properties with the energy-based method in 2D (Xia and Breitkopf 2015) and 3D (Gao et al. 2019), so that the results of the proposed method and the energy-based method in the 2D and 3D examples (Xia and Breitkopf 2015;Gao et al. 2019Gao et al. , 2018 in terms of homogenization, sensitivity analysis, and optimal design are compared and discussed. Thereafter, the extensions of the method to piezocomposite are explored.

2D microstructured materials with negative Poisson's ratio
The base cell is discretized into 100 × 100 elements. The volume constraint is set to 0.5, p = 3 , the filter radius r min = 4 , and the penalty parameter is = 1000 , where the optimization parameter setups are the same as in the literature (Xia and Breitkopf 2015) to compare with the results of the energy-based method. The unit Young's modulus and the Poisson's ratio = 0.3 are used for computing the solid material stiffness tensor C 0 . To construct negative Poisson's ratio materials, the objection function is set as (Xia and Breitkopf 2015): The initial design is shown in Fig. 6, which also includes the homogenized stiffness tensor calculated by the perturbation analysis according to Eq. (11), the differences with the results of the energy-based homogenization approach, and their L 2 -norm. The displacement solutions corresponding to three test strains with penalty methods and the difference with the direct solution are given in Fig. 7. The relative differences of the homogenized stiffness tensor and displacements are all within 10 −5 , indicating the accuracy of the proposed perturbation analysis and penalty methods, respectively.
Based on the proposed method, the sensitivity for the initial design is shown in Fig. 8b. It is compared with the sensitivity from the energy-based method, and the relative difference is also within 10 −5 . The optimal designs of the two methods are very similar, differing only locally, and the difference in homogenized stiffness tensor is also small, as illustrated in Fig. 9. In addition, the method also allows designing materials with negative Poisson's ratio without introducing additional symmetry or isotropy constraints as in the energy-based method (Xia and Breitkopf 2015).
The above comparison shows that the perturbation analysis, penalty methods, and sensitivity analysis in the proposed method are very accurate, especially the sensitivity analysis is consistent without shortcuts (Wang et al. 2021). It should be noted that the Poisson's ratio of the optimized design is − 0.52 based on parameter setups consistent with the literature. Theoretically, Poisson's ratio of the optimal microstructure can approach the limiting value of −1 by adjusting the initial design, optimization algorithm, and parameter setups. However, the core of the proposed method lies in homogenization Table 2 The optimization procedures Homogenized stiffness tensors (2D and 3D) Piezoelectric materials Step 1: Define parameters for topology optimization, i.e., nelx, nely, (nelz), , p, r min Step 1: Define parameters for topology optimization, i.e., nelx, nely, , p, r min Step 2: Step 3: Loop of topology optimization Step 3

3D microstructured materials with maximum bulk modulus
The maximum of the material bulk modulus corresponds to the minimization of the following objective function: The base cell is discretized into 50 × 50 × 50 elements, = 0.5 , p = 5 , r min = 2 , and = 1000 are the solution parameters. Young's modulus and Poisson's ratio of the solid material are also selected as E = 1 and = 0.3 . The initial design is a cube with a matrix density of and a spherical inclusion of density ∕2 , as shown in Fig. 10a. For the initial design, the homogenized stiffness tensor and the (67) c = −(C 11 + C 12 + C 13 + C 21 + C 22 + C 23 + C 31 + C 32 + C 33 ) 1.89e-7 -2.38e-6 1.80e-9 -1.79e-9 1.80e-9 -1.46e-7 ×10 -4 comparison with the energy-based method (Gao et al. 2019(Gao et al. , 2018 are given in Fig. 10b, and the sensitivity c∕ e and its difference are illustrated in slices, as shown in Fig. 10c, d, respectively. All the relative differences are within 10 −5 . The optimal design without symmetry constraints is shown in Fig. 11, which is similar to the results in the literature (Chen et al. 2020;Sivapuram et al. 2018). Compared with the sensitivity analysis of the energy-based method, the vector l (26) needs to be solved in addition to the displacement u. The computation cost for different discretizations can be found in Fig. 12. The logarithm of the computation time is proportional to the mesh density, most of which is spent on solving for u and l. The disadvantage of the proposed method is that the computation cost is about twice the energy-based method for the presented 3D problem.

2D microstructured materials with maximum hydrostatic coupling coefficient
The hydrostatic coupling coefficient ( d h ) is an important performance characteristic for the design of low-frequency 1-3 piezocomposites such as hydrophones (Silva et al. 1999).
For the 2D plane strain microstructures in the 1-3 plane, the hydrostatic coupling coefficient can be written as (Silva et al. 1997): where The base cell is discretized into 100 × 100 elements; we choose the following parameters: = 0.5 , p = 3 , r min = 3 , = 10 5 , and = 10 5 . The initial guess of the material topology layout is a circle with the density of 0.001 in the base cell. The material matrices in (58) are given as follows, and the material properties are listed in Table 3.
Without imposing symmetry constraints, the optimal design, the homogenized material tensor, and the hydrostatic coupling coefficient d h are given in Fig. 13a, where the microstructure is a rotating structure. Due to the asymmetry  of the material model or the objective function, the density consistency ( + = − ) of the optimized microstructure on the boundary cannot be guaranteed, so additional symmetry constraints are sometimes required. When symmetry constraints are applied, the corresponding results are shown in Fig. 13b, where d h is much smaller. The optimal structure is axis-symmetric, and it contains two gray vertical bands due to the existence of gray elements, which can be avoided by the level-set method and needs further study.

Conclusions
This work proposes a new methodology to design material microstructures with perturbation analysis and the penalty methods. For 3D problems, it is more convenient to impose periodic boundary conditions with the penalty methods. Most importantly, this method can be easily and directly extended to multi-field coupling material design issues of first-order and even second-order homogenization. The numerical examples demonstrate the accuracy of the proposed method, where the relative differences of the homogenized effective parameters and sensitivity with the energy-based approach are within 10 −5 . Due to the use of perturbation analysis in homogenization, the sensitivity analysis of unknown freedoms to the design variables requires solving additional equations, increasing the solution time by twice. In addition, this methodology is based on the SIMP method, so there are some shortcomings such as gray elements, and further research is needed to combine the proposed method with the level-set method.