Level set-based isogeometric topology optimization for maximizing fundamental eigenfrequency

Abstract

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.

Key words： topology optimization; level set method; isogeometric analysis; eigenfrequency

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 [¹,²] vibration [³,⁴], and thermal issues [⁵,⁶], after Bendsøe and Kikuchi [⁷] 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 [⁸,⁹], which uses an artificial isotropic material whose physical properties are expressed as a function of continuous penalized material density (design variables). The phase-field method [¹⁰,¹¹], 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) [¹²,¹³] 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 material-based 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 [¹⁴,¹⁵], 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 [¹⁶,¹⁷]. 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 [¹⁸] 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 [¹⁹,²⁰]. 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 [²¹,²²]. Variations of the conventional approach include parameterizing LSF using various basis functions, such as finite element method (FEM) basis functions [²³], radial basis functions (RBFs) [²⁴,²⁵], and spectral parameterization [²⁶], 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) [²⁷], 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 [²⁸]; 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 [²⁹,³⁰] and analysis-aware parameterization methods for single [³¹] and multi-domain geometries [³²], have been proposed for the computational domain. Alternatives to NURBS, such as T-splines [³³,³⁴], polynomial splines over hierarchical T-meshes (PHT-splines) [³⁵–³⁷], and Powell-Sabin splines [³⁸], 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 [³⁹,⁴⁰], and isogeometric collocation method is one of the most important among these methods [³⁹]. With regard to interior discretization obstacles, the isogeometric boundary element method is a suitable candidate [⁴¹] because only boundary data are required for analysis, and it enables stress analysis [⁴²], fracture analysis [⁴³,⁴⁴], acoustic analysis [⁴⁵], and shape optimization [⁴⁶,⁴⁷]. 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. [⁴⁸]. IGA has been applied to a wide range of problems, such as structural vibrations [⁴⁹], fluid-structure interaction [⁵⁰,⁵¹], heat conduction analyses [⁵²], shape optimization [⁵³,⁵⁴], shell analyses [⁵⁵], TO [⁵⁶], and electromagnetics [⁵⁷].

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. [⁵⁸], 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. [⁵⁹], the trimmed spline surface technique was used for spline-based TO. In Refs. [⁶⁰,⁶¹], optimality criteria (OC) and the method of moving asymptotes were implemented in the isogeometric-based SIMP method. In Ref. [⁶²], 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. [⁵⁶], 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 [⁶³] extended TO to eigenfrequency optimization within the framework of the homogenization approach. Tenek and Hagiwara [⁶⁴] introduced the SIMP method to structural shape optimization and TO in eigenfrequency problems. Xie and Steven [⁶⁵] 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 [⁶⁶–⁷⁰]. However, standard density-based TO methods are unsuitable for eigenfrequency optimization due to localized modes in low-density areas [⁷¹]. 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 [³]. 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 [⁷²–⁷⁵].

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 NURBS-parameterized LSM.

NURBS-based IGA

