Eigenfrequency constrained topology optimization of finite strain hyperelastic structures

This paper incorporates hyperelastic materials, nonlinear kinematics, and preloads in eigenfrequency constrained density–based topology optimization. The formulation allows for initial finite deformations and subsequent small harmonic oscillations. The optimization problem is solved by the method of moving asymptotes, and the gradients are calculated using the adjoint method. Both simple and degenerate eigenfrequencies are considered in the sensitivity analysis. A well-posed topology optimization problem is formulated by filtering the volume fraction field. Numerical issues associated with excessive distortion and spurious eigenmodes in void regions are reduced by removing low volume fraction elements. The optimization objective is to maximize stiffness subject to a lower bound on the fundamental eigenfrequency. Numerical examples show that the eigenfrequencies drastically change with the load magnitude, and that the optimization is able to produce designs with the desired fundamental eigenfrequency.


Introduction
Topology optimization is extensively used in industry to guide early design processes. Unfortunately, existing formulations for topology optimization are primarily restricted to stiffness optimization of linear elastic structures. To allow for broader applicability, other objectives and nonlinearities should be considered.
Several studies combine nonlinear elasticity with topology optimization. In the work of Bruns and Tortorelli (2001), finite strain Green-Venant elastic structures and compliant mechanisms were designed. Different structures undergoing nonlinear elastic deformations were optimized for stiffness metrics by Buhl et al. (2000). This subject was revisited by Kemmler et al. (2005) where the contribution from geometrical nonlinearities was investigated, and by Wallin et al. (2018) where a neo-Hookean hyperelastic material model was used. Other examples of nonlinear structural topology optimization include the work of Sigmund (2001a, b) where electromechanical devices with geometrical nonlinearities were designed using multiple materials, and the more recent work of Ivarsson et al. (2018) where the objective was to design impact mitigating structures while considering the effect of transients and finite strain viscoplasticity. Bruns et al. (2002) designed structures exhibiting snap-through behavior by modifying the basic arc-length algorithm. A phase-field-based topology optimization formulation that incorporates multi-material structures and large deformations was developed by Wallin et al. (2015).
Few have investigated the influence of eigenfrequency constraints in the topology optimization of structures undergoing finite deformations. To the authors' knowledge, the work of Yoon (2010) is the only example, wherein the fundamental eigenfrequency optimization of prestrained structures was considered by utilizing an internal element connectivity parameterization in contrast to the standard density-based approach.
Topology optimization of linear elastic structures subject to eigenfrequency constraints was first proposed by Díaz and Kikuchi (1992), where structural reinforcements were introduced to increase the fundamental eigenfrequency. Various optimization formulations that consider eigenfrequencies have subsequently been suggested, e.g. the maximization of the fundamental eigenfrequency as done by Pedersen (2000) or the maximization of the gap between two eigenfrequencies as performed by Du and Olhoff (2007). In the work of Niu et al. (2009), a two scale optimization method was presented to maximize the fundamental eigenfrequency of multiscale structures. A level set-based topology optimization formulation has been proposed by Allaire and Jouve (2005) to maximize the fundamental structural eigenfrequency.
Unlike the aforementioned examples, our topology optimization incorporates frequency constraints of hyperelastic structures that undergo large deformations. The frequencies are obtained by linearizing the governing equations about the deformed state. This results in an eigenvalue problem which is numerically similar to the conventional linear case, the only difference being that the stiffness depends on the deformation. A neo-Hookean model is used to model the material response, which implies that both geometrical and material nonlinearities are included, unlike the formulation proposed by Yoon (2010) which only considers the former. This choice is made to accurately capture the response during finite deformations, since the conventional St. Venant material model is only applicable to large deformation with small strain. For additional discussion on hyperelastic models and topology optimization, see Klarbring and Strömberg (2013).
As pointed out in e.g. Seyranian et al. (1994), Du and Olhoff (2007), and Thore (2016), the occurrence of degenerate eigenfrequencies is rarely avoided in optimization studies that consider eigenfrequencies. The treatment of degenerate eigenfrequencies in topology optimization of linear elastic structures is well documented, where the nondifferentiability of degenerate eigenfrequencies is generally resolved by only considering directional derivatives, cf. Seyranian et al. (1994). This formulation is extended here to incorporate geometrical and material nonlinearities.
An often encountered issue when performing eigenfrequency topology optimization is the presence of spurious eigenmodes localized in nearly void regions. These spurious eigenmodes stem from the ersatz material and SIMP-penalization schemes. To eliminate this issue, Du and Olhoff (2007) suggested a nonlinear interpolation of the volume fraction, such that the mass in low volume fraction elements rapidly tends to zero. Low volume fraction elements are also known to cause convergence problems in the Newton-Raphson iterations when performing nonlinear finite element analyzes as these elements become highly distorted. To combat this problem, Wang et al. (2014) proposed an energy interpolation scheme which models these elements via small strain theory. In this work, we circumvent both of the preceding issues by implementing an element removal method technique proposed by Bruns and Tortorelli (2003), whereby the low volume fraction elements are removed from the finite element discretization.
We formulate the optimization problem as a maximization of the stiffness subject to maximum volume and minimum fundamental eigenfrequency constraints. From an engineer's perspective, this formulation is more relevant than the case where the objective is the fundamental eigenfrequency.
A well-posed topology optimization problem is formulated by using restriction in which a minimum length-scale is imposed on the design via a filter. The method of moving asymptotes (MMA) scheme is used to update the design, where the required gradients are computed using the adjoint sensitivity analysis. Both simple and degenerate eigenfrequencies are considered when performing the sensitivity analysis. The numerical examples illustrate the influence of the eigenfrequency constraints on the finalized topology.

