Level set-based isogeometric topology optimization for maximizing fundamental eigenfrequency

Maximizing the fundamental eigenfrequency is an efficient means for vibrating structures to avoid resonance and noises. In this study, we develop an isogeometric analysis (IGA)-based level set model for the formulation and solution of topology optimization in cases with maximum eigenfrequency. The proposed method is based on a combination of level set method and IGA technique, which uses the non-uniform rational B-spline (NURBS), description of geometry, to perform analysis. The same NURBS is used for geometry representation, but also for IGA-based dynamic analysis and parameterization of the level set surface, that is, the level set function. The method is applied to topology optimization problems of maximizing the fundamental eigenfrequency for a given amount of material. A modal track method, that monitors a single target eigenmode is employed to prevent the exchange of eigenmode order number in eigenfrequency optimization. The validity and efficiency of the proposed method are illustrated by benchmark examples.


Introduction
Topology optimization (TO), which has been extensively studied over the last decades, is a process of determining optimal layout of materials inside a given design domain. TO has been applied to various structural optimization problems, such as minimum compliance [1,2] vibration [3,4], and thermal issues [5,6], after Bendsøe and Kikuchi [7] proposed the homogenization method. Homogenization is a material distribution method in which a design domain is discretized into small rectangular elements, and each element contains an artificial composite material with microscopic voids. The proposal of the homogenization method was followed by a parallel exploration of the solid isotropic material with penalization (SIMP) method [8,9], which uses an artificial isotropic material whose physical properties are expressed as a function of continuous penalized material density (design variables). The phasefield method [10,11], which is also a material distribution method, is based on the theory of phase transitions. A different type of method called evolutionary structural optimization (ESO) [12,13] has also been proposed. This method eliminates elements with the lowest criterion value on the basis of certain heuristic criteria. ESO is computationally expensive because it requires a much larger number of iterations with an enormous number of intuitively generated solutions compared with materialbased methods.
However, these conventional TO methods, which are based on element-wise design variables, suffer from numerical instability problems, such as checkerboards and mesh dependency. Accordingly, several studies have proposed prevention methods. The use of high-order elements has been proven to be an effective means to prevent checkerboards [14,15], but this method entails a considerable increase in computation time. Various filter techniques have been utilized to mitigate checkerboards and mesh dependency because these techniques require only a small amount of extra computation time and are simpler to implement than other methods [16,17]. Notably, filter schemes are purely heuristic. Other prevention schemes, such as perimeter control and gradient constraint, which often make the optimization procedure unstable, are yet to be improved.
A new type of TO approach is the level set-based method (LSM), which was developed by Sethian and Wiegmann [18] to numerically track the motion of structural boundaries. In LSM, boundaries are represented as zero level set contour of an implicit high-dimensional function (level set function or LSF), in which boundary motion, merging, and introduction of new holes are performed. The evolution of a structural interface with respect to time is tracked by solving a Hamilton-Jacobi (H-J) partial differential equation (PDE), where transporting LSF along its outward normal direction is equivalent to moving the boundaries along the descent gradient direction. The conventional level set-based TO approach uses shape derivatives coupled with the original LSM for boundary tracking [19,20]. In this approach, regularization that reinitializes LSF to be signed distance to a zero level set is employed to control the slope of LSF. This conventional approach is updated by solving the H-J equation via an explicit up-wind scheme [21,22]. Variations of the conventional approach include parameterizing LSF using various basis functions, such as finite element method (FEM) basis functions [23], radial basis functions (RBFs) [24,25], and spectral parameterization [26], and corresponding methods for solving the H-J equation. By defining the interfaces between two material phases via the iso-contour of LSF, LSM can handle shape and topology changes during the optimization procedure and provides optimal structures with clear boundaries that are free of checkerboard patterns. Notably, most LSMs rely on finite elements wherein boundaries are still represented by discretized mesh in the analysis field unless alternative techniques are utilized to map the geometry to the analysis model.
Most of these TOs are performed in a fixed domain of finite elements where FEM is used to solve optimization problems. Currently used FEMs are often based on Lagrange polynomials for analysis while the geometrical representation of structures relies on non-uniform rational B-splines (NURBS), which are the criteria in computer aided design (CAD) systems. Thus, conversion of NURBS-based representation into one that is compatible with Lagrange polynomials, that is, mesh generation, is required in structural analysis. The disadvantages of FEM are as follows. First, the geometry approximation inherent in the FEM mesh may generate an approximate error. Second, frequent data interaction between geometry description and the computational mesh, which can be found in several calculations (e.g., fluid, large deformation, and shape optimization problems), is cumbersome and error-prone. An integrating method, namely, isogeometric analysis (IGA) [27], for unifying analysis and CAD processes has been proposed to overcome these disadvantages. This method employs the same basis functions as a technique for describing and analyzing the geometric model, which features the IGA method and CAD-based parameterization of field variables in an isoparametric manner. The first work on isogeometric approximation dates back to 1982 [28]; however, this method is considerably different from the IGA method. Several methods have been devised to help alleviate the difficulties faced by IGA. Special parameterization techniques, such as variational harmonic-based methods [29,30] and analysis-aware parameterization methods for single [31] and multi-domain geometries [32], have been proposed for the computational domain. Alternatives to NURBS, such as T-splines [33,34], polynomial splines over hierarchical T-meshes (PHT-splines) [35][36][37], and Powell-Sabin splines [38], have been studied for local refinement in IGA due to the limitation of the tensor product form of NURBS in computation refinement. Methods of parameterization of the interior domain while retaining the geometry exactness from the CAD model have been devised [39,40], and isogeometric collocation method is one of the most important among these methods [39]. With regard to interior discretization obstacles, the isogeometric boundary element method is a suitable candidate [41] because only boundary data are required for analysis, and it enables stress analysis [42], fracture analysis [43,44], acoustic analysis [45], and shape optimization [46,47]. Considering that the integral efficiency of IGA is limited by the tensor product structure of NURBS, an efficient quadrature rule, which is more suitable for NURBS-based IGA compared with the Gaussian quadrature rule, has been proposed in Ref. [48]. IGA has been applied to a wide range of problems, such as structural vibrations [49], fluid-structure interaction [50,51], heat conduction analyses [52], shape optimization [53,54], shell analyses [55], TO [56], and electromagnetics [57].
IGA-based shape optimization has been extensively investigated because remeshing is eliminated during the optimization process. IGA has also been recently applied to TO. The most commonly used TO approach is material based, and the most commonly solved TO problem is the minimum compliance case. In Ref. [58], TO was proposed based on isogeometric shape optimization. B-spline curves were introduced to represent the material boundary, and the coordinates of their control points were considered design variables. In Ref. [59], the trimmed spline surface technique was used for spline-based TO. In Refs. [60,61], optimality criteria (OC) and the method of moving asymptotes were implemented in the isogeometric-based SIMP method. In Ref. [62], the TO problem was solved by using a phase-field model, and IGA was utilized for the exact representation of the design domain. In Ref. [56], IGA was introduced to a level set TO method, and NURBS basis functions were used for geometry description and LSF parameterization.
Most studies on TO were restricted to compliance optimization, and the number of studies on TO of dynamic problems is limited. Dynamic problems, such as vibration and noise, are critical in many engineering fields, such as aeronautical and automotive industries. As a typical dynamic problem, structure vibration is controlled by the structure's dynamic characteristics, which are usually represented by the eigenfrequencies of the structure. Thus, eigenfrequency optimization plays an important role in improving dynamic characteristics.
As an important research topic, eigenfrequency optimization has been studied by many scholars. Díaaz and Kikuchi [63] extended TO to eigenfrequency optimization within the framework of the homogenization approach.
Tenek and Hagiwara [64] introduced the SIMP method to structural shape optimization and TO in eigenfrequency problems. Xie and Steven [65] applied the ESO scheme to TO, and the specified eigenfrequency of the structure was used as a constraint. Additional research has also been conducted on frequency optimization [66][67][68][69][70]. However, standard density-based TO methods are unsuitable for eigenfrequency optimization due to localized modes in low-density areas [71]. Low-density areas are much more flexible than areas with full densities; hence, they control the lowest eigenmodes of the entire structure. By changing the penalization of stiffness in the SIMP method, a modified algorithm has been proposed and applied to circumvent this numerical instability problem [3]. LSM that employs a crisp description of structure boundaries has advantages over the density-based approach. The LSM approach can avoid artificial modes localized in the weak phase, which makes LSM a choice for eigenfrequency optimization. Many studies have been performed on level set TO for eigenfrequency problems [72][73][74][75].
IGA-based TO has been extensively studied. However, research on the combination of the level set approach and IGA is limited, and that on eigenfrequency problems is absent. In this study, we develop a new optimization method to formulate the TO problem for cases with maximum fundamental eigenfrequency by using LSM under the framework of IGA instead of conventional finite element analysis (FEA). Given that the OC algorithm is particularly efficient for problems with many design variables and few constraints, we consider the OC method for the solution of the optimization problem and conduct a sensitivity analysis. The rest of this paper is organized as follows. In Section 2, IGA and NURBS-based TO are summarized. The eigenvalue optimization problem is described in Section 3, and the TO model is proposed in Section 4. Section 5 presents numerical examples to demonstrate the validity of the proposed approach. The conclusions and discussions are shown in Section 6.

