Stochastic Natural Vibration Analyses of Free-Form Shells

This work investigates the natural vibration characteristics of free-form shells when considering the influence of uncertainties, including initial geometric imperfection, shell thickness deviation, and elastic modulus deviation. Herein, free-form shell models are generated while using a self-coded optimization algorithm. The Latin hypercube sampling (LHS) method is used to draw the samplings of uncertainties with respect to their stochastic probability models. ANSYS finite element (FE) software is adopted to analyze the natural vibration characteristics and compute the natural frequencies. The mean values, standard deviations, and cumulative distributions functions (CDFs) of the first three natural frequencies are obtained. The partial correlation coefficient is adopted to rank the significances of uncertainty factors. The study reveals that, for the free-form shells that were investigated in this study, the natural frequencies is a random quantity with a normal distribution; elastic modulus deviation imposes the greatest effect on natural frequencies; shell thickness ranks the second; geometrical imperfection ranks the last, with a much lower weight than the other two factors, which illustrates that the shape of the studied free-form shells is robust in term of natural vibration characteristics; when the supported edges are fixed during the shape optimization, the stochastic characteristics do not significantly change during the shape optimization process.


Introduction
The term free-form shells refers to shells whose geometric shapes cannot be represented by certain mathematical formulas or their combination [1]. When compared with traditional shells, such as cylindrical shells and spherical shells, the geometrical shape of free-from shells is flexible to satisfy certain structural, architectural, and other requirements, which advances the wide application of free-form shells in civil engineering, aerospace engineering, and so on.
The mechanical properties of free-form shells are strongly dependent on their geometric shapes. Therefore, shape optimization became an indispensable approach for designing a free-form shell. Bletzinger et al. [2] optimized free-form shells for structural stiffness under given loads while using numerical methods, which were emerged with physical experiments. Ohmori and Hamada [3] combined Non-Uniform Rational B-Spline (NURBS) with genetic algorithm (GA), as well as the gradient method to perform shape optimization of free-from shells, in which the objective functions are set as minimum strain energy and minimum geometric deviation from the prescribed shape. Cui and Yan [4] developed a height adjusting method to obtain optimal shape of free-form shells, with respect to minimum strain energy. Feng and Ge [1] described a conjugate gradient optimization method for the shape design of cable-braced free-form grid shells, with strain energy as the optimization objective. Wang and Wu [5] conducted a study on local optimal solutions and a modified optimization method, which was applied on the shape optimization of cable-stiffened latticed free-form shells, with the Appl. Sci. 2019, 9,  objective of the minimization of strain energy. It can be seen that, in the past decades, numerical optimization methods have been developed for free-form shells, many of which made efforts to obtain a geometrical shape with minimum strain energy, i.e., largest structural stiffness. There is a relative dearth of research on mechanical properties of free-form shells despite the fast development of shape optimization of free-form shells. Only minimal literature is concerned with that. For example, Li [6] studied the buckling properties of free-form shells; Cui et al. [7,8] carried out experimental and numerical studies on the static behaviors of the free-form concrete shells. However, the existing work mainly studies static characteristics and it is limited to the deterministic analysis without any consideration of uncertainty factors, which are unavoidable in practical engineering, such as geometrical imperfection and deviations of material properties. The previous stochastic studies on traditional shells have demonstrated the significant effect of uncertain factors on mechanical properties of structures [9][10][11][12][13]. The construction process of free-form shells is more complex when compared with traditional shells, which tends to lead more uncertainty factors. This further motivates this study to consider the effects of uncertainty factors on mechanical properties of free-form shells, particularly on the dynamic characteristics.
In this study, the stochastic natural vibration analyses of free-form shells is carried out while considering three uncertainty factors, including geometric imperfection, shell thickness deviation, and elastic modulus deviation. Free-form shell models are generated by a self-coded numerical optimization algorithm. Latin hypercube sampling (LHS) is used for the random sampling of three uncertainty factors based on their stochastic probability modeling. The global sensitivity study is conducted to rank the importance of these error sources.
The outline of this paper is as follows: in Section 2, the shape generation method of free-form shells is presented and programmed; Section 3 describes the stochastic and sensitivity analysis method; in Section 4, stochastic and sensitivity analyses are carried out on free-form shell, which are designed by the shape generation algorithm presented in Section 2. The results are reported along with discussions; finally, Section 5 provides closing remarks.

