Post-Buckling Analysis of FG-TPMS Shells with Geometric Imperfection and Porosity under Axial Compression

1  Introduction

Light-weight structures, such as shells, play a crucial role in various fields including civil, mechanical, aerospace, and marine engineering. Real-world applications of shells include roofs of football stadiums, factories, railway stations, airports, and submarines, ships, spaceships, and airplanes. In civil engineering, using such shell structures can save materials and costs because of their advantages, such as light-weight and long-span. Recently, bio-inspired structures, such as porous, lattice, and honeycomb structures, have been widely utilized because of their remarkable mechanical properties [1,2]. FG-TPMS structures, which are bio-inspired structures and contain internal pores, have attracted considerable attention from scientists and researchers. There are four types of FG-TPMS architectures: sheet-based, thickened-based, skeletal-based, and network-based structures [3]. The sheet-based FG-TPMS structures investigated in this work possess some remarkable mechanical properties compared with the other types of structures [4,5]. Mathematical formulas can be utilized to express the geometries of FG-TPMS architectures, which possess remarkable advantages including many internal pores, high inter-connectivity, smooth surfaces, and mathematically controllable geometries. FG-TPMS structures can be fabricated numerously and rapidly by using 3D printing technology [6]. FGTPMS-based structures are therefore used in many applications, such as bone implants, batteries, microwave absorption, semipermeable membranes, energy absorbers, and engine brackets [7]. Thus, investigations on the responses of FG-TPMS structures are very necessary and practical. Modeling FG-TPMS structures is challenging and time-consuming because of the complicated geometries and responses. To overcome this difficulty, an effective approach for modeling FG-TPMS structures using equivalent single-layer (ESL) theories was proposed in [8]. The numerical implementation of the approach is accurate, efficient, and simple. Therefore, the approach was extended and investigated for analyses of FG-TPMS plates [9], FG-TPMS plates with graphene platelet reinforcements [10], and FG-TPMS shells [11]. In addition, the free vibration analysis of functionally graded TPMS porous quadrilateral plates reinforced with graphene nanoplatelets was performed in [12], and the wave propagation behaviors of triply periodic minimal surface sandwich shells were investigated using numerical and analytical approaches were performed in [13,14]. In structural engineering, geometric and material imperfections can significantly reduce the load-carrying capacity of structures and lead to unpredictable failures. These imperfections can stem from inaccurate fabrication and erection, and should be taken into account in the analysis and design. For the first time in this paper, a post-buckling analysis of FG-TPMS shells with geometric imperfection and porosity under axial compression is performed using the ESL theory-based approach. The formulation was established by using IGA and FSDT.

There are some theoretical approaches for modeling shells such as unified formulations, ESL theories, and 3D elastic theory [15]. The formulations using ESL theories achieved a balance between the accuracy of solutions and the computational cost, most shells are consequently analyzed using these formulations. We classify thin, thick, or moderately thick shells using their thicknesses, and shell analyses should be performed using the appropriate theories to obtain accurate results. The classical shell theory (CST), which eliminates transverse shear deformations, is appropriate for the analysis of thin shells. For analyses of thicker shells, we can use the FSDT to obtain accurate solutions. Then, higher-order shear deformation shell theories (HSDTs) were proposed to analyze thin, moderately thick, and thick shells and overcome the disadvantages of the CST and FSDT. Several HSDTs can be mentioned, such as the third-order shear deformation theory [16,17], quasi-3D theory [18], the generalized theories [19,20]. In this paper, the formulation is established by using the FSDT and IGA, which is considered a powerful and accurate method. The IGA was proposed by using the core idea of the finite element method (FEM) and NURBS basis functions, which possess higher order derivatives [21]. The NURBS basis functions are utilized in the IGA to approximate solutions and to express exact geometries. One of the benefits of using the IGA is that the geometrical data of structures, which are extracted from computer-aided design (CAD), can be directly used for simulations. As an advantage, the solution error caused by inaccurate mesh generation is completely eliminated in the IGA. Therefore, we can save time on mesh generation, obtain highly accurate solutions, and solve problems with complex geometries quickly by using the IGA. The IGA can be considered a bridge connecting the FEM and CAD. Accordingly, the IGA has great potential for simulations in many engineering fields. A comprehensive review of the IGA can be found in [22]. The IGA was investigated for the analyses of curved structures, such as membrane structures [23], thin shells [24–26], thin and thick shells [27,28], shell optimization [29], and nonlinear analysis of shells [30–33]. Notably, the application of IGA to analyses of FG-TPMS shells is very limited in the literature, except for the linear analyses of FG-TPMS shells in [11]. For the first time, post-buckling behavior of FG-TPMS panels subjected to compression utilizing IGA and FSDT is investigated in this paper. This is the first work which studies on effects of geometric and material imperfections on the nonlinear behavior of FG-TPMS shells. In this paper, the non-linearity of shells is modeled using the von Karman assumption and the modified Riks method is used to iteratively solve the discrete nonlinear equation system.

2  Post-Buckling Analysis of FG-TPMS Shells with Geometric Imperfection and Porosity under Axial Compression

2.1 Mathematical Model of FG-TPMS Materials Considering Porosity

In this work, FG-TPMS shells with the cell geometry I-graph and wrapped package-graph (IWP) are investigated. The shells consider both geometric and material imperfections as shown in Fig. 1. The geometry of sheet-based TPMS can be mathematically expressed as follows [8]:

ψ(x,y,z)=2(cos⁡(χ1x)cos⁡(χ2y)+cos⁡(χ2y)cos⁡(χ3z)+cos⁡(χ3z)cos⁡(χ1x))−−(cos⁡(2χ1x)+cos⁡(2χ2y)+cos⁡(2χ3z))whereχi=2πnili;i=1,2,3(1)

Figure 1: The geometric and material imperfections.