IGA for level set-based TO
In IGA, geometric modeling and analysis are integrated by using NURBS, where the basis function is a bridge of the parameter domain, physical field, and numerical solution.
In the proposed method, we consider the NURBS basis function as a bound between IGA and parameterized LSM. We provide a brief review of IGA and NURBSparameterized LSM.

NURBS-based IGA
We assume that geometrical mapping Ψ maps the parameter domainΩ into the physical domain Ω. Given a knot vector Ξ ¼ ð 1 , 2 ,:::, s Þ with a non-decreasing sequence of values lying in parameter space, the mapping between two domains can be expressed as The NURBS curve is constructed by linear combinations of its basis functions, in which the coefficients are a given set of control points. A NURBS curve of p-degree is defined as where n ¼ sp -1 is the number of control points, P i 2 ℝ d is the ith control points in the physical domain, and R i,p ðÞ is the ith univariate NURBS basis function defined in the parameter spaceΩ as follows: where ω i 2 ð0,1Þ are non-decreasing weights associated with control points and N i,p ðÞ represents the ith B-spline basis functions of p-degree; it is defined by the following Cox-de Boor recursion formula [27].
where basis function N i,p ðÞ has its own support domain ½ i ,:::, iþpþ1 in which N i,p ðÞ is non-zero. A knot vector is deemed open when the knots are repeated p þ 1 times at the ends of the vector. In IGA, the open-knot vector is generally used to satisfy the Kronecker delta property at boundary points. A NURBS surface is a tensor product of bivariate NURBS curves in Ξ and H directions with pand qdegrees, respectively, where a knot vector H ¼ ðη 1 ,η 2 ,:::, η t Þ is given in H direction.
where m ¼ tq -1, P i,j are control points and R p,q i,j are bivariate basis functions of the form: where N i,p ðÞ and N j,q ðηÞ are B-spline basis functions defined on the knot vector Ξ ¼ ð 1 , 2 ,:::, s Þ and H ¼ ðη 1 ,η 2 ,:::,η t Þ, respectively. The interval Ξ Â H forms a patch containing all elements, namely, ½ k , kþ1 Â ½η l ,η lþ1 , 1£k£n þ p, and 1£l£m þ q, which are defined by the two knot vectors. The parameter domain and corresponding physical domain for a surface model are depicted in Fig. 1.
On the basis of the isoparametric concept, the IGA approach utilizes the same parameters for geometry and analysis models, and the basis functions used for geometry representation are also employed to approximate the numerical solution of PDEs. With Eq. (2), the numerical solution u can be expressed as where R i is the ith basis function. u i , which is referred to as a control variable at the ith control point, is the coefficient used to approximate the field variable u, which plays the same role as the nodal value in FEA. For each element, the shape function and strain-displacement matrix can be expressed as : The strain matrix is given by

