Large Deflection Geometrically Nonlinear Bending of Porous Nanocomposite Cylindrical Panels on Elastic Foundation

: Large deflection nonlinear bending of functionally graded (FG) porous cylindrical panels reinforced with graphene platelets (GPLs) on a Pasternak-type elastic foundation is examined by developing a reliable and effective 2D meshfree-based nonlinear numerical method. The large displacement field is express by the first-order shear deformation theory (FSDT) and the von K á rm á n nonlinearity, and approximated by 2D natural element method (NEM) in conjunction with the stabilized MITC3+ shell concept and the shell surface–rectangular grid geometry transformation. The nonlinear simultaneous equations are solved by a load incremental Newton–Raphson scheme. The developed nonlinear numerical method is justified from by comparing with the reference solutions, and the load–deflection and bending moment of FG-GPLRC porous cylindrical panels on elastic foundation are scrutinizingly examined. Four different symmetric GPL distribution patters (except for FG-Λ ) and three different symmetric porosity distributions are considered and their combined effects on the nonlinear bending behavior are investigated, as well as the effects of foundation stiffness and GPL amount. Also, the results are compared with those of FG CNT-reinforced porous cylindrical panels.


Introduction
Composite structures have been widely adopted in various science and engineering fields over several decades because they could maximize the desired performance and reduce the undesired function at the same time.A representative one would be a fiberinserted composite, which has been successfully used in the transportation industry to enhance structural strength while reducing weight.Furthermore, lately, these composite materials have greatly evolved according to the advances in materials science engineering, such as nanocomposite materials reinforced by graphene platelets or carbon nanotubes.Both GPLs and CNTs have attained a great spotlight as greatly promising nanofillers because of their excellent physical properties [1,2].These nanofillers are very strong, so the effective stiffness of the reinforced composite materials with these nanofillers can be surprisingly increased, even when a tiny number of nanofillers are added [3].However, owing to the cost problem, a limited number of nanofillers are not uniformly distributed in the matrix through the thickness with the specific functional gradient [4,5].Reinforced composites (RC) with a functionally graded distribution of nanofillers are called FG-GPLRCs or FG-CNTRCs [5][6][7].
When compared with GPLs, CNTs are in the form of cylinders, so the manufacturing cost is relatively higher, and their mechanical properties are anisotropic, which is dependent on the alignment [8][9][10][11].Furthermore, it has been reported that GPLs exhibit higher elastic modulus, tensile strength and fracture toughness than CNTs, besides the cost-effectiveness [12].For this reason, the mechanical behaviors of FG-GPLRCs have been Symmetry 2024, 16, 224 2 of 20 intensively investigated for developing advanced next-generation functional nanocomposites.Meanwhile, GPL-reinforced composites are characterized by a closed-cell porous form exhibiting low density, high porosity and large specific surface, which provides the possibility for combining the superior properties of GPLs and porous forms [13].Early studies concentrated on the linear responses of FG-GPLRC porous structures in order to investigate the basic responses like static bending, free vibration and buckling.For example, Li et al. [14] numerically investigated these three basic behaviors of FG-GPLRC porous plates by employing an isogeometric analysis.Yang et al. [15] examined the buckling and postbuckling behaviors of FG-GPLRC multilayer beams using the FSDT.Arefi et al. [16] solved the free vibration of FG-GPLRC nanoplates using a two-variable sinusoidal shear deformation theory (SSDT).Jafari and Kiani [17] analyzed the free vibration of FG-GPLRC plates using a quasi-3D shear and normal deformable model.Cho [18] numerically investigated the free vibration of FG-GPLRC porous cylindrical panels by 2D NEM.A more detailed extended survey of the early studies on FG-GPLRC structures may be referred to in review paper [19].
Regarding the large deflection bending of FG GPL-reinforced composites, Shen et al. [20] examined the temperature-dependent nonlinear bending of FG-GPLRC laminated plates on elastic foundations by applying a two-step perturbation scheme to a higher-order shear deformation theory (HSDT) with the von Kármán nonlinearity.Gholami and Ansari [21] analytically investigated a large deflection geometrically nonlinear bending of FG-GPLRC plates by applying the virtual work principle to the SSDT.Sahmani et al. [22] analytically solved the size-dependent nonlinear bending of FG-GPLRC porous micro/nano beams by the Galerkin method together with an improved perturbation scheme.Shen et al. [23] analytically examined the temperature-dependent nonlinear bending behavior of FG-GPLRC cylindrical panels on elastic foundations by using the HSDT.Wang.et al. [24] investigated the nonlinear bending behavior of FG-GPLRC microbeams in thermal environments by applying the Newton-Raphson method to Timoshenko beam theory with the modified couple stress theory.Liu et al. [25] presented an analytical method for solving the nonlinear static bending of FG-GPLRC porous arches by using the virtual work principle based on the Euler-Bernoulli theory.Tam et al. [26] numerically investigated the influence of an open crack on the nonlinear bending of multilayer FG-GPLRC beams elastically restrained at both ends.Wang et al. [27] analytically investigated the nonlinear static response of FG polymer-based circular microarches reinforced by graphenes by developing a modified couple stress-based model based on the exponential FSDT.Anirudh et al. [28] numerically examined the nonlinear bending of porous FG-GRC curved beams using a three-noded curved beam finite element based on a FSDT.Nguyen et al. [29] proposed a numerical model for investigating size-dependent geometrically nonlinear static and dynamic responses of FG microplates reinforced by graphene nanofillers by using four-variable refined plate theory within the framework of isogeometric analysis.Songsuwan and Wattanasakulpong [30] analytically examined the linear and nonlinear bending of FG-GPLRC beams by applying the Gram-Schmidy-Ritz method to the Reddy's HSDT [31] in conjunction of a von Kármán nonlinearity.
The literature survey on the nonlinear deflection FG-GPLRC structures reveals that the studies were mostly focused on beams and plates using analytical methods.In particular, the numerical studies using the meshfree method were rarely reported, and the effects of associated parameters on the nonlinear bending behaviors, for example, the combined effects of GPL and porosity distributions by including the elastic foundation, were not sufficiently investigated.In this context, this study aims at the scrutinizing investigation of the nonlinear bending responses of FG-GPLRC porous cylindrical panels resting on an elastic foundation by developing a reliable and effective meshfree-based 2D nonlinear numerical method.The large deflection of cylindrical panel is expressed by the FSDT together with the von Kármán geometry nonlinearity, but the approximated stresses are enhanced by applying the strain recovery [32] and the Reddy's TSDT.The equivalent elastic properties of FG-GPLRC cylindrical panels are determined by applying the Halpin-Tsai homogenization model, and those are modified by introducing three different porosity distributions.
The nonlinear numerical method is developed in the framework of 2D NEM by introducing a geometry transformation from the shell mid-surface to a 2D planar NEM grid to relax the complex mathematical derivation on the curved shell surface.The 2D NEM is a last introduced meshfree method for the effective analysis of 3D elastic structures and shows a relatively higher accuracy even at coarse grids [18].The concept of the MITC3+ shell element is employed to suppress the troublesome shear membrane locking [33,34], and is further stabilized by introducing two stabilization parameters associated with the porosity and nonlinear iteration number.The constructed nonlinear matrix equations are solved by introducing a load incremental Newton-Raphson iteration scheme, by including the stiffness matrix increments.The benchmark experiment is performed to justify the introduced nonlinear numerical method.And the parametric experiments are performed to profoundly investigate the large deflection bending behaviors of FG-GPLRC porous cylindrical panels on elastic foundation with respect to the associated parameters, particularly to the symmetric GPL and porosity distributions (except for FG-Λ).The comparison with the FG-CNTRC cylindrical panels for the same volume fractions of nanofillers is also presented.