In which, li is the length of a unit cell in the three directions x, y, and z. In addition, ni is the number of unit cells. The popular level-set equation utilized to produce minimal surface lattices is ψ(x,y,z)=c. In which, c is known as the level-set constant and c is equal to zero for TPMS-like surfaces. In this instance, the obtained minimal surfaces separate the space into some sub-domains with equal volumes [3]. Fig. 2 describes a FG-TPMS panel under axial compression as an illustration. A fitting method [8] was employed to calculate the mechanical properties of the FG-TPMS materials. Thus, the relative density at an arbitrary material point inside the structure is computed as ρ=VVs [8]. In this formula, V and Vs are respectively known as volumes of a TPMS cell and its surrounding cube. Density distribution patterns can be classified into two groups as the following:

ρ=ρmin+(ρmax−ρmin)×(zh+12)n Pattern Aρ=ρmin+(ρmax−ρmin)×(1−cos⁡(πzh))nPattern B(2)

Figure 2: A FG-TPMS panel under axial compression.

Based on the values of ρmin, ρmax, and n, the two aforementioned groups (A and B) can be separated into six density distribution patterns, namely, A1, A2, A3, B1, B2, and B3 as given in Table 1. Fig. 3 illustrates a FG-TPMS plate in the case the average relative density ρ=0.35. Besides, illustrations for six density distributions across the thickness are shown in Fig. 4. The mechanical properties of FG-TPMS materials, such as Poisson’s ratio ν, shear modulus G, and elastic modulus E at a material point, can be determined by the relative density ρ at that point as follows [8]:

E=0.597Esρ1.225;ρ≤0.35E=(0.987ρ1.782+0.013)Es;ρ>0.35G=0.529Gsρ1.287;ρ≤0.35G=(0.960ρ2.188+0.040)Gs;ρ>0.35ν=2.597e−0.157ρ−2.244;ρ≤0.13ν=0.201ρ2−0.227ρ+0.326;ρ>0.13(3)

in which, Es and Gs are the mechanical properties of the parent or base material.

Figure 3: A FG-TPMS plate using the unit cell I wrapped package-graph (IWP).

Figure 4: Six density (ρ) distributions cross the shell thickness.

In this work, we establish formulas for FG-TPMS materials considering porosity based on that for functionally graded materials considering porosity. Firstly, the effective properties of functionally graded materials with evenly distributed porosity are expressed as follows [34]:

Pe(z)=(Pc−Pm)(12+zh)n+Pm−α2(Pc+Pm)(4)

The effective properties of functionally graded materials with unevenly distributed porosity are expressed as follows [34]:

Pe(z)=(Pc−Pm)(12+zh)n+Pm−α2(Pc+Pm)(1−2|z|h)(5)

where Pc and Pm are the material properties (E, G, ν) of ceramic and metal, respectively, while n is the power index. α (0≤α≪1) is the porosity volume fraction. The parameter α is equal to zero in the case of no porosity, which is called perfect functionally graded materials. In this work, the present structures use one phase of material, i.e., the parent or base material. Thus, formulas for FG-TPMS materials considering porosity are obtained by simplifying the formulas for functionally graded materials considering porosity (Pc=Pm=P) as follows: The effective properties of FG-TPMS materials with evenly distributed porosity are as follows:

Pe(z)=P(1−α)(6)

The effective properties of FG-TPMS materials with unevenly distributed porosity are obtained as follows:

Pe(z)=P[1−α(1−2|z|h)](7)

Notably, P and Pe(z) are, respectively, the mechanical and effective properties of FG-TPMS materials considering porosity, such as Poisson’s ratio ν, shear modulus G, and elastic modulus E at a material point.

2.2 Post-Buckling Behavior of FG-TPMS Shells Subjected to Compression

Firstly, a shell shown in Fig. 2a is studied. An orthogonal curvilinear coordinate system is introduced as follows: The x-axis and y-axis match the lines of the principal curvatures of the mid-surface, while the z-axis is perpendicular to the mid-surface. Then, physical quantities, such as displacements, stresses, and strains, are expressed based on this orthogonal curvilinear coordinate system. The FSDT displacement field is written as follows [35]:

u(x,y,z)=u0(x,y)+zβx(x,y)v(x,y,z)=v0(x,y)+zβy(x,y)oru¯=u0+zu1w(x,y,z)=w0(x,y)(8)

where

u¯={uvw};u0={u0v0w0};u1={βxβy0}(9)

u0, v0, and w0 respectively are the tangential and radial displacements of a point on the mid-surface. βx and βy respectively denote the rotations in the y-axis and the x-axis. Note that a singly curved shell is obtained if R1=R and R2=∞. We mainly study singly curved shallow shells subjected to compression in this work, the following assumptions can be accepted as: constant radii of curvatures ax,y=ay,x=0 and (1+z/R)≈1. The detailed formulas for nonlinear engineering strain components based on the von Karman assumption in the orthogonal curvilinear coordinate system can be found in [36]. The Sanders’ nonlinear strains are written as [16]

ε={εxxεyyγxy}T=ε0+zκbγ={γxzγyz}T=εs(10)

the membrane strain ε0 taking into account the initial deformation wI is expressed as follows:

ε0=εL+εN+εI;εL={u0,x+w0Rv0,yu0,y+v0,x};εN=12{w0,x2w0,y22w0,xy};εI={w0,xwI,xw0,ywI,yw0,xwI,y+w0,ywI,x};κb={βx,xβy,yβx,y+βy,x};εs={−u0R+w0,x+βxw0,y+βy}(11)

the nonlinear and initial strains are re-expressed as follows:

εN=12Aθθ;Aθ=[w0,x00w0,yw0,yw0,x];θ={w0,xw0,y}εI=12AIθ;  AI=[2wI,x002wI,y2wI,y2wI,x](12)

The virtual work equation based on the total Lagrangian approach is expressed as

∫Ωσ^Tδε^dΩ=∫Ωδu¯TfsdΩ(13)

Ω is the initial configuration of the shell, while fs={fxfyfz}T denotes the external load vector. For the post-buckling analysis, the virtual work equation of the shell subjected to axial compression λNy shown in Fig. 2 is written as follows:

