# Shape Optimization for Natural Frequency with Isogeometric Kirchhoff-Love Shell and Sensitivity Mapping

Authors: Zhen Lei; Frederic Gillot; Louis Jezequel
Journal: Mathematical Problems in Engineering
 (2018)

---

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


---

## Body

## 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 Section2. 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(1)X=Sξ1,ξ2+ξ3N3.The deformed shell body is(2)x=sξ1,ξ2+ξ3n3.Here,S and N3 denote the shell middle surface and the shell normal before deformation and s and n3 for the deformed one. The notation ·,ξi denotes ∂·/∂ξi and is simply noted as ·,i hereinafter. The scripts noted by the Latin letters i, j belong to 1,2,3, and the Greek ones α,β∈1,2.(3)N3=S,1×S,2S,1×S,2,n3=s,1×s,2s,1×s,2The displacement is(4)u=x-X.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:(5)Eαβ=εαβ+ξ3καβ.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(6)εαβ=12S,α·u,β+S,β·u,α.The bending part is(7)καβ=-u,αβ·N3+1S,1×S,2u,1·S,αβ×S,2+u,2·S,1×S,αβ+N3·S,αβS,1×S,2u,1·S,2×N3+u,2·N3×S,1.According to the principle of virtual work, the following equation can be obtained:(8)δΠint=∫Ωnαβδεαβ+mαβδκαβdΩ=∫Ωb·δudΩ+∫ΓNt·δudΓ=δΠext.Here,nαβ=hCαβγδεγδ, mαβ=(h3/12)Cαβγδκγδ. They are effective stress resultants. Ω indicates the shell middle surface, b is the shell middle surface load density, and t is the traction on the boundary ΓN. Cαβγδ is the material tensor component [26] detailed as(9)Cαβγδ=Eν1-ν2S,α·S,βS,γ·S,δ+E21+νS,α·S,γS,δ·S,β+E21+νS,α·S,δS,γ·S,β.E and ν are the material Young’s modulus and the Poisson ratio, respectively. In IGA analysis, the shell middle surfaces (S and s) are expressed by NURBS [1, 27], so is the displacement field u. The structure control equation for static analysis is(10)Ku=F,(11)Klm=∂δΠint∂δul∂um,Fl=∂δΠext∂δul.The structure control equation for modal analysis is(12)Mu¨+Ku=0.M denotes the mass matrix, it is obtained by condensing mass to the shell middle surface and treating them as inertial force.(13)Mij=∫sρNiNjds


### 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 mapT constructed according to different modal synthesis approaches, the system control equation (12) can be simplified as(14)TT-w2M+KTu=0.The transformation matrixT is a product of Ts and Tc, 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 Ts and Tc. We adopt the most often used fixed interface Craig-Bampton method [21] in the present work.The discrete control equation for a substructureΩs with harmonic loads is(15)-w2MiisMibsMbisMbbs+KiisKibsKbisKbbsuisubs=0.Here, the footnotesi and b 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 t is(16)Tst=ΦNt,Φct.Here,ΦNt denotes the selected fixed interface mode, ΦNt=ϕ1t,ϕ2t,⋯,ϕkt. ϕrt indicates the rth eigenvector of substructure t with all interface degrees of freedom constrained. Only the first k eigenvectors are kept. Φct 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.(17)Φct=-Kiis-1KibsIThe condensed substructure equation now can be obtained as(18)TstT-w2M+KTstu=0The 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 areC0 connected, we can assemble the macroelements in the same way as the FEM method. If they are G1 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].


