Porosity-dependent stability analysis of bio-inspired cellular nanocomposite shells

This paper presents an innovative numerical model to investigate buckling behaviour of bio-inspired continuously graded porous (CGP) nanocomposite cylindrical shells. It is postulated that the shell subjected to combined lateral pressure and axial compressive load is constructed from metal foams with closed-cell structures that possess graded internal pores, which exhibit three types of continuously graded porosity profiles based on a power-law distribution. A scaling relation for the effective Young’s modulus of the cellular structure determined by a variational finite element method (FEM) is used. The effective constitutive law of an elastic isotropic metal matrix containing distributed elastic carbon-nanotubes (CNTs) is estimated in consideration of the impact of CNTs agglomeration using a continuum model based on the Eshelby–Mori–Tanaka (EMT) approach. In contrast to conventional approaches, the study employs Euler–Bernoulli beams to model stiffeners within the CGP shells. This choice allows for a more realistic representation of stiffener effects, as opposed to the prevalent approach of uniform smearing across the shell’s surface. The equilibrium equations of the CGP shell, based on the Reddy higher-order shell theory (RHST), is obtained through the application of the Euler equation. Subsequently, the equations for stability are obtained through the utilization of the variational method. This study emphasizes the effects of geometrical parameters, porosity variability, and distribution of CNTs on the buckling performance of the CGP shells. The intricate interplay between CNTs and porosity distributions critically influences the stability behaviour of CGP shells. CNTs arrangement remarkably impacts buckling behaviour at higher length-to-mean radius ratios, while symmetric porosity near the mid-surface significantly enhances stiffness. These findings provide valuable insights for designing closed-cell cellular stiffened shells with optimal porosity to enhance stability.