Preliminaries
In this work, tensors and vectors will be written in bold faced symbols, e.g. C, and all finite element matrices will be denoted in bold sans serif symbols, e.g. K.
Lowercase subscripts e indicate quantities associated with finite element e o = 1, ..., n e , where n e denotes the number of elements in the discretization. The conventional finite element assembly operation is represented by . Since finite deformations are present, we distinguish between the undeformed, reference configuration o , and the deformed, current configuration . It is assumed that a smooth mapping ϕ, between the reference configuration and the current configuration exists, i.e. the position vector X, of each particle in the reference configuration is connected to its position x, in the current configuration via x = ϕ(X, t) = u(X, t) + X, where u denotes the displacement field. The local deformation is defined by the deformation gradient F = ∇ o ϕ = 1 + ∇ o u, where ∇ o is the material gradient operator on o and 1 is the second-order identity tensor. The Jacobian representing the volumetric change is given by J = det(F ) and it is assumed that the mapping ϕ, satisfies J > 0. The local deformation can also be described by the right Cauchy-Green deformation tensor C = F T F and the Green-Lagrange strain tensor E = 1 2 (C − 1). In absence of body forces, the equation of motion is where ρ o denotes the mass density in the reference configuration and the symmetric second Piola-Kirchhoff stress tensor S is related to the Cauchy stress σ , via S = J F −1 σ F −T . The boundary of the body in the reference configuration ∂ o , with the unit normal n o , consists of two complementary surfaces, ∂ ot and ∂ ou , over which the traction t o = F Sn o , and the displacement u = ϕ − X =ū, are prescribed, respectively. The starting point for the finite element formulation is the weak form, which is stated as finding the differentiable u that satisfies the boundary condition u =ū on ∂ ou , such that for all smooth virtual displacements that satisfy δu = 0 on ∂ ou . The virtual Lagrangian strain δE, introduced in Eq. 2, is defined as

Constitutive model
Our constitutive assumption uses the isotropic compressible neo-Hookean strain energy model in which the strain energy U NH is expressed as where K and G in the limit of infinitesimal strain correspond to the bulk and shear moduli. Since we assume a hyperelastic material model, the second Piola-Kirchhoff stress tensor is obtained from the strain energy as The fourth-order stiffness tensor associated with Eq. 4 becomes where a 1 = KJ 2 + 2G 9 J −2/3 tr(C), a 2 = − 2G 3 J −2/3 and a 3 = G 3 J −2/3 tr(C) − K 2 (J 2 − 1), respectively. The dyadic products ⊗ and ⊗ are defined such that (A⊗T ) : H = A · H T · T T and (A⊗T ) : H = A · H · T T , where A, T and H are arbitrary second-order tensors.

