Shape optimization of concrete free-form shells considering material damage

Shape optimization is a pervasive tool for designing concrete free-form shells. However, the existing shape optimization approaches for free-form shells usually assume the material is purely elastic without considering material damage, which may produce undesired results in real-world structural applications, especially in the case of concrete structures. In this article, a framework is presented for performing the shape optimization of concrete free-form shells while considering concrete damage. In the proposed optimization algorithm, the Mazars model is coupled with finite element analysis for modelling structural responses subjected to material damage. Through a series of numerical examples, the proposed algorithm is validated and the effects of material damage on optimal results are investigated. It is demonstrated that the consideration of material damage leads to a more robust design, which highlights the importance of accounting for the material damage during the shape optimization process of concrete free-form shells.


Introduction
Different from regular shells, free-form concrete shells have unique aesthetic configurations that usually cannot be described by analytical formulas (Feng and Ge 2013;Garcia et al. 2019). Over the past few years, the advancement of construction techniques such as 3-D printing motivates their everincreasing application in civil engineering (Su et al. 2018), such as the roofs of long-span structures (e.g. Island City Central Park Gringrin 2005, the Toshima Art Museum 2010 and the Rolex Learning Center 2010) and dry cooling towers. When designing free-form shells, shape optimization is an indispensable step, since their mechanical properties are strongly tied to shape. The pioneering work on the shape optimization of shells started from experimental approaches, such as the inverse hoisting method, the inflatable films method (Isler 1993;Motro and Oliva 2010), etc. However, when the structure is large and geometrically complicated, experimental methods are difficult to implement owing to the constraints of similarity laws and measuring accuracy. In recent decades, numerical shape optimization methods (Gholizadeh and Barzegar 2013;Cao et al. 2017) have been developed and have become beneficial tools for designing free-form structures.
Via sensitivity weighting and a semi-analytical approach, Bletzinger et al. (2005) presented an isogeometric shape optimization method for shells. Cui and Yan (2006) developed a height adjusting method for shape optimization to minimize the strain energy of free-form shells. Ohmori and Hamada (2006) optimized free-form shells for minimum strain energy and minimum geometrical deviation by coupling non-uniform rational B-splines (NURBS) with a gradient method and a genetic algorithm (GA). Winslow, Pellegrino, and Sharma (2010) employed a GA to optimize rod directions with respect to multiple objectives for free-form structures. To achieve target vibration response, Shimoda and  solved the optimal shapes of free-form structures using a non-parametric algorithm. In another of their works (Y. Liu and Shimoda 2015), an effective shape optimization method was proposed for the natural vibration design of stiffened thin-walled or shell structures. Su et al. (2018) integrated a numerical inverse hanging method and multi-objective optimization to generate free-form shapes of shells. Wang and Wu (2018) proposed a shape optimization method for cable-stiffened latticed shells to improve their geometrical smoothness. Tian et al. (2020) proposed a data-driven modelling and optimization framework for undevelopable stiffened curved shells.
Via the existing work, it can be seen that, in most shape optimization approaches and studies, the material is assumed to be purely elastic, without considering material damage. This represents a void within the literature since material damage is inevitable in real engineering, particularly for concrete structures. The neglect of material damage in shape optimization may produce undesired structures that lack robustness. Herein, it would be interesting to point out the relation with topology optimization, which has employed material damage models in the optimization framework. Achtziger and Bendsoe (1995) applied an internal damage parameter to the topology design of truss structures. B. Desmorat and R. Desmorat (2008) maximized structural lifetime subject to fatigue damage by introducing continuum damage to topology optimization. Amir and Sigmund (2013) optimized a reinforced concrete continuum subject to damage. Luo and Kang (2013) developed a method that can automatically generate the optimal layout of reinforced concrete structures, in which the Drucker-Prager yield criterion is applied to predict the failure behaviour of concrete. In other of their work (P. Liu, Luo, and Kang 2016;Luo and Bao 2019;Luo, Xing, and Kang 2020), material damage or separation were also considered in topology optimization. Waisman (2014, 2015) performed topology optimization under damage constraints to achieve the desired damage resistance of the structure. The above studies show the significant effect of material damage on optimal topology, which further motives the current article to consider material damage in the optimal shape design of concrete shells.
In this study, a framework is developed for the shape optimization of free-form concrete shells accounting for material damage. The Mazars model is combined with a finite element (FE) analysis algorithm to model structural response subject to damage. A structural response quantity, strain energy, is determined as the optimization objective. The geometrical shape of free-form shells is represented by NURBS functions and then optimized by a gradient-based method. Three optimization examples are carried out to investigate the effects of material damage on the optimal results and highlight the importance of accounting for damage effects in the shape optimization of concrete shells.