Introduction
Recently, concerns over the prospect of climate change and limited resources of materials and energy have stimulated a pressing demand for the development of resilient and efficient lightweight structures with outstanding multi-functional properties.Moreover, recent advances in high-tech industries require innovative engineering designs to address the tunable multi-functional mechanical properties, including high strength, stiffness, toughness, and energy absorption.Porous cellular materials are one of the emerging lightweight materials that can simultaneously fulfil multiple requirements from structural stiffness to thermal insulation.For instance, a captivating potential of sandwich porous metal foams shines in their role in creating the front structure of high-speed trains, shown in Fig. 1.This application merges lightweight design with structural integrity for pioneering the future of swift and efficient rail travel.Metal foams are porous materials that including the implementation of gradual transitions in either chemical compositions or microstructural characteristics [11].This phenomenon is exemplified in the squid beaks [12], which exhibit a gradual shift from stiff to compliant materials, as well as in the variation of porosity observed between trabecular and cortical bones [13].This has led to the emergence of the novel continuously graded porous (CGP) structures comprised of metal foams.Various theoretical models have been introduced to predict the vibrational and buckling behaviour of CGP beams and plates [14][15][16][17][18][19][20].However, despite the significant importance of cylindrical shells in various engineering disciplines such as mechanical, civil, ocean, and aeronautical engineering, there has been a scarcity of research dedicated to studying the stability of CGP shells.
In recent years, there has been a growing interest in investigating the stability and buckling behaviour of bio-inspired FG shells and plates, driven by their potential applications in diverse engineering fields [23][24][25][26][27][28].Duc et al. [29] conducted an investigation into the nonlinear buckling and post-buckling analyses of stiffened truncated conical sandwich shells.These shells were composed of functionally graded (FG) face sheets and a CGP core, and rested on the Winkler-Pasternak elastic foundation.The stability equations for the shell were derived using the first-order shear deformation theory (FSDT), incorporating a von Karman-Donnell type of kinematic non-linearity.Zhou et al. [30] conducted a nonlinear buckling analysis on cylindrical shells composed of FG porous graphene platelet-reinforced composites, considering pre-buckling effects and in-plane constraints.Their findings indicate that the optimal material distribution for this type of nanocomposite is a symmetrical distribution of both graphene platelets and porosity throughout the thickness of the shell.Nam et al. [31] conducted a study on the nonlinear buckling and post-buckling of cylindrical shells reinforced by orthogonal stiffeners resting on Pasternak elastic foundations, utilizing classical shell theory as the basis for their investigation.An investigation into the dynamic buckling analysis of a truncated conical shell made of FGMs was conducted by Allahkarami et al. [32].The shell was embedded in an elastic foundation and subjected to an excitation compressive load.The study utilized FSDT to model the behaviour of the shell under these conditions.Twinkle and Pitchaimani [33] presented the buckling and vibration characteristics of cylindrical shells reinforced with CGP graphene platelets.The effective mechanical properties of open-cell metal foams were evaluated using a modified Halpin-Tsai micro-mechanics model, in combination with an extended rule of mixture.A finite element model was developed by Zghal and Dammak [34] to demonstrate the impact of porosity parameter on the critical buckling responses of powerbased (P-FGM) and sigmoid (S-FGM) plates and spherical caps in FGMs.This model was developed in the context of FSDT.The researchers utilized this approach to investigate the relationship between porosity parameter and buckling behaviour.Thai et al. [35] presented an analytical method for examining the nonlinear instability of eccentrically stiffened FG sandwich truncated conical shells with porosity in the presence of Pasternak elastic foundations.The approach employed is based on displacement approach.A buckling analysis on a cylindrical shell made of porous nanocomposite reinforced with graphene platelet was conducted by Shahgholian et al. [36] based on the FSDT.The researchers computed various mechanical properties, including the effective modulus of elasticity in the thickness direction, by applying the modified Halpin-Tsai micro-mechanics approach.Additionally, they utilized the rule of mixture to determine density and Poisson's ratio.Tao and Dai [37] carried out a geometrically nonlinear analysis on sandwich FG cylindrical shells with a porous core to investigate their post-buckling behaviour when subjected to central point loads and uniform pressures.They used higher-order shear deformation theories (HSDT) and a non-uniform rational B-spline (NURBS)-based isogeometric analysis (IGA) in combination with the modified arc-length method of Crisfield.The study aimed to obtain load-deflection responses and identify snap-through and snap-back post-buckling behaviours.
The constraints of the classical theories [38] and the FSDT [39-41] stimulated scholars to develop a number of HSDT [42][43][44].Among several kinds of the HSDT, the Reddy higher order shell theory (RHST) [45,46] with five unknowns is the most broadly employed theory in the research of plate and shell structures due to adequate efficiency and straightforwardness without any shear correction coefficient.A concise review on the plate and shell laminate theories have been recently presented by Sobhani Aragh et al. [47].What is more, the shells are commonly strengthened by stiffening components for the purpose of improving load-carrying capacity with a relatively little additional weight.Rings and stringers, which are circumferential and meridional stiffeners, respectively, have widespread usage in a myriad of engineering application [48].When lightweight shells with porosities are subjected to a mechanical loading, they might be buckled despite comparatively increased critical buckling loading as a result of utilizing stiffeners.Therefore, there is a necessity to investigate the buckling behaviour of such structures for multi-purpose and cost-effective design with improved performance.
Studies have shown that the mechanical properties of nanocomposites can be enhanced by integrating carbon-nanotubes (CNTs) into metallic matrices.In the last years, the commonly-used techniques for the fabrication of CNT-metal matrix composites have included spark plasma sintering [49,50], hot isostatic pressing [51], microwave sintering [50], and electrochemical deposition [52].Selective laser melting (SLM) [53] and laser powder bed fusion (LPBF) [54] methods of additive manufacturing (AM) have garnered significant attention due to their ability to produce metallic parts with intricate geometries [55,56].Advanced techniques enable the production of components directly from 3-D CAD models through the selective melting and solidification of thin layers of loose powder using high-energy laser beams in a layer-by-layer manner.Gu et al. [57] utilized a laser-assisted am technique known as SLM to produce nanocomposites composed of Al and CNTs that possess customized microstructures and impressive mechanical properties.Meanwhile, Yu et al. [58] achieved the successful production of AlSi10Mg nanocomposites reinforced with multi-walled carbon nanotubes (MWCNTs) at a mass fraction of 0.5% through the utilization of SLM.In a recent work by Yin et al. [59], CNT/316L stainless steel nanocomposites were fabricated by mixing 316L stainless steel powders with various content ratios of CNT through ball milling followed by LPBF.Analysis of the scanning electron microscopy (SEM) images [60,61] revealed that the resulting nanocomposites exhibited CNT aggregation, which resulted in localized areas of higher CNT concentration than the average volume fraction throughout the material.In recent years, a substantial body of literature has been focused on investigating the effects of nanomaterials, specifically CNTs, on the structural behaviour of engineered systems [62][63][64][65][66][67][68][69][70][71].Sobhani Aragh et al. [72][73][74] conducted a study where they used the Eshelby-Mori-Tanaka (EMT) approach, an effective homogenization technique, to examine how the degree of CNT aggregation within the ceramic matrix phase affects the mechanical response of nanocomposite shells.The interested readers are encouraged to refer to Sobhani Aragh's PhD thesis [75] to find comprehensive details on various homogenization techniques.
Motivated by the lack of research work on mechanical buckling of CGP stiffened cylindrical shells strengthened by agglomerated CNTs, this paper aims to bridge this gap by presenting an effective numerical model to investigate buckling behaviour CGP nanocomposite shells.The main contribution of this work is shed light on the effects of porosities as well as CNTs distributions through the shell's thickness on the stability of the shell.Depending on how eccentrically the stiffeners are positioned in relation to the mid-surface of the shell, two types of stiffeners, i.e. internal and exterior, are considered.The smeared stiffener method is employed under the premise of a tightly spaced configuration of stiffeners.The EMT technique is utilized to develop a continuum model for estimating the effective constitutive law of an elastic isotropic metal matrix containing scattered elastic inhomogeneities in the form of CNTs.The shell under combined lateral pressure and axial compressive load is assumed to be composed of closed-cell metal foams that possess graded internal pores, which exhibit Classic, Symmetric, and Asymmetric distributions of porosity based on a power-law distribution.A scaling correlation for the Young's modulus of closed-cell cellular solids is employed, which is obtained through the application of a variational finite element method (FEM).