Shape Parametrization
The geometric shape of free-form shells cannot be represented by certain mathematical formulas or their combination, different from regular shells, such as spherical shells and cylindrical shells. Therefore, shape parameterization is needed to describe their arbitrary shapes with desired smoothness. NURBS is employed here to represent the shape of structures, which is developed from B-spline and have been widely applied. This section gives a brief introduction on NURBS. For more details, one can refer to [14].
A NURBS surface, as shown in Figure 1, is defined as where S(u, v) is the value of a point on the NURBS surface, P i,j represents the coordinate positions of a set of control points, which form a m × n bi-directional control-point grid, and ω i,j is their respective weights; N i,k (u) and N j,l (v) are B-spline basis functions, defined on the knot vector ; u i and v i are the real numbers representing the coordinates in the parametric space [0, 1]; k and l are degrees.
The aforementioned equations have revealed that the geometry of a NURBS surface is affected by its control points and their weights. However, the work that is presented in [6] shows that the control points play a dominant role in changing shapes rather than the weights. Although setting both locations and weights of control points as design variables could give more flexibility to the shape and enlarge the design space, it would significantly increase the computational cost. For this reason, we only consider the coordinates of control points as the optimization variables while keeping the weights fixed. Particularly, the optimal shape is approached by updating the Z-coordinate of the control points, while setting their X, Y-coordinates as the fixed parameters.
N v is defined in the similar way. The aforementioned equations have revealed that the geometry of a NURBS surface is affected by its control points and their weights. However, the work that is presented in [6] shows that the control points play a dominant role in changing shapes rather than the weights. Although setting both locations and weights of control points as design variables could give more flexibility to the shape and enlarge the design space, it would significantly increase the computational cost. For this reason, we only consider the coordinates of control points as the optimization variables while keeping the weights fixed. Particularly, the optimal shape is approached by updating the Z-coordinate of the control points, while setting their X, Y-coordinates as the fixed parameters.

Optimization Algorithm
In this study, the free-form shape is generated through an optimization approach that minimizes the strain energy. A gradient-based method is employed to solve the problem, in which a gradient descent, or negative gradient direction of the objective function is selected as the search direction of each iteration step, and it gradually approaches the minimum function value. In gradient-based optimization, it is indispensable to solve the derivatives of the objective and/or constraint functions with respect to the design variables, which could be computed while using analytical, semianalytical, or numerical methods, depending on the problem itself and the available resources [15]. Here, a numerical differentiation method is selected to compute derivatives, since we have relatively small number of design variables.
Combining the finite element (FE) forward analysis with the shell shape being parameterized by NURBS, the gradient-based optimization procedure is obtained. The code is programmed while using the FORTRAN language. The FE modeling is based on the Kirchhoff-Love plate theory and the plane triangular shell element with 3 nodes is employed. For better understanding, the algorithm is

Optimization Algorithm
In this study, the free-form shape is generated through an optimization approach that minimizes the strain energy. A gradient-based method is employed to solve the problem, in which a gradient descent, or negative gradient direction of the objective function is selected as the search direction of each iteration step, and it gradually approaches the minimum function value. In gradient-based optimization, it is indispensable to solve the derivatives of the objective and/or constraint functions with respect to the design variables, which could be computed while using analytical, semi-analytical, or numerical methods, depending on the problem itself and the available resources [15]. Here, a numerical differentiation method is selected to compute derivatives, since we have relatively small number of design variables.
Combining the finite element (FE) forward analysis with the shell shape being parameterized by NURBS, the gradient-based optimization procedure is obtained. The code is programmed while using the FORTRAN language. The FE modeling is based on the Kirchhoff-Love plate theory and the plane triangular shell element with 3 nodes is employed. For better understanding, the algorithm is described in Algorithm 1. Given material properties, shell thickness, load, boundary conditions 2.
Set up the location of control points 3.
Assign an initial guess for the design variables, i.e., Z-coordinates of control points 4.
Select the convergence criteria, i.e., the relative change of the objective function is smaller than a given tolerance value 5.
Loop until convergence (1) Implement NURBS to obtain the geometry of the shape (2) Generate FE mesh with surface parameterized by NURBS (3) Solve the FE forward problem and obtain the objective function value (4) Compute derivatives of objective with respect to Z i using a numerical method (5) Compute the step size of Z i by the Golden section method (6) Update the design variables Z i 6.