We assume that geometrical mapping \mathrm{\Psi } maps the parameter domain \widehat{\mathrm{\Omega }} into the physical domain \mathrm{\Omega }. Given a knot vector \mathbf{\Xi }=({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\dots ,{\xi }_s) with a non-decreasing sequence of values lying in parameter space, the mapping between two domains can be expressed as

(1)\mathrm{\Psi }:\widehat{\mathrm{\Omega }}\to \mathrm{\Omega },\mathbf{\Xi }\to \mathrm{\Psi }(\mathbf{\Xi }).

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

(2)\mathrm{\Psi }(\xi )={\displaystyle \sum _{i=\mathrm{1}}^n{R}_{i,p}}(\xi )P_i,

where n=s-p-\mathrm{1} is the number of control points, P_i\in {ℝ}^d is the ith control points in the physical domain, and R_{i,p}(\xi ) is the ith univariate NURBS basis function defined in the parameter space \widehat{\mathrm{\Omega }} as follows:

(3)R_{i,p}(\xi )={\displaystyle \frac{N_{i,p}(\xi ){\omega }_i}{{\displaystyle \sum _{i=\mathrm{1}}^n{N}_{i,p}}(\xi ){\omega }_i}},

where {\omega }_i\in (\mathrm{0},\mathrm{1}) are non-decreasing weights associated with control points and N_{i,p}(\xi ) represents the ith B-spline basis functions of p-degree; it is defined by the following Cox-de Boor recursion formula [²⁷].

(4)

where basis function N_{i,p}(\xi ) has its own support domain [{\xi }_i,\dots ,{\xi }_{i+p+\mathrm{1}}] in which N_{i,p}(\xi ) is non-zero. A knot vector is deemed open when the knots are repeated p+\mathrm{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 \mathbf{\Xi } and \mathrm{H} directions with p- and p-degrees, respectively, where a knot vector \mathrm{H}=({\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\dots ,{\eta }_t) is given in \mathrm{H} direction.

(5)\mathrm{\Psi }(\xi ,\eta )={\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \sum _{j=\mathrm{1}}^m{R}_{i,j}^{p,q}}}(\xi ,\eta )P_{i,j},

where m=t-q-\mathrm{1}, P_{i,j} are control points and R_{i,j}^{p,q} are bivariate basis functions of the form:

(6)R_{i,j}^{p,q}={\displaystyle \frac{N_{i,p}(\xi )N_{j,q}(\eta ){\omega }_{i,j}}{{\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \sum _{j=\mathrm{1}}^m{N}_{i,p}}}(\xi )N_{j,q}(\eta ){\omega }_{i,j}}},

where N_{i,p}(\xi ) and N_{j,q}(\eta ) are B-spline basis functions defined on the knot vector \mathbf{\Xi }=({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\dots ,{\xi }_s) and\mathrm{H}=({\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\dots ,{\eta }_t), respectively. The interval \mathbf{\Xi }\times \mathrm{H} forms a patch containing all elements, namely, [{\xi }_k,{\xi }_{k+\mathrm{1}}]\times [{\eta }_l,{\eta }_{l+\mathrm{1}}], \mathrm{1}\le k\le n+p, and \mathrm{1}\le l\le 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.

Fig.1 Geometrical mapping \mathrm{\Psi }maps the common parameter space (\xi ,\eta ) onto the physical space

Full size|PPT slide

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

(7)u={\displaystyle \sum _{i=\mathrm{1}}^n{R}_i}(\xi )u_i,

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

(8)\mathit{R}=[R_{\mathrm{1}}{R}_{\mathrm{2}}\mathrm{...}{R}_n],\mathit{B}=\left[\begin{array}{lllll}{\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial x}}& \mathrm{0}& \mathrm{...}& {\displaystyle \frac{\partial R_n}{\partial x}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial y}}& \mathrm{...}& \mathrm{0}& {\displaystyle \frac{\partial R_n}{\partial y}}\\ {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial y}}& {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial x}}& \mathrm{...}& {\displaystyle \frac{\partial R_n}{\partial y}}& {\displaystyle \frac{\partial R_n}{\partial x}}\end{array}\right].

The strain matrix is given by

(9)\mathit{\epsilon }=\mathit{B}\mathit{u}.

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

(10)\mathrm{\Phi }(\mathit{x})={\varphi }^{\mathrm{T}}(\mathit{x})\mathit{s}={\displaystyle \sum _i^N{\varphi }_i(\mathit{x})s_i},

where \mathit{x} denotes the spatial coordinate, {\varphi }_i comprises interpolating functions associated with Nspatial 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

(11)\mathrm{\Phi }(x)={\mathit{R}}^{\mathrm{T}}\mathit{s}={\displaystyle \sum _i^N{R}_i(\mathit{x})s_i},

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.

(12){\displaystyle \frac{\partial \mathrm{\Phi }(\mathit{x},t)}{\partial t}}+\nabla \mathrm{\Phi }{\displaystyle \frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}}=\mathrm{0},

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=\mathrm{d}\mathit{x}/\mathrm{d}t. {\mathit{V}}_{\mathrm{n}}=V\cdot \mathit{n} defines the speed of propagation of all level sets along the external normal velocity, where \mathit{n}=-\nabla \mathrm{\Phi }/|\nabla \mathrm{\Phi }|. Therefore, Eq. (12) can be rewritten as

(13){\displaystyle \frac{\partial \mathrm{\Phi }(\mathit{x},t)}{\partial t}}={\mathit{V}}_{\mathrm{n}}|\nabla \mathrm{\Phi }|.

By substituting Eq. (11) into Eq. (13), the H-J equation can be written as