Filtration and penalization
As is usual in topology optimization, the binary valued material indicator function c(X) = {0, 1} is replaced by the continuous varying volume fraction design field z(X) ∈ [0, 1] that is subsequently penalized so that z → c. The use of z instead of c convexifies the design space making the optimization cost and constraint functionals differentiable. Unfortunately, this representation of the design still generates an ill-posed optimization problem resulting in so-called chattering designs 1 . To obtain a well-posed optimization problem, we impose a lengthscale restriction on z via the filter proposed by Bruns and Tortorelli (2003), which is based on studies and proofs presented by Bourdin (2001) and Bruns and Tortorelli (2001). Thus, volume fraction field z is replaced by the filtered volume fraction field where X ∈ e o . In our discretization, we assume z is parameterized to be piecewise uniform over the elements and hence so is ν. The kernel function ω is a cone filter function emanating from the element e o centroid position X e and enveloping the surrounding element s o centroid positions Y s that lie within the filter radius r of X e , cf. Bruns and Tortorelli (2003). The filter relation (7) can ultimately be expressed in discretized format as where Z is the [n e × n e ] filter matrix, z = [z 1 , z 2 , ..., z n e ] T and ν ν ν = [ν 1 , ν 2 , ..., ν n e ] T are the vectors of element volume fractions and filtered volume fractions. Elements possessing ν e (X) = 1 and ν e (X) = 0 indicate elements full and devoid of material, respectively. Unfortunately, the filtering produces gray regions consisting of elements for which ν e (X) ∈ (0, 1). To reduce the extent of these regions, the SIMP scheme-based on the work by Bendsøe (1989) is used, wherein we replace U NH with where Here δ 0 represents a lower threshold in the ersatz material model which ensures the well-posedness of the elasticity problem and p controls the level of penalization. It is clear that as χ(ν) → δ 0 the material mimics void, whereas as χ(ν) → 1 the material realizes the neo-Hookean strain energy function (4).

Total Lagrangian FE formulation
The displacement field u is interpolated in each finite element using the standard finite element polynomial shape functions N, to define the displacement field u, i.e. u(X) = N(X)a e (t) where a e (t) is a vector consisting of the nodal displacements for finite element e at time t. Differentiating u gives the accelerationü(X) = N(X)ä e (t), and virtual displacement δu(X) = N(X)δa e (t), fields. Using these approximations, the discrete virtual Lagrangian strain in Eq. 3 is expressed as where B(a e ) has been decomposed into constant B c , and linear B l (a e ) terms, cf. Crisfield (1993) for details 2 . Using the arbitrariness of the virtual nodal displacement and the finite element interpolations in Eq. 2 yields the finite element equation of motion where the mass matrix M, internal force vector F int , and external force vector F ext , are In Eq. 13, we emphasize that ρ = νρ o .

Displacement controlled Newton-Raphson
In this work, we load the structure quasi-statically by imposing a displacement u =ū on ∂ ou to obtain the deformed configuration . We subsequently evaluate the natural frequencies on . The displacement controlled Newton-Raphson scheme is used to solve the quasi-static residual equilibrium equation obtained by neglecting inertia and equating t o = 0 in Eq. 12, i.e.
where the subscript f refers to the equations associated with the unconstrained degrees of freedom, i.e. the free displacements a f , as opposed to the prescribed displacements a p . To solve (16), we employ the Newton-Raphson scheme where the residual (16) is linearized to obtain the unknown increment da f , in each equilibrium iteration as the solution to where K ff is the partition of the tangent stiffness matrix (18) corresponding to the free degrees of freedom. In Eq. 18, G contains the gradient of the element shape functions, H is a symmetric block matrix that contains the stress, and D is the Voigt notation stiffness matrix that corresponds to D. The explicit format of these matrices appears in e.g. Crisfield (1993).

Eigenfrequency analysis
The conventional eigenvalue problem is formulated for zero initial strain. In contrast to the conventional problem, our eigenvalue problem considers nonzero initial strain wherein the stiffness and eigenfrequencies depend on the displacement field. Our structure is preloaded via the constant displacement a p such that r f = F int,f (a) = 0 where a is the displacement field in the equilibrium configuration. To formulate the eigenvalue problem, we introduce a small time-dependent perturbation δa f from a f so that where δa p = δä p = 0. Insertion of Eqs. 19 into 12 results in A Taylor series expansion of the internal force and ignoring higher-order terms yields which when combined with Eqs. 16 and 18 allows (20) to be expressed as As with the linear problem, we see that δa f is a harmonic the eigenpairs that are obtained from the real, symmetric generalized eigenvalue problem where n denotes the dimension of the problem, i.e. the number of free nodal degrees of freedom in a f , and ω 2 j are the square eigenfrequencies. In our study, the stiffness and mass matrices are symmetric, real, and positive semi-definite, and hence the eigenfrequencies are real and positive, cf. Clough and Penzien (1993). Subsequently, the eigenvalues are placed in ascending order such that 0 < ω 2 1 ≤ ω 2 2 ≤ ... ≤ ω 2 j ≤ ... ≤ ω 2 n , and their associated eigenmodes are mass-normalized such that where δ jk denotes the Kronecker's delta. The formulation of the eigenvalue problem for a nonlinear system (23) is clearly similar to the one obtained in conventional analyzes, cf. e.g. Du and Olhoff (2007); however, it is emphasized that the eigenpairs of the nonlinear system depend on the equilibrium displacement field a since K ff is a function of the displacement a.

