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

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.

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 𝛼₁(𝛼⁰₁≤𝛼₁≤𝛼¹₁), 𝛼₂(𝛼⁰₂≤𝛼₂≤𝛼¹₂) upon the middle surface or reference surface 𝐫(𝛼₁,𝛼₂), and 𝜁 directed along the outward normal 𝐧(𝛼₁,𝛼₂) and measured from the reference surface (−ℎ/2≤𝜁≤ℎ/2), where ℎ(𝛼₁,𝛼₂) 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: 𝐑𝛼₁,𝛼₂,𝜁=𝐫𝛼₁,𝛼₂+ℎ𝛼₁,𝛼₂2𝑧𝐧𝛼₁,𝛼₂, (1) where 𝑧=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 ℎ=^𝑙_𝑘=1ℎ_𝑘 (2) 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 𝐫(𝛼₁,𝛼₂) and the normal vector 𝐧(𝛼₁,𝛼₂) 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 ℎ(𝛼₁,𝛼₂). Hence, writing the position vector of the reference surface, it is possible to define the three components along the three global axes 𝑂𝑥₁𝑥₂𝑥₃ as 𝐫𝛼₁,𝛼₂=𝑟₁𝛼₁,𝛼₂𝐞₁+𝑟₂𝛼₁,𝛼₂𝐞₂+𝑟₃𝛼₁,𝛼₂𝐞₃, (3) where 𝐞₁,𝐞₂, and 𝐞₃ are the unit vector of the global reference system 𝑂𝑥₁𝑥₂𝑥₃. From the definition of the first fundamental form [7, 11, 12, 23] of the reference surface 𝐫(𝛼₁,𝛼₂), the Lamé parameters can be expressed as 𝐴₁𝛼₁,𝛼₂=√𝐫_,1⋅𝐫_,1,𝐴₂𝛼₁,𝛼₂=√𝐫_,2⋅𝐫_,2. (4) The comma in expressions (4) defines the partial derivative with respect to 𝛼₁ or 𝛼₂, respectively. Moreover, by considering an orthogonal curvilinear coordinate system 𝑂^′𝛼₁𝛼₂𝜁 [7, 11, 12, 23, 74] from the position vector (3) the normal vector 𝐧(𝛼₁,𝛼₂) can be written as 𝐧𝛼₁,𝛼₂=𝐫_,1×𝐫_,2𝐴₁𝐴₂. (5) Finally, due to the fact that an orthogonal curvilinear coordinate system 𝑂^′𝛼₁𝛼₂𝜁 is considered and following the definition of the second fundamental form [7, 11, 23, 74] of the reference surface 𝐫(𝛼₁,𝛼₂), the principal radii of curvature can be evaluated as 𝑅₁𝛼₁,𝛼₂𝐫=−^,1⋅𝐫_,1𝐫_,11,𝑅⋅𝐧₂𝛼₁,𝛼₂𝐫=−^,2⋅𝐫_,2𝐫_,22.⋅𝐧 (6) 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: 𝜀_𝑛=𝜀_𝑛(𝛼₁,𝛼₂,𝜁,𝑡)=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: 𝜎_𝑛=𝜎_𝑛(𝛼₁,𝛼₂,𝜁,𝑡)=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: 1ℎ100≤max𝑅_min,ℎ𝐿_min≤1.10 (7) 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) where 𝐻₁𝜁=1+𝑅₁,𝐻₂𝜁=1+𝑅₂, (9) and 𝑢₁,𝑢₂, and 𝑢₃ are the displacement components of points lying on the middle surface (𝜁=0) of the shell, while 𝑡 is the time variable. 𝛽₁ and 𝛽₂ 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 𝑢₃≪ℎ. In-plane displacements 𝑈₁ and 𝑈₂ vary linearly through the thickness, while 𝑈₃ 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: 𝜀⁰₁𝜀⁰₂𝛾⁰₁𝛾⁰₂𝜒⁰₁𝜒⁰₂𝜔⁰₁𝜔⁰₂𝛾⁰_1𝑛𝛾⁰_2𝑛=1𝐴₁𝜕𝜕𝛼₁1𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂1𝑅₁100𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁1𝐴₂𝜕𝜕𝛼₂1𝑅₂−100𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂1𝐴₁𝜕𝜕𝛼₁1000𝐴₂𝜕𝜕𝛼₂−1𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁1000𝐴₁𝑅₁𝜕𝜕𝛼₁−1𝐴₁𝑅²₁𝜕𝑅₁𝜕𝛼₁1𝐴₁𝐴₂𝑅₂𝜕𝐴₁𝜕𝛼₂01𝐴₁𝜕𝜕𝛼₁1𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂1𝐴₁𝐴₂𝑅₁𝜕𝐴₂𝜕𝛼₁1𝐴₂𝑅₂𝜕𝜕𝛼₂−1𝐴₂𝑅²₂𝜕𝑅₂𝜕𝛼₂01𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁1𝐴₂𝜕𝜕𝛼₂−1𝐴₁𝐴₂𝑅₁𝜕𝐴₁𝜕𝛼₂1𝐴₁𝑅₂𝜕𝜕𝛼₁−1𝐴₁𝑅²₂𝜕𝑅₂𝜕𝛼₁10−𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂1𝐴₁𝜕𝜕𝛼₁1𝐴₂𝑅₁𝜕𝜕𝛼₂−1𝐴₂𝑅²₁𝜕𝑅₁𝜕𝛼₂−1𝐴₁𝐴₂𝑅₂𝜕𝐴₂𝜕𝛼₁01𝐴₂𝜕𝜕𝛼₂−1𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁100𝐴₁𝜕𝜕𝛼₁11000𝐴₂𝜕𝜕𝛼₂𝑢01¹𝑢₂𝑢₃𝛽₁𝛽₂ (10)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 𝜀⁰₁,𝜀⁰₂,𝛾⁰₁, and 𝛾⁰₂ are the in-plane components, and 𝜒⁰₁,𝜒⁰₂,𝜔⁰₁, and 𝜔⁰₂ are the analogous curvature changes. The last two components 𝛾⁰_1𝑛,𝛾⁰_2𝑛 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:𝑁₁𝑁₂𝑁₁₂𝑁₂₁𝑀₁𝑀₂𝑀₁₂𝑀₂₁𝑇₁𝑇₂=𝐴₍₎⁽⁰⁾₁₁₂₀𝐴₍₎⁽⁰⁾₁₂₁₁𝐴₍₎⁽⁰⁾₁₆₂₀𝐴₍₎⁽⁰⁾₁₆₁₁𝐴₍₎⁽¹⁾₁₁₂₀𝐴₍₎⁽¹⁾₁₂₁₁𝐴₍₎⁽¹⁾₁₆₂₀𝐴₍₎⁽¹⁾₁₆₁₁𝐴00_()(0)1211𝐴₍₎⁽⁰⁾₂₂₀₂𝐴₍₎⁽⁰⁾₂₆₁₁𝐴₍₎⁽⁰⁾₂₆₀₂𝐴₍₎⁽¹⁾₁₂₁₁𝐴₍₎⁽¹⁾₂₂₀₂𝐴₍₎⁽¹⁾₂₆₁₁𝐴₍₎⁽¹⁾₂₆₀₂𝐴00_()(0)1620𝐴₍₎⁽⁰⁾₂₆₁₁𝐴₍₎⁽⁰⁾₆₆₂₀𝐴₍₎⁽⁰⁾₆₆₁₁𝐴₍₎⁽¹⁾₁₆₂₀𝐴₍₎⁽¹⁾₂₆₁₁𝐴₍₎⁽¹⁾₆₆₂₀𝐴₍₎⁽¹⁾₆₆₁₁𝐴00_()(0)1611𝐴₍₎⁽⁰⁾₂₆₀₂𝐴₍₎⁽⁰⁾₆₆₁₁𝐴₍₎⁽⁰⁾₆₆₀₂𝐴₍₎⁽¹⁾₁₆₁₁𝐴₍₎⁽¹⁾₂₆₀₂𝐴₍₎⁽¹⁾₆₆₁₁𝐴₍₎⁽¹⁾₆₆₀₂𝐴00_()(1)1120𝐴₍₎⁽¹⁾₁₂₁₁𝐴₍₎⁽¹⁾₁₆₂₀𝐴₍₎⁽¹⁾₁₆₁₁𝐴₍₎⁽²⁾₁₁₂₀𝐴₍₎⁽²⁾₁₂₁₁𝐴₍₎⁽²⁾₁₆₂₀𝐴₍₎⁽²⁾₁₆₁₁𝐴00_()(1)1211𝐴₍₎⁽¹⁾₂₂₀₂𝐴₍₎⁽¹⁾₂₆₁₁𝐴₍₎⁽¹⁾₂₆₀₂𝐴₍₎⁽²⁾₁₂₁₁𝐴₍₎⁽²⁾₂₂₀₂𝐴₍₎⁽²⁾₂₆₁₁𝐴₍₎⁽²⁾₂₆₀₂𝐴00_()(1)1620𝐴₍₎⁽¹⁾₂₆₁₁𝐴₍₎⁽¹⁾₆₆₂₀𝐴₍₎⁽¹⁾₆₆₁₁𝐴₍₎⁽²⁾₁₆₂₀𝐴₍₎⁽²⁾₂₆₁₁𝐴₍₎⁽²⁾₆₆₂₀𝐴₍₎⁽²⁾₆₆₁₁𝐴00_()(1)1611𝐴₍₎⁽¹⁾₂₆₀₂𝐴₍₎⁽¹⁾₆₆₁₁𝐴₍₎⁽¹⁾₆₆₀₂𝐴₍₎⁽²⁾₁₆₁₁𝐴₍₎⁽²⁾₂₆₀₂𝐴₍₎⁽²⁾₆₆₁₁𝐴₍₎⁽²⁾₆₆₀₂0000000000𝜅𝐴₍₎⁽⁰⁾₄₄₂₀𝜅𝐴₍₎⁽⁰⁾₄₅₁₁00000000𝜅𝐴₍₎⁽⁰⁾₄₅₁₁𝜅𝐴₍₎⁽⁰⁾₅₅₀₂𝜀⁰₁𝜀⁰₂𝛾⁰₁𝛾⁰₂𝜒⁰₁𝜒⁰₂𝜔⁰₁𝜔⁰₂𝛾⁰_1𝑛𝛾⁰_2𝑛, (11)where the elastic engineering stiffnesses 𝐴^(𝜏)_{𝑖𝑗(𝑝𝑞)} are defined as follows: 𝐴^(𝜏)_{𝑖𝑗𝑝𝑞}=^𝑙_𝑘=1^𝜁𝑘+1_𝜁𝑘𝑄^(𝑘)_𝑖𝑗𝜁^𝜏𝐻₁𝐻₂𝐻^𝑝₁𝐻^𝑞₂𝑑𝜁,𝜏,𝑝,𝑞=0,1,2. (12) Different approaches can be found in the literature to evaluate the engineering elastic constants 𝐴^(𝜏)_{𝑖𝑗(𝑝𝑞)} [7, 17–20, 73, 76]. In the present paper, the relations of the elastic engineering stiffnesses 𝐴^(𝜏)_{𝑖𝑗(𝑝𝑞)} 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 𝐴^(𝜏)_{𝑖𝑗(𝑝𝑞)} depend on curvatures, the corresponding derivatives, with respect to the coordinates along 𝛼₁ and 𝛼₂ 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 𝐴^(𝜏)_{𝑖𝑗(𝑝𝑞)} are numerically evaluated. The corresponding elastic constants 𝑄^(𝑘)_𝑖𝑗 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 𝑁₁,𝑁₂,𝑁₁₂, and 𝑁₂₁ are the in-plane force resultants, and 𝑀₁,𝑀₂,𝑀₁₂, and 𝑀₂₁ are the analogous couples, while 𝑇₁,𝑇₂ are the transverse shear force resultants. In the above definitions (11) the symmetry of shearing force resultants 𝑁₁₂,𝑁₂₁ and torsional couples 𝑀₁₂,𝑀₂₁ 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 𝜁/𝑅₁,𝜁/𝑅₂ 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 𝐸^(𝑘)(𝜁), shear modulus 𝐺^(𝑘)(𝜁), Poisson’s ratio 𝜈^(𝑘)(𝜁), and mass density 𝜌^(𝑘)(𝜁) of the functionally graded shell 𝑘th lamina can be expressed as 𝜌^(𝑘)(𝜁)=𝜌_𝐶^(𝑘)−𝜌_𝑀^(𝑘)𝑉_𝐶^(𝑘)(𝜁)+𝜌_𝑀^(𝑘),𝐸^(𝑘)(𝜁)=𝐸_𝐶^(𝑘)−𝐸_𝑀^(𝑘)𝑉_𝐶^(𝑘)(𝜁)+𝐸_𝑀^(𝑘),𝜈^(𝑘)(𝜁)=𝜈_𝐶^(𝑘)−𝜈_𝑀^(𝑘)𝑉_𝐶^(𝑘)(𝜁)+𝜈_𝑀^(𝑘),𝐺^(𝑘)(𝜁)=𝐸^(𝑘)(𝜁)21+𝑣^(𝑘)(𝜁),for𝜁_𝑘≤𝜁≤𝜁_𝑘+1, (13) 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 𝑉_𝐶^(𝑘)(𝜁) follows two simple four-parameter power-law distributions [49, 51, 53, 70, 75]: FGM_{1(𝑎(𝑘)/𝑏}^(𝑘)_{/𝑐(𝑘)/𝑝(𝑘))}∶𝑉_𝐶^(𝑘)(𝜁)=1−𝑎^(𝑘)𝜁ℎ_𝑘−𝜁_𝑘ℎ_𝑘+𝑏^(𝑘)𝜁ℎ_𝑘−𝜁_𝑘ℎ_𝑘^{𝑐(𝑘)𝑝(𝑘)},FGM_{2(𝑎(𝑘)/𝑏}^(𝑘)_{/𝑐(𝑘)/𝑝(𝑘))}∶𝑉_𝐶^(𝑘)(𝜁)=1−𝑎^(𝑘)𝜁_𝑘+1ℎ_𝑘−𝜁ℎ_𝑘+𝑏^(𝑘)𝜁_𝑘+1ℎ_𝑘−𝜁ℎ_𝑘^{𝑐(𝑘)𝑝(𝑘)}, (14) 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: 𝑉_𝐶^(𝑘)+𝑉_𝑀^(𝑘)=1. (15) 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 𝑝^(𝑘)=0⟶𝑉_𝐶^(𝑘)=1,𝑉_𝑀^(𝑘)=0⟶𝜌^(𝑘)(𝜁)=𝜌_𝐶,𝐸^(𝑘)(𝜁)=𝐸_𝐶,𝜈^(𝑘)(𝜁)=𝜈_𝐶,𝑝^(𝑘)=∞⟶𝑉_𝐶^(𝑘)=0,𝑉_𝑀^(𝑘)=1⟶𝜌^(𝑘)(𝜁)=𝜌_𝑀,𝐸^(𝑘)(𝜁)=𝐸_𝑀,𝜈^(𝑘)(𝜁)=𝜈_𝑀. (16) Some material profiles through the functionally graded shell thickness are illustrated in Figures 2 and 3.