(14){\mathit{R}}^{\mathrm{T}}{\displaystyle \frac{\partial \mathit{s}}{\partial t}}+{\mathit{V}}_{\mathrm{n}}\cdot \left|{(\nabla \mathit{R})}^{\mathrm{T}}\mathit{s}\right|=\mathrm{0.}

The moving speed of the material free boundary during evolution is related to the time derivative of the expansion coefficient as follows:

(15){\mathit{V}}_{\mathrm{n}}=-{\displaystyle \frac{{\mathit{R}}^{\mathrm{T}}}{\left|{(\nabla \mathit{R})}^{\mathrm{T}}\mathit{s}\right|}}{\displaystyle \frac{\partial \mathit{s}}{\partial t}}.

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 \mathrm{\Omega }\subset {ℝ}^d(d=\mathrm{2}\mathrm{o}\mathrm{r}\mathrm{3}) with the boundary \mathrm{\Gamma }=\partial \mathrm{\Omega }. The weak formulation of the undamped free vibration problem can be expressed as

(16)a(\mathit{u},\mathit{v})-\lambda b(\mathit{u},\mathit{v})=\mathrm{0},\mathit{v}\in U,

where eigenfrequency \lambda and corresponding eigenvector \mathit{u}, that is, the displacement subdomain in \mathrm{\Omega }, are solutions of this eigenvalue problem, \mathit{v} is adjoint displacement, which satisfies the kinematic boundary condition, and U is a space of kinematically admissible displacement fields. In LSM, a(\cdot ,\cdot ) and b(\cdot ,\cdot ) are respectively defined as

(17)a(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}}{\epsilon }_{kl}(\mathit{u}){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega },

(18)b(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}\rho }\mathit{u}\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega },

where D_{ijkl} stands for the elasticity tensor component, {\epsilon }_{ij} is the strain tensor component, \rho is the density of the material, and H(\mathit{\Phi }) is the Heaviside function, which takes \mathrm{0} for and \mathrm{1} otherwise.

The eigenvalue problem has a family of solutions {\lambda }_k and u_k, k\ge \mathrm{1}. The first eigenfrequency and its eigenvector are related to each other as

(19){\lambda }_{\mathrm{1}}=\min {\displaystyle \frac{{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}}{\epsilon }_{kl}(\mathit{u}){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}{{\displaystyle {\int }_{\mathrm{\Omega }}\rho }\mathit{u}\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}}.

Optimization model

We consider the TO problem by maximizing the first eigenfrequency under a volume constraint. Under the NURBS-based level set framework, eigenfrequency TO can be expressed as

(20)\begin{array}{l}\mathrm{Maximize}\mathrm{\hspace{1em}\hspace{1em}}J(\mathit{u},\mathrm{\Phi })={\lambda }_{\mathrm{1}}\\ \mathrm{subject}\hspace{0.17em}\mathrm{to}:\hspace{1em}\hspace{0.17em}a(\mathit{u},\mathit{v},\mathrm{\Phi })={\lambda }_{\mathrm{1}}b(\mathit{u},\mathit{v},\mathrm{\Phi }),\\ {\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }\le V_{\max },\\ s_{\min }\le s_i\le s_{\max },\mathrm{\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}} \end{array}

where J(\mathit{u},\mathrm{\Phi }) 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

(21)\mathrm{MAC}({\mathit{u}}_a,{\mathit{u}}_b)={\displaystyle \frac{|{\mathit{u}}_a^{\mathrm{T}}{\mathit{u}}_b{|}^{\mathrm{2}}}{({\mathit{u}}_a^{\mathrm{T}}{\mathit{u}}_a)({\mathit{u}}_b^{\mathrm{T}}{\mathit{u}}_b)}},