FG-GPLRC Porous Cylindrical Panel Resting on Elastic Foundation
A cylindrical panel reinforced with graphene platelets (GPLs) is shown in Figure 1a, where a coordinate system (x, y, z) is introduced on the panel mid-surface ϖ with the relation of x = Rθ.The geometry of the cylindrical panel is governed by radius R, subtended angle θ 0 , length L, and uniform thickness h.The panel is rested on elastic foundation which is assumed to be firmly attached to the panel bottom without causing debonding, and its pressure-deflection relationship p − w is assumed to be where K w is the Winkler stiffness, K s is the shearing layer stiffness and ∇ 2 ϖ is the Laplacian operator defined on the panel mid-surface ϖ, respectively.Graphene platelets are distributed through the thickness in a specific functionally graded pattern.Figure 1b shows four primitive GPL distribution patterns taken for this study, where GPLs are uniformly distributed in FG-U while those are biased to the mid-surface in FG-O, to the upper surface in FG-X, and to the lower surface in FG-Λ, respectively.These four patterns sufficiently represent all the possible GPL distribution patterns, and three patterns, FG-U, FG-O and FG-X, are symmetric with respect to the panel mid-surface.This study assumes that GPLs are uniformly distributed within the matrix and behave as an equivalent rectangular solid fiber with the width GPL w , length GPL l , and thickness GPL t , and the graphene-reinforced composites are assumed to be isotropic.Letting V m (z) and V GPL (z) be the matrix and GPL volume fractions of GPLs in function of thickness co-ordinate z, both must obey the physical relation given by Here, the GPL volume fraction V GPL (z) has one of following thickness functions given by depending on the above GPL distribution pattern, in which the total GPL volume fraction V * GPL is determined by in terms of the GPL mass fraction g GPL and the matrix and GPL densities ρ m and ρ GPL .This study assumes that GPLs are uniformly distributed within the matrix and behave as an equivalent rectangular solid fiber with the width w GPL , length l GPL , and thickness t GPL , and the graphene-reinforced composites are assumed to be isotropic.And the effective elastic modulus E e f f of GPLRC is calculated using the Halphin-Tsai approach [35], which gives with in which E m and E GPL are the elastic moduli of matrix and GPLs, and the geometric parameters ξ L and ξ T are determined by Meanwhile, the effective Poisson's ratio ν e f f and density ρ e f f of GPLRC are calculated as ) using the linear rule of mixture.
Referring to Figure 2, this study considers three symmetric porosity distributions: center-biased (PD1), outer-biased (PD2) and uniform (PD3), which are defined by PD3 : ψ(z) = e 0 (12) with the porosity parameter e 0 (0 ≤ e 0 ≤ 1).Both the magnitude and distribution of porosity affect the magnitudes of the effective elastic modulus E e f f and shear modu-lus G e f f = E e f f /2 1 + ν e f f of GPLRC structures.Letting ℘(z) be the effective elastic properties of porous GPLRC structures, they are determined by should be modified using the relationship of ( ) ( ) between the elastic modulus and the density [36].Letting e be the porosity parameter for the porous CNTRC structures, it is determined by   But, for the effective density ρ e f f , the porosity parameter e 0 in Equations ( 10)-( 12) should be modified using the relationship of between the elastic modulus and the density [36].Letting e be the porosity parameter for the porous CNTRC structures, it is determined by for PD1, for example.The displacement field u = u x , u y , u z T in the first-order shear deformation shell theory is expressed as with the displacement component vector d = u 0 , v 0 , w 0 , ϑ x , ϑ y T at the panel mid-surface.
Here, (u 0 , v 0 , w 0 ) are translations while ϑ x , ϑ y are rotations of the mid-surface.The large deflection bending of the GPLRC cylindrical panel is modeled by a von Kármán-type geometric nonlinear model.Then, the nonlinear strain fields are expressed as with r = R + z ≈ R, where the (3 × 5) and (2 × 5) partial derivative matrices A L , A NL and A s are defined by with A x = ∂/∂x and A y = ∂/∂y.Here, w * ,x = ∂w * /∂x and w * ,y = ∂w * /∂y denote the derivatives of panel deflection which are assumed to be determined previously, as will be given later.Then, the constitutive relations between stresses and strains are expressed as