∫Ωσ^Tδε^dΩ=λ∫Ωδv0NydΩ(14)

λ denotes the load factor, and σ^ represents the resultant stress vector as follows:

σ^={σpσbσs}T(15)

which contains the following components in plane,

σp={NxxNyyNxy}T={∫−h/2h/2(σxσyτxy)dz}T(16)

bending,

σb={MxMyMxy}T={∫−h/2h/2(σxσyτxy)zdz}T(17)

shear,

σs={QxQy}T={∫−h/2h/2(τxzτyz)dz}T(18)

The stress resultant vector σ^ is expressed in terms of the generalized strain vector ε^ and takes into account the initial deformation as follows:

σ^=D^ε^;D^=[AB0BDb000Ds];ε^={εLκbεs}+{εN+εI00}(19)

where

(Aij,Bij,Dijb)=∫−h/2h/2(1,z,z2)Qijdz;i,j=1,2,6Dijs=κ∫−h/2h/2Qijdz;i,j=4,5(20)

and

Q11=Q22=E1−υ2;Q12=Q21=υE1−υ2;Q44=Q55=Q66=G(21)

the shear correction factor is κ=5/6 [16].

3  Post-Buckling Analysis of FG-TPMS Shells Using NURBS Basis Functions and FSDT

3.1 NURBS Basis Functions

The NURBS basis functions are briefly presented in this section. Details and corresponding source codes were presented in [22]. A knot vector Ξ={ξ1,ξ2,...,ξn+p+1}, ξi∈R and i=1,...,n+p+1 is first considered. In the case the first and last knots are repeated p+1 times, the vector is considered open. In addition, a B-spline basis function has continuity C∞ inside a knot span and Cp−1 at each knot. The B-spline basis functions in one-dimensional parametric space Ni,p(ξ) can be expressed as follows [37]:

Ni,p(ξ)=ξ−ξiξi+p−ξiNi,p−1(ξ)+ξi+p+1−ξξi+p+1−ξi+1Ni+1,p−1(ξ) asp=0,Ni,0(ξ)={1  ifξi≤ξ<ξi+10otherwise}(22)

The two knot vectors Ξ={ξ1,ξ2,...,ξn+p+1} and H={η1,η2,...,ηm+q+1} are now considered. B-spline basis functions in the two-dimensional parametric space can be achieved through a tensor product of the basis functions in the dimensions ξ and η as

NAb(ξ,η)=Ni,p(ξ)Mj,q(η)(23)

Taking into account weights w for control points, the non-uniform rational B-spline basic functions (NURBS) are utilized to accurately express complex geometries such as ellipsis, spheres, and circles, etc. Therefore, the NURBS functions are written as follows [21]:

NA(ξ,η)=NAbwA∑I=1m×nNIb(ξ,η)wI(24)

in the case all weights of control points are equal, the NURBS and B-spline functions are equivalent.

3.2 Formulation Utilizing the NURBS Basis Functions and FSDT

In this paper, NURBS basis functions are utilized to interpolate the geometries and displacement fields u of shells as follows:

xh(x,y)=∑A=1m×nNA(ξ,η)PAanduh(ξ,η)=∑A=1m×nNA(ξ,η)qA(25)

NA(ξ,η) represents the NURBS basis function. PA is the control point coordinate. The vector qA={u0Av0Aw0AβxAβyA}Tcontains the degrees of freedom. xh(x,y) includes the coordinates of physical points, while uh denotes the interpolated displacement vector.

By substituting Eqs. (25) into (11), the linear, bending, and shear strains can be expressed as follows:

εL=∑A=1nBALqA;κb=∑A=1nBAbqA;εs=∑A=1nBAsqA(26)

where

BAL=[NA,x01RNA000NA,y000NA,yNA,x000];BAb=[000NA,x00000NA,y000NA,yNA,x]BAs=[−1RNA0NA,xNA000NA,y0NA](27)

By substituting Eqs. (25) into (12), the nonlinear and initial strains can be re-written as follows:

θ=∑A=1nBAgqA; BAg=[00NA,x0000NA,y00]; εN+εI=12∑A=1nBANqA; BAN=(Aθ+AI)BAgBAN={[w0,x00w0,yw0,yw0,x]+[2wI,x002wI,y2wI,y2wI,x]}[00NA,x0000NA,y00](28)

The interpolated strains in Eqs. (26) and (28) can be substituted into Eq. (19). Then, the variation of the generalized strain ε^ can be expressed as follows:

δε^=∑A=1nBAδqA;BA={BALBAbBAs}+{BAN00}(29)

All the terms in Eqs. (25) and (29) are substituted into Eq. (14). By removing the virtual displacement, the nonlinear equations can be achieved.

Ψ(q,λ)=∫ΩBATσ^AdΩ−λ∫ΩNyNAdΩ=0;NA={0NA000}T(30)

The equation system in (30) is nonlinear with respect to unknown variables q and λ. Therefore, the modified Riks method, which includes a series of linear solutions, can be utilized to solve the system. If there is a small change in the external load from λmNy to (λm+Δλm)Ny, a new equilibrium configuration around the old one is achieved as

Ψ(qm+Δqm,λm+Δλm)=0(31)

Utilizing the Taylor series expansion for Eq. (31) and only taking into account the first term, the following system of linear incremental equilibrium equations can be achieved:

KT(qm)Δqm=(λm+Δλm)∫ΩNyNAdΩ−∫ΩBATσ^AdΩ(32)

KT is known as the tangent stiffness matrix estimated at qm as:

KT(qm)=∂Ψ∂q|qm(33)

Employing the modified Riks method for post-buckling and nonlinear analyses, an iterative process is carried out for all load levels [38]. Accordingly, Eqs. (32) and (33) are re-written in a generalized form with the mth load step and the ith iteration as follows:

KT(qm)Δqmi=Fext,mi−Fint,mi(34)

where