NURBS-based parameterized LSM
LSM is a TO method that implicitly defines the interfaces between material and void phases by iso-contours of a high-dimensional LSF. Thus, this method allows a crisp description of the material boundary and helps avoid meshdependent problems. The shape of the interpolating functions of LSF directly influences the smoothness of LSF and the material domain. In its most general form, LSF is described by a summation of interpolating functions scaled by their degree of freedom (DOF).
where x denotes the spatial coordinate, f i comprises interpolating functions associated with N spatial points, and s i are time-dependent optimization variables. The most commonly used interpolation functions in present LSM are FEM shape functions and RBFs, and their corresponding optimization variables are nodal values and expansion coefficients, respectively. We introduce NURBS basis functions, which can be used to approximate a given set of points with smooth polynomial functions, for parameterizing LSF. Thus, the NURBS-based parameterized LSF is constructed as where R i is the ith NURBS basis function and s i is the ith time-dependent expansion coefficient related to the ith control point. The evolution of LSF is governed by the following H-J equation. Fig. 1 Geometrical mapping Ψ maps the common parameter space ð,ηÞ onto the physical space where t denotes pseudo-time, which represents the evolution of the design in the optimization process. The speed of movement of a point on the level set surface can be expressed by V ¼ dx=dt. V n ¼ V ⋅n defines the speed of propagation of all level sets along the external normal velocity, where n ¼ -rΦ=jrΦj. Therefore, Eq. (12) can be rewritten as ∂Φðx,tÞ ∂t ¼ V n jrΦj: By substituting Eq. (11) into Eq. (13), the H-J equation can be written as The moving speed of the material free boundary during evolution is related to the time derivative of the expansion coefficient as follows: 3 Optimization problems of maximizing eigenvalue

