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

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 andmethodologies are implemented.TheParticle SwarmOptimization,Monte Carlo andGenetic 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.


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, twodimensional linear theories of thin shells and plates have been developed including important contributions [1][2][3][4][5][6][7][8][9][10][11][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 Mathematical Problems in Engineering 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  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][100][101][102][103][104][105][106][107][108][109], Differential Evolution [110], Particle Swarm Algorithm [111][112][113][114][115][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][124][125][126][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.

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 2 ) 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 (−ℎ/2 ≤  ≤ ℎ/2), where ℎ( 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: where  = 2/ℎ( 1 ,  2 ) and  ∈ [−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  laminae or plies, the total thickness ℎ is defined as in which ℎ  =  +1 −   is the thickness of the th 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 ℎ( 1 ,  2 ).Hence, writing the position vector of the reference surface, it is possible to define the three components along the three global axes  1  2  3 as where e 1 , e 2 , and e 3 are the unit vector of the global reference system  1  2  3 .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 The comma in expressions (4) defines the partial derivative with respect to  1 or  2 , respectively.Moreover, by considering an orthogonal curvilinear coordinate system    1  2  [7,11,12,23,74] from the position vector (3) the normal vector n( 1 ,  2 ) can be written as Finally, due to the fact that an orthogonal curvilinear coordinate system    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 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:   =   ( 1 ,  2 , , ) = 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:   =   ( 1 ,  2 , , ) = 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: 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: where and  1 ,  2 , and  3 are the displacement components of points lying on the middle surface ( = 0) of the shell, while  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  3 ≪ ℎ.In-plane displacements  1 and  2 vary linearly through the thickness, while  3 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].
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  () (), shear modulus  () (), Poisson's ratio ] () (), and mass density  () () of the functionally graded shell th lamina can be expressed as where  ()  ,  ()  , ] ()  ,  ()  and  ()  ,  ()  , ] ()  ,  ()  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  ()   () Mathematical Problems in Engineering follows two simple four-parameter power-law distributions [49,51,53,70,75]: ) ) where the volume fraction index  () (0 ≤  () ≤ ∞) and the parameters  () ,  () , and  () 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: In order to choose the three parameters  () ,  () , and  ()  suitably, the relation ( 15) must be always satisfied for every volume fraction index  () in each lamina.By considering the relations ( 14), when the power-law exponent is set equal to zero ( () = 0) or equal to infinity ( () = ∞), 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 Some material profiles through the functionally graded shell thickness are illustrated in Figures 2 and 3.
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: where Furthermore, the generalized external actions  1 ,  2 ,   ,  1 , and  2 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: ) , where  + 1 ,  − 1 ,  + 2 ,  − 2 ,  +  , and  −  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 where   , ,  = 1, . . ., 5 are the equilibrium operators and the new mass inertias are defined as follows: Differently from previous works [43, 44, 47, 49-53, 67, 68, 70], the equilibrium operators   , 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 (), the simply supported edge boundary condition (), and the free edge boundary condition ().The equations describing the boundary conditions can be written as follows:

Clamped edge boundary conditions (C):
Simply supported edge boundary conditions (S): Free edge boundary conditions (F): 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  1 2 = 2, it is necessary to implement the kinematic and physical compatibility conditions between the two computational meridians with  0 2 = 0 and with  1 2 = 2:  Figure 2: Variations of the ceramic volume fraction   through the thickness for different values of the three parameters  (1) ,  (1) , and  (1) and the power-law index  (1) for a single-layered shell: (a) FGM

