#### Abstract

A fast shape optimization strategy for free form shell structure design with structural dynamics criteria is proposed in this paper. The structures are modelled with Non-Uniform Rational B-Spline based isogeometric Kirchhoff-Love shell elements. The substitution of the traditional finite elements not only makes the mesh model geometrically exact but also avoids the laborious mesh regeneration during the design update. As for the structural response evaluation, the modal synthesis method is adopted to avoid a repeated evaluation of some substructures where there are no designed variables attached; thus, the model reanalysis is speeded up. A bottom-up strategy for the analytical design sensitivity evaluation is also proposed here; the element-level analytical sensitivity with respect to the inherent shape parameters is firstly calculated from which the design sensitivity is then extracted with the help of a sensitivity map. Finally, gradient based algorithm is used to solve the optimization problem. Several examples show that our approach is flexible and efficient for fast free form shell structure optimization.

#### 1. Introduction

Isogeometric analysis (IGA) [1] is a promising method which can bridge the gap between computer aided design (CAD) and computer aided analysis (CAE). It is a special finite element method (FEM) with an exact geometrical mapping. Unlike traditional FEM where the analysis model (meshes) is a faceted approximation of the CAD model, the IGA analysis model is geometrically exact. The mesh refinements do not affect the physical geometry of the structure [2]. Normally, wherever the traditional FEM can be used, IGA can also be used, e.g., various mechanics problems [1, 3–6], electromagnetic problems [7, 8], and heat transfer problem [9]. The IGA solutions are more efficient and robust for these problems compared with the traditional FEM solutions [10].

In engineering structural optimization, model reanalysis is an indispensable step. It greatly affects the overall efficiency. In traditional FEM based optimization framework, the fixing of CAD model and creating FEM model account for 80% of overall analysis time cost [11]. With IGA analysis, large saving can be expected in this part, since the mesh refinements no longer need tracing back to the CAD model and are easy to implement [2].

IGA provides an explicit mapping between the parametric domain and the physical shape; the element-level analytical sensitivity of the structural equation ingredients, e.g., stiffness matrix and mass matrix, can thus be obtained [12–15]. The mesh creation algorithm also establishes mappings between the different meshes of the same model; the sensitivity of the fine meshed analysis model can be transformed into the coarse meshed model [15]. This allows an exact model analysis as well as a proper design variable definition. In our work, the design variables could be defined as combinations of the bottom mesh-level nodes (control points in IGA) movement. Since any modification of the design variable leads to a model update, it can be finally achieved by a joint movement of a set of nodes. The design sensitivities can thus be obtained by retracing with this map. IGA analysis has been successfully used in shape optimization [16, 17] and topology optimization [18–20] and shown its superiority over the traditional FEM based optimization. The geometrically exact modelling, the high quality of the response evaluation, and the efficient sensitivity analysis constitute the main reasons to choose IGA based optimization.

The modal synthesis method is a traditional method used for analysing the built-up structures concerning dynamics aspects. It is also termed as component modal analysis or branch mode analysis. The full structure is partitioned into some substructures; the substructure control equations are simplified by transforming the physical degrees of freedom into mode coordinates and then coupling with each other according to the interface compatibility. With modal synthesis, the dimensions of the problem will be reduced. In the present work, the Craig-Bampton [21] method is used. A deeper insight into the modal synthesis method can be obtained from the review papers [22–24]. For some optimization problems concerning dynamics aspects, the design variables might only relate to several substructures; it is less computationally costly to only reanalyze the corresponding substructures. Modal synthesis methods are proper to be utilized in such circumstance.

In this paper, we give a brief introduction on the linear IGA Kirchhoff-Love shell and the modal synthesis method in Section 2. In Section 3, we present our bottom-up strategy for the sensitivity evaluation; it enables efficient sensitivity calculation for various design variables. The sensitivity mapping concept and its implementation are also detailed here. In Section 4, we give two examples of shape optimization for natural frequency problem to demonstrate the effectiveness of our approach. The gradient based solver is adopted to make full use of the sensitivity information for fast convergence. In Section 5, we give a short conclusion.