## 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(1)X=Sξ1,ξ2+ξ3N3.The deformed shell body is(2)x=sξ1,ξ2+ξ3n3.Here,S and N3 denote the shell middle surface and the shell normal before deformation and s and n3 for the deformed one. The notation ·,ξi denotes ∂·/∂ξi and is simply noted as ·,i hereinafter. The scripts noted by the Latin letters i, j belong to 1,2,3, and the Greek ones α,β∈1,2.(3)N3=S,1×S,2S,1×S,2,n3=s,1×s,2s,1×s,2The displacement is(4)u=x-X.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:(5)Eαβ=εαβ+ξ3καβ.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(6)εαβ=12S,α·u,β+S,β·u,α.The bending part is(7)καβ=-u,αβ·N3+1S,1×S,2u,1·S,αβ×S,2+u,2·S,1×S,αβ+N3·S,αβS,1×S,2u,1·S,2×N3+u,2·N3×S,1.According to the principle of virtual work, the following equation can be obtained:(8)δΠint=∫Ωnαβδεαβ+mαβδκαβdΩ=∫Ωb·δudΩ+∫ΓNt·δudΓ=δΠext.Here,nαβ=hCαβγδεγδ, mαβ=(h3/12)Cαβγδκγδ. They are effective stress resultants. Ω indicates the shell middle surface, b is the shell middle surface load density, and t is the traction on the boundary ΓN. Cαβγδ is the material tensor component [26] detailed as(9)Cαβγδ=Eν1-ν2S,α·S,βS,γ·S,δ+E21+νS,α·S,γS,δ·S,β+E21+νS,α·S,δS,γ·S,β.E and ν are the material Young’s modulus and the Poisson ratio, respectively. In IGA analysis, the shell middle surfaces (S and s) are expressed by NURBS [1, 27], so is the displacement field u. The structure control equation for static analysis is(10)Ku=F,(11)Klm=∂δΠint∂δul∂um,Fl=∂δΠext∂δul.The structure control equation for modal analysis is(12)Mu¨+Ku=0.M denotes the mass matrix, it is obtained by condensing mass to the shell middle surface and treating them as inertial force.(13)Mij=∫sρNiNjds


## 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 mapT constructed according to different modal synthesis approaches, the system control equation (12) can be simplified as(14)TT-w2M+KTu=0.The transformation matrixT is a product of Ts and Tc, 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 Ts and Tc. We adopt the most often used fixed interface Craig-Bampton method [21] in the present work.The discrete control equation for a substructureΩs with harmonic loads is(15)-w2MiisMibsMbisMbbs+KiisKibsKbisKbbsuisubs=0.Here, the footnotesi and b 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 t is(16)Tst=ΦNt,Φct.Here,ΦNt denotes the selected fixed interface mode, ΦNt=ϕ1t,ϕ2t,⋯,ϕkt. ϕrt indicates the rth eigenvector of substructure t with all interface degrees of freedom constrained. Only the first k eigenvectors are kept. Φct 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.(17)Φct=-Kiis-1KibsIThe condensed substructure equation now can be obtained as(18)TstT-w2M+KTstu=0The 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 areC0 connected, we can assemble the macroelements in the same way as the FEM method. If they are G1 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(19)MaximizeΩ⊂ΩDJΩs.t.hΩ=0gΩ≤0The first equality constrainth is the structural control equation which is usually a partial differential equation (PDE) and g 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(20)MaximizepJΩps.t.HΩp=0GΩp≤0(21)Ωp=Sξ,η=∑i=1n∑j=1mPi,jRi,jp,qξ,η.Here,Pi,j indicates the spatial coordinate of the i,jth control point of the shape to be optimized, Ri,jp,qξ,η denotes the corresponding NURBS basis function, it is nonnegative compactly supported. As a result, all items of the stiffness and mass matrices are positive. p and q denote the orders of the basis function in two parametric directions, respectively.(22)Pi,j=Pi,jx,Pi,jy,Pi,jz(23)p=Pi,j∣i∈1~m,j∈1~nThe discrete structural control equation is(24)H:MpU¨p+KpUp-Fp=0.


