Parametric Structural Optimization of 2D Complex Shape Based on
Isogeometric Analysis

The geometric model and the analysis model can be unified together through the isogeometric analysis method, which has potential to achieve seamless integration of CAD and CAE. Parametric design is a mainstream and successful method in CAD field. This method is not continued in simulation and optimization stage because of the model conversion in conventional optimization method based on the finite element analysis. So integration of the parametric modeling and the structural optimization by using isogeometric analysis is a natural and interesting issue. This paper proposed a method to realize a structural optimization of parametric complex shapes by using isogeometric analysis. By the given feature curves and the constraints, a feature frame model is built. Based on the feature frame model, a parametric representation of complex shape is obtained. After adding some auxiliary curves, the feature frame model is divided into many box-like patches in three dimension or four-sided patches in two dimension. These patches are built into parametric patches by using volume interpolation methods such as Coons method. Based on the parametric patches, isogeometic analysis is applied. Thus, the relationships are constructed among the size parameters, the control points and the physical performance parameters. Then the sensitivity matrix could be derived based on the relationships. The size optimization is carried out in the first stage by taking the size parameters as variables. Based on the result of size optimization, shape optimization with the constraints of stress is carried out in the second stage by taking the control points as variables. Serval planar complex shapes are taken as example to verify our method. The results verify that the parametric modeling and structural optimization can be united together without model conversion. Benefit from this, the optimization design can be executed as a dark box operation without considering the concrete modeling and analysis by input of the sizes, constraints and loads.


Introduction
In the stage of product optimization design, phase, the initial CAD model can be optimized based on the results obtained during the simulation analysis phase. The optimal CAD model improves the performance of the CAD products [1][2][3][4]. Structural optimization has always been an important research area because it usually attempts to integrate structural analysis, optimization algorithms and geometric modeling into an automated design. The main purpose of structural optimization is to improve structural characteristics, such as reducing the stress concentration and weight, and increasing the stiffness by changing the structure boundary geometry. Optimization is divided into three stages: size optimization, shape optimization and topology optimization [5]. Topological optimization is to determine the best product prototype by changing the topological structure according to the best force transmission path in a given material design domain. Shape optimization takes the boundary shape of the structure as the optimization object. It takes the node coordinates as the design variables. Finally it improve the performance of product by changing the boundary shape [6][7]. Size optimization is the optimization of discrete dimensional variables to determine the final accurate dimensional model of the product.
In traditional structural optimization, the design model is converted into the finite element model in the process of applying FEA method. Transformation of design model, analysis model and optimization model leads to inaccuracy and tediousness. In order to describe the boundaries to be optimized, earlier method took node coordinates as design variables [8]. Taking the node coordinates of finite element model has the following disadvantages [11]. Firstly, a large number of design variables are needed to describe the shape of the boundary. Secondly, the independent changes of each node in the optimization iteration can easily lead to grid distortion and distortion. Thirdly, the model of design, analysis and optimization are not unified, so that data transmission can only be unidirectional.
Isogeometric optimization method has received more and more attentions recently [9][10]. Benefiting from the superiority of splines in shape expression, shape optimization based on isogeometric Analysis has achieved great success [12][13]. However, the current isogeometric analysis and optimization methods encounter difficulties when facing the complex shapes or models because of the complexity existing in the stage of construction and analysis.
For the 2D complex shapes, we do not directly treat design model but its feature frame model as the optimization object. In this way, the design variables are greatly reduced. The optimization object is closer to the designer's intent since the feature sizes of the product are mainly mainly the designer's concern. The feature frame model is used to control the shape of overall model and detail. The designer is allowed to edit the product model by a few feature parameters. Thus, the parametric modeling based on feature sizes and the optimization based on the parametric model are unified together. The contributions of this paper as follows: 1) A framework supporting parametric modeling, isogeometric analysis and optimization is proposed. The framework supports to construct the relations existing among the design size variables, control points and physical performance parameters, which will help to derive the sensitivity matrix for optimization.
2) Structural optimization including shape and size optimization is realized with the help of sensitivity matrix to reach the optimization of global and local performance by taking the sizes and control points as the design variables.

