Eigenfrequency constrained topology optimization of finite strain hyperelastic structures

Abstract

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.

1 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.

2 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 ${{\Omega }_{o}^{e}}=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 + ∇_ou, 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 (\boldsymbol {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^TF and the Green-Lagrange strain tensor $\boldsymbol {E} = \frac {1}{2}(\boldsymbol {C} - \boldsymbol {1})$.

In absence of body forces, the equation of motion is

$$ \boldsymbol{\nabla}_{o} \cdot (\boldsymbol{F}\boldsymbol{S}) = \rho^{o} \ddot{\boldsymbol{u}} \quad \text{in} \quad {\Omega}_{o}, $$

(1)

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 = JF^− 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 = FSn^o, and the displacement $\boldsymbol {u} = \boldsymbol {\varphi } - \boldsymbol {X} = \bar {\boldsymbol {u}}$, 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 $\boldsymbol {u} = \bar {\boldsymbol {u}}$ on ∂Ω_ou, such that

$$ {\int}_{{\Omega}_{o}} \rho^{o} \delta\boldsymbol{u}\cdot \ddot{\boldsymbol{u}} dv + {\int}_{{\Omega}_{o}}\delta\boldsymbol{E} : \boldsymbol{S} dv - {\int}_{\partial{\Omega}_{ot}}\delta\boldsymbol{u}\cdot \boldsymbol{t}^{o} ds = 0, $$

(2)

for all smooth virtual displacements that satisfy δu = 0 on ∂Ω_ou. The virtual Lagrangian strain δE, introduced in Eq. 2, is defined as

$$ \delta\boldsymbol{E}(\boldsymbol{\varphi},\delta\boldsymbol{u}) = \frac{1}{2}\left( (\boldsymbol{\nabla}_{o} \delta\boldsymbol{u})^{T} \boldsymbol{F} + \boldsymbol{F}^{T}\boldsymbol{\nabla}_{o} \delta\boldsymbol{u} \right). $$

(3)

3 Constitutive model

Our constitutive assumption uses the isotropic compressible neo-Hookean strain energy model in which the strain energy U^NH is expressed as

$$ U^{NH} = \frac{1}{4}K \left( (J^{2} - 1) - \ln (J^{2}) \right) + \frac{1}{2}G(J^{-2/3}\text{tr}(\boldsymbol{C}) -3), $$

(4)

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

$$ \boldsymbol{S} = 2\frac{\partial U^{NH}}{\partial \boldsymbol{C}} = \frac{K}{2}(J^{2}-1) \boldsymbol{C}^{-1} + GJ^{-2/3}\left( \textsf{\textbf{1}} - \frac{\text{tr}(\boldsymbol{C})}{3}\boldsymbol{C}^{-1}\right). $$

(5)

The fourth-order stiffness tensor associated with Eq. 4 becomes

$$ \begin{array}{@{}rcl@{}} \boldsymbol{D} = 2\frac{\partial \boldsymbol{S}}{\partial \boldsymbol{C}} &= &a_{1} \boldsymbol{C}^{-1} \otimes \boldsymbol{C}^{-1} + a_{2}\left( \boldsymbol{1} \otimes \boldsymbol{C}^{-1} + \boldsymbol{C}^{-1} \otimes \boldsymbol{1}\right)\\ &&+ a_{3}\left( \boldsymbol{C}^{-1}\underline{\otimes}\boldsymbol{C}^{-1} + \boldsymbol{C}^{-1}\overline{\otimes}\boldsymbol{C}^{-1}\right), \end{array} $$

(6)

where $a_{1} = KJ^{2} + \frac {2G}{9}J^{-2/3}\text {tr}(\boldsymbol {C})$, $a_{2} = -\frac {2G}{3}J^{-2/3}$ and $a_{3} = \frac {G}{3}J^{-2/3}\text {tr}(\boldsymbol {C}) - \frac {K}{2}(J^{2} -1)$, respectively. The dyadic products $\underline {\otimes }$ and $\overline {\otimes }$ are defined such that $(\boldsymbol {A}\underline {\otimes }\boldsymbol {T}) : \boldsymbol {H} = \boldsymbol {A}\cdot \boldsymbol {H}^{T} \cdot \boldsymbol {T}^{T}$ and $(\boldsymbol {A}\overline {\otimes }\boldsymbol {T}) : \boldsymbol {H} = \boldsymbol {A}\cdot \boldsymbol {H} \cdot \boldsymbol {T}^{T}$, where A, T and H are arbitrary second-order tensors.

4 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 \rightarrow 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¹. To obtain a well-posed optimization problem, we impose a length-scale 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

$$ \nu(\boldsymbol{T}) = \frac{{\int}_{{\Omega}_{o}} \omega(\boldsymbol{X}- \text{\boldmath \textit{Y}}) z(\text{\boldmath \textit{Y}}) dv}{{\int}_{{\Omega}_{o}} \omega(\boldsymbol{X}- \text{\boldmath \textit{Y}}) dv} \approx \frac{\sum\limits_{s=1}^{n_{e}} \omega(\textsf{\textbf{X}}_{e} - \textsf{\textbf{Y}}_{s}) \textsf{z}_{s}}{\sum\limits_{s=1}^{n_{e}} \omega(\textsf{\textbf{X}}_{e}- \textsf{\textbf{Y}}_{s})}, $$

(7)

where $\boldsymbol {X} \in {{\Omega }_{o}^{e}}$. 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 ${{\Omega }_{o}^{e}}$ centroid position X_e and enveloping the surrounding element ${{\Omega }_{o}^{s}}$ 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

$$ \pmb{\upnu} = \textsf{\textbf{Z}}\ \textsf{\textbf{z}}, $$

(8)

where Z is the [n_e × n_e] filter matrix, $\textsf {\textbf {z}} = [\textsf {z}_{1},\textsf {z}_{2},..., \textsf {z}_{n_{e}}]^{T}$ and $\pmb {\upnu } = [{\upnu }_{1}, {\upnu }_{2}, ..., {\upnu }_{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

$$ U = \chi(\nu) U^{NH}, $$

(9)

where

$$ \chi(\nu) = \nu^{p} (1 -\delta_{0}) + \delta_{0}. $$

(10)

Here δ₀ 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 $\chi (\nu )\rightarrow \delta _{0}$ the material mimics void, whereas as $\chi (\nu )\rightarrow 1$ the material realizes the neo-Hookean strain energy function (4).

5 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 $\ddot {\textsf {\textbf {u}}}(\boldsymbol {X}) = \textsf {\textbf {N}}(\boldsymbol {X})\ddot {\textsf {\textbf {a}}}_{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

$$ \delta\textsf{\textbf{E}}(\textsf{\textbf{a}}_{e},\delta\textsf{\textbf{a}}_{e}) = \textsf{\textbf{B}}(\textsf{\textbf{a}}_{e})\delta\textsf{\textbf{a}}_{e} = (\textsf{\textbf{B}}^{c} + \textsf{\textbf{B}}^{l}(\textsf{\textbf{a}}_{e}))\delta\textsf{\textbf{a}}_{e}, $$

(11)

where B(a_e) has been decomposed into constant B^c, and linear $\textsf {\textbf {B}}^{l}(\textsf {\textbf {a}}_{e})$ terms, cf. Crisfield (1993) for details ². Using the arbitrariness of the virtual nodal displacement and the finite element interpolations in Eq. 2 yields the finite element equation of motion

$$ \textsf{\textbf{M}} \ddot{\textsf{\textbf{a}}} + \textsf{\textbf{F}}_{int}(\textsf{\textbf{a}}) - \textsf{\textbf{F}}_{ext} = \boldsymbol{0}, $$

(12)

where the mass matrix M, internal force vector F_int, and external force vector F_ext, are

(13)

(14)

(15)

In Eq. 13, we emphasize that ρ = νρ^o.

5.1 Displacement controlled Newton-Raphson

In this work, we load the structure quasi-statically by imposing a displacement $\boldsymbol {u} = \bar {\boldsymbol {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.

$$ \textsf{\textbf{r}}_{f}(\textsf{\textbf{a}}) = \textsf{\textbf{F}}_{int,f}(\textsf{\textbf{a}}) = \boldsymbol{0}, $$

(16)

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

$$ \textsf{\textbf{r}}_{f}(\textsf{\textbf{a}}) + \frac{\partial \textsf{\textbf{r}}_{f}(\textsf{\textbf{a}})}{\partial \textsf{\textbf{a}}_{f}}d\textsf{\textbf{a}}_{f} = \boldsymbol{0} \quad \Rightarrow \quad \textsf{\textbf{K}}_{ff}(\textsf{\textbf{a}})d\textsf{\textbf{a}}_{f} = -\textsf{\textbf{r}}_{f}(\textsf{\textbf{a}}), $$

(17)

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).

6 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

$$ \textsf{\textbf{a}} \rightarrow \textsf{\textbf{a}} + \delta\textsf{\textbf{a}} \ \Rightarrow \ \ddot{\textsf{\textbf{a}}} \rightarrow \ddot{\textsf{\textbf{a}}} + \delta\ddot{\textsf{\textbf{a}}} = \delta\ddot{\textsf{\textbf{a}}}, $$

(19)

where $\delta \textsf {\textbf {a}}_{p} = \delta \ddot {\textsf {\textbf {a}}}_{p} = \boldsymbol {0}$. Insertion of Eqs. 19 into 12 results in

$$ \textsf{\textbf{M}}\delta\ddot{\textsf{\textbf{a}}} + \textsf{\textbf{F}}_{int}(\textsf{\textbf{a}} + \delta\textsf{\textbf{a}}) = \boldsymbol{0} \ \Rightarrow \ \textsf{\textbf{M}}_{ff}\delta\ddot{\textsf{\textbf{a}}}_{f} + \textsf{\textbf{F}}_{int,f}(\textsf{\textbf{a}} + \delta\textsf{\textbf{a}}) = \boldsymbol{0}. $$

(20)

A Taylor series expansion of the internal force and ignoring higher-order terms yields

$$ \textsf{\textbf{F}}_{int,f}(\textsf{\textbf{a}} + \delta \textsf{\textbf{a}}) = \textsf{\textbf{F}}_{int,f}(\textsf{\textbf{a}}) + \frac{\partial \textsf{\textbf{F}}_{int,f}(\textsf{\textbf{a}})}{\partial \textsf{\textbf{a}}_{f}}\delta\textsf{\textbf{a}}_{f}, $$

(21)

which when combined with Eqs. 16 and 18 allows (20) to be expressed as

$$ \textsf{\textbf{M}}_{ff}\delta\ddot{\textsf{\textbf{a}}}_{f} + \textsf{\textbf{K}}_{ff} \delta\textsf{\textbf{a}}_{f} = \boldsymbol{0}. $$

(22)

As with the linear problem, we see that δa_f is a harmonic oscillation, i.e. $\delta \textsf {\textbf {a}}_{f} = \sum \limits _{j=1}^{n} \pmb {\upphi }_{j} e^{i{\upomega }_{j} t}$ where $({{\upomega }_{j}^{2}}, \pmb {\upphi }_{j})$ are the eigenpairs that are obtained from the real, symmetric generalized eigenvalue problem

$$ \textsf{\textbf{K}}_{ff} \pmb{\upphi}_{j} = {{\upomega}_{j}^{2}}\textsf{\textbf{M}}_{ff}\pmb{\upphi}_{j} , \quad j = 1,...,n, $$

(23)

where n denotes the dimension of the problem, i.e. the number of free nodal degrees of freedom in a_f, and ${{\upomega }^{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 < {{\upomega }^{2}_{1}} \leq {{\upomega }^{2}_{2}} \leq ... \leq {{\upomega }^{2}_{j}} \leq ... \leq {{\upomega }^{2}_{n}}$, and their associated eigenmodes are mass-normalized such that

$$ \pmb{\upphi}_{j}^{T}\textsf{\textbf{M}}_{ff}\pmb{\upphi}_{k} = \delta_{jk} , \quad j,k = 1,...,n, $$

(24)

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.

7 Consideration of low volume fraction elements

In our formulation, void regions are modeled using the minute volume fraction δ₀ 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 stiffness-to-mass ratio becomes very small in the void regions where $\nu (\boldsymbol {X}) \rightarrow 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 stiffness-to-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 δ₀ to a relatively large value δ₀ = 10^− 3. This choice of δ₀ increases the magnitude of the spurious eigenfrequencies such that they do not obstruct the optimization. For most ersatz material models, a large δ₀ 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.

$$ \begin{array}{@{}rcl@{}} &&E_{pot}= \pmb{\upphi}_{j}^{T}\textsf{\textbf{K}}_{ff}\pmb{\upphi}_{j} \quad = \quad E_{kin} = {{\upomega}_{j}^{2}}\pmb{\upphi}_{j}^{T}\textsf{\textbf{M}}_{ff}\pmb{\upphi}_{j}, \\ &&j=1,...,n, \end{array} $$

(25)

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.

8 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

$$ g_{o} = \textsf{\textbf{F}}_{ext,p}^{T} \textsf{\textbf{a}}_{p}, $$

(26)

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

$$ g_{j} = {\upomega}^{2}_{min} - {{\upomega}_{j}^{2}} \leq 0, \quad j=1,...,n_{c}, $$

(27)

and on the maximum volume $V_{\max \limits }$ as

$$ g_{V} = {\int}_{{\Omega}_{0}} \nu dv - V_{\max} \leq 0. $$

(28)

The optimization problem together with the box constraints on the components of z is formulated as

$$ \mathcal{P} \left\{\begin{array}{ll} \underset{\textsf{\textbf{z}}}{\min } \ -g_{o} \\ \text{s.t} \quad \left\{\begin{array}{ll} g_{j} \leq 0, \quad j=1,...,n_{c} \\ g_{V} \leq 0\\ 0 \leq \textsf{z}_{e} \leq 1, \quad e = 1,...,n_{e} \ \end{array}\right. \end{array}\right. $$

(29)

where we reformulate the objective as a minimization problem and consider the three smallest eigenfrequencies, i.e. n_c = 3.

9 Sensitivity analysis

The optimization problem is solved using the gradient-based 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

$$ \frac{\partial g_{i}}{\partial \textsf{\textbf{z}}} = \left( \frac{\partial\boldsymbol{\upnu}}{\partial \textsf{\textbf{z}}}\right)^{T}\frac{\partial g_{i}}{\partial \boldsymbol{\upnu}} = \textsf{\textbf{Z}}^{T} \frac{\partial g_{i}}{\partial\boldsymbol{\upnu}}, $$

(30)

where we emphasize that Z is the constant filter matrix.

9.1 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

$$ g_{o} = \textsf{\textbf{a}}_{p}^{T} \textsf{\textbf{F}}_{int,p} + \pmb{\upmu}^{T} \textsf{\textbf{F}}_{int,f}, $$

(31)

where we utilize F_int,f = 0 and introduce the adjoint vector μ. Differentiation of Eq. 31 with respect to the filtered volume fraction variables yields

$$ \begin{array}{@{}rcl@{}} \frac{\partial g_{o}}{\partial {\upnu}_{e}} &= &\textsf{\textbf{a}}_{p}^{T} \frac{\partial \textsf{\textbf{F}}_{int,p}}{\partial {\upnu}_{e}} +\pmb{\upmu}^{T} \frac{\partial \textsf{\textbf{F}}_{int,f}}{\partial {\upnu}_{e}} \ + \\ &&+\left( \textsf{\textbf{a}}_{p}^{T} \frac{\partial \textsf{\textbf{F}}_{int,p}}{\partial \textsf{\textbf{a}}_{f}} + \pmb{\upmu}^{T}\frac{\partial \textsf{\textbf{F}}_{int,f}}{\partial \textsf{\textbf{a}}_{f}} \right)\frac{\partial \textsf{\textbf{a}}_{f}}{\partial {\upnu}_{e}}. \end{array} $$

(32)

Assigning the adjoint vector μ such that

$$ \textsf{\textbf{K}}_{ff} \pmb{\upmu} = -\textsf{\textbf{K}}_{fp} \textsf{\textbf{a}}_{p}, $$

(33)

where $\textsf {\textbf {K}}_{fp} = \frac {\partial \textsf {\textbf {F}}_{int,f}}{\partial \textsf {\textbf {a}}_{p}}$ annihilates the $\frac {\partial \textsf {\textbf {a}}_{f}}{\partial {\upnu }_{e}}$ contribution in Eq. 32. The sensitivity is thus reduced to

$$ \frac{\partial g_{o}}{\partial {\upnu}_{e}} = \textsf{\textbf{a}}_{p}^{T}\frac{\partial \textsf{\textbf{F}}_{int,p}}{\partial {\upnu}_{e}} + \pmb{\upmu}^{T} \frac{\partial \textsf{\textbf{F}}_{int,f}}{\partial {\upnu}_{e}} =\pmb{\upalpha}^{T} \frac{\partial \textsf{\textbf{F}}_{int} }{\partial {\upnu}_{e}} , $$

(34)

where

$$ \pmb{\upalpha} = \left[\begin{array}{cc} \pmb{\upmu} \\ \textsf{\textbf{a}}_{p} \end{array}\right]. $$

(35)

9.2 Sensitivity analysis of simple eigenvalues

If eigenfrequency ω_j satisfies ω_j− 1 < ω_j < ω_j+ 1, the eigenpair (${{\upomega }_{j}^{2}},\pmb {\upphi }_{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)

$$ \frac{\partial {{\upomega}_{j}^{2}}}{\partial {\upnu}_{e}} = \pmb{\upphi}_{j}^{T} \left( \frac{\partial \textsf{\textbf{K}}_{ff}}{\partial {\upnu}_{e}} - {{\upomega}_{j}^{2}} \frac{\partial\textsf{\textbf{M}}_{ff}}{\partial {\upnu}_{e}} \right)\pmb{\upphi}_{j} , \quad j = 1,...,n_{c}. $$

(36)

Unfortunately, the corresponding sensitivity analysis for a nonlinear elastic body is less straight forward since the stiffness matrix depends on the displacement, i.e. $\frac {\partial \textsf {\textbf {K}}}{\partial {\upnu }_{e}} \rightarrow \frac {\partial \textsf {\textbf {K}}}{\partial {\upnu }_{e}} + \frac {\partial \textsf {\textbf {K}}}{\partial \textsf {\textbf {a}}}\frac {\partial \textsf {\textbf {a}}}{\partial {\upnu }_{e}}$. To eliminate the implicitly defined sensitivity $\frac {\partial \textsf {\textbf {a}}}{\partial {\upnu }_{e}}$, we utilize the adjoint approach, whereby we manipulate (16) and (23) to obtain

$$ \pmb{\upphi}^{T}_{j}(\textsf{\textbf{K}}_{ff} - {{\upomega}^{2}_{j}}\textsf{\textbf{M}}_{ff})\pmb{\upphi}_{j} - \pmb{\uplambda}_{j}^{T}\textsf{\textbf{F}}_{int,f}= \boldsymbol{0}, \quad j = 1,...,n_{c}, $$

(37)

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

$$ \begin{array}{@{}rcl@{}} \frac{\partial {{\upomega}_{j}^{2}}}{\partial {\upnu}_{e}} &= &\pmb{\upphi}_{j}^{T}\left( \frac{\partial \textsf{\textbf{K}}_{ff}}{\partial {\upnu}_{e}} -{{\upomega}_{j}^{2}} \frac{\partial \textsf{\textbf{M}}_{ff}}{\partial {\upnu}_{e}} \right)\pmb{\upphi}_{j} \ + \\ &&+ \left( \frac{\partial (\hat{\pmb{\upphi}}_{j}^{T} \textsf{\textbf{K}}_{ff} \hat{\pmb{\upphi}}_{j})}{\partial \textsf{\textbf{a}}_{f}} - \pmb{\uplambda}_{j}^{T} \textsf{\textbf{K}}_{ff} \right)\frac{\partial \textsf{\textbf{a}}_{f}}{\partial {\upnu}_{e}} - \pmb{\uplambda}_{j}^{T}\frac{\partial \textsf{\textbf{F}}_{int,f}}{\partial {\upnu}_{e}}, \end{array} $$

(38)

where the mass orthonormalization (24) was utilized and the $\hat {(\cdot )}$ notation indicates quantities to be treated as constant in the differentiation. To annihilate the implicit sensitivities $\frac {\partial \textsf {\textbf {a}}_{f}}{\partial {\upnu }_{e}}$ from Eq. 38, the adjoint vector λ_j is chosen as the solution to

$$ \textsf{\textbf{K}}_{ff}\pmb{\uplambda}_{j} = \left( \frac{\partial (\hat{\pmb{\upphi}}_{j}^{T} \textsf{\textbf{K}}_{ff} \hat{\pmb{\upphi}}_{j})}{\partial \textsf{\textbf{a}}_{f}}\right)^{T}, $$

(39)

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 ${{\upomega }^{2}_{j}}$ with respect to each element volume fraction ν_e as

$$ \frac{\partial {{\upomega}_{j}^{2}}}{\partial {\upnu}_{e}} = \pmb{\upphi}_{j}^{T}\left( \frac{\partial \textsf{\textbf{K}}_{ff}}{\partial {\upnu}_{e}} -{{\upomega}_{j}^{2}} \frac{\partial \textsf{\textbf{M}}_{ff}}{\partial {\upnu}_{e}} \right)\pmb{\upphi}_{j} - \pmb{\uplambda}_{j}^{T}\frac{\partial \textsf{\textbf{F}}_{int,f} }{\partial {\upnu}_{e}}. $$

(40)

9.3 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 ${\upomega }_{\min \limits }$. 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

$$ \bar{\upomega}^{2}={{\upomega}_{j}^{2}} , \quad j = 1,..,N , $$

(41)

and its corresponding M-orthonormalized eigenmodes ϕ span an N-dimensional hyperplane. To begin the sensitivity analysis, we express an eigenmode $\bar {\pmb {\upphi }}$ in this hyperplane as

$$ \bar{\pmb{\upphi}} = \sum\limits_{j=1}^{N} {\upbeta}_{j}\pmb{\upphi}_{j}. $$

(42)

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 $\frac {\partial \textsf {\textbf {a}}}{\partial {\upnu }_{e}}$. The augmented version of the original eigenvalue problem (23) for the degenerate case is cast as

$$ \bar{\pmb{\upphi}}^{T}(\textsf{\textbf{K}}_{ff} - \bar{\upomega}^{2}\textsf{\textbf{M}}_{ff})\bar{\pmb{\upphi}} - \bar{\pmb{\uplambda}}^{T}\textsf{\textbf{F}}_{int,f}= 0 $$

(43)

; however, the adjoint vector $\bar {\pmb {\uplambda }}$ is now ornamented with a $\bar { (\cdot ) }$ to emphasize the degeneracy, cf. the notation in Eq. 41.

Next, we consider the variation of Eq. 43 which yields

$$ \begin{array}{@{}rcl@{}} &&\bar{\pmb{\upphi}}^{T}\left( \delta\textsf{\textbf{K}}_{ff} - \delta\bar{\upomega}^{2}\textsf{\textbf{M}}_{ff} - \bar{\upomega}^{2}\delta\textsf{\textbf{M}}_{ff}\right)\bar{\pmb{\upphi}} - \bar{\pmb{\uplambda}}^{T}\delta\textsf{\textbf{F}}_{int,f}\\ &&+ \frac{\partial \left( \hat{\bar{\pmb{\upphi}}}^{T} \textsf{\textbf{K}}_{ff} \hat{\bar{\pmb{\upphi}}}\right)}{\partial \textsf{\textbf{a}}_{f}}\delta\textsf{\textbf{a}}_{f} - \bar{\pmb{\uplambda}}^{T}\frac{\partial \textsf{\textbf{F}}_{int,f}}{\partial \textsf{\textbf{a}}_{f}}\delta\textsf{\textbf{a}}_{f} = 0, \end{array} $$

(44)

upon utilizing (23). Rearranging the above (44) results in

$$ \begin{array}{@{}rcl@{}} &&\bar{\pmb{\upphi}}^{T}\left( \delta\textsf{\textbf{K}}_{ff} - \bar{\upomega}^{2}\delta\textsf{\textbf{M}}_{ff}\right)\bar{\pmb{\upphi}} - \delta\bar{\upomega}^{2}\bar{\pmb{\upphi}}^{T}\textsf{\textbf{M}}_{ff}\bar{\pmb{\upphi}} - \bar{\pmb{\uplambda}}^{T}\delta\textsf{\textbf{F}}_{int,f}\\ && + \left( \frac{\partial\left( \hat{\bar{\pmb{\upphi}}}^{T} \textsf{\textbf{K}}_{ff} \hat{\bar{\pmb{\upphi}}}\right)}{\partial \textsf{\textbf{a}}_{f}} - \bar{\pmb{\uplambda}}^{T}\textsf{\textbf{K}}_{ff} \right)\delta\textsf{\textbf{a}}_{f} = 0. \end{array} $$

(45)

To annihilate the implicit sensitivities, we require that

$$ \begin{array}{@{}rcl@{}} \textsf{\textbf{K}}_{ff}\bar{\pmb{\uplambda}} &=& \left( \frac{\partial\left( \hat{\bar{\pmb{\upphi}}}^{T} \textsf{\textbf{K}}_{ff} \hat{\bar{\pmb{\upphi}}}\right)}{\partial \textsf{\textbf{a}}_{f}}\right)^{T}\\ & =& \sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N} {\upbeta}_{j}{\upbeta}_{k}\left( \frac{\partial\left( \hat{\pmb{\upphi}}^{T}_{j} \textsf{\textbf{K}}_{ff} \hat{\pmb{\upphi}}_{k}\right)}{\partial \textsf{\textbf{a}}_{f}}\right)^{T}, \end{array} $$

(46)

is fulfilled. In Eq. 46, the definition of the N-fold eigenmode $\bar {\pmb {\upphi }}$ cf. Eq. 42, is used and in a similar fashion, we now express the adjoint vector $\bar {\pmb {\uplambda }}$ as a function of the to-be-determined vectors λ_jk, i.e.

$$ \bar{\pmb{\uplambda}}= \sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N} {\upbeta}_{j}{\upbeta}_{k}\pmb{\uplambda}_{jk} , $$

(47)

which allows (46) to be expressed as

$$ \sum\limits_{j=1}^{N} \sum\limits_{k=1}^{N} {\upbeta}_{j}{\upbeta}_{k}\textsf{\textbf{K}}_{ff}\pmb{\uplambda}_{jk} = \sum\limits_{j=1}^{N} \sum\limits_{k=1}^{N} {\upbeta}_{j}{\upbeta}_{k}\left( \frac{\partial\left( \hat{\pmb{\upphi}}^{T}_{j} \textsf{\textbf{K}}_{ff} \hat{\pmb{\upphi}}_{k}\right)}{\partial \textsf{\textbf{a}}_{f}}\right)^{T}. $$

(48)

Thus, the implicit sensitivity in Eq. 45 is annihilated for any choice of β_j and β_k, j,k = 1,...,N if the λ_jk satisfy

$$ \textsf{\textbf{K}}_{ff}\pmb{\uplambda}_{jk}= \left( \frac{\partial\left( \hat{\pmb{\upphi}}^{T}_{j} \textsf{\textbf{K}}_{ff} \hat{\pmb{\upphi}}_{k}\right)}{\partial \textsf{\textbf{a}}_{f}}\right)^{T}. $$

(49)

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

$$ \bar{\pmb{\upphi}}^{T}\left( \delta\textsf{\textbf{K}}_{ff} - \bar{\upomega}^{2}\delta\textsf{\textbf{M}}_{ff}\right)\bar{\pmb{\upphi}} - \delta\bar{\upomega}^{2}\bar{\pmb{\upphi}}^{T}\textsf{\textbf{M}}_{ff}\bar{\pmb{\upphi}} - \bar{\pmb{\uplambda}}^{T}\delta\textsf{\textbf{F}}_{int,f} = 0. $$

(50)

Use of Eqs. 24, 42, and 47 in Eq. 50 subsequently yields

$$ \begin{array}{@{}rcl@{}} & & \sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N} {\upbeta}_{j} \left[\pmb{\upphi}_{j}^{T} \left( \delta\textsf{\textbf{K}}_{ff} - \bar{\upomega}^{2}\delta\textsf{\textbf{M}}_{ff}\right) \pmb{\upphi}_{k} - \pmb{\uplambda}_{jk}^{T}\delta \textsf{\textbf{F}}_{int,f} \ + \right.\\ &&\left.\qquad\qquad \quad- \delta\bar{\upomega}^{2} \delta_{jk} \right]{\upbeta}_{k} = 0. \end{array} $$

(51)

The above can be arranged as

$$ \pmb{\upbeta}^{T} \left[ \textsf{\textbf{S}} - \delta\bar{\upomega}^{2} \textsf{\textbf{1}} \right]\pmb{\upbeta} = 0, $$

(52)

where β = [β₁,β₂,...,β_N]^T and the components of the N × N matrix S are

$$ \textsf{S}_{jk} = \pmb{\upphi}_{j}^{T} \left( \delta\textsf{\textbf{K}}_{ff} - \bar{\upomega}^{2}\delta\textsf{\textbf{M}}_{ff}\right) \pmb{\upphi}_{k} - \pmb{\uplambda}_{jk}^{T}\delta \textsf{\textbf{F}}_{int,f}. $$

(53)

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 $\textsf {\textbf {S}} - \delta \bar {\upomega }^{2} \boldsymbol {1}$ must be rank-deficient, i.e.

$$ \det \left[ \textsf{\textbf{S}} - \delta\bar{\upomega}^{2} \textsf{\textbf{1}} \right] = 0 , $$

(54)

and thus we see that the eigenvalues of S provide the directional sensitivities $\delta \bar {\upomega }_{j}^{2}$, j = 1,...,N of the degenerate eigenvalue $\bar {\upomega }^{2}$.

The directional derivatives $\delta \bar {\upomega }^{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 $\delta \bar {\upomega }_{j}^{2}$ 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 Krog and Olhoff (1999) and Lund (1994), and require S to be diagonal by imposing the constraints

$$ g_{jk} = \textsf{S}_{jk} = 0 ,\quad j\neq k, \ j,k = 1,...,N, $$

(55)

on our optimization problem. In this way, we can always use Eq. 40 to evaluate the derivatives $\frac {{\partial {\upomega }_{j}^{2}}}{\partial {\upnu }_{e}}$ of simple eigenvalues and directional derivatives $\delta \bar {\upomega }_{j}^{2}$, 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 $\mathcal {P}$ is thus stated as

$$ \mathcal{P} \left\{\begin{array}{ll} \underset{\textsf{\textbf{z}}}{\min } \ -g_{o} \\ \text{s.t} \quad \left\{\begin{array}{ll} g_{j} \leq 0, \quad j=1,...,n_{c} \\ g_{jk} = 0, \quad j\neq k, \ j,k = 1,...,N \\ g_{V} \leq 0\\ 0 \leq \textsf{z}_{e} \leq 1, \quad e = 1,...,n_{e} \ \end{array}\right. \end{array}\right. $$

(56)

10 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 the choices by Du and Olhoff (2007), i.e. density ρ_g = 1.000 kg/m³, Young’s modulus E = 10.00 MPa, and Poisson’s ratio υ = 0.300 which provides the bulk and shear moduli $K = \frac {E}{3(1-2\upsilon )}$ and $G = \frac {E}{2(1+\upsilon )}$. 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 design iteration up to the maximum value p = 3. The maximum allowed volume is set to $V_{\max \limits } = 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.

10.1 Double clamped beam

The design domain and boundary conditions for the double clamped 2D beam are shown in Fig. 1.

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

As reference, we first investigate the design, structural response, and mode shapes for a u_y = 1000 mm load level when omitting the eigenvalue constraints, with the results 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.

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—ϕ₁, c—ϕ₂, and d—ϕ₃ are plotted for u_y = 0 mm. Every 5th iteration is marked in e, where red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, green—${{\upomega }^{2}_{3}}$, and black—$g_{o}/g_{o}^{\text {init}}$

To investigate the influence of the prestrain in conjunction with the eigenfrequency constraints, we provide examples for five load levels u_y with ${\upomega }^{2}_{\min \limits } = 250$ kHz² 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_{o}^{\text {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 ${{\upomega }_{j}^{2}} >{\upomega }^{2}_{\min \limits } = 250$ kHz². au_y ≈ 0 mm, bu_y = 200 mm, cu_y = 500 mm, du_y = 850 mm, and eu_y = 1000 mm. In the right hand side plots, every 5th iteration is marked, where red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, green—${{\upomega }^{2}_{3}}$, and black—$g_{o}/g_{o}^{\text {init}}$

In Fig. 4, the optimized design corresponding to the ${{\upomega }_{j}^{2}} > {\upomega }^{2}_{\min \limits }=250$ kHz² 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² constraint limit when u_y = 1000 mm.

Fig. 4

Optimized Fig. 3 a and e designs in their undeformed configurations

Fig. 5

The smallest eigenvalue as a function of the deflection u_z for the Fig. 3e design

Figure 6 illustrates designs subject to minimum eigenvalue constraints ${\upomega }^{2}_{\min \limits }=400$ kHz² and ${\upomega }^{2}_{\min \limits }=450$ kHz² 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.

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 ${{\upomega }_{j}^{2}} >{\upomega }^{2}_{\min \limits }$. a, b${\upomega }^{2}_{\min \limits } = 400$ kHz² and c, d${\upomega }^{2}_{\min \limits } = 450$ kHz². Every 5th iteration is marked in b and d, where red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, green—${{\upomega }^{2}_{3}}$, and black—$g_{o}/g_{o}^{\text {init}}$

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.

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

In Table 1, the values of g_o, $g_{o}^{\text {init}}$, ${{\upomega }_{1}^{2}}$, ${{\upomega }_{2}^{2}}$, and ${{\upomega }_{3}^{2}}$ 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 Fig. 6d design clearly illustrates the impact of the nonlinear framework.

Table 1 Resulting values for the double clamped beam designs

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 elements that are reintroduced during a design update rarely exceeds 10 elements in a 6250 element mesh.

Fig. 8

a Evolution of the ratio of removed finite elements n_rem with respect to the total number of elements. b The number of reintroduced elements n_rein during a design update. Both (a) and (b) correspond to the Fig. 3b design. Every 5th iteration is marked in (a)

10.2 Simply supported beam

In the last example, we consider the simply supported beam shown in Fig. 9.

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

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 ${{\upomega }_{j}^{2}} > {\upomega }_{\min \limits }^{2} = 150$ kHz² renders the design and its mode shapes illustrated in Fig. 11.

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—ϕ₁, d—ϕ₂, and e—ϕ₃ are plotted for u_y = 0 mm. Every 5th iteration is marked in b where red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, green—${{\upomega }^{2}_{3}}$, and black—$g_{o}/g_{o}^{\text {init}}$

Fig. 11

Optimized design a, evolutions of the eigenvalues and compliance ratio b, and first three eigenmodes c–e of a simply supported beam when u_y ≈ 0 mm subject to the eigenvalue constraints ${{\upomega }_{j}^{2}} >{\upomega }^{2}_{\min \limits } = 150$ kHz². Scaled c—ϕ₁, d—ϕ₂, and e—ϕ₃ are plotted for u_y = 0 mm. Every 5th iteration is marked in (b) where red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, and green—${{\upomega }^{2}_{3}}$

Figure 12 illustrates designs at the u_y = 1000mm prestrain subject to the eigenvalue constraints ${{\upomega }_{j}^{2}} > {\upomega }_{\min \limits }^{2}$ with ${\upomega }^{2}_{\min \limits }=150$ kHz² and ${\upomega }^{2}_{\min \limits }=250$ kHz². 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.

Fig. 12

Optimized designs and evolutions of the eigenvalues and compliance ratio for the geometry in Fig. 9 subject to the eigenvalue constraints ${{\upomega }_{j}^{2}} >{\upomega }^{2}_{\min \limits }$. a${\upomega }^{2}_{\min \limits }=150$ kHz² and b${\upomega }^{2}_{\min \limits }=250$ kHz². In the right hand side plots, every 5th iteration is marked, where red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, green—${{\upomega }^{2}_{3}}$, and black—$g_{o}/g_{o}^{\text {init}}$

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—ϕ₁, second row—ϕ₂, and third row—ϕ₃

The mode shapes associated with the Fig. 12 designs are plotted in Fig. 13. We observe that a local eigenmode highlighted in 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 ${{\upomega }_{j}^{2}}>{\upomega }^{2}_{\min \limits }=150$ kHz² eigenvalue constraint in the undeformed configurations in Fig. 14. The impact of the prestrain on the optimized topology is clear.

Fig. 14

Optimized Figs. 11 and 12a designs in their undeformed configurations.

In Table 2, the values of the response quantities of interest are listed for the various designs. It is concluded that the optimization problem $\mathcal {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.

Table 2 Resulting values in the last design iteration for the simply supported beam

11 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.

Appendix

1.1 Verification of sensitivities

The accuracy of the analytical sensitivity expression is verified by numerical differentiation using the forward finite difference method. Due to the finite precision in the eigenfrequency computations, the degeneracy must be determined to within a given tolerance ε_tol, cf. Seyranian et al. (1994). To set this tolerance, we compare the analytical directional derivative obtained when assuming simple or degenerate eigenfrequencies, i.e. $\overset {\text {S}}{{\delta {\upomega }_{j}^{2}}}$ when using Eq. 40 or $\overset {\text {D}}{{\delta {\upomega }_{j}^{2}}}$ when using Eq. 54, to the forward finite difference derivative approximation

$$ \overset{\text{FD}}{{\delta{\upomega}_{j}^{2}}} \approx \frac{{{\upomega}_{j}^{2}}(\textsf{z}_{1},...,\textsf{z}_{e} +\upvarepsilon,...,\textsf{z}_{n_{e}}) - {{\upomega}_{j}^{2}}(\textsf{z}_{1},...,\textsf{z}_{n_{e}})}{\upvarepsilon}, $$

(57)

where ε denotes the perturbation of variable z_e. Three different scenarios are then possible. Firstly, if $\overset {\text {S}}{{\delta {\upomega }_{j}^{2}}} \approx \overset {\text {FD}}{{\delta {\upomega }_{j}^{2}}}$ and $\overset {\text {D}}{{\delta {\upomega }_{j}^{2}}} \neq \overset {\text {FD}}{{\delta {\upomega }_{j}^{2}}}$, the eigenfrequency is treated as simple. Secondly, if $\overset {\text {D}}{{\delta {\upomega }_{j}^{2}}} \approx \overset {\text {FD}}{{\delta {\upomega }_{j}^{2}}}$ and $\overset {\text {S}}{{\delta {\upomega }_{j}^{2}}} \neq \overset {\text {FD}}{{\delta {\upomega }_{j}^{2}}}$, the eigenfrequency is treated as degenerate. Finally, if both $\overset {\text {S}}{{\delta {\upomega }_{j}^{2}}} \approx \overset {\text {D}}{{\delta {\upomega }_{j}^{2}}} \approx \overset {\text {FD}}{{\delta {\upomega }_{j}^{2}}}$, S is diagonal and the eigenvalue is degenerate, but Eqs. 40 and 54 yield identical result, cf. Seyranian et al. (1994).

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.

Fig. 15

Simply supported 3D plate. The prescribed displacement u_z is applied to the top-center surface in the negative z-direction

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 $\left ({\int \limits }_{{\Omega }_{0}} \nu dv\right )/V =V^{*}/V$ are denoted in Table 3.

Fig. 16

Optimized design a and first three eigenmodes b–d of simply supported 3D plate when u_z = 280 mm without eigenvalue constraints. Scaled b—ϕ₁, c—ϕ₂, and d—ϕ₃ are plotted for u_z = 0 mm

Fig. 17

Evolutions of the eigenvalues and compliance ratio value for the Fig. 16 design. Every 5th iteration is marked. Red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, green—${{\upomega }^{2}_{3}}$, and black—$g_{o}/g_{o}^{\text {init}}$

Table 3 Resulting values for the 3D plate designs

Next, we repeat the previous problem, but now we enforce the eigenvalue constraints ${{\upomega }^{2}_{j}} > {\upomega }^{2}_{\min \limits } = 900$ kHz². The resulting design is depicted in Fig. 18, where the mode shapes also are seen.

Fig. 18

Optimized design a and first three eigenmodes b–d of simply supported 3D plate when u_z = 280 mm subject to the eigenvalue constraints ${{\upomega }_{j}^{2}} >{\upomega }^{2}_{\min \limits } = 900$ kHz². Scaled b—ϕ₁, c—ϕ₂, and d—ϕ₃ are plotted for u_z = 0 mm

As is shown in Fig. 19, the twofold degenerate fundamental eigenvalue is increased to fulfill the lower bound ${\upomega }^{2}_{\min \limits } = 900$ kHz² 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.

Fig. 19

Evolutions of the eigenvalues and compliance ratio value for the Fig. 18 design. Every 5th iteration is marked. Red—${{\upomega }^{2}_{1}}$, blue—${{\upomega }^{2}_{2}}$, green—${{\upomega }^{2}_{3}}$, and black—$g_{o}/g_{o}^{\text {init}}$

Fig. 20

Evolution of the total volume fraction V^∗/V for the Fig. 18 design. Every 5th iteration is marked.