Stochastic and Sensitivity Analysis Method
The natural vibration analyses are coupled with LHS to investigate the stochastic characteristics of natural frequencies of free-form shells and to explore the sensitivity of the uncertainty factors.
LHS [16][17][18] is a widely-used sampling method, which allows for the extraction of a large amount of stochastic and sensitivity information with a relatively small sample size. It can ensure that the entire distribution of each input random variable is covered and good at estimating the mean values of output variables [19]. LHS consists of two main steps, sampling and permutation. Sampling is implemented to produce representative samples to describe the distribution of each input random variable, while permutation aims at reducing the correlations between the samples of different random variables.
Herein, LHS is used for sampling. e 1 , e 2 , · · · , e n e denote the input uncertainty factors. Their respective probability cumulative function is expressed as F i = F i (e i ). f represents the output variable, i.e., the natural frequency in the current study. For a sample size n s , the range of F i is divided into n s equal interval zones, which are non-overlapping. One value is randomly selected for each intervals and the sampling. A row of sampled values e i1 , e i2 , · · · , e in s is obtained. A n e × n s sampling matrix can be generated after all the input random variables are sampled. Subsequently, random permutation is performed, in which the sampled value of every input variable is randomly chosen out of the sample values from the corresponding row in sampling matrix without replacement.
To check the sample size, the convergence estimation formulas are given, as follows where n s is the number of samplings; ch is added number of samplings; f (n s − ch), f (n s ) is mean value of output results with n s − ch samplings and n s samplings, respectively; σ(n s − ch) and σ(n s ) is standard deviation of output results with n s − ch samplings and n s samplings, respectively; referring the existing work [19][20][21], the tolerances ξ 1 and ξ 2 are determined as values that are less than or equal to 0.01. Based on LHS, the natural frequencies of free-form shells, which involve inherent randomness, is investigated by a partial correlation-based method [17][18][19]22,23]. The partial correlation coefficient was proposed by Iman and Helton [22], which is used to measure the degree of the linear correlation between the input variable and the model output after making an adjustment to remove the linear effect of all the remaining variables. The derivation of the partial correlation coefficient can be referred to [23]. Herein, the formulation of the partial correlation coefficients between input uncertainty factors e i and output variables f (i.e., natural frequencies in the current study) is given, as follows [22,23] where c i f , c ii , and c f f are the elements of matrix C, written as follows c 11 c 12 · · · c 1n e c 1 f · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · c n e 1 c n e 2 · · · c n e n e c n e f c f 1 c f 2 · · · c f n e c f f r 11 r 12 · · · r 1n e r 1 f · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · r n e 1 r n e 2 · · · r n e n e r n e f r f 1 r f 2 · · · r f n e r f f where r ij is linear correlation coefficient between input variables e i and e j ; r i f is the linear correlation coefficient between input variable e i and output variable f , see Equation (6) Partial correlation coefficients S i , are dimensionless parameters with a range between −1.0 and +1.0 and reveal the degree of association between the input and output variables. Herein, partial correlation coefficients are applied to rank the importance of each input variables.
In this study, three uncertainty sources e 1 , e 2 , e 3 are considered, which are initial geometric imperfection of shells, shell thickness deviation, and elastic modulus deviation. The output variable f , i.e., natural frequency, is computed by software ANSYS while using the Block lanczos method.

Results and Discussion
As presented above, three main stochastic factors that may occur on shells are explored, including geometric imperfection of shells, shell thickness deviation, and elastic modulus deviation. The distribution of initial geometric imperfection throughout the shell is determined by the consistent mode imperfection method [24], while the shell thickness deviation and elastic modulus deviation are assumed to be uniformly distributed over the shell. Since each of uncertainty error sources is determined by several independent factors, they are assumed to follow normal distribution [25]. Referring to the existing work and code [26,27], their respective mean values and standard deviations are determined, as shown in Table 1, where initial geometric imperfection is described by the maximum height deviation of the shell; the thickness and elastic modulus deviation is represented by the relative error between actual and design values, with "+" and "−" meaning the actual value larger and smaller than the design value, respectively. The presented stochastic analysis method is implemented on two widely used free-form shells. One is a free-form shell with negative Gaussian curvature and the other one is a combination of two positive-Gaussian domes. Both of the models are generated by the optimization approach, which is presented in Section 2. The natural vibration analysis is implement by ANSYS finite element software.

