Numerical Study on the Buckling Behavior of FG Porous Spherical Caps Reinforced by Graphene Platelets

The buckling response of functionally graded (FG) porous spherical caps reinforced by graphene platelets (GPLs) is assessed here, including both symmetric and uniform porosity patterns in the metal matrix, together with five different GPL distributions. The Halpin–Tsai model is here applied, together with an extended rule of mixture to determine the elastic properties and mass density of the selected shells, respectively. The equilibrium equations of the pre-buckling state are here determined according to a linear three-dimensional (3D) elasticity basics and principle of virtual work, whose solution is determined from classical finite elements. The buckling load is, thus, obtained based on the nonlinear Green strain field and generalized geometric stiffness concept. A large parametric investigation studies the sensitivity of the natural frequencies of FG porous spherical caps reinforced by GPLs to different parameters, namely, the porosity coefficients and distributions, together with different polar angles and stiffness coefficients of the elastic foundation, but also different GPL patterns and weight fractions of graphene nanofillers. Results denote that the maximum and minimum buckling loads are reached for GPL-X and GPL-O distributions, respectively. Additionally, the difference between the maximum and minimum critical buckling loads for different porosity distributions is approximately equal to 90%, which belong to symmetric distributions. It is also found that a high weight fraction of GPLs and a high porosity coefficient yield the highest and lowest effects of the structure on the buckling loads of the structure for an amount of 100% and 12.5%, respectively.


Introduction
Nowadays, there is a high demand for materials with a low weight and high strength, for many engineering applications. Among them, FG-GPL porous materials have attracted the interest of many researchers due to their mechanical potentials in aerospace and marine industries. A large variety of works from the scientific literature have focused on the static and/or dynamic behavior of different structural members, such as beams, plates, shells, with arbitrary shapes and made of composite materials [1][2][3][4]. For example, Zhang et al. [5] applied the DSC-regularized Dirac-delta method using the Timoshenko theory to explore the dynamics of FG-GPL porous beams resting on elastic foundations and subjected to a moving load. Based on a shear and normal deformation theory and by employing the Ritz approach, Priyanka et al. [6] investigated the stability and dynamic responses of porous beams made of FG-GPLs. Moreover, the free vibrations of rotating, FG-GPL, porous Timoshenko beams were studied by Binh et al. [7], using the generalized differential quadrature method (GDQM). Xu et al. [8] adopted the differential transformation method to investigate the free vibration behavior of FG-GPL porous beams based on the on the conveying fluid flow. Among the recent literature, Zhou et al. [32] examined the flutter and vibration properties of FG-GPL, porous cylindrical panels under a supersonic flow. At the same time, the vibration of FG-GPL porous shells was analytically investigated by Ebrahimi et al. [33]. Pourjabari et al. [34] analytically investigated the effect of porosity on the free and forced-vibration characteristics of GPL-reinforcement composite cylindrical shells in a nonlocal sense, based on a modified strain gradient theory (MSGT). In line with the previous works, a limited attention has been paid to the buckling response of FG-GPL porous materials and structures. Among the available literature, Zhou et al. [35] presented an accurate nonlinear buckling study of FG-GPL, porous, composite cylindrical shells based on Donnell's theory and HSDT. Shahgholian-Ghahfarokhi et al. [36,37] investigated the torsional buckling behavior of FG-GPL, porous cylindrical shells, according to a FSDT and Rayleigh-Ritz method. Similarly, Yang [38] applied the Chebyshev polynomials-based Ritz method to study the natural frequencies and buckling response of FG-GPL porous rectangular plates, using the FSDT approach. Dong [39] investigated the buckling behavior of spinning cylindrical shells made of FG-GPL porous nanocomposites, while applying a FSDT and Galerkin approach. A novel numerical DQ-FEM solution to investigating the buckling and post-buckling of FG-GPL porous plates with different shapes and boundary conditions was applied by Ansari et al. [40]. Kitipornchai [41] analyzed the natural frequencies and elastic buckling of FG-GPL porous beams using the Timoshenko beam approach and the Ritz method. Twinkle et al. [42] focused on the impacts of grading, porosity and edge loads on the natural frequency and buckling problems of porous cylindrical panels made of FG-GPLs. Nguyen [43] investigated the buckling, instability and natural-frequency response of FG porous plates reinforced by GPLs using three-variable higher order isogeometric analysis (IGA). Rafiei Anamagh and Bediz [44], instead, applied the FSDT to study the free vibration and buckling behavior of porous plates made of FG-GPLs using a spectral Chebyshev approach.
In the available literature, it seems that the static, buckling and dynamic behavior of porous spherical shells made of FG-GPLs has not been surveyed so far, despite their geometry being of great interest in various engineering applications, such as heat exchangers or energy absorbers, among other applications in the areas of aerospace, mechanical engineering and marine engineering. Among the different shell geometries, a spherical shell structures, indeed, features a high strength with a simple geometry, even compared to a cylindrical structure. The design of such structural members considering only static loading conditions may fail in dynamic situations. In such context, we focus on the buckling capacities of spherical shells made of porous FG nanocomposites reinforced by graphene, due to their exceptional flexibility and enhanced physical features. It is well known from the literature, indeed, that porous ceramic nanocomposites can ensure different beneficial effects, such as a reduced electrical and thermal conductivity; low weight; reasonable hardness; and resistance to wear, corrosion and high-temperature applications [45]. Among the few works on spherical shell dynamics available in the literature, we cite Refs. [46,47], where a Ritz-Galerking procedure was proposed to solve a dynamic buckling problem for clumped spherical members. A finite difference method was applied, instead, in [48][49][50], to check for the sensitivity of the dynamic buckling response of spherical caps to some initial manufacturing imperfections. Novel theoretical shear deformation theories were applied in Refs. [51,52] to treat the buckling response of isotropic and orthotropic shallow spherical caps, whose problem was solved analytically by means of Chebychev series [51], or numerically according to classical finite elements [52]. At the present state, however, there is a general lack of works from the literature focusing on the dynamic buckling of GPL-reinforced porous nanocomposite spherical shells, whose aspects are explored here according to the 3D elasticity basics and Green deformation nonlinearities, rather than common shell theories and Von-Karman nonlinearities, as proposed in [53].
The equilibrium equations of a pre-buckling state are determined from the principle of virtual work, whose solution is found according to classical finite elements. The buckling loads are obtained according to the nonlinear Green strain field and the generalized geo-metric stiffness concept, for spherical caps featuring a uniform and non-uniform pattern of GPLs in the metallic matrix, including open-cell internal pores and for various porosity distributions along the shell's thickness with uniform and symmetric FG patterns. More specifically, five different patterns of GPL dispersion pattern are assumed throughout the shell's thickness, namely, a FG GPL-X, A, V, UD and O patterns. A systematic investigation checks for the effects of various porosity distributions and GPL patterns, along with the weight fractions and porosity coefficients of nano-fillers and different polar angles, on the buckling behavior of FG-GPL, porous spherical caps.