Natural Element Approximation of Large Deflection Bending
In order to compute the geometrically nonlinear deformation of an FG-GPLRC porous cylindrical panel using 2D NEM, the panel mid-surface ϖ is divided into a number of nodes and Delaunay triangles, as shown in Figure 3.And the displacement field u(x, y, z) is approximated as using Laplace interpolation (L/I) functions ψ I (x, y) [37,38] and the nodal displacement  A and 2 A : Then, the natural element approximation of the in-plane strains ε in Equation (11) and the transverse shear (T/S) strains γ in Equation ( 18) arrives at Then, L/I functions ψ I (x, s) are transformed to φ I (ζ 1 , ζ 2 ) on the rectangular NEM grid, and the relations of x = R • ζ 1 and y = ζ 2 introduce the inverse Jacobi matrix J −1 , defined by Furthermore, the partial derivatives A x and A y in Equations ( 13) and ( 14) on the panel mid-surface are switched to On the 2D planar NEM grid, according to the chain rule.
Introducing Equation (27) into Equations ( 19)-( 21) leads to ÂL , ÂNL and Âs in which A x and A y are replaced with A 1 and A 2 : Then, the natural element approximation of the in-plane strains ε in Equation ( 11) and the transverse shear (T/S) strains γ in Equation ( 18) arrives at However, the standard NE approximation (30) of the T/S strain γ h using C 0 -L/I functions φ I may frequently suffer from shear-membrane locking [22,23].One way to overcome this problem is to indirectly interpolate the T/S strains by employing the notion of the MITC3+ shell element [39], as described in the previous work [40].According to this concept, the NE approximation u h is re-interpolated using three-node bilinear triangular basis functions [41], from which the element-wise locking-free T/S strains γ h e are evaluated at six tying points [39] within each triangular element.And the resulting re-interpolated T/S strains γ h e are expressed as with the (2 × 15) element-wise matrices Be (ξ, η, z, R) and the (15 × 1) element-wise nodal vectors b e = {b e 1 , b e 2 , b e 3 }.Meanwhile, the Galerkin variational form for the geometrically nonlinear bending analysis of FG-GPLRC porous cylindrical panels on resting on elastic foundation is derived from the energy principle [42], given by Introducing Equations ( 29) and (31) into Equation ( 29), together with the stress-strain relations (22) and (23), one can obtain the nonlinear simultaneous equations given by to solve the nonlinear nodal displacement d I N