KT=∫Ω[{BALBAbBAs}+{BAN00}]T[AB0BDb000Ds][{BALBAbBAs}+{BAN00}]dΩ++∫Ω(BAg)T[NxxNxyNxyNyy]BAgdΩ(35)

the external load vector

Fext,mi=(λmi+Δλmi)∫ΩNy{0NA000}TdΩ=(λmi+Δλmi)F0(36)

where F0 denotes the reference load. The internal load is determined as follows:

Fint,mi=∫Ω(Bmi)T(σ^mi)dΩ(37)

Replacing Eqs. (19), (26), (28), and (29) into Eq. (37), we obtain

Fint,mi=Kmiqmi(38)

where

Kmi=∫Ω[{BALBAbBAs}+{BAN00}]T[AB0BDb000Ds][{BALBAbBAs}+0.5{BAN00}]dΩ(39)

The iterative process at each load level is conducted till the convergence criterion for the residual load is matched as follows:

‖λmiF0−Fint,mi‖‖(λmi+Δλmi)F0‖<10−8(40)

Notably, incremental solutions can be obtained by solving Eq. (34). Thereafter, the load factor and displacement vector are determined as follows [38]:

λmi+1=λmi+Δλmiqmi+1=qmi+ΔqmiΔqmi=ΔqR,mi+ΔλmiqF,m(41)

ΔqR,mi and qF,m are the displacement vectors caused by the residual and reference load vectors, respectively.

ΔqR,mi=[KT(qm)]−1(λmiF0−Fint,mi)qF,m=[KT(qm)]−1F0(42)

4  Results and Discussions

Steel-SUS304 was used as the parent (or base) material of FG-TPMS structures for all the problems in this work. Its mechanical characteristics are as follows [8]: Es=200 GPa and υs=0.3. For all the following FG-TPMS structures, the average relative density ρ=0.35 was fixed [8]. The present formulation was verified through several benchmark problems of an isotropic plate considering geometric imperfection, and a perfect panel under axial compression, and perfect FG-TPMS panels under uniform pressure. This is the first study on FG-TPMS panels considering material porosity, so the results of these structures are not available in the literature for reference. Thus, we verified the proposed formulation for FG-TPMS panels without material porosity. Notably, in the present formulation, FG-TPMS panels without material porosity are special cases of that considering material porosity when the porosity volume fraction α=0. All the problems in this work were solved utilizing 14 × 14 cubic NURBS elements. For numerical integration, we used 4 × 4 Gauss points per element. The boundaries of a plate and a shell described in Figs. 2 and 5 are mathematically expressed as follows

•   The fully simply supported boundary for the plate [39]

u0=w0=βx=0at the loaded edgesw0=βy=0at the unloaded edges(43)

•   The simply supported (S) boundary for the panel [39]

w0=βy=0atx=0,αR(44)

•   The clamped (C) boundary for the panel [39]

u0=v0=w0=βx=βy=0aty=0;u0=w0=βx=βy=0aty=L(45)

Figure 5: A simply supported isotropic square plate subjected to the axial compression [39].

In this paper, the normalized parameters and geometrical imperfections for plates and panels are mathematically expressed as follows:

•   Normalized central deflection [39]

w^=w/h:plates, panels(46)

•   Load parameters [40,41]

Ncr=Nya2π2Eh3for isotropic plateNb=NyL2Esh3for panels(47)

•   The initially geometrical imperfection of plate [40]

wI=W0sin⁡(πxa)sin⁡(πya)(48)

•   The initially geometrical imperfection of panels [42]

wI=W0sin⁡(xR)sin⁡(πyL)(49)

where W0 denotes the maximum magnitude of the initial deflection due to geometrical imperfection. R and L respectively are the radius and the length of the panel while a denotes the width of the square plate. Ny denotes the compressive load depicted in Figs. 2 and 5. In this paper, the geometrical imperfections given in Eqs. (48) and (49) were assumed to be initial sinusoidal deformation. In addition, the material imperfections given in Eqs. (6) and (7) were assumed to be evenly or unevenly distributed pores. These assumptions are simplifications of geometrical and material imperfections. As a matter of fact, FG-TPMS structures are usually manufactured by 3D printing technology. The material and geometrical imperfections can be caused by the manufacturing process. Thus, uncertainties in material and geometry need to be considered in future works.

4.1 Verification

4.1.1 Plate with Geometric Imperfection

The square plates are particular circumstances of singly curved shells if R=∞. We first consider a simply supported square plate with a/h=120 and subjected to the compression depicted in Fig. 5. The mechanical characteristics of material are as follows: Poisson’s ratio v=0.3 and Young’s modulus E=300 GPa. The initial deflection due to geometrical imperfection is determined by Eq. (48), where the maximum magnitude is W0=0.1h. The normalized central deflection w^ is determined by Eq. (46) and the load parameter Ncr is defined by Eq. (47). The obtained result using the proposed formulation is presented in Fig. 6 in the comparison with the analytical solution of Yamaki [40]. In [40], post-buckling analysis of functionally graded material sandwich plates subjected to in-plane mechanical compression was performed using a mesh-free method and a higher-order shear deformation plate theory. Then, the formulation was used to analyze an isotropic plate considering the geometric imperfection, as mentioned above. Note that the formulation in this paper was established to analyze the post-buckling behavior of FG-TPMS shells with geometric imperfection and porosity under axial compression using FSDT and IGA. Interestingly, both approaches (in [40] and in this paper) produced accurate results compared with the analytical solution of Yamaki. The present formulation for isotropic plates considering geometric imperfections is verified.

Figure 6: Normalized central deflection and load parameter of the plate.

4.1.2 Isotropic Perfect Panel

