Abstract
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.
Keywords:
shells; NURBS; dynamic; optimization
1. Introduction
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 NonUniform Rational BSplines (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)KIENDL, J. et al. Isogeometric shell analysis with KirchhoffLove elements. Computer Methods in Applied Mechanics and Engineering, v. 198, p. 39023914, 2009. also demonstrate shells modeled using NURBS, employing a KirchhoffLove 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)BREITENBERGER, M. et al. Analysis in computer aided design: nonlinear isogeometric BRep analysis of shell structures. Computer Methods in Applied Mechanics and Engineering, v. 284, p. 401457, 2015. present a nonlinear analysis of a computeraided shell design modeled with NURBS via boundary representation, also using KirchhoffLove shells and isogeometric analysis. Lastly, Zhang et al. (2017)ZHANG, X. et al. NURBS modeling and isogeometric shell analysis for complex tubular engineering structures. Computational and Applied Mathematics, v. 36, n. 4, p. 16591679, 2017. show an analysis of complex tubular shells modeled by isogeometric shell analysis using degenerated NURBS elements and the MindlinReissner theory. To verify the efficiency of this method, they constructed a ComputerAided 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)FALCO, S. A.; AFONSO, S. M. B.; VAZ, L. E. Analysis and optimal design of plates and shells under dynamic loadsI: finite element and sensitivity analysis. Structural and Multidisciplinary Optimization, v. 27, n. 3, p. 189196, 2004a. present the development and application of a computational tool for geometry modeling, mesh generation, structural analysis (a HuangHinton element, based on a MindlinReissner 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)FALCO, S. A.; AFONSO, S. M. B.; VAZ, L. E. Analysis and optimal design of plates and shells under dynamic loadsII: optimization. Structural and Multidisciplinary Optimization, v. 27, n. 3, p. 197209, 2004b.. Park and Dang (2010)PARK, H. S.; DANG, X. P. Structural optimization based on CADCAE integration and metamodeling techniques. ComputerAided Design, v. 42, n. 10, p. 889902, 2010. present a structural optimization method based on ComputerAidedDesign and ComputerAidedEngineering integration, using metamodeling techniques to reduce the time taken for solving the problem and to automatize the process. Alves and Vaz (2013)ALVES, E. C.; VAZ, L. E. Optimum design of plates structures under random loadings. REM  Revista Escola de Minas, v. 66, n. 1, p. 4147, 2013. present the optimization of plates subject to random dynamic loading. Espath et al. (2011)ESPATH, L. F. R.; LINN, R. V.; AWRUCH, A. M. Shape optimization of shell structures based on NURBS description using automatic differentiation. International Journal for Numerical Methods in Engineering, v. 88, n. 7, p. 613636, 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)KANG, P.; YOUN, S. K. Isogeometric topology optimization of shell structures using trimmed NURBS surfaces. Finite Elements in Analysis and Design, v. 120, p. 1840, 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.
2. 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 ijth 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 Bspline basis functions of p and q degrees, respectively, and w_{ij} is the weight given to each basis function product.
The Bspline basis functions are defined recursively by using Equations (3) and (4):
where U=[u _{0},u _{1},... ,u_{n} ] is a knot vector with nondecreasing 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.
The transformation of the parametric domain into the geometric domain. a) Parametric domain. b) Geometric domain.
It is not recommended to use NURBS surfaces with noninjective geometric functions, since these functions cause the parametric into geometric 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)PIEGL, L.; TILLER, W. The NURBS book. [S.l.]: Springer Science & Bussiness Media. 1996. 646p..
3. Dynamic analysis
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:
To solve the equation of motion, Newmark's method, with constant acceleration, was used. Details about Newmark's method can be found in Newmark (1959)NEWMARK, N. A Method of Computation of Structural Dynamics. ASCE Journal of the Engineering Mechanics Division, v. 85, p. 6794, 1959., while details about dynamic analysis can be found in Clough and Penzien (2003)CLOUGH, R. W.; PENZIEN, J. Dynamics of structures. [S.l.]: Computer & Structures, Inc. 2003. 730p..
4. Thickness optimization
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 threedimensional components) or the vertical displacement, as explicit on the example, of the jth node at time t calculated by Equation (5), u_{max} is the maximum allowable displacement, ω_{min} is the minimum vibration frequency, ω_{k} is the kth vibration mode, and lb and ub are the vectors that contain the minimum and maximum thicknesses of each element, respectively.
The solution to the optimization problem was obtained by using Sequential Quadratic Programming (SQP) (Horowitz, 1999HOROWITZ, B. Implementation Considerations in SQP ALgorithms for LargeScale Structural Optimization. In: IBERIAN LATINAMERICAN CONGRESS ON COMPUTATIONAL METHODS IN ENGINEERING (CILAMCE), 20., 1999, São Paulo. Proceedings [...]. São Paulo: Departamento de Engenharia de Estruturas e Fundações  EPUSP, 1999. p. 93.0193.10.).
5. Results
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)CLOUGH, R. W.; JOHNSON, C. P. A finite element approximation for the analysis of thin shells. International Journal of Solids and Structures, v. 4, p. 4360, 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)LOGAN, D. L. A first course in the finite element method. [S.l.]: Cengage Learning, 2011. 808p., whereas the stiffness and mass matrices of the DKT element can be found in Bathe (1981)BATHE, K. J.; HO, L. W. A simple and effective element for analysis of general shell structures. Computers & Structures, v. 13, p. 673681, 1981. and Luo and Hutton (2002)LUO, Z.; HUTTON, S. G. Formulation of a threenode traveling triangular plate element subjected to gyroscopic and inplane forces. Computers & Structures, v. 80, n. 26, p. 19351944, 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 256element 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 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 selfweight. 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 selfweight 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 selfweight 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.
Example 2  Igrejinha da Pampulha: The "Igrejinha da Pampulha" (São Francisco de Assis Church), located in Belo Horizonte, in the state of Minas Gerais, Brazil, was designed by Oscar Niemeyer and inaugurated in 1943. The building distinguishes itself by the parabolic shells that comprise its architecture, as can be seen in Figure 5a.
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 and Figure 5b. Gomes Filho et. al. (2018)GOMES FILHO, H.; ALVES, E. C.; GONÇALVES JÚNIOR, E. Otimização de estruturas de cascas com domínio modelado por meio de NURBS. In: SIMPÓSIO DE MECÂNICA COMPUTACIONAL (SIMMEC), 13., 2018. Vitória. Anais [...]. Vitória: SIMMEC, 2018. present a similar example, which employs only a static analysis.
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 selfweight 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.
a) The horizontal load applied to the model. b) Absolute displacement obtained at the central node of the highest point of the shell.
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 frequencyrelated 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.
Thicknesses calculated in the optimization process  groups based on geometry and aesthetic aspect.
6. Discussion
In Example 1, when the selfweight was not applied, a reduction of 3.76% was obtained when comparing the use of uniform and variable thicknesses. When the selfweight 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 selfweight 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 selfweight 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.
7. Conclusions
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 selfweight 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 architectural 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.
References
 ALVES, E. C.; VAZ, L. E. Optimum design of plates structures under random loadings. REM  Revista Escola de Minas, v. 66, n. 1, p. 4147, 2013.
 ALVES, E. C.; VAZ, L. E. Optimization of structures subject to dynamic load: deterministic and probabilistic methods. REM  Revista Escola de Minas, v. 69, n. 3, p.281286, 2016.
 BATHE, K. J.; HO, L. W. A simple and effective element for analysis of general shell structures. Computers & Structures, v. 13, p. 673681, 1981.
 BREITENBERGER, M. et al. Analysis in computer aided design: nonlinear isogeometric BRep analysis of shell structures. Computer Methods in Applied Mechanics and Engineering, v. 284, p. 401457, 2015.
 CLOUGH, R. W.; JOHNSON, C. P. A finite element approximation for the analysis of thin shells. International Journal of Solids and Structures, v. 4, p. 4360, 1968.
 CLOUGH, R. W.; PENZIEN, J. Dynamics of structures [S.l.]: Computer & Structures, Inc. 2003. 730p.
 ESPATH, L. F. R.; LINN, R. V.; AWRUCH, A. M. Shape optimization of shell structures based on NURBS description using automatic differentiation. International Journal for Numerical Methods in Engineering, v. 88, n. 7, p. 613636, 2011.
 FALCO, S. A.; AFONSO, S. M. B.; VAZ, L. E. Analysis and optimal design of plates and shells under dynamic loadsI: finite element and sensitivity analysis. Structural and Multidisciplinary Optimization, v. 27, n. 3, p. 189196, 2004a.
 FALCO, S. A.; AFONSO, S. M. B.; VAZ, L. E. Analysis and optimal design of plates and shells under dynamic loadsII: optimization. Structural and Multidisciplinary Optimization, v. 27, n. 3, p. 197209, 2004b.
 GOMES FILHO, H.; ALVES, E. C.; GONÇALVES JÚNIOR, E. Otimização de estruturas de cascas com domínio modelado por meio de NURBS. In: SIMPÓSIO DE MECÂNICA COMPUTACIONAL (SIMMEC), 13., 2018. Vitória. Anais [...] Vitória: SIMMEC, 2018.
 HOROWITZ, B. Implementation Considerations in SQP ALgorithms for LargeScale Structural Optimization. In: IBERIAN LATINAMERICAN CONGRESS ON COMPUTATIONAL METHODS IN ENGINEERING (CILAMCE), 20., 1999, São Paulo. Proceedings [...] São Paulo: Departamento de Engenharia de Estruturas e Fundações  EPUSP, 1999. p. 93.0193.10.
 KANG, P.; YOUN, S. K. Isogeometric topology optimization of shell structures using trimmed NURBS surfaces. Finite Elements in Analysis and Design, v. 120, p. 1840, 2016.
 KIENDL, J. et al. Isogeometric shell analysis with KirchhoffLove elements. Computer Methods in Applied Mechanics and Engineering, v. 198, p. 39023914, 2009.
 LOGAN, D. L. A first course in the finite element method [S.l.]: Cengage Learning, 2011. 808p.
 LUO, Z.; HUTTON, S. G. Formulation of a threenode traveling triangular plate element subjected to gyroscopic and inplane forces. Computers & Structures, v. 80, n. 26, p. 19351944, 2002.
 NEWMARK, N. A Method of Computation of Structural Dynamics. ASCE Journal of the Engineering Mechanics Division, v. 85, p. 6794, 1959.
 PARK, H. S.; DANG, X. P. Structural optimization based on CADCAE integration and metamodeling techniques. ComputerAided Design, v. 42, n. 10, p. 889902, 2010.
 PIEGL, L.; TILLER, W. The NURBS book [S.l.]: Springer Science & Bussiness Media. 1996. 646p.
 ZHANG, X. et al. NURBS modeling and isogeometric shell analysis for complex tubular engineering structures. Computational and Applied Mathematics, v. 36, n. 4, p. 16591679, 2017.
Publication Dates

Publication in this collection
20 Dec 2019 
Date of issue
JanMar 2020
History

Received
23 Mar 2019 
Accepted
12 Sept 2019