### 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. Figure1 shows an example where the kth control point of the parametric surface is taken as design variable; it indicates an optimization of the z coordinate component of that control point. The design sensitivity can be obtained as the perturbation of design variable tends to zero.(25)pkzδ=pkz+δ(26)JSpkz,pkz=limδ→0⁡Jpkzδ-Jpkzδ⁡Figure 1
Variation of a control point of a NURBS surface.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.(27)Mij,pkz=∂∫SpkzρNiNjdSpkz∂pkz(28)Klm,pkz=∂δΠintpkz∂δul∂um∂pkzDepending 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 asz; its design sensitivity with respect the kth design variable can be formulated as follows:(29)zkδ=zk+δ(30)JSzk,zk=limδ→0⁡JSzkδ-JSzkδ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 z dimension of the design surface as an example, the following formulations can be obtained:(31)Δ=0,0,δT(32)Szkδ-Szk=Sfzkδ-Sfzk=∑j=1mΔTk,jRjfξ,ηHere,Sf denotes the meshed design surface. T indicates the linear relationship between the initial control point coordinates and meshed control points coordinates. Rjf denotes the basis function after mesh creation. As a result, (30) can be transformed as(33)JSzk,zk=⁡limδ→0⁡JSfzk+∑j=1mΔTk,jRjf-JSfzkδ=∑j=1mTk,jJSf,pj.Here,m 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 Tk,j by FDM, a selected step δ should be added to the design variables zk to create the perturbed model, then both the initial model and the perturbed model are meshed, their corresponding coordinates are rearranged into vectors A and B, respectively, then subtract A from B, and divide the result by δ; the vector Tk can be obtained. Tk,j indicates the jth component of vector Tk. 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.Figure 2
The feature variation and the movement of the relevant control points.


## 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. Figure1 shows an example where the kth control point of the parametric surface is taken as design variable; it indicates an optimization of the z coordinate component of that control point. The design sensitivity can be obtained as the perturbation of design variable tends to zero.(25)pkzδ=pkz+δ(26)JSpkz,pkz=limδ→0⁡Jpkzδ-Jpkzδ⁡Figure 1
Variation of a control point of a NURBS surface.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.(27)Mij,pkz=∂∫SpkzρNiNjdSpkz∂pkz(28)Klm,pkz=∂δΠintpkz∂δul∂um∂pkzDepending 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 asz; its design sensitivity with respect the kth design variable can be formulated as follows:(29)zkδ=zk+δ(30)JSzk,zk=limδ→0⁡JSzkδ-JSzkδ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 z dimension of the design surface as an example, the following formulations can be obtained:(31)Δ=0,0,δT(32)Szkδ-Szk=Sfzkδ-Sfzk=∑j=1mΔTk,jRjfξ,ηHere,Sf denotes the meshed design surface. T indicates the linear relationship between the initial control point coordinates and meshed control points coordinates. Rjf denotes the basis function after mesh creation. As a result, (30) can be transformed as(33)JSzk,zk=⁡limδ→0⁡JSfzk+∑j=1mΔTk,jRjf-JSfzkδ=∑j=1mTk,jJSf,pj.Here,m 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 Tk,j by FDM, a selected step δ should be added to the design variables zk to create the perturbed model, then both the initial model and the perturbed model are meshed, their corresponding coordinates are rearranged into vectors A and B, respectively, then subtract A from B, and divide the result by δ; the vector Tk can be obtained. Tk,j indicates the jth component of vector Tk. 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.Figure 2
The feature variation and the movement of the relevant control points.