Parameterization of free-form shapes
In the shape optimization of free-form shells, shape parameterization is essential for describing arbitrary shell shapes. NURBS are known as one of widely-used parameterization approaches and have been applied in FE discretizations, computer aided design and so on (Kapoor and Kapania 2012;Shojaee et al. 2012;Y. Xia, Wu, and Hendriks 2019). Therefore, the current study selects NURBS to parameterize free-form shapes.
For brevity, only major NURBS equations are presented in this section. The reader can refer to Piegl and Tiller (2012) and Azegami, Fukumoto, and Aoyama (2013) for more details.
NURBS surfaces (see Figure 1) are defined as where S (u, v) is the value of the point (u, v) on the NURBS surface, u and v represent the parametric coordinates in the space [0, 1], P i,j are the coordinates of control points, which are located along the directions u and v with the respective numbers m and n, w i,j are the weights of control points, and N i,k (u) and N j,g (v) represent the B-spline basis functions with degrees k and g, respectively. N i,k (u) is given in Equation (2) and N i,g (v) is defined in a similar way.

Coupled material damage model and FE analysis
Damage may refer to the formation of material voids or cracks at the microscopic level. In the early work of Mazzucco et al. (2016), an internal damage parameter D was employed to represent the extent of local damage and model the process at the macro-scale. D ranges between zero and one. The value zero represents no damage and one represents completely damaged. Phenomenological damage laws are adopted to model the damage field and predict the damage initiation and propagation (James and Waisman 2015). In this framework, constitutive relationship involving material damaged can be written as where C ijkl represents the material 4th-order elasticity tensor, and σ ij and ε kl represent the 2nd-order stress and strain tensors, respectively. From Equation (3), it can be seen that the damage is isotropic and described by scalar values, which leads to an easy combination of damage and the FE model. Herein, the Mazars damage law (Mazars and Pijaudier-Cabot 1989) is employed, which was developed to describe concrete brittle damage and is relatively easy to introduce into the FE code. Using the Mazars damage law, the damage parameter D is defined in terms of equivalent strain:ε where where ε i denotes the principle strains. The equivalent strain function indicates that only tensile strains can initiate damage since all negative principle strains are neglected. It agrees with the fact that the compression capacity of concrete is much higher than its tensile capacity. Herein, the damage parameter is defined for three different stress states as follows.
In the case of uniaxial tensile stress, the damage parameter function is denoted as D T and written as where the equivalent strainε is obtained by Equation (4), ε f is the damage threshold, i.e., when local strains exceed ε f , damage is initiated. The damage threshold and constants A T and B T are coefficients that are determined by experiments. The damage parameter for the uniaxial compression state, denoted as D c , is written similarly to Equation (6) except for the values of coefficients: where A C and B C are material coefficients determined through experimental data. For complex stress states, the coupling effect between tensile and compressive damage is considered in the following way (Kotronis et al. 2005): where the damage parameter D is the total weight of D T and D C (see Equations (6) and (7)), with respect to weight coefficients α T and α C , which are solved as where α T + α C = 1, H i is the Heaviside function (Ventura and Benvenuti 2015), and ε Ti and ε Ci are the principle tensile strain and the principle compressive strain, respectively. Via the above equations, the damage parameter can be solved and then coupled with the FE model to simulate material damage and solve the mechanical response.
In this study, FE modelling is based on Kirchhoff-Love plate theory (Reddy 2007). Considering the compromise between surface smoothness and computational cost, a plane triangular shell element with three nodes is selected, which can be deemed as the combination of a plane stress element and a plate bending element. The elements are assumed to be plane concrete. The damage parameters of the upper and lower surfaces of the element are solved with respect to their respective strain states. Then the maximum one is selected to represent the damage of the element and to modify the constitutive relation in Equation (3). The introduction of material damage leads to material nonlinearity in the FE analysis, which is solved by the Newton-Raphson repeat method. The FE algorithm is written in Fortran.