where {\mathit{u}}_a and {\mathit{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 {\mathit{u}}_{\mathrm{obj}}^{}, a new objective eigenvector is obtained in each iteration step of the optimization process, and it can be expressed as

(22){\mathit{u}}_{\mathrm{obj}}^n={\mathit{u}}_k^n\hspace{1em}\mathrm{that} \underset{k}{\max }[\mathrm{MAC}({\mathit{u}}_{\mathrm{obj}}^{n-\mathrm{1}},{\mathit{u}}_k^n)], k=\mathrm{1},\mathrm{2},\dots ,N_m,

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 {\mathit{u}}_{\mathrm{obj}}^{n-\mathrm{1}} 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

(23)L(\mathit{u},\mathrm{\Phi })=J(\mathit{u},\mathrm{\Phi })+a(\mathit{u},\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b(\mathit{u},\mathit{v},\mathrm{\Phi })+\mathrm{\Lambda }\left[{\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }-V_{\max }\right],

where \mathrm{\Lambda } is the Lagrangian multiplier. Assuming that V(\mathit{u},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }-V_{\max } is the volume constraint, the shape derivative of Lagrangian function L(\mathit{u},\mathrm{\Phi }) is

(24){\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}={{\lambda }_{\mathrm{1}}}^{\prime }+a^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}{b}^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })-{{\lambda }_{\mathrm{1}}}^{\prime }b(\mathit{u},\mathit{v},\mathrm{\Phi })+\mathrm{\Lambda }{V}^{\prime }(\mathit{u},\mathrm{\Phi }) ,

where the material derivatives of a(\mathit{u},\mathit{v},\mathrm{\Phi }) and b(\mathit{u},\mathit{v},\mathrm{\Phi }) are respectively given by

(25)\begin{array}{l}{a}^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}\\ +{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}({\mathit{v}}^{\prime })H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}\\ +{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}(\mathit{v})\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega ，}}\end{array}

(26)b^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}\rho {\mathit{u}}^{\prime }\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\mathrm{\Omega }}\rho \mathit{u}{\mathit{v}}^{\prime }H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\mathrm{\Omega }}\rho \mathit{u}\mathit{v}\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega }},

where {\mathit{u}}^{\prime } and {\mathit{v}}^{\prime } are partial derivatives of \mathit{u} and \mathit{v}, respectively, with respect to pseudo-time. \delta (x) is the Dirac function.

The adjoint state equation can be obtained by the Kuhn-Tucker condition.

(27)a({\mathit{u}}^{\prime },\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b({\mathit{u}}^{\prime },\mathit{v},\mathrm{\Phi })=\mathrm{0},

(28)a(\mathit{u},{\mathit{v}}^{\prime },\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b(\mathit{u},{\mathit{v}}^{\prime },\mathrm{\Phi })=\mathrm{0},

(29)\mathrm{1}-b(\mathit{u},\mathit{v},\mathrm{\Phi })=\mathrm{0.}

Given that the real mode \mathit{u} is equal to the adjoint mode \mathit{v}, Eq. (29) is a normalization condition for the eigenvector. In this case, Eq. (24) can be simplified as

(30){\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}={\displaystyle {\int }_{\mathrm{\Omega }}(F(\mathit{u})+\mathrm{\Lambda })\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega }},

where F(\mathit{u})=D_{ijkl}\mathit{\epsilon }(\mathit{u})\mathit{\epsilon }(\mathit{u})-{\lambda }_{\mathrm{1}}\rho \mathit{u}\mathit{u}. By substituting Eqs. (11) and (15) into the preceding shape derivative equation, we obtain

(31){\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}=-{\displaystyle {\int }_{\mathrm{\Omega }}(F(\mathit{u})+\Lambda )\delta (\mathrm{\Phi })\mathit{R}{\displaystyle \frac{\partial s}{\partial t}}\mathrm{d}\mathrm{\Omega }}.

Given that

(32){\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial s}}={\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s}}+\mathrm{\Lambda }{\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s}}={\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial s}}{\displaystyle \frac{\partial s}{\partial t}},

the sensitivity of the objective function and volume constraint with respect to the design variables is respectively obtained as follows:

(33){\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s}}=-{\displaystyle {\int }_{\mathrm{\Omega }}F(\mathit{u})\delta (\mathrm{\Phi })\mathit{R}\mathrm{d}\mathrm{\Omega }},

(34){\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s}}=-{\displaystyle {\int }_{\mathrm{\Omega }}\delta (\mathrm{\Phi })\mathit{R}\mathrm{d}\mathrm{\Omega }}.

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

(35)s_i^{(k+\mathrm{1})}=c_i^{(k)}{s}_i^{(k)}.

Theoretically, the iteration coefficients c_i^{(k)} are obtained by setting Eq. (32) equal to 0, which can be written as

