# Mixed Static and Dynamic Optimization of Four-Parameter Functionally Graded Completely Doubly Curved and Degenerate Shells and Panels Using GDQ Method

Authors: Francesco Tornabene; Alessandro Ceruti
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This study deals with a mixed static and dynamic optimization of four-parameter functionally graded material (FGM) doubly curved shells and panels. The two constituent functionally graded shell consists of ceramic and metal, and the volume fraction profile of each lamina varies through the thickness of the shell according to a generalized power-law distribution. The Generalized Differential Quadrature (GDQ) method is applied to determine the static and dynamic responses for various FGM shell and panel structures. The mechanical model is based on the so-called First-order Shear Deformation Theory (FSDT). Three different optimization schemes and methodologies are implemented. The Particle Swarm Optimization, Monte Carlo and Genetic Algorithm approaches have been applied to define the optimum volume fraction profile for optimizing the first natural frequency and the maximum static deflection of the considered shell structure. The optimization aim is in fact to reach the frequency and the static deflection targets defined by the designer of the structure: the complete four-dimensional search space is considered for the optimization process. The optimized material profile obtained with the three methodologies is presented as a result of the optimization problem solved for each shell or panel structure.


---

## Body

## 1. Introduction

Shell structures are widely used in many fields of engineering thanks to their optimum dynamic behavior, strength, and stability guaranteed by the curvature effect. The dynamic and static deflection of these structures, caused by different external forces, can have serious consequences for their strength and safety, like resonance. Therefore, an accurate static deflection and frequency determination are of paramount importance for the technical design of these structural elements. One of the aims of this work is to study the static and dynamic behavior of completely doubly curved shell structures. During the last sixty years, two-dimensional linear theories of thin shells and plates have been developed including important contributions [1–12]. The transverse shear deformation has been incorporated into shell theories by applying the theory of Reissner-Mindlin [13, 14], also named First-order Shear Deformation Theory (FSDT). Abandoning the assumption of the preservation of the normals to the shell middle surface after the deformation, a comprehensive analysis for elastic isotropic shells and plates was made by Kraus [7] and Gould [15, 16]. Indeed, the present work is just based on the FSDT. In order to include the effect of the initial curvature in the evaluation of the stress resultants a generalization of the classical Reissner-Mindlin theory (CRMT) has been proposed in the literature by Kraus [7], by Leissa and Chang [17], by Qatu [18, 19] and by Toorani and Lakis [20]. As a consequence of these contributions the stress resultants directly depend on the geometry of the structure in terms of the curvature coefficients, and the hypothesis of the symmetry of the in-plane shearing force resultants and the torsional couples declines. A further improvement of the previous theories of shells has been proposed by Toorani and Lakis [21]. In the present work the kinematical model is generalized in order to include the effect of the curvatures from the beginning of the shell formulation. In this way, the strain relationships have been changed, and the equilibrium equations in terms of displacements have to be modified. The General First-order Shear Deformation Theory (GFSDT) is herein considered. It is worth noting that no results are available in the literature about this general theory for doubly curved shells. Thus, the motivation of the present work is based on the lack of results about completely doubly curved shells and panels.Due to the significant developments that have taken place in composite materials [22], the laminated composite doubly curved shells and panels are considered in this work. As for the static and dynamic analysis of shells, several studies have been presented earlier [23]. Furthermore, some complicated effects in shell structures have been considered in the recent years [24, 25]. Referring to the formulation of the static and dynamic equilibrium in terms of midsurface displacements and rotations, in this paper the system of second-order linear partial differential equations is solved. The static and dynamic solutions are obtained by using the numerical technique named Generalized Differential Quadrature (GDQ) method. The mathematical fundamentals and recent developments of the GDQ method as well as its major applications in engineering are discussed in detail in the book by Shu [26]. In the GDQ method the governing differential equations of equilibrium are directly transformed in one step to obtain the final algebraic form. The interest in researches dealing with this procedure is increasing due to its great simplicity and versatility. As shown in the literature [23, 27], GDQ technique is a global method, which can obtain very accurate numerical results by using a considerably small number of grid points. Therefore, this simple direct procedure has been applied in a large number of cases [28–77] to circumvent the difficulties of programming complex algorithms for the computer, as well as to avoid the excessive use of storage and computing time.In the last decades, the increased use of functionally graded materials (FGMs) [39, 49–55, 57–62, 64, 65, 70, 75, 78–91] in engineering structures calls for improved analysis and tailored design tools. Thus, in the present paper, functionally graded shells are considered. Typically, FGMs consist of a mixture of ceramic and metal, or a combination of different materials. In this study, ceramic-metal graded shells and panels with two different volume fraction power-law variations of the constituents in the thickness direction are considered. Two different four-parameter power-law distributions, proposed by Tornabene [49], are used for the ceramic volume fraction. Various material profiles through the functionally graded lamina thickness are chosen by varying the four parameters of the power-law distributions.The need for optimization lies in the mathematical formulation of FGMs, based on the four coefficients [49]: it is in fact quite difficult to handle these materials from the perspective of a designer. Small changes in parameters can lead in fact to strong changes in the distribution of base materials mixing in the thickness. Moreover, the domain of the four parameters is not continuous since some combinations of the parameters lead to distribution in which the mix of ceramics and metal resulting from the mathematical formula is negative or a complex number. From an operative perspective, a typical structural design scenario can include the need for minimizing the weight, the request of avoiding a resonance frequency, the control of the maximum displacement, or other constraints: in this case, there is no way to analytically relate the four parameter values to these important design requirements. As an answer, several authors proposed to apply optimization methods to achieve a feasible solution. For instance, Yas et al. [92] propose to minimize the weight of FGMs by using Imperialist Competitive Algorithm and Artificial Neural Networks. A very similar approach has been followed by Jam et al. [93] to optimize FGM conical shells. According to these examples, heuristic or semi heuristic methods are suitable in such a case since they consider the function to optimize as a “black box”, in which only inputs and outputs are considered. A lot of heuristic (or semi heuristic) optimization techniques have been proposed in the literature. A comprehensive but not complete list includes Tabu Search [94], Simulated Annealing [95, 96], Ant Colonies Optimization [97, 98], Genetic Algorithms [99–109], Differential Evolution [110], Particle Swarm Algorithm [111–116], Immune Systems [117, 118], Gravity Optimization [119], Imperialist Competitive Algorithm [120], and Intelligent Water Drop [121, 122]. In this paper, Genetic Algorithm (GA) and Particle Swarm Optimization (PSO) will be used since they are tested, and widespread methods and other techniques could have been applied. Also an optimization with Monte Carlo (MC) technique [115, 123–127] will be presented to check the good performance of GA and PSO and to evaluate how much these methods are efficient and reliable with respect to a completely random approach. From the point of view of the software implementation, the PSO method is more simple than the GA, but it requires more attention to boundary conditions. MC is considered as a support to the GA and PSO techniques: in this paper, it is used to check the functionality of GA and PSO and to confirm that heuristic methods are to be preferred to random approaches. A well-designed algorithm implementing optimization techniques like GA or PSO can be useful to obtain a better solution than MC, requiring a shorter time to solve the optimization problem. MC is in fact a random technique, in which large number of tries is necessary to obtain a good result: moreover, the search is blind and without memory of the past, so that the previous tests of the optimization domain are not exploited by the following computations. The above considerations suggested the application of heuristic methods to the optimization of FGMs; a more detailed description of the implementation of these methods will be provided in the next sections of the present work.The structure of this article is as follows: after this introduction, the second section describes the mathematical framework necessary to model and study FGM doubly curved shell structures; the third section shows the GDQ numerical implementation for the static and dynamic analyses; the fourth section describes the implementation of GA, PSO, and MC methods with respect to peculiarities of the case study, while the fifth one presents the results of the FGM optimization procedure applied to a set of six structures. A final section provides some conclusive comments regarding the problem description, the approach followed to solve it, and results obtained. The main novelty of the papers lies in the proposal of a design methodology which can be used by the designer involved in FGM applications to achieve a short time particular performances or characteristics for the FGM doubly curved shell structures.


## 2. Geometry Description and Shell Fundamental Systems of Equations