Definition of the eigenvalue problem
We describe the eigenvalue problems in the linear elastic continuum to facilitate the computation of vibration frequencies and modes. A linear elastic continuum structure with a constant mass density is defined in domain Ω & ℝ d ðd ¼ 2 or 3Þ with the boundary Γ ¼ ∂Ω. The weak formulation of the undamped free vibration problem can be expressed as where eigenfrequency l and corresponding eigenvector u, that is, the displacement subdomain in Ω, are solutions of this eigenvalue problem, v is adjoint displacement, which satisfies the kinematic boundary condition, and U is a space of kinematically admissible displacement fields. In LSM, að⋅,⋅Þ and bð⋅,⋅Þ are respectively defined as where D ijkl stands for the elasticity tensor component, ε ij is the strain tensor component, is the density of the material, and HðΦÞ is the Heaviside function, which takes 0 for Φ < 0 and 1 otherwise. The eigenvalue problem has a family of solutions l k and u k , k³1. The first eigenfrequency and its eigenvector are related to each other as where J ðu,ΦÞ is the objective function, V max represents the maximum admissible volume of the design domain, and s min and s max stand for the lower and upper bounds of the design variables, respectively. However, in the eigenfrequency optimization process, the value of higher-order eigenfrequency may decrease whereas that of lower-order target eigenfrequency may increase, which may possibly lead to the repetition and exchange of mode order number. Given that the objective and constraint functions are typically defined based on a fixed modal order, the sensitivities of these functions are discontinuous in the repeated eigenfrequency case. Approaches are often used to maintain the simplicity of eigenfrequency during the entire optimization process and overcome this ill-posed problem. The modal assurance criterion (MAC) method, an efficient and accurate strategy, is introduced to monitor a single target mode, which is the first mode in this study. The definition of MAC is where u a and u b represent two eigenvectors: One is the reference eigenvector of the current optimization cycle and the other is the objective eigenvector of the previous cycle in the eigenvalue optimization process. The value of MAC varies between 0 and 1. Theoretically, a MAC value of 1 means that the two eigenvectors representing modal shapes are exactly the same. However, this condition is impossible because the structural configuration changes in each iteration, and the modal shapes of adjacent iterations are not orthogonal to each other. By comparing a few reference eigenvectors with an objective eigenvector u obj , a new objective eigenvector is obtained in each iteration step of the optimization process, and it can be expressed as u n obj ¼ u n k that max k ½MACðu n-1 obj ,u n k Þ, k ¼ 1,2,:::,N m , (22) where N m is the number of modal eigenvectors that need to be checked, superscript of displacement indicates the number of iteration step, that is, in the nth iteration, objective eigenvector of previous iteration step u n-1 obj is used to calculate the MAC value.