Methodology
The geometrical configuration of a CGP open cylindrical shell with mean radius , thickness h and length , is shown in Fig. 2. The three orthogonal displacements of an arbitrary point on the shell middle surface in  1 ,  2 , and  3 directions are represented as  1 ,  2 and  3 , respectively.The CGP nanocomposite shell is reinforced by eccentrically closely spaced circular rings and longitudinal stringers with the width and thickness of circular and longitudinal stiffeners   , ℎ  and   , ℎ  , respectively.It is worth highlighting that the investigated shell is subjected to a concurrent application of lateral pressure and axial compressive load.In the present work, the shell is made of closed-cell porous nanocomposite where the porosity coefficient and total volume fraction of the CNTs are constant throughout the shell, however, they smoothly vary along the thickness direction.

Micro-mechanical modelling of accumulated CNTs
As reported in several SEM images [61], CNTs highly tend to agglomerate in the matrix phase, which can be attributed to their low bending stiffness and high aspect ratio.In this part, the effective mechanical characteristics of the nanocomposite shells are determined using a two-parameter EMT model [72,73].The primary limitation of the original elastic inclusion theory proposed by Eshelby is that it assumes the consideration of a solitary inclusion within a semi-infinite elastic medium that is isotropic and homogeneous.This can be alleviated by developing the Mori-Tanaka approach to consider the presence of several inhomogeneities within in a finite region [77].In other words, the EMT method is primarily based on combining Eshelby's concept of equivalent elastic inclusions [78] with Mori-Tanaka's idea of average stress across the matrix phase [79].Inclusions possessing distinct material properties from the surrounding regions are identified in areas exhibiting a high concentration of CNTs.The arrangement of CNTs within a metal matrix phase is determined by their location relative to the inclusions, either inside or outside.The overall volume fraction of the reinforcing phase consisting of CNTs is divided into two components where    and    denote the volume fraction of CNTs situated in the inclusions and distributed through the metal matrix phase, respectively.The presence of agglomerated CNTs in a metal matrix results in the deterioration of elastic properties within the resulting nanocomposite.This phenomenon can be qualitatively explained by proposing two accumulation parameters, denoted as  and , given by where   denotes the effective volume fraction of the inclusions, and the parameter  is the volume fraction of inclusions in relation to the overall volume fraction V of the representation volume element (RVE).Provided that parameter  is equal to unity, the CNTs are not accumulating in the metal matrix phase.In other words, there is only one spherical inclusion, which practically coincides with the total volume fraction.Moreover, an increase of the parameter  leads to the decline of the accumulation degree of CNTs.Moreover, the parameter  specifies the ratio of CNTs volume fraction situated in the inclusions to the total volume fraction of the CNTs.It is apparent that, when  is equal to one, all CNTs are contained within the inclusions.It is noteworthy to mention that a specific constraint must be satisfied in order for the accumulation of CNTs to occur as Mori and Tanaka [79] have developed an approach to evaluate the effective stiffness tensor.The basic idea is to consider the interaction between the inhomogeneities through their effect on the mean values of the fields in the matrix phase.Benveniste's revisitation [80] provides the expression for the effective elastic tensor, which is as follows [81] where where I denotes the fourth-order unit tensor, S represents the fourthorder Eshelby tensor [78], which is specialized to the inhomogeneities with cylindrical geometry, representative of the straight and long CNTs, A  is a fourth-order tensor referred to as concentration factor, and the brackets denote an average.The subscripts m and c stand for the quantities of the matrix and the CNT reinforcing phase, respectively,   and   denote the volume fractions of the respective phases.
The effective bulk modulus   ( 3 ) and shear modulus   ( 3 ) of spherical inclusions containing randomly arranged transversely isotropic CNTs are calculated using the EMT method as Similarly, the effective bulk modulus   ( 3 ) and shear modulus   ( 3 ) of the metal matrix phase are elaborated as ) where   and   represent the bulk and shear moduli of the metal matrix, respectively, and   ,   ,   , and   are defined by ) where   ,   ,   and   denote the Hill's elastic moduli for the reinforcing phase.Ultimately, the effective bulk modulus and shear modulus of the nanocomposite are obtained by where in which   ( 3 ) represents the Poisson's ratio of the metal matrix phase as In above-mentioned formulations,   ( 3 ) and   ( 3 ) are the volume fraction of the CNTs and the metal phase, respectively, which fulfil the expression of   ( 3 ) +   ( 3 ) = 1.The CNTs volume fraction is defined by a gradual change from the inner to the outer surface of the shell, proposing that where in which   denotes the mass fraction of CNTs,   and   are density of CNTs and matrix phase, respectively.  ( 3 ) is the distribution of CNTs along the thickness direction of the shell.This study considers the various distributions of CNTs through the shell's thickness, as classified by Jalali et al. [82].The profiles examined are Pyramid, Inverted Pyramid, Sandglass, and Uniform.