Model Generation
The shell model in this example is as shown in Figure 2, which cover a rectangle plane section with a dimension of 8 × 12m. Both long-span edges are simply supported, while the two short-span edges are free. The material properties and the design shell thickness are as shown in Table 2. The material is assumed to be linear elastic.
The shape is generated by the proposed optimization approach. NURBS functions are adopted to represent the geometric shape, with control points that are shown in Figure 2a, which forms a  53 grid. As aforementioned, the weight factors have small effects on the NURBS surface [6]. Therefore, the weight factors of all control points are chosen randomly from the range 0~1.0 as 0.5 and not changed during the optimization process. Both the degrees of the basis functions in the two directions are selected to be 2. The optimization aims at minimizing the strain energy under a given uniformly distributed vertical load  2 3.5kN/m q . The strain energy is solved by a self-coded program, in which triangular shell element is used with 768 elements and 425 nodes.
The gradient-based algorithm is used to solve the optimization problem. A cylindrical shell with the height of 4 m is given as an initial guess, as shown in Figure 2b. In this example, the control point on the long-span edges are fixed meaning the location and shape of the two edges are not changing during the optimization process. Consequently, the total number of optimization variables is 5.    The shape is generated by the proposed optimization approach. NURBS functions are adopted to represent the geometric shape, with control points that are shown in Figure 2a, which forms a 5 × 3 grid. As aforementioned, the weight factors have small effects on the NURBS surface [6]. Therefore, the weight factors of all control points are chosen randomly from the range 0~1.0 as 0.5 and not changed during the optimization process. Both the degrees of the basis functions in the two directions are selected to be 2. The optimization aims at minimizing the strain energy under a given uniformly distributed vertical load q = 3.5kN/m 2 . The strain energy is solved by a self-coded program, in which triangular shell element is used with 768 elements and 425 nodes.
The gradient-based algorithm is used to solve the optimization problem. A cylindrical shell with the height of 4 m is given as an initial guess, as shown in Figure 2b. In this example, the control point on the long-span edges are fixed meaning the location and shape of the two edges are not changing during the optimization process. Consequently, the total number of optimization variables is 5.
The variation of the objective function during the optimization process is depicted in Figure 3 to show the convergence process, in which the vertical coordinate indicates the ratio of the strain energy of the ith optimization step to that of the initial guess. It can be seen that the strain energy decreases significantly during the first 50 optimization steps. After the 50th step, the strain energy decreases slowly and then converges at the 400th step. Through the optimization, the strain energy decreases to 40% of that of the original design.
In Figure 4, the initial, 20th, 200th, and 400th (optimal) shapes are selected to show the shape variation during optimization process. It is observed that, via the first 20 optimization steps, the shape changes from a cylinder to a saddle to earn larger stiffness. Since the 20th step, the shape remains The variation of the objective function during the optimization process is depicted in Figure 3 to show the convergence process, in which the vertical coordinate indicates the ratio of the strain energy of the ith optimization step to that of the initial guess. It can be seen that the strain energy decreases significantly during the first 50 optimization steps. After the 50th step, the strain energy decreases slowly and then converges at the 400th step. Through the optimization, the strain energy decreases to 40% of that of the original design.  In Figure 4, the initial, 20th, 200th, and 400th (optimal) shapes are selected to show the shape variation during optimization process. It is observed that, via the first 20 optimization steps, the shape changes from a cylinder to a saddle to earn larger stiffness. Since the 20th step, the shape remains saddle-shaped, and the surface curvature increases overall. The shape change is too slight to be observed after the 200th optimization step.  Prior to the stochastic analysis, the deterministic analysis of free vibration is implemented on shells with respect to several selected optimization steps, and the first three natural frequencies are plotted in Figure 5. It can be seen that, all the first three natural frequencies significantly change during the optimization process and increase by 158.0% (1st), 61.1% (2nd), and 111.5% (3rd), respectively.   slowly and then converges at the 400th step. Through the optimization, the strain energy decreases to 40% of that of the original design.  In Figure 4, the initial, 20th, 200th, and 400th (optimal) shapes are selected to show the shape variation during optimization process. It is observed that, via the first 20 optimization steps, the shape changes from a cylinder to a saddle to earn larger stiffness. Since the 20th step, the shape remains saddle-shaped, and the surface curvature increases overall. The shape change is too slight to be observed after the 200th optimization step. Prior to the stochastic analysis, the deterministic analysis of free vibration is implemented on shells with respect to several selected optimization steps, and the first three natural frequencies are plotted in Figure 5. It can be seen that, all the first three natural frequencies significantly change during the optimization process and increase by 158.0% (1st), 61.1% (2nd), and 111.5% (3rd), respectively.   Prior to the stochastic analysis, the deterministic analysis of free vibration is implemented on shells with respect to several selected optimization steps, and the first three natural frequencies are plotted in Figure 5. It can be seen that, all the first three natural frequencies significantly change during the optimization process and increase by 158.0% (1st), 61.1% (2nd), and 111.5% (3rd), respectively. slowly and then converges at the 400th step. Through the optimization, the strain energy decreases to 40% of that of the original design.  In Figure 4, the initial, 20th, 200th, and 400th (optimal) shapes are selected to show the shape variation during optimization process. It is observed that, via the first 20 optimization steps, the shape changes from a cylinder to a saddle to earn larger stiffness. Since the 20th step, the shape remains saddle-shaped, and the surface curvature increases overall. The shape change is too slight to be observed after the 200th optimization step. Prior to the stochastic analysis, the deterministic analysis of free vibration is implemented on shells with respect to several selected optimization steps, and the first three natural frequencies are plotted in Figure 5. It can be seen that, all the first three natural frequencies significantly change during the optimization process and increase by 158.0% (1st), 61.1% (2nd), and 111.5% (3rd), respectively.