I=1
. The linear and nonlinear stiffness matrices K L and K NL and the load vector F in Equation ( 33) are defined in Appendix A. Here, K e L,s is the linear stiffness matrix derived using the notion of MITC3+ shell element [39].
To solve the nonlinear matrix Equation ( 33), the incremental loading scheme and the Newton-Raphson method are combined.In the uniform incremental loading, the load vector F at load step k is calculated by with F 0 being a unit-step load vector.By simplifying Then, the linearization of Taylor series expansion of ℵ(u) becomes Here, the gradient ∂K NL /∂u of the nonlinear stiffness matrix K NL is obtained as follows: By enforcing that ℵ(u + ∆u) in Equation ( 35) be zero, it is not hard to obtain the following iterative matrix equations: for step n + 1.The Newton-Raphson iteration (37) terminates when the following stop criterion is satisfied: which is defined in the Euclidean norm.And the solution u n+1 at the iteration stage n is updated as u n+1 = u n + ∆u n+1 .

Benchmark Test
The developed nonlinear numerical method is verified with the simply supported isotropic and FGM cylindrical panels under the uniform distributed load q 0 .Referring to Figure 3, the simply-supported condition (SSSS) in FSDT is enforced as w 0 = v 0 = ϑ y = 0 for sides 2  ⃝ and 4 ⃝ at θ = 0 and θ 0 and w 0 = u 0 = ϑ x = 0 for sides 1 ⃝ and 3 ⃝ at x = 0 and L, respectively.The panel mid-surface is uniformly divided into an 11 × 11 NEM grid, and the stiffness matrices and load vector in Equations (A5)-(A9) are numerically integrated using seven Gauss points.The load q 0 is uniformly divided into 400 increments and the convergence tolerance ε T in Equation ( 38) is set by 1.0 × 10 −2 , unless otherwise stated.
The first example is an isotropic cylindrical panel with the elastic modulus E m of 70.0 GPa and the Poisson's ratio ν of 0.3.The geometry dimensions are R = 1.0 m, b = L = 0.2 m and h = 0.01 m, and the porosity and the elastic foundation are not considered.The central deflection w c and the uniform load q 0 are calibrated as w c = w c /h and q 0 = q 0 b 4 / E m h 4 , respectively.The convergence test was performed to the density of NEM grid, as given in Table 1.The uniform 15 × 15 NEM grid shows a relative difference w re f c /h less than 5.0%, so this grid density is chosen for all of the numerical experiments in this paper.The load-deflection plots are comparatively represented in Figure 4a with Shen et al. [23], Maleki and Tahani [43] and Thang et al. [44].An excellent agreement between the present method and three reference solutions is shown up to the nondimensional load q 0 = 45.A second example is a ZrO 2 /Al FGM cylindrical panel, and the zirconia volume fraction f c (z) is governed by (0.5 − z/h) N .The geometry dimensions and the properties of Al are same with the previous isotropic cylindrical panel, but the material of ZrO 2 is E c = 151.0GPa and ν = 0.3.The effective elastic properties of ZrO 2 /Al FGM cylindrical panel are estimated using the Mori-Tanaka approach given in Appendix B. The load-deflection curves for two different ceramic power-law indices N = 0.2 and 2.0 are comparatively represented in Figure 4b with Shen and Wang [45] and Zhao and Liew [46].One can clearly see that the present method fairly agrees with two reference solutions up to q 0 = 70 for N = 2.0.70 for N = 0.2. .A good agreement with Shen and Wang [45] is seen for the load-central deflection plots, but the present method produces the load-bending moment curves lower than those of Shen and Wang [45].[45] is seen for the load-central deflection plots, but the present method produces the load-bending moment curves lower than those of Shen and Wang [45].
to the Mori-Tanaka approach.The bending moment is defined by .A good agreement with Shen and Wang [45] is seen for the load-central deflection plots, but the present method produces the load-bending moment curves lower than those of Shen and Wang [45].To obtain better numerical results, the iteration compensation parameter β(n) = √ 1 + 4n/400 was introduced to the shear elastic modulus G e f f in Equation ( 23) (see also Equation (A10) in Appendix A) and the strains and stresses are enhanced by applying the strain recovery scheme [32] and Reddy and Liu's HSDT, given in Appendix A. The introduction of β(n) was motivated that the variation in stiffness matrix K NL in Equation (A8) along the nonlinear iteration requires the adaption of K e L,s in Equation (A6), which was derived according to the MITC3+ shell element concept.The improvement results are given in Figure 6a, where the combined use of β(n) and the higher-order shell theory provides the load-bending moment curve, which is closer to the reference solution.Furthermore, from Figure 6b, it is seen that the introduction of β(n) slightly reduces the magnitude of nondimensional central deflection such that the curve moves towards the reference solution.So, all of the load-deflection and load-bending moment plots in this paper are obtained by applying the above-mentioned combined use.