Related Works
The isogeometric analysis method based on volume parametric model was first proposed by Hughes in 2005 [14]. Compared with other methods, he application of isogeometric analysis to product design has more advantages. The consistency of the design model and the analysis model are ensured, and the seamless integration between the models is achieved. It will obtain more precise result than the traditional finite element method which uses geometry approximation. The current research topics of isogeometric analysis include: 1) the basic theory of isogeometric analysis, including accuracy calculation, convergence, continuity and singularity, and so on [15][16]; 2) improvement of Isogeometric analysis method, such as boundary element method, collocation method, etc. [17][18]; 3) extension of isogeometric analysis based on various spline basis functions [19][20]; 4) application of geometric analysis for various practical applications [21][22]. Before widely exploited in engineering filed, isogeometric analysis of complex shapes is an important and necessary issue. Due to the inherent attribute of volume parameterization presented with tensor product, parameterization and analysis of complex shapes is a bottleneck problem. There have been some works that try to solve this problem, such as wear coupling and trimmed multipatch merging or tearing [35][36].
We call the optimization method based on isogeometric analysis as isogeometric optimization, which has been carried out in the field of size optimization, shape optimization and topology optimization. Great research progress has been made in shape optimization based on isogeometric optimization benefiting from the superiority of splines in expressing shape [23][24][25][26][27][28]. There have some tentative studies on topology optimization, but the current research has not made much progress due to the limitations of the spline function in expressing topology [29][30][31][32]. There has few studies on size optimization based on isogeometric optimization [33,[42][43][44][45][46], because it is not an intuitive way to present the curved shapes expressed as splines into sizes which are mostly horizontal, vertical, and angular. At the same time, the shapes expressed by splines are mostly free shapes, and parametric modeling of these shapes is not a trivial work. Since isogeometric optimization uses NURBS splines as basis function, which improves the continuity of the computational domain, sensitivity analysis becomes an important issue in isogeometric optimization problems [34].

Overview of Our Algorithm
In our method, a feature frame model composed of feature curves is firstly constructed. The feature curves can be acquired through given sizes or extracting from a point cloud model. A parametric feature frame model can be driven to modify by adding some constraints. Then the feature frame model is divided into many box-like subdomains by adding some curves. Volume interpolation such as Coons method and the optimization methods applied on the inner control points are used to create volumetric solid. By using the IGA method, displacement, stress and other physical parameters of each control point can be obtained after applying boundary conditions and constraints on to the analysis model [37]. The relationships between physical parameters and characteristic parameters are established. According to the above established relationships, the sensitivity matrix can be obtained by the derivatives of physical parameters with respect to feature parameters. Through the optimization algorithm and the sensitivity matrix, the feature parameters and the corresponding feature frame model are updated until the given iteration termination conditions are satisfied. The work process of product optimization design is shown in Fig. 1.

.1 Domain Composition
The first step of the isogeometric analysis is creating the parametric computational domain. Computational domain with complex boundary or topology is supposed to be divided into multi parametric quadrilateral computational domains. This problem is similar to the domain decomposition method by dividing the solving region into a number of sub-systems or sub-regions according to some certain rules. The methods can be divided into two categories, non-overlapping subdomain decomposition and overlapping domain decomposition, as shown in Fig. 2. In the figure, (a) denotes the Noneoverlapping domain decomposition, (b) denotes the overlapping domain decomposition.
Assuming the computational domain is divided into a series of subdomain { i }, the overall domain with overlapping domains can be expressed as follows: And the non-overlapping domain is,

Parametric Modeling Method
In this paper, all the subdomains denote NURBS surfaces, so the objective of modeling is to construct these subdomains. A feature frame model is consisted of feature points, feature curves or feature surfaces. The feature points mainly emphasize the morphological location of details but they are inadequate to describe the whole shape. The shape of the surface has a specific function, which is synthesized by the feature curves. Though the whole content of the shape can be expressed, it is too over complicated and not easy to be operated. The feature curves contain important modeling information that can be linked to the feature points and the feature surfaces. They have no cumbersome operations caused by excessive number of points. They also have no redundancy of surfaces, so the processing of the feature curves is more superior to that of the feature points and feature surfaces.
Here, feature frame model is only composed of feature curves. The feature curves can be generated by curve parameterization method by giving the sizes. The feature curves can also be extracted from point cloud models. All the feature curves are linked by all kinds of geometric constraints. The whole process is depicted in Fig. 3.
When constructing the feature curve frame, the following basic principles should be considered when designing the characteristic curve: 1) Feature curves have semantic meaning in order to reflect the shape, size and position of the design object; 2) The feature curves are easy to edit in order to control the shape or size during the interactive operation; 3) The feature curves are simple in order to reduce the difficulty of modeling. Following the three principles, the curves of profiles and sections are selected as feature curves, as well as the curves with design significance. The feature curves are stratified into geometric layer, constraint layer and semantic layer. The geometric layer is generally composed of control points, knot vectors or interpolation points. The layer constitutes the lowest level of the feature curve. The constraint layer is composed of the constraints of geometry, dimension and topology, which are generated during the process of generation and modification of feature curves. The semantic layer describes the semantic parameters and functions of the feature curves. The semantic parameters are defined for the control feature curves according to the actual requirements, such as height, length, width, girth.
The mathematic form of the feature curve can be defined as F = {D; R; S}. F denotes the feature curve. D is the set of control points. R is the constraint type illustrated in Tab. 1. S is the semantic parameter of the feature curve.
After generation, the feature frame can be modified to get a new model. There are two ways to realize the modifications: size-driving and sketch-driving. Size-driving means to drive the target feature curve by modifying its sizes. Sketch-driving refers to driving the target feature curve by sketch interaction. During the process of driving, it is necessary to use appropriate constraint solving methods to maintain the constraints between the characteristic curves. Resultantly, a new feature frame can be generated.  Based on the feature frame, some key points are generated according to the given sizes and the parametric relations among them. On the basis of these key points, all the boundary curves of the foursided subdomain can be obtained. If we only discuss the case in two dimension, all the subdomains denote the NURBS surfaces. The feature frame is divided into multi subdomains by adding some curves. Then all the control points of each subdomain can be acquired using the interpolation method such as Coons method. Accordingly, the parametric representation of the whole model can be finished. Thus, the relation between the size parameters and the control points can be built.
If we discuss the case in three dimension, the thing gets complicated. The first difficulty exists in the segmentation of the feature frame. Dividing the feature frame into many cube-like subdomains is far from a simple task. The second difficulty exists in the volume parameterization. Construction of volume parameterization for each subdomain is the core task of the research of parameterization and has attracted many scholars' interests. Here we only focus on the two dimensional case.

