Modern software capabilities for shape optimization of shells

. The article is devoted to the shape optimization of shell structures in Comsol Multiphysics using three gradient-based methods: IPOPT (Interior Point OPTimizer), SNOPT (Sparse Nonlinear OPTimizer) and MMA or GCMMA (the Method of Moving Asymptotes). Two types of complex shapes, such as right helicoid and developable helicoid are taken for the computational experiment. The task is to investigate the initial design and optimization process of two helicoids. To obtain a more accurate result and an interesting design solution, the calculation is carried out using three physics-controlled mesh sizes: extra coarse, fine and extra fine with varying values of special optimization settings, such as maximum displacement ( d max ) and filter radius ( R min ). The results obtained using these three methods allow to conclude that the studied parameters d max and R min have a significant impact on the final optimization result.


INTRODUCTION
Shape optimization is a powerful tool for efficient structural design of buildings and structures.Providing designers with a universal set of tools for modeling and calculation contributes to solving a wide range of optimization problems of shells, especially complex shapes.Currently, in many calculation programs, there is a special module for optimizing the shape and it is possible to quickly get the best design option.Among the most used shape optimization modules are Workbench-Structural optimization (Ansys Mechanical) [1], Simulia Tosca Structure Shape (Simulia Abaqus) [2] and Optimization-Shape Optimization (Comsol Multi-physics) [3].Also, optimization is possible in Rhino 3D modeling program in connection with the Grasshopper plugin and special applications, such as Galapagos, Millipede and so on [4].In this work all necessary calculations are carried out in the Comsol Multiphysics software package, providing wide opportunities for modeling and subsequent calculation of shells of simple and complex shapes.
Program modules contain a certain set of methods for solving any optimization problem, thereby expanding the choice of the most appropriate solution.In addition to the traditional preparatory stage with material, load and boundary condition settings, special attention is also paid to constructing an appropriate mesh, which can significantly affect the calculation result.For example, in articles [5 -8] the effect of different mesh sizes during the computational process is demonstrated.
After the structural analysis, the optimization settings are installed directly and the most interesting section in the Shape optimization module in Comsol is the control variable settingmaximum displacement (d max ) and filtering settingfilter radius (R min ), which are situated in Free Shape Boundary feature and described in the present work.All these settings are specified by the user and are responsible for changing the shape during the optimization process.The control variable (d max ) controls the amount of optimization displacement and filter radius (R min ) controls the smoothness of the computational mesh.The filtering method is widely used in topological optimization [9 -10], but it can also be applied to shape optimization tasks [11 -12].It is necessary to mention that no works have been found in the public domain that describe the control variable application and the general process of using d max and R min together in the Comsol Multiphysics environment, which confirms the novelty of the current study.Also, the significant advantage of the work is the proposed recommendations for the selection of suitable values d max and R min for shells, using the example of the considered helicoids.
The purpose of this work is to explore the capabilities of Comsol Multiphysics to provide optimization of two complex shapes with different values of maximum displacement (d max ) and filter radius (R min ) using multiple mesh sizes and optimization methods.Two types of helicoid shells are taken as the test models for calculation.Thus, it is necessary to solve the following tasks: -to give a brief description of shells under study and to provide a structural analysis to analyze the stress-strain state of the shells with helical surfaces; -to perform optimization calculations with d max and R min using three gradient-based methods and three mesh sizes from extra coarse to extra fine; -to compare the obtained results.