The clamped-simply-clamped-simply cylindrical shell (CCSS) subjected to compression described in Fig. 7 is next considered. The material properties are as follows [41]: Poisson’s ratio v=0.33 and Young’s modulus E=10 Msi. The geometric properties are as follows: length L=14.75 in., radius R=60 in., thickness h=0.13 in., and b=14.5 in. As a result, the response of the panel is obtained in Fig. 8, where the compression Pb=Nyb. The present solution is compared with that obtained using the mesh-free kp-Ritz method (MKR) and FSDT [41]. The limit points are successfully detected on the equilibrium path. The two solutions well agree when the radial displacement is smaller than 0.3 in. In the case the radial displacement is greater than 0.3 in., two solutions do not well agree. This discrepancy can be explained as follows: For a relatively large number of load steps, the accumulated geometrical description error of IGA is eliminated by utilizing NURBS but that of MKR method becomes significant. Again, the present formulation for isotropic perfect panels under compression is confirmed.

Figure 7: A clamped-simply-clamped-simply (CCSS) supported panel subjected to the axial compression [39].

Figure 8: The central deflection and compression load of the CCSS isotropic panel.

4.1.3 FG-TPMS Perfect Panels

The present formulation for simulating FG-TPMS materials is verified in this section. However, the results of FG-TPMS shells subjected to compression are not available in the literature. Thus, we verify the proposed formulation for the linear analysis of FG-TPMS panels under uniform pressure. An FG-TPMS panel under uniform pressure q=0.2 MPa is now considered. The panel is fully clamped, and its geometry is described as following: h/R=0.002, α=0.1 rad., R=1 m, and L/R=0.1. Table 2 shows the obtained results using the present formulation in the comparison with that using linear static analysis based on FSDT and IGA [11]. An excellent agreement is found. It is concluded that the present formulation possesses high accuracy in predicting the nonlinear responses of isotropic plates considering geometric imperfection, isotropic perfect panels under axial compression, and in modeling FG-TPMS materials.

4.2 Parametric Study

The CCSS FG-TPMS panel depicted in Fig. 7 is continued investigating in this section. The compressive load is applied with the load step ΔNy=1 kN/m. The panel geometry is provided as R=2540 mm, h=6.35 mm, and L=B=508 mm. For panels, the geometrical imperfection is determined by Eq. (49) [42], and the load parameter Nb is determined by Eq. (47) [41].

4.2.1 FG-TPMS Perfect Panels without Geometric Imperfection and Porosity

FG-TPMS perfect panels without geometric imperfection and porosity under axial compression are studied in this section. Fig. 9 shows the load-deflection curves of panels with six density distribution patterns. As seen, the panels using patterns B (including B1, B2, and B3) have higher post-buckling strengths than those using patterns A (including A1, A2, and A3). We can conclude that the panels with patterns B have higher bending and membrane stiffness than those using patterns A. Hence, we recommend to utilize FG-TPMS panels with patterns B, which possess two advantages as depicted in Fig. 4b, as follows: 1) Density is rich at both top and bottom surfaces; 2) Symmetrical density distributions with respect to the mid-plane. As seen in Fig. 9, the panel with pattern B3 produces the highest post-buckling strength and the highest bending and membrane stiffness compared with the panels utilizing the remaining patterns A1, A2, A3, B1, and B2. Fig. 10 shows the load-deflection curves of panels with various curvatures κ=L/R. As seen the higher the curvature of the FG-TPMS panel, the higher the post-buckling strength.

Figure 9: Load-deflection paths of CCSS FG-TPMS perfect panels subjected to compression, R=2540 mm, L=508 mm, B = L, h=6.35 mm, W0=0, α=0.

Figure 10: Load-deflection paths of CCSS FG-TPMS perfect panels subjected to compression, pattern B3, L=508 mm (fixed), B = L, h=6.35 mm, W0=0, α=0, various curvatures κ=L/R.

4.2.2 FG-TPMS Panels with Geometric Imperfection but without Porosity

FG-TPMS panels with geometric imperfection and without porosity (α=0) under axial compression are studied in this section. Figs. 11 and 12 show the load-deflection curves of panels with various geometric imperfections for six density distribution patterns. As seen, the panels with patterns B2 and B3 have the highest post-buckling strengths. This is because of the advantages of FG-TPMS panels with patterns B, as explained in Section 4.2.1. All panels have similar critical buckling loads for an arbitrary density distribution pattern. After the load-deflection curve passes the second limit point, the larger the geometric imperfection, the lower the post-buckling strength.

Figure 11: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, patterns A1, A2, A3, α=0, various geometric imperfections.

Figure 12: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, patterns B1, B2, B3, α=0, various geometric imperfections.

4.2.3 FG-TPMS Panels without Geometric Imperfection but with Porosity

FG-TPMS panels without geometric imperfection (W0=0) and with porosity under axial compression are studied in this section. Figs. 13 and 14 show the load-deflection curves of panels with various porosity volume fractions for even porosity distribution. In the meanwhile, Figs. 15 and 16 show the load-deflection curves of panels with various porosity volume fractions for uneven porosity distribution. Clearly, the porosity volume fraction significantly affects the post-buckling strength of the FG-TPMS cylindrical shell. The higher the porosity volume fraction, the lower the post-buckling strength of FG-TPMS panel. We can conclude that porosity significantly reduces nonlinear strength of FG-TPMS cylindrical shells. Interestingly, the effect of the porosity volume fraction on post-buckling strength of the FG-TPMS panel with even porosity distribution is more significant than that of the panel with uneven porosity distribution. Especially in the case of uneven porosity distribution, the mentioned effect with respect to panels using patterns A is more significant than panels using patterns B. This is because of the advantages of FG-TPMS panels with patterns B as explained in Section 4.2.1. These advantages results in the panels with patterns B having higher bending and membrane stiffness than the panels utilizing patterns A. Hence, the effect of porosity volume fraction on the post-buckling strength of the FG-TPMS panel with patterns B is smaller than that with patterns A. Fig. 17 shows the load-deflection curves of panels with various porosity distributions and volume fractions. Clearly, effect of porosity distribution on post-buckling strength of FG-TPMS panel is significant. In the case of uneven porosity distribution, pores are concentrated near the neutral axis, where its contribution to bending stiffness is minimal. As a result for the same porosity volume fraction α, the panel with uneven porosity distribution has a higher post-buckling strength than the panel with even porosity distribution. It is recommended to use FG-TPMS panels with uneven porosity distribution in analysis and design.