Basic Formulation of IGA
If given a set of knot vectors ½n 1 ; n 2 ; n 3 ; ……; n nþpþ1 , the B-spline curve can be defined as [34], where P i is the set of control points. N i;p is the B-spline basis function with the order p. Similarly, the B-spline surface can be defined as, B-spline solid is defined as, In this paper, the linear elastic problem is taken as an example for study. The principle of minimum potential energy is applied to obtain the equation of equilibrium according to relative mechanics principles and variation principles.
where u is the array of unknowns for the control points. f is the loads applied on the model. K is the stiffness matrix presented as follows, To calculate the stiffness matrix K, three domains are defined. They are the physical domain and its subdomain or element e , the parameter domain and its subdomain or element e , and the parent domaiñ , respectively. The relative mapping relations are defined as follows.f e : ! e is the mapping from the parent domain to the parameter subdomain. S : e ! e is the mapping from the parameter subdomain to the physical subdomain. Therefore, the function is S f e for mapping X e : ! e from the parent domain to the physical subdomain. For the three dimensional case, an element in the parameter space is Any point in parameter space n; g; f ð Þ can be expressed as follows with the pointñ;g;f by the mappingf e : ! e n g f The Jacobian transformation matrix of mapping S : e ! e is J 1 ¼ @x @n @x @g @x @f @y @n @y @g @y @f @z @n @z @g @z @f The Jacobian transformation determinant for the transformation from the parameter coordinate to the parent unit coordinate is The Jacobian determinant of mapping X e : ! e is The integral form of the stiffness matrix in the whole physical space is as follows.
where n e is the number of the elements or subdomains. The item B is defined as follows.
R s is defined as follows, Then, the matrixes of strain and stress can be computed as follows, where, D is a constant matrix.
For the plane stress problem, D if defined as follows, For the plane strain problem, D if defined as follows. Where E represents the Elastic Modulus, l represents Poisson's ratio.