(a)

(b)

(c)

(d)

(e)

(f)

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) FGM1(𝑎(1)=1/𝑏(1)=0/𝑐(1)/𝑝(1)), (b) FGM2(𝑎(1)=1/𝑏(1)=0/𝑐(1)/𝑝(1)), (c) FGM1(𝑎(1)=1/𝑏(1)=1/𝑐(1)=2/𝑝(1)), (d) FGM2(𝑎(1)=1/𝑏(1)=1/𝑐(1)=2/𝑝(1)), (e) FGM1(𝑎=1/𝑏=0.5/𝑐=2/𝑝), and (f) FGM2(𝑎(1)=1/𝑏(1)=0.5/𝑐(1)=2/𝑝(1)).

(a)

(b)

(c)

(d)

Figure 3 

Variations of the ceramic volume fraction 𝑉𝐶 through the thickness for two-layered and three-layered laminated shells for different values of the power-law index 𝑝=𝑝(1)=𝑝(2) and 𝑝=𝑝(1)=𝑝(3): (a) FGM1(𝑎(1)=1/𝑏(1)=0/𝑐(1)/𝑝(1))/FGM2(𝑎(2)=1/𝑏(2)=0/𝑐(2)/𝑝(2)), (b) FGM2(𝑎(1)=1/𝑏(1)=0/𝑐(1)/𝑝(1))/FGM1(𝑎(2)=1/𝑏(2)=0/𝑐(2)/𝑝(2)), (c) FGM2(𝑎(1)=1/𝑏(1)=0/𝑐(1)/𝑝(1))/FGM𝐶/FGM1(𝑎(3)=1/𝑏(3)=0/𝑐(3)/𝑝(3)), and (d) FGM2(𝑎(1)=1/𝑏(1)=0.5/𝑐(1)=2/𝑝(1))/FGM𝐶/FGM1(𝑎(3)=1/𝑏(3)=0.5/𝑐(3)=2/𝑝(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: 1𝐴₁𝐴₂𝜕𝑁₁𝐴₂𝜕𝛼₁+1𝐴₁𝐴₂𝜕𝑁₂₁𝐴₁𝜕𝛼₂+𝑁₁₂𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂−𝑁₂𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁+1𝑅₁1𝐴₁𝐴₂𝜕𝑀₁𝐴₂𝜕𝛼₁+1𝐴₁𝐴₂𝜕𝑀₂₁𝐴₁𝜕𝛼₂+𝑀₁₂𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂−𝑀₂𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁+𝑞₁=𝐼₀+2𝐼₁𝑅₁+𝐼₂𝑅²₁𝑢̈₁+𝐼₁+𝐼₂𝑅₁𝛽̈₁,1𝐴₁𝐴₂𝜕𝑁₁₂𝐴₂𝜕𝛼₁+1𝐴₁𝐴₂𝜕𝑁₂𝐴₁𝜕𝛼₂+𝑁₂₁𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁−𝑁₁𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂+1𝑅₂1𝐴₁𝐴₂𝜕𝑀₁₂𝐴₂𝜕𝛼₁+1𝐴₁𝐴₂𝜕𝑀₂𝐴₁𝜕𝛼₂+𝑀₂₁𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁−𝑀₁𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂+𝑞₂=𝐼₀+2𝐼₁𝑅₂+𝐼₂𝑅²₂𝑢̈₂+𝐼₁+𝐼₂𝑅₂𝛽̈₂,1𝐴₁𝐴₂𝜕𝑇₁𝐴₂𝜕𝛼₁+1𝐴₁𝐴₂𝜕𝑇₂𝐴₁𝜕𝛼₂−𝑁₁𝑅₁+𝑁₂𝑅₂+𝑞_𝑛=𝐼₀𝑢̈₃,1𝐴₁𝐴₂𝜕𝑀₁𝐴₂𝜕𝛼₁+1𝐴₁𝐴₂𝜕𝑀₂₁𝐴₁𝜕𝛼₂+𝑀₁₂𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂−𝑀₂𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁−𝑇₁+𝑚₁=𝐼₁+𝐼₂𝑅₁𝑢̈₁+𝐼₂𝛽̈₁,1𝐴₁𝐴₂𝜕𝑀₁₂𝐴₂𝜕𝛼₁+1𝐴₁𝐴₂𝜕𝑀₂𝐴₁𝜕𝛼₂+𝑀₂₁𝐴₁𝐴₂𝜕𝐴₂𝜕𝛼₁−𝑀₁𝐴₁𝐴₂𝜕𝐴₁𝜕𝛼₂−𝑇₂+𝑚₂=𝐼₁+𝐼₂𝑅₂𝑢̈₂+𝐼₂𝛽̈₂, (17) where 𝐼_𝜏=^𝑙_𝑘=1^𝜁𝑘+1_𝜁𝑘𝜌^(𝑘)𝜁^𝜏𝐻₁𝐻₂𝑑𝜁,𝜏=0,1,2. (18) Furthermore, the generalized external actions 𝑞₁,𝑞₂,𝑞_𝑛,𝑚₁, and 𝑚₂ 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: 𝑞₁=𝑞⁺₁ℎ1+2𝑅₁²ℎ1+2𝑅₂+𝑞⁻₁ℎ1−2𝑅₁²ℎ1−2𝑅₂,𝑞₂=𝑞⁺₂ℎ1+2𝑅₁ℎ1+2𝑅₂²+𝑞⁻₂ℎ1−2𝑅₁ℎ1−2𝑅₂²,𝑞_𝑛=𝑞⁺_𝑛ℎ1+2𝑅₁ℎ1+2𝑅₂+𝑞⁻_𝑛ℎ1−2𝑅₁ℎ1−2𝑅₂,𝑚₁=𝑞⁺₁ℎ2ℎ1+2𝑅₁ℎ1+2𝑅₂−𝑞⁻₁ℎ2ℎ1−2𝑅₁ℎ1−2𝑅₂,𝑚₂=𝑞⁺₂ℎ2ℎ1+2𝑅₁ℎ1+2𝑅₂−𝑞⁻₂ℎ2ℎ1−2𝑅₁ℎ1−2𝑅₂, (19) where 𝑞⁺₁,𝑞⁻₁,𝑞⁺₂,𝑞⁻₂,𝑞⁺_𝑛, and 𝑞⁻_𝑛 are the external forces in the three principal directions 𝛼₁,𝛼₂, 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𝐿₁₁𝐿₁₂𝐿₁₃𝐿₁₄𝐿₁₅𝐿₂₁𝐿₂₂𝐿₂₃𝐿₂₄𝐿₂₅𝐿₃₁𝐿₃₂𝐿₃₃𝐿₃₄𝐿₃₅𝐿₄₁𝐿₄₂𝐿₄₃𝐿₄₄𝐿₄₅𝐿₅₁𝐿₅₂𝐿₅₃𝐿₅₄𝐿₅₅𝑢₁𝑢₂𝑢₃𝛽₁𝛽₂+𝑞₁𝑞₂𝑞_𝑛𝑚₁𝑚₂=𝐼₀₁00𝐼₁₁00𝐼₀₂00𝐼₁₂00𝐼₀𝐼00₁₁00𝐼₂00𝐼₁₂00𝐼₂𝑢̈₁𝑢̈₂𝑢̈₃𝛽̈₁𝛽̈₂, (20)where 𝐿_𝑖𝑗,𝑖,𝑗=1,…,5 are the equilibrium operators and the new mass inertias are defined as follows: 𝐼₀₁=𝐼₀+2𝐼₁𝑅₁+𝐼₂𝑅²₁,𝐼₁₁=𝐼₁+𝐼₂𝑅₁,𝐼₀₂=𝐼₀+2𝐼₁𝑅₂+𝐼₂𝑅²₂,𝐼₁₂=𝐼₁+𝐼₂𝑅₂. (21) 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):𝑢₁=𝑢₂=𝑢₃=𝛽₁=𝛽₂=0at𝛼₁=𝛼⁰₁or𝛼₁=𝛼¹₁,𝛼⁰₂≤𝛼₂≤𝛼¹₂,𝑢₁=𝑢₂=𝑢₃=𝛽₁=𝛽₂=0at𝛼₂=𝛼⁰₂or𝛼₂=𝛼¹₂,𝛼⁰₁≤𝛼₁≤𝛼¹₁, (22)

Simply supported edge boundary conditions (S):𝑢₁=𝑢₂=𝑢₃=𝛽₂=0,𝑀₁=0at𝛼₁=𝛼⁰₁or𝛼₁=𝛼¹₁,𝛼⁰₂≤𝛼₂≤𝛼¹₂,𝑢₁=𝑢₂=𝑢₃=𝛽₂=0,𝑀₂=0at𝛼₂=𝛼⁰₂or𝛼₂=𝛼¹₂,𝛼⁰₁≤𝛼₁≤𝛼¹₁, (23)

Free edge boundary conditions (F):𝑁₁+𝑀₁𝑅₁=0,𝑁₁₂+𝑀₁₂𝑅₂𝑇=0,₁=0,𝑀₁=𝑀₁₂=0at𝛼₁=𝛼⁰₁or𝛼₁=𝛼¹₁,𝛼⁰₂≤𝛼₂≤𝛼¹₂,𝑁₂+𝑀₂𝑅₂=0,𝑁₂₁+𝑀₂₁𝑅₁𝑇=0,₂=0,𝑀₂=𝑀₂₁=0at𝛼₂=𝛼⁰₂or𝛼₂=𝛼¹₂,𝛼⁰₁≤𝛼₁≤𝛼¹₁. (24)(25) In addition to the external boundary conditions (22)–(25), the kinematic and physical compatibility conditions should be satisfied at the common closing meridians with 𝛼₂=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 𝛼¹₂=2𝜋, it is necessary to implement the kinematic and physical compatibility conditions between the two computational meridians with 𝛼⁰₂=0 and with 𝛼¹₂=2𝜋:

Kinematic compatibility conditions along the closing meridian (𝛼₂=0,2𝜋)𝑢₁𝛼₁,0,𝑡=𝑢₁𝛼₁,2𝜋,𝑡,𝑢₂𝛼₁,0,𝑡=𝑢₂𝛼₁,𝑢,2𝜋,𝑡₃𝛼₁,0,𝑡=𝑢₃𝛼₁,𝛽,2𝜋,𝑡₁𝛼₁,0,𝑡=𝛽₁𝛼₁,2𝜋,𝑡,𝛽₂𝛼₁,0,𝑡=𝛽₂𝛼₁𝛼,2𝜋,𝑡₀₁≤𝛼₁≤𝛼¹₁, (26)Physical compatibility conditions along the closing meridian (𝛼₂=0,2𝜋)𝑁₂𝛼₁+𝑀,0,𝑡²𝛼₁,0,𝑡𝑅₂=𝑁₂𝛼₁+𝑀,2𝜋,𝑡²𝛼₁,2𝜋,𝑡𝑅₂,𝑁₂₁𝛼₁+𝑀,0,𝑡²¹𝛼₁,0,𝑡𝑅₁=𝑁₂₁𝛼₁+𝑀,2𝜋,𝑡²¹𝛼₁,2𝜋,𝑡𝑅₁,𝑇₂𝛼₁,0,𝑡=𝑇₂𝛼₁,𝑀,2𝜋,𝑡₂𝛼₁,0,𝑡=𝑀₂𝛼₁,𝑀,2𝜋,𝑡₂₁𝛼₁,0,𝑡=𝑀₂₁𝛼₁,𝛼,2𝜋,𝑡₀₁≤𝛼₁≤𝛼¹₁. (27) 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 and with 𝛼¹₁=2𝜋:

Kinematic compatibility conditions along the closing parallel (𝛼₁=0,2𝜋)𝑢₁0,𝛼₂,𝑡=𝑢₁2𝜋,𝛼₂,𝑡,𝑢₂0,𝛼₂,𝑡=𝑢₂2𝜋,𝛼₂,𝑢,𝑡₃0,𝛼₂,𝑡=𝑢₃2𝜋,𝛼₂,𝛽,𝑡₁0,𝛼₂,𝑡=𝛽₁2𝜋,𝛼₂,𝑡,𝛽₂0,𝛼₂,𝑡=𝛽₂2𝜋,𝛼₂,𝛼,𝑡₀₂≤𝛼₂≤𝛼¹₂, (28)

Physical compatibility conditions along the closing parallel (𝛼₁=0,2𝜋)𝑁₁0,𝛼₂+𝑀,𝑡¹0,𝛼₂,𝑡𝑅₁=𝑁₁2𝜋,𝛼₂+𝑀,𝑡¹2𝜋,𝛼₂,𝑡𝑅₁,𝑁₁₂0,𝛼₂+𝑀,𝑡¹²0,𝛼₂,𝑡𝑅₂=𝑁₁₂2𝜋,𝛼₂+𝑀,𝑡¹²2𝜋,𝛼₂,𝑡𝑅₂,𝑇₁0,𝛼₂,𝑡=𝑇₁2𝜋,𝛼₂,𝑀,𝑡₁0,𝛼₂,𝑡=𝑀₁2𝜋,𝛼₂,𝑀,𝑡₁₂0,𝛼₂,𝑡=𝑀₁₂2𝜋,𝛼₂𝛼,𝑡₀₂≤𝛼₂≤𝛼¹₂. (29)

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 (𝛼_1𝑖,𝛼_2𝑗) along the reference surface are in the discrete form: 𝛼_1𝑖=1−cos𝑖−1𝜋𝛼𝑁−1¹¹−𝛼⁰₁2+𝛼⁰₁,𝑖=1,2,…,𝑁,for𝛼₁∈𝛼⁰₁,𝛼¹₁,𝛼_2𝑗=1−cos𝑗−1𝜋𝛼𝑀−1¹²−𝛼⁰₂2+𝛼⁰₂,𝑗=1,2,…,𝑀,for𝛼₂∈𝛼⁰₂,𝛼¹₂, (30) where 𝑁,𝑀 are the total number of sampling points used to discretize the domain in 𝛼₁ and 𝛼₂ 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]: 𝜕^𝑛𝑓(𝑥)𝜕𝑥^𝑛_{𝑥=𝑥𝑚}=^𝑇_𝑘=1𝜍^(𝑛)_𝑚𝑘𝑓𝑥_𝑘,𝑚=1,2,…,𝑇. (31) 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) 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: 𝐊_𝑑𝑑−𝐊_𝑑𝑏𝐊⁻¹_𝑏𝑏𝐊_𝑏𝑑𝜹_𝑑=𝐟_𝑑−𝐊_𝑑𝑏𝐊⁻¹_𝑏𝑏𝐟_𝑏. (33) 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 𝑞⁺₁,𝑞⁻₁,𝑞⁺₂,𝑞⁻₂,𝑞⁺_𝑛, 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: 𝑢₁𝛼₁,𝛼₂,𝑡=𝑈¹𝛼₁,𝛼₂𝑒^𝑖𝜔𝑡,𝑢₂𝛼₁,𝛼₂,𝑡=𝑈²𝛼₁,𝛼₂𝑒^𝑖𝜔𝑡,𝑢₃𝛼₁,𝛼₂,𝑡=𝑈³𝛼₁,𝛼₂𝑒^𝑖𝜔𝑡,𝛽₁𝛼₁,𝛼₂,𝑡=𝐵¹𝛼₁,𝛼₂𝑒^𝑖𝜔𝑡,𝛽₂𝛼₁,𝛼₂,𝑡=𝐵²𝛼₁,𝛼₂𝑒^𝑖𝜔𝑡, (34) where the vibration spatial amplitude values 𝑈¹,𝑈²,𝑈³,𝐵¹, and 𝐵² 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) 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.

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<𝑉_𝐶^(𝑘)<1). The pseudocode of the GA implemented is shown in the following Pseudocode 1.

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

