Accurate numerical solutions for the static analysis of functionally graded porous deep curved beams under various boundary conditions

Abstract

This article presents a numerical solution for the static analysis of functionally graded porous deep curved beams with various boundary conditions. The main objective is to develop an accurate reference model for structures with high curvature and arbitrary support conditions, where simplified theories—such as higher-order shear deformation models—exhibit limited accuracy. The strong form of the governing equations is derived from two-dimensional elasticity theory for curved beams and solved using the Differential Quadrature Method, which ensures interlaminar continuity of displacements and stresses as well as boundary condition equations. Numerical validation is performed by comparison against results reported in the literature and those obtained using commercial finite element software. A parametric study is performed to investigate the influence of geometric curvature, material gradation, and porosity on the structural response. The results demonstrate the model’s capability to yield accurate solutions applicable to the design of curved structural components made of functionally graded materials.

Introduction

The structural analysis of functionally graded materials (FGMs) represents a key challenge in the design of advanced components. These structures are widely employed in sectors such as aerospace, biomedical, and automotive engineering, where high strength, complex geometries, and spatially varying mechanical properties are required. ¹ For curved or porous configurations, accurately modeling their structural response demands robust numerical formulations capable of capturing the continuous variation of properties and the actual boundary conditions, beyond the simplified cases traditionally assumed.

Several researchers have developed analytical and numerical models to address the complex behavior of functionally graded structures. Filippi et al. ² applied the one-dimensional Carrera Unified Formulation (CUF) to the static analysis of FG structures, employing various displacement functions and solving the governing equations using the finite element method. Mantari and Granados ³ proposed a first-order shear deformation theory (FSDT) with only four unknowns for analyzing FG plates, obtaining closed-form solutions using the Navier series. Subsequently, Mantari et al. ⁴ extended this approach through the Carrera Unified Formulation (CUF) and new non-polynomial displacement fields, employing trigonometric, exponential, and hyperbolic functions to analyze simple and sandwich plates under static loading. Brischetto ⁵ developed an exact three-dimensional model for the static analysis of simply supported structures made of FGM layers and subjected to harmonic normal loads. In a subsequent extension, Brischetto ⁶ formulated a general three-dimensional model based on layer-wise elasticity theory, employing a closed-form solution using the exponential matrix method. This model enables the application of normal and shear loads at various positions, providing precise solutions for different geometries, material configurations, and loading conditions. Sankar ⁷ proposed an exact elasticity solution for beams with an exponentially varying Young’s modulus through the thickness, complemented by an Euler–Bernoulli-type formulation for slender beams. He found that stress concentrations depend on the material distribution, decreasing when the more compliant side is subjected to loading. Zhong and Yu ⁸ developed a general two-dimensional solution based on the Airy stress function for beams with thickness-wise varying properties, yielding analytical expressions that are useful for analyzing various boundary and loading conditions. Kadoli et al. ⁹ employed a higher-order shear deformation theory and the finite element method to investigate metal-ceramic FG beams, revealing material gradation to be a significant parameter for static analysis. Thai and Vo ¹⁰ formulated several higher-order shear deformation models for bending and free vibration analysis without the need for a shear correction factor. Nguyen et al. ¹¹ developed a first-order theory with an analytically derived shear correction factor, providing analytical solutions for simply supported beams under axial loading. Vo et al. ¹² extended this approach by incorporating axial–bending coupling and implementing cubic Hermite elements, evaluating the influence of the power-law exponent and material anisotropy under various boundary conditions. Wang and Liu ¹³ applied the Airy stress function to curved beams composed of orthotropic FG layers, solving the problem using generalized hypergeometric functions for various boundary conditions. Pydah and Sabale ¹⁴ investigated the bending behavior of bidirectional FG circular beams using analytical solutions based on Euler–Bernoulli beam theory. Yarasca et al. ¹⁵ analyzed the static behavior of FG beams using a quasi-3D hybrid theory with seven degrees of freedom, solving the equations through finite element analysis with cubic and linear Hermite interpolation. Sayyad and Ghugal ¹⁶ analyzed FG curved sandwich beams using a sinusoidal theory to obtain analytical solutions under simply supported conditions. Garg et al. ¹⁷ evaluated the behavior of FG sandwich beams under combined hygrothermal and mechanical loading, considering temperature- and moisture-dependent properties through a higher-order zigzag theory (HOZT) finite element model. Belarbi et al. ¹⁸ proposed a new refined shear deformation theory that accounts for parabolic shear stress distribution in FG curved sandwich beams, which was solved numerically using a two-node finite element. Excellent reviews on the modeling and analysis of FG structures, including the effects of porosity, were presented by Sayyad and Ghugal, ¹⁹ Swaminathan and Sangeetha, ²⁰ and Garg et al. ²¹ .

Functionally graded porous materials have emerged as an efficient alternative in advanced structural design due to their ability to combine continuous property gradation with density control through the distribution of porosity. This feature allows for the tailoring of mechanical and thermal properties, making them suitable for applications in lightweight, biomedical, and aerospace structures. Chen et al. ²² presented a recent review of this type of material, addressing advances in mechanical analysis techniques, manufacturing methods, and the emerging use of artificial intelligence for modeling and experimental evaluation. In the field of structural analysis, several studies have investigated the impact of porosity on the mechanical behavior of FG structures. Chen et al. ²³ investigated the elastic buckling and static bending of porous FG beams using the Timoshenko beam theory and the Ritz method, considering various porosity patterns and boundary conditions. Demirhan and Taskin ²⁴ investigated porous FG plates using a four-variable theory, which accounts for the effect of porosity on bending and free vibration behavior. Masjedi et al. ²⁵ analyzed the nonlinear behavior of porous FG beams subjected to large deflections using an intrinsic formulation and the Chebyshev collocation method. Their study included various porosity distributions (transverse and longitudinal) and demonstrated the effectiveness of the proposed model in capturing large deformations. Similarly, Zghal et al. ²⁶ employed a refined mixed finite element model to investigate the effect of porosity on the static bending of FG beams, considering different porosity distributions, power-law indices, and boundary conditions. Wang et al. ²⁷ developed a higher-order beam model for the dynamic analysis of porous FG sandwich beams under the influence of a moving distributed mass. Their results showed that the mass location has a significant effect on dynamic deflections. Sayyad et al. ²⁸ investigated the static deformation and free vibration of porous FG circular beams using a higher-order hyperbolic theory and Navier-type solutions. This study incorporated both shear deformation and rotary inertia effects, considering uniform and non-uniform porosity distributions.

The differential quadrature method (DQM), introduced by Bellman et al., ²⁹ emerged as an efficient numerical technique for solving nonlinear partial differential equations. DQM has been widely adopted in structural analysis, particularly for composite material structures, where obtaining closed-form solutions is challenging due to through-thickness property variations and geometric complexity. Brischetto and Tornabene ³⁰ applied the Generalized Differential Quadrature (GDQ) method in a comparative study between 2D and 3D models for the static analysis of multilayered and sandwich plates and shells subjected to transverse loads. While the 3D models were solved using exact solutions, the 2D models employed the GDQ method, demonstrating accuracy for structures with complex geometry, including double-curvature shells. In a complementary work, Tornabene and Brischetto ³¹ extended this approach by developing both closed-form and numerical 3D solutions using the exponential matrix and GDQ methods. The 2D solutions were based on the Unified Formulation and incorporated both classical and refined theories. Monge and Mantari ³² developed a three-dimensional numerical solution for the static analysis of FG shells. They employed the Navier method and the DQM in curvilinear coordinates to accurately capture interlaminar continuity effects and transverse stresses. Subsequently, Padilla et al. ³³ extended the semi-analytical formulation and the numerical solution using DQM presented in ³² to the static analysis of laminated arch beams in elevation with deep curvature, employing a two-dimensional DQM solution.

In recent years, several studies have addressed the modeling of advanced composite materials with functionally graded porous structures. Among them, Tounsi et al. ³⁴ and Belabed et al.^35,36 have employed higher-order shear deformation theories (HSDT) and quasi-3D theories, together with finite element methods, to analyze the static and dynamic behavior of porous nanocomposite beams reinforced with carbon nanotubes (CNTs) on elastic foundations. For their part, Mohamed et al. ³⁷ investigated the postbuckling behavior and snap-through instability in bidirectional functionally graded (BDFG) porous plates resting on elastic foundations, using a parabolic shear deformation theory and solving the system through quadrature methods.