Nonlinear Behavior
Next, the large deflection bending of FG-GPLRC porous cylindrical panels on elastic foundation is scrutinizingly examined by changing the associated parameters.The simply supported cylindrical panels subject to uniform pressure 0 q is taken, and its geometric dimensions are , and the effects of temperature and humidity are not considered.According to Yasmin and Daniel [47] and Rafiee et al. [12],

Nonlinear Behavior
Next, the large deflection bending of FG-GPLRC porous cylindrical panels on elastic foundation is scrutinizingly examined by changing the associated parameters.The simply supported cylindrical panels subject to uniform pressure q 0 is taken, and its geometric dimensions are R = 10.0 m, b = L = 0.2 m and h = 0.01 m.The matrix is epoxy with the elastic properties of E m = 3.0 GPa and v m = 0.34, and the effects of temperature and humidity are not considered.According to Yasmin and Daniel [47] and Rafiee et al. [12], the geometry dimensions of GPLs are l GPL = 2.5 µm, t GPL = 1.5 nm and w GPL = 1.5 µm and the elastic properties are E GPL = 1.01 TPa and v GPL = 0.186, respectively.
Figure 7a,b represent the load-central deflection and bending moment plots of a nonporous FG-U GPLRC cylindrical panel with the GPL mass fraction g GPL of 0.4% for three different values of elastic foundation.Both the magnitude and nonlinearity of the central deflection and bending moment uniformly decrease proportionally to the foundation stiffness.The reason is because the deflection is suppressed when the panel is rested on the elastic foundation and the suppression intensity becomes larger as the foundation stiffness increases.Figure 8a,b comparatively represent the load-central deflection and load-bending moment plots for four different GPL distribution patterns when the foundation stiffness is K w , K s = 10 2 , 0 .In the load-deflection curves, FG-O produces the highest level while FG-X exhibits the lowest one, because GPLs are concentrated at the mid-surface in FG-O, while they are at the top and bottom surfaces in FG-X.The panel stiffness becomes larger as GPLs are biased far from the mid-surface.Meanwhile, the order in the magnitude of bending moment is almost reversed, because the bending moment magnitude is proportional to the panel stiffness, contrary to the central deflection.
The influence of GPL total mass on FG-GPLRC cylindrical panel on elastic foundation was analyzed and the plots are comparatively represented in Figure 9a,b.The central deflection and its nonlinearity decrease in proportion to g GPL , because the panel stiffness becomes larger as the GPL amount increases.On the hand, the bending moment uniformly increases with the increase in the GPL amount, since the elastic modulus increase of the panel is superior to the central deflection decrease according to the increase in GPL amount.
is ( ) ( ) In the load-deflection curves, FG-O produces the highest level while FG-X exhibits the lowest one, because GPLs are concentrated at the mid-surface in FG-O, while they are at the top and bottom surfaces in FG-X.The panel stiffness becomes larger as GPLs are biased far from the mid-surface.Meanwhile, the order in the magnitude of bending moment is almost reversed, because the bending moment magnitude is proportional to the panel stiffness, contrary to the central deflection.The influence of GPL total mass on FG-GPLRC cylindrical panel on elastic foundation was analyzed and the plots are comparatively represented in Figure 9a,b.The central deflection and its nonlinearity decrease in proportion to GPL g , because the panel stiffness becomes larger as the GPL amount increases.On the hand, the bending moment uniformly increases with the increase in the GPL amount, since the elastic modulus increase of the panel is superior to the central deflection decrease according to the increase in GPL amount.Next, the nonlinear bending behaviors of FG-GPLRC cylindrical panels resting on elastic foundation are analyzed.Figure 10a represents the effect of porosity distribution on the load-deflection curve of FG-Λ GPLRC cylindrical panel.Referring to the previous work [18], the porosity stabilization parameter α in Equation (A9) is used as times (by vertical integration of Equations ( 10)-( 12)).On the other hand, this trend is completely reversed for the load-bending moment curves, as shown in Figure 10b.The reason is the same with the previous Figure 9b for the influence of the GPL total mass.Next, the nonlinear bending behaviors of FG-GPLRC cylindrical panels resting on elastic foundation are analyzed.Figure 10a represents the effect of porosity distribution on the load-deflection curve of FG-Λ GPLRC cylindrical panel.Referring to the previous work [18], the porosity stabilization parameter α in Equation (A9) is used as α = 1 for PD1, for PD3, respectively.The occurrence of porosity increases the central deflection, and porosity distribution 3 (i.e., uniform distribution) produces larger central deflection than the center-(PD1) and outer-biased (PD2) porosity distributions.This is because the porosity decreases the panel stiffness and the total amounts of pores in PD1 and PD2 are less than that of PD3 by 2/π times (by vertical integration of Equations ( 10)-( 12)).On the other hand, this trend is completely reversed for the load-bending moment curves, as shown in Figure 10b.The reason is the same with the previous Figure 9b for the influence of the GPL total mass.distribution) produces larger central deflection than the center-(PD1) and outer-biased (PD2) porosity distributions.This is because the porosity decreases the panel stiffness and the total amounts of pores in PD1 and PD2 are less than that of PD3 by π / 2 times (by vertical integration of Equations ( 10)-( 12)).On the other hand, this trend is completely reversed for the load-bending moment curves, as shown in Figure 10b.The reason is the same with the previous Figure 9b for the influence of the GPL total mass.The influence of porosity parameter 0 e (i.e., the porosity intensity) is also analyzed, and the results are given in Figure 11a,b for four different porosity parameters.As is expected, the central deflection uniformly increases proportionally to the value of 0 e , while The influence of porosity parameter e 0 (i.e., the porosity intensity) is also analyzed, and the results are given in Figure 11a,b for four different porosity parameters.As is expected, the central deflection uniformly increases proportionally to the value of e 0 , while the bending moment monotonically decreases proportional to e 0 .It is, of course, the panel stiffness is reversely proportion to the porosity intensity.The axial stress σ xx is extracted from the center of the panel at the final load step q 0 = 800 and calibrated as σ xx = σ xx h 2 / q 0 b 2 .Figure 12a shows the thicknesswise axial stress distributions for four different GPL distributions when the PD3 with e 0 = 0.6 is used.The linear variation is seen for FG-U because GPLs are uniformly distributed, while FG-O and FG-X show almost antisymmetric stress distributions to each other because GPLs are center-and outer-biased in FG-O and FG-X.Meanwhile, FG-Λ leads to a quadratic variation with the peak at the bottom because GPLs in FG-Λ have a linear variation through the thickness with the peak at the bottom.The combination of a linear variation in axial strain and a linear variation in elastic modulus leads to a quadratic variation in σ xx .Figure 12b shows the influence of porosity distribution on the axial stress distribution, where PD1 and PD2 show an almost antisymmetric distribution because pores are center-and outer-biased in PD1 and PD2, respectively, differing from the uniform distribution in PD3.The combination of a linear variation in axial strain and center-and outer-biased porosity distributions leads to outerand center-biased thicknesswise stress distributions, respectively.
Next, the nonlinear bending of FG-GPLRC porous cylindrical panels on elastic foundation is analyzed when GPLs are replaced with CNTs by keeping the matrix and the simulation parameters unchanged.The CNTs are (10,10) single-walled [48] and their orthotropic properties are given in Table 2. Using the modified rule of mixtures, the effective elastic properties FG-CNTRC cylindrical panels are estimated as (1, 2, 3 = x, y, z) with the CNT volume fraction V CNT (= 1 − V m ) and the CNT efficiency parameters [49]  for reflecting the macro-micro interaction effect between the matrix and CNTs.13a,b compare the load-central deflection and load-bending moment plots between GPL and CNT for the same volume fractions 0.12 and 0.28, respectively.The axis of CNTs is aligned in the x-direction.The foundation stiffness is set by 10 2 , 0 and the GPL and porosity distributions are FG-O and PD1 with e 0 of 0.3.The central deflection and bending moment of GPLs are seen to be smaller than CNTs by at least 3.5 and 2.1 times, respectively.Thus, it is found that GPLs are more effective than CNTs in aspects of both the deflection and bending moment.Figure 14a compares the thicknesswise axial stress distributions between GPL and CNT at the half-iteration stage of 400 0 = q (i.e., iteration number n is 200) where the sine-like stress distributions are produced by the combination of the center-biased GPL and porosity distributions in FG-O and PD1.It is clearly seen that the stress level of GPL is significantly lower than that of GPL because the strength per unit volume of GPL is higher than that of CNT.Meanwhile, Figure 14b compares the iteration histories of axial stress distributions between GPL and CNT.It is clearly seen that the axial stress of CNT is increased larger than GPL along the nonlinear iteration.Figure 15a  Figure 14a compares the thicknesswise axial stress distributions between GPL and CNT at the half-iteration stage of q 0 = 400 (i.e., iteration number n is 200) where the sine-like stress distributions are produced by the combination of the center-biased GPL and porosity distributions in FG-O and PD1.It is clearly seen that the stress level of GPL is significantly lower than that of GPL because the strength per unit volume of GPL is higher than that of CNT.Meanwhile, Figure 14b compares the iteration histories of axial stress distributions between GPL and CNT.It is clearly seen that the axial stress of CNT is increased larger than GPL along the nonlinear iteration.Figure 15a,b show the dependence of the load-bending moment plot and the axial stress distribution on the CNT alignment, respectively.The bending moment M x and the axial stress σ xx become significantly smaller when the CNT axis is changed from 0 0 to 90 0 (i.e., from the circumferential-direction alignment to the shell-axis alignment).However, it was seen that the variation in central deflection to the change in CNT axis alignment is not remarkable, even though not represented.Meanwhile, in the shell-axis alignment, the bending moment M y and the axial stress σ yy show almost the same levels with M x and σ xx in the circumferential-direction alignment.