Modelling of porosity dispersion
To improve the strength-to-weight ratio and energy absorbing capability of the structure, porous materials, such as metal foams, has provided unique potential for a number of engineering applications.In these state-of-the-art materials, embedding pores can significantly decrease the structural weight whereas a loss of structural stiffness is partially compensated by including a certain amount of CNTs.Currently, a new paradigm shift is taking place in designing of porous composite materials from uniform dispersion to continuously graded distribution of porosity.Consequently, tailoring the desired properties of the CGP shell is achievable by employing an appropriate pattern for the distribution of porosity through the shell's thickness, in addition to Eq. ( 20) for the smooth variation of CNTs.This study concentrates on the closed-cell porous materials whose theoretical formulation for deriving the mechanical properties are different from those of opencell materials.Nevertheless, this point has not been carefully notified in some literature [83][84][85].The elastic moduli and density of the shell are estimated as where where ) ]   (25) in which  denotes either the effective Young's modulus or shear modulus of the nanocomposite shell and  represents the effective mass density.Subscripts 0 refers to the material properties of a perfect nanocomposite without any voids.  and   known as porosity coefficient and mass density coefficient, respectively.The function   governs through-the-thickness variation of Young's modulus, which have values that range from 0 to 1.In other words, it dominates the graded porosity distribution in the radial direction of the shell using a power-law distribution.For example, three key material profiles through the radial ( 3a-3c.The classical profile through the radial direction is presented as a case of the power-law distributions by setting  = 1 and  = 0.In Fig. 3a, the function   increases through the thickness from 0 at   3 = −0.5 to 1 at   3 = 0.5 for different values of the parameter   .By appropriate choice of the parameter   the porosity distribution can be regulated such that at   3 close to   3 = −0.5, by rising the parameter   the rate of increase in Young's modulus grows while this trend is reversed at   3 close to   3 = 0.5.By adjusting the values of a, b, and , it is feasible to generate both symmetric and asymmetric profiles, as illustrated in Figs.3b and 3c.Symmetric porosity distribution is achieved by setting  = 1,  = 1, and  = 2 in Eq. ( 23).As can be seen in Fig. 3b, the function   on the inner surface is the same as that on the outer surface (  = 1), however, the Young's modulus drops with an increase in the parameter   and tends to   = 0 for higher values of parameter   .Fig. 3c demonstrates asymmetric profiles obtained by setting  = 1,  = 0.5, and  = 2.This figure shows a specific case of the power-law distribution obtained by modifying the parameters , , and .This profile is characterized by the fact that on the inner surface of shell we have a fixed value of   = 1 while on the outer surface the value of the function   is restrained by the parameter   .Moreover, the rate of growth in Young's modulus grows by increasing the parameter   at   3 near the inner surface.
A statistical closed-cell model based on Gaussian random fields (GRFs), which generates a substantial variation in cell sizes and shapes, has been established, as shown in Fig. 4 [86].To predict the Young's modulus, a variational FEM was used.Based on the results obtained, a scaling relation can describe the Young's modulus of the model as which is employed to obtain mass density coefficient,   , as Note that   = 0 corresponds to a case when no pore exists in the nanocomposite while   = 1 is physically meaningless.With an increase in   , density and size of the pores in the material rise and, as a result, modulus of elasticity and mass density of the nanocomposite reduce.Furthermore, Poisson's ratio of the CGP shell is expressed as where  0 = 1 − ( 3 )∕ 0 .