Optimization formulation and algorithm
In this study, the authors optimize the shell shape for minimum strain energy under the given loads. From Equations (1) and (2), it can be seen that NURBS surfaces are determined by both the locations of the control points and their weights. However, considering the tradeoff between computational cost and accuracy, only the parameters that have significant effect on strain energy would be determined as optimization variables in this study. The exiting research (Li 2011) showed that, for free-form shells as shown in Figure 1, the Z-coordinates of control points play a more dominant role in the strain energy compared with their X-, Y-coordinates and weights. In addition, it was found that the update of X-, Y-coordinates during the optimization tends to cause mesh distortion. Therefore, the Z-coordinates of the control points are selected as the optimization variables while keeping the other parameters constant during the optimization. Different from the previous work, the strain energy is solved by a nonlinear FE model accounting for material damage as presented in Section 3. In this way, the effect of material damage is involved in the optimization process. Combining with the NURBS representation, the proposed optimization problem is formulated as searching for the optimal locations of the control points and is written as where the objective function f, i.e. the total structural strain energy, is equal to the sum of the external work of each load increment in the Newton-Raphson approach, P i is the external nodal force vector, U i is the ith increment of the nodal displacement vector subject to material damage parameter D, and nl is the total number of load increments; Z 1 , Z 2 , . . . , Z n×m are the global Z-coordinates of the NURBS control points (see Figure 1); and denotes the feasible space of control points.
The gradient-based method (San, Zhao, and Qiu 2019;Zhang et al. 2019;San, Waisman, and Harari 2020) is employed to solve the proposed optimization problem. It has been proved to be an effective method for shape optimization problems. In the gradient-based approach, sensitivity is needed for determining the searching direction. In the current problem, the presence of material nonlinearity makes it difficult to derive an analytical or semi-analytical expression for the sensitivities. Therefore, the numerical differentiation method (Hay, Borggaard, and Pelletier 2009;Akhtar, Borggaard, and Hay 2010;Noel, Van-Miegroet, and Duysinx 2016) is adopted as an alternative.
Combining the damage-considered FE analysis with the shell shape parameterized by NURBS, the gradient-based optimization procedure is obtained and described in detail in Algorithm 1. The whole algorithm is self-coded using the Fortran language. The subroutine of gradient-based optimization in the algorithm has been verified in the authors' previous work (San, Waisman, and Harari 2020). The validation of the proposed damage-considered FE model by comparing the proposed results with the experimental results in Cui, Cui, and Jiang (2014) can be found in the online Supplemental Data for this article, which can be accessed at https://doi.org/10.1080/0305215X.2021.1968854.