Description of Geometry and Mechanical Properties
Let assume a spherical cap with uniform thickness h and mean radius of R. The outer r out = b and inner r in = a radii of the spherical cap are denoted as reported in Figure 1. The spherical cap is defined using the spherical coordinates that can be defined as follows: of virtual work, whose solution is found according to classical finite elemen buckling loads are obtained according to the nonlinear Green strain field a generalized geometric stiffness concept, for spherical caps featuring a uniform a uniform pattern of GPLs in the metallic matrix, including open-cell internal pores various porosity distributions along the shell's thickness with uniform and symm patterns. More specifically, five different patterns of GPL dispersion pattern are a throughout the shell's thickness, namely, a FG GPL-X, A, V, UD and O patt systematic investigation checks for the effects of various porosity distributions a patterns, along with the weight fractions and porosity coefficients of nano-fill different polar angles, on the buckling behavior of FG-GPL, porous spherical caps

Description of Geometry and Mechanical Properties
Let assume a spherical cap with uniform thickness h and mean radius of outer out r b = and inner in r a = radii of the spherical cap are denoted as repo  As can be seen in Figure 2, two types of non-uniform symmetric distribution uniform one are assumed in the present work, such that three different porosity  As can be seen in Figure 2, two types of non-uniform symmetric distributions and a uniform one are assumed in the present work, such that three different porosity profiles are here considered throughout the thickness of the spherical cap. In distribution 1, the porosity is nonlinear and symmetric. Furthermore, the distribution around the mid-radius is larger than the corresponding one around the external surfaces of the structure. In distribution 2, the porosity is also nonlinear and symmetric, but the porosity near the inner and outer surfaces of the spherical cap is higher than that one around the mid-radius. Equations (1)-(3) define mathematically the distributions of the material properties considering the effect of porosity, for the three selected distributions. At the same time, Figure 2 reports the five GPL distribution profiles throughout the spherical cap, thickness-wise, which are defined next [25,26]. More specifically, the mechanical properties refer to the mass density ρ(r), Young's modulus E(r) and shear modulus G(r) of porous nanocomposite spherical caps [54][55][56][57][58].
where ρ * , E * and G * refer to the mass density, Young's modulus and shear modulus of the GPL spherical cap without interior cavities, respectively. In addition, e 0 and e * 0 (0 ≤ e 0 (e * 0 ) < 1) refer to the coefficients of porosity for the first two profiles, respectively; e m and e * m stand for the mass density coefficients for these two distributions, respectively; α and α are two parameters referring to a uniform porosity profile. For an increased size and density of the internal cavities, the porosity increases, with a subsequent reduction of the mechanical properties.
The relation between the elasticity modulus and density for an open-cell metal foam is assumed as [58,59] which is adopted to derive the relation between porosity and mass density coefficients for various porosity patterns; i.e., Here, we assume that the masses of spherical caps with various porosity patterns and GPL dispersions are identical. To compare the stiffness of different distributions, indeed, the analyses should be implemented for shells with equal masses. Hence, the values of e * 0 and α can be evaluated for a fixed value of e 0 [38,39], as According to Equation (6), the values of e * 0 and α can be estimated for a fixed value of e 0 , as shown in Table 1. It can be seen that e * 0 increases as the value of e 0 increases. When e 0 equals 0.6, e * 0 becomes equal to 0.9612, which is near to the upper bound. Hereafter, e 0 ∈ [0, 0.6] is used within the numerical investigation. According to the Halpin-Tsai micromechanics model [60], the elasticity modulus of nanocomposites without internal cavities is defined as where indices m and GPL stand for properties of the metallic matrix and graphene platelets, respectively; V GPL is the volumetric content of GPLs; and l GPL , w GPL and t GPL refer to the length, width and thickness of the nano-filler platelets, respectively. Based on the rule of mixtures, the mass density and Poisson's ratio of nanocomposite materials are defined as [61][62][63] whereas refers to the associated shear modulus. The volumetric content of GPLs, V GPL , is assumed to vary throughout the spherical cap's thickness, having five different dispersion patterns (see Figure 2): where S i1 , S i2 , S i3 , S i4 and S i5 denote the upper limits of V GPL ; and subscript i = 1, 2 or 3 refers to various porosity distributions within each pattern. Moreover, V T GPL stands for the total volumetric content of GPLs, which is defined in terms of the nanofiller weight fraction ∆ GPL in the following form: This is, in turn, used to derive S i1 , S i2 , S i3 , S i4 and S i5 as