Consideration of low volume fraction elements
In our formulation, void regions are modeled using the minute volume fraction δ 0 which is assumed to have a negligible influence on the structural response. Unfortunately, for eigenfrequency constrained finite deformation problems, the existence of low volume fraction elements is problematic.
Firstly, low volume fraction elements yield so-called spurious, localized eigenmodes. The issue arises when solving (23) since the mass scales linearly, i.e. p = 1, with the volume fraction, whereas the stiffness scales at p > 1, cf. Eqs. 10 and 13, which implies that the stiffnessto-mass ratio becomes very small in the void regions where ν(X) → 0. Several ways to suppress these nonphysical eigenmodes have been suggested in the literature. By far the most common method, as proposed by Du and Olhoff (2007), weights the volume fraction such that the stiffnessto-mass ratio, i.e. the spurious eigenfrequency, becomes large.
Secondly, regions consisting of low volume fraction elements often become highly distorted in our finite deformation analysis. This causes convergence issues in the Newton-Raphson equilibrium iterations. Different means for addressing this issue have also been proposed, e.g. by introducing an interpolation of the strain energy such that void regions are modeled using small deformation theory as done by Wang et al. (2014), or by simply relaxing the convergence criterion of the Newton-Raphson iterative scheme as done by Pedersen et al. (2001).
To eliminate the aforementioned problems associated with low volume fraction elements, we employ the element removal method proposed by Bruns and Tortorelli (2003). This procedure is able to both remove and reintroduce the finite elements. Elements are removed from the analysis if their filtered volume fraction is smaller than a prescribed small threshold ε r > 0, i.e. if ν e < ε r ; they are reintroduced if ν e > ε r . The threshold ε r must be chosen such that the discontinuities in the gradient information due to the nonzero tolerance ε r do not greatly affect the optimization. We did not experience any numerical difficulties when using this method. More details regarding the implementation of this method are discussed in Bruns and Tortorelli (2003).
Despite using the element removal scheme, the remaining small regions with intermediate volume fractions ν(x) ∈ (ε r , 1) gave rise to spurious eigenmodes due to the ersatz material model. To combat this, we follow Ferrari et al. (2018) and set the lower bound δ 0 to a relatively large value δ 0 = 10 −3 . This choice of δ 0 increases the magnitude of the spurious eigenfrequencies such that they do not obstruct the optimization. For most ersatz material models, a large δ 0 leads to significant load-carrying capabilities of void regions which might adversely affect the resulting design. However, since we employ element removal, this problem is avoided as the effected regions are small.
By examining the evolution of the potential energy or the kinetic energy of the spurious modeshapes, i.e.
we verify the absence of low energy spurious eigenmodes in the nearly void interface regions. Indeed, such modes would be readily identified by low E pot values.

Topology optimization
The objective of our topology optimization is to distribute material to build a stiff structure with fundamental eigenfrequency greater than ω min . Following Niu et al. (2011), the stiffness of a displacement loaded structure is defined by where F ext,p contains the reaction forces of the degrees of freedom for which the displacements a p are imposed.
Constraints on the n c lowest eigenvalues are cast as and on the maximum volume V max as The optimization problem together with the box constraints on the components of z is formulated as where we reformulate the objective as a minimization problem and consider the three smallest eigenfrequencies, i.e. n c = 3.

Sensitivity analysis
The optimization problem is solved using the gradientbased MMA, where the sensitivities of the objective and constraint functions are computed using the adjoint approach. An example verifying the analytical sensitivities of the eigenfrequencies is found in the Appendix. In the subsequent sensitivity analysis, the nondifferentiable effect of the element removal scheme for elements with ν e = ε r is neglected. Fortunately, for small choices of ε r , we did not experience any convergence issues during the optimization. Since the problem is regularized using a filter, the sensitivities are obtained using the chain rule and Eq. 8 as where we emphasize that Z is the constant filter matrix.

Sensitivity analysis of objective
The sensitivity of the objective function as defined in Eq. 26 has previously been outlined by Luo et al. (2017), but for completeness, it is summarized below. Using the adjoint method, we augment g o and write it equivalently as where we utilize F int,f = 0 and introduce the adjoint vector μ μ μ.
Assigning the adjoint vector μ μ μ such that where