Sensitivity analysis
Establishing the relationship between the optimization function and design variables by using a sensitivity analysis approach is necessary to solve the optimization problem. According to the material derivative and the adjoint method, the Lagrangian function can be defined as where Λ is the Lagrangian multiplier. Assuming that V ðu,ΦÞ ¼ ! Ω HðΦÞdΦ -V max is the volume constraint, the shape derivative of Lagrangian function Lðu,ΦÞ is where the material derivatives of aðu,v,ΦÞ and bðu,v,ΦÞ are respectively given by where u # and v # are partial derivatives of u and v, respectively, with respect to pseudo-time. δðxÞ is the Dirac function.
The adjoint state equation can be obtained by the Kuhn-Tucker condition.

Numerical implementation
Many design variables, which correspond to large-scale nonlinear equations in the eigenvalue problem, exist in continuum structural TO. Thus, OC is introduced to solve this eigenvalue TO problem. By properly iterating and updating the design variables, this optimization problem is guaranteed to converge to a final solution. Starting from an initial value, the iterative formula of the design variables is expressed as Theoretically, the iteration coefficients c where is a very small number that can avoid singularity when the sensitivity of the volume constraint with respect to the design variables is equal to 0. The Lagrangian multiplier Λ is calculated by the bisection method [76].
A flowchart of the structural TO for maximization of the first eigenfrequency problem is depicted in Fig. 2. Given the condition of constraints, when the relative difference value of the objective function in the current and previous iterations is less than 10 -3 , this optimization process is considered convergent, and the current optimization process is terminated.

Numerical examples
In this section, the proposed IGA-based level set TO framework is applied to two 2D optimization problems.
For all examples, the properties of the isotropic material are set as follows: Young's modulus E ¼ 210 GPa and mass density ¼ 7:8 Â 10 3 kg=m 3 . The properties of the artificial weak material are E 0 ¼ 210 Â 10 -3 GPa and mass density ¼ 7:8 Â 10 -3 kg=m 3 . Poisson's ratio ¼ 0:3 and plane stress state are assumed for all the materials.
The two examples are TO of maximizing the fundamental eigenfrequency of the plane structure with a unit thickness of 0:001 m and a prescribed material volume fraction of α ¼ 50%. In the initial design, the available material is uniformly distributed over the entire admissible design domain. In the following examples, the boundary conditions are imposed by the collocation method, which enforces these conditions to be satisfied pointwise. Given that NURBS basis functions associated with the interior control points vanish at the structural boundary when open-knot vectors are employed, the displacement boundary condition applied on the left and right sides of the beam is imposed by setting the displacement values at left and right boundary control points to zero. All results are produced with programs developed in the MATLAB R2018a environment on a computer with an Intel Core i3-3240 CPU, 3.4 GHz clock speed, and 6 GB RAM. Additional details on the implementation of IGA on the MATLAB platform were presented in Ref. [77].