#### 2. Preknowledge

##### 2.1. Isogeometric Kirchhoff-Love Shell Formulation

The linear Kirchhoff-Love shell formulation is adopted here. It is a linearization of the version proposed by Kiendl [25]. Such shell model does not account for the transversal shear deformation; thus, it is only applicable for thin shells. On the side, such simplification makes the node only consist of three translational degrees of freedom, reducing the dimensions of the problem for the same mesh size. The initial shell body is expressed as

The deformed shell body is

Here, and denote the shell middle surface and the shell normal before deformation and and for the deformed one. The notation denotes and is simply noted as hereinafter. The scripts noted by the Latin letters , belong to , and the Greek ones .

The displacement is

The components of the Lagrange-Green strain tensor are obtained with a constant part which indicates the membrane effect and a linear part which denotes the curvature change as follows:

For linear analysis, the deformed and unreformed configurations are approximately the same and the higher order item in the strain formulation is omitted. The linear membrane strain is obtained as

The bending part is

According to the principle of virtual work, the following equation can be obtained:

Here, , . They are effective stress resultants. indicates the shell middle surface, is the shell middle surface load density, and is the traction on the boundary . is the material tensor component [26] detailed as

and are the material Young’s modulus and the Poisson ratio, respectively. In IGA analysis, the shell middle surfaces ( and ) are expressed by NURBS [1, 27], so is the displacement field . The structure control equation for static analysis is

The structure control equation for modal analysis is

denotes the mass matrix, it is obtained by condensing mass to the shell middle surface and treating them as inertial force.

##### 2.2. Modal Synthesis

Modal synthesis method is traditionally used for the analysis of built-up structures concerning dynamics aspects. With the help of these methods, the control equation of large structures consisted by multiple parts can be simplified to reduce the problem dimensions. With a linear map constructed according to different modal synthesis approaches, the system control equation (12) can be simplified as

The transformation matrix is a product of and , they are used for separated substructure reduction and the reduced substructures coupling, respectively. There are many modal synthesis methods available which can be classified into three categories, the fixed interface method [21], the free interface method [28], and the mixed interface method [29, 30]. They differ from each other in the construction of and . We adopt the most often used fixed interface Craig-Bampton method [21] in the present work.

The discrete control equation for a substructure with harmonic loads is

Here, the footnotes and indicate the inner control points and the interface control points, respectively. The Craig-Bampton method adopts the selected fixed interface modes and the constraints modes to simplify system equation (12). The transformation matrix for substructure is

Here, denotes the selected fixed interface mode, . indicates the th eigenvector of substructure with all interface degrees of freedom constrained. Only the first eigenvectors are kept. denotes the constraint mode which is a set of displacement vectors obtained by sequentially applying a unite displacement on each interface degree of freedom and keeping others fully constrained.

The condensed substructure equation now can be obtained as

The simplified substructure FEM model can be viewed as a macroelement which only possesses interface degrees of freedom and the inner hanging nodes’ degrees of freedom. The inner nodes’ dynamics is now approximately considered by the interface degrees of freedom. When applying Craig-Bampton method to static problem, it works the same as the static condensation method (but with more degrees of freedom) which is usually used for the creation of super element. After model reduction, the interface degrees of freedom are kept. A connection between the independent macroelements should be established. If the substructures are connected, we can assemble the macroelements in the same way as the FEM method. If they are connected, it can be achieved with the bending strip method [31]. The correctness and effectiveness of this strategy for the prediction of structural dynamic properties of built-up shell structures have been testified in paper [32].

#### 3. Shape Optimization

Structural shape optimization is to adjust the shape of a continuous material region to meet specific criteria under structural control equation constraint and other user specified constraints. A general shape optimization problem can be formulated as