Pseudocode 1 

Genetic Algorithm pseudocode.

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) can be evaluated by the generation of random samples 𝑥₁,𝑥₂,𝑥₃,𝑥₄,…,𝑥_𝑆 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 1𝐼≈𝑆^𝑆_𝑠=1𝑔𝑥_𝑠. (37) 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.

Step (1) Set a max and min value for each parameter Step (2) Set the max number of population (𝑛𝑔) members, the convergence tolerance (𝜀) and    the number of consecutive no improvements loops (𝑛𝜀) Step (3) while (𝑖<𝑛𝑔) or (𝜀 for 𝑛𝜀 loops false) Step (4)    Randomly generate a member of population 𝑃𝑖, following a probability       distribution (normal, exponential, Weibull [113]). Step (5)    Compute volume fraction distribution (𝑉𝐶) in the thickness for 𝑃𝑖 Step (6)    if 0<𝑉𝐶<1 Step (7)     𝑖=𝑖+1 Step (8)     Compute Fitness Step (9)     if Fitness (𝑖) > Best Fitness Step (10)     Best Fitness = Fitness (𝑖) Step (11)     Best member = 𝑖 Step (12)    end if Step (12)    end if Step (13) end while

Pseudocode 2 

Monte Carlo pseudocode.

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 𝑛-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 𝑣_𝑖^(𝑚)=𝜙^(𝑚)𝑣_𝑖^(𝑚−1)+𝑐₁𝑎₁𝑃_𝑏𝑖−𝑃_𝑖^(𝑚−1)+𝑐₂𝑎₂𝑃_𝑏𝐺𝑏−𝑃_𝑖^(𝑚−1), (38) where 𝑚 is the algorithm step (𝑛_𝐺 generations); 𝑖 is the index of parameter of the single particle; 𝜙^(𝑚) is the inertia function; 𝑣_𝑖^(𝑚) 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); 𝑎₁,𝑎₂ are the acceleration constants; 𝑐₁,𝑐₂ are random numbers in the interval [0,1].