(36)c_i^{(k)}=-{\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s_i^{(k)}}}/\max \left\{\mu ,{\mathrm{\Lambda }}^{(k)}{\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s_i^{(k)}}}\right\},

where \mu 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 \mathrm{\Lambda } is calculated by the bisection method [⁷⁶].

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⁻³, this optimization process is considered convergent, and the current optimization process is terminated.

Fig.2 Flowchart of the optimization procedure

Full size|PPT slide

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=\mathrm{210}\mathrm{GPa}and mass density \rho =\mathrm{7.8}\times {\mathrm{10}}^{\mathrm{3}}{\mathrm{kg/m}}^3. The properties of the artificial weak material are E_{\mathrm{0}}=\mathrm{210}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{GPa} and mass density \rho =\mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{kg/m}}^3. Poisson’s ratio \nu =\mathrm{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 \mathrm{0.001}\mathrm{m} and a prescribed material volume fraction of \alpha =\mathrm{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. [⁷⁷].

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=\mathrm{15.6} kg that is one-tenth 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.

Fig.3 Design domain of a cantilever beam structure

Full size|PPT slide

In this example, IGA- and FEA-based TO approaches are applied to a similar problem for comparison. For this 2D example, equally spaced open-knot vectors are used for x and y directions, and the degrees of NURBS basis functions on the two directions are the same, i.e., p=q=\mathrm{2}. Nine-node quadratic rectangle elements are used for FEA, and the element number of both methods is the same, which facilitates a comparison under identical conditions.

The computation consumption times of the two methods are shown in Table 1. Notably, for simplicity in the proceeding table and figures, IGA-L and FEM-L represent IGA- and FEM-based LSM, respectively. The solution time of the system equation includes the time spent on assembling the stiffness matrix and solving the equation to obtain an objective function. In this case, when the numbers of Lagrange and NURBS elements are a\times b (a and b represent the number of elements in x and y directions, respectively), the number of DOFs of the IGA-based method is (a+\mathrm{2})(b+\mathrm{2})=ab+\mathrm{2}a+\mathrm{2}b+\mathrm{4} and that of FEM is (\mathrm{2}a+\mathrm{1})(\mathrm{2}b+\mathrm{1})=\mathrm{4}ab+\mathrm{2}a+\mathrm{2}b+\mathrm{1}. 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.

Tab.1 Comparison of IGA- and FEA-based LSM TO

Method	Number of elements	Number of DOFs	Time of each iteration/s	Time of solution of system equation/s
IGA-L	256×128	67080	873.26	713.68
FEA-L	256×128	263682	1580.42	1237.16
IGA-L	128×64	17160	30.06	25.31
FEA-L	128×64	66306	54.17	43.97
IGA-L	64×32	4488	6.41	5.18
FEA-L	64×32	16770	8.55	6.85

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.

Fig.4 Optimized results obtained by using IGA and FEA methods with 128 × 64 elements

Full size|PPT slide

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 {\omega }_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.

Fig.5 Comparison of IGA and FEA in terms of convergence history

Full size|PPT slide

Fig.6 Comparison of IGA and FEA in terms of eigenfrequency history

Full size|PPT slide

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.

Fig.7 Design domain of a clamped beam structure

Full size|PPT slide

Fig.8 Optimal layouts obtained by using (a) 64×16 meshes, (b) 128×32 meshes, and (c) 256×64 meshes

Full size|PPT slide

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 eigenfrequency gradually increases after the fifth iteration, whereas that of the second and third eigenfrequencies decline with the increase in iterations.

Fig.9 Convergence history

Full size|PPT slide

Fig.10 Iteration history of the first three eigenfrequencies

Full size|PPT slide

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.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (Grant No. 51675197).

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://creativecommons.org/licenses/by/4.0/.

References

1 Suzuki K, Kikuchi N. A homogenization method for shape and topology optimization. Computer Methods in Applied Mechanics and Engineering, 1991, 93(3): 291–318

2 Sigmund O. A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 2001, 21(2): 120–127

3 Pedersen N L. Maximization of eigenvalues using topology optimization. Structural and Multidisciplinary Optimization, 2000, 20(1): 2–11 doi:10.1007/s001580050130

4 Du J, Olhoff N. Minimization of sound radiation from vibrating bi-material structures using topology optimization. Structural and Multidisciplinary Optimization, 2007, 33(4‒5): 305–321

5 Iga A, Nishiwaki S, Izui K, Topology optimization for thermal conductors considering design-dependent effects, including heat conduction and convection. International Journal of Heat and Mass Transfer, 2009, 52(11‒12): 2721–2732

6 Yamada T, Izui K, Nishiwaki S. A level set-based topology optimization method for maximizing thermal diffusivity in problems including design-dependent effects. Journal of Mechanical Design, 2011, 133(3): 031011

7 Bendsøe M P, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197–224

8 Bendsøe M P. Optimal shape design as a material distribution problem. Structural Optimization, 1989, 1(4): 193–202

9 Rozvany G I N, Zhou M, Birker T. Generalized shape optimization without homogenization. Structural Optimization, 1992, 4(3–4): 250–252

10 Bourdin B, Chambolle A. Design-dependent loads in topology optimization. ESAIM. Control, Optimisation and Calculus of Variations, 2003, 9: 19–48

11 Wang M Y, Zhou S. Phase field: A variational method for structural topology optimization. Computer Modeling in Engineering & Sciences, 2004, 6(6): 547–566

12 Xie Y M, Steven G P. A simple evolutionary procedure for structural optimization. Computers & Structures, 1993, 49(5): 885–896

13 Querin O M, Steven G P, Xie Y M. Evolutionary structural optimisation (ESO) using a bidirectional algorithm. Engineering Computations, 1998, 15(8): 1031–1048

14 Díaz A, Sigmund O. Checkerboard patterns in layout optimization. Structural Optimization, 1995, 10(1): 40–45

15 Jog C S, Haber R B. Stability of finite element models for distributed-parameter optimization and topology design. Computer Methods in Applied Mechanics and Engineering, 1996, 130(3‒4): 203–226

16 Sigmund O. Materials with prescribed constitutive parameters: An inverse homogenization problem. International Journal of Solids and Structures, 1994, 31(17): 2313–2329

17 Petersson J, Sigmund O. Slope constrained topology optimization. International Journal for Numerical Methods in Engineering, 1998, 41(8): 1417–1434

18 Sethian J A, Wiegmann A. Structural boundary design via level set and immersed interface methods. Journal of Computational Physics, 2000, 163(2): 489–528

19 Osher S J, Santosa F. Level set methods for optimization problems involving geometry and constraints: I. Frequencies of a two-density inhomogeneous drum. Journal of Computational Physics, 2001, 171(1): 272–288

20 Allaire G, Jouve F, Toader A M. A level-set method for shape optimization. Mathematical Rendering, 2002, 334(12): 1125–1130 (in French)

21 Wang M Y, Wang X, Guo D. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1‒2): 227–246

22 Allaire G, Jouve F, Toader A M. Structural optimization using sensitivity analysis and a level-set method. Journal of Computational Physics, 2004, 194(1): 363–393

23 Dijk N P, Langelaar M, Keulen F. Explicit level-set-based topology optimization using an exact Heaviside function and consistent sensitivity analysis. International Journal for Numerical Methods in Engineering, 2012, 91(1): 67–97

24 Wang S, Wang M Y. Radial basis functions and level set method for structural topology optimization. International Journal for Numerical Methods in Engineering, 2006, 65(12): 2060–2090

25 Luo Z, Tong L, Wang M Y, et al. Shape and topology optimization of compliant mechanisms using a parameterization level set method. Journal of Computational Physics, 2007, 227(1): 680–705

26 Gomes A A, Suleman A. Application of spectral level set methodology in topology optimization. Structural and Multidisciplinary Optimization, 2006, 31(6): 430–443

27 Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

28 Grebennikov A I. Isogeometric approximation of functions of one variable. USSR Computational Mathematics and Mathematical Physics, 1982, 22(6): 42–50 doi:10.1016/0041-5553(82)90095-7

29 Nguyen T, Jüttler B. Parameterization of contractible domains using sequences of harmonic maps. In: Proceedings of International Conference on Curves and Surfaces. Berlin: Springer, 2010, 501–514

30 Xu G, Mourrain B, Duvigneau R, Constructing analysis-suitable parameterization of computational domain from CAD boundary by variational harmonic method. Journal of Mathematical Physics, 2013, 252: 275–289

31 Xu G, Mourrain B, Duvigneau R, Optimal analysis-aware parameterization of computational domain in isogeometric analysis. In: Proceedings of International Conference on Geometric Modeling and Processing. Berlin: Springer, 2010, 236–254

32 Xu G, Mourrain B, Duvigneau R, Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications. Computer-Aided Design, 2013, 45(2): 395–404

33 Bazilevs Y, Calo V M, Cottrell J A, Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5‒8): 229–263