A 2D Equivalent Single Layer (ESL) model is proposed to study generic doubly curved shells and panels. The position of an arbitrary point within the shell medium is defined by coordinatesα1(α10≤α1≤α11), α2(α20≤α2≤α21) upon the middle surface or reference surface r(α1,α2), and ζ directed along the outward normal n(α1,α2) and measured from the reference surface (-h/2≤ζ≤h/2), where h(α1,α2) is the total thickness of the shell (Figure 1). The differential geometry [7, 11, 12, 23, 74] is used to describe the shell structure, starting with the position vector written in the global reference system:
(1)R(α1,α2,ζ)=r(α1,α2)+h(α1,α2)2zn(α1,α2),
where z=2ζ/h(α1,α2) and z∈[-1,1]. The basic configuration of the problem considered is a laminated composite doubly curved shell, as shown in Figure 1. For a laminated composite shell made of l laminae or plies, the total thickness h is defined as
(2)h=∑k=1lhk
in which hk=ζk+1-ζk is the thickness of the kth lamina or ply. In this study, doubly curved shells and degenerate shells such as plates are considered. It should be noted from (1) that the location of each point of the 3D shell is a function of the location of the point on the reference surface r(α1,α2) and the normal vector n(α1,α2) to the reference surface at the given point (Figure 1). Moreover, the position of the generic point of the shell volume is also a function of the shell thickness h(α1,α2). Hence, writing the position vector of the reference surface, it is possible to define the three components along the three global axes Ox1x2x3 as
(3)r(α1,α2)=r1(α1,α2)e1+r2(α1,α2)e2+r3(α1,α2)e3,
where e1,e2, and e3 are the unit vector of the global reference system Ox1x2x3. From the definition of the first fundamental form [7, 11, 12, 23]  of the reference surface r(α1,α2), the Lamé parameters can be expressed as
(4)A1(α1,α2)=r,1·r,1,A2(α1,α2)=r,2·r,2.
The comma in expressions (4) defines the partial derivative with respect to α1 or α2, respectively. Moreover, by considering an orthogonal curvilinear coordinate system O′α1α2ζ [7, 11, 12, 23, 74] from the position vector (3) the normal vector n(α1,α2) can be written as
(5)n(α1,α2)=r,1×r,2A1A2.
Finally, due to the fact that an orthogonal curvilinear coordinate system O′α1α2ζ is considered and following the definition of the second fundamental form [7, 11, 23, 74] of the reference surface r(α1,α2), the principal radii of curvature can be evaluated as
(6)R1(α1,α2)=-r,1·r,1r,11·n,R2(α1,α2)=-r,2·r,2r,22·n.
As concerns the shell theory, the present study is based on the following assumptions: (1) the transverse normal is inextensible so that the normal strain is equal to zero: εn=εn(α1,α2,ζ,t)=0; (2) the transverse shear deformation is considered to influence the governing equations so that normal lines to the reference surface of the shell before deformation remain straight but not necessarily normal after deformation (a relaxed Kirchhoff-Love hypothesis is considered); (3) the shell deflections are small and the strains are infinitesimal; (4) the shell is moderately thick, therefore it is possible to assume that the thickness-direction normal stress is negligible so that the plane assumption can be invoked: σn=σn(α1,α2,ζ,t)=0; (5) the linear elastic behavior of anisotropic materials is assumed; and (6) the rotary inertias and the initial curvatures are also taken into account. The theory under development is meant to be in the group of moderately thick shells, in which the following ratios of thickness over curvature and curve length are valid:
(7)1100≤max(hRmin,hLmin)≤110.
Consistent with the assumptions of a moderately thick shell theory reported above, the displacement field considered in the present study is that of the First-order Shear Deformation Theory (FSDT) and can be put in the following form:
(8)U1(α1,α2,ζ,t)=H1u1(α1,α2,t)+ζβ1(α1,α2,t),U2(α1,α2,ζ,t)=H2u2(α1,α2,t)+ζβ2(α1,α2,t),U3(α1,α2,ζ,t)=u3(α1,α2,t),
where
(9)H1=1+ζR1,H2=1+ζR2,
and u1,u2, and u3 are the displacement components of points lying on the middle surface (ζ=0) of the shell, while t is the time variable. β1 and β2 are normal-to-mid-surface rotations, respectively. The kinematic hypothesis expressed by relations (8) should be supplemented by the statement that the shell deflections are small and strains are infinitesimal, that is u3≪h. In-plane displacements U1 and U2 vary linearly through the thickness, while U3 remains independent of ζ. Differently from the previous works [43, 44, 47, 49–53, 67, 68, 70], the displacement field has been improved taking into account the effective geometry of the shell and in particular the curvature effect has been directly introduced into the kinematical model, as proposed by Toorani and Lakis [21].Figure 1
Geometry description and coordinate system of a doubly curved shell.Due to the change of the kinematical model, the relationships between strains and generalized displacements along the shell reference surface(ζ=0) become the following: (10)[ε10ε20γ10γ20χ10χ20ω10ω20γ1n0γ2n0]=[1A1∂∂α11A1A2∂A1∂α21R1001A1A2∂A2∂α11A2∂∂α21R200-1A1A2∂A1∂α21A1∂∂α10001A2∂∂α2-1A1A2∂A2∂α10001A1R1∂∂α1-1A1R12∂R1∂α11A1A2R2∂A1∂α201A1∂∂α11A1A2∂A1∂α21A1A2R1∂A2∂α11A2R2∂∂α2-1A2R22∂R2∂α201A1A2∂A2∂α11A2∂∂α2-1A1A2R1∂A1∂α21A1R2∂∂α1-1A1R22∂R2∂α10-1A1A2∂A1∂α21A1∂∂α11A2R1∂∂α2-1A2R12∂R1∂α2-1A1A2R2∂A2∂α101A2∂∂α2-1A1A2∂A2∂α1001A1∂∂α110001A2∂∂α201][u1u2u3β1β2]that are different from those presented in previous papers [43, 44, 47, 49–53, 67, 68, 70].  In the above (10), the first four strains ε10,ε20,γ10, and γ20 are the in-plane components, and χ10,χ20,ω10, and ω20 are the analogous curvature changes. The last two components γ1n0,γ2n0 are the transverse shearing strains. The shell is assumed to be made of a linear elastic composite material. Accordingly, the following constitutive equations relate internal stress resultants and internal couples with generalized strain components (10) on the middle surface:(11)[N1N2N12N21M1M2M12M21T1T2]=[A11(20)(0)A12(11)(0)A16(20)(0)A16(11)(0)A11(20)(1)A12(11)(1)A16(20)(1)A16(11)(1)00A12(11)(0)A22(02)(0)A26(11)(0)A26(02)(0)A12(11)(1)A22(02)(1)A26(11)(1)A26(02)(1)00A16(20)(0)A26(11)(0)A66(20)(0)A66(11)(0)A16(20)(1)A26(11)(1)A66(20)(1)A66(11)(1)00A16(11)(0)A26(02)(0)A66(11)(0)A66(02)(0)A16(11)(1)A26(02)(1)A66(11)(1)A66(02)(1)00A11(20)(1)A12(11)(1)A16(20)(1)A16(11)(1)A11(20)(2)A12(11)(2)A16(20)(2)A16(11)(2)00A12(11)(1)A22(02)(1)A26(11)(1)A26(02)(1)A12(11)(2)A22(02)(2)A26(11)(2)A26(02)(2)00A16(20)(1)A26(11)(1)A66(20)(1)A66(11)(1)A16(20)(2)A26(11)(2)A66(20)(2)A66(11)(2)00A16(11)(1)A26(02)(1)A66(11)(1)A66(02)(1)A16(11)(2)A26(02)(2)A66(11)(2)A66(02)(2)0000000000κA44(20)(0)κA45(11)(0)00000000κA45(11)(0)κA55(02)(0)][ε10ε20γ10γ20χ10χ20ω10ω20γ1n0γ2n0],where the elastic engineering stiffnesses Aij(pq)(τ) are defined as follows:
(12)Aij(pq)(τ)=∑k=1l∫ζkζk+1Q-ij(k)ζτH1H2H1pH2qdζ,τ,p,q=0,1,2.
Different approaches can be found in the literature to evaluate the engineering elastic constants Aij(pq)(τ) [7, 17–20, 73, 76]. In the present paper, the relations of the elastic engineering stiffnesses Aij(pq)(τ) are numerically evaluated using an integral function in order to avoid numerical instabilities. It is worth noting that due to the fact that the elastic engineering stiffnesses Aij(pq)(τ) depend on curvatures, the corresponding derivatives, with respect to the coordinates along α1 and α2 directions of the reference surface, have to be evaluated. In order to perform this operation the Differential Quadrature rule [26] is used. Thus, the derivatives of the elastic engineering stiffnesses Aij(pq)(τ) are numerically evaluated. The corresponding elastic constants Q-ij(k) can be found in the article by Tornabene et al. [70, 73], in which all the constants above introduced are explicitly defined for laminated composite and functionally graded shells and panels. κ is the shear correction factor, which is usually taken as κ=5/6, such as in the present work. In particular, the determination of shear correction factors for composite laminated structures is still an unresolved issue, since these factors depend on various parameters [20]. In (11), the four components N1,N2,N12, and N21 are the in-plane force resultants, and M1,M2,M12, and M21 are the analogous couples, while T1,T2 are the transverse shear force resultants. In the above definitions (11) the symmetry of shearing force resultants N12,N21 and torsional couples M12,M21 is not assumed as a further hypothesis, as done in Reissner-Mindlin theory [7, 12, 22]. This hypothesis is satisfied only in the case of spherical shells and flat plates [7]. The assumption under discussion is derived from the consideration that ratios ζ/R1,ζ/R2 cannot be neglected with respect to unity.Typically, the functionally graded materials are made of a mixture of two constituents. In the present work, it is assumed that the functionally graded material lamina is made of a mixture of ceramic and metal constituents: Silicon Nitride and Stainless Steel. The material properties of the functionally graded shell vary continuously and smoothly in the thickness directionζ of each lamina and are functions of volume fractions of two constituent materials. The Young’s modulus E(k)(ζ), shear modulus G(k)(ζ), Poisson’s ratio ν(k)(ζ), and mass density ρ(k)(ζ) of the functionally graded shell kth lamina can be expressed as
(13)ρ(k)(ζ)=(ρC(k)-ρM(k))VC(k)(ζ)+ρM(k),E(k)(ζ)=(EC(k)-EM(k))VC(k)(ζ)+EM(k),ν(k)(ζ)=(νC(k)-νM(k))VC(k)(ζ)+νM(k),G(k)(ζ)=E(k)(ζ)2(1+v(k)(ζ)),forζk≤ζ≤ζk+1,
where  ρC(k),EC(k),νC(k),VC(k) and ρM(k),EM(k),νM(k),VM(k) represent mass density, Young’s modulus, Poisson’s ratio, and volume fraction of the ceramic and metal constituent materials, respectively. In this work, the ceramic volume fraction VC(k)(ζ) follows two simple four-parameter power-law distributions [49, 51, 53, 70, 75]:
(14)FGM1(a(k)/b(k)/c(k)/p(k)):VC(k)(ζ)=(1-a(k)(ζhk-ζkhk)+b(k)(ζhk-ζkhk)c(k))p(k),FGM2(a(k)/b(k)/c(k)/p(k)):VC(k)(ζ)=(1-a(k)(ζk+1hk-ζhk)+b(k)(ζk+1hk-ζhk)c(k))p(k),
where the volume fraction index p(k)(0≤p(k)≤∞) and the parameters a(k),b(k), and c(k) dictate the material variation profile through the functionally graded shell lamina thickness. It is important to remark that the volume fractions of all the constituent materials should add up to unity:
(15)VC(k)+VM(k)=1.
In order to choose the three parameters a(k),b(k), and c(k) suitably, the relation (15) must be always satisfied for every volume fraction index p(k) in each lamina. By considering the relations (14), when the power-law exponent is set equal to zero (p(k)=0) or equal to infinity (p(k)=∞), the homogeneous isotropic material is obtained as a special case of functionally graded material. In fact, from (15), (14), and (13) it is possible to obtain
(16)p(k)=0⟶VC(k)=1,VM(k)=0⟶ρ(k)(ζ)=ρC,E(k)(ζ)=EC,ν(k)(ζ)=νC,p(k)=∞⟶VC(k)=0,VM(k)=1⟶ρ(k)(ζ)=ρM,E(k)(ζ)=EM,ν(k)(ζ)=νM.
Some material profiles through the functionally graded shell thickness are illustrated in Figures 2 and 3.Variations of the ceramic volume fractionVC through the thickness for different values of the three parameters a(1),b(1), and c(1) and the power-law index p(1) for a single-layered shell: (a) FGM1(a(1)=1/b(1)=0/c(1)/p(1)), (b) FGM2(a(1)=1/b(1)=0/c(1)/p(1)), (c) FGM1(a(1)=1/b(1)=1/c(1)=2/p(1)), (d) FGM2(a(1)=1/b(1)=1/c(1)=2/p(1)), (e) FGM1(a=1/b=0.5/c=2/p), and (f) FGM2(a(1)=1/b(1)=0.5/c(1)=2/p(1)).
(a)
(b)
(c)
(d)
(e)
(f)Variations of the ceramic volume fractionVC through the thickness for two-layered and three-layered laminated shells for different values of the power-law index p=p(1)=p(2) and p=p(1)=p(3): (a) FGM1(a(1)=1/b(1)=0/c(1)/p(1))/FGM2(a(2)=1/b(2)=0/c(2)/p(2)), (b) FGM2(a(1)=1/b(1)=0/c(1)/p(1))/FGM1(a(2)=1/b(2)=0/c(2)/p(2)), (c) FGM2(a(1)=1/b(1)=0/c(1)/p(1))/FGMC/FGM1(a(3)=1/b(3)=0/c(3)/p(3)),  and (d) FGM2(a(1)=1/b(1)=0.5/c(1)=2/p(1))/FGMC/FGM1(a(3)=1/b(3)=0.5/c(3)=2/p(3)).
(a)
(b)
(c)
(d)Following the Hamilton’s principle [7, 12, 19, 22, 23], the five governing equations in terms of internal actions can be written for the shell element:
(17)1A1A2∂(N1A2)∂α1+1A1A2∂(N21A1)∂α2+N12A1A2∂A1∂α2-N2A1A2∂A2∂α1+1R1(1A1A2∂(M1A2)∂α1+1A1A2∂(M21A1)∂α2+M12A1A2∂A1∂α2-M2A1A2∂A2∂α11A1A2∂(M1A2)∂α1)+q1=(I0+2I1R1+I2R12)u¨1+(I1+I2R1)β¨1,1A1A2∂(N12A2)∂α1+1A1A2∂(N2A1)∂α2+N21A1A2∂A2∂α1-N1A1A2∂A1∂α2+1R2(1A1A2∂(M12A2)∂α1+1A1A2∂(M2A1)∂α2+M21A1A2∂A2∂α1-M1A1A2∂A1∂α21A1A2∂(M12A2)∂α1)+q2=(I0+2I1R2+I2R22)u¨2+(I1+I2R2)β¨2,1A1A2∂(T1A2)∂α1+1A1A2∂(T2A1)∂α2-(N1R1+N2R2)+qn=I0u¨3,1A1A2∂(M1A2)∂α1+1A1A2∂(M21A1)∂α2+M12A1A2∂A1∂α2-M2A1A2∂A2∂α1-T1+m1=(I1+I2R1)u¨1+I2β¨1,1A1A2∂(M12A2)∂α1+1A1A2∂(M2A1)∂α2+M21A1A2∂A2∂α1-M1A1A2∂A1∂α2-T2+m2=(I1+I2R2)u¨2+I2β¨2,
where
(18)Iτ=∑k=1l∫ζkζk+1ρ(k)ζτH1H2dζ,τ=0,1,2.
Furthermore, the generalized external actions q1,q2,qn,m1, and m2 due to the external forces, acting on the top and bottom surfaces of the shell, can be evaluated using the static equivalence principle [7, 23] and can be written on the reference surface of the doubly curved shell as follows:
(19)q1=q1+(1+h2R1)2(1+h2R2)+q1-(1-h2R1)2(1-h2R2),q2=q2+(1+h2R1)(1+h2R2)2+q2-(1-h2R1)(1-h2R2)2,qn=qn+(1+h2R1)(1+h2R2)+qn-(1-h2R1)(1-h2R2),m1=q1+h2(1+h2R1)(1+h2R2)-q1-h2(1-h2R1)(1-h2R2),m2=q2+h2(1+h2R1)(1+h2R2)-q2-h2(1-h2R1)(1-h2R2),
where q1+,q1-,q2+,q2-,qn+, and qn- are the external forces in the three principal directions α1,α2, and ζ at the top and the bottom surface of the shell, respectively. The three basic sets of equations, namely, the kinematic (10), constitutive (11), and motion (17) equations, may be combined to give the fundamental system of equations, also known as the governing system of equations. By replacing the kinematic equations (10) into the constitutive equations (11) and the result of this substitution into the motion equations (17), the complete equations of motion in terms of displacement and rotational components can be written as(20)[L11L12L13L14L15L21L22L23L24L25L31L32L33L34L35L41L42L43L44L45L51L52L53L54L55][u1u2u3β1β2]+[q1q2qnm1m2]=[I0100I1100I0200I1200I000I1100I200I1200I2][u¨1u¨2u¨3β¨1β¨2],where Lij,i,j=1,…,5 are the equilibrium operators and the new mass inertias are defined as follows:
(21)I01=I0+2I1R1+I2R12,I11=I1+I2R1,I02=I0+2I1R2+I2R22,I12=I1+I2R2.
Differently from previous works [43, 44, 47, 49–53, 67, 68, 70], the equilibrium operators Lij, introduced in (20), have been changed due to the choice of using the kinematical model (8).Three kinds of boundary conditions are considered, namely, the fully clamped edge boundary condition(C), the simply supported edge boundary condition (S), and the free edge boundary condition (F). The equations describing the boundary conditions can be written as follows:Clamped edge boundary conditions (C):
(22)
u
1
=
u
2
=
u
3
=
β
1
=
β
2
=
0
at
α
1
=
α
1
0
or
α
1
=
α
1
1
,
α
2
0
≤
α
2
≤
α
2
1
,
u
1
=
u
2
=
u
3
=
β
1
=
β
2
=
0
at
α
2
=
α
2
0
or
α
2
=
α
2
1
,
α
1
0
≤
α
1
≤
α
1
1
,Simply supported edge boundary conditions (S):
(23)
u
1
=
u
2
=
u
3
=
β
2
=
0
,
M
1
=
0
at
α
1
=
α
1
0
or
α
1
=
α
1
1
,
α
2
0
≤
α
2
≤
α
2
1
,
u
1
=
u
2
=
u
3
=
β
2
=
0
,
M
2
=
0
at
α
2
=
α
2
0
or
α
2
=
α
2
1
,
α
1
0
≤
α
1
≤
α
1
1
,Free edge boundary conditions (F):
(24)
N
1
+
M
1
R
1
=
0
,
N
12
+
M
12
R
2
=
0
,
T
1
=
0
,
M
1
=
M
12
=
0
at
α
1
=
α
1
0
or
α
1
=
α
1
1
,
α
2
0
≤
α
2
≤
α
2
1
,
(25)
N
2
+
M
2
R
2
=
0
,
N
21
+
M
21
R
1
=
0
,
T
2
=
0
,
M
2
=
M
21
=
0
at
α
2
=
α
2
0
or
α
2
=
α
2
1
,
α
1
0
≤
α
1
≤
α
1
1
.
In addition to the external boundary conditions (22)–(25), the kinematic and physical compatibility conditions should be satisfied at the common closing meridians with α2=0,2π, if a complete shell of revolution has to be considered. The kinematic compatibility conditions include the continuity of displacements. The physical compatibility conditions can only be represented by the five continuous conditions for the generalized stress resultants. To consider complete revolute shells characterized by α21=2π, it is necessary to implement the kinematic and physical compatibility conditions between the two computational meridians with α20=0 and with α21=2π:Kinematic compatibility conditions along the closing meridian
(
α
2
=
0,2
π
)
(26)
u
1
(
α
1
,
0
,
t
)
=
u
1
(
α
1
,
2
π
,
t
)
,
u
2
(
α
1
,
0
,
t
)
=
u
2
(
α
1
,
2
π
,
t
)
,
u
3
(
α
1
,
0
,
t
)
=
u
3
(
α
1
,
2
π
,
t
)
,
β
1
(
α
1
,
0
,
t
)
=
β
1
(
α
1
,
2
π
,
t
)
,
β
2
(
α
1
,
0
,
t
)
=
β
2
(
α
1
,
2
π
,
t
)
α
1
0
≤
α
1
≤
α
1
1
,
Physical compatibility conditions along the closing meridian
(
α
2
=
0,2
π
)
(27)
N
2
(
α
1
,
0
,
t
)
+
M
2
(
α
1
,
0
,
t
)
R
2
=
N
2
(
α
1
,
2
π
,
t
)
+
M
2
(
α
1
,
2
π
,
t
)
R
2
,
N
21
(
α
1
,
0
,
t
)
+
M
21
(
α
1
,
0
,
t
)
R
1
=
N
21
(
α
1
,
2
π
,
t
)
+
M
21
(
α
1
,
2
π
,
t
)
R
1
,
T
2
(
α
1
,
0
,
t
)
=
T
2
(
α
1
,
2
π
,
t
)
,
M
2
(
α
1
,
0
,
t
)
=
M
2
(
α
1
,
2
π
,
t
)
,
M
21
(
α
1
,
0
,
t
)
=
M
21
(
α
1
,
2
π
,
t
)
,
α
1
0
≤
α
1
≤
α
1
1
.
In analogous way, in order to consider a toroidal shell of revolution it is necessary to implement the kinematic and physical compatibility conditions between the two computational parallels with α10=0 and with α11=2π:Kinematic compatibility conditions along the closing parallel
(
α
1
=
0,2
π
)
(28)
u
1
(
0
,
α
2
,
t
)
=
u
1
(
2
π
,
α
2
,
t
)
,
u
2
(
0
,
α
2
,
t
)
=
u
2
(
2
π
,
α
2
,
t
)
,
u
3
(
0
,
α
2
,
t
)
=
u
3
(
2
π
,
α
2
,
t
)
,
β
1
(
0
,
α
2
,
t
)
=
β
1
(
2
π
,
α
2
,
t
)
,
β
2
(
0
,
α
2
,
t
)
=
β
2
(
2
π
,
α
2
,
t
)
,
α
2
0
≤
α
2
≤
α
2
1
,Physical compatibility conditions along the closing parallel
(
α
1
=
0,2
π
)
(29)
N
1
(
0
,
α
2
,
t
)
+
M
1
(
0
,
α
2
,
t
)
R
1
=
N
1
(
2
π
,
α
2
,
t
)
+
M
1
(
2
π
,
α
2
,
t
)
R
1
,
N
12
(
0
,
α
2
,
t
)
+
M
12
(
0
,
α
2
,
t
)
R
2
=
N
12
(
2
π
,
α
2
,
t
)
+
M
12
(
2
π
,
α
2
,
t
)
R
2
,
T
1
(
0
,
α
2
,
t
)
=
T
1
(
2
π
,
α
2
,
t
)
,
M
1
(
0
,
α
2
,
t
)
=
M
1
(
2
π
,
α
2
,
t
)
,
M
12
(
0
,
α
2
,
t
)
=
M
12
(
2
π
,
α
2
,
t
)
α
2
0
≤
α
2
≤
α
2
1
.