Numerical examples
In the context of concrete roof applications, three shape optimization examples are presented to investigate the effect of material damage and highlight the importance of considering damage in the shape Algorithm 1 Shape optimization algorithm considering material damage parameterized by NURBS functions and using nonlinear FEM (1) Input material property parameters, damage coefficients, shell thickness, load and boundary conditions (2) Assign initial guesses for design variables: Z i = Z i0 (i = 1, 2, . . . , n × m) (3) Set convergence criteria, i.e. the relative change of the objective function and constraint is smaller than a given tolerance (e.g. = 0.1%) (4) Loop until converging (a) Generate the NURBS surface (i) Obtain the knot vectors: = {u 1 , u 2 , . . . , u m+k+1 }, = {v 1 , v 2 , . . . , v n+g+1 } (ii) Input the number of interpolation points. In this work, finite element nodes are set up as interpolation points (iii) Via Equations (1) and (2), the global coordinates of interpolation points (X, Y, Z) are solved (b) Generate FE mesh with surface parameterized by NURBS (c) Solve the nonlinear FE problem accounting for material damage, as presented in Section 3 (d) Compute the strain energy as objective function f in Equation (10) (e) Compute the sensitivity of objective f with respect to Z i where ε zi is a given small value (f) Compute the step size of Z i at the lth load step, denoted as δ (l) Z i , by the golden section method (Koupaei, Hosseini, and Ghaini 2016) (g) Modify the optimization variables: optimization. In the first example, a free-form shell with a continuous surface is optimized with respect to different magnitudes of the applied load; in the second example, a free-form shell with an opening is optimized; in the third example, a free-form shell with weaker boundary constraints than those of Examples 5.1 and 5.2 is optimized. Through the examples, various cases are discussed: different magnitudes of the applied load, the presence or absence of an opening in the shell, and strong or weak boundary conditions. As described in Section 4, the objective of shape optimization is to minimize the strain energy. Algorithm 1 is implemented to solve the optimization problems with the consideration of damage. In all the examples, the material is assumed to be plain concrete with elastic modulus E c = 3.00 × 10 4 MPa and Poisson's ratio μ = 0.167. The concrete damage coefficients in Equations (6) and (7) are determined as A T = 0.75, B T = 1.0 × 10 4 , A C = 1.25, B C = 1.0 × 10 3 and ε f = 0.5 × 10 −4 (X. Xia, Zhang, and Tang 2007). The shell thickness is determined as t = 140 mm, which is uniform over the shells and fixed in the optimization process. The convergence criterion in Algorithm 1, i.e. the tolerance value of the relative change of objective, is set to be = 0.1%.
To study the effect of material damage on the optimal result, shape optimization is repeated for the same shells except for purely elastic material. The optimization is performed using the Fortran code of Algorithm 1. In the case of purely elastic material, the damage parameter D is assumed to be zero and linear FE analysis is implemented. The computer specs are as follows: the operating system is Windows 7; the CPU is an Intel Core tm i7-6700 CPU@3.40 Hz (8CPUs), 3.4 GHz; the memory is 12,288 MB RAM.

Example 5.1 (A free-form shell with a continuous surface):
In this example, a free-form shell with a continuous surface is optimized. The geometry and boundary conditions as well as the mesh are given in Figure 2. The plane configuration of the shell is a sector. Two edges are simply supported, which are depicted as heavy lines in Figure 2, while the other two are free edges. Using the proposed algorithm, the shell is optimized for minimum strain energy under a uniformly distributed vertical load q. Considering that the material damage is related to the magnitude of the load, the optimization is repeated under q = 18.0 kN/m 2 and q = 8.0 kN/m 2 , respectively. NURBS functions are used to generate the free-form shape. The control points are located as shown in Figure 2(a) (represented by dots), with a total number of 20, which is determined by the tradeoff between accuracy and efficiency. Both degrees of the basis function in two directions are given as two. All control-point weights are determined as 0.5, which are constant in the optimization.
As shown in Algorithm 1, the optimization problem is solved by the gradient-based method. The initial Z-coordinates of control points are given randomly, by which the initial guess shape is generated, as shown in Figure 2(b). Since the Z-coordinates of control points are various, the initial shape shows a free-form feature. It is noted that the NURBS surface of the initial guess has good smoothness, which means the degree of two is acceptable. A triangular mesh is generated, as Figure 2 shows.
First, the optimization is performed under q = 18.0 kN/m 2 accounting for material damage. The Newton-Raphson repeat method is used to solve the nonlinear FE problem. The load increments are set to be the same and equal to 0.1q. In each load increment, the iterative step is set as 10. It is noted that, in the following examples, load increments and iterative steps are set to be the same as in Example 5.1. To ensure numerical accuracy, the optimization is repeated with the number of elements, i.e. 384, 600, 720 and 864, respectively. All the optimizations converge within 100 steps. Their respective evolution curves are illustrated in Figure 3, in which the vertical axis, denoted as C i /C 0 , indicates the ratio of the strain energy of the ith optimization step to that of the initial guess. It can be seen that, when the number of elements is increased from 384 to 600, the evolution curve changes noticeably; when the number is larger than 600, the evolution processes are very close. Therefore, the number of this model is determined as 600. The shape variation process with respect to the element number of 600 is also presented in Figure 3. It is observed that, through the optimization, the wavy surface of the shell becomes smooth and finally converges to a cylinder-like shape, which is supported on the two opposite edges.
The damage variation during the optimization is depicted in Figure 4. Note that this section selects the maximum damage value of the upper and lower surface to represent its distribution. It is the same in the case of stress. It is observed that, on the initial shell (Figure 4(a)), the damage occurs in a large area with a large damage parameter close to 1.0, while at the 5th optimization step, the damage area has been narrowed significantly (Figure 4(b)). After the 20th optimization step, the damage parameter over all the shell has decreased to zero (Figure 4(c)), meaning that material damage has been avoided.
To investigate the effect of concrete damage on the optimal shape, optimization is also conducted for the same model except with linear elastic material. The optimization process is similar to the damage-considered case except that it converges at the 60th step. For comparison, the optimal surfaces of damage-considered and elastic designs are illustrated in Figure 5 with Gaussian curvature distributions. It can be seen that the overall shapes of the two optimal designs are similar. However, compared with the elastic design, the curvature of the damage-considered design is larger in some regions. In addition, the computational costs are also compared. The damage-considered optimization costs 2.50 h, while the optimization with respect to linear elastic material costs 0.17 h. The consideration of damage increases the computational cost. However, the computational cost is still acceptable.
Moreover, the displacement and principle stress of the two optimal designs are evaluated using damage-considered nonlinear FE analysis under q = 18.0 kN/m 2 , as presented in Figures 6. It can be seen that the distribution patterns of displacement and stress of the two designs are similar, which is due to the fact that their overall shapes are similar. However, the maximum displacement and stress of the damage-considered design are smaller than those of the elastic design by 12.50% and 33.00%, respectively. That is to say, the damage-considered design has improved the mechanical performance compared with the elastic design.
Furthermore, to investigate the effect of initial guess shapes on optimization results, the damageconsidered optimization is re-performed with two other initial guess shapes. It is observed that the  effect of initial guess shape is relatively small in the case studied. The details can be found in the online Supplemental Data for this article.
The optimization and analysis is repeated under load value q = 8.0 kN/m 2 . The results show that the optimal shape and structural response of the damage-considered design are very close to those of the elastic design when q = 8.0 kN/m 2 . The reason is that, when the surface is continuous, without any sharp corner or opening, the material damage is not serious under small loads. However, when the load is large, such as q = 18.0 kN/m 2 , the effect of material damage needs to be considered as presented above.