The first equality constraint is the structural control equation which is usually a partial differential equation (PDE) and represents the designer imposed other constraints. A general description of such type of optimization problem can be found in literature [33]. There are mainly two subcategories of shape optimization problem, the boundary optimization and the region optimization. For the first one, the design variables are the boundary curves by which the two-dimensional design area is bounded. Such topic has been intensively studied [34–36]. For the second one, the design variable is the region itself; it is the free form surface design problem [15]. The key concern is the parametrization of the design domain. In IGA analysis, both the shape and the approximate solution are expressed by NURBS; consequently, the design can be reformulated into the discrete form

Here, indicates the spatial coordinate of the th control point of the shape to be optimized, denotes the corresponding NURBS basis function, it is nonnegative compactly supported. As a result, all items of the stiffness and mass matrices are positive. and denote the orders of the basis function in two parametric directions, respectively.

The discrete structural control equation is

##### 3.1. Element-Level Sensitivity Analysis

Until now, the design variable has not been defined. However, in the above formulations there exists inherent shape design variables; they are the spatial coordinates of the control points of the IGA meshed model. Figure 1 shows an example where the th control point of the parametric surface is taken as design variable; it indicates an optimization of the coordinate component of that control point. The design sensitivity can be obtained as the perturbation of design variable tends to zero.

For the structural optimization problem, the sensitivity of the control equation ingredients, e.g., the mass and the stiffness matrix, with respect to the inherent shape parameters should be evaluated. Wall et al. [12] and Qian [13] gave a general description on this topic where they use IGA plane stress elements to model the structures. Nagy et al. [15] adopted the IGA Kirchhoff-Love shell elements and the analytical sensitivity. The similar work can also be found in literature [34, 35]. A thorough study on the sensitivity analysis under the IGA frame is given by [37]. Despite the large variety of sensitivity derivation, the key of its calculation is the chain differentiation of the structural ingredients.

Depending on the types of the elements adopted, the evaluation of (27) and (28) might be quite involved, since the variation of the shape parameters will change both the integrand and the integral domain. Kiendl et al. [38] used numerical differentiation to approximately evaluate them. It is less accurate and will be slow when there are many design variables. Here, we compute them analytically. The structural ingredients and their sensitivities are computed at the element level simultaneously with a single part of code. It is more efficient, since they use many intermediate variables in common.

##### 3.2. Sensitivity Mapping

In general shape design, although the initial CAD model is usually expressed with NURBS, the design variables which are defined on the CAD model may be of great diversity; for example, they could be the radius of a cylinder or some concerned control points through which the designer could manipulate the shape in their way. The variation of these features will guide the shape change and our ultimate goal is to find their sensitivities to support general fast shape optimization.

In traditional FEM based optimization, the design variation cannot directly drive the CAE model update; it should firstly update the CAD model and then remesh it to obtain the structural response; it is not a trivial work. Alternatively, some researchers studied the parameter guided mesh morphing technique, e.g., FFD, to directly update the meshes [39–41]. The key step is the construction of a design parameter relevant parametric domain of the meshes. It is not so easy, and the more severe problem is that it is difficult to maintain the mechanical features during the mesh deformation.

In IGA analysis, the CAD models and the meshes have the same parametric domain. The mesh creation only adds more control points without changing the model parametrization. Denote the CAD model based designed variables as ; its design sensitivity with respect the k design variable can be formulated as follows:

To obtain a reasonable result, the CAD model should be meshed with more control points by spline knot insertion or order elevation, and they do not change the geometry [1]. If the design variation does not change the weight of the meshes, taking dimension of the design surface as an example, the following formulations can be obtained:

Here, denotes the meshed design surface. indicates the linear relationship between the initial control point coordinates and meshed control points coordinates. denotes the basis function after mesh creation. As a result, (30) can be transformed as