Physical compatibility conditions along the closing meridian
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  0 1 = 0 and with  1 1 = 2: Kinematic compatibility conditions along the closing parallel  (1) =  (2) and  =  (1) =

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 ( 1 ,  2 ) along the reference surface are in the discrete form: where ,  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]: 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: In the above mentioned matrices and vectors, the partitioning is set forth by subscripts  and , referring to the system degrees of freedom and standing for boundary and domain, respectively.In this sense, -equations represent the discrete boundary conditions, which are valid only for the points lying on the constrained edges of the shell, while -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: 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 forces  + 1 ,  − 1 ,  + 2 ,  − 2 ,  +  , and  −  (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  = /2.The generalized displacements can be written as follows: where the vibration spatial amplitude values  1 ,  2 ,  3 ,  1 , and  2 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 The natural frequencies of the structure   =   /2, for  = 1, 2, . . ., 5 ( − 2) × ( − 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  and , 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.

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 populationbased 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 <  ()  < 1).The pseudocode of the GA implemented is shown in the following Pseudocode 1.

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 can be evaluated by the generation of random samples  1 ,  2 ,  3 ,  4 , . . .,   of , when () is an appropriate probability distribution and () is the function of which the integral have to be computed.The integral can be so computed by the expression 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][125][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.

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 -dimensional solution space can be considered a solution to the problem in the -dimensional space of the parameters; the -dimensional speed which is computed in the algorithm represents the direction towards a new position in the -dimensional space.The velocity can be obtained by where  is the algorithm step (  generations);  is the index of parameter of the single particle;  () is the inertia function; V ()  is the velocity of the th particle at the -step;   is the best position found by the th particle;   is the global best position (it is the best position found by the whole swarm);  1 ,  2 are the acceleration constants;  1 ,  2 are random numbers in the interval [0, 1].
The following formulation can be applied to obtain the new position: where  ()  is the position of the th particle at the -step and V ()  is the velocity of the th particle at the -step.One of the problems of the FGM profile shape in the thickness is that a small variation in  () ,  () ,  () , and  () power-law parameters can lead to a solution showing a volume fraction inconsistent ( ()   < 0 or  ()  > 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   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.

Numerical Applications and Results
In the present section, some results and considerations about the mixed static and dynamic optimization problem of fourparameter 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 Figure 4.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  =  = 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][50][51][52][53][68][69][70][71][72][73][74][75][76] for shell and panel structures.For a rectangular flat plate, the position vector [74,78] can be written as In Figure 4(a) the considered square plate has the sides  =  = 1m, thickness ℎ = 0.1 m, and it is subjected to a normal load  +  = −10000 Pa 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 FGM 1( (1) / (1) / (1) / (1) ) four-parameter powerlaw distribution.In this case, the four parameters  =  (1) ,  =  (1) ,  =  (1) , and  =  (1) have to be determined by the optimization procedure.
The position vector for the cylindrical panel and conical shell of Figures 4(b) and 4(c) can be obtained from the conical shell formula, already presented in [74,78] and reported here Step (1) Set a max and min value for each gene Step (2) Set the max number of generations (), the convergence tolerance () and the number of consecutive no improvements loops (  ) after which algorithm ends, the number of genes (  ), and the number of population ( pop ) members Step (3) Initialize the first generation of population () by randomly set  pop chromosomes, each one made by   genes Step (4) while ( < ) or ( for   loops false) Step (5) Select   ⊂  (mating pool), initialize   = 0 (set of children) Step (6) for  = 1 to  Step (7) Randomly select individuals (chromosomes   and   ) from   Step (8) Obtain  child by applying crossover to   and   (with probability  cross ) Step (9) Mutate produced child  child to   and   (with probability  mut ) Step (10) Apply elitarism (if set) Step (11) Update population   =   ∪  child Step (12) end for Step (13) P = survival respect to the fitness * (  ,   ) 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 [ () ,  () ,  () ,  () ] 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).
As already reported in [74,78], the position vector of a toroidal shell panel (Figure 4(d)) can be written as where   is the shift of the circular meridian curve with respect to the axis of revolution [23,43,44,47,68,69] and  is the radius of the circular curve section of the toroidal shell panel.The toroidal structure is characterized by   = 9 m,  = 3 m, meridian angle  1 =  ∈ [0, 360 ∘ ], circumferential angle  2 =  ∈ [0, 120 ∘ ], and a thickness ℎ = 0.6 m.Also the toroidal shell has CF boundary conditions, and the normal load  +  = −10000 Pa 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 where  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 ) .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  =  (1) =  (3) ,  =  (1) =  (3) ,  =  (1) =  (3) , and  =  (1) =  (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]: If the generatrix parabola needs to be characterized (the other parabola can be defined analogously), it can be described by the following equation: where   1 = ( 2 −  2 )/ is a characteristic parameter of the parabolic curve,   1 0 is the abscissa of a point of the parabola, and   1 3 is its ordinate in the generatrix plane of the parabolic curve.The abscissa   1 0 ( 1 ) of the parabolic curve assumes the form The parameters describing the two parabolas of the elliptic paraboloid under consideration are  = 3 m,  = −3 m,  = 0 m, and  = 0.8 m, and the thickness is ℎ = 0.4 m.The elliptic paraboloid is completely clamped at its four edges (CCCC), and it is subjected to a normal load  +  = −10000 Pa at the top surface.The elliptic paraboloid is a single-layered structure with FGM 2( (1) / (1) / (1) / (1) ) four-parameter powerlaw distribution.In this case, the four parameters  =  (1) ,  =  (1) ,  =  (1) , and  =  (1) have to be determined by the optimization procedure.
In Table 1 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 FGM  , the singlelayered structure is made of ceramic isotropic material, while for the second one, named FGM  , 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   = 322.27GPa, ] C = 0.24, and   = 2370 kg/m 3 , and for the Stainless Steel are   = 207.78GPa, ]  = 0.3177, and  M = 8166 kg/m 3 , respectively.As 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 function  has been modelled to obtain a defined first resonance frequency target ( 1 ) set by user, a maximum displacement target ( max  =  3 max  ) 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  is still monodimensional but keeps into account two aspects, ranging  from zero to one.According to the above considerations, the fitness expression  can be so represented by the formula It is worth to note that the minimization of the fitness  is equivalent to obtain values of frequency ( =  1 ) and of displacement ( =  max =  3 max ) as close as possible to the desired ones ( 1 ,  max  =  3 max  ), 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.Figures 5,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  and  parameter for each try; the best solution is represented by a red circle.In the upper centre there are the values of  and  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  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 Table 2: 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 , , , and .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 Mathematical Problems in Engineering 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 , , , and  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.

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 powerlaw.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.

Figure 1 :
Figure 1: Geometry description and coordinate system of a doubly curved shell.

Table 1 :
First ten frequencies and maximum static deflections for the six shell structures of Figure4made of isotropic materials: FGM : full ceramic material and FGM  : full metal material using a 31 × 31 Chebyshev-Gauss-Lobatto (C-G-L) grid distribution with different boundary conditions.