Conclusions
The geometrically nonlinear bending of FG-GPLRC porous cylindrical panels resting on an elastic foundation was scrutinizingly investigated by developing a reliable and effective nonlinear numerical method.The numerical method was based on 2D NEM and a geometry transformation between the panel mid-surface, and a 2D regular planar NEM grid and a stabilized locking-free scheme were introduced.The nonlinear displacement field was expressed by a FSDT with the von Kármán nonlinearity.The nonlinear numerical results were enhanced by combining the stabilization parameter and the HSDT.The proposed numerical method was justified from the benchmark experiments.The nonlinear bending behaviors of FG-GPLRC cylindrical panels were parametrically investigated in terms of the load-deflection and bending moment plots and the thicknesswise axial stress distributions.The numerical results provide the following major findings:

•
The reliability of the developed nonlinear numerical method was justified through isotropic and two metal-ceramic FGM cylindrical panels such that both the load-

Conclusions
The geometrically nonlinear bending of FG-GPLRC porous cylindrical panels resting on an elastic foundation was scrutinizingly investigated by developing a reliable and effective nonlinear numerical method.The numerical method was based on 2D NEM and a geometry transformation between the panel mid-surface, and a 2D regular planar NEM grid and a stabilized locking-free scheme were introduced.The nonlinear displacement field was expressed by a FSDT with the von Kármán nonlinearity.The nonlinear numerical results were enhanced by combining the stabilization parameter and the HSDT.The proposed numerical method was justified from the benchmark experiments.The nonlinear bending behaviors of FG-GPLRC cylindrical panels were parametrically investigated in terms of the load-deflection and bending moment plots and the thicknesswise axial stress distributions.The numerical results provide the following major findings: • The reliability of the developed nonlinear numerical method was justified through isotropic and two metal-ceramic FGM cylindrical panels such that both the loaddeflection and bending moment plots agree well with the representative reference solutions, even at coarse grids.• The introduction of β(n) and the higher-order shell theory enhance the nonlinear bending results such that both the load-central deflection and bending moment plots become closer to the reference solutions.• The large deflection of FG-GPLRC cylindrical panels on elastic foundation uniformly decreases proportionally to both the foundation stiffness and the GPL amount.Meanwhile, the nonlinear bending moment does also decrease with the former, but increases with the latter.