Here, is the number of degrees of freedom of the meshed model. Equation (33) indicates that the design sensitivity can be finally expressed by a linear combination of mesh level sensitivities. Figure 2 shows how the change of the radius of the cylinder is achieved by the variation of mesh level control points. This is the advantage of using IGA, there is no such natural link in the FEM based optimization. This characteristic is used to obtain the design sensitivity from the mesh level sensitivities. The mapping matrix can be calculated analytically according to the NURBS operation [27] or with the finite difference method (FDM). To numerically evaluate by FDM, a selected step should be added to the design variables to create the perturbed model, then both the initial model and the perturbed model are meshed, their corresponding coordinates are rearranged into vectors and , respectively, then subtract from , and divide the result by ; the vector can be obtained. indicates the j component of vector . With this mapping, the final design sensitivities can be obtained by computing the mesh level sensitivities firstly and then mapping them onto the design variables.

#### 4. Numerical Results and Discussions

##### 4.1. Single Curved Free Form Surface Design

Here, a free from surface design example is given to test the sensitivity mapping strategy. Figure 3 shows the initial structure; it is extruded by the boundary curve. There are totally 15 control points. Its boundary condition and material properties are also shown on the picture. The structure is modelled with IGA Kirchhoff-Love shell elements. The coordinates of the profile, denoted on Figure 3, are optimized to obtain the maximum fundamental natural frequency under mass constraint.

The optimization problem can be formulated as follows:

Here, is the design vector. is the first eigenvalue; it relates to the fundamental natural frequency with . is the corresponding eigenvector. denotes the structure mass, . is the user defined maximum mass; here we set it as .

The analysis model is shown in Figure 4, there are 450 quadratic elements and 544 control points, and the total degrees of freedom are 1632. At the element level, all the movable control point positions are considered as designs variables. Let denote a coordinate component of a control point of the analysis model. The element level sensitivities of the objective function and constraint are

Here, and are obtained with (27) and (28). Figure 5 shows the sensitivity of the natural frequency with respect to the coordinates of the movable control points of the fine meshed model, it is denoted as , and there are 480 mesh level design variables.

For the initial model shown in Figure 3, there are 15 movable control points, and their design sensitivity can be obtained by the sensitivity mapping as follows:

As has been stated in (33), is can be obtained by NURBS operation algorithm [27]. indicates that a perturbation of the th design variable will lead to movement of the th fine mesh control point coordinate component. Here, denotes the sensitivities if we take the coordinates of the 15 movable control points of the initial model as design variables. In the problem set-up, the surface is obtained by extrusion; thus, the 15 movable control points are categorized into 5 groups, and each group has its own coordinate that corresponds to the initial 5 design variables defined on the boundary curve. The final profile design sensitivities can be obtained by a secondary mapping from to . Specifically speaking, is obtained by summing up all its corresponding group member sensitivities , since the model is extruded by the profile. The overall design sensitivity is thus obtained analytically.

Figure 6 shows the design sensitivity with respect to the initial model control points (15 design variables), obtained with (37). Figure 7 shows the final design sensitivity with respect to the user defined variables (5 design variables).

According to (33), the mapping from mesh level sensitivity to the final design sensitivity can also be achieved by one step with FDM method; that is,

Here, is obtained by directly calculating (33) with a small . In this example, is used. The procedure is semianalytical, which means that the mapping itself is obtained by FDM, i.e., geometrical perturbation, and the mesh level sensitivity is calculated analytically.

In order to make comparisons, the design sensitivity (five design variables) is also calculated by FDM method (step ), the design variables are firstly perturbed with the previous value in turns, then the mesh is created, and the perturbed structure responses are evaluated and compared to the initial responses to get the design sensitivity. The sensitivities obtained from different methods are shown in Table 1. It can be seen that the two analytical methods obtain the same results, and they agree with the results from the full FDM method. It shows the correctness of the sensitivity mapping procedure.

The optimization is solved with function in MATLAB platform. The “active set" algorithm method is adopted. The sensitivities are evaluated with full analytical approach. Figure 8 shows the iteration history. The optimal solution is achieved with merely 10 iterations. It is very fast. Table 2 gives the optimum geometrical coefficients. Figure 9 shows the optimum structure. The structure gets the maximum allowable mass; it is quite reasonable since the natural frequency decreases as the mass increases. It is also noteworthy that the structure is symmetric, and no such constraint is prescribed in prior; it indicates that the sensitivity information obtained in each iteration by our approach is quite accurate.

