Asymmetric Design Sensitivity and Isogeometric Shape Optimization Subject to Deformation-Dependent Loads

: We present a design sensitivity analysis and isogeometric shape optimization with path-dependent loads belonging to non-conservative loads under the assumption of elastic bodies. Path-dependent loads are sometimes expressed as the follower forces, and these loads have characteristics that depend not only on the design area of the structure but also on the deformation. When such a deformation-dependent load is considered, an asymmetric load stiffness matrix (tangential operator) in the response region appears. In this paper, the load stiffness matrix is derived by linearizing the non-linear non-conservative load, and the geometrical non-linear structure is optimally designed in the total Lagrangian formulation using the isogeometric framework. In particular, since the deformation-dependent load changes according to the change and displacement of the design area, the isogeometric analysis has a signiﬁcant inﬂuence on the accuracy of the sensitivity analysis and optimization results. Through several numerical examples, the applicability and superiority of the isogeometric analysis method were veriﬁed in optimizing the shape of the problem subject to deformation-dependent loads.


Introduction
Special attention is required in the structural analysis of the pressure loading, as its loading direction could change with the deformation of the structure. Engineering problems with deformation-dependent loading are frequently encountered, for example, in the design of dams, tires, airbags, and pressurized vessels. When such a deformationdependent load is considered, an asymmetric tangential operator expressed as a load stiffness matrix appears in the governing equation. Hibbit [1] first mentioned the concept of load stiffness, and in the case of a surface traction load, the load stiffness is symmetric if it is applied to the fully enclosed volume. If it is not the case, the equation has an asymmetric term which is to be solved using the approximation method, assuming it is symmetric while refining the mesh.
In the previous study from Mok et al. [2], the pressure load was decomposed to the body-attached pressure load and the space-attached pressure load. They derived the load stiffness term by linearizing the nonlinear pressure load term. Simo [3] expressed the pressure-loaded moving boundary in an isoparametric presentation and obtained the tangent stiffness through linearization of the loading boundary. In addition, the convergence rate dependent on tangent stiffness was discussed in numerical examples for an axisymmetric case.
If the pressure loading is considered in the design optimization, it becomes a function of both the deformation by loading and the design change on boundaries. Many previous The fixed load, p 0 , in Figure 1a is a distributed load defined at initial configuration without any changes, and the design-dependent load, p(d), in Figure 1b is the pressure load with the changed loading direction according to perturbed design (d). The symmetric shape sensitivity analysis for the design-dependent load in Figure 1b with the assumption of infinitesimal strain was previously derived and utilized in several optimization problems in the isogeometric framework [12,13].
On the other hand, the design-and deformation-dependent load, p(d, u), in Figure 1c changes according to the deformation (u) and the variation of design (d), which leads to an asymmetric load stiffness matrix by linearization in the response and design sensitivity analysis. In this paper, the isogeometric continuum-based sensitivity analysis suggested by Cho and Ha [14] is performed for accurate sensitivity analysis of the problems subject to a design-and deformation-dependent load. The design boundary change expressed by the NURBS basis function leads to a natural shape change, which leads to accurate design shape sensitivity. The design sensitivity analysis based on the geometric analysis method was induced by the design sensitivity theory [15,16], and its accuracy is verified for the response field. The asymmetric sensitivity equation for the deformation-dependent loadings is expressed as the dependent term according to the displacement and the design change of the structure in Figure 1c, including the geometric higher-order terms such as normal vector and curvature. Compared to the traditional finite element-based sensitivity analysis, NURBS-based design sensitivity analysis results in accurate sensitivity values because the load-stiffness matrix by the deformation-and design-dependent loads, in general taking an asymmetric form, is incorporated into the sensitivity equation.
The remainder of this paper is organized as follows. Section 2 of this paper briefly explains the B-spline basis function and provides a brief explanation of the features of the geometric analysis based on it. We formulate a load representation that describes the characteristics of a pressure load in a continuum form and develop a governing equation and a weak formulation for a geometric nonlinear structure in a Lagrangian formulation, which incorporates the asymmetric load stiffness terms. The details on our formulation of the design sensitivity equations for geometric nonlinear structures subjected to deformationdependent pressure loads are also described. In Section 3, numerical examples showing the accuracy of isogeometric continuum-based sensitivity compared based on the isogeometric approach to that of finite difference for the structures subjected to deformation-dependent loads are presented. Through the shape optimization examples, the efficiency and applicability of the proposed method will also be checked. Finally, in Section 4, we summarize the conclusions, explaining the importance of the isogeometric analysis method in nonlinear structures subjected to pressure loads.