Cantilever beam
The first numerical example of a short cantilever beam for maximizing the first eigenfrequency optimization problem is shown in Fig. 3. The entire design domain is a rectangle with a size of 0.2 mÂ0.1 m, a Dirichlet boundary, and fixed displacement at the left edge of the design domain. A concentrated nonstructural mass M ¼ 15:6 kg that is onetenth of the total structural mass of the plate is placed at the center of the right side. Notably, the structure disappears without the nonstructural mass because no structure leads to the highest eigenfrequencies.
In this example, IGA-and FEA-based TO approaches The DOF of IGA is much less than that of FEM (nearly a quarter of FEM) when the element number is sufficiently large. Table 1 shows that the computation efficiency of IGA-based LSM is higher than that of FEM-based LSM in TO due to the fewer DOFs or smaller size of equations in the IGA-based approach. However, because of the extra calculations of the basis functions and their derivative in the IGA-based approach, the ratio between the iteration time and DOFs of the IGA-based optimization method is higher than that of the FEM-based method.
The optimal layouts obtained by using IGA and FEA with 128Â64 elements are shown in Fig. 4. In this case, the results are similar, and the crisp boundary is obtained due to the level set method. Adopting B-spline basis functions to parameterize the level set function together with IGA for calculation instead of FEA leads to the rapid convergence of the IGA-based optimization process.
The convergence history of optimization using IGA and FEA with 128 Â 64 elements is shown in Figs. 5 and 6. Figure 5 shows the convergence history of the first eigenfrequency and volume ratio of the structure by using IGA-and FEA-based optimization frameworks, respectively. The initial designs and resultant structures of both methods are basically the same. In optimization with the IGA method, the ω i of the initial design and the resultant optimum are 173.5 and 188.7, respectively. The volume ratio of the initial design and the resultant optimum are 0.79 and 0.5, respectively. Thus, the fundamental eigenfrequency increased by 8.8% and the volume decreased by 36.7%. Figure 6 shows the iteration history of the first three eigenfrequencies by using both methods. The first eigenfrequency always remains simple, whereas the second and third eigenfrequencies tend to oscillate and  overlap. By using the MAC method, problems regarding the orders of eigenfrequency exchange during the optimization process are avoided. The first eigenfrequency generally increases, and the second and third eigenfrequencies decrease as the volume ratio decreases.

Beam with clamped ends
In this section, we present an example of maximizing the fundamental frequency of a clamped beam structure shown in Fig. 7. The working domain has a size of 0.4 mÂ0.1 m.
A fixed displacement boundary condition is imposed on both sides, and a concentrated nonstructural mass M = 31.2 kg is placed at the center of structure. In this example, the capability of the proposed method to capture the optimum topology and the effect of the number of elements are studied. For this purpose, a clamped beam is solved with three mesh sizes, namely, 64Â16, 128Â32, and 256Â64.
The degree of NURBS basis function is 2 in both directions. The resulting layouts shown in Fig. 8 indicate that the mesh dependency problem is avoided in the proposed method. This figure also shows that an inaccurate optimum topology with a rough boundary is obtained as a consequence of coarse meshes. By refining the mesh, the smoothness of the boundaries of the optimal layout is improved. However, when the number of elements is larger than 128Â32, the resulting layout slightly changes. Accounting for the computation cost and precision of results, 128Â32 meshes are used in this example.
The convergence history of the objective function and volume ratio is given in Fig. 9. The first eigenfrequency of the initial design and the resultant topology are 244.3 and 257.4, respectively. The volume ratio of the initial design and the resultant topology are 0.82 and 0.5, respectively. Thus, the fundamental frequency increases by 5.4%, and the volume decreases by 39%. Figure 10 shows the convergence history of the first three eigenfrequencies. The fundamental frequency remains simple throughout the entire optimization process. The value of first eigenfre-

Conclusions
We solved maximum fundamental eigenfrequency TO problems with a level set model based on the IGA technique. IGA combines the fundamental idea of FEM with the spline technique from a computer-aided geometry design for the integration of CAD and CAE. The IGA method was also introduced to TO due to its superiority over currently used FEM in terms of accuracy and efficiency. The feature of the proposed method is the combination of IGA and LSM in eigenfrequency optimization where the same basis functions (NURBS) are used for geometry representation, dynamic analysis, and parameterization of the implicit LSF. High accuracy and smoothness of LSF were achieved by using smooth NURBS basis functions to approximate LSF. In the case of maximizing the fundamental eigenfrequency, a regularization method for the repetition or exchange of eigenfrequencies was employed to guarantee the simple behavior of structural eigenfrequency.
Two benchmark numerical examples of TO for dynamic problems were applied to verify the validity of the proposed approach. The results obtained from the comparison of FEA-and IGA-based level set TO methods demonstrated that the proposed IGA-based optimization method has better computational efficiency and converges faster than the traditional FEA-based optimization method. The results also showed that solving dynamic TO problems by using IGA together with the level set method is possible. Although we only presented examples of 2D structures, no theoretical difficulties will be encountered if the proposed is extended to the optimization of 3D structures.