##### 4.2. Built-Up Shell Structure Design

A built-up shell structure design example is given here. It is composed of two parts depicted in Figure 10. They are all modelled with Kirchhoff-Love IGA shell elements. The fundamental natural frequency of the full structure will be maximized under the mass and strain energy constraints. In this design, only the shape of the curved part is changeable, and the flat plate shape is kept unchanged during the optimization. The prementioned modal synthesis method is adopted to reduce the computational burden. Figure 10 shows the initial shape of the structure. Figure 11 shows the fine meshed analysis model the boundary condition and the martial properties. The design variables are the coordinates of some control points shown in Figure 12.

The optimization problem is formulated as follows:

In addition to the constraints used in the previous example, here we also introduce static analysis equation and strain energy constraints. The substructure is simplified with the prementioned modal synthesis method with only the first 10 fixed interface modes. In this case, works as a macroelement, its initial 621 degrees of freedom are condensed to 136 degrees of freedom, 10 of which are the fixed interface modes, and the rest are the interface degrees of freedom, which is decided by the number of control points on the interfaces. The gain increases as the meshes become finer, since after modal reduction only the interface degrees of freedom and the fixed interface modes are kept. On the other side, the flat part of the structure is unchanging during the iteration; it means the model reduction process will be conducted only one time and its reduced model will be used in all iterations; thus, it can greatly reduce the computational cost.

The strain energy sensitivity is

In the present work, and are set as and , respectively. The function with “active set” algorithm in MATLAB platform is also used here. All the sensitivities are obtained with our full analytical sensitivity mapping approach.

Figures 13, 14, and 15 show the iteration history of the cost functions. Figure 16 shows the optimal shape of the structure. Table 3 shows the coefficients of the initial shape and optimal shape and their corresponding cost function values. A symmetric structure is obtained again, it is noteworthy that such feature is not prescribed before, and it proves the effectiveness of our approach.

##### 4.3. Discussion

The above optimization framework is flexible in treating various types of design variables. The bottom structure ingredient sensitivities are calculated simultaneously with the structure ingredients, and they share a lot of intermediate variables, so it will save some computational resources. After this step, the structure ingredients sensitivity with respect to all degrees of freedom is obtained. Based on them, the sensitivity mapping transforms the bottom sensitivity to the user defined variables. The users only need to define the design variables by defining sensitivity mapping matrix. This framework is promising for the commercial software.

Another aspect that should be stated is the efficiency. In the present example, the time cost in a single step with sensitivity evaluation embedded is about 23 times as a pure cost functions evaluation without sensitivity evaluation. It indicates that only when the number of design variables is larger than 22, the optimization based on analytical sensitivity is more economical in terms of the time costs in a single iteration. However, the exactness of the analytical sensitivity leads to less iteration steps, and this will save some time. The incorporation of modal synthesis method in the built-up structure shape optimization also makes it more efficient, since the unchanging parts are only calculated one time and their degrees of freedom are greatly reduced with modal synthesis method.

#### 5. Conclusion

In this paper, a framework for shell shape optimization considering structural dynamic criteria is proposed. It is based on the IGA Kirchhoff-Love shell element. The mesh level sensitivities are calculated firstly; then sensitivity mapping is used to extract the user defined design sensitivity. It makes the program general and flexible. Modal synthesis method is also adopted to avoid a repeated evaluation of the unchanging parts of the structure. Only the parts where the design variables are defined are reanalyzed, the full structure responses are based on the synthesis of the substructures, and it is more time efficient. Examples show that our method is promising in dealing with free form shape optimization problem of large number of design variables considering dynamics criteria.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: this work was supported by the French ANR Project PEPITO (ANR-14-CE23-0011), the Fundamental Research Funds for the Central Universities (Grants nos. 300102258108 and 300102258201), the National Natural Science Foundation of China (Grant no. 51509006), and the Natural Science Basic Research Plan in Shaanxi Province of China (Grant no. 2018JQ5038).