Figure 13: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, patterns A1, A2, A3, W0=0, various porosity volume fractions, and even porosity distribution.

Figure 14: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, patterns B1, B2, B3, W0=0, various porosity volume fractions, and even porosity distribution.

Figure 15: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, patterns A1, A2, A3, W0=0, various porosity volume fractions, and uneven porosity distribution.

Figure 16: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, patterns B1, B2, B3, W0=0, various porosity volume fractions, and uneven porosity distribution.

Figure 17: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, pattern B3, W0=0, various porosity distributions and volume fractions.

4.2.4 FG-TPMS Panels with Geometric Imperfection and Porosity

In the last section, FG-TPMS panels with both geometric imperfection and porosity under axial compression are examined. Without loss of generality, we consider FG-TPMS panels with pattern B3. Fig. 18 shows the load-deflection curves of panels with various geometric imperfections and porosity volume fractions. As seen, the greater the geometric and material imperfections, the lower the post-buckling strength of FG-TPMS panel. If we simultaneously increase the geometric imperfection and porosity volume fraction, the post-buckling strength of the panel significantly decreases. Again, the simultaneous effect of porosity volume fraction and geometric imperfection on the post-buckling strength of the FG-TPMS panel with even porosity distribution is more significant than that of the panel with uneven porosity distribution. Finally, Fig. 19 shows the load-deflection curves of panels with various porosity distributions. Again, it is seen that the effect of porosity distribution on the post-buckling strength of FG-TPMS panel is significant. In addition, Fig. 20 shows the panel deformation during the buckling process.

Figure 18: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, pattern B3, various geometric imperfections and porosity volume fractions.

Figure 19: Load-deflection paths of CCSS FG-TPMS panels subjected to compression, pattern B3, various porosity distributions, geometric imperfections, and porosity volume fractions.

Figure 20: Structural deformation during the buckling process.

5  Conclusions

An investigation into the post-buckling behavior of FG-TPMS shells under axial compression was conducted for the first time in this paper. The proposed formulation considers geometric imperfections and porosity inside the panels. Two types of porosity distributions were considered: even and uneven. The post-buckling behavior of FG-TPMS shells was investigated considering six density distribution patterns along the thickness. The mechanical properties of the FG-TPMS materials were determined via a fitting technique. The present formulation was established by using IGA and FSDT. The NURBS basis functions were used to simulate exact geometries and approximate displacements. The non-linearity of shells was determined using the von Karman assumption and the total Lagrangian approach. Then, a modified Riks method was used to iteratively solve the discrete nonlinear equation system. The present results well agreed with those available in the literature. From the parametric study, the notable conclusions are drawn as follows:

1) For FG-TPMS perfect panels without geometric imperfection and porosity under axial compression:

•   The panels using patterns B (including B1, B2, and B3) have higher post-buckling strengths than those utilizing patterns A (including A1, A2, and A3). We can conclude that the panels with patterns B have higher bending and membrane stiffness than those utilizing patterns A. The panel with pattern B3 produces the highest post-buckling strength and the highest bending and membrane stiffness compared with the panels utilizing the remaining patterns A1, A2, A3, B1, and B2.

•   The higher the curvature (κ=L/R) of FG-TPMS panel, the higher the post-buckling strength.

2) For FG-TPMS panels with geometric imperfection but without porosity under axial compression:

•   All panels have similar critical buckling loads for an arbitrary density distribution pattern. After the load-deflection curve passes the second limit point the larger the geometric imperfection, the lower the post-buckling strength.

3) For FG-TPMS panels without geometric imperfection but with porosity under axial compression:

•   The porosity volume fraction significantly affects the post-buckling strength of FG-TPMS panel. The higher the porosity volume fraction, the lower the post-buckling strength of the FG-TPMS panel. We conclude that porosity significantly reduces the nonlinear strength of the FG-TPMS panels.

•   The effect of the porosity volume fraction on the post-buckling strength of the FG-TPMS panel with even porosity distribution is more significant than that of the panel with uneven porosity distribution. In the case of uneven porosity distribution, the mentioned effect with respect to panels using patterns A is more significant than panels using patterns B.

4) For FG-TPMS panels with both geometric imperfection and porosity under axial compression:

•   The greater the geometric and material imperfections, the lower the post-buckling strength of the FG-TPMS panel. If we simultaneously increase the geometric imperfection and porosity volume fraction, the post-buckling strength of the panel significantly decreases.

•   The simultaneous effect of porosity volume fraction and geometric imperfection on the post-buckling strength of the FG-TPMS panel with even porosity distribution is more significant than that of the panel with uneven porosity distribution.

5) The effect of porosity distribution on the post-buckling strength of the FG-TPMS panel is significant. For the same porosity volume fraction α, the panel with uneven porosity distribution has a higher post-buckling strength than the panel with even porosity distribution. It is recommended to use FG-TPMS panels with uneven porosity distribution in analysis and design.

Although the proposed approach and structure possess some advantages, there are several limitations which should be overcome as follows: 1) The proposed approach can be applied to singly or doubly curved shells and plates. It should be improved for analyses of shells with arbitrary geometries; 2) There is no experimental validation.

Acknowledgement: None.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: Study conception and design: Tan N. Nguyen, Aman Garg; Literature review and data collection: Tan Khoa Nguyen, Aman Garg; Analysis and interpretation of literature: Tan N. Nguyen, Canh V. Le, Mohamed-Ouejdi Belarbi; Visualization and graphical representation: Tan Khoa Nguyen, Canh V. Le; Draft manuscript preparation: Tan N. Nguyen, Mohamed-Ouejdi Belarbi; Critical revision of the manuscript: Tan N. Nguyen, Aman Garg, Tan Khoa Nguyen. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data will be made available on request.