## 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. Figure3 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 z coordinates of the profile, denoted on Figure 3, are optimized to obtain the maximum fundamental natural frequency under mass constraint.Figure 3
The problem setup.The optimization problem can be formulated as follows:(34)findz=z1,z2,…,zimaximizeλ1zs.t.mz≤mdKzu1=λ1Mzu1u1TMu1=1u1TKu1=λ1Here,z is the design vector. λ1 is the first eigenvalue; it relates to the fundamental natural frequency with λ1=2πf12. u1 is the corresponding eigenvector. m denotes the structure mass, m=(1/3)sumM. md is the user defined maximum mass; here we set it as 3000.The analysis model is shown in Figure4, 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 pk denote a coordinate component of a control point of the analysis model. The element level sensitivities of the objective function and constraint are(35)λ1,pk=u1TM,pk-λ1K,pku1(36)m,pk=13sumM,pk.Figure 4
The fine meshed model.Here,K,pk and M,pk are obtained with (27) and (28). Figure 5 shows the sensitivity of the natural frequency with respect to the z coordinates of the movable control points of the fine meshed model, it is denoted as λ1,pk, and there are 480 mesh level design variables.Figure 5
The natural frequency sensitivity with respect toz coordinates of the movable control points of the analysis model.For the initial model shown in Figure3, there are 15 movable control points, and their design sensitivity can be obtained by the sensitivity mapping as follows:(37)λ1,pck=Taλ1,pk.As has been stated in (33), Ta is can be obtained by NURBS operation algorithm [27]. Taij indicates that a perturbation of the ith design variable will lead to Taij movement of the jth fine mesh control point coordinate component. Here, λ1,pck 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 z 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 λ1,pc to λ1,z. Specifically speaking, λ1,zk is obtained by summing up all its corresponding group member sensitivities λ1,pcl, since the model is extruded by the profile. The overall design sensitivity is thus obtained analytically.Figure6 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).Figure 6
Design sensitivity with respect to the initial model control points.Figure 7
Design variables sensitivity.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,(38)λ1,z=Tfλ1,pk.Here,Tf is obtained by directly calculating (33) with a small δ. In this example, δ=10-6 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 (step10-5), 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.Table 1
Design sensitivities obtained from the fullFDM method, full analytical method (Anal1), and FDM-analytical method (Anal2).
z 1 z 2 z 3 z 4 z 5 FDM 16.51388 5.690376 -4.635494 5.690507 16.51380 Anal1 16.51392 5.690476 -4.635433 5.690476 16.51392 Anal2 16.51392 5.690476 -4.635433 5.690476 16.51392The optimization is solved withfmincon 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.Table 2
Free form surface design results: the initial and optimal values and the structural response.
design variables frequency mass initial 1, 2, 3, 2, 1 3.2752 2135.7 optimum 6.8403, 5.9856, 6.2363, 5.9856, 6.8403 3.8891 3000.0Figure 8
Shape design iteration history.Figure 9
Optimum shape from the free form surface design.


### 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 Figure10. 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 z coordinates of some control points shown in Figure 12.Figure 10
The geometrical model and the substructures.Figure 11
The analysis model and the simulation setup.Figure 12
The design part and the design variables for the built-up shell design.The optimization problem is formulated as follows:(39)findz=z1,z2,…,z25maximizeλ1zs.t.mz≤mdW=uTKzu≤WdKzu=FKzu1=λ1Mzu1u1TMu1=1u1TKu1=λ1In addition to the constraints used in the previous example, here we also introduce static analysis equation and strain energy constraints. The substructureΩ2 is simplified with the prementioned modal synthesis method with only the first 10 fixed interface modes. In this case, Ω2 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(40)W,=-uTK,u.In the present work,Wd and md are set as =8×105 and 2700, respectively. The fmincon function with “active set” algorithm in MATLAB platform is also used here. All the sensitivities are obtained with our full analytical sensitivity mapping approach.Figures13, 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.Table 3
Built-up shell design, the initial and optimum design variable values, and the structural response.
design variables frequency mass strain energy initial -1,-2,-4,-2,-1,-1,-2,-4,-2,-1,-1,-2,-4 7.7210 2477.0 1.0546×106 -2,-1,-1,-2,-4,-2,-1,-1,-2,-4,-2,-1 optimum 3.1609, 7.8490, 10.6005, 7.8490, 3.1609 11.5760 2700.0 8.0×105 2.0118, 5.4873, 5.5142, 5.4873, 2.0118 1.3932, 5.0478, 6.4851, 5.0478, 1.3932 2.0118, 5.4873, 5.5142, 5.4873, 2.0118 3.1609, 7.8490, 10.6005, 7.8490, 3.1609Figure 13
The objective function.Figure 14
The strain energy constraint.Figure 15
The mass constraint.Figure 16
The optimum shape.


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


## 4.1. Single Curved Free Form Surface Design