Fundamental formulations
Under the framework of the RHST, the displacements of a generic point at distance  3 from the middle surface of the shell elaborated by ( 1 ,  2 ,  3 ) are related to the middle surface displacement by ) ) where  1 and  2 are rotations at mid-surface of the shell with respect to  1 and  2 .The nonlinear kinematic relations of the shell are defined by (1)   11  (1)   22 11  (3)   22 13 where 11 (1)   11  (1)   22  (1)   12 13 13  (1)   23 where  (33) in which   are the elastic coefficients.
In most literature, the stiffener impact has been smeared over the stiffener spacing, which implies that the impact of discrete ring and stringers have not been accounted.Motivated by this, this paper considers ring and stringers as Euler-Bernoulli beams.Moreover, taking the Kirchhoff-Love assumptions into account, it is assumed that the straightness of the normal to the ring and stringers is preserved even after deformation.Consequently, the displacement fields of the stringer and ring are given, respectively, by where   1 ,   2 , and   3 represent the components of displacement over the thickness of the shell at the th nodal position.Hence, the kinematic relations for stiffeners are obtained ) where  1 and  2 denote the axial and circumferential strains.Consequently, the strain energy for circumferential and axial stiffeners are given by: where   and   represent ring and stringer spacing, respectively, and  denotes shell angle.E is the effective Young's modulus, and the subscripts r and s stand for ring and stringer stiffeners, respectively.Based on the stationary potential energy criterion, the equilibrium equations of the shell under mechanical loadings can be derived.To this end, first, the expression for the overall potential energy of the shell, which is strengthened by ring and stringer stiffeners, is provided as follows  =  +   +   +   (39) where  denotes the potential energy of the mechanical loads, which is defined as the negative value of the work performed by mechanical loading during the deformation of the shell.Considering axial compressive edge load  0 1 and lateral pressure  0 2 , the potential energy of the loading is expressed by It is worthwhile noting that the impact of quadratic terms in Eq. ( 40) on the buckling of short-or medium-length shells is trivial.In particular, the reduction in length resulting from circumferential displacement can be disregarded for stiffened shells having a length-to-radius ratio ranging from 0.7 to 3 [87].Furthermore,   represents the total strain energy of the shell provided by The equilibrium equations of shell may be derived by employing the Euler equations.It can be stated that a loaded shell attains equilibrium when its total potential energy, denoted by , remains constant.This state of equilibrium is achieved when the integrand in the expression for  follows the Euler equations of the calculus of variations.The equilibrium equations of the shell are given in Appendix A.
The equations governing the stability of shells are derived by applying the principle of minimum potential energy.The transition from stable to neutral equilibrium occurs when the total potential energy expression no longer represents a relative minimum.The criterion for stability loss is met when the integrand in the second variation expression of  satisfies the Euler equations of calculus of variations.This procedure will be explained in detail in the following.Employing Taylor series, the expansion of  about the equilibrium state is given by The first variation  is related to the equilibrium state.The stability condition of the original configuration of the shell in the neighbourhood of the equilibrium state may be obtained by the sign of second variation  2 .The condition  2  = 0 can be applied to determine the stability equations of the buckling problem.To obtain the corresponding expression for the second variation of the total potential energy, let where  1 0 ,  2 0 ,  3 0 ,  1 0 , and  2 0 denote the displacement fields of the equilibrium state. 1 1 ,  2 1 ,  3 1 ,  1 1 , and  2 1 are a virtual increment in the displacement fields, i.e. at the neighbourhood of the stable equilibrium.The stress resultants are comprised of two components that express the stable equilibrium and the neighbouring condition.Additionally, the critical buckling is determined by mechanical loading applied to the original shell configuration, provided that the variation equation  2  = 0 is met.By calculating the second variation of the potential and utilizing the Euler equations, the stability equations can be derived, as presented in Appendix B.

Numerical solution
Since the stability equations of the shell, Eqs.(A.1)-(A.6),cannot be solved analytically, consequently, an efficient combined numerical technique is employed to determine the critical buckling load.To this end, first, the Levy method is employed to seek a solution that satisfies the simply supported boundary conditions, and thereby reduces the 2-D problem to a 1-D problem with respect to the coordinate  1 .Afterwards, the generalized differential quadrature (GDQ) technique is used to discretize the governing equations.It has been deemed that the forces  0 in  1 and  2 directions, respectively, may are given by With consideration of the simply supported boundary conditions at the two edges of the shells as: In the Levy procedure, the following representation of the displacement are assumed where   = , and m =1, 2, … denotes the number of half waves in  2 -direction.The partial differential equations (PDEs) obtained by substituting Eq. ( 46) into Eqs.(A.1)-(A.6)are discretized using the GDQ method.The comprehensive overview of the GDQ approach can be found in [75,89].In compact vector notation, the discretized form of PDEs obtained together with related boundary conditions elaborated in Eq. ( 45) is expressed as where and where In order to convert Eq. ( 47) into the standard eigenvalue equation, it is necessary to eliminate the vector {   } .The vector {   } can be obtained from Eq. (47) using the following expression By back-substituting Eq. ( 51) into Eq.( 47), the following expression is derived Eq. ( 52) can be transformed into an eigenvalue problem in the following manner Eventually, the critical buckling load, denoted as F, is determined by solving Eq. ( 53).

Validation
Due to lack of appropriate results for mechanical buckling of the shell for direct comparison, validation of the present methodology is carried out in two ways.Initially, in the context of continuously graded closed cylindrical shells, the precision of the current numerical work utilizing two-parameter EMT is evaluated against that presented in Ref. [88] which employed an extended rule of mixture and graded aligned, straight CNTs.In Table 1, Poly (methyl methacrylate), named as PMMA, has been selected for the matrix phase while the (10, 10) SWCNTs have been considered as reinforcements.Moreover, at a proposed temperature of  = 300 K, which is equivalent to room temperature, there are no thermal strains.As can be seen in this table, there is an excellent agreement between the results verifying the accuracy of the present numerical work.It is interesting to note that for two types of CNTs profiles, the critical mechanical buckling load determined using the EMT is higher than that determined using the extended rule of mixture for various CNTs volume fractions.Secondly, the results are contrasted with those obtained from isotropic closed cylindrical shells.This study centres on the critical buckling load of closed cylindrical shells made of alumina-based isotropic material, both stiffened and unstiffened, as discussed in Ref. [90].As shown in Table 2, the results presented in this comparison for different values of mid-radius to thickness are in good agreement.