In parallel, other recent works have focused on the analysis of functionally graded (FG) sandwich beams using higher-order shear deformation theories and quasi-3D theories. Draiche et al. ³⁸ proposed a high-order model for the static analysis of FG sandwich curved beams, accurately capturing shear stress without the need for correction factors. Addou et al. ³⁹ developed a refined hyperbolic theory for the static analysis of FG sandwich beams under different boundary conditions, whereas Mohamed et al.^40,41 employed quasi-3D theories to study the static and dynamic behavior of FG sandwich beams with porous core on elastic foundations, considering different porosity patterns and gradations of material properties.

Although formulations based on HSDT or quasi-3D theories, widely reported in the literature, have proven efficient in various applications, their accuracy decreases when applied to beams with deep curvatures, materials with high gradation, or complex boundary conditions. These limitations motivate the development of a two-dimensional elasticity-based numerical formulation for the static analysis of functionally graded porous curved beams, capable of delivering high-accuracy solutions for beams with deep curvatures and arbitrary boundary conditions. Unlike simplified approaches, the proposed model provides highly accurate solutions that can serve both as reference benchmarks for validating other models and as a basis for the structural design of curved components made of functionally graded materials.

To this end, the approach is based on a two-dimensional elasticity formulation derived through the principle of virtual displacements, following the methodology reported in. ³³ The governing equations are numerically solved using the two-dimensional Differential Quadrature Method, which allows for the consideration of arbitrary boundary conditions. The results obtained enable the assessment of the influence of key parameters, such as curvature, material gradation, and porosity fraction, on the structural response of functionally graded porous curved beams.

The present article is organized as follows: Section 2 provides a detailed description of the equilibrium equations and their associated boundary conditions; Section 3 formulates the numerical solution using the DQM; Section 4 discusses the numerical results obtained for functionally graded porous curved beams under various boundary conditions; and finally, Section 5 presents the main conclusions of the study.

Two-dimensional curved beam model

Geometry

Consider a functionally graded curved beam with length $L$ , thickness $h$ , and radius of curvature $R$ , as illustrated in Figure 1. The beam is described using the orthogonal coordinates ( $\alpha$ , $z$ ), where $\alpha$ is the curvilinear coordinate along the length of the beam, and $z$ is taken along the thickness direction. The mathematical domain for the curved beam is defined as $0\le \alpha \le L$ and $-h/2\le z\le h/2$ . The total number of layers is denoted by $\mathrm{NL}$ , and $z_{k-1}$ and $z_k$ represent the lower and upper thickness coordinates of the $k$ -th layer, respectively.

Figure 1.

Geometry and reference coordinate system of a single-layer functionally graded metal–ceramic curved beam with uniform porosity.

Functionally graded porous material

In this work, single-layer functionally graded porous metal–ceramic beams with continuous variation of material properties along the thickness direction are analyzed (see Figure 1). The FGM varying properties are described by the distribution law $P\left(z\right)$ which changes through the thickness. Nikbakht et al. ⁴² presented a review of the most common distribution laws reported in the literature. In the present study, a power-law distribution that accounts for porosity is employed.

The material properties $P\left(z\right)$ for a metal–ceramic FGM following a power-law distribution and assuming a uniform porosity distribution are defined as ²⁸ :

$$\begin{array}{c}P\left(z\right)=(P_c-P_m)V_c\left(z\right)+P_m-\frac{a}{2}(P_c+P_m),\hfill \\ V_c\left(z\right)={\left(\frac{2z+h}{2h}\right)}^p\hfill \end{array}$$

(1)

where $V_c\left(z\right)$ is the volume fraction of the ceramic phase, $p$ is the power-law index that defines the property gradation across the functionally graded structure, and $a$ denotes the porosity volume fraction ⁴² $P_c$ and $P_m$ are the material properties of the ceramic and metallic phases, respectively. Figure 2 shows the influence of the power-law index on the distribution of Young’s modulus through the thickness of the beam.

Figure 2.

Distribution of Young’s modulus through the thickness of the beam for various values of the power-law index $p$ ( $a=0$ ).

Stress-displacement relations

For a general case of a multilayer curved beam with deep curvature, the two-dimensional strain–displacement relations are expressed as follows:

$$\begin{array}{c}{\epsilon }_{\alpha }^k=\frac{1}{H^k}\frac{\partial u^k}{\partial \alpha }+\frac{1}{R{H}^k}{w}^k{\epsilon }_z^k=\frac{\partial w^k}{\partial z}{\gamma }_{\alpha z}^k=\frac{1}{H^k}\frac{\partial w^k}{\partial \alpha }\hfill \\ +\frac{\partial u^k}{\partial z}-\frac{1}{R{H}^k}{u}^k\hfill \end{array}$$

(2)

where $k$ denotes the layer index, and $H$ is a parametric coefficient defined by the following equation: $H=1+z/R$ . It is important to note that Eq. (2) is derived from the three-dimensional strain–displacement relations for shells, ⁴³ preserving the kinematic terms associated with displacements along the $\alpha$ and $z$ coordinates.

For multilayer two-dimensional beams, the stress–strain relations for the functionally graded layer $k$ are determined by the generalized Hooke’s law:

$$\left\{\begin{array}{c}{\sigma }_{\alpha }\hfill \\ \multicolumn{1}{c}{{\sigma }_z}\\ \multicolumn{1}{c}{{\tau }_{\alpha z}}\end{array}\right\}{}^{\left(k\right)}={\left[\begin{array}{ccc}{C}_{11}\left(z\right)\hfill & C_{13}\left(z\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{13}\left(z\right)}& C_{33}\left(z\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & C_{55}\left(z\right)\hfill \end{array}\right]}^{\left(k\right)}\left\{\begin{array}{c}{\epsilon }_{\alpha }\hfill \\ \multicolumn{1}{c}{{\epsilon }_z}\\ \multicolumn{1}{c}{{\gamma }_{\alpha z}}\end{array}\right\}{}^{\left(k\right)}$$

(3)

where $C_{\mathrm{ij}}$ represents the transformed elastic coefficients, which are derived from the general three-dimensional state by considering either a plane stress ( ${\sigma }_y={\tau }_{\mathrm{yz}}={\tau }_{\alpha y}=0$ ) or plane strain ( ${\epsilon }_y={\gamma }_{\mathrm{yz}}={\gamma }_{\alpha y}=0$ ) condition. ⁴⁴

For the plane stress condition, the coefficients $C_{\mathrm{ij}}$ are defined as follows:

$$\begin{array}{c}{C}_{11}\left(z\right)=C_{33}\left(z\right)=\frac{E\left(z\right)}{1-v^2\left(z\right)},C_{13}\left(z\right)=\frac{E\left(z\right)v\left(z\right)}{1-v^2\left(z\right)},\hfill \\ C_{55}\left(z\right)=\frac{E\left(z\right)}{2(1+v\left(z\right))}\hfill \end{array}$$

(4)

whereas, for a plane strain condition the coefficients $C_{\mathrm{ij}}$ are expressed as:

$$\begin{array}{c}{C}_{11}\left(z\right)=C_{33}\left(z\right)=\frac{E\left(z\right)(1-v\left(z\right))}{(1-2v\left(z\right))(1+v\left(z\right))},\hfill \\ C_{13}\left(z\right)=\frac{E\left(z\right)v\left(z\right)}{(1-2v\left(z\right))(1+v\left(z\right))},C_{55}=\frac{E\left(z\right)}{2(1+v\left(z\right))}\hfill \end{array}$$

(5)

The stress–displacement relationships for the curved beam are obtained by substituting equation (2) into equation (3):

$$\begin{array}{c}{\sigma }_{\alpha }^k={D_1}^k\frac{\partial u^k}{\partial \alpha }+{D_2}^k{w}^k+{D_3}^k\frac{\partial w^k}{\partial z}\hfill \\ {\sigma }_z^k={F_1}^k\frac{\partial u^k}{\partial \alpha }+{F_2}^k{w}^k+{F_3}^k\frac{\partial w^k}{\partial z}\hfill \\ {\tau }_{\alpha z}^k={G_1}^k{u}^k+{G_2}^k\frac{\partial u^k}{\partial z}+{G_3}^k\frac{\partial w^k}{\partial \alpha }\hfill \end{array}$$

(6)

where:

$$\begin{array}{c}({D_1}^k,{D_2}^k,{D_3}^k)=(\frac{C_{11}^k}{H^k},\frac{C_{11}^k}{R{H}^k},C_{13}^k)\hfill \\ ({F_1}^k,{F_2}^k,{F_3}^k)=(\frac{C_{13}^k}{H^k},\frac{C_{13}^k}{R{H}^k},C_{33}^k)\hfill \\ ({G_1}^k,{G_2}^k,{G_3}^k)=(-\frac{C_{55}^k}{R{H}^k},C_{55}^k,\frac{C_{55}^k}{H^k})\hfill \end{array}$$

(7)

Two-dimensional governing equations

The governing equations for multilayer curved beams are derived from the principle of virtual displacements, as presented in Padilla et al. ³³ . The variational form of the governing equations is expressed as follows:

$$\begin{array}{c}\sum _{k=1}^{\mathrm{NL}}\underset{0}{\overset{L}{\int }}\underset{z_{k-1}}{\overset{z_k}{\int }}\delta u^k(-\frac{\partial \left({\sigma }_{\alpha }^k\right)}{\partial \alpha }-\frac{\partial \left({\tau }_{\alpha z}^k{H}^k\right)}{\partial z}-\frac{{\tau }_{\alpha z}^k}{R})\hfill \\ +\delta w^k(\frac{{\sigma }_{\alpha }^k}{R}-\frac{\partial \left({\sigma }_z^k{H}^k\right)}{\partial z}-\frac{\partial \left({\tau }_{\alpha z}^k\right)}{\partial \alpha })\mathrm{dz}d\alpha =0\hfill \end{array}$$

(8)

The governing equations in terms of stresses are then defined by collecting the coefficients of the virtual displacements:

$$\begin{array}{c}\delta u^k:\frac{\partial \left({\sigma }_{\alpha }^k\right)}{\partial \alpha }+\frac{\partial \left({\tau }_{\alpha z}^k{H}^k\right)}{\partial z}+\frac{{\tau }_{\alpha z}^k}{R}=0\hfill \\ \delta w^k:\frac{\partial \left({\tau }_{\alpha z}^k\right)}{\partial \alpha }+\frac{\partial \left({\sigma }_z^k{H}^k\right)}{\partial z}-\frac{{\sigma }_{\alpha }^k}{R}=0\hfill \end{array}$$

(9)

The first equation enforces axial equilibrium, where the variation of the axial normal stress ${\sigma }_{\alpha }^k$ along the arc, the gradient of the shear stress ${\tau }_{\alpha z}^k$ through the thickness (scaled by $H^k$ ), and the curvature effect ${\tau }_{\alpha z}^k/R$ must balance. The second equation enforces transverse equilibrium, combining the variation of the shear stress along the arc, the gradient of the transverse normal stress ${\sigma }_z^k$ through the thickness (scaled by $H^k$ ), and the curvature contribution $-{\sigma }_{\alpha }^k/R$ . Together, these equations ensure that the axial and transverse stress components interact consistently under the combined influence of material gradation and beam curvature.

Finally, the governing equations in terms of displacements are obtained by substituting equation (6) into equation (9):

$$\begin{array}{c}{A_1}^k{u}^k+{A_2}^k\frac{\partial u^k}{\partial z}+{A_3}^k\frac{{\partial }^2{u}^k}{\partial {\alpha }^2}+{A_4}^k\frac{{\partial }^2{u}^k}{\partial z^2}+{A_5}^k\frac{\partial w^k}{\partial \alpha }\hfill \\ +{A_6}^k\frac{{\partial }^2{w}^k}{\partial \alpha \partial z}=0\hfill \\ {B_1}^k\frac{\partial u^k}{\partial \alpha }+{B_2}^k\frac{{\partial }^2{u}^k}{\partial \alpha \partial z}+{B_3}^k{w}^k+{B_4}^k\frac{\partial w^k}{\partial z}+{B_5}^k\frac{{\partial }^2{w}^k}{\partial {\alpha }^2}\hfill \\ +{B_6}^k\frac{{\partial }^2{w}^k}{\partial z^2}=0\hfill \end{array}$$

(10)

The coefficients ${A_i}^k$ and ${B_i}^k$ do not directly represent physical quantities; instead, they result from the mathematical formulation obtained by combining the stress–displacement relations with the governing equations.

where:

$$\begin{array}{l}\left({A_1}^k,{A_2}^k,{A_3}^k,{A_4}^k,{A_5}^k,{A_6}^k\right)=\left(-\frac{1}{R}\frac{\partial C_{55}^k}{\partial z}-\frac{C_{55}^k}{R^2{H}^k},H^k\frac{\partial C_{55}^k}{\partial z}+\frac{C_{55}^k}{R},\frac{C_{11}^k}{H^k},H^k{C}_{55}^k,\frac{C_{11}^k}{R{H}^k}+\frac{\partial C_{55}^k}{\partial z}+\frac{C_{55}^k}{R{H}^k},C_{13}^k+C_{55}^k\right)\\ \left({B_1}^k,{B_2}^k,{B_3}^k,{B_4}^k,{B_5}^k,{B_6}^k\right)=\left(-\frac{C_{55}^k}{R{H}^k}+\frac{\partial C_{13}^k}{\partial z}-\frac{C_{11}^k}{R{H}^k},C_{55}^k+C_{13}^k,\frac{1}{R}\frac{\partial C_{13}^k}{\partial z}-\frac{C_{11}^k}{R^2{H}^k},\frac{C_{33}^k}{R}+H^k\frac{\partial C_{33}^k}{\partial z},\frac{C_{55}^k}{H^k},H^k{C}_{33}^k\right)\end{array}$$

(11)

Boundary conditions

Boundary conditions equations for the beam ends, upper and lower surfaces, and interlaminar continuity conditions for multilayered beams are defined.

The boundary conditions at the beam ends ( $\tilde{\alpha }=0$ and $\tilde{\alpha }=L$ ) are expressed as:

Free (F):

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\alpha }^k(\tilde{\alpha },z)=\left({D_1}^k\frac{\partial }{\partial \alpha }\right)u^k(\tilde{\alpha },z)+({D_2}^k+{D_3}^k\frac{\partial }{\partial z})w^k(\tilde{\alpha },z)=0}\\ {\tau }_{\alpha z}^k(\tilde{\alpha },z)=({G_1}^k+{G_2}^k\frac{\partial }{\partial z})u^k(\tilde{\alpha },z)+\left({G_3}^k\frac{\partial }{\partial \alpha }\right)w^k(\tilde{\alpha },z)=0\end{array}$$

(12)

Simply supported (S):

$$\begin{array}{c}{\sigma }_{\alpha }^k(\tilde{\alpha },z)=\left({D_1}^k\frac{\partial }{\partial \alpha }\right)u^k(\tilde{\alpha },z)\hfill \\ +({D_2}^k+{D_3}^k\frac{\partial }{\partial z})w^k(\tilde{\alpha },z)=0w^k(\tilde{\alpha },z)=0\hfill \end{array}$$

(13)

Clamped (C):

$$\begin{array}{c}{u}^k(\tilde{\alpha },z)=0\hfill \\ w^k(\tilde{\alpha },z)=0\hfill \end{array}$$

(14)

The boundary conditions on the lower ( $z=z_0$ ) and upper ( $z=z_{\mathrm{NL}}$ ) surfaces for transversal loads are defined as:

$$\begin{array}{c}{\sigma }_z^1(\alpha ,z_0)=\left({F_1}^1\frac{\partial }{\partial \alpha }\right)u^1(\alpha ,z_0)+({F_2}^1+{F_3}^1\frac{\partial }{\partial z})w^1(\alpha ,z_0)=P_0\left(\alpha \right)\hfill \\ {\tau }_{\alpha z}^1(\alpha ,z_0)=({G_1}^1+{G_2}^1\frac{\partial }{\partial z})u^1(\alpha ,z_0)+\left({G_3}^1\frac{\partial }{\partial \alpha }\right)w^1(\alpha ,z_0)=0\hfill \end{array}$$

(15)