Stochastic and Sensitivity Analysis
Three representative shapes generated in the optimization, i.e., the initial shape, the 20th-step shape, the 400th-step optimal shape, see Figure 4, are selected for the stochastic analysis of natural frequencies.
LHS is adopted to generate samples. To determine the sample size, sampling is repeated on the optimal shape with size of 350, 400, 450, 500, 550, and 600, respectively. The mean values and standard deviations of first three frequencies are solved and checked by Equation (3), with ξ 1 , ξ 2 set up as 0.01 and 0.001, respectively. It is found that, when the sampling size is 400, Equation (3) is satisfied for the first three frequency of all the selected shapes. For brevity, only the check result of the first natural frequency of the optimal shape is given in Figure 6. Herein, the size of samples is determined as 400 and Figure 7 depicts the samples.
shape, the 400th-step optimal shape, see Figure 4, are selected for the stochastic analysis of natural frequencies.
LHS is adopted to generate samples. To determine the sample size, sampling is repeated on the optimal shape with size of 350, 400, 450, 500, 550, and 600, respectively. The mean values and standard deviations of first three frequencies are solved and checked by Equation (3), with  1 ,  2 set up as 0.01 and 0.001, respectively. It is found that, when the sampling size is 400, Equation (3) is satisfied for the first three frequency of all the selected shapes. For brevity, only the check result of the first natural frequency of the optimal shape is given in Figure 6. Herein, the size of samples is determined as 400 and Figure 7 depicts the samples. The stochastic analysis is implemented on the initial shell model, the 20th-step shell model and the 400th-step (optimal) shell model, respectively. The natural frequencies of all samples are solved by FEM using ANSYS. We first focus on the study of first natural frequencies, which are shown in Figure 8. Mean values and standard deviations of the first natural frequency are obtained as 4.56Hz and 0.15 Hz (initial), 6.13 Hz and 0.20 Hz (20th step), and 11.76 Hz and 0.36 Hz (optimal shape), respectively. LHS is adopted to generate samples. To determine the sample size, sampling is repeated on the optimal shape with size of 350, 400, 450, 500, 550, and 600, respectively. The mean values and standard deviations of first three frequencies are solved and checked by Equation (3), with  1 ,  2 set up as 0.01 and 0.001, respectively. It is found that, when the sampling size is 400, Equation (3) is satisfied for the first three frequency of all the selected shapes. For brevity, only the check result of the first natural frequency of the optimal shape is given in Figure 6. Herein, the size of samples is determined as 400 and Figure 7 depicts the samples. The stochastic analysis is implemented on the initial shell model, the 20th-step shell model and the 400th-step (optimal) shell model, respectively. The natural frequencies of all samples are solved by FEM using ANSYS. We first focus on the study of first natural frequencies, which are shown in Figure 8. Mean values and standard deviations of the first natural frequency are obtained as 4.56Hz and 0.15 Hz (initial), 6.13 Hz and 0.20 Hz (20th step), and 11.76 Hz and 0.36 Hz (optimal shape), respectively. The stochastic analysis is implemented on the initial shell model, the 20th-step shell model and the 400th-step (optimal) shell model, respectively. The natural frequencies of all samples are solved by FEM using ANSYS. We first focus on the study of first natural frequencies, which are shown in Figure 8. Mean values and standard deviations of the first natural frequency are obtained as 4.56Hz and 0.15 Hz (initial), 6.13 Hz and 0.20 Hz (20th step), and 11.76 Hz and 0.36 Hz (optimal shape), respectively. Initial guess shape 20th-step shape 400th-step (optimal )shape Via the stochastic analysis, it is found that the first natural frequency follows a normal distribution. The cumulative distributions of first natural frequency of the three models are depicted and fitted in Figure 9 for clear analysis and comparison. As presented in Figure 9a, it can be seen that the value of first natural frequency of the initial guess ranges from 4.05 to 5.00 Hz and the cumulative Via the stochastic analysis, it is found that the first natural frequency follows a normal distribution. The cumulative distributions of first natural frequency of the three models are depicted and fitted in Figure 9 for clear analysis and comparison. As presented in Figure 9a, it can be seen that the value of first natural frequency of the initial guess ranges from 4.05 to 5.00 Hz and the cumulative distributions function (CDF) is fitted as where x denotes the first natural frequency. The determination coefficient is 0.9883, which shows a good approximation.  (7) where x denotes the first natural frequency. The determination coefficient is 0.9883, which shows a good approximation. It is also observed that the CDFs of 20th-step shape and optimal shape is similar to Equation (7) and fitted as Equations (8) and (9)  Partial correlation coefficients with respect to the first natural frequency are calculated and subsequently normalized to rank the significance of each factor, as shown in Figure 10. It is illustrated that, for initial guess, shell thickness and elastic modulus deviations impose relatively significant effects with similar weights, while the initial geometric imperfection is the least significant factor. In the case of the 20th-step shell, the significance ranking is similar to the initial guess, while the case of the optimal (400th-step) shape shows a difference that the effect of elastic modulus is remarkably higher than that of shell thickness deviation. The reason is that the increase (or decrease) of shell thickness leads to the increase (or decrease) of both stiffness and mass, which affects the frequencies in the opposite ways; in contrast, the increase (or decrease) of elastic modulus only affects the stiffness. Therefore, in most cases, the impact of the elastic modulus is larger than shell thickness. The difference between their significance depends on the geometrical shape to some extent.
The stochastic and sensitivity analysis is repeated for the second and third natural frequency. Similar results are obtained. 100 Initial geometric imperfection It is also observed that the CDFs of 20th-step shape and optimal shape is similar to Equation (7) and fitted as Equations (8) and (9), respectively Partial correlation coefficients with respect to the first natural frequency are calculated and subsequently normalized to rank the significance of each factor, as shown in Figure 10. It is illustrated that, for initial guess, shell thickness and elastic modulus deviations impose relatively significant effects with similar weights, while the initial geometric imperfection is the least significant factor. In the case of the 20th-step shell, the significance ranking is similar to the initial guess, while the case of the optimal (400th-step) shape shows a difference that the effect of elastic modulus is remarkably higher than that of shell thickness deviation. The reason is that the increase (or decrease) of shell thickness leads to the increase (or decrease) of both stiffness and mass, which affects the frequencies in the opposite ways; in contrast, the increase (or decrease) of elastic modulus only affects the stiffness. Therefore, in most cases, the impact of the elastic modulus is larger than shell thickness. The difference between their significance depends on the geometrical shape to some extent. Partial correlation coefficients with respect to the first natural frequency are calculated and subsequently normalized to rank the significance of each factor, as shown in Figure 10. It is illustrated that, for initial guess, shell thickness and elastic modulus deviations impose relatively significant effects with similar weights, while the initial geometric imperfection is the least significant factor. In the case of the 20th-step shell, the significance ranking is similar to the initial guess, while the case of the optimal (400th-step) shape shows a difference that the effect of elastic modulus is remarkably higher than that of shell thickness deviation. The reason is that the increase (or decrease) of shell thickness leads to the increase (or decrease) of both stiffness and mass, which affects the frequencies in the opposite ways; in contrast, the increase (or decrease) of elastic modulus only affects the stiffness. Therefore, in most cases, the impact of the elastic modulus is larger than shell thickness. The difference between their significance depends on the geometrical shape to some extent.
The stochastic and sensitivity analysis is repeated for the second and third natural frequency. Similar results are obtained.  Figure 11 shows the shell of this example. The two long-span edges are simply supported and the other edges are free. The material properties are same to that in Table 2. The shape is generated by the proposed optimization approach. NURBS functions are adopted to represent the geometric shape. The control points are located, as shown in Figure 11a, which forms a  53 grid. Weight factors and the degrees of the basis functions are determined as similar to that in Example 1. The The stochastic and sensitivity analysis is repeated for the second and third natural frequency. Similar results are obtained. Figure 11 shows the shell of this example. The two long-span edges are simply supported and the other edges are free. The material properties are same to that in Table 2. The shape is generated by the proposed optimization approach. NURBS functions are adopted to represent the geometric shape. The control points are located, as shown in Figure 11a, which forms a 5 × 3 grid. Weight factors and the degrees of the basis functions are determined as similar to that in Example 1. The optimization aims at minimizing the strain energy under applied load, which is also determined as the uniformly distributed vertical load of q = 3.5 kN/m 2 . The triangular shell element is used with 1024 elements and 561 nodes. As presented in Section 2, the gradient method is used to solve the optimization problem. The initial shape of this model is set up as a free-form shape, which consists of two connected positive -Gaussian domes. Similar to Example 1, the control point on the long-span edges is fixed, meaning that the location and shape of these two edges are not changing during the optimization process. The total number of optimization variables is 5, which is same to Example 1.