IGA for Multi-Patches
One model can be divided into many subdomains. The whole interface curve is shared by two subdomains. The equation of equilibrium for each subdomain is integrated into the global equation of equilibrium. By solving the equation, the analysis of the whole structure is realized. If the two subdomains don't share the same interface curve, a connectivity matrix is imported to establish a correspondence between the two subdomains. The following derivation is in accordance with the subdomains with the same interface curve [38].
where u m denotes the unknown displacement field. r m denotes the stress tensor. b m denotes volume force. displacement of boundary S Ã . n 1 and n 2 are the outward unit normal of 1 and 2 respectively. According to the method [38], the stiffness matrix in multiple patches is expressed as: The relevant stiffness matrix is given by: R m i;j n; g ð Þ is a binary NURBS basis function supporting the mth patcheŝ The force vector f is expressed as Substitute Eqs. (19) and (24) into Eq. (4), we can get the equation of equilibrium for the whole domain. By solving the equation, we can get the matrix of displacement u for each patch. If the number of subdomain is more than two, the relevant matrices of stiffness and load are listed here.
where ph denotes the index of the subdomain, nph is the total number of all the subdomains.

Optimization Principle
For the linear elastic problem, the single object of structure optimization is usually to minimize the volume or the area of the structure [39][40][41]. There are two optimization function. The objective function w 1 is represented for the flexibility of the model by the formula w 1 b i ½ ð Þ ¼ f T u. The objective function w 2 is represented for the von Mises stress by the formula w 2 a j Â Ã À Á ¼ r von . The design variables b i is the feature sizes, and a j is the coordinate component of control points to be optimized. In two dimensional case, there are two coordinate components denoting x and y. All the coordinate components of the optimized control points are arranged into a linear set a j È É .