$$\begin{array}{c}{\sigma }_z^{\mathrm{NL}}(\alpha ,z_{\mathrm{NL}})=\left({F_1}^{\mathrm{NL}}\frac{\partial }{\partial \alpha }\right)u^{\mathrm{NL}}(\alpha ,z_{\mathrm{NL}})\hfill \\ +({F_2}^{\mathrm{NL}}+{F_3}^{\mathrm{NL}}\frac{\partial }{\partial z})w^{\mathrm{NL}}(\alpha ,z_{\mathrm{NL}})=P_{\mathrm{NL}}\left(\alpha \right)\hfill \\ {\tau }_{\alpha z}^{\mathrm{NL}}(\alpha ,z_{\mathrm{NL}})=({G_1}^{\mathrm{NL}}+{G_2}^{\mathrm{NL}}\frac{\partial }{\partial z})u^{\mathrm{NL}}(\alpha ,z_{\mathrm{NL}})\hfill \\ +\left({G_3}^{\mathrm{NL}}\frac{\partial }{\partial \alpha }\right)w^{\mathrm{NL}}(\alpha ,z_{\mathrm{NL}})=0\hfill \end{array}$$

(16)

where $P_0$ and $P_{\mathrm{NL}}$ represent the transversal load applied on the lower and upper surfaces of the beam, respectively.

In the case of multilayer beams, continuity of displacements and stresses in the $z$ -direction must be satisfied. The interlaminar continuity equations for displacements are mathematically defined as:

$$\begin{array}{c}{u}^k(\alpha ,z_k)=u^{k+1}(\alpha ,z_k)\hfill \\ w^k(\alpha ,z_k)=w^{k+1}(\alpha ,z_k)\hfill \end{array}$$

(17)

Whereas the interlaminar continuity equations for stresses are defined as:

$$\begin{array}{c}{\sigma }_z^k(\alpha ,z_k)={\sigma }_z^{k+1}(\alpha ,z_k)\hfill \\ {\tau }_{\alpha z}^k(\alpha ,z_k)={\tau }_{\alpha z}^{k+1}(\alpha ,z_k)\hfill \end{array}$$

(18)

Substituting equation (6) into equation (18), the interlaminar stress continuity equations in the $z$ -direction can be expressed in terms of displacements as follows:

$$\begin{array}{c}\left({F_1}^k\frac{\partial }{\partial \alpha }\right)u^k(\alpha ,z_k)+({F_2}^k+{F_3}^k\frac{\partial }{\partial z})w^k(\alpha ,z_k)\hfill \\ =\left({F_1}^{k+1}\frac{\partial }{\partial \alpha }\right)u^{k+1}(\alpha ,z_k)+({F_2}^{k+1}+{F_3}^{k+1}\frac{\partial }{\partial z})w^{k+1}(\alpha ,z_k)\hfill \\ ({G_1}^k+{G_2}^k\frac{\partial }{\partial z})u^k(\alpha ,z_k)+\left({G_3}^k\frac{\partial }{\partial \alpha }\right)w^k(\alpha ,z_k)\hfill \\ =({G_1}^{k+1}+{G_2}^{k+1}\frac{\partial }{\partial z})u^{k+1}(\alpha ,z_k)+\left({G_3}^{k+1}\frac{\partial }{\partial \alpha }\right)w^{k+1}(\alpha ,z_k)\hfill \end{array}$$

(19)

Numerical solution using the Differential Quadrature Method

The governing equations (see equation (10)), together with the boundary conditions presented in Section Boundary conditions constitute a system of partial differential equations in which the displacements $u^k(\alpha ,z)$ and $w^k(\alpha ,z)$ are the unknown functions, and $\alpha$ and $z$ are the independent variables. This system can be fully transformed into a linear system of algebraic equations by means of the DQM. This method involves discretizing the problem domain using a finite set of points, known as a grid, and approximating the derivatives of the unknown function at these points by means of a weighted linear combination of the function values at all the grid points. ⁴⁵

DQM in a two-dimensional domain

Let $f(\alpha ,z)$ be a continuous function defined over the domain $[{\alpha }_a,{\alpha }_b]\times [z_a,z_b]$ . The domain is discretized by dividing it into a set of points $({\propto }_i,{\zeta }_j)$ , with $\{{\propto }_i|i=1,2,\dots ,N\}$ in the $\alpha$ -direction and $\{{\zeta }_j|j=1,2,\dots ,M\}$ in the $z$ -direction, where $N$ and $M$ are the total number of divisions along the $\alpha$ and $z$ directions, respectively. The total number of nodes in the grid is thus given by $N\times M$ .

The function $f(\alpha ,z)$ and its partial derivatives with respect to $\alpha$ and $z$ , evaluated at a point $({\propto }_i,{\zeta }_j)$ , are approximated as follows:

$$\begin{array}{c}f({\propto }_i,{\zeta }_j)=\sum _{n=1}^N\sum _{m=1}^M{\lambda }_{\mathrm{in}}^{\left(\alpha 0\right)}{\lambda }_{\mathrm{jm}}^{\left(z0\right)}f({\propto }_n,{\zeta }_m)\hfill \\ \frac{{\partial }^{n_t}f({\propto }_i,{\zeta }_j)}{\partial {\alpha }^{n_t}}=\sum _{n=1}^N\sum _{m=1}^M{\lambda }_{\mathrm{in}}^{\left(\alpha n_t\right)}{\lambda }_{\mathrm{jm}}^{\left(z0\right)}f({\propto }_n,{\zeta }_m)\hfill \\ \frac{{\partial }^{n_t}f({\propto }_i,{\zeta }_j)}{\partial z^{n_t}}=\sum _{n=1}^N\sum _{m=1}^M{\lambda }_{\mathrm{in}}^{\left(\alpha 0\right)}{\lambda }_{\mathrm{jm}}^{\left(z{n}_t\right)}f({\propto }_n,{\zeta }_m)\hfill \\ \frac{{\partial }^{n_t}f({\propto }_i,{\zeta }_j)}{\partial {\alpha }^{n_1}\partial z^{n_2}}=\frac{{\partial }^{n_t}f({\propto }_i,{\zeta }_j)}{\partial z^{n_2}\partial {\alpha }^{n_1}}\hfill \\ =\sum _{n=1}^N\sum _{m=1}^M{\lambda }_{\mathrm{in}}^{\left(\alpha n_1\right)}{\lambda }_{\mathrm{jm}}^{\left(z{n}_2\right)}f({\propto }_n,{\zeta }_m);n_1+n_2=n_t\hfill \end{array}$$

(20)

where ${\lambda }_{\mathrm{in}}^{\left(\alpha n_t\right)}$ and ${\lambda }_{\mathrm{jm}}^{\left(z{n}_t\right)}$ are the weighting coefficients associated with the discretization of the domain in the $\alpha$ - and $z$ -directions, respectively; and ${\lambda }_{\mathrm{in}}^{\left(\alpha 0\right)}=0$ for $n\ne i$ , ${\lambda }_{\mathrm{in}}^{\left(\alpha 0\right)}=1$ for $n=i$ , ${\lambda }_{\mathrm{jm}}^{\left(z0\right)}=0$ for $m\ne j$ and ${\lambda }_{\mathrm{jm}}^{\left(z0\right)}=1$ for $m=j$ .

Domain discretization

In the present study, the Chebyshev-Gauss-Lobatto distribution was employed to discretize the problem domain. The coordinates of the grid distribution along the $z$ -axis are obtained using the following expression:

$${\zeta }_j={\zeta }_1+\frac{{\zeta }_M-{\zeta }_1}{2}(1-\cos \left(\frac{j-1}{M-1}\pi \right))j=1,2,\dots ,M;{\zeta }_1=z_a;{\zeta }_M=z_b$$

(21)

Similarly, the grid coordinates along the $\alpha$ -axis are obtained.

Weighting coefficients

The selection of the basis functions is a crucial step in the numerical implementation of the DQM. In this work, Lagrange interpolation polynomials are used as the set of basis functions. The use of Lagrange polynomials allows for the computation of the weighting coefficients for any derivative of order $n_t$ and grid distribution. The weighting coefficients associated with the discretization of the domain in the $z$ -direction are obtained by:

$$\begin{array}{c}{\lambda }_{\mathrm{jm}}^{\left(n_t\right)}=n_t({\lambda }_{\mathrm{jm}}^{\left(1\right)}{\lambda }_{\mathrm{jj}}^{(n_t-1)}-\frac{{\lambda }_{\mathrm{jm}}^{(n_t-1)}}{{\zeta }_j-{\zeta }_m})\hfill \\ \mathrm{para}j\ne m,j,m=1,2,\dots ,M,n_t=2,3,\dots ,M-1\hfill \\ {\lambda }_{\mathrm{jj}}^{\left(n_t\right)}=-\sum _{m=1,m\ne j}^M{\lambda }_{\mathrm{jm}}^{\left(n_t\right)}\hfill \\ \mathrm{para}j=1,2,\dots ,M,n_t=1,2,\dots ,M-1\hfill \\ {\lambda }_{\mathrm{jm}}^{\left(1\right)}=\frac{L^{\left(1\right)}\left({\zeta }_j\right)}{({\zeta }_j-{\zeta }_m)L^{\left(1\right)}\left({\zeta }_m\right)}\hfill \\ \mathrm{para}j\ne m,j,m=1,2,\dots ,M\hfill \\ L^{\left(1\right)}({\zeta }_j)=\underset{m=1,m\ne j}{\overset{M}{\Pi }}({\zeta }_j-{\zeta }_m)\hfill \end{array}$$

(22)

The weighting coefficients related to the discretization of the domain in the $\alpha$ -direction are obtained in a similar manner.

Finally, the linear system of algebraic equations for curved beams with arbitrary boundary conditions is obtained by applying equations (20), (21), and (22) to the governing equations (equation (10)), the boundary conditions at $\alpha =0$ and $\alpha =L$ (equations (12), (13), or (14)), the boundary conditions at $z=z_0$ and $z=z_{\mathrm{NL}}$ (equations (15) and (16)), and the interlaminar continuity equations (equations (17) and (19)). The resulting algebraic equations are presented in Appendix A.

Numerical results and discussion

In this section, the static analysis of functionally graded porous curved beams with simply supported (S-S), clamped-clamped (C-C), and clamped-free (C-F) boundary conditions is presented. The analysis was conducted for deep curved beams considering various values of the power-law index $p$ and porosity volume fraction $a$ .

The following assumptions were considered for all cases:

• The FGM consists of aluminum (Al: $E_m=70$ GPa, ${\nu }_m=0.3$ ) as the metallic phase and alumina (Al₂O₃: $E_c=380$ GPa, ${\nu }_c=0.3$ ) as the ceramic phase.

• In the DQM solutions, the domain of each layer was discretized along the $\alpha$ and $z$ axes using a 35 × 35 grid for all cases (S-S, C-C, and C-F).

• The two-dimensional finite element method (2D FEM) was employed to obtain reference results via the commercial software ANSYS. The numerical model was built using PLANE182 elements. The beam thickness was discretized into 250 layers for $p=2$ and 110 layers for $p=5$ , such that the layer thickness ensured a constant variation of Young’s modulus between adjacent layers.

• The dimensionless results presented use the following normalization formulas:

$$\begin{array}{c}\overline{u}=\frac{100E_m{h}^3}{L^4\left|P_0\right|}u;\overline{w}=\frac{100E_m{h}^3}{L^4\left|P_0\right|}w;\hfill \end{array}$$

$${\overline{\sigma }}_{\alpha }=\frac{h}{L\left|P_0\right|}{\sigma }_{\alpha };{\overline{\sigma }}_z=\frac{h}{L\left|P_0\right|}{\sigma }_z;{\overline{\tau }}_{\alpha z}=\frac{h}{L\left|P_0\right|}{\tau }_{\alpha z}$$

Convergence and computational efficiency

The convergence of the proposed model was evaluated in curved beams with symmetric (C-C) and asymmetric (C-F) support conditions, considering deep curvature ( $L/R=1$ ) and different material gradations ( $p=5,7,10$ ). Grid configurations of $N$ × $M$ were analyzed with an equal number of divisions in both directions ( $N=M$ ) and with higher refinement in the axial direction ( $N>M$ ).

The results presented in Table 1 confirm the rapid convergence of the model. For all support conditions and values of $p$ , the solutions stabilize from grids of approximately 35 × 35 ( $N=M$ ) and 75 × 25 ( $N>M$ ) discretization points, with differences smaller than ${10}^{-3}$ observed in subsequent refinements. This behavior reflects the high efficiency of the method, capable of achieving stable solutions with moderate discretizations and avoiding the use of excessively dense grids that would increase the computational cost.

Table 1.

Convergence of the DQM model for functionally graded porous curved beams ( $L/R=1$ , $L/h=5$ , $a=0.2$ ): axial stress ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ .

Grid	Number of nodes	Symmetric beam (C-C)			Asymmetric beam (C-F)
		$p=5$	$p=7$	$p=10$	$p=5$	$p=7$	$p=10$
11 × 11	121	−4.3264	−4.9226	−5.7286	10.265	11.285	12.485
19 × 19	361	−4.3364	−4.9116	−5.5928	10.503	11.263	12.578
27 × 27	729	−4.3353	−4.9091	−5.5915	10.531	11.295	12.619
35 × 35	1225	−4.3351	−4.9087	−5.5906	10.537	11.303	12.631
43 × 43	1849	−4.3350	−4.9085	−5.5903	10.539	11.306	12.634
51 × 51	2601	−4.3350	−4.9085	−5.5902	10.539	11.307	12.636
35 × 25	875	−4.3373	−4.9113	−5.5936	10.535	11.301	12.627
45 × 25	1125	−4.3328	−4.9057	−5.5867	10.555	11.325	12.657
55 × 25	1375	−4.3366	−4.9102	−5.5921	10.535	11.302	12.631
65 × 25	1625	−4.3344	−4.9076	−5.5890	10.544	11.312	12.642
75 × 25	1875	−4.3354	−4.9088	−5.5904	10.538	11.306	12.635
85 × 25	2125	−4.3350	−4.9084	−5.5899	10.540	11.308	12.637

The computational efficiency of the model was evaluated through comparison with a finite element method (FEM) analysis in commercial software. In the case of FEM, representing the functional gradation of the material requires discretizing the thickness into numerous layers, which considerably increases the number of elements and, consequently, the computational cost. In contrast, the DQM directly captures the continuous variation of properties without the need for additional subdivisions. In the tests performed, the MATLAB implementation of the DQM model, using a 35 × 35 grid, achieved the same accuracy as a 2D FEM model in ANSYS APDL with a 250-layer discretization and approximately 1,000,000 elements, but with a computation time about 93% lower. The computation times of both models are presented in Table 2.

Table 2.

Computation time comparison between DQM and FEM (ANSYS APDL) for functionally graded porous curved beams.

Method	Average computation time (s)
2D FEM	197
Present (DQM)	13.2

The remarkable efficiency of the DQM lies in its direct solution of the original differential equations of the problem at each discretization point, using high-order basis functions. This formulation allows for a continuous representation of the spatial variation of the mechanical properties of functionally graded materials, avoiding the artificial interfaces associated with the layer-wise discretization of conventional FEM methods. Moreover, the use of optimal distributions of discretization points promotes rapid convergence with a reduced number of degrees of freedom.

Nevertheless, the method presents certain limitations that must be considered. The analysis of highly irregular geometries or complex boundary conditions may require further developments, and numerical stability largely depends on an appropriate selection of discretization points. In addition, due to the inherent formulation of the DQM, the resulting stiffness matrix is not positive definite, which necessitates the implementation of specific techniques to ensure the robustness of the computation.

Despite these considerations, the DQM represents a highly efficient alternative to conventional FEM-based methods for the static analysis of functionally graded curved beams, offering a favorable balance between numerical stability and computational cost.

Validation

A static analysis was performed on functionally graded porous curved beams subjected to a uniformly distributed load applied on their upper surface to validate the accuracy of the proposed DQM-based model. The analysis considered various boundary conditions as well as different values of the power-law index $p$ and porosity volume fraction $a$ .

The results obtained using the proposed model were compared with those of the 2D FEM numerical model. Additionally, for S-S beams, the results were also compared with those reported by Sayyad and Ghugal, ¹⁶ Belarbi et al., ¹⁸ and Sayyad et al., ²⁸ which are based on higher-order shear deformation theories (HSDTs). Table 3 presents the displacement and stress results for shallow curved beams with S-S supports, a length-to-radius ratio of $L/R=0.2$ , and a length-to-thickness ratio of $L/h=5$ . The results obtained for different values of $p$ and $a$ show high accuracy when compared to those from the 2D FEM model and exhibit improved precision over the models presented in refs.,^16,18,28 both in terms of displacement and stress evaluation.