Ethics Approval: No applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Sun Y, Li QM. Dynamic compressive behaviour of cellular materials: a review of phenomenon, mechanism and modelling. Int J Impact Eng. 2018;112:74–115. [Google Scholar]

2. Zheng Z, Yu J, Li J. Dynamic crushing of 2D cellular structures: a finite element study. Fifth Intl Symp Impact Eng. 2005;32(1):650–64. doi:10.1016/j.ijimpeng.2005.05.007. [Google Scholar] [CrossRef]

3. Al-Ketan O, Abu Al-Rub RK. Multifunctional mechanical metamaterials based on triply periodic minimal surface lattices. Adv Eng Mater. 2019;21(10):1900524. doi:10.1002/adem.201900524. [Google Scholar] [CrossRef]

4. Krishnan K, Lee DW, Al Teneji M, Abu Al-Rub RK. Effective stiffness, strength, buckling and anisotropy of foams based on nine unique triple periodic minimal surfaces. Int J Solids Struct. 2022;238(5):111418. doi:10.1016/j.ijsolstr.2021.111418. [Google Scholar] [CrossRef]

5. Abueidda DW, Abu Al-Rub RK, Dalaq AS, Lee DW, Khan KA, Jasiuk I. Effective conductivities and elastic moduli of novel foams with triply periodic minimal surfaces. Mech Mater. 2016;95(1838):102–15. doi:10.1016/j.mechmat.2016.01.004. [Google Scholar] [CrossRef]

6. Gado MG, Al-Ketan O, Aziz M, Al-Rub RA, Ookawara S. Triply periodic minimal surface structures: design, fabrication, 3D printing techniques, state-of-the-art studies, and prospective thermal applications for efficient energy utilization. Energy Technol. 2024;12(5):2301287. [Google Scholar]

7. Feng J, Fu J, Yao X, He Y. Triply periodic minimal surface (TPMS) porous structures: from multi-scale design, precise additive manufacturing to multidisciplinary applications. Intl J Extreme Manuf. 2022;4(2):022001. doi:10.1088/2631-7990/ac5be6. [Google Scholar] [CrossRef]

8. Nguyen-Xuan H, Tran KQ, Thai CH, Lee J. Modelling of functionally graded triply periodic minimal surface FG-TPMS plates. Compos Struct. 2023;315(10):116981. doi:10.1016/j.compstruct.2023.116981. [Google Scholar] [CrossRef]

9. Thai CH, Hung P, Nguyen-Xuan H, Phung-Van P. A meshfree method for functionally graded triply periodic minimal surface plates. Compos Struct. 2024;332(1):117913. doi:10.1016/j.compstruct.2024.117913. [Google Scholar] [CrossRef]

10. Tran KQ, Hoang TD, Lee J, Nguyen-Xuan H. Three novel computational modeling frameworks of 3D-printed graphene platelets reinforced functionally graded triply periodic minimal surface (GPLR-FG-TPMS) plates. Appl Math Model. 2024;126:667–97. doi:10.1016/j.apm.2023.10.043. [Google Scholar] [CrossRef]

11. Nguyen TN, Wattanasakulpong N, Nguyen NP, Fakharian P, Eiadtrong S. Isogeometric analysis of functionally graded triply periodic minimal surface shells. Mech Adv Mater Struct. 2025;32(21):5447–60. doi:10.1080/15376494.2024.2423278. [Google Scholar] [CrossRef]

12. Wang S, Song M, Yang J, Zhu W, Kitipornchai S. Free vibration of functionally graded TPMS porous quadrilateral plates reinforced with graphene nanoplatelets. Thin-Walled Struct. 2025;215(1):113440. doi:10.1016/j.tws.2025.113440. [Google Scholar] [CrossRef]

13. Abdoli E, Zarastvand MR, Talebitooti R. Wave propagation numerical simulation approach based on a novel TPMS-based lattice metamaterial for improved vibration transmission of doubly curved sandwich systems. Eng Comput. 2025;41(5):3737–54. doi:10.1007/s00366-025-02181-5. [Google Scholar] [CrossRef]

14. Zarastvand MR, Abdoli E, Talebitooti R. A lattice metamaterial-based sandwich cylindrical system for numerical simulation approach of vibroacoustic transmission considering triply periodic minimal surface. Appl Math Mech. 2025;46(11):2035–54. doi:10.1007/s10483-025-3314-9. [Google Scholar] [CrossRef]

15. Thai HT, Kim SE. A review of theories for the modeling and analysis of functionally graded plates and shells. Compos Struct. 2015;128(1):70–86. doi:10.1016/j.compstruct.2015.03.010. [Google Scholar] [CrossRef]

16. Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. Boca Raton, FL, USA: CRC Press; 2003. [Google Scholar]

17. Duc ND, Khoa ND, Thiem HT. Nonlinear thermo-mechanical response of eccentrically stiffened Sigmoid FGM circular cylindrical shells subjected to compressive and uniform radial loads using the Reddy’s third-order shear deformation shell theory. Mech Adv Mater Struct. 2018;25(13):1156–67. doi:10.1080/15376494.2017.1341581. [Google Scholar] [CrossRef]

18. Mantari JL. Refined and generalized hybrid type quasi-3D shear deformation theory for the bending analysis of functionally graded shells. Comp B Eng. 2015;83(7):142–52. doi:10.1016/j.compositesb.2015.08.048. [Google Scholar] [CrossRef]

19. Mantari JL, Guedes Soares C. Analysis of isotropic and multilayered plates and shells by using a generalized higher-order shear deformation theory. Compos Struct. 2012;94(8):2640–56. doi:10.1016/j.compstruct.2012.03.018. [Google Scholar] [CrossRef]

20. Viola E, Tornabene F, Fantuzzi N. General higher-order shear deformation theories for the free vibration analysis of completely doubly-curved laminated shells and panels. Compos Struct. 2013;95:639–66. doi:10.1016/j.compstruct.2012.08.005. [Google Scholar] [CrossRef]

21. Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng. 2005;194(39):4135–95. [Google Scholar]