Derivation of Sensitivity Matrix
The derivative of objective function in the formula (26) with respect to a i needs to be calculated.
Since the load vector is independent of b i , so df T db i ¼ 0. Thus the first item of the formula (27) is eliminated and the following formula comes out.
To calculate the formula (28), the item du db i is required to compute. From the equation of equilibrium Ku ¼ f, we can get the formula (29) by deriving the equation on both sides of it. The formula (30) is also obtained by transposing the item.
Taking the formula (30) into the formula (28), we can get the following formula since the matrix K is symmetric.
According to the chain rule of derivation, we can get, According to the definition of a j , the item da j db i can be written as d P j :x P j :y ½ db i . Since we have constructed the relation between the control points and the feature sizes in Section 4.2, the item da j db i can be derived from the parametric representation of each subdomain.
For another object function, If defining r ¼ ½r xx r yy r xx À r yy ffiffi ffi 6 p r xy , we can get r von ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r r T =2 q and then get the sensitivity matrix according to the partial derivative of r von with respect to a j .
If introducing a matrix A, r can also be written as the following matrix.
If defining C ¼ A T A, we can get following formula since r ¼ ½r xx r yy r xy ¼ DBu.
Taking the formula into the Eq. (33), we can get Calculating both dw 1 da j and dw 2 da j requires to solve dK da j .
The problem is further deduced as @K b e @a j ; @K n e @a j ; @K s e @a j since K b ; K n ; K s are all assembled by K b e ; K n e ; K s e . The design variable a j is not included in K s e , so @K s e @a j will vanishes.
@K n e @a j ¼ @K n;11 e @a j @K n;12 e @a j @K n;21 e @a j @K n;22 e @a j 2 6 6 4 3 When calculating @K n e @a j ; @K s e @a j , we need to integrate along the curve S. Under the premise of ensuring the continuity of the sub-patches, the Gaussian integral of the linear element can be used to solve the interface integral, that is to calculate the sum of the integrand which is transformed by the curve S Ã . @k n;11 In the formula, ngp denotes the number of Gaussian integration points, and w i is the weight.
In the above formula, the items of @B m e @a j and @ J m j j @a j should be derived. In order to calculate these two items, two matrices are defined in the first place.
Derivative on both sides of Eq. (11): We can get, If the formula (3) is written into a matrix formation, we can get the following expression for one calculating element.
where R m e ¼ R 1 R 2 Á Á Á R s Á Á Á R n e ½ T and R s have been defined in the formula (13). P m e is the coordinate matrix of all the control points in a calculating element.
On derivate of the formula (47) on both sides, we can get the formula (49).
The second item after the equal-sign of the formula (49) is equal to zero. According to the definition of T , the following formula can be derived.
On the derivate of the formula (46) and the formula (50) on both sides respectively, we can get the formula (51).
According to the definitions of the matrices B and T, which are defined in the formula (11) and (43) respectively, we can get @B m e @a j by rearranging the non-zero items of @T m e @a j .
Next, @ J m e @a j is calculated. To avoid confusion with the parameter of J 2 , the parameter of a j is written into the form a p since there is no relationship between them.
. . . @x ne @a j @y ne @a j According to formula (13), formula (53) is obtained: And we can get the following equation based on Eq. (15): Taking the derivative of both sides of formula (53) to get formula (55): where tr (Z) denotes the sum of the diagonal elements of Z. The load applied on the model is independent of the design variable, so the derivative of the element load array F to the design variables is 0, i.e., @f e @a j = 0. By the above Eqs. (51) and (55) Bringing these two data into Eqs. (36), (37) and (39), we are able to respectively calculate @k b e @a j which is the derivative of element volume matrix to design variable, @k n e @a j and @k s e @a j are the derivatives of the element coupling matrix to design variables.
Optimization algorithm is the key technology that influences the quality and efficiency of shape optimization. This paper employs the method of moving asymptotes called MMA [36]. The MMA algorithm is used to update the design variables in each iterative step by using the sensitivity matrix.

Examples and Discussions
We select two planar complex models to verify our method. The first is a connecting-rod and the second is a load-bearing part. In the two examples, the elasticity modulus E is equal to 1e5 and the Poisson's ratio ν is equal to 0.3. Since the constraints of each model are not same, a general feature frame model cannot be given. The size optimization is performed directly on the feature frame model with the objective of flexibility by taking the size parameters as design variables. Based on the result of the size optimization, shape optimization is realized with the objective of equivalent stress by taking the control points as the design variables.