Governing Equations of the Problem
The stress-strain relations are defined in matrix form as where the stress and strain field, together with the elasticity matrix D, read as follows: where υ * is the Poisson's ratio and E * denotes the Young's modulus that depends on the r coordinate. Based on the linear elasticity theory, the strain field in spherical coordinate is defined as Nanomaterials 2023, 13, 1205 and In addition, u, v and w define the kinematic components along r, φ and θ directions, respectively. According to the above-mentioned relations, the linear strain relation can be rewritten as ε L = LQ (25) where Q is the displacements vector and L is an operator matrix involving the partial derivatives of a function

Finite Element Modeling
A FEM-based approach is now adopted to solve the governing equations of the problem, where the spherical cap is divided into 8-node linear brick elements. For element (e), the 3D kinematic field is approximated as where Φ is the matrix of linear shape functions in spherical coordinates, whereas Λ (e) refers to the nodal displacement vector of the element, which is defined as The components of Φ are V being the volume of each element; i.e., 1 ξ 5 η 5 ζ 5 ξ 5 η 5 ξ 5 ζ 5 η 5 ζ 5 ξ 5 η 5 ζ 5 1 ξ 6 η 6 ζ 6 ξ 6 η 6 ξ 6 ζ 6 η 6 ζ 6 ξ 6 η 6 ζ 6 1 ξ 7 η 7 ζ 7 ξ 7 η 7 ξ 7 ζ 7 η 7 ζ 7 ξ 7 η 7 ζ 7 and It is also ξ = r cos θ sin φ, η = r sin θ sin φ, ζ = r cos φ (35) where r i , θ i and φ i are the nodal coordinates and A ij is obtained by elimination of the ith row and jth column from V. Substituting Equation (26) into Equation (23) gives the strain matrix of element (e) as where The FEM-based governing equations are determined from the principle of virtual work, where the potential energy U and virtual work of external loads δW are defined as A being the area under the external radial load σ rr = 1, which is subjected to the external surface of a spherical cap. In a pre-buckling state, the displacement field can be considered to be small, and the nonlinear terms of strain-displacement relations vanish. Therefore, one may write