## 3. Discretized Equations and Numerical Implementation

The Generalized Differential Quadrature method is used to discretize the spatial derivatives in the governing equations in terms of generalized displacements, as well as boundary conditions (see Tornabene [49] for a brief review). Throughout the paper, the Chebyshev-Gauss-Lobatto (C-G-L) grid distribution is assumed, for which the coordinates of grid points (α1i,α2j) along the reference surface are in the discrete form:
(30)α1i=(1-cos(i-1N-1π))(α11-α10)2+α10,i=1,2,…,N,forα1∈[α10,α11],α2j=(1-cos(j-1M-1π))(α21-α20)2+α20,j=1,2,…,M,forα2∈[α20,α21],
where N,M are the total number of sampling points used to discretize the domain in α1 and α2 directions, respectively, of the doubly curved shell. It has been proven that, for the Lagrange interpolating polynomials, the Chebyshev-Gauss-Lobatto sampling points rule guarantees convergence and efficiency to the GDQ technique [23, 43–45, 67, 68]. For the static analysis, when the inertias (21) are set to zero, the GDQ procedure enables to write the governing equations (20) and the boundary and compatibility conditions (22)–(29) in discrete form, transforming each space derivative into a weighted sum of node values of independent variables using the Differential Quadrature rule [26, 49]:
(31)∂nf(x)∂xn|x=xm=∑k=1Tςmk(n)f(xk),m=1,2,…,T.
Each approximate equation is valid in a single sampling point. Thus, the whole system of differential equations has been discretized and the global assembling leads to the following set of linear algebraic equations:
(32)[KbbKbdKdbKdd][δbδd]=[fbfd].
In the above mentioned matrices and vectors, the partitioning is set forth by subscripts b and d, referring to the system degrees of freedom and standing for boundary and domain, respectively. In this sense, b-equations represent the discrete boundary conditions, which are valid only for the points lying on the constrained edges of the shell, while d-equations are the equilibrium equations, assigned on the interior nodes. In order to make the computation more efficient, static condensation of nondomain degrees of freedom is performed:
(33)(Kdd-KdbKbb-1Kbd)δd=fd-KdbKbb-1fb.
The deflection of the considered structures can be determined by solving the linear algebraic problem (33). In particular, the solution procedure by means of the GDQ technique has been implemented in a personal code.Differently from the static case, when the external forcesq1+,q1-,q2+,q2-,qn+, and qn- (19) are set to zero, the free vibration of laminated composite doubly curved shells and panels can be studied. Using the method of variable separation, it is possible to seek solutions that are harmonic in time and whose frequency is f=ω/2π. The generalized displacements can be written as follows:
(34)u1(α1,α2,t)=U1(α1,α2)eiωt,u2(α1,α2,t)=U2(α1,α2)eiωt,u3(α1,α2,t)=U3(α1,α2)eiωt,β1(α1,α2,t)=B1(α1,α2)eiωt,β2(α1,α2,t)=B2(α1,α2)eiωt,
where the vibration spatial amplitude values U1,U2,U3,B1, and B2 fulfil the fundamental differential system (20). Each approximate equation is valid in a single sampling point. Thus, the whole system of differential equations can be discretized and the global assembling leads to a set of linear eigenvalue problem. When kinematic condensation of nondomain degrees of freedom is performed, one gets
(35)(Kdd-Kdb(Kbb)-1Kbd)δd=ω2Mddδd.
The natural frequencies of the structure fr=ωr/2π, for r=1,2,…,5(N-2)×(M-2), can be determined by solving the standard eigenvalue problem (35). In particular, the solution procedure by means of the GDQ technique has been implemented in a personal code. The above partitioning (35) is set forth by subscripts b and d, referring to the system degrees of freedom and standing for boundary and domain, respectively. Finally, the results in terms of frequencies are obtained using an eigenvalue function. With the present approach, differing from the finite element method, no integration occurs prior to the global assembly of the linear system, and this represents a further computational cost saving in favour of the Differential Quadrature technique.