Model Generation
(a) (b) Figure 11. Layout of control points, mesh and the initial geometry of Model 2: (a) Layout of control points and mesh; and, (b) Initial guess shape. Figure 12 shows the convergence process of Model 2, which is slower than Model 1, since the shape is relatively complicated. The optimization converges at the 1000th step with a 65% decrease of strain energy. Figure 13 shows the changes of geometric shape of Model 2 during the optimization process. It can be seen that, since the control points on the edges are fixed, the general shape remain Figure 11. Layout of control points, mesh and the initial geometry of Model 2: (a) Layout of control points and mesh; and, (b) Initial guess shape.
As presented in Section 2, the gradient method is used to solve the optimization problem. The initial shape of this model is set up as a free-form shape, which consists of two connected positive -Gaussian domes. Similar to Example 1, the control point on the long-span edges is fixed, meaning that the location and shape of these two edges are not changing during the optimization process. The total number of optimization variables is 5, which is same to Example 1. Figure 12 shows the convergence process of Model 2, which is slower than Model 1, since the shape is relatively complicated. The optimization converges at the 1000th step with a 65% decrease of strain energy. Figure 13 shows the changes of geometric shape of Model 2 during the optimization process. It can be seen that, since the control points on the edges are fixed, the general shape remain similar to the initial shape, which consists of two connected positive-Gaussian surfaces. However, the curvature increases overall.
points and mesh; and, (b) Initial guess shape. Figure 12 shows the convergence process of Model 2, which is slower than Model 1, since the shape is relatively complicated. The optimization converges at the 1000th step with a 65% decrease of strain energy. Figure 13 shows the changes of geometric shape of Model 2 during the optimization process. It can be seen that, since the control points on the edges are fixed, the general shape remain similar to the initial shape, which consists of two connected positive-Gaussian surfaces. However, the curvature increases overall.
Prior to the stochastic analysis, the deterministic analysis of free vibration is implemented on shell models with respect to several selected optimization steps, and Figure 14 depicts their first three natural frequencies. It is observed that, different from Model 1, both the first and third natural frequencies decrease through the optimization, while only the second natural frequency increases. It is also noted that the change of all the frequencies are slight. For example, the relative change of the first natural frequency is only 11.9%, which is significantly smaller than Model 1. This can be explained by the fact that, in this case, the stiffness is mainly increased by the increase of curvature, which results in the increase of mass. The natural frequencies are affected by both stiffness and mass, and they consequently show a slight change. The natural frequencies decrease when the increase of mass plays a more significant role than stiffness.