D(BΛ (e) )dV = δΛ (e) T KΛ (e)
Therefore, based on the principle of virtual work, the static balance equation of the problem for each element in a pre-buckling state takes the following form Equation (43) in compact form can be written as where K (e) is the linear stiffness matrix and F (e) is the force matrix for each element, defined as By assembling each element matrix, the equilibrium equation of the spherical cap in the pre-buckling state is obtained as whose solution is determined in terms of the strain field in pre-buckling state for σ rr = 1. Afterward, the stress field due to these deformations is used in the geometric stiffness matrix, as detailed in the following. Finally, in order to determine the governing equations of an instability problem, the following equation can be used: Therefore, based on Equations (44) and (48), it is In line with Equation (22), the linear and nonlinear terms of the kinematic relations have been considered in the strain energy of the shell. In the pre-buckling state, the radial displacement components or large deformations can be assumed to be small, whereas only the linear strain terms appear. At the buckling state, instead, the nonlinear kinematic relations have to be considered. Therefore, the following relation can be obtained: Π Ext. can be rewritten as In the last relation, S 0 refers to the stresses obtained in pre-buckling state. By substituting Q (e) = ΦΛ (e) , we obtain where More in detail, it is Thus, according to Equation (50), we have δ 2 Π = δ (Λ (e)T )K (e) δ (Λ (e) ) + δ (Λ (e)T )( Equation (50) can be redefined in the following form: After the assembly of the element matrices, the following determinant should be assumed as null for the structure.
Note that K and K G refer to the linear stiffness matrix and geometric stiffness matrix, respectively, which are computed using the Gauss 8-point numerical integration rules.
Hereafter, we assume the following clamped boundary conditions:

Numerical Results and Discussion
In this section, we discuss the numerical results in terms of buckling loads of an FG-GPL, porous spherical cap with clamped boundary conditions, for various volume or weight fractions of GPL, and for different GPL distribution patterns, porosity distributions and coefficients, along with two polar angles of the FG-GPL, porous spherical shell.

Validation
In order to verify our results, we started the analysis with a comparative evaluation of the buckling predictions using the commercial Ansys Workbench code, for an isotropic homogeneous spherical cap. Hence, the following changes were considered in our study: e 0 = 0, γ GPL = 0. As far as the mechanical properties are concerned, we assumed E m = 130 GPa, ρ m = 8960 kg/m 3 , ν m = 0.34 for the copper material. As geometrical dimensions, we assumed: a = 0.225 m, b = 0.25 m, θ = 180 • , φ = 180 • , 90 • . In this way, the FG-GPL porous structure changes to an isotropic homogenous structure. The comparison between our results and predictions from Ansys Workbench is shown in Table 2, with an excellent agreement among them.

Parametric Analysis of the Buckling Load
We now study the effects of two polar angles, porosity coefficient, porosity distribution, GPL patterns and the weight fraction of GPL nanofillers (for the first time) on the buckling load of an FG-GPL, porous spherical cap with clamped boundaries. Hence, the following geometrical properties are considered: a = 0.225 m, b = 0.25 m, θ = 180 • , φ = 180 • , 90 • . The materials is characterized by the following mechanical properties: E m = 130 GPa, ρ m = 8960 kg/m 3 and ν m = 0.34 for the copper material [28]; and E GPL = 1.01 TPa, ρ GPL = 1062.5 kg/m 3 , ν GPL = 0.186, w GPL = 1.5 µm, l GPL = 2.5 µm and t GPL = 1.5 nm for GPLs. Table 3 indicates the influences of two types of polar angles and various GPL patterns on the buckling load of the FG-GPL-reinforced, porous spherical structure, under the assumptions PD3, e 0 = 0.2 and γ = 0.01 wt%. The extreme values of buckling load are related to GPLX and GPL-O distributions, respectively, which means that, when GPLs accumulate around the inner and outer surfaces of the shell, the stiffness reaches its highest value. Moreover, when GPLs are sparser around the outer and inner surfaces of the shell, the minimum buckling load is obtained. Note also that the critical buckling loads for a GPLA and GPL-V distributions are approximately the same. The results also show that the surface area of the shell increases by increasing the polar angle, and there are consecutive increases in the structural stiffness and buckling load. The influences of two polar angles and various porosity distributions are reported in Table 4 (GPLX, e 0 = 0.4, γ = 0.01 wt%), which shows that the maximum and minimum buckling loads belong to PD1 and PD2 distributions, respectively. PD1 provides higher structural stiffness, and PD2 gives the minimum stiffness of the spherical cap shell. The main difference between the extreme values of critical buckling load is approximately 90% for different porosity patterns. This means that the porosity distribution has a considerable effect on the buckling loads of FG-GPL, porous, spherical cap shells. The influences of two polar angles and weight fractions of nanofillers on the buckling loads of FG-GPL, porous spherical shell (PD1, e 0 = 0.5, GPLX) are given in Table 5. Note that, by increasing the weight fraction of GPLs, the critical buckling loads of shell increases significantly (approximately 100%), along with a small variation of the structural mass. This issue can be useful for aerospace structures where the high stiffness and low density are extremely important.  Table 6 shows the effect of the porosity coefficient on the critical buckling loads of FG-GPL, porous spherical shells (PD1, γ = 0.01 wt%, GPLX). When the porosity of the structure increases, the critical buckling load of FG-GPL porous spherical shells decreases, because of the decreased structural stiffness. The comparative evaluation of Tables 3-6 denotes that the influence of the porosity coefficient on the critical buckling load is lower than the GPL pattern and weight fraction of the nanofiller (its impact is approximately equal to 12.5%). On the other hand, the effects of the GPL pattern, porosity distribution and weight fraction of nanofillers on the critical buckling loads of FG-GPL, porous spherical shells are more pronounced than the porosity coefficient. The first six buckling mode shapes are shown in Figures 3 and 4 for the clamped spherical caps with ϑ = 180 • , φ = 180 • and 90 • , respectively. The figures clearly show that the first two buckling mode shapes and loads for each polar angle of spherical cap are the same. It is also observable that the number of buckling waves increases for higher buckling modes.

Concluding Remarks
The present work has studied the buckling responses of spherical caps made of FG porous materials reinforced by GPLs. Three different porosity distributions and five GPL patterns have been considered, along with the shell thickness. The equilibrium equations for the pre-buckling state have been determined according to the linear 3D elasticity theory and the principle of virtual work, whereas the buckling load associated with the problem has been computed according to the nonlinear Green strain field and generalized geometric stiffness concepts. We have studied the influence of the GPL pattern, weight fraction of nanofillers, porosity coefficient, porosity distribution and polar angles on the buckling loads of porous spherical cap made of FG-GPLs. Based on a large systematic investigation, the final remarks can be summarized as follows: (a) The maximum and minimum buckling loads seem to be reached for GPL-X and GPL-O distributions, respectively. (b) The maximum and minimum buckling loads belong to the PD1 and PD2 cases, respectively. (c) The difference between the maximum and minimum critical buckling loads for different porosity distributions is approximately equal to 90%, and the buckling loads of the selected structure increase considerably (approximately of 100%) with an increase in the weight fraction of GPLs. (d) The effect of the porosity coefficient on the critical buckling load for porous spherical cap shells made of FG-GPLs is lower than the weight fraction of the nanofillers, being approximately equal to 12.5%.
Such results could be useful for designing similar shell members with optimized mechanical properties and structural performances, as required by various engineering applications.