Sensitivity analysis of simple eigenvalues
If eigenfrequency ω j satisfies ω j −1 < ω j < ω j +1 , the eigenpair (ω 2 j , φ φ φ j ) is said to be simple and the corresponding normalized eigenmode φ φ φ j is unique. The sensitivity of a simple eigenvalue with respect to the design variables of a linear elastic structure is trivially obtained as, cf. e.g. Haftka and Gürdal (1992) where as in Eq. 31 λ λ λ j is the adjoint vector and F int,f = 0. Differentiation of Eq. 37 with respect to ν e and some manipulations yield where the mass orthonormalization (24) was utilized and thê ( · ) notation indicates quantities to be treated as constant in the differentiation. To annihilate the implicit sensitivities ∂a f ∂ν e from Eq. 38, the adjoint vector λ λ λ j is chosen as the solution to where the right hand side has been computed by Wallin et al. (2018). Insertion of λ λ λ j obtained from Eq. 39 into 38 yields the sensitivity of the simple eigenvalue ω 2 j with respect to each element volume fraction ν e as ∂ω 2 (40)

Sensitivity analysis of degenerate eigenvalues
Frequency constrained topology optimization invariably produces designs with repeated, i.e. degenerate, eigenfrequencies. Indeed, N-fold degenerate eigenfrequencies of the initial design often exist due to domain symmetry of the body. Additionally, degenerate eigenfrequencies arise during the optimization as multiple simple frequencies ω j that violate the constraint approach the limiting value ω min . In either case, the subsequent sensitivity analysis is affected since the eigenmodes corresponding to the N-fold degenerate eigenfrequency are not unique, i.e. any set of vectors on the N-dimensioned hyperplane will satisfy the original eigenvalue problem (23), and consequently degenerate eigenfrequencies are not differentiable in the common (Fréchet) sense. The phenomena of degenerate eigenfrequencies in structural optimization have previously been studied in detail by e.g. Haftka and Gürdal (1992), using the concept of directional derivatives. Later, Seyranian et al. (1994) calculated the directional sensitivities of degenerate eigenfrequencies using a mathematical perturbation method. This formulation has thereafter been implemented in topology optimization by e.g. Du and Olhoff (2007) and Pedersen and Nielsen (2003). The mathematical perturbation-based sensitivity analysis by Seyranian et al. (1994) is limited to linear elasticity and extended here to incorporate nonlinear hyperelasticity.
By definition, the N-fold degenerate eigenvalue satisfies and its corresponding M-orthonormalized eigenmodes φ φ φ span an N-dimensional hyperplane. To begin the sensitivity analysis, we express an eigenmodeφ φ φ in this hyperplane as The choice β j of Eq. 42 is made to obtain eigenmodes that remain continuous with respect to design changes in the subsequent sensitivity analysis, cf. Courant and Hilbert (1953). As done above, we invoke the adjoint method to eliminate the implicit displacement derivatives ∂a ∂ν e . The augmented version of the original eigenvalue problem (23) for the degenerate case is cast as ; however, the adjoint vectorλ λ λ is now ornamented with ā ( · ) to emphasize the degeneracy, cf. the notation in Eq. 41. Next, we consider the variation of Eq. 43 which yields upon utilizing (23). Rearranging the above (44) results in To annihilate the implicit sensitivities, we require that is fulfilled. In Eq. 46, the definition of the N-fold eigenmodeφ φ φ cf. Eq. 42, is used and in a similar fashion, we now express the adjoint vectorλ λ λ as a function of the to-be-determined vectors λ λ λ jk , i.e.
which allows (46) to be expressed as Thus, the implicit sensitivity in Eq. 45 is annihilated for any choice of β j and β k , j, k = 1, ..., N if the λ λ λ jk satisfy We solve these j, k = 1, ..., N adjoint equations for λ λ λ jk , j, k = 1, ..., N. So in all we have to solve N × N adjoint problems, but due to the symmetry of the right-hand side of Eq. 49, we note that λ λ λ jk = λ λ λ kj and hence we only solve N(N + 1)/2 adjoint problems. Also notable is that the computational cost for solving the N(N + 1)/2 adjoint problems is negligible when using direct solvers since K ff is identical in all k, j -combinations. Upon inserting (49) into (45), we obtain (50) Use of Eqs. 24, 42, and 47 in Eq. 50 subsequently yields The above can be arranged as where β β β = [β 1 , β 2 , ..., β N ] T and the components of the Due to symmetries K T = K, M T = M and λ λ λ jk = λ λ λ kj , we see that S is real and symmetric. To satisfy (52), the matrix S − δω 2 1 must be rank-deficient, i.e.
and thus we see that the eigenvalues of S provide the directional sensitivities δω 2 j , j = 1, ..., N of the degenerate eigenvalueω 2 .
The directional derivatives δω 2 j , j = 1, ..., N of the degenerate eigenvalues are nonlinear functions of the design perturbations δM ff , δK ff , and δF int,f and do not admit a usual linearization in contrast to simple eigenvalues. However, if S is diagonal, the directional derivative δω 2 j can be computed using the simple eigenvalue (37).
Standard gradient-based methods for optimization are not directly applicable for the solution of optimization problems in which degenerate eigenvalues are present due to the lack of differentiability. To combat this, we follow Fig. 1 Double clamped beam. The prescribed displacement u y is uniformly applied to 14 nodes on the center-top surface in the negative y-direction Krog and Olhoff (1999) and Lund (1994), and require S to be diagonal by imposing the constraints on our optimization problem. In this way, we can always use Eq. 40 to evaluate the derivatives ∂ω 2 j ∂ν e of simple eigenvalues and directional derivatives δω 2 j , j = 1, ..., N of degenerate eigenvalues. In the MMA, this constraint is easily accommodated in the subproblem wherein M, K, r and hence S are linear in the update z. As such the g jk sensitivities are trivially evaluated. The final optimization problem P is thus stated as