Description of the shells under study
The surface of the helicoid is obtained by rotating a straight line with a constant angular velocity around a vertical fixed axis (a helical surface), intersecting the axis at a constant angle, and simultaneously moving translationally with a constant velocity along this axis.The general parametric equations for the helicoid can be written in the following form [13]: ( , ) cos sin sin ( , ) sin sin cos ( , ) cos where parameters u, v are the curvilinear coordinates and define the position of an arbitrary point of the helicoid; a-the least distance between the straight line and the helical axis Oz;the angle of the generatrix straight line with the helical axis Oz; pthe pitch of the screw.
Helical surfaces are characterized by unusual architectural properties, since this type of shell can be of various shapes, forming the final volume at the base.Helicoids are widely used in the design of suspended structures, small architectural forms, spiral ramps, staircases, and so on (Figure 1).The considered types of helicoids can be obtained, depending on the shape of the generatrix and the position of the axis of rotation.The right helicoid is considered a minimal surface, it has mean curvature of zero [14].The one loop of the right helicoid has parametric equations of the following form: 12 ( , ) cos( ), ( , ) sin( ), ( , ) 02 Parameters u and v are the curvilinear coordinates of an arbitrary point of the helicoid, c is the vertical displacement of the generatrix, which is the segment from u 1 to u 2 on the axis upon its rotation by 1 radian and is related to the height of 1 loop of the helicoid H (with the height of the surface, the generatrix of which has rotated 360°) by the following equation: c = H/2π [15] (Fig. 2).
A developable helicoid is a ruled surface formed by the tangent lines to the helix of constant slope on a circular cylinder with radius a [14].The one loop of the developable helicoid has parametric equations in the following form: sin cos ( , ) cos , ( , ) sin , ( , ) 02 Here, the curvilinear coordinate u (or v = const) is a rectilinear generatrix, tangent to the helical cuspidal edge, but lines v (u = const) are helixes [16].The surface of the developable helicoid is bounded by the cuspidal edge √ , b is the lead of a helix u = 0, v is an angle measured from an axis Ox, a is a radius of a cylinder on which the helical cuspidal edge is laid [13] (Figure 2).All above-mentioned shells are created in the Comsol Multiphysics program with the following values: inner radius -3.1 m, outer radius -7.1 m, and the height of one loop -3 m.The width of the one helicoid loop is 4 m and the total thickness of the shell is 0.16 m.The selected material for the shells is concrete B25 with a density (D) -2300 kg/m 3 , Young's modulus (Е) -30 × 10 9 Pa and Poisson's ratio (μ) -0.2.
Before performing optimization, it is necessary to analyze the design by the finite element method to determine the deformation of the structure.At the preparatory stage, the models are divided into triangle elements and mesh with the following selected element sizes -«extremely fine», «fine» and «extremely coarse».These types of meshes are physics-controlled since the geometry does not have complex outlines and there is no need to rebuild manually.In a right helicoid, «extremely fine» mesh consists of 1324 triangle elements, «fine» -264 triangle elements, and «extremely coarse» -30 triangle elements.In a developable helicoid, «extremely fine» mesh consists of 1360 triangles, «fine» -259 triangles, and «extremely coarse» -32 triangles.Also, the necessary restrictions are setinitial, final faces, and inner radius of one helicoid loop are fixed (Figure 3).A vertical uniform distributed force F = 10000 N/m 2 acts on the surface of the shell.Then the structural analysis is carried out for two shells.The optimization process after structural analysis is performed in the Shape optimization module with the use of three gradient-based methods -IPOPT (Interior Point OPTimizer), SNOPT (Sparse Nonlinear OPTimizer), and MMA or GCMMA (the Method of Moving Asymptotes), which are widely used for structural optimization problems.These methods contain three basic steps: finding the direction of movement to achieve the optimal result, determining the step dimension, and checking convergence.

IPOPT (Interior Point OPTimizer
) is an open-source solver, which implements an interiorpoint line-search filter method [17].The approach is based on choosing a starting point, considering the curvature of the objective function in the vicinity of this starting point, and tracking the curvature until they find a point of local optimum.
MMA or GCMMA (the Method of Moving Asymptotes) is a general method, which is based on the problem of approximation of convex functions and, at each iterative stage, a strictly convex approximating sub-problem is solved [3].This sub-problem is controlled by adaptive mechanism, which is called «moving asymptotes», which move in the direction of the search and change with each iteration, gradually approaching the optimal solution.SNOPT (Sparse Nonlinear OPTimizer) is a solver, which is based on the use of the sequential quadratic programming (SQP) method [18].This means that SNOPT solves a sequence of approximations to the original problem, where the objective function is quadratic with several variables under linear constraints on these variables [3].The steps in this sequence are called major and minor iterations.
As the objective function criterion of the system, the total elastic strain energy is selected, and the task is to minimize it.Fifty iterations are set for each method.

Structural analysis of the right helicoid with three types of meshes
After the preparatory analysis, the stress-strain state of the model of the right helicoid is calculated with the use of three types of meshes: «extra fine», «fine», and «extra coarse».The graphical results are shown in Figure 4. «Extra fine» mesh: The maximum total displacement value is 31.7 mm.The maximum value of the first principal stress is 4.75 MPa.The maximum value of the second principal stress is 0.18 MPa.The value of the third principal stress is -25.2MPa (min).Total elastic strain energy before the optimization is 7.74 kJ.
«Fine» mesh: The maximum total displacement value is 31.8mm.For the first principal stress, the maximum value is 5.02 MPa (max).The second principal stress is 0.35 MPa (max).The value of the third principal stress is -23.9MPa (min).Total elastic strain energy before the optimization is 7.78 kJ.

«Extra coarse» mesh:
The maximum total displacement value is 31.2mm.The maximum value of the first principal stress is 7.31 MPa (max).The second principal stress is 1.56 MPa (max).The value of the third principal stress is -25.4MPa (min).Total elastic strain energy before the optimization is 7.07 kJ.

Structural analysis of the developable helicoid
Furthermore, the stress-strain state analysis is performed for the developable helicoid shell.The obtained results are demonstrated in Figure 5. «Extra fine» mesh: The maximum total displacement value is 10.5 mm.For the first principal stress, the maximum value is 4.28 MPa (max).The second principal stress is 0.37 MPa (max).For the third principal stress, the minimum value is -14.9MPa.The value of the total elastic strain energy is 2.55 kJ.
«Fine» mesh: The maximum total displacement value is 10.3 mm.For the first principal stress, the maximum value is 3.66 MPa (max).The second principal stress is 0.268 MPa (max).
The value of the third principal stress is -17.7 MPa (min).Total elastic strain energy before optimization is 2.52 kJ.