## 4. Optimization Algorithms

### 4.1. Genetic Algorithm (GA) Optimization Method

The GA approach to optimization was probably introduced at first by Holland [107], while a comprehensive reference including implementation procedures and application notes can be found in later work, like, for instance, the seminal book by Goldberg [99]. GA tries to implement and imitate in a mathematical framework the law of evolution, which Darwin introduced to explain the changes in nature towards individuals better suited to the environment in which they live: GA can be in fact classified as a typical population-based optimization algorithm. Generation by generation, populations evolve improving the fitness function which represents the ability to survive in a defined environment; the individuals less adapting to the surrounding environment do not mate and their genetic set of chromosomes is lost. A solution is represented by a chromosome, constituted by a set of genes representing the parameters of the solution. The mathematical implementation follows with the coding of a chromosome in binary and with the application of some computational functions like mutations, crossover, and elitarism. When elitarism is set, some of the best individuals are replicated in the following generation without any change in a perfect replication from father to son. The concept of crossover implies a change between genes of two solutions and imitates the reproduction in which the son possesses a part of genes from his mother and the remaining from the father. Aim of mutations, as it happens in nature, is to randomly change some genes of the individual to test new configurations: from an optimization point of view, mutations help in exploring new zones of the space and are useful to avoid the problem of “local minimum” capturing. According to the work by Konak et al. [106], the fitness is the main driver of the capability to survive and to pass genes to the next generation: the chromosome is in fact decoded from binary to decimal and tested in the fitness only after the application of mating functions requiring a binary coding. Modern approaches to the application of GA lie in new formulations and in the introduction of hybridization with other optimization strategies [108]; also Pareto-based analysis for multiobjective optimization [103] has been evaluated, and improvement of GA by the application of fuzzy sets and neural networks [109] has been proposed to solve complex tasks. The algorithm implemented in this work follows the procedure proposed by Goldberg [99], and the fitness has been defined by authors considering the closeness of the solution found to the value set by the designer and the feasibility of the volume fraction distribution in the thickness (0<VC(k)<1). The pseudocode of the GA implemented is shown in the following  Pseudocode 1.Pseudocode 1:Genetic Algorithm pseudocode.
Step(1) Set a max and min value for each gene
Step(2)  Set the max number of generations (T), the convergence tolerance (ε) and the number of consecutive no improvements
loops (nε) after which algorithm ends, the number of genes (ng), and the number of population (npop) members
Step(3) Initialize the first generation of population (P) by randomly set npop chromosomes, each one made by ng genes
Step(4) while (i<T) or (ε for nε loops false)
Step(5)    Select P′⊂P (mating pool), initialize P′′=0 (set of children)
Step(6)     for j=1 to n
Step(7)     Randomly select individuals (chromosomes xa and xb) from P′
Step(8)      Obtain xchild by applying crossover to xa and xb (with probability pcross)
Step(9)      Mutate produced child xchild to xa and xb (with probability pmut)
Step(10)      Apply elitarism (if set)
Step(11)      Update population P′′=P′∪xchild
Step(12)     end for
Step(13)P = survival respect to the fitness* (P′, P′′)
Step(14)  end while
Step(15) best chromosome detection
*In this case the fitness is multiplied by a penalty term if one of the sets [a(k),b(k),c(k),p(k)] leads to volume fractions which are
not feasible (e.g. percentage of one constituent in the thickness larger than 1, less than 0, imaginary number).


### 4.2. Monte Carlo (MC) Optimization Method

The Monte Carlo technique has been developed to support early studies related to the nuclear physics; the idea is to find the best approximation of a constant or to solve a problem through a statistical way, obtaining the results from a very large set of random inputs. From a mathematical and formal point of view, for instance, the paper by Mosegaard and Sambridge [123] introduces the way by which an integral in the form
(36)I=∫Xg(x)f(x)dx
can be evaluated by the generation of random samples x1,x2,x3,x4,…,xS of x, when f(x) is an appropriate  probability distribution and g(x) is the function of which the integral have to be computed. The integral can be so computed by the expression
(37)I≈1S∑s=1Sg(xs).
The MC methods are very simple, but studies are focused on the probability distribution shapes providing best results and on the software methods to generate random numbers; this is not a trivial issue, since it can be complex to generate a set of random numbers which are not dependent on the clock of the processor or on other hardware timers. Obviously, considering computational efforts, the MC method is very time expensive; however, due to the increasing computations capabilities of personal computers and the improving of the capabilities which are a constant trend in years, this approach has been recently reconsidered and still applied to a wide range of applications. The continuous improvement in randomizer algorithms helps in the gain of good performances, which are often obtained exploiting the capabilities of clusters of computers. The main critical aspect related to the MC methods application in engineering is that it requires a lot of iteration to obtain a solution; following this approach, in fact, the space is explored in a completely random way and no attention is focused on zones in which solutions are better than the average. The literature presents a lot of interesting applications of MC to engineering problems, like the works [124–126] show; MC is also widely applied to games and strategy since it allows keeping in consideration all the possible scenarios evolving from a situation: some papers dealing with this issue can be found in [127]. A pseudocode of the MC algorithm is included in Pseudocode 2.Pseudocode 2:Monte Carlo pseudocode.
Step(1) Set a max and min value for each parameter
Step(2) Set the max number of population (ng) members, the convergence tolerance (ε) and
the number of consecutive no improvements loops (nε)
Step(3) while (i<ng) or (ε for nε loops false)
Step(4)   Randomly generate a member of population Pi, following a probability
distribution (normal, exponential, Weibull [113]).
Step(5)   Compute volume fraction distribution (VC) in the thickness for Pi
Step(6)   if 0<VC<1
Step(7)   i=i+1
Step(8)     Compute Fitness
Step(9)     if Fitness (i) > Best Fitness
Step(10)     Best Fitness = Fitness (i)
Step(11)     Best member = i
Step(12)    end if
Step(12)    end if
Step(13) end while


