Optimization of shells modeled using NURBS subjected to a dynamic loading condition

This study presents an optimization method for the thickness of shell elements subjected to dynamic loads. Structural analyses were carried out using the finite element method. The analyzed domains were modeled using NURBS, and meshes were generated using a transformation of the parametric domain into the geometric domain via a geometric function. The shell element considered is a combination of a CST membrane element and a DKT plate element, forming an element with 15 degrees of freedom. The Newmark method with constant acceleration was applied to solve the equation of motion. The optimization approach was considered in two different ways: with uniform and variable thickness throughout the shell element. Due to the nonlinear constraints of the problem, the Sequential Quadratic Programming (SQP) method was employed. SQP routines are available in MATLAB, in which this study was performed. Useful examples were detailed in this study to demonstrate the applicability of the optimization method to real structures, such as the Igrejinha da Pampulha, a church located in Belo Horizonte, Brazil, whose modeling was performed.


Optimization of shells modeled using NURBS subjected to a dynamic loading condition
The evolution of computers over the past decades has allowed for the development of new technologies in the field of structural engineering. In turn, this enabled the construction of more slender and economical buildings. One of the main factors that influenced this change is the possibility of combining routines for solving optimization problems with others in order to analyze and design structures. One example is the fact that not long ago, several abacuses were used in the process of designing a structure, whereas today it is very unlikely that a structural design office will not employ computer software for the same purpose.
Despite this evolution, the analysis and design of shells is still not a trivial task, with its study usually limited to specific courses or graduate classes. The objective of this study is to combine a few of the theories that aim to facilitate the modeling of shells (and, consequentially, plates) using thickness analyses and optimization processes. With this approach, it is possible to perform a preanalysis to provide a sketch of the thickness of the designed structure, which nonetheless does not make the design process unnecessary. The Non-Uniform Rational B-Splines (NURBS) representation was used to model the domains by projecting a plane mesh over a surface. The following studies support this methodology: Kiendl et al. (2009) also demonstrate shells modeled using NURBS, employing a Kirchhoff-Love shell element for a nonlinear isogeometric analysis, and present examples that show the applicability of this modeling method to the integration of design and analysis . Breitenberger et al. (2015) present a nonlinear analysis of a computer-aided shell design modeled with NURBS via boundary representation, also using Kirchhoff-Love shells and isogeometric analysis. Lastly, Zhang et al. (2017) show an analysis of complex tubular shells modeled by isogeometric shell analysis using degenerated NURBS elements and the Mindlin-Reissner theory. To verify the efficiency of this method, they constructed a Computer-Aided Design (CAD) model of a car body and simulated its deformation.
For the analysis performed in this study, a combination of a membrane element and a plate element was used to form a shell element. A dynamic analysis method was adopted due to the complexity of shells. In the optimization process, a deterministic method was applied to problems with nonlinear constraints. The studies below support this analysis technique: Falco et al. (2004a) present the development and application of a computational tool for geometry modeling, mesh generation, structural analysis (a Huang-Hinton element, based on a Mindlin-Reissner shell, was used) and sensitivity analysis of plates and shells under dynamic loads. Falco et. al (2004b) used structural sizing and shape optimization procedures to obtain optimum designs for plates and shells under dynamics loads, and the Sequential Quadratic Programming algorithm was used in this process (this article is a continuation of Falco et al. (2004a). Park and Dang (2010) present a structural optimization method based on Computer-Aided-Design and Computer-Aided-Engineering integration, using metamodeling techniques to reduce the time taken for solving the problem and to automatize the process. Alves and Vaz (2013) present the optimization of plates subject to random dynamic loading. Espath et al. (2011) demonstrates the shape optimization of shells (combining NURBSbased surfaces, mathematical optimization, the finite element method in a structural analysis, and automatic differentiation) by applying varied techniques to modify the shell geometry while conserving the same parameterization without generating a new finite element mesh. Kang and Young (2016) present the topology optimization of shells using isogeometric analysis and trimmed NURBS surfaces to handle complex domains that cannot be analyzed with a conventional isogeometric analysis. Alvez and Vaz (2016) carry out a comparative analysis of the thickness optimization of plates under dynamic loads in time domain and its equivalent optimization in frequency domain.
The objective of this article is to present a formulation of the optimization problem of a shell structure under dynamic loads, using NURBS to create a real model of this structure. It also demonstrates the application of this process to real problems such as that of the Igrejinha da Pampulha (a church), located in Belo Horizonte, Brazil, in a hypothetical situation involving dynamic loads. As constraints of the optimization problem were considered to be the maximum displacements suffered by the shell and the minimum frequency for the modal analysis. Strength and buckling constraints are not considered, since this was not the focus of the problem.

NURBS and mesh creation
The surfaces (shells) analyzed in this study can be modeled with NURBS by using the geometric function F:[0,1] 2 → Ω, as in Equation (1): where R ij (u,v) is the ij-th NURBS basis function, and P ij is a control point in R 3 . This geometric function transforms the unitary square domain, called parametric domain, into the geometric domain Ω.
The NURBS basis functions are defined as in Equation (2): where N i,p and M j,q are B-spline basis functions of p and q degrees, respectively, and w ij is the weight given to each basis function product.
The B-spline basis functions are defined recursively by using Equations (3) and (4):

Introduction
It is not recommended to use NURBS surfaces with non-injective geometric functions, since these functions cause the parametric into geo-metric domain transformation to generate elements with null areas, making the finite element analysis impossible. For more details about NURBS and data regarding the construction of the main surface geometries, see Piegl and Tiller (1996).

Mu (t) + Cu (t) + Ku(t) = P(t)
The dynamic analysis consists of solving Equation (5), called the equation of motion, where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, P is the matrix of applied forces and u is the vector of nodal displacements.
Rayleigh's damping was employed, where the damping matrix C is a linear combination of the stiffness and mass matrices: For the thickness optimization of the shell elements, the optimization problem represented by Equation (6) was considered: The objective of this optimization is to minimize the volume of the structure and, consequently, minimize its weight and value. The design variables are the thicknesses of element groups (they are combined into groups to facilitate the optimization process). The maximum displacement and the minimum frequency of the analysis are considered as constraints. In Equation (6), V is the volume, h is the vector of design variables that contains the thicknesses of all the elements (or group of elements), A is a vector that contains the areas of all the elements, u j is the absolute displacement (i.e. considering the three-dimensional components) or the vertical displacement, as explicit on the example, of the j-th node at time t calculated by Equation (5), u max is the maximum allowable displacement, ω min is the minimum vibration frequency, ω k is the k-th vibration mode, and lb and ub are the vectors that contain the To solve the equation of motion, Newmark's method, with constant acceler-ation, was used. Details about Newmark's method can be found in Newmark (1959), while details about dynamic analysis can be found in Clough and Penzien (2003).  where U=[u 0 ,u 1 ,... ,u n ] is a knot vector with non-decreasing numbers from zero to one. The repetition of a value in the knot vector reduces the derivative of the function in one degree at that point. However, the first and last knots are repeated p+1 times to ensure the partition of unity property. The mesh of the shell can be generated through the image of the mesh of the parametric domain by using the F function, as shown in Figure 1. The mesh of the parametric domain can be easily created by dividing the edges into n by m parts, and each resulting rectangle is subdivided into eight triangles.

Thickness optimization
(4) Optimization of shells modeled using NURBS subjected to a dynamic loading condition  The computational routines were developed in MATLAB and the solution to the optimization problem is obtained using the Matlab function fmincon. The sensitivity analysis is performed by the application of a finite difference method. For the analysis of the shell elements, a combination of a plate element and a membrane element was used, according to Clough and Johnson (1968). The Constant Strain Triangle (CST) membrane element and the Discrete Kirchhoff Theory (DKT) plate element were selected, because they show good results when paired with a moderately large number of elements, and they have an explicit formulation with a low computation cost when used in the optimization problem. The stiffness and mass matrices of the CST element are easily found in introductory books on finite elements, such as Logan (2011), whereas the stiffness and mass matrices of the DKT element can be found in Bathe (1981) and Luo and Hutton (2002), respectively. The dynamic analysis routines were validated by a comparison with the results obtained by ANSYS Workbench 18.2, with which the Shell 181 element was used.
The use of NURBS can reduce considerably the number of elements necessary to perform the analysis in comparison to the traditional way, due to its exact modeling of the analyzed domain, without affecting the results.

Example 1 -Bus stop shelter:
Consider the shell presented in Figure  2a, which simulates a bus stop shelter. The data required for its NURBS modeling is presented in Table 1. The entire shelter was made of steel, with a modulus of elasticity of 200GPa, Poisson's ratio of 0.3 and mass density of 7850kg/m 3 . The bottom of the shelter was fixed, and an applied force F=1000 sin(t) was employed, as shown in Figure  2a. The force was, divided into 100 steps of 0.1 second each. Damping was not used. Figure 2b illustrates the 256-element finite element model that was generated from a parametric domain with 8×4 rectangles subdivided into eight triangles each. Prior to the optimization process, a dynamic analysis was performed considering a thickness of 2cm, and the vertical displacement in one of the applied force nodes was obtained, comparing the author routines (i.e. those of this study) and the ANSYS results, as presented in Figure  3. Based on the displacement obtained, a maximum vertical displacement equal to 0.02m was taken as the reference in the optimization process. Furthermore, the limiting thicknesses of 1cm and 4cm

Results
minimum and maximum thicknesses of each element, respectively.
The solution to the optimization problem was obtained by using Se-quential Quadratic Programming (SQP) (Horowitz, 1999). were adopted. The minimum frequency obtained in the dynamic analysis was 1.36Hz, and it was employed as a limit; in addition, two examples were performed: with and without a self-weight. The minimum allowable thickness established in Section 4 was considered as the initial value in the optimization process with constant thickness, and the calculated optimum thickness was used as the initial value in the optimization process with variable thickness.  In the optimization process without the self-weight and with only one design variable (i.e. uniform thickness), the optimum thickness of 2.05cm was calculated, corresponding to a volume of 0.2446m 3 . When the thicknesses were allowed to vary in horizontal strips, the minimum thickness, maximum thicknesses and volume of 1.19cm, 2.36cm, and 0.2354cm 3 were calculated, respectively. The thickness results are shown in Figure 4a.
When the self-weight was applied, the optimum uniform thickness of 3.03cm was calculated, corresponding to a volume of 0.3604. When allowing the thicknesses to vary, the minimum thickness, maximum thickness and volume of 1.44 cm, 2.95cm, and 0.2579cm 3 were calculated, respectively. The thickness results are shown in Figure 4b.  Optimization of shells modeled using NURBS subjected to a dynamic loading condition The finite element model was created with the aim of constructing a format similar to the actual structure while not necessarily possessing the exact dimensions, as this would require a deeper study of the geometry of the shell. There were 704 elements used in the model. It is worth pointing out the repetition of a few knots in the knot vector to create the discontinuity of the surface derivative, as presented in Table 2   [0; 0; 0;---; ; ; ; ; ; ; ; ; ; 1;1;1] 1 2 3 4 5 6 7 8 9 10 11 11 11 11 11 11 11 11 11 11 The shells were made of concrete with an f ck of 25MPa, modulus of elasticity of 23.8GPa, Poisson ratio of 0.2 and specific mass of 2500kg/m 3 . The bottom ends of the shell were fixed, while the rest of the structure remained free to deform.
In addition to the self-weight of the structure, a dynamic load was distributed over the entire shell in the horizontal direction (i.e. the direction where the perpendicular projection of the area of the shell is not null) with intensity, following Figure 6a.
The load was, divided into 50 steps of 0.2 second each. The horizontal loading condition sought to simulate a strong and short wind gust over the structure. For this example, a damping matrix C =0.05M + 0.05K was also used. To obtain the limiting values for displacement, a uniform thickness of 20cm was first adopted. A dynamic analysis was then performed, and the absolute displacement for the central node of the highest point of the shell was acquired. The displacement over time is shown in Figure  6b, where the maximum value was slightly short of 0.1m (0.09862m, to be exact). Therefore, the maximum displacement restriction was considered to be 0.1m. A frequency-related limit was not used.
For the optimization process, 10cm and 25cm were taken as the minimum and maximum allowable thicknesses, respectively. In the problem involving uniform thickness, this minimum thickness was employed as the initial value, whereas for the problem involving variable thickness, the optimum uniform thickness calculated in the previous process was used as the initial value. Furthermore, to reduce the number of design variables, the thickness variation was considered to be the combination of four parallel horizontal strips of the elements, i.e. there was one different thickness variation variable for every four strips (which differs from the previous example where each horizontal strip was one different design variable). The optimum uniform thickness calculated was to be 23.17cm, which represents a volume of 60.3177m 3 . When the thickness of parallel strips was allowed to vary, the resulting minimum and maximum thicknesses were imposed as constraints and the minimum volume was 47.3094m 3 . The thicknesses calculated are illustrated in Figure 7.  A different approach can be taken, for aesthetic reasons, by defining the thicknesses to allow each "curve" to have its own thickness. This way, there are only four design variables being considered and the minimum volume calculated is 56.9057m 3 .
The resulting thicknesses were, from left to right, 25.00cm, 24.20cm, 20.87cm and 18.87cm, as shown in Figure 8.

References
In Example 1, when the self-weight was not applied, a reduction of 3.76% was obtained when comparing the use of uniform and variable thicknesses. When the self-weight was applied, this reduction was much higher, by 28.44%. In Example 2, reductions of 21.57% and 5.67% were observed in the first and second cases, respectively. The first is more economical, while the second is more interesting from a constructive and aesthetic point of view.
It is evident that as the number of design variables increases, so does the computational cost; however, the optimum calculated volume decreases. It is common, for structural optimization problems using prefabricated parts, to use discrete variables, but in this work because it is a problem of shell optimization, especially in the second example that is a reinforced concrete shell, we preferred to use continuous variables.
When the self-weight is applied, due to the displacement constraints, the optimization process places the thickest (and heaviest) elements close to the constraints and the thinnest (and lightest) close to the cantilever, as observed in both examples. When the self-weight is not used, the distribution of the elements is a function of the stiffness of the structure as a whole, since the variation of thickness does not influence the loading.
The use of NURBS in the modeling of the domain makes it unnecessary to use software to create the mesh, in addition to making the modeling itself directly connected to the discretization of the domain.
The verification via dynamic analysis ensures the reliability of the structure, mainly for those with slim thicknesses. In such cases, the variation of the wind can be more important than the self-weight and should be considered.
It was observed that using the optimization process is an advantage, since it is possible to use a routine to obtain the optimum thickness of the structure, making the structural design process -which is not trivial for shells -easier. Despite the reductions resulting from the optimization process with variable thicknesses, it is not always feasible due to cost issues related to constructive processes and aesthetics. However, it is always possible to adopt an approach similar to the one presented in Example 2 and select which elements or strips should have equal thicknesses, in accordance with the architec-tural project. This approach would also result in savings when compared to the use of a constant thickness throughout the entire structure. As shown above, the minimum volume calculated decreases (or does not change) as the number of design variables increases.
Finally, it is noteworthy that if considered were the constraints that take into account the strength analysis and a buckling analysis, the results of the problem could be different from those obtained here.