The following formulation can be applied to obtain the new position: 𝑃_𝑖^(𝑚)=𝑃_𝑖^(𝑚−1)+𝑣_𝑖^(𝑚), (39) where 𝑃_𝑖^(𝑚) is the position of the 𝑖th particle at the 𝑚-step and 𝑣_𝑖^(𝑚) 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.

Step (1) Set number of particles (𝑛𝑃), number of parameters for each particle (𝑛param), max    number of generations (𝑛𝐺), convergence tolerance (𝜀) and the number of    consecutive no improvements loops (𝑛𝜀), 𝛿 (ratio between position and speed) Step (2) Set particles boundary values: max1∶𝑛param, min1∶𝑛param Step (3) Randomly generate particles position 𝑃𝑖(0)∣min(𝑖)<𝑃𝑖(0)<max(𝑖) and velocity 𝑣𝑖(0) Step (4) Evaluate fitness function 𝐹𝑖=𝑓𝑃𝑖(0),𝑖=1,…,𝑛𝑃 Step (5) Set 𝑃𝑏=𝑃(0), 𝐹best=𝐹, 𝐺𝑏=min𝐹𝑖, 𝑔𝑏=indexofthemin𝐹𝑖value Step (6) for 𝑚=1 to 𝑛𝐺 Step (7)   for 𝑖=1 to 𝑛𝑃 Step (8)     Δ𝑣𝑖=𝑐1𝑎1𝑃𝑏𝑖−𝑃𝑖(𝑚−1)+𝑐2𝑎2(𝑃𝑏𝐺𝑏−𝑃𝑖(𝑚−1)) Step (9)     𝑣𝑖(𝑚)=𝜙(𝑚)𝑣(𝑚−1)+Δ𝑣𝑖 Step (10)   𝑃𝑖(𝑚)=𝑃𝑖(𝑚−1)+𝑣𝑖(𝑚) Step (11)    Compute Volume fraction distribution (𝑉𝐶) in the thickness for 𝑃𝑖(𝑚) Step (12)     while 𝑃𝑖(𝑚)>max(𝑖)or𝑃𝑖(𝑚)<min(𝑖)or𝑉𝐶<0or𝑉𝐶>1 Step (13)      𝑣𝑖(𝑚)=𝑣(𝑚−1)+𝛿⋅(max(𝑖)−min(𝑖))⋅random(0÷1) Step (14)      𝑃𝑖(𝑚)=𝑃𝑖(𝑚−1)+𝑣𝑖(𝑚) Step (15)      end while Step (16)       Evaluate fitness 𝐹𝑖=𝑓𝑃𝑖(𝑚),𝑖=1,…,𝑛𝑃 Step (17)       for 𝑖=1,…,𝑛𝑃 Step (18)      if 𝐹𝑖<𝐹𝑏𝑖 then 𝑃𝑏𝑖=𝑃𝑖(𝑚) and 𝐹𝑏𝑖=𝐹𝑖 Step (19)      if 𝐹𝑖<𝐺𝑏 then 𝐺𝑏=𝐹𝑖 and 𝑔𝑏=𝑖 Step (20)       end for Step (21)  end for Step (22) end for