Example 5.2 (A free-form shell with an opening):
A concrete free-form shell with an opening in the centre is considered in this example. Such openings are usually designed in civil structures for skylight or other non-load-bearing components. The model is described in Figure 7. Both the outer and inner edges are square, with respective side lengths of 36 and 12 m. All four outer edges are simply supported and the inner edges are free. The shell is prescribed to be symmetrical about central lines. Therefore, only a quarter of the shell is modelled, as shown in Figure 7(c).
Similar to Example 5.1, shape optimization is performed accounting for material damage under a uniformly distributed vertical load to obtain minimum strain energy. And the element number of the quarter model is determined as 1600 by a mesh sensitivity study. Considering that the opening tends to cause stress concentrations, a relatively small load is selected: q = 8.0 kN/m 2 . The locations of control points are shown in Figure 7(c), with a total number of 15. The weight factors and basis degree are taken to be the same as those in Example 5.1.
The optimization starts with a given initial guess shape (Figure 7(b)). The initial Z-coordinates of control points are given randomly, leading to a free-form shape. The evolution of objective value is plotted in Figure 8, which shows good convergence. At the final step, i.e. the 88th step, the strain energy decreases to 6.83% of the initial strain energy. The initial shape, the 10th step shape and the optimal (88th step) shape are also presented in Figure 8 to show the shape evolution. It can be observed that the shape changes significantly during the optimization: the surface around the square opening is warped and waved at the initial stage and then gradually becomes lower and smooth; by contrast, the central part of the quarter model changes from a flat surface to a surface with large positive Gaussian curvature.
The damage parameter contours are shown in Figure 9. Large material damage is observed in the initial shape, which is particularly large around the opening and near the outer edges. This can be explained by the stress concentrations. Through the first five optimization steps, material damage is reduced dramatically and only appears in a small region near the opening. After the 10th optimization step, material damage disappears over all the shell.
Again, the elastic optimization is conducted for comparison. In this example, the computational costs of the damage-considered and purely elastic cases are 3.20 and 0.31 h, respectively. Figures 10(a) and 10(b) show the damage-considered and purely elastic optimal designs, respectively. Compared with Example 5.1, the curvature of the damage-considered design is much larger than that of the elastic design. It is demonstrated that the effect of material damage on shape is more significant because the opening leads to stress concentrations and consequent large material damage. Figure 11 shows the displacement and principle stress of both the damage-considered and elastic optimal designs when evaluated using damage-considered nonlinear FE analysis. It can be seen that their displacement and stress patterns are similar and not influenced significantly by the shape difference. The reason could be that the boundary conditions play a more important role on the distribution of displacement and stress than the geometrical shape does. However, the shape difference leads to a diverse magnitude of structural responses. The maximum displacement and principle stress of the elastic design are larger than those of the damage-considered design by 18.22% and 20.12%, respectively.
It can be concluded from this example that the opening is a source of stress concentration and therefore of material damage, which results in the significant effect of material damage on the optimal design.