Stochastic and Sensitivity Analysis
The initial, 100th, and 1000th optimized shape of Model 2 in Figure 13 are selected to investigate the influence of uncertainty factors on free frequencies. The samples are generated by LHS. The Prior to the stochastic analysis, the deterministic analysis of free vibration is implemented on shell models with respect to several selected optimization steps, and Figure 14 depicts their first three natural frequencies. It is observed that, different from Model 1, both the first and third natural frequencies decrease through the optimization, while only the second natural frequency increases. It is also noted that the change of all the frequencies are slight. For example, the relative change of the first natural frequency is only 11.9%, which is significantly smaller than Model 1. This can be explained by the fact that, in this case, the stiffness is mainly increased by the increase of curvature, which results in the increase of mass. The natural frequencies are affected by both stiffness and mass, and they consequently show a slight change. The natural frequencies decrease when the increase of mass plays a more significant role than stiffness.

Stochastic and Sensitivity Analysis
The initial, 100th, and 1000th optimized shape of Model 2 in Figure 13 are selected to investigate the influence of uncertainty factors on free frequencies. The samples are generated by LHS. The sample size is also determined as 400, which has been checked for the current model by Equation (3), with  1 2 , as 0.01 and 0.001, respectively.
Firstly, the first natural frequency is studied, which is depicted in Figure 15. The mean values of the initial guess, 100th-step shape, 1000th-step shape is 11.91 Hz, 11.58 Hz, and 10.58 Hz, respectively. Their standard deviation is 0.28 Hz, 0.32 Hz, and 0.28 Hz, respectively. The stochastic analysis shows that the first frequency of this shell model always follows the normal distribution during the optimization process. Their CDFs are depicted in Figure 16 and fitted, which are similar to Model 1. Initial guess shape 100th-step shape 1000th-step (optimal )shape

