12892 FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 23, No 4, 2025, pp. 827 - 860 https://doi.org/10.22190/FUME240920015M © 2025 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper POROSITY EFFECTS ON NONLINEAR STATIC PERFORMANCES OF FUNCTIONALLY GRADED SHELLS CONSIDERING THICKNESS STRETCHING Hanen Mallek1,2, Hana Mellouli2,3, Lotfi Ben Said4, Mondher Wali2,5, Fakhreddine Dammak2, Muapper Alhadri4 1Higher Institute of Applied Sciences and Technology of Sousse, ISSAT Sousse, University of Sousse, Tunisia 2Laboratory of Electrochemistry and Environment (LEE), National Engineering School of Sfax, (ENIS), University of Sfax, Tunisia 3Gafsa Preparatory Engineering Institute, IPEI Gafsa, University of Gafsa, Tunisia 4Department of Mechanical Engineering, College of Engineering, University of Ha’il, Ha’il City, Saudi Arabia 5Higher School of Sciences and Technology of Hammam Sousse, ESSTHS, University of Sousse, Tunisia ORCID iDs: Hanen Mallek https://orcid.org/0000-0003-1681-9558 Hana Mellouli https://orcid.org/0000-0001-6739-1550 Lotfi Ben Said https://orcid.org/0000-0001-7652-3068 Mondher Wali https://orcid.org/0000-0001-8384-8936 Fakhreddine Dammak https://orcid.org/0000-0002-9629-671X Muapper Alhadri https://orcid.org/0000-0001-5848-5803 Abstract. This study delves into the intriguing realm of nonlinear responses exhibited by porous functionally graded (FG) 3D shell structure. A power law approach is formulated to simulate the mechanical behavior of FG materials, considering two distinct porosity distributions. This approach provides a comprehensive exploration of porosity phenomena in FG materials. The finite element (FE) formulation is based on an improved first-order shear deformation theory (FSDT) with the inclusion of thickness stretching parameters. The theoretical framework improved in this manner provides a more accurate representation of the transverse shear stress distribution within the structure, capturing the complexities of its behavior under loading conditions. This research enriches understanding of FG materials by integrating porosity that may be distributed evenly or unevenly, thus contributing to the advancements in the field. Key words: Porous FGM, Thickness stretching, FE analyses, Porosity, Shell structure Received: September 20, 2024 / Accepted April 05, 2025 Corresponding author: Hanen Mallek Laboratory of Electrochemistry and Environment (LEE), National Engineering School of Sfax, (ENIS), University of Sfax, Route de Soukra km 4, 3038 Sfax, Tunisia E-mail: hanen.mallek@enis.tn 828 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI 1. INTRODUCTION Functionally Graded Materials (FGMs) have garnered increasing interest in modern technologies over the last decades due to their unique composition, synthesized from a blend of two distinct material phases in specific volume ratios. These materials exhibit gradual variations in properties along specified directions, typically through-thickness, thereby mitigating abrupt changes in thermal or mechanical properties and reducing interlaminar stresses. Leveraging the remarkable advantages of FGMs in industrial and engineering applications, numerous researchers have focused on analyzing the behaviors of FGM plates and shell structures using finite element (FE) methods [1]. Recent studies have also explored the integration of FGM materials with piezoelectric materials to create smart structural components. For instance, Mallek et al. [2] investigated the piezoelastic response of FGM shells embedded with intelligent materials under mechanical and electrical loads, highlighting the potential for advanced multifunctional applications in structural engineering. The materials industry has witnessed a significant evolution, with porous materials gaining prevalence owing to their distinct characteristics, including lightweight nature, excellent energy absorption, and heat resistance [3]. These attributes have motivated the emergence of a novel category of lightweight materials known as Functionally Graded Porous (FGP) materials, characterized by a graded distribution of internal pores within the microstructural section [4]. FGP materials, offering unique mechanical properties, have recently found applications in various domains, including static, free vibration, and dynamic analysis. Numerous studies have investigated both linear and nonlinear behaviors of FGP beams [5-7]. For instance, Wang and Wu [8] conducted free vibration analysis of cylindrical porous FG shells, combining sinusoidal shear deformation theory with the Rayleigh-Ritz method to discern the effect of graded porosity on structural characteristics. Additionally, Nguyen et al. [9] proposed an efficient approach based on the polygonal finite element (FE) method to predict static and free vibration responses of FGP plates reinforced with graphene platelets (GPLs), considering various porosity distributions. Subsequent works by other authors [10-12] have explored the bending, buckling, and vibration behaviors of FGP structures, highlighting the influence of porosity distributions and mechanical properties. Notably, Zghal and Dammak [13, 14] investigated buckling responses and vibrational behavior of FGP plates and shells with varying mechanical properties using a mixed FE formulation, emphasizing the critical role of porosity percentages and distributions. The literature review underscores the importance of considering porosity effects and distributions in predicting the mechanical behavior of FGP materials. Indeed, Wattanasakulpong et al. [5] demonstrated the emergence of porosity and micro-voids during the multi-step infiltration technique in layered FGMs, attributed to solidification temperature differences among material phases. Two models for studying the mechanical behaviors of FGP materials have been proposed: one based on mixture rules, assuming continuous variation of effective material properties and even or uneven dispersion of porosities [15], and another employing a specific continuous function to determine porous material properties throughout the structure's thickness [16]. While considerable attention has been devoted to analyzing the mechanical behavior of FGP beams and plates, there is a notable dearth of investigations concerning porous FG shell structures, as noted by Wang and Wu [8]. Thus, it is imperative to explore the behaviors of porous FG shells to unveil the effects of porosity. Shell structures find widespread applications [17] in various industrial fields, including pressure vessels, roll- forming, and deep drawing, often undergoing large deformations [18] and finite rotations Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 829 where linear models prove inefficient. Consequently, accounting for geometric nonlinearities becomes crucial in shell modeling, given the significant influence of nonlinear parameters in numerous analysis fields, particularly in porous FG structural design. Within this framework, Abid et al. [19] examined the geometrically nonlinear static responses of porous FG ring and cylindrical panels, considering even and uneven porosity distributions along the thickness direction. The High-order Shear Deformation Theory (HSDT) employing double directors' vectors was utilized with an isoparametric concept, highlighting the disparities between linear and nonlinear deflections of FGP shell structures. It was deduced that under high loading conditions, linear shell modeling may yield inaccurate predictions, emphasizing the importance of geometric nonlinearity in FE formulations. Moreover, when modeling moderate or thick shell structures under the Plane Stress Assumption (PSA), neglecting the thickness stretching effect of normal deformation may lead to inaccuracies, especially for thick structures. In this regard, evolving an efficient FE formulation capable of accommodating arbitrary three-dimensional constitutive laws without condensation is crucial. To achieve this objective, the standard FSDT shell formulation with five parameters is extended to a 6- parameter theory, including the full set of stress and strain fields. This extension involves considering an extensible director vector, acknowledging the sixth degree of freedom as a scalar stretching component. The exploration of mechanical behaviors using the 3D theory of elasticity was pioneered by Simo and Rifai [20] employing a 6-parameter shell model, which avoids displacement formulation errors and thickness locking phenomena. Various methods have been proposed to address thickness locking, including quadratic displacement distribution through the thickness direction and enriching the thickness strain via methods like the Enhanced Assumed Strain (EAS) method. Recent efforts by Beheshti and Ramezani [21] introduced a nonlinear extensible shell model to investigate the large deformation of FG shell structures, incorporating seven parameters to account for thickness stretch, while employing the EAS concept to mitigate numerical instability issues and enhance strain components. Triply Periodic Minimal Surface (TPMS) structures are currently renowned for their exceptional energy absorption capabilities, which make them highly suitable for applications requiring efficient mechanical performance. In this regard, several recent studies have explored advanced approaches such as nonlocal strain gradient-based isogeometric analysis for functionally graded TPMS nanoplates reinforced with graphene platelets. These studies aim to improve the accuracy of mechanical behavior predictions for such complex materials. Furthermore, new computational frameworks [22, 23] have been proposed for modeling 3D- printed graphene-reinforced functionally graded TPMS plates, emphasizing multiscale techniques to better capture the mechanical properties and performance in additive manufacturing contexts. These advancements offer a deeper understanding of the behavior of graphene-based composites and present innovative computational strategies for their analysis. From the literature review, it is evident that previous studies primarily focused on analyzing the responses of FG plates and shell structures with internal pores, employing the Plane Stress Assumption (PSA) model, which neglects thickness stretch. However, a comprehensive investigation into the large elastic deformations and finite rotations of FG shell structures with internal pores, considering three-dimensional constitutive equations and through-the-thickness stretching, remains unexplored. This paper aims to address this research gap. Our work introduces a 3D extensible shell model utilizing six independent kinematic parameters, including thickness stretching, and incorporates the EAS method. Considering the limitations of conventional theories like FSDT and HSDT, an efficient finite element shell model is proposed based on an improved FSDT theory, ensuring a realistic parabolic distribution of Transverse Shear Strains (TSS) across the structure's 830 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI thickness. Additionally, the importance of thickness stretching is highlighted by incorporating linear variation through the thickness via the EAS concept to analyze large elastic deformations. Moreover, locking problems inherent in extensible models are addressed by applying the Assumed Natural Strain (ANS) method to improve membrane behavior. The subsequent sections of this paper delve into the materials modeling and geometrically nonlinear analysis methods considered herein. The proposed finite element method necessitates an effective elastic part, including offsetting the neutral surface, accounting for geometric nonlinearity in stiffness matrices, enhancing quantities for bending and membrane deformations using the EAS method, and improving TSS using a modified FSDT theory and ANS method. 2. GENERALIZED FORMULATION AND FE MODELING OF FGP SHELLS The study of shell structures is significant in several fields of structural mechanics. With the expansion of FE methods, various theories have been proposed to examine shell issues. In this work, the extensible director shell approach is used to develop the governing mathematical expressions for the elastic case to evaluate the large deflection performance of porous functionally graded structures. Three essential concepts underpin the extensible FSDT shell theory. The first concept illustrates the position of the director dispersion in the deformed state. The second concept, called parabolic TSS, assumes that the distribution of TSS across the structure's thickness is formulated as polynomial of degree two. The inclusion of the thickness stretching parameter in shell elements is established in accordance with the third concept. 2.1. FGP Properties’ Estimation The increasing application of functionally graded materials in thin-walled structures necessitates the development of accurate predictive models to characterize their material behavior. In this investigation, the rule of mixture (ROM) is employed to analyze the performance of these advanced materials. FGMs can be categorized into two classes: perfect and imperfect. Perfect FGMs are composed of two phases: ceramic and metal where the composition and properties change gradually within the material. Imperfect FGMs, on the other hand, are characterized by the presence of internal pores within their microstructure (see Fig.1). Fig. 1 Perfect and imperfect FGM Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 831 2.1.1. Perfect FGM Material Properties FGMs are composite materials composed of distinct phases, typically metals and ceramics, which vary gradually in volumetric fraction along the thickness direction. This gradual variation results in material properties that smoothly transition from one constituent to another. Typically, the top surface consists of ceramic while the bottom surface is made of metal. In this study, the mechanical properties of the material evolve according to a power law distribution (where n denotes the power low index), considering the specific material properties of the ceramic and metal phases (as illustrated in Fig. 2 with the porosity volume fraction α set to 0). It is worth mentioning that when n equals zero, the material is entirely ceramic, and as n tends to infinity, the material transitions to being entirely metal. In situations where there are asymmetrical variations in material properties, particularly in terms of E, the offset 𝑧̅ of the neutral surface from the discretized mid-plane is defined as detailed in reference [24]. /2 /21 /2 1 /2 ( ) ; ( ) h h h h E E z dzE z E E E z z dz − −  =  =   =    (1) where �̅� and �̅�1 are the effective moduli for membrane and bending, respectively. Subsequently, the analytical representation of 𝑧̅ /h is approximated as [24]: ( 1)( 1) / ; 2( 1)( 2)( ) c m En k n z h k n n n k E − + = = + + + (2) In this study, the modified rule of mixture (MROM) refers to an enhancement of the traditional rule of mixture used in composite materials analysis. The rule of mixture (Fig. 2), also entitled the FGM homogenization schema, is used to predict the overall properties of an FGM based on the properties of its individual constituents and their volume fractions. In the context of Wattanasakulpong and Ungbhakorn's work [25], the modified rule of mixture extends this concept by incorporating additional factors or adjustments to more accurately predict the behavior of functionally graded materials (FGMs). These factors may firstly refer to the porosity volume fraction α that defines the proportion of the FGM's volume that is occupied by voids or pores. These porosities can be distributed either evenly or unevenly (as shown in Fig. 2) across the thickness of the structure. It is assumed that the top surface is predominantly ceramic, while the bottom surface is primarily metal, with a continuous variation of the material composition throughout the thickness. The material properties of the porous FGM follow the rule of mixtures, where the porosities are distributed uniformly across the thickness, as described by: ( ) ( ) ( ) 1 2 2 c m c m c m z P z P P V P P P b h    = − + − + −    (3) where α (0≤ α ≤1) represents the porosity volume fraction, 𝑃𝑚 and 𝑃𝑐 denote the material properties of the metal and ceramic respectively such as the Young’s modulus E, the Poisson’s ratio υ and the density ρ. Further, 𝑉𝑚 and 𝑉𝑐 indicate the volume fractions of metal and ceramic respectively assuming the condition of 𝑉𝑚 + 𝑉𝑐 = 1 and h denotes the structure’s thickness. 832 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI Fig. 2 FGM homogenization scheme When α is zero, it indicates that the FGM is completely dense without any voids or pores. Increasing this factor results in more voids along the FGM’s volume. Secondly, the factor b represents the uneven parameter. Setting b to zero means that the porosity is evenly distributed throughout the material. Conversely, setting b to unity indicates an uneven distribution of porosity within the FGM (Fig. 2). Not forgetting the power law index n of the FGM that determines the rate at which the material composition and properties change along the material's gradient. In fact, the volume fraction of the ceramic phase varies along the structure's thickness, denoted as:   1 ( ) ; / 2, / 2 2 n c z V z z h h h   = +  −    (4) As the power law index n increases, the FGM structure becomes richer in metal phase, while lower values of n result in a structure that is predominantly ceramic. When n is set to zero, the structure is entirely ceramic, and as n increases, the structure shifts towards a higher metal content. 2.1.2. Imperfect FGM Material Properties Porosities arise as a consequence of the manufacturing process, manifesting as defects within the microstructure. It is imperative to account for these porosities when studying the mechanical behavior of FGMs. In this study, two distinct porosity distributions along the thickness direction are examined: even and uneven, as depicted in Fig. 1. In contrast to the even distribution of porosity, where the distribution is uniform across the thickness, the uneven distribution reveals a concentration of porosities primarily within the middle zone. Fig. 3 illustrates the evaluation of effective FGP material properties along the thickness, considering both even and uneven distributions of porosities. The exponent n represents the material power-law index describing the material profile gradation. When n is set to zero, a single-phase ceramic material is obtained. The porosity volume fraction α varies from zero to 1, Setting α to zero implies a perfect FGM shell, while increasing introduces porosity and micro-voids, rendering the FGM imperfect. Further, b represents the uneven parameter. When it is set to zero, an evenly distributed porosity is achieved, whereas setting it to unity indicates an uneven porosity distribution. Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 833 It is important to acknowledge that the variations of porous FG material properties along the thickness direction are asymmetric. In mixed finite element formulations, this asymmetry can pose challenges during simulation due to the coupled nature of bending and membrane effects. To address this issue, the concept of a neutral surface can be employed. All effective quantities are then evaluated relative to this neutral surface. The analytical expression for the offset of the neutral surface from the discretized mid-plane is given as follows [26]: ( 1)( 1) / ; ( 1)( 2)(2( ) ( 1)( 1)(1 0.5 )) c m En k n z h k n n n k k n b E − + = = + + + − + + − (5) The distribution of the offset 𝑧̅/ℎ along with the variation of the material exponent n and the porosity volume fraction α is plotted in Fig. 3. a) b) Fig. 3 Variation of the offset 𝑧̅/ℎ versus different values of the material exponent n and porosity volume fraction α for the (a) Even; (b) Uneven configuration 834 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI 2.2. Extensible Shell Kinematics and Fundamental Equations This section provides a brief overview of the kinematics underlying the nonlinear extensible director shell model. The initial and deformed configurations are denoted as C0 and Ct, respectively. The present research utilizes curvilinear dimensions (𝜉1, 𝜉2, 𝜉3 = 𝑧) to characterize the location and placement of physical parameters. 2.2.1. Extensible Director Shell Kinematics In the following, the FSDT model, also sometimes denoted as 5-parameter model, is here enhanced with an extensible director, thus being able to depict the shell thickness stretch during its deformation. According to the cited theory, the displacement field of an arbitrary point (q) on the shell, in both its initial state and its current configuration can be expressed using curvilinear coordinates:   1 2 1 2 1 2 0 01 2 1 2 1 2 ˆ( , , ) ( , ) ( , ) ˆ; ; / 2, / 2 ˆ( , , ) ( , ) ( , ) q p q p z z z z z z h h z z              = + = −  − = + X X D x x d (6) in which the in-plane curvilinear coordinates (𝜉1, 𝜉2) are tangent to the mid-surface, whereas the thickness coordinate 𝜉3 is perpendicular to this mid-surface, as depicted in Fig. 4. Xp and xp represent the position vector located on the mid-surface in C0 and Ct, respectively. ho refers to the initial thickness of the structure and D designates the director shell vectors in C0. Fig. 4 Geometry and kinematics of the shell element The initial improvement to the shell's kinematic properties involves incorporating an extensible director vector using an extra variable known as the shell thickness stretch. Consequently, the deformation of the director within Ct is characterized by both an inextensible unit vector t and a thickness stretch parameter represented as λ: Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 835 1 2 1 2 1 2 1 2( , ) 1 ( , ) ( , ); ( , ) 1         = − = d t t (7) According to Eq. (4), the vectors xq and d are treated as independent quantities, with xq expressed by three components, while d is described by only two. The shell thickness can be determined in relation to the thickness stretching parameter λ (𝜉1, 𝜉2) as follows [27]: 0 01 / , (1 )h h h h − = = − (8) It is worth mentioning that for the reference configuration C0 where there is no thickness stretch, the scalar function λ is initially defined as zero. The derivative of the director vector d is expressed as below: , , ,(1 ) , 1,2k k k k = − − =d t t (9) For the ease of defining the shell deformation gradient, the covariant basis vectors are determined in the initial and actual configurations as: , 3 , 3 ; , 1,2 ˆ ; k k k k k k z k z = + = = = + = G A D G D g a d g d (10) where (𝑨1, 𝑨2, 𝑨3) denote the tangent basis vectors of the mid-surface in the undeformed reference configuration and (𝒂1, 𝒂2, 𝒂3) are their image at the actual configuration. The Riemannian metric at the initial configuration C0 is designated by G, that of the current configuration Ct by g, having the following expressions: . , , 1,2,3 . i j i j i j   =   =  =   G G G g g g (11) Introducing the expressions of the covariant basis vectors and the Riemannian metric (Eqs. (10) and (11)) yields the covariant components of the 3D metric tensor in Ct as: 2 0 1 ij i j ˆ ( ) ˆ. , + , , 1,2,3, , 1,2 . kβ kβ k k3 k k 33 g a zb O z g g c zc i j k g    = + +  = = = =  = g g d d (12) where 𝑎𝑘𝛽, 𝑏𝑘𝛽, 𝑐𝑘 0 and 𝑐𝑘 1 represent the kinematic variables, which are expressed as: , , 0 1 , . , . . , , 1,2 . , . kβ k β k k β β k k k k k a b k c c    = = + = = = a a a d a d a d d d (13) For later use, geometrical variables are also added in the initial state C0 as: 1 2 1 2d d , , d d dξ ξA A A A A ξ ξ= =  =A A (14) 836 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI 2.2.2. Green Compatible Strain Tensor A suitable strain measure, accommodating arbitrarily large rotations and displacements, is provided by the Green compatible strain tensor 𝑬𝑐, given as: ( ) 2 c 0 1 ij ij 3 c 2 33 33 ˆ ( ) 1 1 ˆ( ), ( ), 2 +z , , 1,2 2 2 (1 ) 1 / 2 c kβ kβ k c c ij k k k E e zχ O z E g G E k E e       = + +  = − = − = =  = = − − E g G (15) The components of strain field, including membrane 𝑒𝑘𝛽, transverse normal 𝑒33, bending 𝜒𝑘𝛽 and shear (𝛾𝛼 0, 𝛾𝛼 1 ) strains, are defined as energetically conjugate parts to the kinematic variables. 0 0 0 1 , ( ) / 2 , ( ) / 2 , , 1,2 , ( 1) kβ kβ kβ k k k k k k k k e a A χ b B k γ c C γ        = − = − = = − = − (16) By utilizing these relations, the variation of the corresponding strains can be evaluated as: , , 33 , , , , , , 0 1 , , , ( . . ) / 2 , ( 1) ( . . . . ) / 2, , 1,2 . . , ( 1) kβ k β β k k k β β k k β β k k k k k k k δe δ δ δe δχ δ δ k δγ δ δ δγ            = + = −  = + + + =  = + = + − a x a x a d a d x d x d a d x d (17) Whereby, the variation of the extensible shell director is computed as follows: , , , , ,(1 ) , (1 )k k k k k          = − − = − − − −d t t d t t t t (18) The cited quantities are condensed obtaining the following matrix form: 11 11 0 1 22 22 0 11 1 0 1 33 2 2 12 12 , , , 0 2 2 e e e e                       = = = =                    e χ γ γ (19) Next, these strain vectors are organized into a generalized strain 𝞢𝑐 as follows: 0 1 12 1 c       =       e χ Σ γ γ (20) 2.2.3. Improved FSDT Theory The Mindlin's First-Order Shear Deformation Theory (FSDT) is based on the assumption of a linear shear strain distribution 𝐸𝑘3 𝑐 across the thickness of the shell. However, it is well Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 837 established that shear strain follows a parabolic distribution across the thickness, vanishing at the top and bottom surfaces of the shell. To correct this discrepancy, several studies have introduced constant shear correction coefficients. Additionally, Tanov and Tabiei [28] incorporated a shear function in the FSDT formulation, resulting in a parabolic distribution of transverse shear stresses. Building on the work of Shi [29], the present study improves the FSDT theory by assuming a parabolic shear strain distribution through the thickness, using a quadratic function f(z) as detailed in the following: 2 0 1 3 2 5 ˆ2 ( )( ), ( ) 1 4 , 1,2 4 c k z z E f z z f z k h   = + = − =    γ γ (21) 2.3. Shell Variational Formulation The Enhanced Assumed Strain (EAS) method is employed for defining the strain tensor, aiming to prevent numerical locking problems and enhance the deformation gradient. According to this technique, the Green strain field E comprises the compatible strain field 𝑬𝐶 and the enhanced incompatible strain field �̃�: c= + E E E (22) The EAS formulation is employed to represent the three-dimensional equilibrium equations in terms of three independent variables: the Green strain tensor E, the 2nd Piola- Kirchhoff stress tensor S, and the deformation field Φ. As demonstrated in the work of Simo and Rifai [20], the following variational formulation is articulated: int int0 ; : ( )c ext V W W W W dV = − = = + S E E (23) where 𝑊𝑒𝑥𝑡 represents the virtual work of the external loading and dV expresses the shell initial volume element. Note that the stress tensor S is considered to be orthogonal to the enhanced strain tensor �̃� over the element as outlined by Simo et al. [30]. This assumption of orthogonality is expressed as: 0 V dV = E S (24) Moreover, no additional shear enhancement components are included in �̃�. Thus, the internal virtual work, upon integration across the thickness of the shell, can be expressed as: ( )0 1 int 0 1( ) ( ) z A W dA     = + + + + + N e e M χ χ T γ T γ (25) where the terms (N, M, 𝑻𝑘(k=0,1)) represent the membrane, bending, and shear stress resultants expressed in vector form as: 11 11 1 122 22 0 1 0 12 233 33 0 1 12 12 N M T TN M , , T TN M N M               = = = =                     N M T T, (26) 838 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI where its constituents are defined as following: h /2 h /2 33 33 -h /2 -h /2 h /2 h /2 33 33 -h /2 -h /2 h /2 h /2 0 1 -h /2 -h /2 / d , / d ˆˆ / d , z d , , 1,2 ˆ( ) / d , / d kβ kβ k kβ k k3 k3 N S G A z N S G A z M zS G A z M S z k T f z S G A z T zS G A z     = =   = = =   = =        (27) To simplify the expression of the internal work as shown in Eq. (23), the stress and strain resultants are organized in a generalized vector form, yielding: 0 0 1 1 12 1 , ,c                           = = =                     N e e M χ χ R Σ Σ T γ 0 T γ 0 (28) Thus, Eq. (25) can be expressed as follows: ( )c T int A W dA = + Σ Σ R (29) 2.4. Governing Equations of FGP Shell A linear relationship is assumed between the generalized stress resultant R and the generalized strain 𝞢 as follows: T .=R H Σ (30) where 𝑯𝑇 represents the material tangent modulus, described in terms of the linear elastic in-plane H and out-of-plane 𝑯𝜏 sub-matrices: ( ) 11 22 33 34 44 2 2 11 12 22 2 2 2 2 2 2 33 34 44 2 ( ) (1 ) ( ) ( ) ( ) ( ) T h/ h / h / h / Sym ˆ ˆ, , ,z,z dz ˆ ˆ, , f z ,zf z ,z f z dz −  −      =        =   =    H 0 0 0 H 0 0 H H H H H H H H H H H H (31) For isotropic materials, the matrices H and 𝑯𝜏 are specified as follows: 1 0 1 0 1 0( ) 1 0 0 1(1 )(1 2 ) 2(1 ) 0 0 0 1 2 ) 2 E z E , /  −       −     = =     −+  −  +      −   H H (32) Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 839 2.5. FE Modelling of FGP Shell In this section, the FE implementation of the current theoretical formulation based on a nonlinear four-node shell element is established according to the extensible director theory to predict the static FGP responses of plate and shell structures. Thus, the development of an isoparametric FE, which involves rotation-displacement-thickness-stretch, can be achieved. This extensible shell element according to the director model encompasses a total of six degrees of freedom at each node I: including 3 displacements, 2 rotational degrees and the stretching parameter presented in the generalized displacement vector as 𝜱𝑰 = (𝒖𝑰, 𝒕𝑰, 𝜆𝑰). The displacement vector, written in the form u=x-X and its variation are approximated using the standard finite element bilinear shape functions 𝑁𝐼, as: 4 1 ( ) ( )I I I I I , , N , , =  = X u u X u u (33) The unit director vector t and the λ parameter are discretized using the same functions: 4 1 ( ) ( )I I I I I , , N , , =   =  t t t t (34) As the approximation of the different fields is assessed in the Local Cartesian system, Eq. (35) defines the Jacobian matrix J to facilitate the transformation from curvilinear coordinate vectors to the Local Cartesian system {𝒏1 0, 𝒏2 0, 𝒏𝟎}. Here, 𝒏𝟎 represents the normal field to the mid-surface in the initial configuration, computed as 𝒏𝟎 = 𝑨𝟏 ∧ 𝑨𝟐/‖𝑨𝟏 ∧ 𝑨𝟐‖: 0 0 1 1 2 1 0 0 1 2 2 2 . . . .   =      n A n A J n A n A (35) Subsequently, the derivatives of the Cartesian standard shape functions defined at node I are converted to the global Cartesian basis as follows:   1,1 ,1 ,2 ,2 I I I I N N N N −       =           J (36) 2.5.1. Compatible Membrane and Transverse Normal Strain As detailed in the study by Bischoff and Ramm [31], to mitigate the occurrence of parasitic strain distribution, the Assumed Natural Strain (ANS) method is employed to interpolate the compatible thickness strain defined at the sampling points instead of the Gauss points (refer to Fig. 5). 4 33 33 33 1 , ( 1)I I I I I I e N e e     = = = − (37) 840 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI The interpolated compatible membrane strains (𝑒11, 𝑒22 and 𝑒33) and transverse normal strain 𝑒33 are connected to the node’s kinematic variables 𝛿𝜱𝑛 through the discrete operator matrix 𝑩𝑚, as outlined in Appendix A .m n =e B Φ (38) Fig. 5 Presentation of transverse normal and shear strain interpolation points 2.5.2. Compatible Bending Strain The linearized compatible bending strain is written in term of the node’s kinematic variables 𝛿𝜱𝑛 and the discrete bending operator matrix 𝑩𝑏 (detailed in Appendix A), as follows: .b n =χ B Φ (39) 2.5.3. Compatible First Term of the Transverse Shear Strain It is hypothesized that the shear locking issue, which hampers the effectiveness of the displacement Finite Element (FE) formulation, particularly for thin shells, can be effectively mitigated by the Assumed Natural Strain (ANS) method [32, 33]. The compatible first term of the TSS 𝜸0 is computed within the element at sampling points (A, B, C and D) as follows: ( ) ( ) ( ) ( ) 2 0 2 00 1 10 1 1 0 1 00 2 22 (1 ) (1 )1 (1 ) (1 )2 B D A C              − + + = =    − + +    γ (40) Here, (A, B, C and D) represent a set of mid-points along the boundaries of the element, as illustrated in Fig. 5. This interpolation technique safeguards against unreasonable modes of compatible shear strains, enabling deformations without parasitic shear strains during pure bending. Next, the variation of the constant TSS is interpolated within the element as: 0 0 .s n =γ B Φ (41) Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 841 in which, 𝑩𝑠 0 is the discrete first shear strain operator matrix, illustrated in Appendix A. 2.5.4. Compatible Second Term of the Transverse Shear Strain The calculation of the compatible second TSS follows the same procedure as the first shear strain, utilizing the ANS method. 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2 2 2 (1 ) ( ) (1 ) ( )1 2 (1 ) ( ) (1 ) ( ) B D A D                − + + = =    − + +    γ (42) The variation of the second TSS strain is presented in the matrix form as follows: 1 1.s n =γ B Φ (43) Here, the discrete operator matrix for the second shear strain-displacement 𝑩𝑠 1 is depicted in Appendix A. Ultimately, the virtual compatible strain tensor is obtained using Eqs. (39), (40), (42), and (43) as follows: 0 1 . T c n      = = Σ e χ γ γ B Φ (44) 2.5.5. Enhanced Transverse Normal Strains Approximation The enhanced strain field �̃� according to the EAS method is expressed as following: 0 0 E G G =E T Mα (45) where α refers to the internal variable vector. �̃� denotes the interpolation function matrix and the index ‘0’ represents its evaluation at the center of the element. Four configurations of enhancement for the �̃� matrix are developed with one, four, five, or eight parameters. In Eq. (18), a linear distribution is assumed along the shell thickness for the in-plane normal strains 𝐸𝛼𝛽 𝑐 and 𝐸𝛼3 𝑐 , while the transverse normal strain �̃�33 is presumed constant with respect to the thickness coordinate z. This constant distribution may result in artificial stresses in the thickness direction. To address this issue, the transverse normal strain 𝐸33 𝑐 should be augmented by an additional linear strain �̃�33. This enhancement component is introduced discretely using the Enhanced Assumed Strain (EAS) concept. According to the EAS approach, the interpolation of �̃�33 is treated independently for each element and is introduced via additional parameters, which are condensed within the element. The following enhancement is achieved using either one or four additional parameters. The expression of the interpolation matrix when considering one additional parameter, �̃�1 is introduced as follows:  1 2 0 0 0 0 0 T z , h  = =M (46) Next, the variation of the generalized incompatible strain 𝛅�̃� is determined using the interpolation matrix �̃�1 defined in the curvilinear basis (refer to Appendix B), as follows: 842 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI 1 =Σ B  (47) If four strain-supplemented parameters are taken into account, the expression of the interpolation matrix �̃�4 can be described as outlined in Appendix B. 2.5.6. Enhanced Membrane Strains Approximation To tackle the membrane locking, the EAS method is employed. Within this approach, the membrane strains are enhanced through incompatible strains. Hence, selecting the additional enhanced strain modes is crucial. Four bilinear enhanced strain components can be incorporated to form a complete polynomial, resulting in: 1 11 2 22 1 2 12 0 0 0 span 0 0 0 2 0 0 e e e                                                           (48) This technique was initially explored by Bischoff and Ramm [31] to improve the compatible membrane strains. The interpolation for the enhanced strains can now be articulated as follows: 1 1 1 2 0 2 4 4 4 4 2 0 1 2 3 3 0 0 0 , , 0 0 0 0 0 T A A     −            = = =               G G G e G α G G T 0 G G  (49) where 𝑇𝟎 is the following matrix evaluated at the center of the element as: 2 2 2 2 11 21 11 21 2 2 2 2 12 22 22 12 2 2 2 2 2 2 2 2 11 12 22 21 11 22 12 21 2 2 J J J J J J J J J J J J J J J J     =    +  T (50) Here, 𝐽𝛼𝛽 represent the Jacobian components expressed in Eq. (36). The fifth and eighth strain parameters are then derived, considering not only the enhanced thickness strain component but also the membrane strains, allowing for the matrices �̃�5 and �̃�8, as illustrated in Appendix B. It is now appropriate to summarize the relations between strain and kinematic variables at the node for the compatible and enhanced strains in matrix notation as follows: . , .n   = =Σ B Φ Σ B α (51) 2.5.7. Linearization Eq. (24) delineates the nonlinear shell model, which is solved using the Newton iterative algorithm. The consistent tangent operator for the Newton procedure can be derived from the directional derivative of the weak form in the direction of the increment ΔΦ: Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 843   0 T int ext T n nW DW         −  + = + =             ΦK L f f Φ α αL H h (52) where K, L and H are given by: , , T D G D T A T T T T A A dA dA dA = + = = =    K K K K B H B L B H .B H B H B (53) 𝑲𝐷 and 𝑲𝐺 are the material and geometric tangent matrices. The components of the geometric tangent matrix 𝑲𝐺 are given in the Appendix C. In Eq. (50), 𝒇𝑖𝑛𝑡 and h, are given by the following expressions ,int T T A A dA dA= =  f B R h B R (54) Removing the internal variable vector α from Eq. (53) at the element level yields the element tangent matrix 𝑲𝑇 , expressed as: 1T T D G −= + −K K K L H L (55) and the residual vector R defined as: 1T ext int−= −R L H h+ f f (56) At the element level, transitioning from the generalized displacement vector 𝜱𝑛 = (𝒖, 𝒕, 𝜆)𝑛 to the nodal variable vector 𝜞𝑛 = (𝒖, 𝜣, 𝜆)𝑛, where 𝜣 represents the two rotations vector in material description, is significant. At node I, the relationship between the generalized displacement vector 𝜱𝑛 = (𝒖, 𝒕, 𝜆)𝑛 and the nodal variable vector 𝜞𝑛 = (𝒖, 𝜣, 𝜆)𝑛 is elucidated by: 7 6 1 4I I I I I 0 , I ,.., , 0 0      0 1         = = =       _ I 0 Φ 0 Λ   (57) For all element nodes, the transformation matrix Π is expressed as: ( )1 2 3 4diag , , ,    = (58) For further insights into nodal transformation, refer to Mallek et al. [34]. 3. NUMERICAL TESTS AND DISCUSSION This work investigates geometrically nonlinear performances of FGP shells via the formulated extensible director shell approach. The 3D shell structure is examined for its nonlinear responses taking into account the thickness-stretch. Existing literature lacks such results, necessitating model verification. Through numerical examples, the approach 844 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI proves versatile for analyzing nonlinear structural responses of 3D shell structures, including Porous Functionally Graded materials. Verification tests and simulations validate the methodology's reliability and accuracy. New numerical results are presented showing the influence of porosity dispersion on the large deformation responses of the semi- cylindrical shell structure made by FGP. This research provides crucial insights into the geometrically nonlinear behavior of 3D shell structures with the phenomenon of the thickness-stretch included. 3.1. Linear Behavior of Simply Supported Rectangular Plate The first example is considered as dimensionless. It deals with static bending of a simply supported rectangular plate with the edges a = 7.5 and b = 5 and with the thickness of h = 0.2. The material properties are taken as follows: the elastic modulus E = 1 and the Poisson’s ratio υ = 0.3. The plate is subjected to a uniformly distributed transverse load q = 1, as shown in Fig. 6. The symmetry allows to model only one quarter of the plate. Refinement of the FE mesh was done in order to study the convergence of the proposed model. Fig. 6 Simply supported rectangular plate Table 1 gives a comparison of the central transverse deflection w obtained with the current element and the same results reported in [35]. In [35], a mixed enhanced finite plate element was used. The reference solution used for comparison is 6624.4 which represents an exact 3-D solution given by Piltner [36]. As can be seen, by increasing the number of elements, the results converges and the results by the 32×32 mesh is taken as representative. Table 1 Central deflection results of simply supported plate Mesh 2x2 4x4 8x8 16x16 32 x 32 Piltner and Joseph [35] 6503.6 6591.2 6612.3 6617.9 6619.4 Present 6488.8 6587.7 6613.6 6620.1 6621.8 3.2. Nonlinear Behavior of FGP Square Plate under Uniformly Distributed Load The second example examines deformation of a simply supported square plate under a uniform load. The plate has a width of L=100 mm and a thickness of h=1mm. The material properties include the Poisson’s ratio υ = 0.3 and the Young modulus of aluminum and Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 845 Zr/O2: 𝐸1= 70 GPa and 𝐸2= 151 GPa, respectively. To simplify analysis, the symmetry is used and one quarter of the plate is discretized using 16×16 enhanced four-node elements. The plate surface is subjected to a uniform load q. The parameters considered in this test are the normalized load P= 𝑞𝐿4 𝐸1ℎ4 and the normalized central deflection �̅� = 𝑤/ℎ. Fig. 7 illustrates the normalized load - central deflection curves for various porosity volume fractions, considering a material exponent of n=0.1. The response of the FGP plate closely aligns with findings from Nguyen et al. [37]. Notably, the configurations with evenly distributed pores exhibit higher deflections compared to the uneven configurations due to reduced bending stiffness. Additionally, an increase in porosity volume fraction, particularly with even distribution, results in higher deflections, whereas this effect is less pronounced with uneven distribution due to pores vanishing at the plate's surfaces. a) b) Fig. 7 Results of normalized central deflections of the porous FG plate versus the load parameter using the material exponent n=0.1 at the (a): even; (b) uneven configuration 3.3. Nonlinear Behavior of FGP Cantilever Plate under Uniformly Distributed Load The next test investigates tip deformation of a cantilever plate subjected to a shear line load at the free end. This example shows the significance of the enhanced thickness -5 -4 -3 -2 -1 0 -0.25 -0.2 -0.15 -0.1 -0.05 0 Load parameter N o rm a liz e d c e n tr a l d e fl e c ti o n Present: Perfect Present: Even (=0.2) Present: Even (=0.6) Nguyen et al.[37]: Perfect Nguyen et al.[37]: Even (=0.2) Nguyen et al.[37]: Even (=0.6) -5 -4 -3 -2 -1 0 -0.14 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0 Load parameter N o rm a liz e d c e n tr a l d e fl e c ti o n Present: Perfect Present: Uneven (=0.2) Present: Uneven (=0.6) Nguyen et al.[37]: Perfect Nguyen et al.[37]: Uneven (=0.2) Nguyen et al.[37]: Uneven (=0.6) 846 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI strains �̃�33. The plate with dimensions L=10 mm, b=1 mm, and h=0.1mm is clamped at one end and subjected to a distributed transverse line load of the magnitude -1N/mm at the free end, as shown in Fig. 8. a) b) Fig. 8 Plate’s geometry: a) underformed; b) deformed In this problem, the enhanced four-node element is employed for modeling. The functionally graded material plate is fabricated using transition alloys of metals and ceramics to enhance the ductility of the structure. It is assumed that the top surface is rich in ceramics, while the bottom surface is rich in metals. The material properties of the FGM vary continuously in the thickness direction following a power law distribution. The Poisson’s ratio of the FG material is υ = 0.3 and the Young’s modulus of the ceramic and metal phases are 𝐸𝑐= 3.8 × 105 GPa and 𝐸𝑚= 2.1 × 105 GPa, respectively, as reported by Beheshti and Ramezani [21]. In Fig. 9, the nonlinear tip displacement of the FGM plate is presented using the extensible director shell model, with enhancements of membrane and transverse strains, across various values of the material exponent. It is observed that as the material exponent Fig. 9 Load-deflection curves of the FGM cantilever plate for various values of the material exponent Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 847 n transitions from ceramic to metal, the tip deflection increases, with metallic plates exhibiting greater deflection due to the higher bending stiffness of ceramic plates. Furthermore, the nonlinear response of the FGM cantilever plate closely aligns with the findings reported by Beheshti and Ramezani [21], who employed a seven-parameter shell model to account for the thickness change and utilized the Enhanced Assumed Strain (EAS) method to mitigate the locking phenomena. Hence, the current nonlinear model demonstrates the capability to yield convergent results and avoid locking. Fig. 10 depicts a comparison between the linear and nonlinear deflections of cantilever plates made of ceramic and metal. It is evident that the displacement magnitudes are consistently overestimated in linear analysis. This discrepancy arises because the result of the linear analysis is dominated by the bending stiffness, while the nonlinear behavior correctly accounts for the influence of the membrane stiffness. Furthermore, the nonlinear effect is more pronounced for metallic plates, which exhibit lower stiffness compared to ceramic plates. Fig. 10 Linear and nonlinear deflection curves of cantilever plate made of ceramic and metal 3.4. Nonlinear Behavior of Hinged FGM Cylindrical Panel The next example deals with the cylindrical panel, often referred to as the roof-like panel, depicted in Fig. 11. The example aims to assess the accuracy of the presented enhanced element in nonlinear analysis. The geometrical characteristics of the panel include a length of L=508 mm, a radius of R=2540 mm, a thickness of h=12.7 mm and a span angle of θ = 0.1 rad. The shell is free at its straight edges, hinged over the two curved edges, and subjected to a concentrated normal force F acting in the middle of the panel. The panel is made of a mixture of two material phases with the following mechanical properties: elastic modulus 𝐸𝑚=3102.75 MPa and 𝐸𝑐=3615.77 MPa, for metal and ceramic, respectively, and a Poisson’s ratio υ =0.3. Due to the structure's symmetry, one quarter of the panel is discretized by 4×4 enhanced elements. 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Displacement [mm] L o a d q [ N /m m ] Nonlinear deflection: ceramic Linear deflection: ceramic Nonlinear deflection: metal Linear deflection: metal 848 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI a) b) Fig. 11 (a) Undeformed and (b) deformed shapes of the FGM hinged cylinder The nonlinear deflection at the loaded point, obtained using the developed extensible 3D director shell model, is depicted in Fig. 12 for various values of the material exponent n, and are compared with the results by the models available in the literature [21]. It is evident that the obtained results closely align with the numerical solutions presented by Beheshti and Ramezani [21], particularly in the case of the fully metallic shell, displaying good agreement across different material exponent n. This underscores the high accuracy of the developed FE shell model in predicting large displacements. Furthermore, as illustrated in Fig. 11, the load required to deform a metallic panel is significantly lower than that required for a ceramic panel due to the lower rigidity of the metal. Fig. 12 Load-deflection diagram of the hinged FGM cylindrical panel for various values of the material exponent at loaded point Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 849 3.5. Nonlinear Behavior of FGM Ring Shell under Shear Force at the End The slit annular plate serves as a sensitive benchmark to assess the capability and accuracy of the finite element formulation in handling large deformations. As illustrated in Fig. 13 (a), the structure is subjected to a distributed vertical shear force with a magnitude of q=2120 N/m over the free edge, while a clamped boundary condition is applied at the other edge. The geometric parameters of the plate include the inner radius of r = 6 mm, the outer radius of R =10 m, and the thickness of h = 0.03 m. Regarding material properties, the Young's modulus of the metal and ceramic phases are 𝐸𝑚= 21 × 103MPa and 𝐸𝑐= 38 × 103 MPa, respectively, with a Poisson coefficient ν assumed to be zero. The ring is discretized by 8×46 enhanced elements. Upon application of predefined line load over the free edge, the initially planar structure undergoes arge displacements and rotations, thus deforming into a doubly crooked structure, as depicted in Fig. 13 (b). a) b) Fig. 13 Geometry of the ring plate: a) undeformed shape; b) deformed shape The nonlinear deflection at point A on the outer circle of the functionally graded (FG) ring obtained with the developed model is compared to the results of the models developed by Beheshti and Ramezani [21]. The first model is a displacement-based four-node model, while the second model is an enhanced model incorporating the enhanced assumed strain formulation and accounting for the thickness change. Fig. 14 displays the load–deflection curves at point A Fig. 14 Load-deflection diagram of the FGM ring shell for various values of the material exponent at point A 850 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI for different material exponents n. It is evident from the plotted curves that the present nonlinear deflection of the studied structure closely aligns with the results by the enhanced model of Beheshti and Ramezani [21], although the authors employed a completely different numerical approach based on the strain enhancements and thickness stretch. Hence, this test serves as a highly sensitive benchmark for assessing the accuracy of finite rotation analysis. 3.6. Parametric Study of a Porous Semi-Cylinder In this example, a semi-cylindrical shell structure (Fig. 15 (a)) is considered to analyze its large deflection response subjected to an end pinching force at point A equal to F=2000. Due to the symmetry of the structure, a half of the cylinder is modelled using 40×32 finite elements. The geometrical characteristics of the cylinder are given dimensionless: length L=3.048, radius R=1.016 and thickness h=0.03. In the first step is to validate the nonlinear response of the isotropic semi-cylinder by comparing the result with the solutions by Mallek et al. [2]. The isotropic material considered for the validation has the following material properties: Young’s modulus E =2.068×107 and Poisson’s coefficient υ=0.3. Fig. 16 depicts results of large deflections at the point A of the isotropic semi-cylinder versus the load parameter. As may be seen in this figure, a good accuracy is associated for the present model where present results coincide with the referenced solutions. a) b) Fig. 15 (a) Undeformed and (b) deformed shapes of the porous FGM semi-cylinder Once the nonlinear response of the isotropic semi-cylinder is validated, the subsequent phase involves applying porous Functionally Graded Material (FGM) to the analyzed structure. This FGM comprises a blend of two materials characterized by mechanical properties, 𝐸𝑚= 21 × 109 Pa, 𝜐𝑚=0.3 and 𝐸𝑐= 38 × 109 Pa, 𝜐𝑐=0.3. Fig. 17 illustrates deflection outcomes at point A of the porous semi-cylinder, considering an FGM material exponent (n=0.1) for various porosity distributions: perfect, even, and uneven configurations. This figure demonstrates that the increasing porosity volume fraction α results in significant deflections of the semi-cylinder for both even and uneven configurations. This is due to the reduced bending stiffness caused by a higher density of pores within the structure. Moreover, the even distribution of porosity exhibits greater deflections compared to other configurations with the same porosity volume fraction, especially at higher values of α due to its impact on the bending stiffness of the structure. In an even distribution, pores are uniformly spread throughout the material, affecting its mechanical properties consistently across its volume. This results in a more pronounced reduction in bending stiffness compared to uneven distributions where pores may cluster or be concentrated in specific areas. When the porosity Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 851 is evenly distributed, a larger portion of the material experiences a reduction in stiffness, leading to more significant deformations under loading conditions. Essentially, the even distribution creates more consistent reduction of the structure's bending stiffness, allowing for greater overall deflections. Conversely, in uneven distributions, the concentration of pores in certain regions may lead to localized reduction in stiffness, but other areas may remain relatively unaffected. This can result in a less uniform reduction in bending stiffness throughout the structure, leading to comparatively lower deflections than those observed with even distributions for the same overall porosity volume fraction. Fig. 16 Load-deflection diagram of the isotropic semi-cylinder at point A Furthermore, the influence of the FGM material exponent n on the large deflection response at point A, considering the porosity volume fraction α, is highlighted in Fig. 18 for even and uneven configurations of pores. This figure indicates that the increasing FGM material exponent n leads to higher deflections due to a decrease in bending stiffness resulting from a greater proportion of the metal phase in the structure. Additionally, Fig. 19 depicts the impact of the ratio R/h on the nonlinear bending response at point A of the cylindrical structure, considering the FGM material exponent (n=0.1) for various pore configurations. High values of the ratio R/h correspond to increased deflection of the porous semi-cylinder. This is expected from the geometric influence, flexibility, and reduced bending stiffness associated with the higher R/h ratios. With a larger radius relative to thickness, the shell becomes geometrically less stiff, allowing for greater deformation under applied loads. Additionally, the thinner profile in relation to the curvature enhances flexibility, thus making the structure more susceptible to bending. Consequently, the shell experiences reduced resistance to bending moments, resulting in more significant deflections under the same loading conditions. 0 0.5 1 1.5 2 0 500 1000 1500 2000 Deflection w A L o a d Present Mallek et al. [2] 852 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI a) b) Fig. 17 Results of deflections at point A of the porous semi-cylinder for the material exponent of the FGM n=0.1 for the (a): even; (b) uneven configuration 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 100 200 300 400 500 Deflection w A [m] Lo ad ( × 10 3 [ N ]) Perfect Even (=0.2) Even (=0.6) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 100 200 300 400 500 Deflection w A [m] L o a d ( × 1 0 3 [ N ]) Perfect Uneven (=0.2) Uneven (=0.6) Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 853 a) b) Fig. 18 Results of deflections at point A versus different material exponent n of the porous semi-cylinder for the porosity volume fraction at the (a): even; (b) uneven configuration 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 100 200 300 400 500 Deflection w A [m] L o a d ( × 1 0 3 [ N ]) n=0 n=1 n=3 n=6 0 0.1 0.2 0.3 0.4 0.5 0 100 200 300 400 500 Deflection w A [m] L o a d ( × 1 0 3 [ N ]) n=0 n=1 n=3 n=6 854 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI Fig. 19 Influence of the ration R/h on the large deflections at point A of the porous semi- cylinder considering the material exponent n=0.1 of the FGM for the different configurations of porosity with various values of α Finally, Fig. 20 illustrates the distributions of normal 𝜎𝑥𝑥 and transverse shear stresses 𝜎𝑥𝑧 at point A for different configurations of the porous FGM, considering the material exponent (n=0.1). These stress distributions underscore the influence of porosity volume fractions, as an increase in pore density reduces the structure's stiffness and leads to higher stress values. As evident from Fig. 20 (b), a notable parabolic variation of transverse shear stress throughout the thickness of the structure is observed, owing to the utilization of an enhanced FSDT. This improved theoretical framework provides a more accurate representation of the transverse shear stress distribution within the structure, capturing the complexities of its behavior under loading conditions. 0 10 20 30 40 50 60 0 0.5 1 1.5 2 D e fl e c ti o n w A [ m ] R/h Perfect Even(=0.2) Even(=0.6) Uneven(=0.2) Uneven(=0.6)) Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 855 a) b) Fig. 20 Stress distributions at point A along the thickness of the porous semi-cylinder with a material exponent n=0.1 for (a) 𝜎𝑥𝑥; (b) 𝜎𝑥𝑧 -2 -1 0 1 2 3 -0.5 0 0.5 Stress  xx [×109Pa] z /h Perfect Even(=0.2) Even(=0.6) Uneven(=0.2) Uneven(=0.6) 0 0.5 1 1.5 2 2.5 3 3.5 -0.5 0 0.5 Stress  xz [×108Pa] z /h Perfect Even(=0.2) Even(=0.6) Uneven(=0.2) Uneven(=0.6) 856 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI 4. CONCLUSIONS This study introduces an innovative finite element-based model, incorporating enhancements to the standard FSDT, to investigate large deformation behaviors in porosity- containing Functionally Graded (FG) shell structures. Emphasizing thickness stretching through the enhancement of the third bending strain, with additional parameters facilitated by the EAS method, and augmentation of the membrane strains with enriched parameters, a comprehensive six-degree-of-freedom shell model is developed. The investigation evaluates the model's performance through simulations of various boundary conditions and different types of perfect FGMs, demonstrating convergent and accurate results compared to existing literature. A parametric study explores the impact of porosities on FGM characteristics, particularly focusing on cylindrical shell structures. The presence of pores reduces the structure's stiffness, resulting in significant deformations, particularly notable in unevenly distributed porosity configurations. Moreover, an increase in porosity volume fraction exacerbates the deflection of porous FGM shell structures, underscoring the structural implications of porosity in FGM materials. This research significantly advances the understanding of PFG shell structures, offering valuable insights into their nonlinear behavior under diverse loading conditions. Acknowledgement: This research is carried out under Young Researchers Encouragement Program, funded by the Ministry of Higher Education and Scientific Research (Tunisia) through project number PEJC2023-D6P19. REFERENCES 1. Hajlaoui, A., Dammak, F., 2022, A modified first shear deformation theory for three-dimensional thermal post- buckling analysis of FGM plates, Meccanica, 57(2), pp. 337-353. 2. Mallek, H., Jrad, H., Algahtani, A., Wali, M., Dammak, F., 2019, Geometrically non-linear analysis of FG-CNTRC shell structures with surface-bonded piezoelectric layers, Computer Methods in Applied Mechanics and Engineering, 347, pp. 679-699. 3. Mota, A. F., Loja, M. A. R., 2019, Mechanical behavior of porous functionally graded nanocomposite materials, Journal of Carbon Research C, 5(2), 34. 4. Zghal, S., Dammak, F., 2021, Buckling responses of porous structural components with gradient power-based and sigmoid material variations under different types of compression loads, Composite Structures, 273, 114313. 5. Wattanasakulpong, N., Prusty, B.G., Kelly, D.W., Hoffman, M., 2012, Free vibration analysis of layered functionally graded beams with experimental validation, Materials & Design (1980-2015), 36, pp. 182-190. 6. Wattanasakulpong, N., Ungbhakorn, V., 2014, Linear and nonlinear vibration analysis of elastically restrained ends FGM beams with porosities, Aerospace Science and Technology, 32(1), pp. 111-120. 7. Wattanasakulpong, N., Chaikittiratana, A., 2015, Flexural vibration of imperfect functionally graded beams based on Timoshenko beam theory: Chebyshev collocation method, Meccanica, 50, pp. 1331-1342. 8. Wang, Y., Wu, D., 2017, Free vibration of functionally graded porous cylindrical shell using a sinusoidal shear deformation theory, Aerospace Science and Technology, 66, pp. 83-91. 9. Nguyen, N. V., Nguyen-Xuan, H., Lee, D., Lee, J., 2020, A novel computational approach to functionally graded porous plates with graphene platelets reinforcement, Thin-Walled Structures, 150, pp. 106684. 10. Nguyen, N. V., Nguyen-Xuan, H., Lee, J., 2022, A quasi-three-dimensional isogeometric model for porous sandwich functionally graded plates reinforced with graphene nanoplatelets, Journal of Sandwich Structures & Materials, 24(2), pp. 825-859. 11. Thang, P. T., Nguyen-Thoi, T., Lee, D., Kang, J., Lee, J., 2018, Elastic buckling and free vibration analyses of porous- cellular plates with uniform and non-uniform porosity distributions, Aerospace science and technology, 79, pp. 278-287. 12. Amir, M., Talha, M., 2019, Nonlinear vibration characteristics of shear deformable functionally graded curved panels with porosity including temperature effects, International Journal of Pressure Vessels and Piping, 172, pp. 28-41. Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 857 13. Zghal, S., Dammak, F., 2021, Buckling responses of porous structural components with gradient power-based and sigmoid material variations under different types of compression loads, Composite Structures, 273, 114313. 14. Zghal, S., Dammak, F., 2021, Vibration characteristics of plates and shells with functionally graded pores imperfections using an enhanced finite shell element, Computers & Mathematics with Applications, 99, pp. 52-72. 15. Phuong, N.T.B., Tu, T.M., Phuong, H. T., Van Long, N., 2019, Bending analysis of functionally graded beam with porosities resting on elastic foundation based on neutral surface position, Journal of Science and Technology in Civil Engineering (JSTCE)-HUCE, 13(1), pp. 33-45. 16. Binh, C. T., Van Long, N., Tu, T. M., Minh, P. Q., 2021, Nonlinear vibration of functionally graded porous variable thickness toroidal shell segments surrounded by elastic medium including the thermal effect, Composite Structures, 255, 112891. 17. Milić P., Marinković D., Klinge S., Ćojbašić Ž., 2023, Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells, Tehnicki Vjesnik, 30(2), pp. 416-425. 18. Rama G., Marinkovic D., Zehn M., 2018, High performance 3-node shell element for linear and geometrically nonlinear analysis of composite laminates, Composites Part B: Engineering, 151, pp. 118-126. 19. Abid, M., Koubaa, S., Abdelkefi, A., Frikha, A., Dammak, F., 2021, Numerical modeling of porous functionally graded shells response in large deflection, Computers & Mathematics with Applications, 104, pp. 59-70. 20. Simo, J. C., Rifai, M., 1990, A class of mixed assumed strain methods and the method of incompatible modes, International journal for numerical methods in engineering, 29(8), pp. 1595-1638. 21. Beheshti, A., Ramezani, S., 2015, Nonlinear finite element analysis of functionally graded structures by enhanced assumed strain shell elements, Applied Mathematical Modelling, 39(13), pp. 3690-3703. 22. Nguyen, N. V., Tran, K. Q., Lee, J., Nguyen-Xuan, H., 2024, Nonlocal strain gradient-based isogeometric analysis of graphene platelets-reinforced functionally graded triply periodic minimal surface nanoplates, Applied Mathematics and Computation, 466, 128461. 23. Tran, K. Q., Hoang, T. D., Lee, J., Nguyen-Xuan, H., 2024, Three novel computational modeling frameworks of 3D- printed graphene platelets reinforced functionally graded triply periodic minimal surface (GPLR-FG-TPMS) plates, Applied Mathematical Modelling, 126, pp. 667-697. 24. Fotiu, P. A., Murin, J., 2021, A novel GBT-formulation for thin-walled FGM-beam-structures based on a reference beam problem, Composite Structures, 257, 113158. 25. Wattanasakulpong, N., Ungbhakorn, V., 2013, Linear and nonlinear vibration analysis of elastically restrained ends FGM beams with porosities, Aero. Sci. Technol. 32, pp. 111–120. 26. Kugler, P. Fotiu, J. Murin, 2013, The numerical analysis of FME shells with enhanced finite elements, Engineering Structures, 49, pp. 920–935. 27. Mellouli, H., Mallek, H., Wali, M., Dammak, F., Gamaoun, F., Abdulrahman, A., 2022, An extensible double director 3D shell formulation for FGM-CNTRC shell bending analysis, Engineering Analysis with Boundary Elements, 145, pp. 258-270. 28. Tanov, R., Tabiei, A., 2000, A simple correction to the first-order shear deformation shell finite element formulations, Finite elements in analysis and design, 35(2), pp. 189-197. 29. Shi, G. , 2007, A new simple third-order shear deformation theory of plates, Int J Solids Struct, 44, pp. 4399-4417. 30. Simo, J.C., Rifai, M.S., Fox, D., 1990, On a stress resultant geometrically exact shell model. Part IV: Variable thickness shells with through-the-thickness stretching, Computer methods in applied mechanics and engineering, 81(1), pp. 91-126. 31. Bischoff, M., Ramm, E., 1997, Shear deformable shell elements for large strains and rotations, International Journal for Numerical Methods in Engineering, 40(23), pp. 4427-4449. 32. Hughes, T. J., Tezduyar, T., 1981, Finite elements based upon Mindlin plate theory with particular reference to the four-node bilinear isoparametric element, Journal of Applied Mechanics, pp. 587-596. 33. Dvorkin, E. N., Bathe, K. J., 1984, A continuum mechanics based four‐node shell element for general non‐linear analysis, Engineering computations, 1(1), pp. 77-88. 34. Mallek, H., Mellouli, H., Said, L.B., Wali, M., Dammak, F., Boujelbene, M., 2023, Bending and free vibration analyses of CNTRC shell structures considering agglomeration effects with through-the-thickness stretch, Thin- Walled Structures, 191, 111036. 35. Piltner, R., Joseph, D.S., 2001, A mixed finite element for plate bending with eight enhanced strain modes, Communications in numerical methods in engineering, 17(7), pp. 443-454. 36. Piltner, R., 1988, The application of a complex 3-dimensional elasticity solution representation for the analysis of a thick rectangular plate, Acta Mechanica, 75(1), pp. 77-91. 37. Nguyen, N.V., Nguyen, H.X., Lee, S., Nguyen-Xuan, H., 2018, Geometrically nonlinear polygonal finite element analysis of functionally graded porous plates, Advances in Engineering software, 126, pp. 110-126. 858 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI APPENDIX A: STRAIN-DISPLACEMENT INTERPOLATION Membrane : .m n =e B Φ I I I m mu m  =  B B 0 B T 1 ,1 T 2 ,2 T T 1 ,2 2 ,1 I I I mu I I N N N N+      =        n n B 0 n n 0 0 ( 1) 0 I m I I N        =  −     B Bending : .b n =χ B Φ I I I I b bu bt b  =  B B B B ,1 ,1 ,2 ,2 ,1 ,2 ,2 ,1 T I T I I bu T I T I N N N N+      =        d d B 0 d d 1 1 2 2 2 1 1 2 I T I T I bt I T I T L L L L+      =        n n B 0 n n 1 1 2 2 1 2 2 1 . . 0 . . T I T I I b T I T I  +      =        n m n m B n m n m , , ,(1 )   = − −d t t , ,(1 )I I IL N N   = − − , ,( )I I IN N  = − +m t t Transverse Shear, First Terme : 0 0 .s n =γ B Φ 0 1 0 s  −=B J B 1 2 2 2 2 2 ,1 1 1 ,1 1 ,1 1 1 ,1 10 1 4 4 2 3 3 ,2 2 2 ,2 2 ,2 2 2 ,2 2 . . . . T T T T T T B B B B B B B B B B T T T T T T A A A A A C C C C C N P N N P N N P N N P N   − − =  − − d a a t d a a t B d a a t d a a t 3 3 3 4 3 3 ,1 1 1 ,1 1 ,1 1 1 ,1 1 3 3 3 4 4 4 ,2 2 2 ,2 2 ,2 2 2 ,2 2 . . . . T T T T T T D D D D D D D D D D T T T T T T C C C C C A A A A A N P N N P N N P N N P N − −  − −  d a a t d a a t d a a t d a a t , (1 )I I M MP N  = − Transverse Shear, second Terme : 1 1.s n =γ B Φ 1 1 1 s  −=B J B 1 1I I    =  B 0 0 B 2 2 3 3 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,11 4 3 3 4 ,2 ,2 ,2 ,2 ,2 ,2 ,2 ,2 B B D D A C C A N N N N N N N N            =      B 1 2 3 4 ,1 ,1 ,1 ,1 1 2 3 4 ,2 ,2 ,2 ,2 ( 1) ( 1) ( 1) ( 1) ( 1) ( 1) ( 1) ( 1) B B D D A C C A N N N N N N N N          − − − − +   − − − −   Porosity Effects on Nonlinear Static Performances of FG Shells Considering Thickness Stretching 859 APPENDIX B: THE INTERPOLATION MATRIX FOR ENHANCEMENT OF MEMBRANE AND TRANSVERSE NORMAL STRAIN Transverse normal strain with one additional parameter : 6 1 0 5 12 1 / 1A A      =       0 B 0 Transverse normal strain with four additional parameters : 6 4 4 1 4 5 4 12 4 T         =       0 B V 0 1 1 4 0 2 1 2 1 /A A           =       V Enhanced membrane and transverse normal strain with one additional parameter : 4 4 5 1 4 4 4 4 12 5         =         G B B 0 0 Enhanced membrane and transverse normal strain with four additional parameters : 4 4 8 4 4 4 4 4 12 8         =         G B B 0 0 APPENDIX C: GEOMETRIC TANGENT MATRIX COMPONENTS 6 6 ( ) ( ) ( ) IJ IJ IJ IJ IJ GIJ JI JI IJ IJ IJ IJ T T T T JI JI JI JI IJ UUM UBF UBC UBF UBC BBF BBC LL  + +    = + + +    + +  I I ULF ULC K I I BLF BLC ULF ULC BLF BLC Components of GIJK : ( )11 12 12 22 ,1 ,1 ,2 ,2 ,1 ,2( ) ( )I J J I J J IJ A UUM N N N N N N N N N N dA= + + + ( )11 12 12 22 ,1 1 2 ,2 1 2( ) ( )I J J I J J IJ A UBF N M L M L N M L M L dA= + + + ( )11 12 12 22 ,1 1 2 ,2 1 2 0 ( ) ( ) IJ I I II I A BBF if I J BBF L M M L M M dA =   = − + + +  n n n n t ( )11 12 12 22 ,1 1 2 ,2 1 2( ) ( )I J J I J J IJ A N M M N M M dA= + + +ULF m m m m 860 H. MALLEK, H. MELLOULI, L. BEN SAID, M. WALI, F. DAMMAK, M. ALHADRI ( )11 12 12 22 1 1 2 2 1 1 2 2( ) ( )IJ IJ IJ IJ IJ A M n M n M n M n dA= − + + +BLF n n n n 33 0IJ I II A LL if I J LL N N dA =   =  , ,(1 )I I IL N N   = − − , ,( )I I IN N  = − +m t t , , IJ I J I Jn N N N N  = + 0 01 08 0 A B B A B B C C C C D D A D A D − − − −    − −  =  +   − −  UBC 1 8 − − − −    − −  =  +   − −  a b b 0 a b b c c 0 ULC 0 c c d d a 0 d a d 1 8 +    +  = −  +   +  α β β 0 α β β γ γ 0 BLC 0 γ γ δ δ α 0 δ α δ  1 1 2 2 3 3 4 4 1 Diag . . . . 4 = −BBC g d I g d I g d I g d I 2 0(1 ) A a T dA= − 1 0(1 ) A b T dA= − 2 0(1 ) A c T dA= + 1 0(1 ) A d T dA= + (1 )AA a= − (1 )BB b= − (1 )CC c= − (1 )DD d= − Aa=a t Bb=b t Cc=c t Dd=d t 2(1 )A A A a = −g a 1(1 )B B B b = −g a 2(1 )C C C c = −g a 1(1 )D D D d = −g a 1 A B= +g g g 2 B C= +g g g 3 C D= +g g g 4 D A= +g g g 2 Aa=α a 1 Bb=β a 2 Cc=γ a 1 Dd=δ a