«Extra coarse» mesh:
The maximum total displacement value is 9.96 mm.The value of the first principal stress is 3.51 MPa (max).The second principal stress is 0.678 MPa (max).For the third principal stress, the minimum value is -23.7 MPa.Total elastic strain energy before optimization is 2.41 kJ.

The optimization process after structural analysis
After the stress-strain analysis is carried out, the Shape optimization module is launched and its settings are switched to manage the existing control variable settingsmaximum displacement (d max ) and filteringfilter radius (R min ), which are part of the Free Shape Boundary feature.The SI unit for these settings is m.By default, the following values are set in the program: 5 % of the geometry bounding box (BBox) for maximum displacement and 10 % of the geometry bounding box (BBox) for filter radius [19].For two helicoids, the program automatically installed d max = 1 m and R min = 2 m, which are added to the computational process.

Initial design of two shells after structural analysis
The obtained results show that the largest displacements in right and developable helicoid shells are marked along the outer radius.The maximum first principal stress values here are concentrated closer to the base of the helicoid.The lowest values of the second stress and total elastic strain energy are located along the inner radius of the two shells.The lowest distribution of the third principal stress in the developable shell is concentrated at the beginning and end of the loop.In the right helicoid shell, the lowest values are marked mainly along the inner radius.It should also be noted that with «extra fine» and «fine» meshes a smoother distribution of stresses and total elastic strain energy can be observed.This is especially evident on the shell of a right helicoid.Furthermore, in the developable helicoid with «extra coarse» mesh, the stress distribution is more uniform than in the right helicoid.
In addition, the values obtained with «extra fine» and «fine» mesh sizes have an almost insignificant difference in both the right helicoid and the developable one.Also, when analyzing the stresses arising in the right helicoid, it can be noticed that as the mesh size decreases, the stress at the same location increases.Considering all of the above, a mesh with extra fine element size is taken for further optimization calculations.

Optimized design of two shells
Total displacement: During the optimization process with the use of extra fine mesh size, the smallest displacement values are obtained with d max = 1 m, R min = 2 m and d max = 0.36 m, R min = 0.71 m.The displacement values are reduced by 94 -98 % for the right helicoid shell and 88 -92 % for the developable helicoid shell.The results obtained using the three methods are almost the same for the developable helicoid, while the right one has a noticeable difference between MMA and the remaining two methods.In the right helicoid the smallest values are obtained using the MMA method.In the case of the developable helicoid, the smallest results are obtained using the SNOPT method.The graphical results are demonstrated in Figure 6.   ), the helicoid changes within acceptable limits, the structure is not deformed, and the mesh is smoother and has no irregularities.When setting these parameters, the upper part of the right helicoid is raised and does not bend down as with the other values.
Total elastic strain energy: After optimization, a significant decrease in the objective function can be observed in both shells.The maximum decrease in the optimized total elastic strain energy function is noted during setting d max = 1 m and R min = 2 m.The function value is reduced by 90 -96 % compared to the original result before optimization.In other cases, the decrease in the function values is 70 -80 %.

Figure 1 .
Figure 1.Application of helical surfaces in architecture.

Figure 3 .
Figure 3. Preparatory stage for the calculation.

Figure 4 .
Figure 4. Structural analysis of the right helicoid shell.
Also, the values with the above-mentioned percentage and intermediate values for more extended results are taken: d max = 0.15 m and R min = 0.3 m; d max = 0.36 m and R min = 0.71 m; d max = 0.5 m and R min = 5 m; d max = 1 m and R min = 5 m.

Figure 6 .
Figure 6.Total displacement and Total elastic strain energy for different values of d max and R min (right and developable helicoid shells).Principal stress:The smallest values of the first principal stress are obtained using the MMA method.When analyzing the second principal stress, a slight increase in values can be seen, especially if the parameters, such as d max = 0.15 m, R min =0.3 m, d max = 0.36 m, R min = 0.71 m and d max = 1 m, R min = 2 m are set.As for the third principal stress, the values are negligibly small (Figure7).

Figure 7 .
Figure 7. First, Second, Third stress for different values of d max and R min (right and developable helicoid shells).

Figure 8 .
Figure 8. Normal boundary displacement of the right and developable helicoid shells.Normal boundary displacement: The highest percentage of the normal boundary displacement in the shells can be observed with d max = 0.15 m, R min = 0.3 m and d max = 0.36 m, R min = 0.71 m.As can be seen from Figure 8, the shape of the helicoid becomes rough and bumpy.The same situation is formed with d max = 1 m, R min = 2 m, where the shell of the right helicoid begins to bend strongly.With the remaining combinations (d max = 0.5 m, R min = 5 m and d max = 1 m, R min = 5 m), the helicoid changes within acceptable limits, the structure is not deformed, and the mesh is smoother and has no irregularities.When setting these parameters, the upper part of the right helicoid is raised and does not bend down as with the other values.