Mechanical buckling analysis of porous nanocomposite shells
Having assessed the validity of the current methodology, a number of test cases including different class of structures, ranging from closed to open CGP shells, with or without exterior and interior stiffening elements, different geometrical parameters, and various types of material profiles, are presented.The investigated shell experiences a concurrent application of lateral pressure and axial compressive load.The methodology developed here, indeed, has a broad variety of applications and may be effectively used to do exploratory research on many sorts of design.Aside from showcasing the potential application to a range of design possibilities, this section intends to presenting the high efficiency and fast rate of convergence of the method by presenting the results to establish its very high accuracy and versatility.
For all findings discussed in this section, the buckling load intensity factor (BLIF) is represented by where denotes the rigidity modulus of the shell.The (10,10) SWCNTs as the reinforcement phase and copper as the metal matrix [91,92] are selected.The material properties of metal matrix (copper) are   = 130 GPa and   = 0.34.According to the molecular dynamics simulation results predicted by Shen and Zhang [93], the material properties of the (10,10) SWCNTs are   11 = 5.6466 TPa,   22 = 7.06 TPa,   12 = 1.9445TPa, and   22 = 0.175 at room temperature (300 K) and the effective wall thickness is 0.067 nm satisfying the Vodenitcharova-Zhang criterion [94].All the computations reported in this section are carried out by considering  *  = 0.28 unless the otherwise mentioned.The purpose of this part is to demonstrate the precision of the approach via comparison with both benchmark data and FEM computations.In this work, ABAQUS FE code version 6.17 [95] is used.The FE models use the first-order shell element S4R as a general-purpose shell element with six degrees of freedom per node (three displacements and three rotations).The mesh density is determined through a preliminary convergence analysis, with an average component size of 5 mm.It is of significance to mention that FE analyses have been conducted employing both the shell element S4R and the solid continuum element C3D8R.However, for the purpose of conciseness, the outcomes of the continuum element have not been included in this work.It is crucial to highlight that upon comparing the results of both element types, the S4R element exhibits superior computational efficiency, but the differences in outcomes are minimal.This determination is underpinned by the understanding that using solid continuum elements becomes progressively uneconomical when the R/h ratio surpasses 25 [96,97].Table 3 presents a summary of the results obtained for the CGP shells using the methodology proposed in this study and FEM.In contrast to FEM, the discrepancies observed in the comparison remain negligible, with a maximum deviation of approximately 1%.Figs. 5 and 6 show different buckling mode shapes of the CGP shells with different number of shell angle and stiffeners, respectively.
To attain an accurate BLIF of the CGP, a series of calculations are conducted to ascertain the necessary quantity of grid points along the longitudinal axis of the shell.A convergence study of the numerical results for CGP shells with Sandglass CNTs profile is shown in Figs.7  and 8 for different values of the mid-radius to thickness ratio, ∕ℎ, and the aggregation parameter, , respectively.The results are obtained by setting  = 1,  = 1,  = 2, and   = 1.As observed in these figures, not only the numerical method converges very fast as the number of grid points in the longitudinal direction rises, but it also shows that using only 9 number of grid points one can gain accurate BLIF for the shell.It is worthwhile noting that with the escalation of the aggregation parameter, , resulting in a reduction of the CNTs agglomeration within the metal matrix, there is an observable enhancement in the convergence behaviour.
Figs. 9a and 9b show the variation of the BLIF of the CGP shell with the length-to-mean radius ratio, L/R, different types of CNTs distribution profiles, and porosity dispersion.It can be observed from these figures that, on the one hand, the BLIF is significantly influenced in that the BLIF increases steadily as L/R ratio becomes larger, manifesting that a longer shell possesses a higher critical buckling load.On the other hand, it is interesting to note that the through-thickness variation of CNTs distribution according to the Sandglass type leads to a higher BLIF among other CNTs dispersion types.Comparing Symmetric and Asymmetric porosity distribution in Figs.9a and 9b, respectively, shows that the discrepancy between the Uniform and Pyramid types of CNTs dispersion becomes more tangible for lower values of the L/R ratio.Moreover, the type of CNTs distribution within the matrix through the thickness of the shell has a noticeable impact of the BLIF in higher values of the L/R ratio, irrespective of how the porosity is distributed through the shell's radial direction.
The effect of the parameter   on the buckling behaviour of the CGP stiffened shell is compared in Figs.10a and 10b for different porosity distribution profiles.The porosity distribution governed by  Eq. ( 23) can be influenced through the appropriate choice of the parameter   .It has been explained in Ref. [74] that the tendency of CNTs to aggregate more in the inclusions has a detrimental impact on the mechanical buckling behaviour of cylindrical shells.This trend has been also observed in this figure.Comparing Figs.10a and 10b reveals that for different agglomeration degrees of CNTs within the metal matrix, the porosity distributions symmetrically with respect to the mid-surface of the shell has a remarkable impact on the buckling behaviour of the CGP shell.However, when it comes to Classic porosity distribution, the parameter   has a negligible effect on the BLIF of the CGP shell.Another point that can stand out from Fig. 10b is that the change in the parameter   and, as a result, the symmetric distribution of the porosity is much more appreciable on the buckling behaviour of the CGP shell when the agglomeration degree of the CNTs is more sever.This will be more investigated in the next tables.Effect of number of ring and stringer stiffeners on the buckling load intensity factor of CGP shells with varying the parameter   and the length-to-mean radius ratio for Classic and Symmetric porosity Profiles is investigated in Tables 4 and 5. Certain significant findings evident from these tables are as follows: (i) For various types of porosity distribution profiles, the buckling load intensity factor rises with the number of stiffeners.(ii) By increasing the parameter   for Classic porosity profile, the BLIF slightly drops, whilst for Symmetric porosity profile it is the other way round.This is attributed to different distributions of porosities and, as a result, the Young's modulus of the shell in relation to the mid-surface in Classic and Symmetric Profiles.(iii) In comparison to shells stiffened by interior stiffeners, those stiffened by exterior stiffeners have greater critical buckling loads.
Fig. 11 represents the comparison of the buckling behaviour of continuously graded Symmetric, Asymmetric, and Classic porous shells with three different CNTs distributions.It can be inferred from this figure that the impact of different CNTs distribution profiles on the BILF of the stiffened shell is noticeably higher than those of various porosity dispersion profiles.This highlights the proper selection of the grade CNTs distributions through the shell's thickness.Furthermore, as can be seen from this figure, symmetrical distribution of porosities leads to an improvement in the buckling behaviour of the shell, followed by the Asymmetric profile.In other words, porosities distributed close to the mid-surface are more effective in increasing the stiffness of the structure than other non-uniform porosities distributions, Classic or Asymmetric profiles.Therefore, designers can obtain desired stiffness Table 4 Effect of number of ring and stringer stiffeners on the buckling load intensity factor of CGP shells with varying the parameter   and the length-to-mean radius ratio for Classic porosity Profile.(∕ = 2,  = 2∕3,  =  = 0.5,  = 1,  = 0).a i = r or s .