34 Wu Z, Huang Z, Liu Q, A local solution approach for adaptive hierarchical refinement in isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2015, 283: 1467–1492

35 Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Rotation free isogeometric thin shell analysis using PHT-splines. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47‒48): 3410–3424

36 Nguyen-Thanh N, Nguyen-Xuan H, Bordas S P A, Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids. Computer Methods in Applied Mechanics and Engineering, 2011, 200(21–22): 1892–1908

37 Atroshchenko E, Tomar S, Xu G, Weakening the tight coupling between geometry and simulation in isogeometric analysis: From sub- and super-geometric analysis to geometry-independent field approximation (GIFT). International Journal for Numerical Methods in Engineering, 2018, 114(10): 1131–1159

38 Speleers H, Manni C, Pelosi F, Isogeometric analysis with Powell-Sabin splines for advection-diffusion-reaction problems. Computer Methods in Applied Mechanics and Engineering, 2012, 221 ‒ 222: 132–148

39 Auricchio F, Da Veiga L B, Hughes T J R, Isogeometric collocation methods. Mathematical Models and Methods in Applied Sciences, 2010, 20(11): 2075–2107

40 Xu G, Li M, Mourrain B, Constructing IGA-suitable planar parameterization from complex CAD boundary by domain partition and global/local optimization. Computer Methods in Applied Mechanics and Engineering, 2018, 328: 175–200