### 4.3. Particle Swarm Optimization (PSO) Method

The PSO algorithm used in this paper is similar to the one proposed by Birge [115], with some variations due to the particular application to the FGM optimization problem. According to the general implementation of PSO, the position of a particle in the n-dimensional solution space can be considered a solution to the problem in the n-dimensional space of the parameters; the n-dimensional speed which is computed in the algorithm represents the direction towards a new position in the n-dimensional space. The velocity can be obtained by
(38)vi(m)=ϕ(m)vi(m-1)+c1a1(Pbi-Pi(m-1))+c2a2(PbGb-Pi(m-1)),
where m is the algorithm step (nG generations); i is the index of parameter of the single particle; ϕ(m) is the inertia function; vi(m) is the velocity of the ith particle at the m-step; Pbi is the best position found by the ith particle; Gb is the global best position (it is the best position found by the whole swarm); a1,a2 are the acceleration constants; c1,c2 are random numbers in the interval [0,1].The following formulation can be applied to obtain the new position:(39)Pi(m)=Pi(m-1)+vi(m),
where Pi(m) is the position of the ith particle at the m-step and vi(m) is the velocity of the ith particle at the m-step. One of the problems of the FGM profile shape in the thickness is that a small variation in a(k),b(k),c(k), and p(k) power-law parameters can lead to a solution showing a volume fraction inconsistent (VC(k)<0 or VC(k)>1). The PSO algorithm, in fact, belongs to the so-called “trajectory based methods” and sweeps the space following a path in which the new position is equal to the previous one plus a constant segment (the speed). For this reason, sometimes the updated position of a particle lies in a zone in which the volume fraction is out of limit; in this case, the classical implementation of the algorithm stops since the trajectory enters in a trap from which it is impossible to escape. The new position, in fact, can lead to inconsistent volume fractions and the algorithm stalls. In order to solve this problem, a check has been introduced in the algorithm (step 11 and step 12 in Pseudocode 3), so that if an unfeasible position is found by the algorithm, the velocity is rejected and changed until a new valid position is found. By this way, a forecast of the new position is computed and the PSO velocity is accepted if leading to a valid solution, randomized if it does not. The end criterion is due to one of these conditions: the achievement of the maximum number of iterations, or a condition in which after nε generations the solution does not improve. The PSO seems to be very attractive for optimization since the studies by Hu et al. [116] and by Ceruti et al. [102] present advantages with respect to GA: easiness of software code implementation, need for the definition of few setting parameters, and good fitting to engineering application. The PSO implementation and the modifications introduced for this application are described in Pseudocode 3, where an illustrative pseudocode is presented.Pseudocode 3:Particle Swarm Optimization pseudocode.
Step(1) Set number of particles (nP), number of parameters for each particle (nparam), max
number of generations (nG), convergence tolerance (ε) and the number of
consecutive no improvements loops (nε), δ (ratio between position and speed)
Step(2) Set particles boundary values: max(1:nparam), min(1:nparam)
Step(3) Randomly generate particles position (Pi(0))∣min(i)<Pi(0)<max(i) and velocity vi(0)
Step(4) Evaluate fitness function Fi=f(Pi(0)),i=1,…,nP
Step(5) Set Pb=P(0), Fbest=F, Gb=min(Fi), gb=indexofthe(min(Fi))value
Step(6) for m=1 to nG
Step(7)  for i=1 to nP
Step(8)   Δvi=c1a1(Pbi-Pi(m-1))+c2a2(PbGb-Pi(m-1))
Step(9)   vi(m)=ϕ(m)v(m-1)+Δvi
Step(10)   Pi(m)=Pi(m-1)+vi(m)
Step(11)    Compute Volume fraction distribution (VC) in the thickness for Pi(m)
Step(12)     while (Pi(m)>max(i))or(Pi(m)<min(i))or(VC<0)or(VC>1)
Step(13)      vi(m)=v(m-1)+δ·(max(i)-min(i))·random(0÷1)
Step(14)      Pi(m)=Pi(m-1)+vi(m)
Step(15)      end while
Step(16)    Evaluate fitness  Fi=f(Pi(m)),i=1,…,nP
Step(17)      for i=1,…,nP
Step(18)     if Fi<Fbi then Pbi=Pi(m) and Fbi=Fi
Step(19)     if Fi<Gb then Gb=Fi and gb=i
Step(20)      end for
Step(21)  end for
Step(22) end for


## 4.1. Genetic Algorithm (GA) Optimization Method

The GA approach to optimization was probably introduced at first by Holland [107], while a comprehensive reference including implementation procedures and application notes can be found in later work, like, for instance, the seminal book by Goldberg [99]. GA tries to implement and imitate in a mathematical framework the law of evolution, which Darwin introduced to explain the changes in nature towards individuals better suited to the environment in which they live: GA can be in fact classified as a typical population-based optimization algorithm. Generation by generation, populations evolve improving the fitness function which represents the ability to survive in a defined environment; the individuals less adapting to the surrounding environment do not mate and their genetic set of chromosomes is lost. A solution is represented by a chromosome, constituted by a set of genes representing the parameters of the solution. The mathematical implementation follows with the coding of a chromosome in binary and with the application of some computational functions like mutations, crossover, and elitarism. When elitarism is set, some of the best individuals are replicated in the following generation without any change in a perfect replication from father to son. The concept of crossover implies a change between genes of two solutions and imitates the reproduction in which the son possesses a part of genes from his mother and the remaining from the father. Aim of mutations, as it happens in nature, is to randomly change some genes of the individual to test new configurations: from an optimization point of view, mutations help in exploring new zones of the space and are useful to avoid the problem of “local minimum” capturing. According to the work by Konak et al. [106], the fitness is the main driver of the capability to survive and to pass genes to the next generation: the chromosome is in fact decoded from binary to decimal and tested in the fitness only after the application of mating functions requiring a binary coding. Modern approaches to the application of GA lie in new formulations and in the introduction of hybridization with other optimization strategies [108]; also Pareto-based analysis for multiobjective optimization [103] has been evaluated, and improvement of GA by the application of fuzzy sets and neural networks [109] has been proposed to solve complex tasks. The algorithm implemented in this work follows the procedure proposed by Goldberg [99], and the fitness has been defined by authors considering the closeness of the solution found to the value set by the designer and the feasibility of the volume fraction distribution in the thickness (0<VC(k)<1). The pseudocode of the GA implemented is shown in the following  Pseudocode 1.Pseudocode 1:Genetic Algorithm pseudocode.
Step(1) Set a max and min value for each gene
Step(2)  Set the max number of generations (T), the convergence tolerance (ε) and the number of consecutive no improvements
loops (nε) after which algorithm ends, the number of genes (ng), and the number of population (npop) members
Step(3) Initialize the first generation of population (P) by randomly set npop chromosomes, each one made by ng genes
Step(4) while (i<T) or (ε for nε loops false)
Step(5)    Select P′⊂P (mating pool), initialize P′′=0 (set of children)
Step(6)     for j=1 to n
Step(7)     Randomly select individuals (chromosomes xa and xb) from P′
Step(8)      Obtain xchild by applying crossover to xa and xb (with probability pcross)
Step(9)      Mutate produced child xchild to xa and xb (with probability pmut)
Step(10)      Apply elitarism (if set)
Step(11)      Update population P′′=P′∪xchild
Step(12)     end for
Step(13)P = survival respect to the fitness* (P′, P′′)
Step(14)  end while
Step(15) best chromosome detection
*In this case the fitness is multiplied by a penalty term if one of the sets [a(k),b(k),c(k),p(k)] leads to volume fractions which are
not feasible (e.g. percentage of one constituent in the thickness larger than 1, less than 0, imaginary number).


## 4.2. Monte Carlo (MC) Optimization Method

The Monte Carlo technique has been developed to support early studies related to the nuclear physics; the idea is to find the best approximation of a constant or to solve a problem through a statistical way, obtaining the results from a very large set of random inputs. From a mathematical and formal point of view, for instance, the paper by Mosegaard and Sambridge [123] introduces the way by which an integral in the form
(36)I=∫Xg(x)f(x)dx
can be evaluated by the generation of random samples x1,x2,x3,x4,…,xS of x, when f(x) is an appropriate  probability distribution and g(x) is the function of which the integral have to be computed. The integral can be so computed by the expression
(37)I≈1S∑s=1Sg(xs).
The MC methods are very simple, but studies are focused on the probability distribution shapes providing best results and on the software methods to generate random numbers; this is not a trivial issue, since it can be complex to generate a set of random numbers which are not dependent on the clock of the processor or on other hardware timers. Obviously, considering computational efforts, the MC method is very time expensive; however, due to the increasing computations capabilities of personal computers and the improving of the capabilities which are a constant trend in years, this approach has been recently reconsidered and still applied to a wide range of applications. The continuous improvement in randomizer algorithms helps in the gain of good performances, which are often obtained exploiting the capabilities of clusters of computers. The main critical aspect related to the MC methods application in engineering is that it requires a lot of iteration to obtain a solution; following this approach, in fact, the space is explored in a completely random way and no attention is focused on zones in which solutions are better than the average. The literature presents a lot of interesting applications of MC to engineering problems, like the works [124–126] show; MC is also widely applied to games and strategy since it allows keeping in consideration all the possible scenarios evolving from a situation: some papers dealing with this issue can be found in [127]. A pseudocode of the MC algorithm is included in Pseudocode 2.Pseudocode 2:Monte Carlo pseudocode.
Step(1) Set a max and min value for each parameter
Step(2) Set the max number of population (ng) members, the convergence tolerance (ε) and
the number of consecutive no improvements loops (nε)
Step(3) while (i<ng) or (ε for nε loops false)
Step(4)   Randomly generate a member of population Pi, following a probability
distribution (normal, exponential, Weibull [113]).
Step(5)   Compute volume fraction distribution (VC) in the thickness for Pi
Step(6)   if 0<VC<1
Step(7)   i=i+1
Step(8)     Compute Fitness
Step(9)     if Fitness (i) > Best Fitness
Step(10)     Best Fitness = Fitness (i)
Step(11)     Best member = i
Step(12)    end if
Step(12)    end if
Step(13) end while