Numerical examples
To illustrate the influence of the eigenfrequency constraints, we design 2D beam-like structures. To connect with previous related work, the material parameters are based on  Fig. 2 Optimized design a, first three eigenmodes b-d, and evolutions of the eigenvalues and compliance ratio e for the geometry in Fig. 1 when u y = 1000 mm without eigenvalue constraints. Scaled bφ φ φ 1 , c-φ φ φ 2 , and d-φ φ φ 3 are plotted for u y = 0 mm. Every 5th iteration is marked in e, where red-ω 2 1 , blue-ω 2 2 , green-ω 2 3 , and black-g o /g init  Fig. 3 Optimized designs and evolutions of the three smallest eigenvalues and compliance ratio for the geometry in Fig. 1 subject to the eigenvalue constraints ω 2 j > ω 2 min = 250 kHz 2 . a u y ≈ 0 mm, b u y = 200 mm, c u y = 500 mm, d u y = 850 mm, and e u y = 1000 mm. In the right hand side plots, every 5th iteration is marked, where red-ω 2 1 , blue-ω 2 2 , green-ω 2 3 , and black-g o /g init o the choices by Du and Olhoff (2007), i.e. density ρ g = 1.000 kg/m 3 , Young's modulus E = 10.00 MPa, and Poisson's ratio υ = 0.300 which provides the bulk and shear moduli K = E 3(1−2υ) and G = E 2(1+υ) . A continuation scheme is utilized for the penalty factor p, cf. Eq. 10, where the initial value p = 2 is increased by 0.05 every fifth Fig. 4 Optimized Fig. 3 a and e designs in their undeformed configurations design iteration up to the maximum value p = 3. The maximum allowed volume is set to V max = 0.5V where V denotes the volume of the design domain in the reference configuration o . The initial volume fraction distribution is uniform such that ν e = 0.5, e = 1, ..., n e . The threshold ε r for determining when to remove an element is set to ε r = 0.01, i.e. the same choice as Bruns and Tortorelli (2003). The numerical parameters used in the MMA algorithm are those proposed by Svanberg (1987).
Two different boundary conditions are considered, namely a double clamped beam and a simply supported beam. The length, width, and out-of-plane thickness of the beam are 10000 mm, 1000 mm, and 100 mm, respectively. The beam is discretized using 250 × 25 × 1 = 6250 eight node fully integrated 3D brick elements, i.e. one element in thickness. Plane strain is assumed and thus the out-of-plane displacement is zero. The filter radius of Eq. 7 is r = 60 mm, i.e. 1.5 times the element side length.