Overview of NURBS Basis Function and Isogeometric Analysis
In isogeometric analysis (IGA), the solution space adopts the same NURBS basis function to describe the geometry. Over the conventional finite element method, IGA has several advantages: exact geometry representation and simple refinement process non-interacting with CAD system due to the use of NURBS basis functions based on B-splines [9]. With the set of knots ξ i in parametric space, a knot vector Ξ is defined in a one-dimensional space, where p is the order of the basis function, and n is the number of control points. The B-spline basis functions are defined, recursively, as and Using the B-spline basis function N i,p (ξ) and weight w i , the NURBS basis function R i,p (ξ) is defined as In general, the higher-order basis functions in IGA offer higher regularity than the conventional finite element approach. The NURBS possesses the following desirable properties as a basis function.
For the given n pairs of the p-th order NURBS basis function R i,p (ξ) and the corresponding projected weighted control point B i ≡ B(x i ), a NURBS curve C is obtained by In the same manner, NURBS surface is defined as a product of NURBS basis functions, In isogeometric analysis, geometric points and responses are expressed using the relation (6): and z(Ξ) = ∑ W p where B I (x) = B(x 1 , x 2 , x 3 ) are the control points, and u I = u(x 1 , x 2 , x 3 ) are called the response coefficients at the control points. Note that the NURBS basis functions have a non-interpolatory property. The equilibrium of a deformable body at current configuration (n + 1) under deformationand design-dependent load can be expressed by using the principle of virtual work in total Lagrangian formulation as a n+1 z, z = n+1 z, z (9) where the strain energy and load form are defined, respectively, as

Isogeometric
and and c ijkl , z, z and Z are material response tensor, displacement, virtual displacement, and variational space, respectively. Generally, the pressure boundary at the current configuration (n + 1) can be represented as shown in Figure 2: where F and J are deformation gradient and Jacobian, respectively. The Green-Lagrange and virtual strain tensors are defined by and Since the strain energy form and load form are nonlinear in their arguments, Equation (9) cannot be solved directly. To solve a nonlinear system of equations, an incrementaliterative scheme is adopted. The solution is decomposed into two terms, the solution at configuration (n) and the increment as Using Equation (15), the Green-Lagrange tensor is decomposed into where Additionally, load form can be linearized by taking Taylor expansion with respect to the structural response and dropping the higher-order terms n+1 z, z = n z, z + sti f f n z; ∆z, z , where the load stiffness term is defined as Substituting Equations (16), (18), and (20) into Equation (9) and linearizing by neglecting the higher-order terms, Equation (9) can be rewritten, in incremental form, as a * * n z; ∆z, z ≡ a * n z; ∆z, z − sti f f n z; ∆z, z = * * n z; z for all z ∈ Z. (22) The linearized strain energy form and load linear form at the configuration (n) are defined as where the second Piola-Kirchhoff stress tensor S ij is defined as Since the load stiffness in Equation (21) has asymmetric property, the left-hand side of the final system equation, Equation (22), is also non-symmetric through linearization. The response solution is obtained by repeating iterative steps in which ∆z is computed in Equation (25) until updated n+1 z is converged within a certain tolerance.
Using the relation (8), the incremental variational Equation (22) can be rewritten as where the linearized energy form is discretized as and the load form as The p, CP, and W denote the order of basis function, the number of control points, and the NURBS basis function for the boundary integral, respectively. Note that the domain and boundary in Equations (27) and (28) are maintained exactly with the discretization, different from the conventional finite element discretization.