Example 5.3 (A free-form shell with weak boundary constraints):
The shell model shown in Figure 12 is optimized in this example. Differing from Examples 5.1 and 5.2, the boundary constraint of this model is relatively weak: only two edges are simply supported, which are denoted by heavy  lines in Figure 12, while the other four edges are free. The shell is subject to uniformly distributed vertical load q = 8.0 kN/m 2 . Following the previous two examples, optimization and comparative analysis are performed.
The control points of the NURBS surface are located as shown in Figure 12(a), with the total number being 21. The weight factors and basis degree are taken to be the same as those in Examples 5.1 and 5.2. The initial guess shape is given as shown in Figure 12(b), which is generated by randomly given Zcoordinates of control points. The element number of the quarter model is determined to be 600 by a mesh sensitivity study. The evolutions of objective and shape are presented in Figure 13, which show a slower convergence than the previous two examples. In the first 100 steps, the strain energy decreases sharply. After the 100th step, the decrease of strain energy becomes gradual and finally converges at the 440th step. The relatively slow convergence leads to a large computational cost, which is 11.00 h. The strain energy decreases to 2.14% of the initial strain energy. During the optimization process, the shape of the shell changes from a simple surface to a complex waved surface. This significantly differs from Examples 5.1 and 5.2. The reason is that, in this example, the boundary constraint is weak (only two of six edges are supported), therefore the shell needs more waves to obtain large stiffness. Similar to Examples 5.1 and 5.2, material damage occurs in most regions of the initial shell, while with the evolution of the shape, the material damage decreases to zero over all the shell. For brevity, the figures for damage variation are not given.
Elastic optimization is also conducted, which costs 0.97 h. Figure 14 shows both the damageconsidered optimal design and the purely elastic optimal design. It can be seen that the shape difference between the two designs is much more significant than those in Examples 5.1 and 5.2. Consequently, their displacement and stress responses are also significantly different when evaluated using damage-considered nonlinear FE analysis. The differences in displacement and principle stress are 14.89% and 22.40%, respectively. Through this example, it is concluded that, when the boundary constraint is weak, the shape of the shells plays a more significant role, and therefore considering material damage or not leads to diverse optimization results. The damage-considered design is superior in mechanical performance.

Conclusions
A damage-considered shape optimization algorithm was developed for concrete free-form shells. In the algorithm, the Mazars model was combined with an FE analysis code to model the initiation and propagation of damage. NURBS were employed to represent the free-form shapes. The locations of NURBS control points were updated by the gradient-based method to minimize the structural strain energy. The algorithm was self-coded and several numerical examples were presented. Via the study, the following conclusions were drawn.
(1) Through the numerical examples, the proposed damage-considered optimization algorithm was proved to be accurate and efficient. Differently from a purely elastic shape optimization  algorithm, material damage was modelled and decreased significantly during the optimization process.
(2) The consideration of material damage affects the resulting optimal shapes. When the load is large or openings exist, the effect of material damage is large; when the boundary constraint is strong, material damage mainly influences the curvatures of the shell; however, when the boundary constraint is relatively weak, not only the curvatures but also the overall optimal shape is changed. It is demonstrated that it is necessary to consider material damage in shape optimization, especially in the aforementioned cases. (3) As a consequence of shape diversity, the structural performance of the damage-considered design is vastly superior to that of the purely elastic design. That is to say, the consideration of material damage leads to more robust designs of concrete shells. This further highlights the importance of accounting for material damage in the shape optimization of concrete free-form shells.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Data availability statement
The data that support the findings of this study are available from the corresponding author, Ye Qiu, upon reasonable request.