Table 5
The buckling load intensity factor of CGP shells with varying number of ring and stringer stiffeners, the parameter   and buckling behaviour of CGP shells by regulating distributions of CNTs as well as porosities through the thickness of structure.
In order to scrutinize the impact of the agglomeration parameter  used in the EMT approach on the buckling behaviour of CGP stiffened shells, variation of the BLIF with the parameter  for different porosities and CNTs distribution profiles is shown in Fig. 12.As explained in Eq. ( 2), by increasing the parameter  the agglomeration degree of CNTs within the metal matrix drops and, as a result, the BLIF of the shell gradually grows, followed by higher growth rate for greater values of the parameter .It is interesting to note that the appropriate choice  of CNTs and porosities distribution profiles through the thickness plays a crucial role in the buckling behaviour of the shell with the lower degree of CNTs agglomeration.The BLIF of the CGP shell with Symmetric and Asymmetric porosity profiles increases by 32% and 18%, respectively, than that of Classic profile.
Fig. 13 shows the variation of the BLIF with the shell angle for two values of the mid-radius to thickness ratio, ∕ℎ = 10 and ∕ℎ = 50 and different porosities distribution profiles.As can be seen in this figure, overall, the BLIF of both R/h ratios declines by increasing the shell angle irrespective of the type of porosity profile.However, the BLIF of shells with lower values of the R/h ratio slumps as  < 100, whilst the downturn for thicker shells follows a mild trend and then remains almost unaltered for closed shells.The impacts of the agglomeration parameter  and CNTs profiles on the variation of the BLIF with the shell angle are studied in Fig. 14.It is found that the discrepancy in the buckling behaviour of CGP shells with different CNTs profiles substantially increases by reducing the agglomeration degree of CNTs dispersed within the metal matrix phase.This trend is more appreciable in smaller values of the shell angle.

Conclusion
This study has introduced a numerical model to investigate the buckling behaviour of bio-inspired CGP nanocomposite cylindrical shells.These shells were constructed using metal foams with closedcell structures characterized by power-law distributed porosity profiles.By incorporating an elastic isotropic metal matrix containing dispersed elastic CNTs, the study considered the impact of CNTs agglomeration using the EMT approach to derive the effective constitutive law.The focus was on understanding the effects of geometrical parameters, porosity variability, and CNTs distribution on the buckling performance of CGP shells, thereby offering insights for designing enhanced stability in closed-cell cellular stiffened shells.Through the application of the adjacent equilibrium criterion, the stability equations were solved  employing the GDQ and Levy techniques.The results have contributed substantially to the field by providing fundamental knowledge that serves as a basis for developing cellular stiffened shells exhibiting optimal porosity and CNTs distributions, consequently leading to improved stability.
The primary findings underscored the pivotal role of CNTs distribution within the matrix on the enhancement of the BLIF for higher L/R ratios.Furthermore, the symmetric distribution of porosity around the mid-surface significantly influenced the buckling behaviour, particularly when dealing with varying degrees of CNTs agglomeration.Notably, the controlling porosity parameter (  ) displayed minimal Fig. 13.BLIF versus shell angle for different values of mid-radius to thickness ratio and porosities distribution profiles.As evident from the this figure, the BLIF for both R/h ratios decreases with an increasing shell angle, regardless of the porosity profile.However, for shells with lower R/h ratios and  < 100, the BLIF experiences a significant decline.In contrast, the BLIF for thicker shells exhibits a gradual decline, stabilizing for closed shells.(∕ = 2,  = 0.2,  = 0.7,   =   = 10).impact on buckling behaviour when employing a Classic porosity distribution.In addition, the strategic allocation of porosities closer to the mid-surface, in contrast to alternative non-uniform distributions, effectively heightened the stiffness of the structure.The study further disclosed that the deviation in buckling characteristics became more pronounced as the degree of CNTs agglomeration decreased within the metallic matrix phase, a trend more noticeable for lower shell angles.This work has contributed valuable insights towards the design of CGP shells with amplified stability and structural resilience.The findings laid out in this study set the stage for future endeavours in optimizing lightweight and robust composite structures that harness the synergistic potential of both CNTs and porosity distributions.))] = 0 where N 1  1 and N 2  2 are the mechanical forces applied to the ring and stringer stiffeners expressed by:

Fig. 1 .
Fig. 1.(a) Prototype a 25-millimeter-thick layer made of porous foam.(b) Highly porous aluminum alloy foam used in the fabrication of the frontal structure of the Intercity-Express-Train (ICE) train, jointly developed by Voith Engineering and the Fraunhofer Institute for Machine Tools and Forming Technology (IWU).[21,22].

Fig. 2 .
Fig. 2. (a) Schematic configuration of a CGP shell strengthened by rings and stringers.(b) Cross-sectional image of the shell showing the varying pore size along the thickness direction [76].(c) Visual representation illustrating the arrangement of CNTs within inclusions and their dispersion throughout the metal matrix phase at the meso-scale.

Fig. 7 .
Fig. 7. Convergence of the buckling load intensity factor versus the mid-radius to thickness ratio, ∕ℎ, of CGP shells with Sandglass CNTs profile.As evident from these figures, the numerical approach not only exhibits rapid convergence as the longitudinal grid points increase in number but also underscores the possibility of achieving accurate BLIF for the shell using a minimal count of 9 grid points.(∕ = 2,  = 2∕ = 2,   =   = 10,  = 1,  = 1,  = 2,   = 1).

Fig. 8 .
Fig. 8. Convergence of the buckling load intensity factor versus the aggregation parameter, , of CGP shells with Sandglass CNTs profile.As depicted in this figure, with the increase of the aggregation parameter denoted as , leading to a reduction in the agglomeration of CNTs within the metal matrix, there is an observable improvement in the convergence behaviour.(∕ = 2, ∕ℎ = 20,  = ∕2,   =   = 15,  = 1,  = 1,  = 2,   = 1).

Fig. 9 .
Fig. 9. Variation of the buckling load intensity factor of CGP shells with the length-to-mean radius ratio, L/R, for different types of CNTs distribution profiles.As discerned from the figure, the manner in which CNTs are distributed within the matrix across the shell's thickness exerts a notable influence on the BLIF for elevated values of the L/R ratio, independent of the mode of porosity distribution along the radial direction of the shell.(∕ℎ = 10,  = 2,   =   = 20,  =  = 0.5,  = 2,   = 1).(a) Symmetric profile ( = 1,  = 1).(b) Asymmetric profile ( = 1,  = 0.5).

Fig. 11 .
Fig. 11.Comparison of the buckling behaviour of continuously graded Symmetric, Asymmetric, and Classic porous shells with three different CNTs distributions.It can be discerned from this figure that the influence of distinct CNTs distribution profiles on the BLIF of the stiffened shell is markedly more pronounced compared to the effects stemming from diverse porosity dispersion profiles.(∕ = 2,  = 2∕3,   =   = 10,   = 1,  =  = 0.5).

Fig. 12 .
Fig. 12. Variation of BLIF with agglomeration parameter  for different porosities and CNTs distribution profiles.The appropriate selection of distribution profiles for CNTs and porosities throughout the thickness assumes a pivotal role in influencing the buckling behaviour of the shell, particularly when dealing with a lower degree of CNTs agglomeration.In comparison to the Classic profile, the BLIF of the CGP shell experiences a notable increase of 32% with the adoption of the Symmetric porosity profile and an 18% increase with the utilization of the Asymmetric porosity profile.(∕ℎ = 20, ∕ = 2,  = ∕2,   =   = 15,   = 1,  = 0.95).

Fig. 14 .
Fig. 14.Variation of BLIF with shell angle for different values of agglomeration parameter  and CNTs profiles.It can be observed from this figure that the disparity in the buckling behaviour of CGP shells with varying CNTs profiles is notably amplified by diminishing the degree of CNTs agglomeration dispersed within the metal matrix phase.(∕ℎ = 40, ∕ = 2,  = 0.2,   =   = 20,  = 1,  = 1,  = 2).

Table 2
Comparison of the critical buckling loads of unstiffened and stiffened isotropic closed cylindrical shell ( =  = 0.3).