Double clamped beam
The design domain and boundary conditions for the double clamped 2D beam are shown in Fig. 1.  Fig. 6 Optimized designs a, c and evolutions of the eigenvalues and compliance ratio b, d for the geometry in Fig. 1 subject to the eigenvalue constraints ω 2 j > ω 2 min . a, b ω 2 min = 400 kHz 2 and c, d ω 2 min = 450 kHz 2 . Every 5th iteration is marked in b and d, where red-ω 2 1 , blue-ω 2 2 , green-ω 2 3 , and black-g o /g init o Fig. 7 Scaled first three eigenmodes of the Figs. 3e and 6 designs plotted for u y = 0 mm. a-c Fig. 3e design, d-f Fig. 6a design, and g-i Fig. 6c design depicted in Fig. 2. Note that the φ φ φ j (X) = 0, j = 1, 2, 3 over the 240 mm top center surface over which the displacement is prescribed.
To investigate the influence of the prestrain in conjunction with the eigenfrequency constraints, we provide examples for five load levels u y with ω 2 min = 250 kHz 2 fixed. The optimized designs along with the evolutions of the eigenvalues are shown in Fig. 3, where we also plot the compliance ratio g o /g init o . In Fig. 4, the optimized design corresponding to the ω 2 j > ω 2 min = 250 kHz 2 eigenvalue constraint for the linear u y ≈ 0 mm and nonlinear u y = 1000 mm cases is shown in their undeformed configurations. The difference between the designs demonstrates the impact of the prestrain on the structural response. To further illustrate this effect, we compute the smallest eigenvalue of the Fig. 3e design corresponding to different load levels, cf. Fig. 5. As expected, the structure's stiffness increases under load, causing the lowest eigenvalue to increase until it ultimately equals the 250 kHz 2 constraint limit when u y = 1000 mm. Figure 6 illustrates designs subject to minimum eigenvalue constraints ω 2 min = 400 kHz 2 and ω 2 min = 450 kHz 2 at the load level u y = 1000 mm. The mode shapes associated with the designs in Figs. 3e and 6 are depicted in Fig. 7. From Figs. 2, 3e, and 6, we conclude that the designs are highly influenced by the eigenvalue constraints.
A comparison of the modal shapes in Figs. 2 and 7 indicates that mode switching occurs during the optimization. Our experience is that this phenomena did not systematically affect the convergence of the optimizer, cf. the evolutions of the eigenvalues and the compliance ratio in Figs. 2, 3, and 6.
In Table 1, the values of g o , g init o , ω 2 1 , ω 2 2 , and ω 2 3 for each design are given. An interesting observation is the absence of degenerate eigenfrequencies. The smallest observed relative difference between two eigenfrequencies during all design iterations is ≈ 10 −5 , i.e. the eigenfrequencies are always treated as simple. The last row in Table 1 shows the performance for the Fig. 6a design which is optimized for u y ≈ 0 mm when it is subjected to the u y = 1000 mm displacement. A comparison with the performance of the  The evolution of the number of removed and reintroduced finite elements is depicted for the Fig. 3b design in Fig. 8. We note that the percentage of removed elements reaches a terminal value of ≈ 30%, and that the number of Fig. 9 Simply supported beam. The prescribed displacement u y is uniformly applied to 14 nodes on the center-top surface in the negative y-direction elements that are reintroduced during a design update rarely exceeds 10 elements in a 6250 element mesh.

Simply supported beam
In the last example, we consider the simply supported beam shown in Fig. 9.
Again, we initially apply the u y = 1000 mm load level and omit the eigenvalue constraints which yields the design and mode shapes depicted in Fig. 10. Thereafter, we investigate the influence of the prestrain and use the load u y ≈ 0 mm to mimic a linear design. Imposing the constraint ω 2 j > ω 2 min = 150 kHz 2 renders the design and its mode shapes illustrated in Fig. 11. Figure 12 illustrates designs at the u y = 1000mm prestrain subject to the eigenvalue constraints ω 2 j > ω 2 min with ω 2 min = 150 kHz 2 and ω 2 min = 250 kHz 2 . Similar to the previous example, we observe great influences of the eigenvalue constraints on the final designs. We again emphasize the lack of degenerate eigenfrequencies.
The mode shapes associated with the Fig. 12 designs are plotted in Fig. 13. We observe that a local eigenmode  Fig. 10 Optimized design a, evolutions of the eigenvalues and compliance ratio b, and first three eigenmodes c-e for the geometry in Fig. 9 when u y = 1000 mm without eigenvalue constraints. Scaled c-φ φ φ 1 , d-φ φ φ 2 , and e-φ φ φ 3 are plotted for u y = 0 mm. Every 5th iteration is marked in b where red-ω 2 1 , blue-ω 2 2 , green-ω 2 3 , and black-g o /g init Optimized designs and evolutions of the eigenvalues and compliance ratio for the geometry in Fig. 9 subject to the eigenvalue constraints ω 2 j > ω 2 min . a ω 2 min = 150 kHz 2 and b ω 2 min = 250 kHz 2 . In the right hand side plots, every 5th iteration is marked, where red-ω 2 1 , blue-ω 2 2 , green-ω 2 3 , and black-g o /g init o Fig. 13 First three eigenmodes of the Fig. 12 designs in their undeformed configurations. a-c Fig. 12a and d-g Fig. 12b. f, g The third eigenmode for two different design iterations. First row-φ φ φ 1 , second row-φ φ φ 2 , and third row-φ φ φ 3   Fig. 13f has emerged, whereas for a previous design iteration, the global mode displayed in Fig. 13g was found instead. This phenomena can be identified by the small drops in the evolution of the third eigenvalue in Fig. 12. We compare the optimized designs corresponding to u y ≈ 0 mm and u y = 1000 mm subject to the ω 2 j > ω 2 min = 150 kHz 2 eigenvalue constraint in the undeformed configurations in Fig. 14. The impact of the prestrain on the optimized topology is clear.
In Table 2, the values of the response quantities of interest are listed for the various designs. It is concluded that the optimization problem P succeeds in drastically increasing the fundamental eigenfrequency. Again, we evaluate the impact of the nonlinear framework by analyzing the performance of the Fig. 11 design obtained for u y ≈ 0 mm at a u y = 1000 mm load level, cf. the last row in Table 2. As expected, the Fig. 12a design outperforms its linear counterpart at this load level.