Table 3.

Dimensionless results for $\overline{u}(0,-h/2)$ , $\overline{w}(L/2,0)$ , ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ , and ${\overline{\tau }}_{\alpha z}(0,0)$ for S-S beams ( $L/h=5$ , $L/R=0.2$ ). ^a

$p$	$a$	Model	$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$
2	0	Sayyad and Ghugal 16	5.1816	7.8773	6.9448	0.6832
		Belarbi et al. 18	-	8.0513	6.9000	0.6902
		Sayyad et al. 28	5.3034	8.0675	6.8835	0.6661
		2D FEM	3.6830	8.3619	6.9423	0.6893
		Present	3.6833	8.3624	6.9511	0.6893
	0.1	Sayyad et al. 28	6.7135	10.008	7.5691	0.6552
		2D FEM	4.7213	10.384	7.6427	0.6798
		Present	4.7218	10.385	7.6535	0.6798
	0.2	Sayyad et al. 28	4.4048	13.541	8.8032	0.6390
		2D FEM	6.6845	14.077	8.9014	0.6628
		Presente	6.6856	14.078	8.9154	0.6627
5	0	Sayyad and Ghugal 16	6.2523	9.6298	8.1850	0.6074
		Belarbi et al. 18	-	9.8153	8.1387	0.6107
		Sayyad et al. 28	6.3831	9.8268	8.1116	0.5862
		2D FEM	4.4133	10.173	8.1430	0.6586
		Present	4.4228	10.189	8.1792	0.6591
	0.1	Sayyad et al. 28	8.5075	12.916	9.0725	0.5457
		2D FEM	6.0553	13.397	9.1313	0.6412
		Present	6.0736	13.426	9.1791	0.6419
	0.2	Sayyad et al. 28	14.033	20.089	11.062	0.4698
		2D FEM	10.089	20.955	11.185	0.6049
		Present	10.141	21.035	11.263	0.6061

a

All results are in absolute values.

Tables 4 to 6 present displacement and stress results for deep curved beams ( $L/R=1$ and $L/h=5$ ) with S-S, C-C, and C-F boundary conditions, respectively. In all cases, excellent agreement is observed between the proposed model and the 2D FEM results. These findings demonstrate the capability of the DQM-based model to yield accurate solutions in the analysis of functionally graded porous curved beams under various boundary conditions.

Table 4.

Dimensionless results for $\overline{u}(0,-h/2)$ , $\overline{w}(L/2,0)$ , ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ , and ${\overline{\tau }}_{\alpha z}(0,0)$ for S-S beams ( $L/h=5$ , $L/R=1$ ).

$p$	$a$	Model	$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$
2	0	2D FEM	−7.4118	−11.413	−7.8060	−0.8268
		Present	−7.4124	−11.414	−7.8167	−0.8268
	0.1	2D FEM	−9.4420	−14.281	−8.6266	−0.8189
		Present	−9.4430	−14.282	−8.6396	−0.8189
	0.2	2D FEM	−13.281	−19.599	−10.119	−0.8045
		Present	−13.283	−19.601	−10.136	−0.8044
5	0	2D FEM	−8.9047	−13.836	−9.1445	−0.7876
		Present	−8.9235	−13.861	−9.1880	−0.7883
	0.1	2D FEM	−12.085	−18.355	−10.277	−0.7696
		Present	−12.121	−18.401	−10.335	−0.7705
	0.2	2D FEM	−19.831	−29.079	−12.642	−0.7312
		Present	−19.934	−29.207	−12.738	−0.7328

Table 5.

Dimensionless results for $\overline{u}(L/4,-h/2)$ , $\overline{w}(L/2,0)$ , ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ , and ${\overline{\tau }}_{\alpha z}(L/4,0)$ for C-C beams ( $L/h=5$ , $L/R=1$ ).

$p$	$a$	Model	$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$
2	0	2D FEM	−0.4550	−1.3918	−2.3810	−0.2234
		Present	−0.4551	−1.3918	−2.3847	−0.2234
	0.1	2D FEM	−0.5682	−1.6696	−2.6064	−0.2168
		Present	−0.5683	−1.6696	−2.6103	−0.2168
	0.2	2D FEM	−0.7630	−2.1133	−2.9648	−0.2034
		Present	−0.7632	−2.1134	−2.9694	−0.2034
5	0	2D FEM	−0.5781	−1.8734	−3.0352	−0.2208
		Present	−0.5791	−1.8744	−3.0481	−0.2210
	0.1	2D FEM	−0.7711	−2.4109	−3.4802	−0.2118
		Present	−0.7731	−2.4127	−3.4962	−0.2120
	0.2	2D FEM	−1.1805	−3.4792	−4.3120	−0.1890
		Present	−1.1857	−3.4831	−4.3351	−0.1893

Table 6.

Dimensionless results for $\overline{u}(L,h/2)$ , $\overline{w}(L,0)$ , ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ , and ${\overline{\tau }}_{\alpha z}(L/4,0)$ for C-F beams ( $L/h=5$ , $L/R=1$ ).

$p$	$a$	Model	$\overline{u}$ .	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$
2	0	2D FEM	37.568	−70.677	6.6102	−1.0848
		Present	37.566	−70.676	6.6194	−1.0849
	0.1	2D FEM	46.536	−88.835	7.3028	−1.0753
		Present	46.532	−88.830	7.3151	−1.0754
	0.2	2D FEM	63.055	−122.83	5847	−1.0573
		Present	63.050	−122.82	8.6014	−1.0573
5	0	2D FEM	44.866	−84.045	7.7132	−1.0380
		Present	44.958	−84.242	7.7559	−1.0391
	0.1	2D FEM	58.330	−111.53	8.5965	−1.0179
		Present	58.503	−111.91	8.6558	−1.0194
	0.2	2D FEM	89.632	−177.27	10.433	−0.9733
		Present	90.121	−178.33	10.537	−0.9758

Figure 3 shows through-the-thickness results, including the deflection $\overline{w}(L/2,z)$ , the normal stresses ${\overline{\sigma }}_{\alpha }(L/2,z)$ and ${\overline{\sigma }}_z(L/2,z)$ , and the shear stress ${\overline{\tau }}_{\alpha z}(L/4,z)$ . The DQM model exhibits an accurate distribution of results in comparison with those obtained from the 2D FEM model for all three types of support conditions analyzed. In particular, the stress distributions satisfy the expected physical conditions: the axial normal stress ( ${\overline{\sigma }}_{\alpha }$ ) is predominantly concentrated on the upper surface of the beam; the transverse normal stress ( ${\overline{\sigma }}_z$ ) is zero at the lower surface and equal to the applied load at the upper surface; and the shear stress ( ${\overline{\tau }}_{\alpha z}$ ) vanishes at both the upper and lower surfaces, in accordance with theoretical expectations. These results confirm the validity and accuracy of the DQM-based model for the static analysis of functionally graded porous deep curved beams with various power-law index $p$ and porosity volume fraction $a$ .

Figure 3.

Comparison of dimensionless deflection and stress distributions through the thickness of the beam between the DQM and 2D FEM models ( $L/h=5$ , $L/R=1$ , $p=5$ , $a=0$ ): (a) $\overline{w}(L/2,z)$ , (b) ${\overline{\sigma }}_{\alpha }(L/2,z)$ , (c) ${\overline{\sigma }}_z(L/2,z)$ , and (d) ${\overline{\tau }}_{\alpha z}(L/4,z)$ .

Analysis of curved beams with S-S, C-C, and C-F support conditions

This section presents the analysis of functionally graded porous curved beams under different boundary conditions: S-S, C-C, and C-F. In all cases, a uniformly distributed load is applied to the upper surface of the beam. The study evaluates the influence of curvature ( $L/R$ ), the power-law index $p$ , and the porosity volume fraction $a$ on the structural behavior. Tables 7 to 9 present the results obtained using the DQM for each type of support.

Table 7.

Dimensionless results for $\overline{u}(0,-h/2)$ , $\overline{w}(L/2,0)$ , ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ , and ${\overline{\tau }}_{\alpha z}(0,0)$ for S-S beams ( $L/h=5$ ).