22. Nguyen VP, Anitescu C, Bordas SPA, Rabczuk T. Isogeometric analysis: an overview and computer implementation aspects. Math Comput Simul. 2015;117:89–116. [Google Scholar]

23. Nguyen TN, Hien TD, Nguyen-Thoi T, Lee J. A unified adaptive approach for membrane structures: form finding and large deflection isogeometric analysis. Comput Methods Appl Mech Eng. 2020;369:113239. [Google Scholar]

24. Kiendl J, Bletzinger KU, Linhard J, Wüchner R. Isogeometric shell analysis with Kirchhoff-Love elements. Comput Methods Appl Mech Eng. 2009;198(49):3902–14. doi:10.1016/j.cma.2009.08.013. [Google Scholar] [CrossRef]

25. Guo Y, Ruess M. Nitsche’s method for a coupling of isogeometric thin shells and blended shell structures. Comput Methods Appl Mech Eng. 2015;284(5–8):881–905. doi:10.1016/j.cma.2014.11.014. [Google Scholar] [CrossRef]

26. Kiendl J, Bazilevs Y, Hsu MC, Wüchner R, Bletzinger KU. The bending strip method for isogeometric analysis of Kirchhoff-Love shell structures comprised of multiple patches. Comput Methods Appl Mech Eng. 2010;199(37):2403–16. doi:10.1016/j.cma.2010.03.029. [Google Scholar] [CrossRef]

27. Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. Isogeometric shell analysis: the Reissner-Mindlin shell. Comput Methods Appl Mech Eng. 2010;199(5):276–89. doi:10.1016/j.cma.2009.05.011. [Google Scholar] [CrossRef]

28. Casanova CF, Gallego A. NURBS-based analysis of higher-order composite shells. Compos Struct. 2013;104(10):125–33. doi:10.1016/j.compstruct.2013.04.024. [Google Scholar] [CrossRef]

29. Kiendl J, Schmidt R, Wüchner R, Bletzinger KU. Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting. Comput Methods Appl Mech Eng. 2014;274(2):148–67. doi:10.1016/j.cma.2014.02.001. [Google Scholar] [CrossRef]

30. Hosseini S, Remmers JJC, Verhoosel CV, De Borst R. An isogeometric continuum shell element for non-linear analysis. Comput Methods Appl Mech Eng. 2014;271(3):1–22. doi:10.1016/j.cma.2013.11.023. [Google Scholar] [CrossRef]

31. Leonetti L, Liguori F, Magisano D, Garcea G. An efficient isogeometric solidshell formulation for geometrically nonlinear analysis of elastic shells. Comput Methods Appl Mech Eng. 2018;331(7):159–83. doi:10.1016/j.cma.2017.11.025. [Google Scholar] [CrossRef]

32. Nguyen TN, Nguyen-Xuan H, Lee J. A novel data-driven nonlinear solver for solid mechanics using time series forecasting. Finite Elem Anal Des. 2020;171(3):103377. doi:10.1016/j.finel.2019.103377. [Google Scholar] [CrossRef]

33. Nguyen TN, Dang LM, Lee J, Nguyen PV. Load-carrying capacity of ultra-thin shells with and without CNTs reinforcement. Mathematics. 2022;10(9):1481. doi:10.3390/math10091481. [Google Scholar] [CrossRef]

34. Nguyen NV, Nguyen HX, Lee S, Nguyen-Xuan H. Geometrically nonlinear polygonal finite element analysis of functionally graded porous plates. Adv Eng Softw. 2018;126:110–26. doi:10.1016/j.advengsoft.2018.11.005. [Google Scholar] [CrossRef]

35. Nguyen TN, Thai CH, Nguyen-Xuan H, Lee J. NURBS-based analyses of functionally graded carbon nanotube-reinforced composite shells. Compos Struct. 2018;203(5):349–60. doi:10.1016/j.compstruct.2018.06.017. [Google Scholar] [CrossRef]

36. Nguyen TN, Thai CH, Luu AT, Nguyen-Xuan H, Lee J. NURBS-based postbuckling analysis of functionally graded carbon nanotube-reinforced composite shells. Comput Methods Appl Mech Eng. 2019;347(6348):983–1003. doi:10.1016/j.cma.2019.01.011. [Google Scholar] [CrossRef]

37. Piegl L, Tiller W. The NURBS book. Berlin, Germany: Springer; 1997. [Google Scholar]

38. Crisfield MA. A fast incremental/iterative solution procedure that handles “snap-through”. Comput Struct. 1981;13(1):55–62. doi:10.1016/0045-7949(81)90108-5. [Google Scholar] [CrossRef]

39. Nguyen TN, Lee S, Nguyen PC, Nguyen-Xuan H, Lee J. Geometrically nonlinear postbuckling behavior of imperfect FG-CNTRC shells under axial compression using isogeometric analysis. Eur J Mech-A/Solids. 2020;84(12):104066. doi:10.1016/j.euromechsol.2020.104066. [Google Scholar] [CrossRef]

40. Do VNV, Lee CH. Numerical investigation on post-buckling behavior of FGM sandwich plates subjected to in-plane mechanical compression. Ocean Eng. 2018;170(3):20–42. doi:10.1016/j.oceaneng.2018.10.007. [Google Scholar] [CrossRef]

41. Liew KM, Lei ZX, Yu JL, Zhang LW. Postbuckling of carbon nanotubereinforced functionally graded cylindrical panels under axial compression using a meshless approach. Comput Methods Appl Mech Eng. 2014;268:1–17. doi:10.1016/j.cma.2013.09.001. [Google Scholar] [CrossRef]

42. Duc ND, Thang PT. Nonlinear response of imperfect eccentrically stiffened ceramic–metal–ceramic FGM thin circular cylindrical shells surrounded on elastic foundations and subjected to axial compression. Compos Struct. 2014;110(3):200–6. doi:10.1016/j.compstruct.2013.11.015. [Google Scholar] [CrossRef]