## 4.3. Particle Swarm Optimization (PSO) Method

The PSO algorithm used in this paper is similar to the one proposed by Birge [115], with some variations due to the particular application to the FGM optimization problem. According to the general implementation of PSO, the position of a particle in the n-dimensional solution space can be considered a solution to the problem in the n-dimensional space of the parameters; the n-dimensional speed which is computed in the algorithm represents the direction towards a new position in the n-dimensional space. The velocity can be obtained by
(38)vi(m)=ϕ(m)vi(m-1)+c1a1(Pbi-Pi(m-1))+c2a2(PbGb-Pi(m-1)),
where m is the algorithm step (nG generations); i is the index of parameter of the single particle; ϕ(m) is the inertia function; vi(m) is the velocity of the ith particle at the m-step; Pbi is the best position found by the ith particle; Gb is the global best position (it is the best position found by the whole swarm); a1,a2 are the acceleration constants; c1,c2 are random numbers in the interval [0,1].The following formulation can be applied to obtain the new position:(39)Pi(m)=Pi(m-1)+vi(m),
where Pi(m) is the position of the ith particle at the m-step and vi(m) is the velocity of the ith particle at the m-step. One of the problems of the FGM profile shape in the thickness is that a small variation in a(k),b(k),c(k), and p(k) power-law parameters can lead to a solution showing a volume fraction inconsistent (VC(k)<0 or VC(k)>1). The PSO algorithm, in fact, belongs to the so-called “trajectory based methods” and sweeps the space following a path in which the new position is equal to the previous one plus a constant segment (the speed). For this reason, sometimes the updated position of a particle lies in a zone in which the volume fraction is out of limit; in this case, the classical implementation of the algorithm stops since the trajectory enters in a trap from which it is impossible to escape. The new position, in fact, can lead to inconsistent volume fractions and the algorithm stalls. In order to solve this problem, a check has been introduced in the algorithm (step 11 and step 12 in Pseudocode 3), so that if an unfeasible position is found by the algorithm, the velocity is rejected and changed until a new valid position is found. By this way, a forecast of the new position is computed and the PSO velocity is accepted if leading to a valid solution, randomized if it does not. The end criterion is due to one of these conditions: the achievement of the maximum number of iterations, or a condition in which after nε generations the solution does not improve. The PSO seems to be very attractive for optimization since the studies by Hu et al. [116] and by Ceruti et al. [102] present advantages with respect to GA: easiness of software code implementation, need for the definition of few setting parameters, and good fitting to engineering application. The PSO implementation and the modifications introduced for this application are described in Pseudocode 3, where an illustrative pseudocode is presented.Pseudocode 3:Particle Swarm Optimization pseudocode.
Step(1) Set number of particles (nP), number of parameters for each particle (nparam), max
number of generations (nG), convergence tolerance (ε) and the number of
consecutive no improvements loops (nε), δ (ratio between position and speed)
Step(2) Set particles boundary values: max(1:nparam), min(1:nparam)
Step(3) Randomly generate particles position (Pi(0))∣min(i)<Pi(0)<max(i) and velocity vi(0)
Step(4) Evaluate fitness function Fi=f(Pi(0)),i=1,…,nP
Step(5) Set Pb=P(0), Fbest=F, Gb=min(Fi), gb=indexofthe(min(Fi))value
Step(6) for m=1 to nG
Step(7)  for i=1 to nP
Step(8)   Δvi=c1a1(Pbi-Pi(m-1))+c2a2(PbGb-Pi(m-1))
Step(9)   vi(m)=ϕ(m)v(m-1)+Δvi
Step(10)   Pi(m)=Pi(m-1)+vi(m)
Step(11)    Compute Volume fraction distribution (VC) in the thickness for Pi(m)
Step(12)     while (Pi(m)>max(i))or(Pi(m)<min(i))or(VC<0)or(VC>1)
Step(13)      vi(m)=v(m-1)+δ·(max(i)-min(i))·random(0÷1)
Step(14)      Pi(m)=Pi(m-1)+vi(m)
Step(15)      end while
Step(16)    Evaluate fitness  Fi=f(Pi(m)),i=1,…,nP
Step(17)      for i=1,…,nP
Step(18)     if Fi<Fbi then Pbi=Pi(m) and Fbi=Fi
Step(19)     if Fi<Gb then Gb=Fi and gb=i
Step(20)      end for
Step(21)  end for
Step(22) end for


## 5. Numerical Applications and Results