$p$	$a$	$L/R=0.5$				$L/R=1$
		$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$	$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$
1	0	−3.5823	−7.0154	−6.0864	−0.7713	−5.6011	−8.7808	−6.6555	−0.8744
	0.1	−4.2269	−8.1424	−6.4447	−0.7724	−6.5988	−10.224	−7.0600	−0.8772
	0.2	−5.1807	−9.7757	−6.9648	−0.7741	−8.0752	−12.328	−7.6501	−0.8815
2	0	−4.7693	−9.0950	−7.1437	−0.7283	−7.4124	−11.414	−7.8167	−0.8268
	0.1	−6.0924	−11.327	−7.8769	−0.7193	−9.4430	−14.282	−8.6396	−0.8189
	0.2	−8.5915	−15.427	−9.2002	−0.7033	−13.283	−19.601	−10.136	−0.8044
3	0	−5.3108	−10.115	−7.6756	−0.7060	−8.2444	−12.688	−8.3946	−0.8010
	0.1	−7.0591	−13.052	−8.5886	−0.6902	−10.916	−16.453	−9.4126	−0.7853
	0.2	−10.874	−19.260	−10.430	−0.6583	−16.745	−24.485	−11.480	−0.7529
5	0	−5.7349	−11.067	−8.4019	−0.6956	−8.9235	−13.861	−9.1880	−0.7883
	0.1	−7.8303	−14.624	−9.4371	−0.6784	−12.121	−18.401	−10.335	−0.7705
	0.2	−12.970	−23.025	−11.599	−0.6422	−19.934	−29.207	−12.738	−0.7328
7	0	−5.9320	−11.618	−9.0609	−0.6957	−9.2617	−14.534	−9.9121	−0.7876
	0.1	−8.1483	−15.439	−10.187	−0.6800	−12.651	−19.399	−11.158	−0.7716
	0.2	−13.766	−24.715	−12.426	−0.6496	−21.188	−31.293	−13.640	−0.7399
10	0	−6.1389	−12.248	−9.9870	−0.6976	−9.6269	−15.301	−10.929	−0.7892
	0.1	−8.4467	−16.319	−11.297	−0.6831	−13.172	−20.473	−12.380	−0.7743
	0.2	−14.347	−26.225	−13.712	−0.6577	−22.160	−33.141	−15.057	−0.7480

Table 8.

Dimensionless results for $\overline{u}(L/4,-h/2)$ , $\overline{w}(L/2,0)$ , ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ , and ${\overline{\tau }}_{\alpha z}(L/4,0)$ for C-C beams ( $L/h=5$ ).

$p$	$a$	$L/R=0.5$				$L/R=1$
		$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$	$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$
1	0	−0.4310	−1.4561	−2.2932	−0.3328	−0.3441	−1.0570	−1.9428	−0.2366
	0.1	−0.5088	−1.6621	−2.4197	−0.3326	−0.4034	−1.2027	−2.0461	−0.2355
	0.2	−0.6226	−1.9492	−2.6003	−0.3317	−0.4877	−1.4004	−2.1872	−0.2328
2	0	−0.5721	−1.9168	−2.7767	−0.3150	−0.4551	−1.3918	−2.3847	−0.2234
	0.1	−0.7257	−2.3218	−3.0557	−0.3092	−0.5683	−1.6696	−2.6103	−0.2168
	0.2	−1.0060	−3.0073	−3.5343	−0.2977	−0.7632	−2.1134	−2.9694	−0.2034
3	0	−0.6427	−2.2033	−3.0614	−0.3074	−0.5156	−1.6097	−2.6687	−0.2193
	0.1	−0.8461	−2.7648	−3.4425	−0.2987	−0.6664	−1.9967	−2.9881	−0.2101
	0.2	−1.2662	−3.8373	−4.1671	−0.2793	−0.9555	−2.6834	−3.5428	−0.1894
5	0	−0.7088	−2.5351	−3.4282	−0.3058	−0.5791	−1.8744	−3.0481	−0.2210
	0.1	−0.9620	−3.2958	−3.9246	−0.2972	−0.7731	−2.4127	−3.4962	−0.2120
	0.2	−1.5462	−4.9312	−4.9392	−0.2764	−1.1857	−3.4831	−4.3351	−0.1893
7	0	−0.7429	−2.7299	−3.7026	−0.3068	−0.6146	−2.0342	−3.3218	−0.2238
	0.1	−1.0199	−3.6089	−4.2751	−0.2995	−0.8331	−2.6708	−3.8644	−0.2164
	0.2	−1.6871	−5.5987	−5.4557	−0.2823	−1.3207	−4.0123	−4.9087	−0.1965
10	0	−0.7755	−2.9194	−4.0410	−0.3079	−0.6487	−2.1895	−3.6395	−0.2266
	0.1	−1.0729	−3.9111	−4.7116	−0.3017	−0.8906	−2.9253	−4.2986	−0.2209
	0.2	−1.8078	−6.2365	−6.0787	−0.2878	−1.4507	−4.5508	−5.5906	−0.2046

Table 9.

Dimensionless results for $\overline{u}(L,h/2)$ , $\overline{w}(L,0)$ , ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ , and ${\overline{\tau }}_{\alpha z}(L/4,0)$ for C-F beams ( $L/h=5$ ).

$p$	$a$	$L/R=0.5$				$L/R=1$
		$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$	$\overline{u}$	$\overline{w}$	${\overline{\sigma }}_{\alpha }$	${\overline{\tau }}_{\alpha z}$
1	0	18.124	−58.579	5.7335	−1.1649	29.633	−54.741	5.6850	−1.1430
	0.1	20.840	−68.183	6.0756	−1.1668	34.316	−63.936	6.0350	−1.1471
	0.2	24.729	−82.159	6.5751	−1.1698	41.103	−77.403	6.5486	−1.1531
2	0	22.650	−75.454	6.6687	−1.1045	37.566	−70.676	6.6194	−1.0849
	0.1	27.722	−94.358	7.3525	−1.0918	46.532	−88.830	7.3151	−1.0754
	0.2	36.901	−129.41	8.6086	−1.0684	63.050	−122.82	8.6014	−1.0573
3	0	24.827	−83.266	7.1236	−1.0744	41.309	−77.925	7.0720	−1.0544
	0.1	31.283	−107.82	7.9457	−1.0525	52.805	−101.43	7.9041	−1.0357
	0.2	44.658	−160.33	9.6386	−1.0065	77.118	−152.18	9.6258	−0.9955
5	0	27.080	−90.186	7.8065	−1.0603	44.958	−84.242	7.7559	−1.0391
	0.1	34.666	−119.23	8.6968	−1.0380	58.503	−111.91	8.6558	−1.0194
	0.2	51.947	−188.43	10.555	−0.9895	90.121	−178.33	10.537	−0.9758
7	0	28.675	−94.318	8.4842	−1.0582	47.380	−87.992	8.4332	−1.0362
	0.1	36.796	−125.08	9.4485	−1.0386	61.823	−117.22	9.4096	−1.0190
	0.2	55.503	−200.02	11.317	−1.0004	96.079	−188.90	11.303	−0.9848
10	0	30.732	−99.374	9.4590	−1.0572	50.418	−92.580	9.4019	−1.0347
	0.1	39.509	−131.86	10.610	−1.0388	65.882	−123.37	10.568	−1.0184
	0.2	59.425	−210.62	12.639	−1.0077	102.19	−198.51	12.631	−0.9907

In S-S beams, the displacements $\overline{u}(0,-h/2)$ and $\overline{w}(L/2,0)$ increase with $L/R$ , $p$ , and $a$ . This behavior is attributed to the fact that greater curvature reduces geometric stiffness, while high values of $p$ and $a$ decrease the elastic modulus and the effective load-carrying area, leading to larger deformations. The axial stress ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ also increases with $p$ and $a$ due to stress concentration in regions of lower stiffness, whereas the shear stress ${\overline{\tau }}_{\alpha z}(0,0)$ increases primarily with curvature because of shear flow redistribution in more highly arched sections.

In C-C beams, the displacements $\overline{u}(L/4,-h/2)$ and $\overline{w}(L/2,0)$ decrease with curvature because the end constraints partially compensate for the loss of stiffness. However, these displacements increase with $p$ and $a$ as a result of the degradation of mechanical properties. The axial stress ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ also exhibits an increasing trend with $p$ and $a$ , whereas the shear stress ${\overline{\tau }}_{\alpha z}(L/4,0)$ decreases with greater curvature due to the structural constraint imposed at the ends, which limits transverse deformation in that region.