•
In the load-deflection curve, FG-O and PD3 show the highest level, while FG-X and PD1 display the lowest one.But this relative order is reversed in the load-bending moment plot.

•
Regarding the porosity, the large deflection uniformly increases but the resulting bending moment monotonically decreases as the porosity parameter increases.And the thicknesswise axial stress distribution is diversely affected by the combination of GPL and porosity distributions.• The nonlinear central deflection, axial stress and bending moment of FG-GPLRC porous cylindrical panels on elastic foundation remarkably increase when GPLs are replaced with CNTs, for the same volume fraction.
The current study was performed with the temperature-independent material properties and the spatially varying porosity distributions, but the temperature dependence and the porosity variation during operation would be worthwhile, and which represents a subject that deserves future work.
ϑ y T I .Here, the subscripts I indicate the I-th node within the NEM grid ℑ C constructed with N nodes and M Delaunay triangles.Symmetry 2024, 16, x FOR PEER REVIEW 8 of 22

Figure 3 .
Figure 3. L/I functions φ I (ζ 1 , ζ 2 ) on 2D planar NEM grid and their mapping to ψ I (x, y) on the panel mid-surface.The derivation of the L/I function and its calculation on the curved surface are difficult and painstaking.This problem could be effectively relaxed by introducing a geometry transformation T C from the physical NEM grid ℑ C = [0, Rθ 0 ] × [0, ↕] on the panel midsurface to a computational NEM grid ℑ R = [0, θ 0 ] × [0, ↕] on the rectangular plane the local co-ordinates ζ 1 and ζ 2 such that