Example 1: Optimization of Connecting-Rod
The connecting-rod model has 9 size parameters. The feature frame model is constructed according to these size parameters, as shown in Fig. 5a. After adding some auxiliary curves, the feature frame model is subdivide into 12 sub-patches. Each sub-patch is interpolated as NURBS surfaces, as shown in Fig. 5b. The degree of the 12 sub-patches is equal to 2 in both parameter directions of U and V, and the knot vectors in the two parameter directions are {0,0,0,1,1,1}. Giving all the parametric equations is a tedious work, so we select some key points and give the equations of them. All the parameters in the equations are denoted in Fig. 5a. Then, all the control points of all the subdomains can be obtained using the Coons method on the basis of these key points. The constraints listed in the Tab. 1 are applied on the parametric model, such as the geometric constraint F1 and F11, the topologic constraint F3 and F14, the size constraint F4. Since parametric design for 2D model is not a complex work, we will not give the detail here.
The two distributed load P1 = 100 N and P2 = 100 N are applied to the center of the right hole. Full constraint is applied to the left hole. The objective of optimization is to minimize the global flexibility and local equivalent stress of the structure. Optimization functions consists of two functions. The first is . During the optimization process, the area of the model is less than a given value and the rules of the Elasticity.
A global flexibility optimization is performed, then an optimization of the local equivalent stress is applied. As shown in Figs. 6a and 6b, in the case of V ≤ V 0 = 371, the flexibility is 610 and 507 before and after optimization. As shown in Fig. 6c, the maximum equivalent stress before optimization is at the slot. The slot is in a state of stress concentration with the value 4300. While in Fig. 6e, the maximum equivalent stress is cut down to 3700. The optimization efficiency reaches 16%, with the area decreasing and the stress concentration disappearing. The iteration process of objective functions is shown in Fig. 7. The size parameters before and after optimization are shown in Tab. 2.

Example 2: Optimization of Load-Bearing Part
The second example is a load-bearing part from a planar flexible mechanism. As shown in Fig. 8a, the constraints and loads are imposed the parametric model, and the layout of the sub-patches that amount to 34 is displayed in Fig. 8b. The degree of these sub-patches is 2 in both U and V parameter direction. The knot vectors in the two parameter directions are {0, 0, 0, 1, 1, 1}. Similar with the first example, we select some key points to give the formulation of them and then get the parametric representation of all the subdomains.
The constraints is also applied on the model, such as F11, F14 and F4. Fixed displacement constraints are applied in three places depicted in the figure. The applied point load F is equate to 1000N. The target of optimization is minimization of the flexibility w 1 P i ½ ð Þ ¼ f T u and the equivalent stress with the constraint of the area and the rules of the Elasticity.
From the result presented in the Fig. 9, the flexibility is 4721 and 3800 before and after optimization with the area constraint V V 0 = 862. The optimized rate was 19%. The optimization of equivalent stress follows the size optimization. As seen in the Fig. 9c, the ring groove notch has the most severe stress concentration with the value 32483. In Fig. 9e, the maximum equivalent stress is reduced after optimization with the

Conclusions
Our method has serval advantages. The method combines the parametric design method with the isogeometric analysis. By constructing the relationship chain between design size parameters, control points and physical performance parameters, the feature frame model is used to control the overall and details of the shape. This allows the user to edit the product model through a few feature size parameters. This also avoids the complex modeling operations caused by traditional shape optimization using finite element nodes as design variables. Thus the parametric modeling and optimization are unified together.
Benefit from the parametric design, our method lightens the burden on the designer and transfer the designer's special attention from mesh operation during the stage of analysis and optimization to the creative design itself. The designers will put their attentions on optimization of the key part sizes, without taking into account nonsignificant sizes. It will be easier and more accurate to capture the designer's design intent, since it is a very important thing for computer-aided design software or computer-aided engineering software.
Thirdly, the computation efficiency is significantly improved for the variables of design and optimization are greatly reduced. In the design process, the interactions among the model of design, analysis and optimization are more convenient. This greatly reduces the product research period. Due to the B-spline basis function, the sensitivity matrix is easily acquired because the entire model has complete continuity.
It must be mentioned that the feature frame models presented in this paper is simple, so the constraints between the feature curves are single and simple as well. Meanwhile, the Coons does not always work well if the subdomains are not divided well when the model is complex. Henceforth, the next work is to construct the feature frame model with complex feature elements and constraints. More volume parameterization methods can also be applied onto the model. Thus, our method will be able to treat the more complex optimization problem.
Acknowledgement: Thank all the anonymous reviewers for their valuable comments and suggestions.
Funding Statement: This work is supported by the National Nature Science Foundation of China (No. 51475309).

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.