In C-F beams, the displacement $\overline{u}(L,h/2)$ increases significantly with $L/R$ , $p$ , and $a$ , due to the asymmetry in the support conditions allowing greater free deformation at the unconstrained end. The deflection $\overline{w}(L,0)$ also increases with $p$ and $a$ , but shows a slight reduction at very high curvatures because of a combined effect of geometric stiffness loss and partial constraints. The axial stress ${\overline{\sigma }}_{\alpha }(L/2,h/2)$ maintains the increasing trend with $p$ and $a$ , whereas the shear stress ${\overline{\tau }}_{\alpha z}(L/4,0)$ exhibits minimal variations with respect to the analyzed parameters.

The trends described are illustrated in Figures 4 to 6, which show the variation of displacements and stresses as a function of curvature ( $L/R$ ), the power-law index $p$ , and porosity fraction $a$ for the three boundary conditions analyzed.

Figure 4.

Distribution of dimensionless deflection and stresses through the thickness of the beam for various values of $L/R$ ( $L/h=5$ , $p=5$ , $a=0$ ): (a) $\overline{w}(L/2,z)$ , (b) ${\overline{\sigma }}_{\alpha }(L/2,z)$ , (c) ${\overline{\sigma }}_z(L/2,z)$ , and (d) ${\overline{\tau }}_{\alpha z}(L/4,z)$ .

Figure 5.

Distribution of dimensionless deflection and stresses through the thickness of the beam for various values of $p$ ( $L/h=5$ , $L/R=1$ , $a=0$ ): (a) $\overline{w}(L/2,z)$ , (b) ${\overline{\sigma }}_{\alpha }(L/2,z)$ , (c) ${\overline{\sigma }}_z(L/2,z)$ , and (d) ${\overline{\tau }}_{\alpha z}(L/4,z)$ .

Figure 6.

Distribution of dimensionless deflection and stresses through the thickness of the beam for various values of $a$ ( $L/h=5$ , $L/R=1$ , $p=5$ ): (a) $\overline{w}(L/2,z)$ , (b) ${\overline{\sigma }}_{\alpha }(L/2,z)$ , (c) ${\overline{\sigma }}_z(L/2,z)$ , and (d) ${\overline{\tau }}_{\alpha z}(L/4,z)$ .

Figure 4 presents the through-thickness distributions of displacements and stresses for different curvatures, where the case $L/R=0$ corresponds to a straight beam. In Figure 4(a), it is observed that, as curvature increases, the deflection $\overline{w}(L/2,z)$ increases in S-S beams, decreases in C-C beams, and exhibits a non-monotonic behavior in C-F beams, with an initial increase at $L/R=0.5$ followed by a reduction at $L/R=1$ . These variations are attributed to the effect of curvature on geometric stiffness: in S-S beams, the reduction in bending stiffness results in larger deformations, whereas in C-C and C-F beams the end constraints partially mitigate this loss of stiffness. The stress distribution behavior shown in Figures 4b, 4c, and 4d is primarily governed by the combined effects of axial and transverse bending induced by curvature, as well as the redistribution of shear stresses. In S-S beams, rotational freedom allows for greater deformation, which in turn leads to stress increases in certain regions. In C-C beams, the constraints at both ends limit deflection and redistribute internal loads, promoting axial compression and reducing shear. In C-F beams, the asymmetry in support conditions results in a non-uniform redistribution of stresses and deformations along the arc. In Figure 4(b), the distribution of axial stress ${\overline{\sigma }}_{\alpha }(L/2,z)$ reveals that, in S-S and C-F beams, the curves intersect at a common point associated with zero axial stress. In S-S beams, axial stresses increase with curvature, primarily in tension, whereas in C-F beams, axial stresses remain nearly constant in tension but slightly intensify in compression. In C-C beams, the stress curve shifts to the left with increasing curvature, indicating lower tensile stresses and higher compressive stresses. Figure 4(c) shows that, for straight beams, the three configurations share identical distributions of transverse stress ${\overline{\sigma }}_z(L/2,z)$ , which are entirely compressive. However, curvature introduces differences: in S-S beams, tensile stresses appear, increasing primarily in the lower region; in C-F beams, compression intensifies in the upper region, whereas in C-C beams, the variation with curvature is minimal. Finally, Figure 4d shows that, in straight beams, the shear stress distributions ${\overline{\tau }}_{\alpha z}(L/4,z)$ are identical for S-S and C-C beams. With increasing curvature, shear stress increases in S-S beams and decreases in C-C beams, whereas in C-F beams it remains practically constant.

Figure 5 shows the through-thickness distributions of displacements and stresses for different values of the power-law index $p$ . In Figure 5(a), the deflection $\overline{w}(L/2,z)$ increases with $p$ for all support types, reflecting the loss of stiffness associated with steeper gradients in the elastic modulus. In Figure 5(b), the axial stress ${\overline{\sigma }}_{\alpha }(L/2,z)$ progressively tends to concentrate near the upper and lower regions as $p$ increases, particularly in S-S and C-F beams. For high values of $p$ , the distribution in the lower region approaches a linear profile typical of isotropic materials, consistent with the variation in Young’s modulus shown in Figure 2. On the other hand, Figure 5(c) and (d) show that the parameter $p$ exerts a minimal influence on the transverse stress ${\overline{\sigma }}_z(L/2,z)$ and the shear stress ${\overline{\tau }}_{\alpha z}(L/4,z)$ , indicating that the gradient of mechanical properties mainly affects the axial stresses and plays a limited role in these stress components.

Figure 6 shows the through-thickness distributions of displacements and stresses for different values of the porosity volume fraction $a$ . In Figure 6a, it is observed that the deflection $\overline{w}(L/2,z)$ increases significantly with $a$ , reflecting the progressive loss of structural stiffness caused by the reduction of the material’s effective load-bearing area. In Figure 6(b), the axial stress ${\overline{\sigma }}_{\alpha }(L/2,z)$ tends to concentrate in the upper region as porosity increases, while remaining practically constant in the lower region. Figure 6(c) and (d) show that the parameter $a$ has a limited influence on the transverse stress ${\overline{\sigma }}_z(L/2,z)$ and the shear stress ${\overline{\tau }}_{\alpha z}(L/4,z)$ , with only a slight increase recorded in the upper region of the beam. This pattern is consistent with the mechanical effect of porosity, which has a stronger impact on regions subjected to higher normal stresses.

Conclusions

A two-dimensional numerical solution for the static analysis of functionally graded porous curved beams under various boundary conditions has been presented in this work. The proposed model is based on the general 2D elasticity theory for curved beams. The strong form of the governing equations is solved via the Differential Quadrature Method (DQM), thus ensuring interlaminar continuity and boundary conditions. Validation against results obtained by commercial finite element software, higher-order shear deformation theories (HSDT) reported in the literature, assesses the model’s accuracy. The following conclusions can be drawn:

1. The developed model provides accurate distributions of displacements and stresses in functionally graded porous curved beams under different boundary conditions and material configurations. It outperforms HSDT-based approaches in terms of accuracy, particularly in cases with high curvature, and shows excellent agreement with FEM solutions.

2. From a computational standpoint, the DQM approach is more efficient than conventional finite element methods. Whereas FEM models require multilayer discretization to capture the continuous variation of material properties, the proposed approach directly captures the continuous variation of properties without the need for additional subdivisions, achieving comparable accuracy at a lower computational cost.

3. Geometric curvature, the material gradation ( $p$ ), and the porosity fraction ( $a$ ) have a significant influence on the structural response. The proposed model accurately predicts the redistributions of axial, transverse, and shear stresses induced by these parameters, which is valuable for the optimal design of curved structural components in composite materials.

The solutions serve as valuable benchmarks for validating reduced-order models (e.g., HSDT) and as design tools for advanced functionally graded components. Limitations include the linear elasticity assumption, excluding large deformation and nonlinear effects, and the current omission of multilayered or sandwich configurations, as well as thermal and dynamic phenomena.

Future research should extend this framework to dynamic and thermal analyses, laminated or sandwich structures, and multiphysics loading scenarios.