Isogeometric Shape Sensitivity Analysis
Taking the material derivative of Equation (9) with respect to the shape design parameter τ, the following is obtained: and where a ex ( n+1 z, z) and ex ( n+1 z; z) denote the explicit variation terms in which the dependence of their arguments on the shape design parameter is suppressed. Using the relation n+1 . z = n+1 z + ∇ n+1 z T V (design velocity V), Equations (29) and (30) are rewritten as and Additionally, using the fact that a n+1 z, z = n+1 z, z for all z ∈ Z and a n+1 z, z for all . z ∈ Z, the shape sensitivity equation is finally derived as where and e ij n+1 z; z , V = 1 2 Using the isogeometric discretization, the variational Equation (33) can be rewritten as where a * * ( n+1 u; and Note that the normal vector in the boundary integrals is crucial for an accurate evaluation of shape sensitivity. Due to the exact consideration of higher-order terms, the isogeometric DSA gives accurate values over that of the finite element, which eventually could lead to the precise results of shape design optimization. NURBS basis function requires a recursive evaluation as opposed to a fixed polynomial function in traditional FEM, which is a computational burden if a large number of degrees of freedom (DOF) is used. The asymmetric load stiffness term in Equations (34) and (35) is identified in the continuum-based shape sensitivity formulation to problems subject to deformationdependent loads.

Results and Discussion-Numerical Examples
In the shape optimization using the traditional finite element method, design parametrization and regularity issues occur. In the present study, the shape design method using the B-spline curve (Bribant and Fluery [17]) and the smoothing gradient method (Azegami et al. [18]) were utilized to avoid these problems. Note that the free shape and continuity adjustment of the geometries using the NURBS basis function has an advantage in shape optimization. Due to the flexible nature of the NURBS basis function, engineering meaningful solutions can be easily obtained without any modification to the sensitivity.
In this section, the derived shape sensitivity using the isogeometric approach is compared with FEM-based sensitivity and also with the exact analytic solutions for verification. The obtained isogeometric sensitivity means a gradient of the objective function, which is utilized as important information in the gradient-based optimization algorithm, the modified method of feasible direction (MMFD) in this study. MMFD is a gradient-based optimization solver for constrained problems. The DOT package by Vanderplaats Research and Development [19] is used to apply the MMFD algorithm for the optimization problems. Additionally, the quadratic NURBS basis function is used in all numerical models for both IGA and FEM. Figure 3 shows the model description and its deformed shape for the structure subjected to pressure load. Pressure intensity (p) is 5.0 × 10 9 N/m 2 , and the steel material is used with a Young's modulus of E = 2.1 × 10 11 Pa and a Poisson's ratio of ν = 0.3. The termination criteria of Newton-Raphson analysis for the nonlinear system is defined by ∑ Total DOFs j ∆u j 2 / ∑ Total DOFs j u j 2 with the tolerance of 1.0 × 10 −8 . As shown in Figure 3, the quarter model is used for the investigation with an account of the symmetric boundary condition. To verify the derived shape sensitivity equation, an irregular shape perturbation as shown in Figure 3b is considered. In Table 1, excellent agreement is observed between the finite difference sensitivity and the analytical one at sampled points using 0.01% perturbation. Additionally, the dependency on the magnitude of pressure load and perturbation amount is investigated. The numerical model with a height of 2 m and a width of 1 m is presented as in Figure 4, and pressure loading is applied on the left wall of the simplified dam. The material property and the tolerance of the iterative solution are the same as in the previous numerical example in Figure 3. Considering the geometrical nonlinearity, von Mises stress results and red-colored design velocity vectors are presented in Figure 4b,c. The performance measure is the x-directional displacement of the right upper corner of the structure. Table 2 shows the analytical sensitivity and finite difference sensitivity, and its agreements are presented in Figure 5. It can be noted that the accuracy of shape design sensitivity decreases as the magnitude of pressure increases, and the accuracy increases as the perturbation amount is further reduced under the same loading.

Superior Convergence of Isogeometric Sensitivity (FEM vs. IGA)
A hollow cylinder subjected to inner pressure is considered in Figure 6a. For the comparison with the exact solution, the linear elastic structure problem is considered. On the inner boundary, pressure which has 10 N/m 2 is applied. The displacement and traction are free on the outer boundary. The Young's modulus and Poisson's ratio are 1.0 × 10 5 (MPa) and 0.3, respectively. Inner and outer radii are set as 1 and 7 m. Results of the quarter model for both FEM and IGA are shown in Figure 6b,c using 81 elements. The analytic (exact) solutions to this problem are as follows: The error norm between the numerical results (FEM and IGA) and the analytic solution is considered to assess the convergence characteristics. The L 2 error norm is defined on the internal domain: (43) Figure 7 compares the error norm of FEM and IGA with the number of elements. Because of the exact geometry and higher continuity between elements in IGA, the convergence of the isogeometric method is faster than that of FEM. Next, to investigate the convergence of the isogeometric shape design sensitivity, the previous example is considered again. The inner radius is chosen as a design variable, and the analytic exact design sensitivity is derived as follows: The convergence of the two methods is shown in Figure 8, and the isogeometric approach converges faster than FEM. Such results can be attributed the higher-order geometric term (normal vector), which can be obtained in the isogeometric method even with using fewer DOFs.

Shape Design Optimization of Structures Subject to Pressure Load
Utilizing accurate and efficient design sensitivity, shape optimization is performed for the various problems. Under the constraint, the objective function such as compliance is minimized using the MMFD algorithm. For the comparison, FEM, IGA, and linear elasticity-based optimization results are presented. The material properties of the steel were the same used in the previous problem.

L-Shape Structure
The shape design optimization to minimize compliance at the (n+1) configuration is stated as follows: where d i is a shape design variable, i.e., control points along the design boundary. The Lshaped structure with a height of 20 m and thickness of 5 m was selected for the numerical example as shown in Figure 9a. Pressure is applied on the outer boundary, and the symmetry condition is imposed, which is represented by the circular shape 'o' in Figure 9a. Design variables are on the inner and outer boundary as blue-colored control points in Figure 9b. The material used is steel, which is defined in Section 3.1. The shape sensitivity of Equations (45) and (46) can be derived using n+1 . u from Equation (38) as follows: where V i is a shape design velocity.  Figure 10 shows the stress distribution of the initial and optimized model. In the initial model, stress concentration occurs at the inner corner. On the contrary, the optimal shape has a round one, and therefore uniform stress distribution is observed. The optimization history of compliance is shown in Figure 11. The compliance is decreased by about 32% from the initial compliance of 2.011 × 10 8 to 1.368 × 10 8 after 15 design iterations.

Rectangular Plate
In this example, the shape design optimization of a rectangular structure with a width of 20 m and a height of 10 m considering linear and nonlinear elasticity is performed. The numerical model and design boundary are described in Figure 12. The pressure loading is applied on the lower boundary, and both ends are simply supported, which is represented by the triangle in Figure 12. The material used is steel, which is defined in Section 3.1. The objective function is compliance, and the allowable volume is 60% of the initial one.  Figure 13 shows the optimal shapes for both problems, which have a round shape to resist pressure loads. The optimized shape in Figure 13a has a wiggly lower boundary, which is obtained with consideration of a design-dependent load based on a linear elastic formulation [12,13]. However, the optimized shape considering the design-and deformation-dependent load in the geometrically nonlinear formulation shows a smoother lower boundary where the pressure loading is applied. Optimization histories for both problems are presented in Figure 14. Since the MMFD algorithm meets the constraint first, it is observed that the objective function goes up first in the second step. From the next step, 60% volume is maintained, satisfying the constraint. Although only 60% of the initial mass can be used, it can be seen that it has similar compliance compared to the initial one.

Irregular Shape Plate
In this example, FEM-and IGA-based shape design optimization is performed. A quarter model of a 5 m × 5 m plate with an irregular inner hole and design boundary is described in Figure 15. The pressure loading is applied on the irregular boundary, which is the design boundary, and the roller boundary condition is imposed on the lower and right sides. The material used is steel, which is defined in Section 3.1. The objective function is compliance, and the allowable volume is the initial one. As expected, the optimal boundary shape subjected to pressure loading is round, and Figure 16 shows a round optimal shape in both methods. However, the IGA-based optimal shape shows a smoother boundary represented by NURBS. Optimization histories for both problems are shown in Figure 17. Compliance of the IGA-based optimal shape is slightly lower compared to that of FEM.

Conclusions
In the present study, design sensitivity analysis and shape optimization were performed on nonlinear structures subjected to loads that depend on design change and structural deformation. General non-conservative loads have path-dependent properties, and in this paper, the pressure loads are expressed in a total Lagrangian formulation. The resultant load stiffness matrix is derived by linearization, and a tangential stiffness matrix with an asymmetric structure is also derived in the present study. The asymmetric load stiffness term is identified in the continuum-based shape sensitivity formulation to the problems subject to deformation-dependent loads. In addition, since the total Lagrangian formulation was used, no additional iterations are required in the design sensitivity analysis using the final tangential stiffness matrix converged in the nonlinear analysis. In terms of optimization, since normal vectors, which are high-order geometrical terms, can be accurately obtained based on geometrical geometry, it shows an excellent convergence rate even for the problem of a structure under pressure loading, which is expressed with normal vectors of boundaries. Finally, there was no need for design parameterization and smooth design changes by the design variables, as the NURBS basis function was used in optimizing geometry. Analyses in the present study suggest the superiority of the isogeometric optimization technique, which demonstrates the fast convergence with more accurate calculations of the load stiffness terms as well as the normal vector of the pressure load.