Pseudocode 3 

Particle Swarm Optimization pseudocode.

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 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–53, 68–76] for shell and panel structures.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

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

For a rectangular flat plate, the position vector [74, 78] can be written as 𝐫𝛼₁,𝛼₂=−𝛼₂𝐞₂+𝛼₁𝐞₃. (40) In Figure 4(a) the considered square plate has the sides 𝑎=𝑏=1m, thickness ℎ=0.1m, and it is subjected to a normal load 𝑞⁺_𝑛=−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 FGM_{1(𝑎(1)/𝑏}⁽¹⁾_{/𝑐(1)/𝑝(1))} four-parameter power-law distribution. In this case, the four parameters 𝑎=𝑎⁽¹⁾,𝑏=𝑏⁽¹⁾,𝑐=𝑐⁽¹⁾, and 𝑝=𝑝⁽¹⁾ 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 for the sake of completeness. The position vector for the conical shell can be written as follows: 𝐫𝛼₁,𝛼₂=𝑅_𝑏+𝛼₁sin𝛼cos𝛼₂𝐞₁−𝑅_𝑏+𝛼₁sin𝛼sin𝛼₂𝐞₂+𝛼₁cos𝛼𝐞₃. (41) 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 𝐫𝛼₁,𝛼₂=𝑅_𝑏cos𝛼₂𝐞₁−𝑅_𝑏sin𝛼₂𝐞₂+𝛼₁𝐞₃. (42) The meridian abscissa 𝛼₁ is defined as 𝛼₁∈[0,𝐿], where 𝐿 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 𝑅=𝑅_𝑏=2m, 𝐿=3m, circumferential angle 𝛼₂=𝜗∈[0,120^∘], and a thickness ℎ=0.2m. It is subjected to a normal load 𝑞⁺_𝑛=−10000Pa at the top surface, and its boundary conditions are SFSF. On the contrary, the conical shell (Figure 4(c)) has 𝑅_𝑏=1m, 𝐿=3m, circumferential angle 𝛼₂=𝜗∈[0,360^∘], 𝛼=30^∘, and a thickness ℎ=0.3m. Also the conical shell is subjected to a normal load 𝑞⁺_𝑛=−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 FGM_{2(𝑎(1)/𝑏}⁽¹⁾_{/𝑐(1)/𝑝(1))}/FGM_{1(𝑎(2)/𝑏}⁽²⁾_{/𝑐(2)/𝑝(2))} lamination scheme, while the conical shell presents FGM_{1(𝑎(1)/𝑏}⁽¹⁾_{/𝑐(1)/𝑝(1))}/FGM_{2(𝑎(2)/𝑏}⁽²⁾_{/𝑐(2)/𝑝(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 𝑎=𝑎⁽¹⁾=𝑎⁽²⁾,𝑏=𝑏⁽¹⁾=𝑏⁽²⁾,𝑐=𝑐⁽¹⁾=𝑐⁽²⁾, and 𝑝=𝑝⁽¹⁾=𝑝⁽²⁾ 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 𝐫𝛼₁,𝛼₂=𝑅_𝑏+𝑅sin𝛼₁cos𝛼₂𝐞₁−𝑅_𝑏+𝑅sin𝛼₁sin𝛼₂𝐞₂+𝑅cos𝛼₁𝐞₃, (43) 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 𝑅_𝑏=9m, 𝑅=3m, meridian angle 𝛼₁=𝜑∈[0,360^∘], circumferential angle 𝛼₂=𝜗∈[0,120^∘], and a thickness ℎ=0.6m. Also the toroidal shell has CF boundary conditions, and the normal load 𝑞⁺_𝑛=−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 𝐫𝛼₁,𝛼₂=1𝐴cosharcsinhtan𝛼₁+𝑅_𝑏cos𝛼₂𝐞₁−1𝐴cosharcsinhtan𝛼₁+𝑅_𝑏sin𝛼₂𝐞₂1+𝐴arcsinhtan𝛼₁𝐞₃, (44) 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(e)) its geometrical properties are 𝐴=2m, 𝐵=2m, 𝐷=1m, 𝑅_𝑏=0m, circumferential angle 𝛼₂=𝜗∈[0,90^∘], and a thickness ℎ=0.5m. The catenoidal panel has CFCF boundary conditions, and it is subjected to a normal load 𝑞⁺_𝑛=−10000Pa at the top surface. The toroidal shell panel and catenary panel of revolution have the same three-layered lamination scheme FGM_{2(𝑎(1)/𝑏}⁽¹⁾_{/𝑐(1)/𝑝(1))}/FGM_𝐶/FGM_{1(𝑎(3)/𝑏}⁽³⁾_{/𝑐(3)/𝑝(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 𝑎=𝑎⁽¹⁾=𝑎⁽³⁾,𝑏=𝑏⁽¹⁾=𝑏⁽³⁾,𝑐=𝑐⁽¹⁾=𝑐⁽³⁾, and 𝑝=𝑝⁽¹⁾=𝑝⁽³⁾ 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]: 𝐫𝛼₁,𝛼₂=𝑅^𝛼1₀𝛼₁−𝑥^𝛼2₃𝛼₂sin𝛼₁cos𝛼₂𝐞₁−𝑅^𝛼2₀𝛼₂𝐞₂+𝑥^𝛼1₃𝛼₁+𝑥^𝛼2₃𝛼₂cos𝛼₁𝐞₃. (45) If the generatrix parabola needs to be characterized (the other parabola can be defined analogously), it can be described by the following equation: 𝑅^𝛼1₀²−𝜅^𝛼1𝑥^𝛼1₃=0, (46) where 𝜅^𝛼1=(𝐴²−𝐷²)/𝐵 is a characteristic parameter of the parabolic curve, 𝑅^𝛼1₀ is the abscissa of a point of the parabola, and 𝑥^𝛼1₃ is its ordinate in the generatrix plane of the parabolic curve. The abscissa 𝑅^𝛼1₀(𝛼₁) of the parabolic curve assumes the form 𝑅^𝛼1₀𝛼₁=𝜅^𝛼1tan𝛼₁2. (47) The parameters describing the two parabolas of the elliptic paraboloid under consideration are 𝐴=3m, 𝐶=−3m, 𝐷=0m, and 𝐵=0.8m, and the thickness is ℎ=0.4m. The elliptic paraboloid is completely clamped at its four edges (CCCC), and it is subjected to a normal load 𝑞⁺_𝑛=−10000Pa at the top surface. The elliptic paraboloid is a single-layered structure with FGM_{2(𝑎(1)/𝑏}⁽¹⁾_{/𝑐(1)/𝑝(1))} four-parameter power-law distribution. In this case, the four parameters 𝑎=𝑎⁽¹⁾,𝑏=𝑏⁽¹⁾,𝑐=𝑐⁽¹⁾, and 𝑝=𝑝⁽¹⁾ 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 single-layered 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 𝜌_𝐶=2370kg/m³, and for the Stainless Steel are 𝐸_𝑀=207.78GPa, 𝜈_𝑀=0.3177, and 𝜌_M=8166kg/m³, respectively.

Table 1 

First ten frequencies and maximum static deflections for the six shell structures of Figure 4 made 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.

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𝑇=𝑢_3max𝑇) 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 𝐹=𝜂𝑓₁−𝑓_1𝑇𝑓_1𝑇²+1−𝜂𝑤_max−𝑤_max𝑇𝑤_max𝑇². (48) It is worth to note that the minimization of the fitness 𝐹 is equivalent to obtain values of frequency (𝑓=𝑓₁) and of displacement (𝑤=𝑤_max=𝑢_3max) as close as possible to the desired ones (𝑓_1𝑇,𝑤_max𝑇=𝑢_3max𝑇), 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 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 𝑓=265.364Hz.

Figure 6 

Square plate: Genetic Algorithm optimization with 𝜂=1; final optimized frequency 𝑓=268.816Hz.

Figure 7 

Square plate: Particle Swarm Optimization with 𝜂=1; final optimized frequency 𝑓=268.710Hz.

Figure 8 

Cylindrical panel: Monte Carlo optimization with 𝜂=1; final optimized displacement 𝑤=1.1596⋅10−5m.

Figure 9 

Cylindrical panel: Genetic Algorithm optimization with 𝜂=1; final optimized displacement 𝑤=1.1668⋅10−5m.

Figure 10 

Cylindrical panel: Particle Swarm Optimization with 𝜂=1; final optimized displacement 𝑤=1.1653⋅10−5m.

Figure 11 

Conical shell: Monte Carlo optimization with 𝜂=0.5; final optimized frequency 𝑓=330.946Hz and final optimized displacement 𝑤=5.2073⋅10−7m.

Figure 12 

Conical shell: Genetic Algorithm optimization with 𝜂=0.5; final optimized frequency 𝑓=321.890Hz and final optimized displacement 𝑤=5.2624⋅10−7m.

Figure 13 

Conical shell: optimization by Particle Swarm Algorithm with 𝜂=0.5; final optimized frequency 𝑓=322.049Hz and final optimized displacement 𝑤=5.2615⋅10−7m.

Figure 14 

Toroidal shell panel: Monte Carlo optimization with 𝜂=0; final optimized displacement 𝑤=3.9727⋅10−5m.

Figure 15 

Toroidal shell panel: Genetic Algorithm optimization with 𝜂=0; final optimized displacement 𝑤=3.9495⋅10−5m.

Figure 16 

Toroidal shell panel: Particle Swarm Optimization with 𝜂=0; final optimized displacement 𝑤=3.9519⋅10−5m.

Figure 17 

Catenoidal panel: Monte Carlo optimization with 𝜂=1; final optimized frequency 𝑓=595.868Hz.

Figure 18 

Catenoidal panel: Genetic Algorithm optimization with 𝜂=1; final optimized frequency 𝑓=590.800Hz.

Figure 19 

Catenoidal panel: Particle Swarm optimization with 𝜂=1; final optimized frequency 𝑓=590.119Hz.

Figure 20 

Elliptic paraboloid: Monte Carlo optimization with 𝜂=0.5; final optimized frequency 𝑓=332.688Hz and final optimized displacement 𝑤=1.8923⋅10−6m.

Figure 21 

Elliptic paraboloid: Genetic Algorithm optimization with 𝜂=0.5; final optimized frequency 𝑓=327.600Hz and final optimized displacement 𝑤=1.9075⋅10−6m.

Figure 22 

Elliptic paraboloid: Particle Swarm Optimization with 𝜂=0.5; final optimized frequency 𝑓=327.398Hz and final optimized displacement 𝑤=1.9065⋅10−6m.

Figures 5–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 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 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.

Table 2 

Optimization of the six different structures of Figure 4: (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 𝑓₁; 𝜂=0: optimization of maximum bending displacement 𝑤_max=𝑢_3max; 𝜂=0.5: optimization of the first frequency 𝑓₁ and maximum bending displacement 𝑤_max=𝑢_3max.

Table 3 

Perceptual relative errors in first frequtency 𝑒_𝑓=((𝑓₁−𝑓_1𝑇)/𝑓_1𝑇)⋅100 and in maximum bending displacement 𝑒_𝑤=((𝑤_max−𝑤_max𝑇)/𝑤_max𝑇)⋅100 with Monte Carlo, Particle Swarm Optimization, and Genetic Algorithm methods.

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

Acknowledgments

This research was supported by the Italian Ministry Education, Universities and Research (MIUR).