Here, a free from surface design example is given to test the sensitivity mapping strategy. Figure3 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 z coordinates of the profile, denoted on Figure 3, are optimized to obtain the maximum fundamental natural frequency under mass constraint.Figure 3
The problem setup.The optimization problem can be formulated as follows:(34)findz=z1,z2,…,zimaximizeλ1zs.t.mz≤mdKzu1=λ1Mzu1u1TMu1=1u1TKu1=λ1Here,z is the design vector. λ1 is the first eigenvalue; it relates to the fundamental natural frequency with λ1=2πf12. u1 is the corresponding eigenvector. m denotes the structure mass, m=(1/3)sumM. md is the user defined maximum mass; here we set it as 3000.The analysis model is shown in Figure4, 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 pk denote a coordinate component of a control point of the analysis model. The element level sensitivities of the objective function and constraint are(35)λ1,pk=u1TM,pk-λ1K,pku1(36)m,pk=13sumM,pk.Figure 4
The fine meshed model.Here,K,pk and M,pk are obtained with (27) and (28). Figure 5 shows the sensitivity of the natural frequency with respect to the z coordinates of the movable control points of the fine meshed model, it is denoted as λ1,pk, and there are 480 mesh level design variables.Figure 5
The natural frequency sensitivity with respect toz coordinates of the movable control points of the analysis model.For the initial model shown in Figure3, there are 15 movable control points, and their design sensitivity can be obtained by the sensitivity mapping as follows:(37)λ1,pck=Taλ1,pk.As has been stated in (33), Ta is can be obtained by NURBS operation algorithm [27]. Taij indicates that a perturbation of the ith design variable will lead to Taij movement of the jth fine mesh control point coordinate component. Here, λ1,pck 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 z 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 λ1,pc to λ1,z. Specifically speaking, λ1,zk is obtained by summing up all its corresponding group member sensitivities λ1,pcl, since the model is extruded by the profile. The overall design sensitivity is thus obtained analytically.Figure6 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).Figure 6
Design sensitivity with respect to the initial model control points.Figure 7
Design variables sensitivity.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,(38)λ1,z=Tfλ1,pk.Here,Tf is obtained by directly calculating (33) with a small δ. In this example, δ=10-6 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 (step10-5), 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.Table 1
Design sensitivities obtained from the fullFDM method, full analytical method (Anal1), and FDM-analytical method (Anal2).
z 1 z 2 z 3 z 4 z 5 FDM 16.51388 5.690376 -4.635494 5.690507 16.51380 Anal1 16.51392 5.690476 -4.635433 5.690476 16.51392 Anal2 16.51392 5.690476 -4.635433 5.690476 16.51392The optimization is solved withfmincon 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.Table 2
Free form surface design results: the initial and optimal values and the structural response.
design variables frequency mass initial 1, 2, 3, 2, 1 3.2752 2135.7 optimum 6.8403, 5.9856, 6.2363, 5.9856, 6.8403 3.8891 3000.0Figure 8
Shape design iteration history.Figure 9
Optimum shape from the free form surface design.


## 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 Figure10. 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 z coordinates of some control points shown in Figure 12.Figure 10
The geometrical model and the substructures.Figure 11
The analysis model and the simulation setup.Figure 12
The design part and the design variables for the built-up shell design.The optimization problem is formulated as follows:(39)findz=z1,z2,…,z25maximizeλ1zs.t.mz≤mdW=uTKzu≤WdKzu=FKzu1=λ1Mzu1u1TMu1=1u1TKu1=λ1In addition to the constraints used in the previous example, here we also introduce static analysis equation and strain energy constraints. The substructureΩ2 is simplified with the prementioned modal synthesis method with only the first 10 fixed interface modes. In this case, Ω2 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(40)W,=-uTK,u.In the present work,Wd and md are set as =8×105 and 2700, respectively. The fmincon function with “active set” algorithm in MATLAB platform is also used here. All the sensitivities are obtained with our full analytical sensitivity mapping approach.Figures13, 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.Table 3
Built-up shell design, the initial and optimum design variable values, and the structural response.
design variables frequency mass strain energy initial -1,-2,-4,-2,-1,-1,-2,-4,-2,-1,-1,-2,-4 7.7210 2477.0 1.0546×106 -2,-1,-1,-2,-4,-2,-1,-1,-2,-4,-2,-1 optimum 3.1609, 7.8490, 10.6005, 7.8490, 3.1609 11.5760 2700.0 8.0×105 2.0118, 5.4873, 5.5142, 5.4873, 2.0118 1.3932, 5.0478, 6.4851, 5.0478, 1.3932 2.0118, 5.4873, 5.5142, 5.4873, 2.0118 3.1609, 7.8490, 10.6005, 7.8490, 3.1609Figure 13
The objective function.Figure 14
The strain energy constraint.Figure 15
The mass constraint.Figure 16
The optimum shape.


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



---

*Source: 9531651-2018-06-26.xml*