In the present section, some results and considerations about the mixed static and dynamic optimization problem of four-parameter functionally graded doubly curved and degenerate shells and panels are presented. The analysis has been carried out by means of numerical procedures illustrated above. Different types of structures are considered in the present paper. One of the aims of this study is to show some numerical examples about flat plates, singly curved, and doubly curved shells and panels made of FGMs. The six considered structures are depicted in Figure4. In order to describe the middle surface of the given structures, theoretical formulae of differential geometry [7, 23, 43, 74, 78] are used. The mathematical development of the differential geometry applied to doubly curved shells was deeply explained in [23, 43, 74, 78]. So, in the following, only a few formulae are reported. Furthermore, it is worthwhile noting that the GDQ procedure enables to evaluate the parameters concerning the shell geometry as reported in [74, 78]. For all the GDQ results presented below, the Chebyshev-Gauss-Lobatto grid distributions (30) with N=M=31 along the reference surface have been assumed, and the geometrical parameters with their derivatives are numerically evaluated using the GDQ method [26]. The FGM nomenclature used in the present work is identified by the same convention presented in the previous work by Tornabene et al. [70]. In an analogous way, the geometrical boundary conditions are defined considering similar convention used in the previous works [23, 43, 44, 47, 49–53, 68–76] for shell and panel structures.Six different types of FGM shells and panels: (a) rectangular plate (degenerate panel), (b) cylindrical panel (singly curved panel), (c) conical shell (singly curved shell), (d) toroidal shell panel (doubly curved shell panel of revolution), (e) catenoidal panel (doubly curved panel of revolution),and (f) elliptic paraboloid (completely doubly curved panel).
(a)
(b)
(c)
(d)
(e)
(f)For a rectangular flat plate, the position vector [74, 78] can be written as
(40)r(α1,α2)=-α2e2+α1e3.
In Figure 4(a) the considered square plate has the sides a=b=1m, thickness h=0.1m, and it is subjected to a normal load qn+=-10000Pa at the top surface. The plate is free on two adjacent sides and clamped on the other two. The boundary conditions are indicated by the current nomenclature as CCFF. The square plate is a single-layered structure with FGM1(a(1)/b(1)/c(1)/p(1)) four-parameter power-law distribution. In this case, the four parameters a=a(1),b=b(1),c=c(1), and p=p(1) have to be determined by the optimization procedure.The position vector for the cylindrical panel and conical shell of Figures4(b) and 4(c) can be obtained from the conical shell formula, already presented in [74, 78] and reported here for the sake of completeness. The position vector for the conical shell can be written as follows:
(41)r(α1,α2)=(Rb+α1sinα)cosα2e1-(Rb+α1sinα)sinα2e2+α1cosαe3.
As it is well known [7, 74, 78], a conical shell [23, 42, 43, 49, 52] is a 2D structure having the middle surface generated by the rotation, about a fixed vertical axis, of an inclined straight line. The angle α is the top vertex angle of the cone, where the meridian angle φ can be written as a function of it: φ=π/2-α. Thus, for cylindrical panel, considering α=0, the position vector (41) takes the form
(42)r(α1,α2)=Rbcosα2e1-Rbsinα2e2+α1e3.
The meridian abscissa α1 is defined as α1∈[0,L], where L is the length of the meridian. The presented parameters are graphically shown in the works [23, 42, 43, 49, 52]. The cylindrical panel (Figure 4(b)) is defined by R=Rb=2m, L=3m, circumferential angle α2=ϑ∈[0,120°], and a thickness h=0.2m. It is subjected to a normal load qn+=-10000Pa at the top surface, and its boundary conditions are SFSF. On the contrary, the conical shell (Figure 4(c)) has Rb=1m, L=3m, circumferential angle α2=ϑ∈[0,360°], α=30°, and a thickness h=0.3m. Also the conical shell is subjected to a normal load qn+=-10000Pa at the top surface. The conical shell is clamped at the bottom and free at the top, so the boundary conditions are indicated by CF. The cylindrical panel is a two-layered structure with FGM2(a(1)/b(1)/c(1)/p(1))/FGM1(a(2)/b(2)/c(2)/p(2)) lamination scheme, while the conical shell presents FGM1(a(1)/b(1)/c(1)/p(1))/FGM2(a(2)/b(2)/c(2)/p(2)) lamination scheme. In these two cases the FGM four parameters of the two different power-law distributions are assumed to be equal to the two different laminae of the structure. Thus, only the four parameters a=a(1)=a(2),b=b(1)=b(2),c=c(1)=c(2), and p=p(1)=p(2) have to be determined by the optimization procedure.As already reported in [74, 78], the position vector of a toroidal shell panel (Figure 4(d)) can be written as
(43)r(α1,α2)=(Rb+Rsinα1)cosα2e1-(Rb+Rsinα1)sinα2e2+Rcosα1e3,
where Rb is the shift of the circular meridian curve with respect to the axis of revolution [23, 43, 44, 47, 68, 69] and R is the radius of the circular curve section of the toroidal shell panel. The toroidal structure is characterized by Rb=9m, R=3m, meridian angle α1=φ∈[0,360°], circumferential angle α2=ϑ∈[0,120°], and a thickness h=0.6m. Also the toroidal shell has CF boundary conditions, and the normal load qn+=-10000Pa is applied at its top surface. When a catenary curve is considered as a meridian curve of a revolution shell, the position vector of the catenoidal panel (Figure 4(e)) assumes the aspect
(44)r(α1,α2)=(Acosh(arcsinh(1tanα1))+Rb)cosα2e1-(Acosh(arcsinh(1tanα1))+Rb)sinα2e2+Aarcsinh(1tanα1)e3,
where A is the distance of the throat apex of the catenary curve. For further details about the geometry definition of the catenary curve, the reader might refer to [23, 43, 74, 78]. Regarding the catenoidal panel (Figure 4(e)) its geometrical properties are A=2m, B=2m, D=1m, Rb=0m, circumferential angle α2=ϑ∈[0,90°], and a thickness h=0.5m. The catenoidal panel has CFCF boundary conditions, and it is subjected to a normal load qn+=-10000Pa at the top surface. The toroidal shell panel and catenary panel of revolution have the same three-layered lamination scheme FGM2(a(1)/b(1)/c(1)/p(1))/FGMC/FGM1(a(3)/b(3)/c(3)/p(3)). The lamination scheme presents the middle lamina of the structure made of ceramic isotropic material. Also, in these two cases the FGM parameters of the two different power-law distribution are assumed to be equal to the first and the third laminae of the structure. In this way, only the four parameters a=a(1)=a(3),b=b(1)=b(3),c=c(1)=c(3), and p=p(1)=p(3) have to be determined by the optimization procedure.Finally, the position vector of the reference surface of the elliptic paraboloid is given in [74, 78]:
(45)r(α1,α2)=(R0α1(α1)-x3α2(α2)sinα1)cosα2e1-R0α2(α2)e2+(x3α1(α1)+x3α2(α2)cosα1)e3.
If the generatrix parabola needs to be characterized (the other parabola can be defined analogously), it can be described by the following equation:
(46)(R0α1)2-κα1x3α1=0,
where κα1=(A2-D2)/B is a characteristic parameter of the parabolic curve, R0α1 is the abscissa of a point of the parabola, and x3α1 is its ordinate in the generatrix plane of the parabolic curve. The abscissa R0α1(α1) of the parabolic curve assumes the form
(47)R0α1(α1)=κα1tanα12.
The parameters describing the two parabolas of the elliptic paraboloid under consideration are A=3m, C=-3m, D=0m, and B=0.8m, and the thickness is h=0.4m. The elliptic paraboloid is completely clamped at its four edges (CCCC), and it is subjected to a normal load qn+=-10000Pa at the top surface. The elliptic paraboloid is a single-layered structure with FGM2(a(1)/b(1)/c(1)/p(1)) four-parameter power-law distribution. In this case, the four parameters a=a(1),b=b(1),c=c(1), and p=p(1) have to be determined by the optimization procedure.In Table1 the first ten frequencies and maximum static deflections for the six structures of Figure 4 are presented. In particular, the two limit cases of functionally graded materials are considered. For the first case, named FGMC, the single-layered structure is made of ceramic isotropic material, while for the second one, named FGMM, the single-layered structure is made of metal isotropic material. In the present work, it is assumed that the functionally graded material lamina is made of a mixture of ceramic and metal constituents: Silicon Nitride and Stainless Steel. Young’s modulus, Poisson’s ratio, and mass density for the Silicon Nitride are EC=322.27GPa, νC=0.24, and ρC=2370kg/m3, and for the Stainless Steel are EM=207.78GPa, νM=0.3177, and ρM=8166kg/m3, respectively.Table 1
First ten frequencies and maximum static deflections for the six shell structures of Figure4 made of isotropic materials: FGMC: full ceramic material and FGMM: full metal material using a 31×31 Chebyshev-Gauss-Lobatto (C-G-L) grid distribution with different boundary conditions.
Square plate CCFF
Cylindrical panel SFSF
Conical shell CF
Toroidal shell panel CF
Catenoidal panel CFCF
Elliptic paraboloid CCCC
Mode [Hz]
FG
M
C
FG
M
M
FG
M
C
FG
M
M
FG
M
C
FG
M
M
FG
M
C
FG
M
M
FG
M
C
FG
M
M
FG
M
C
FG
M
M
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
GDQ (31×31)
f
1
375.598
162.089
271.645
115.701
449.402
194.323
28.665
12.452
828.011
359.407
451.863
201.635
f
2
1252.975
536.895
282.473
120.922
449.402
194.323
29.385
12.679
846.255
365.918
485.074
211.268
f
3
1373.081
602.852
546.801
234.445
585.799
247.965
80.330
34.397
1042.714
451.381
487.409
212.299
f
4
2367.239
1019.777
549.102
239.032
585.799
247.965
80.416
34.706
1136.811
486.278
632.857
276.825
f
5
2868.597
1235.869
642.802
276.494
680.296
298.116
113.948
49.115
1147.395
492.825
719.612
314.418
f
6
3084.992
1334.067
648.822
279.014
680.296
298.116
121.533
52.684
1212.360
518.534
747.215
328.314
f
7
3194.352
1399.648
847.611
371.648
860.268
373.956
133.553
57.452
1476.266
632.458
864.072
377.506
f
8
3629.667
1590.719
933.540
404.836
860.268
373.956
160.575
70.062
1566.446
679.318
878.506
383.706
f
9
4018.936
1736.048
973.695
423.554
930.739
411.564
167.115
71.979
1890.222
821.620
1019.096
443.351
f
10
4156.483
1800.873
1143.969
483.642
966.531
418.914
197.566
85.108
1941.698
833.159
1034.972
450.206
u
3
max
1.565687
E
-
05
2.456358
E
-
05
1.012876
E
-
05
1.623087
E
-
05
4.773661
E
-
07
7.045025
E
-
07
3.440893
E
-
05
5.244087
E
-
05
3.511216
E
-
07
5.452771
E
-
07
1.819132
E
-
06
2.611289
E
-
06As described before, one of the goals of the present paper is to develop and test an operative methodology to support the design of FG materials. Several problems can be faced by the designer; in particular some of the most commonly occurred engineering needs are to avoid resonance frequency, to reduce or control the deformations, to increase the strength to weight ratio, to reduce stresses due to different dilatation properties of two materials, and so on. In this study, frequency and deformation and a combination of the two will be considered. The fitness functionF has been modelled to obtain a defined first resonance frequency target (f1T) set by user, a maximum displacement target (wmaxT=u3maxT) due to bending deformation, or to try to satisfy at the best both the requirements. In this latter case, the designer can decide to stress more the importance on one of the two requirements by using the coefficient η which acts as a weight: in this way the fitness F is still monodimensional but keeps into account two aspects, ranging η from zero to one. According to the above considerations, the fitness expression F can be so represented by the formula
(48)F=η(f1-f1Tf1T)2+(1-η)(wmax-wmaxTwmaxT)2.
It is worth to note that the minimization of the fitness F is equivalent to obtain values of frequency (f=f1) and of displacement (w=wmax=u3max) as close as possible to the desired ones (f1T,wmaxT=u3maxT), and the ratio between the difference from the obtained value minus the target and the target value itself is used to make dimensionless frequencies and displacements. By this way, parameters with different magnitudes can be compared, allowing to equate items  with unit of measurement defined by the experimenter depending on his practice. In the paper, this simple monoobjective fitness has been considered to privilege the simplicity and to suggest the use of this approach also in industrial applications, but more complex multiobjective fitness functions can be implemented depending on the specific case study and on the requirements and needs expressed by the designer.Figures5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, and 22 present the results of the optimization of the six structures described above. A mix of values of η parameter has been chosen (η=0,0.5,1), keeping it constant for each structure as to make possible the comparison of the optimization methods. Following formula (48), if η=1 only the frequency is optimized, while for η=0 only the displacement is considered; an intermediate value η=0.5 can be set to find the better compromise between the guess of both frequency and displacement. Each figure presents in the upper left corner the values of a and b parameter for each try; the best solution is represented by a red circle. In the upper centre there are the values of c and p parameters for each try with the best solution in red; in the upper left corner there is the material distribution in the thickness of the structure; an FGM made by one, two, or three layers can be visualized, depending on the structure definition (e.g. rectangular plate presents only one layer; cylindrical panel presents two symmetric layers, while, for instance, the conical shell is made by a sandwich between an homogeneous material and two symmetric about the midplane layers of FG materials): the final optimized distribution is in red. The lower part of the figures from left to right presents the trend of frequency, displacement, and fitness in time. In case of Monte Carlo frequency and displacement of the single run are shown (with in blue the trend of the best solution found), while in case of Particle Swarm and Genetic Algorithms the frequency and the displacement referred to the best individual found at the moment are shown. A green line represents the target values for the frequency and the displacement, while in case of the fitness, the ideal value is obviously equal to zero.Figure 5
Square plate: Monte Carlo optimization withη=1; final optimized frequency f=265.364Hz.Figure 6
Square plate: Genetic Algorithm optimization withη=1; final optimized frequency f=268.816Hz.Figure 7
Square plate: Particle Swarm Optimization withη=1; final optimized frequency f=268.710Hz.Figure 8
Cylindrical panel: Monte Carlo optimization withη=1; final optimized displacement w=1.1596·10-5m.Figure 9
Cylindrical panel: Genetic Algorithm optimization  withη=1; final optimized displacement w=1.1668·10-5m.Figure 10
Cylindrical panel: Particle Swarm Optimization withη=1; final optimized displacement w=1.1653·10-5m.Figure 11
Conical shell: Monte Carlo optimization withη=0.5; final optimized frequency f=330.946Hz and final optimized displacement w=5.2073·10-7m.Figure 12
Conical shell: Genetic Algorithm optimization withη=0.5; final optimized frequency f=321.890Hz and final optimized displacement w=5.2624·10-7m.Figure 13
Conical shell: optimization by Particle Swarm Algorithm withη=0.5; final optimized frequency f=322.049Hz and final optimized displacement w=5.2615·10-7m.Figure 14
Toroidal shell panel: Monte Carlo optimization withη=0; final optimized displacement w=3.9727·10-5m.Figure 15
Toroidal shell panel: Genetic Algorithm optimization withη=0; final optimized displacement w=3.9495·10-5m.Figure 16
Toroidal shell panel: Particle Swarm Optimization withη=0; final optimized displacement w=3.9519·10-5m.Figure 17
Catenoidal panel: Monte Carlo optimization withη=1; final optimized frequency f=595.868Hz.Figure 18
Catenoidal panel: Genetic Algorithm optimization withη=1; final optimized frequency f=590.800Hz.Figure 19
Catenoidal panel: Particle Swarm optimization withη=1; final optimized frequency f=590.119Hz.Figure 20
Elliptic paraboloid: Monte Carlo optimization withη=0.5; final optimized frequency f=332.688Hz and final optimized displacement w=1.8923·10-6m.Figure 21
Elliptic paraboloid: Genetic Algorithm optimization withη=0.5; final optimized frequency f=327.600Hz and final optimized displacement w=1.9075·10-6m.Figure 22
Elliptic paraboloid: Particle Swarm Optimization withη=0.5; final optimized frequency f=327.398Hz and final optimized displacement w=1.9065·10-6m.Figures5–7 present the optimization of a rectangular plate when the designer is interested in setting the first resonance frequency of the structure. The plate is made by only one FGM lamina, and a value of η=1 is set. Figures 8–10 depict the result of the optimization with MC, PSO, and GA of a cylindrical panel made by two symmetrical FGM distributions about the midplane; in this case η=0 is set since the need for a required maximum displacement is simulated. Figures 11–13 refer to a singly curved conical shell made by two FGM laminae, with a distribution symmetric about the midplane of the whole structure; in this case η=0.5 is set to consider an interest both in frequency and in maximum displacement. Figures 14–16 present a toroidal shell panel which is a doubly curved shell panel of revolution, made by three laminae: a homogeneous core and two external symmetrical FGM laminae as to simulate a typical sandwich structure. The interest is in maximum displacement, so that η is set equal to zero. Figures 17–19 list the optimization trend of a catenoidal panel structure which can be defined as a doubly curved panel of revolution; the laminate is obtained by three layers also in this case, representing a typical sandwich structure; in this case the attention is focused on frequency optimization, so that η=1 is set. Finally, Figures 20–22 present an elliptic paraboloid which is a completely doubly curved panel structure made by only one FGM lamina; considerations will be addressed both on frequency and displacement since η=0.5.As to help the interpretation of the optimization process, the final results are collected in Table2: the final value of the fitness provides an idea of the goodness of the results found. It is important to note that heuristic methods do not provide the best solution, but a suboptimal solution; in other words, sometimes, similar results of the fitness function can be found with several different set of parameters a,b,c, and p. The Genetic Algorithm has been run with a population of 30 individuals and the Particle Swarm with a number of particles of 30; these algorithms stop when a maximum number of runs are reached or no improvements are found in solution (or alternatively when the percentage error is less than the 0.5%). The Monte Carlo simulation has been run for 200 times for each optimization problem, in order to obtain a comparable number of runs with Genetic Algorithm and Particle Swarm Optimizers. A first comparison between methods can be performed based on Figure 23, in which there is a graph collecting final fitness versus computational time for each optimization performed. Values for Monte Carlo are in red triangles, Particle Swarm in blue squares, while Genetic Algorithm in red circles. As expected GA and PSO results are better than MC ones and provide data with good approximation in a relative short time, which is also compatible with the needs of conceptual design or Multidisciplinary Optimization. Simulations have been carried out on a Notebook equipped with an 8 GB RAM, a 2 GHz core processor, and Windows 7 operating system. Table 3 lists the errors in frequency and in displacement (depending on the case study) for each simulation run; as can be seen by data provided, both GA and PSO present a similar efficiency. The two methods exploit different strategies to solve the optimization problem in this complex task: in this case, in fact, the domain is not continuous and some combinations of the four parameters a,b,c, and p can lead to a distribution of metal and ceramics which is not consistent (percentage represented by a negative or complex or major than one number, which is physically unfeasible). As well described before in the paper, one of the challenges is that slight changes in one of the four parameters can lead to inconsistent solutions. When a particle of the PSO arrives to the border of the domain with a speed pointing to an inconsistent zone of the domain, it is needed to randomly redirect the particle towards an allowed zone, and this reduces the efficiency of the method. On the other hand, GA solves in a better way this problem since mutations and crossover help in exploring new zones of the domain in which the solution is feasible and the distribution percentage of the FGM distributions (ceramics and metal) is always a positive number between zero and one in all the thickness of the structure. Only a small part of the individuals of a generation should mutate because it slows down the convergence. After tests performed it seems that these two effects compensate and both PSO and GA present a very good behavior in optimization. The comparison between random optimization (MC) with heuristic methods (GA, PSO) shows a different relative efficiency of MC when applied to cases in which only frequency or displacement is studied or when both are considered. In the first case the difference in precision obtained by MC with respect to GA and PSO is quite small, while this difference increases dramatically when both frequency and displacement are considered: it can be inferred that when the fitness function deals with only one parameter also the MC alone can be used, while in case of fitness involving more parameters an heuristic or semi heuristic method provides better results.Table 2
Optimization of the six different structures of Figure4: (a) square plate, (b) cylindrical panel, (c) conical shell, (d) toroidal shell panel, (e) catenoidal panel, and (f) elliptic paraboloid. η=1: optimization of the first frequency f1;  η=0: optimization of maximum bending displacement wmax=u3max; η=0.5: optimization of the first frequency f1 and maximum bending displacement wmax=u3max.
Method
a
b
c
p
Final frequency
Final displacement
Final fitness
Computational time (min)
Rectangular plate
η
=
1: target frequency f1T=268.753 Hz
MC
0.05582
−0.56831
17.86455
9.14279
265.36499
1.8270
E
-
05
1.2606
E
-
02
25
GA
0.08800
−0.06160
5.02440
9.53080
268.81630
1.7493
E
-
05
2.3583
E
-
04
42
PSO
−0.25829
−0.13011
0.67170
3.43176
268.71000
1.7466
E
-
05
1.5970
E
-
04
31
Cylindrical panel
η
=
0: target displacement wmaxT=u3maxT=1.1646·10-5m
MC
1.11895
1.10588
1.09493
13.33798
186.78559
1.1596
E
-
05
4.3053
E
-
03
32
GA
0.12710
−0.33040
6.11930
3.03030
202.41150
1.1668
E
-
05
1.8770
E
-
03
14
PSO
0.35190
−0.26277
6.90518
1.75727
197.75230
1.1653
E
-
05
5.8904
E
-
04
20
Conical shell
η
=
0.5: target frequency f1T=268.753Hz and target displacement wmaxT=u3maxT=5.2639·10-7m
MC
0.00984
−0.74086
14.77422
27.83370
330.94622
5.2073
E
-
07
1.9799
E
-
02
35
GA
0.08800
−0.32650
16.65690
7.72240
321.89420
5.2624
E
-
07
4.9477
E
-
04
25
PSO
0.05769
−0.49331
9.85753
8.35549
322.04900
5.2615
E
-
07
8.2088
E
-
04
43
Toroidal shell panel
η
=
0: target displacement wmaxT=u3maxT=3.9516·10-5m
MC
0.05350
−0.54251
17.60320
20.64836
20.21863
3.9727
E
-
05
5.3396
E
-
03
21
GA
0.02930
−0.78200
12.72730
25.90420
20.56750
3.9495
E
-
05
5.3143
E
-
04
3
PSO
0.10430
−0.58284
8.53859
8.22430
20.48775
3.9519
E
-
05
7.5919
E
-
05
6
Catenoidal panel
η
=
1: target frequency f1T=590.557Hz
MC
1.06122
0.90588
2.06831
2.50895
595.86858
3.9753
E
-
07
8.9937
E
-
03
55
GA
0.08800
−0.48390
13.11830
13.97850
590.77480
4.0135
E
-
07
3.6832
E
-
04
45
PSO
0.09089
−0.32061
18.33638
14.94431
590.11938
4.0129
E
-
07
7.4152
E
-
04
43
Elliptic paraboloid
η
=
0.5: target frequency f1T=327.239Hz and target displacement wmaxT=u3maxT=1.9073·10-6m
MC
0.29917
0.02997
19.54311
2.65521
332.68799
1.8923
E
-
06
1.2258
E
-
02
43
GA
0.08800
−0.29910
3.01080
4.93650
327.62800
1.9075
E
-
06
6.4729
E
-
04
22
PSO
0.25269
−0.65120
13.49764
2.90612
327.39815
1.9065
E
-
06
4.5339
E
-
04
25Table 3
Perceptual relative errors in first frequtencyef=((f1-f1T)/f1T)·100 and in maximum bending displacement ew=((wmax-wmaxT)/wmaxT)·100 with Monte Carlo, Particle Swarm Optimization, and Genetic Algorithm methods.
Optimization method
Square plate
Cylindrical panel
Conical shell
Toroidal shell panel
Catenoidal panel
Elliptic paraboloid
e
f
e
w
e
f
e
w
e
w
e
f
e
f
e
w
MC
−1.261
−0.431
2.884
−1.076
0.534
0.899
1.665
−0.786
GA
0.024
0.188
0.070
−0.029
−0.053
0.037
0.119
0.010
PSO
−0.016
0.059
0.118
−0.046
0.008
−0.074
0.049
−0.042Figure 23
Computational optimization time [in minutes] versus fitness for MC (triangle), GA (circle), and PSO (square).