Figure 4 .
Figure 4. Comparison of the load-deflection plots of cylindrical panels without elastic foundation (Ref_1: [23], Ref_2: [43], Ref_3: [44], Ref_4: [45], Ref_5: [46]): (a) isotropic; (b) zirconia/Al FGM.Third example is a Si 3 N 4 /SUS304 FGM cylindrical panel on elastic foundation with the ceramic power-law index of N = 5.0.The load-deflection and bending moment plots in Figure 5 are computed for three different foundation stiffnesses: K w , K s = (0, 0), 10 2 , 0 and 10 2 , 10 .The geometry dimensions are R = 0.2 m, b = L = 0.04 m and h = 0.002 m, and the porosity is not considered.The properties of base materials are as follows: E c = 322.27GPa and ν c = 0.24 for Si 3 N 4 and E m = 207.79GPa and ν m = 0.318 for SUS304, and the effective elastic properties are also determined according to the Mori-Tanaka approach.The bending moment is defined by M x = h/2 −h/2 zσ xx dz and extracted at the center of the panel, while two foundation stiffnesses are calibrated as K w = K w b 4 /D m and K s = K s b 2 /D m with D m = E m h 3 .A good agreement with Shen and Wang[45] is seen for the load-central deflection plots, but the present method produces the load-bending moment curves lower than those of Shen and Wang[45].

2 σ
and extracted at the center of the panel, while two foundation stiffnesses are calibrated as

Figure 7 .
Figure 7. Variation of the nonlinear bending of FG-U cylindrical panel on an elastic foundation: (a) load-deflection; (b) load-bending moment.

Figure.Figure 8 .
Figure 8a,b comparatively represent the load-central deflection and load-bending moment plots for four different GPL distribution patterns when the foundation stiffness is (

Figure 7 .
Figure 7. Variation of the nonlinear bending of FG-U cylindrical panel on an elastic foundation: (a) load-deflection; (b) load-bending moment.

Figure 8 .
Figure 8. Influence of the GPL distribution pattern on FG-CNTRC cylindrical panels on an elastic foundation: (a) load-deflection, (b) load-bending moment.

Figure 8 .Figure 9 .
Figure 8. Influence of the GPL distribution pattern on FG-CNTRC cylindrical panels on an elastic foundation: (a) load-deflection, (b) load-bending moment.Symmetry 2024, 16, x FOR PEER REVIEW 14 of 22 respectively.The occurrence of porosity increases the central deflection, and porosity distribution 3 (i.e., uniform distribution) produces larger central deflection than the center-(PD1) and outer-biased (PD2) porosity distributions.This is because the porosity decreases the panel stiffness and the total amounts of pores in PD1 and PD2 are less than that of PD3 by π / 2

Figure 9 .
Figure 9. Influence of the GPL total mass on FG-X cylindrical panels on an elastic foundation: (a) load-deflection, (b) load-bending moment.

Figure 10 .
Figure 10.Influence of the porosity distribution on FG-Λ porous cylindrical panels on an elastic foundation: (a) load-deflection, (b) load-bending moment.

Figure 10 .
Figure 10.Influence of the porosity distribution on FG-Λ porous cylindrical panels on an elastic foundation: (a) load-deflection, (b) load-bending moment.
,b show the dependence of the load-bending moment plot and the axial stress distribution on the CNT alignment, respectively.The bending moment x M and the axial stress xx σ become significantly smaller when the CNT axis is changed from 0 0 to 0 90(i.e., from the circumferential-direction alignment to the shell-axis alignment).However, it was seen that the variation in central deflection to the change in CNT axis alignment is not remarkable, even though not represented.Meanwhile, in the shell-axis alignment, the bending moment y M and the axial stress yy σ show almost the same levels with x M and xx σ in the circumferential-direction alignment.

Figure 14 .
Figure 14.Thicknesswise distributions of σ xx of FG-O porous cylindrical panel on elastic foundation: (a) comparison between GPLRC and CNTRC; (b) variation to the nonlinear iteration number.

Figure 15 .
Figure 15.Influence of CNT alignment on the nonlinear bending of FG-O porous cylindrical panel on elastic foundation: (a) load-bending moment; (b) nondimensional normal stresses.

Table 1 .
Convergence of the calibrated central deflection of isotropic cylindrical panel without elastic foundation at the final step q 0 = 45 (R = 1.0 m, b = L = 0.2 m, h = 0.01 m, SSSS).