41 Simpson R N, Bordas S P, Trevelyan J, A two-dimensional isogeometric boundary element method for elastostatic analysis. Computer Methods in Applied Mechanics and Engineering, 2012, 209–212: 87–100

42 Lian H, Simpson R N, Bordas S. Stress analysis without meshing: Isogeometric boundary-element method. Proceedings of the Institution of Civil Engineers: Engineering and Computational Mechanics, 2013, 166(2): 88–99

43 Peng X, Atroshchenko E, Kerfriden P, Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 151–185

44 Peng X, Atroshchenko E, Kerfriden P, Linear elastic fracture simulation directly from CAD: 2D NURBS-based implementation and role of tip enrichment. International Journal of Fracture, 2017, 204(1): 55–78

45 Simpson R N, Scott M A, Taus M, Acoustic isogeometric boundary element analysis. Computer Methods in Applied Mechanics and Engineering, 2014, 269: 265–290

46 Lian H, Kerfriden P, Bordas S P A. Shape optimization directly from CAD: An isogeometric boundary element approach using T-splines. Computer Methods in Applied Mechanics and Engineering, 2017, 317: 1–41

47 Lian H, Kerfriden P, Bordas S P. Implementation of regularized isogeometric boundary element methods for gradient-based shape optimization in two-dimensional linear elasticity. International Journal for Numerical Methods in Engineering, 2016, 106(12): 972–1017

48 Hughes T J R, Reali A, Sangalli G. Efficient quadrature for NURBS-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5‒8): 301–313

49 Reali A.An isogeometric analysis approach for the study of structural vibrations. Journal of Earthquake Engineering, 2006, 10(spec01): 1–30

50 Bazilevs Y, Calo V M, Zhang Y, Isogeometric fluid-structure interaction analysis with applications to arterial blood flow. Computational Mechanics, 2006, 38(4‒5): 310–322

51 Bazilevs Y, Calo V M, Hughes T J R, Isogeometric fluid-structure interaction: Theory, algorithms, and computations. Computational Mechanics, 2008, 43(1): 3–37

52 Xu G, Mourrain B, Duvigneau R, Parameterization of computational domain in isogeometric analysis: Methods and comparison. Computer Methods in Applied Mechanics and Engineering, 2011, 200(23‒24): 2021–2031

53 Wall W A, Frenzel M A, Cyron C. Isogeometric structural shape optimization. Computer Methods in Applied Mechanics and Engineering, 2008, 197(33–40): 2976–2988