## 6. Conclusion Remarks and Summary

Functionally graded material applications are increasing due to their potentialities in problems related to aggressive environments, thermal properties, and special structural needs. Theoretical frameworks have been developed to solve the static and dynamic analysis of this type of structures. In the present paper, the Generalized Differential Quadrature method has been presented as a mean to investigate the static and dynamic analysis of functionally graded and laminated composite doubly curved shells and panels. The adopted shell theory is the First-order Shear Deformation Theory. In particular, the Toorani-Lakis theory has been used as a starting point to obtain the governing equations for shells. The 2D equilibrium equations have been discretized with the GDQ method through standard linear algebraic and eigenvalue problems. Thanks to the mathematical approach, the analysis of an FGM shell or panel can be provided. The material distribution can be expressed by a four-parameter power-law. However, the main drawback of this approach from a designer perspective is that it is impossible to analytically relate the set of the four parameters to static or dynamic performances of the structures. An approach based on the application of optimization algorithms to the problem has been carried out to evaluate its suitability to the topic. Due to the domain of the four parameters which includes zones in which no solutions can be found, the optimization algorithms have been modified and tested. Genetic Algorithm, Particle Swarm Optimization, and Monte Carlo algorithms have been selected for tests and applied to six different geometrical structures optimizing the first frequency, maximum bending displacement, or a mix of the two. Results obtained show the convergence of the optimization for GA and PSO, which provides better results than MC for similar computational effort. In particular, GA and PSO are more precise when both frequency and displacement are considered at the same time. The herein developed methodology can dramatically help the designer in the definition of the mix ceramics/metal which can present the requested properties.



---

*Source: 867079-2013-05-22.xml*