Conclusions
The mathematical framework for performing topology optimization considering eigenfrequencies is extended to finite strain hyperelastic structures subjected to small harmonic oscillations about a prestrained equilibrium configuration. A stiff structure is obtained, while constraining the lowest eigenfrequencies and the maximum allowable volume. The low volume fraction element removal method proposed by Bruns and Tortorelli (2003) is implemented to mitigate problems associated with nearly void regions in finite deformation eigenfrequency topology optimization, i.e. gross distortion and spurious eigenmodes. The problem in the sensitivity analysis associated with the nondifferentiability of degenerate eigenfrequencies is treated by extending the perturbation method of Seyranian et al. (1994) to include both geometrical and material nonlinearities.
The numerical examples show that the fundamental eigenfrequency of prestrained structures can be increased significantly using the proposed framework. The sensitivity analysis implementation works for both simple and degenerate eigenfrequencies. As expected, the magnitude of the prestrain significantly affects the optimized designs.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.
Based on the analysis above, we conclude that N > 1 eigenvalue is N-fold degenerate when they agree to within the tolerance of ε tol = 10 −8 . This choice is smaller than the tolerance used by e.g. Seyranian et al. (1994). Surely, this is a case-by-case choice. For example, here we are dealing with nonlinear boundary value problem, and thus structural equilibrium is only ensured to within a given tolerance. Another factor affecting the value of ε tol is the eigensolver; We use the FEAST Eigenvalue Solver 2.0, cf. Polizzi (2009).
The problem considered for verification of the sensitivity computation is a simply supported 3D plate-like structure cf. the design domain in Fig. 15. The length, width, and thickness of the plate are, 2000 mm, 2000 mm, and 100 mm, respectively. We use a coarse mesh that consists of 20 × 20 × 1 = 400 fully integrated 3D brick elements. The filter radius, cf. Eq. 7, is r = 150 mm, i.e. 1.5 times the element side length.
This geometry is chosen due to its inherent domain symmetry which is ideal for verifying the sensitivities of degenerate eigenfrequencies. First, we obtain an optimized design without constraining the eigenvalue. The resulting design and associated mode shapes appear in Fig. 16, the evolutions of the eigenvalue and the compliance ratio are plotted in Fig. 17, and the compliance and the total volume fraction 0 ν dv /V = V * /V are denoted in Table 3. Next, we repeat the previous problem, but now we enforce the eigenvalue constraints ω 2 j > ω 2 min = 900 kHz 2 . The resulting design is depicted in Fig. 18, where the mode shapes also are seen. As is shown in Fig. 19, the twofold degenerate fundamental eigenvalue is increased to fulfill the lower bound ω 2 min = 900 kHz 2 constraint. Since a coarse mesh is used to discretize the domain, the optimization algorithm is working in a relatively small design space, which might explain such small differences between the Figs. 16 and 18 designs. An interesting phenomenon is noted as the Fig. 18 design does not utilize all available material, cf. Fig. 20. Similar observations are made in Stolpe (2010), who discusses the nonnecessity of an active volume constraint at an optimal solution for continuous and discrete design problems. The terminal values of the constrained eigenvalues and g o are denoted in Table 3.