54 Cho S, Ha S H. Isogeometric shape design optimization: Exact geometry and enhanced sensitivity. Structural and Multidisciplinary Optimization, 2009, 38(1): 53–70

55 Kiendl J, Bletzinger K U, Linhard J, Isogeometric shell analysis with Kirchhoff-Love elements. Computer Methods in Applied Mechanics and Engineering, 2009, 198(49‒52): 3902–3914

56 Wang Y, Benson D J. Isogeometric analysis for parameterized LSM-based structural topology optimization. Computational Mechanics, 2016, 57(1): 19–35

57 Buffa A, Sangalli G, Vázquez R. Isogeometric analysis in electromagnetics: B-splines approximation. Computer Methods in Applied Mechanics and Engineering, 2010, 199(17–20): 1143–1152

58 Lee S, Kwak B M, Kim I Y. Smooth boundary topology optimization using B-spline and hole generation. International Journal of CAD/CAM, 2007, 7(1): 11–20

59 Seo Y D, Kim H J, Youn S K. Isogeometric topology optimization using trimmed spline surfaces. Computer Methods in Applied Mechanics and Engineering, 2010, 199(49–52): 3270–3296

60 Hassani B, Khanzadi M, Tavakkoli S M. An isogeometrical approach to structural topology optimization by optimality criteria. Structural and Multidisciplinary Optimization, 2012, 45(2): 223–233

61 Tavakkoli S M, Hassani B, Ghasemnejad H. Isogeometric topology optimization of structures by using MMA. International Journal of Optimization in Civil Engineering, 2013, 3(2): 313–326

62 Dedè L, Borden M J, Hughes T J R. Isogeometric analysis for topology optimization with a phase field model. Archives of Computational Methods in Engineering, 2012, 19(3): 427–465

63 Díaaz A R, Kikuchi N. Solutions to shape and topology eigenvalue optimization problems using a homogenization method. International Journal for Numerical Methods in Engineering, 1992, 35(7): 1487–1502

64 Tenek L H, Hagiwara I. Static and vibrational shape and topology optimization using homogenization and mathematical programming. Computer Methods in Applied Mechanics and Engineering, 1993, 109(1–2): 143–154

65 Xie Y M, Steven G P. Evolutionary structural optimization for dynamic problems. Computers & Structures, 1996, 58(6): 1067–1073

66 Ma Z D, Kikuchi N, Hagiwara I. Structural topology and shape optimization for a frequency response problem. Computational Mechanics, 1993, 13(3): 157–174

67 Ma Z D, Cheng H C, Kikuchi N. Structural design for obtaining desired eigenfrequencies by using the topology and shape optimization method. Computing Systems in Engineering, 1994, 5(1): 77–89

68 Kim T S, Kim Y. Mac-based mode-tracking in structural topology optimization. Computers & Structures, 2000, 74(3): 375–383

69 Sigmund O, Jensen J S. Systematic design of phononic band-gap materials and structures by topology optimization. Philosophical Transactions: Mathematical, Physical and Engineering Sciences, 2003, 361(1806): 1001–1019

70 Kosaka I, Swan C C. A symmetry reduction method for continuum structural topology optimization. Computers & Structures, 1999, 70(1): 47–61

71 Neves M M, Rodrigues H, Guedes J M. Generalized topology design of structures with a buckling load criterion. Structural Optimization, 1995, 10(2): 71–78

72 Allaire G, Jouve F. A level-set method for vibration and multiple loads structural optimization. Computer Methods in Applied Mechanics and Engineering, 2005, 194(30–33): 3269–3290

73 Shu L, Wang M Y, Fang Z, Level set based structural topology optimization for minimizing frequency response. Journal of Sound and Vibration, 2011, 330(24): 5820–5834

74 Xia Q, Shi T, Wang M Y. A level set based shape and topology optimization method for maximizing the simple or repeated first eigenvalue of structure vibration. Structural and Multidisciplinary Optimization, 2011, 43(4): 473–485

75 Liu T, Li B, Wang S, Eigenvalue topology optimization of structures using a parameterized level set method. Structural and Multidisciplinary Optimization, 2014, 50(4): 573–591

76 Hassani B, Hinton E. A review of homogenization and topology optimization III—Topology optimization using optimality criteria. Computers & Structures, 1998, 69(6): 739–756

77 Nguyen V P, Anitescu C, Bordas S P A, Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117: 89–116