Stochastic and Sensitivity Analysis
The initial, 100th, and 1000th optimized shape of Model 2 in Figure 13 are selected to investigate the influence of uncertainty factors on free frequencies. The samples are generated by LHS. The sample size is also determined as 400, which has been checked for the current model by Equation (3), with ξ 1 , ξ 2 as 0.01 and 0.001, respectively.
Firstly, the first natural frequency is studied, which is depicted in Figure 15. The mean values of the initial guess, 100th-step shape, 1000th-step shape is 11.91 Hz, 11.58 Hz, and 10.58 Hz, respectively. Their standard deviation is 0.28 Hz, 0.32 Hz, and 0.28 Hz, respectively. The stochastic analysis shows that the first frequency of this shell model always follows the normal distribution during the optimization process. Their CDFs are depicted in Figure 16 and fitted, which are similar to Model 1.

Stochastic and Sensitivity Analysis
The initial, 100th, and 1000th optimized shape of Model 2 in Figure 13 are selected to investigate the influence of uncertainty factors on free frequencies. The samples are generated by LHS. The sample size is also determined as 400, which has been checked for the current model by Equation (3), with  1 2 , as 0.01 and 0.001, respectively.
Firstly, the first natural frequency is studied, which is depicted in Figure 15. The mean values of the initial guess, 100th-step shape, 1000th-step shape is 11.91 Hz, 11.58 Hz, and 10.58 Hz, respectively. Their standard deviation is 0.28 Hz, 0.32 Hz, and 0.28 Hz, respectively. The stochastic analysis shows that the first frequency of this shell model always follows the normal distribution during the optimization process. Their CDFs are depicted in Figure 16 and fitted, which are similar to Model 1.   Through the sensitivity analysis, the ranking of the three factors with respect to the first frequency is obtained, as shown in Figure 17. It can be seen that the elastic modulus deviation is the most important uncertainty, with a weight of around 60%; the shell thickness deviation has relatively small influence with a high weight of around 30%; and, the initial geometric imperfection imposes negligible effects with a weight of less than 10%. Generally, the ranking is similar to Example 1, except for that the sensitivity value of elastic modulus deviation is much higher. This could be attributed to the change of the mass, as discussed above. Through the sensitivity analysis, the ranking of the three factors with respect to the first frequency is obtained, as shown in Figure 17. It can be seen that the elastic modulus deviation is the most important uncertainty, with a weight of around 60%; the shell thickness deviation has relatively small influence with a high weight of around 30%; and, the initial geometric imperfection imposes negligible effects with a weight of less than 10%. Generally, the ranking is similar to Example 1, except for that the sensitivity value of elastic modulus deviation is much higher. This could be attributed to the change of the mass, as discussed above.   Through the sensitivity analysis, the ranking of the three factors with respect to the first frequency is obtained, as shown in Figure 17. It can be seen that the elastic modulus deviation is the most important uncertainty, with a weight of around 60%; the shell thickness deviation has relatively small influence with a high weight of around 30%; and, the initial geometric imperfection imposes negligible effects with a weight of less than 10%. Generally, the ranking is similar to Example 1, except for that the sensitivity value of elastic modulus deviation is much higher. This could be attributed to the change of the mass, as discussed above.

Conclusions
In this work, an optimization algorithm generates free-form shells. The stochastic analysis is carried out to investigate the stochastic characteristics of natural frequencies of free-form shells, in which three uncertainty factors are involved, including initial geometric imperfection, shell thickness deviation, and elastic modulus deviation. Furthermore, the global sensitivity analysis method is employed to rank the effect of three uncertainty factors. The main conclusions are drawn, as follows: (i) for the free-form shells investigated in this study, the natural frequencies of free-form shells follow a normal distribution in term of the three uncertainty factors. The elastic modulus deviation imposes the greatest effect; shell thickness ranks the second; geometrical imperfection ranks the last with a much lower weight than the previous two factors. It is demonstrated that the free-form shells studied in this paper have good robustness on geometrical shape error; and, (ii) in the current study, the supported edges are fixed during the shape optimization, which lead to a limited shape change. Therefore, the stochastic characteristics do not significantly change during the optimization process. The case with variable edges needs further study.

Conflicts of Interest:
The authors declare no conflict of interest.