Evaluation of Stress Distribution of Isotropic, Composite, and FG Beams with Different Geometries in Nonlinear Regime via Carrera-Unified Formulation and Lagrange Polynomial Expansions

Abstract

In this study, the geometrically nonlinear behaviour caused by large displacements and rotations in the cross sections of thin-walled composite beams subjected to axial loading is investigated. Newton–Raphson scheme and an arc length method are used in the solution of nonlinear equations by finite element method to determine the mechanical effect. The Carrera-Unified formulation (CUF) is used to solve nonlinear, low or high order kinematic refined structure theories for finite beam elements. In the study, displacement area and stress distributions of composite structures with different angles and functionally graded (FG) structures are presented for Lagrange polynomial expansions. The results show the accuracy and computational efficiency of the method used and give confidence for new research.

1. Introduction

In the last decade, engineering materials have been classified into composite materials, smart materials, micro and nanomaterials, and gossamer space structures. In this classification, materials are distinguished by their tailoring capabilities. These are mechanical, electrical, magnetic, dimensional changes, and reduction in mass density. One-dimensional (1D) structures that exist due to the use of these new engineering materials include cables and beams, and two-dimensional (2D) structures include many linear and nonlinear structural theories have been and continue to be proposed for membranes, plates, and shells [1]. The first study was conducted by Euler in 1744 to evaluate the elastic limits in determining the nonlinear post-buckling and bending behaviour of 1D and 2D structures [2]. In this study, instead of the small deflection approximation, the exact expression for curvature is used. Subsequent studies are based on the Timoshenko beam theory. This exact solution analysis resulted in a post-buckling curve that rises slowly enough that there is no significant increase in load-carrying capacity until the deformations increase significantly. Then, studies with the exact solution, Barsoum and Gallagher [3] presented the buckling and post-buckling stability of two symmetrical I-beams subjected to conservative loads and Woolcock and Trahair [4] studied theoretically and experimentally the post-buckling behaviour of thin-walled I-beams for different boundary conditions. The proposed nonlinear theories of beam structures are based on the Timoshenko beam theory, which assumes a uniform shear distribution with the effects of rotatory inertia across the beam’s cross section, see, for example, the works in [5,6,7,8]. In general, an approximate solution has been proposed for complex structures of exact solutions, since analysis is difficult and time-consuming. The large displacement and rotation of beams are investigated using the finite element method and a current Lagrangian procedure [9,10,11,12,13,14]. Later, the nonlinear behaviour of mono- and bi-symmetric beams with Lagrange formulation is presented by Hsiao and Lin [15,16,17]. They demonstrated the effect of third-order terms of the nodal forces on the nonlinear buckling and post-buckling behaviour of three-dimensional (3D) beams through different models they created. Rizzetto et al. [18] adopted the finite element method for the analysis of buckling and post-buckling behaviour of isotropic and composite cylindrical shells. They employed a numerical time integration, as pulsating axial loading was considered. Gonçalves and Del Prado [19] studied the nonlinear vibration and dynamic instability of circular cylindrical shells employing the Galerkin and Ringe–Kutta methods. Moreover, geometrically nonlinear forced vibrations of circular cylindrical shells were investigated in detail by Popov [20]. Yu et al. [21,22] proposed the variational asymptotic method to separate the 3D nonlinear elasticity problem into the 2D linear cross-sectional analysis and the 1D nonlinear beam problem. With this method, cross-section and cross section stiffness constants are calculated under the effect of buckling and bending of beams with different geometry and material properties. Many refined beam theories have been proposed to determine nonlinear behaviour in more complex structures such as thin-walled beams and other higher-order phenomena that are subject to local effects [23,24,25,26,27]. One of these refined theories is the CUF [28,29]. CUF is a hierarchical formulation that sees the order of the structural model as input to the analysis and therefore does not need special formulations to obtain any refined models. Any higher-order theory can be used, thanks to fundamental nuclei that are independent of the order of expansion used to express unknown displacements in cross section. Thus, it provides accuracy and calculation efficiency in obtaining linear and nonlinear structural analysis. In the literature, CUF is used to solve many structural problems with different linear and nonlinear material properties [30,31,32,33,34,35,36,37,38,39,40,41,42].

CUF formulation has been proposed for geometric nonlinear problems by Pagani and Carrera [32]. Here, the combined formulation is first used the total Lagrangian approach to develop a nonlinear model. For the solution of this model, a Newton–Raphson linearization scheme and a method based on the arc length constraint are used together. In recent years, many studies have been done for large displacements with this approach. For example, the vibration of isotopic and anisotropic structures in studies [43,44,45], and the nonlinear bending and buckling behaviour of anisotropic beams and crust structures in studies [46,47,48] are investigated.

In this study, nonlinear displacement and stress distributions in the cross-section effect of the isotropic beam, composite beam placed at different angles, and FG beams are investigated by CUF theory. The results obtained are confirmed by previous studies [47,49].

2. Preliminaries

2.1. Constitutive Relations

In this study, the cross section of each beam is located in the xz plane of the Cartesian reference system and has a length L along the y-axis. Displacement vector is given below:

$$\mathbf{u}={\left\{\begin{array}{ccc}{u}_x& u_y& u_z\end{array}\right\}}^T$$

(1)

The strain $\epsilon$ and stress $\sigma$ tensors are expressed:

$$\epsilon ={\left\{\begin{array}{cccccc}{\epsilon }_{xx}& {\epsilon }_{yy}& {\epsilon }_{zz}& {\epsilon }_{xz}& {\epsilon }_{yz}& {\epsilon }_{xy}\end{array}\right\}}^T$$

(2)

and stress tensors are expressed, respectively,

$$\sigma ={\left\{\begin{array}{cccccc}{\sigma }_{xx}& {\sigma }_{yy}& {\sigma }_{zz}& {\sigma }_{xz}& {\sigma }_{yz}& {\sigma }_{xy}\end{array}\right\}}^T$$

(3)

After that, using the Green–Lagrange strain component and hypotheses of small deformations, the geometric relationships are described with the following equations:

$$\epsilon ={\epsilon }_l+{\epsilon }_{nl}=\left({\mathbf{b}}_l+{\mathbf{b}}_{nl}\right)u$$

(4)

where ${\mathbf{b}}_l$ and ${\mathbf{b}}_{nl}$ are linear and nonlinear differential operators, see the work in [32] for details on these operators. The material properties of composite structures with different angles and FG layered structures are taken into account with constitutive equations. Generalized Hooke’s law for homogeneous, linear, and elastic materials is given below:

$$\sigma =\mathbf{C}\epsilon$$

(5)

where C is a material matrix, and its open form [28] and examination with finite elements can be found in [50,51].

2.2. Adopted Refined Beam Theory

CUF assumes that the 3D displacement field ($\mathbf{u}(\mathit{x},\mathit{y},\mathit{z})$) can be expressed as the generalized displacement ${\mathbf{u}}_s$ (depending on the y coordinate) and the arbitrary functions of the section $F_s(x,z)$ (depending on the x and z coordinates). Accordingly, the displacement field for 1D theories is

$$\mathbf{u}(x,y,z)=F_s(x,z){\mathbf{u}}_s\left(y\right)s=1,2,\dots ,M$$

(6)

where M and s represent the number of terms used in the expansion and the sum of these terms, respectively. The $F_s$ section function determines the class of the 1D CUF model. With CUF, Taylor-like (TE) and Lagrange expansion (LE) section functions can be used. There have been many studies in the literature using both extensions with CUF. For TE models, M determines the polynomial degree, and if M is equal to 1, this model calculates with first-order shear theory. As M increases, the number of terms used increases, and the effects included in the calculation increase. For LE models, the number of points determines the order of the polynomial. For example, these polynomials can be used as linear three (L3), four (L4) points, quadratic six (L6), nine (L9), and cubic sixteen (16) points in the CUF framework. The LE model has been used specifically to determine the behaviour of layered composite and FG materials [29,39,47,50,52,53,54]. In this study, L4 and L9 polynomials are used. In the finite element analysis, classical four-node (B4) beam elements are used along the beam axis, providing a cubic approach. The choice of this element in CUF theory does not depend on the choice of sectional functions. FEM is used to separate the beam axis along y.

Accordingly, the generalized displacement vector ${\mathbf{u}}_s\left(y\right)$. $N_j$ and p represent the jth shape function and the order of the shape functions, respectively, where j represents the sum.

$${\mathbf{u}}_s\left(y\right)=N_j\left(y\right){\mathbf{q}}_{sj}j=1,2,\dots ,p+1$$

(7)

The vector of FE node parameters is given below with ${\mathbf{q}}_{\mathbf{sj}}$.

$${\mathbf{q}}_{\mathbf{sj}}={\left\{\begin{array}{ccc}{q}_{u_{x_{sj}}}& q_{u_{y_{sj}}}& q_{u_{z_{sj}}}\end{array}\right\}}^T$$

(8)

Details on the $N_j$ shape functions can be found in [29,50].

2.3. Nonlinear FE Equations

In an elastic system in equilibrium, the sum of the virtual changes of the strain energy caused by any arbitrary infinitesimal virtual displacements under the influence of internal and external forces is zero.

$$\delta L_{int}-\delta L_{ext}=0$$

(9)

$L_{int}$ represents work done by deformations and $L_{ext}$ represents work done by external forces. The work done by deformations (strain vector ($\epsilon$)) can be written in terms of stress and strain.

$$\delta L_{int}={\int }_V\delta {\epsilon }^T\sigma dV$$

(10)

Here, V is the initial body volume. Equation (4) can be written in terms of generalized node unknowns ${\mathbf{q}}_{sj}$ using Equations (6) and (7):

$$\epsilon =({\mathbf{B}}_l^{sj}+{\mathbf{B}}_{nl}^{sj}){\mathbf{q}}_{sj}$$

(11)

When the relevant equations are written in place, ${\mathbf{B}}_l$ and ${\mathbf{B}}_{nl}$ show matrices consisting of displacement, section function, and shape function. The name for brevity is not given, details can be found in [32]. The virtual variations of the strain tensor component can be written using the Green-Lagrange strain component and the small deformation hypothesis.

$$\delta \epsilon =\delta ({\mathbf{B}}_l^{\tau i}+{\mathbf{B}}_{nl}^{\tau i}){\mathbf{q}}_{\tau i}=({\mathbf{B}}_l^{\tau i}+2{\mathbf{B}}_{nl}^{\tau i})\delta {\mathbf{q}}_{\tau i}$$

(12)

where the transpose of the $\delta \epsilon$ tensor is taken,

$$\delta {\epsilon }^T=\delta {\mathbf{q}}_{\tau i}^T{({\mathbf{B}}_l^{\tau i}+2{\mathbf{B}}_{nl}^{\tau i})}^T$$

(13)

Here, for the sake of convenience, the indexes of the shape and cross-section functions have been expanded as follows.

$$\tau ,s=1,2,\dots ,Mi,j=1,2,\dots ,p+1$$

(14)

Substituting Equations (5) and (13) into Equation (10) yields

$$\delta L_{int}=\delta {\mathbf{q}}_{sj}^T\left({\int }_V{({\mathbf{B}}_l^{sj}+2{\mathbf{B}}_{nl}^{sj})}^{T}C({\mathbf{B}}_l^{\tau i}+{\mathbf{B}}_{nl}^{\tau i})dV\right){\mathbf{q}}_{\tau i}$$

(15)

where,

$${\mathbf{K}}_{\mathbf{S}}^{\mathbf{ij}\tau \mathbf{s}}={\int }_V{({\mathbf{B}}_l^{sj}+2{\mathbf{B}}_{nl}^{sj})}^{T}C({\mathbf{B}}_l^{\tau i}+{\mathbf{B}}_{nl}^{\tau i})dV$$

(16)

where ${\mathbf{K}}_{\mathbf{S}}^{\mathbf{ij}\tau \mathbf{s}}$ is the secant stiffness matrix, and the first term of this matrix represents the linear component, the next two terms represent first-order nonlinear components, and the last term represents the second-order nonlinear component. The secant stiffness matrix $\mathbf{K}$ is not symmetric and the mathematical and practical disadvantages of this situation are detailed in [32]. In addition, in the same study, the asymmetric form of the secant stiffness matrix used in the linearization of geometric stiffness terms is also included. In the solution of this nonlinear system under the influence of external loads, it is necessary to linearize the virtual change of the strain energy.

$$\delta {\mathbf{q}}_{sj}{F}_{ext}^{js}=\delta {\mathbf{q}}_{sj}^T{\mathbf{K}}_{\mathbf{S}}^{\mathbf{ij}\tau \mathbf{s}}{\mathbf{q}}_{\tau i}$$

(17)

$$F_{ext}^{js}={\mathbf{K}}_{\mathbf{S}}^{\mathbf{ij}\tau \mathbf{s}}{\mathbf{q}}_{\tau i}$$

(18)

For this purpose, a Newton–Raphson linearized incremental scheme and an arc–length constraint relationship are utilized.

$$\begin{array}{ccc}\delta \left(\delta L_{int}\right)\hfill & =\hfill & {\int }_V\left(\delta \left(\delta {\epsilon }^T\right)\sigma +\left(\delta {\epsilon }^T\delta \sigma \right)\right)dV\hfill \\ \\ \hfill & =\hfill & \delta \delta {\mathbf{q}}_{\tau i}^T\left(K_0^{ij\tau s}+K_{T_1}^{ij\tau s}+K_{\sigma }^{ij\tau s}\right)\delta {\mathbf{q}}_{sj}\hfill & =\hfill & \delta \delta {\mathbf{q}}_{\tau i}^T{K}_T^{ij\tau s}\delta {\mathbf{q}}_{sj}\hfill \end{array}$$

(19)

Here, ${\mathbf{K}}_{\mathbf{0}}^{\mathbf{ij}\tau \mathbf{s}}$ is the linear stiffness matrix and ${\mathbf{K}}_{\sigma }^{ij\tau s}$ is the geometric stiffness matrix. ${\mathbf{K}}_{\mathbf{T}}^{\mathbf{ij}\tau \mathbf{s}}$ is the fundamental nucleus of the tangent stiffness matrix that including the nonlinear contribution derived from the linearization of Hooke’s law. More details on this linearization can be found in [32]. The equation after linearization (18) is obtained:

$$F_{ext}^{i\tau }={\mathbf{K}}_{\mathbf{T}}^{\mathbf{ij}\tau \mathbf{s}}{\mathbf{q}}_{sj}$$

(20)

3. Numerical Results

In the numerical analysis section of this study, several problems are addressed to evaluate the nonlinear buckling analysis behaviour of beams with CUF theory: (1) thin-walled single-cell isotropic box beam, (2) thin-walled single-cell composite box beam, (3) thin-walled single-cell FG box beam, (4) thin-walled two-cell composite box beam, and (5) thin-walled two-cell FG box beam. Here, particular focus is on the displacements and stress behaviour of these structures under clamped-free boundary conditions. The properties of orthotropic composite materials with different fiber angles are taken as $E1=155$ GPa, $E2=E3=15.5$ GPa, Poisson’s ratio $\nu =0.25$ and shear modulus $G12=G13=9.3$ GPa, $G23=7.75$ GPa. The material properties of FG for ceramic (c) and metal (m) components are $E_c=380$ GPa, $E_m=70$ GPa, ${\rho }_c=3950$ kg/m${}^3$, ${\rho }_m=2700$ kg/m${}^3$, and $\nu =0.3$. In this study, the beam with a cross section of the cartesian reference system in the xz plane is located along the y-axis and has a length of $\mathrm{L}=8$ m. Furthermore, 12-B4 elements along the beam axis with CUF formulation and L4 and L9 point extensions are used in LE model. Studies have been validated with the literature.

3.1. Thin-Walled Single-Cell Isotopic Box Beam

The geometric properties of the cantilever thin-walled single-cell isotropic box beam are given in Figure 1, where $a_1=0.2$ m, $a_2=0.1$ m, and $t=0.005$ m. The material of the beam is ceramic and the material properties are given above. In this section, when we examine the isotropic nonlinear box beam model presented by Kvaternik [49] for the clamped-free boundary condition in order to verify our analysis, we obtain similar results as given in

Table 1. The critical buckling load (MN) of the isotropic box beam in the nonlinear (NL) regime.

	Ref. [49]	$\mathbf{CUF},\mathbf{LE}$
critical buckling load	$0.85837$	$0.85837$

. Here, the beam material is 100% ceramic.

Large Deflection of Cantilever Thin-Walled Single-Cell Isotopic Box Beam for Post-Buckling

In this section, the cantilever thin-walled cross section isotropic beam that is subject to large deflection due to an axial loading, is examined. Figure 2 shows the equilibrium curves of the box beam subjected to axial loading at $(-a_2/2,L,0)$ for models with polynomial degrees L4 and L9. The displacement distributions of L4 and L9 in the NL regime is similar. In Figure 3, the distribution of axial and shear stress components for two different ${\mathrm{P}}^{*}$ values of L4 and L9 of an isotropic box beam subjected to axial loading for the NL regime is shown at $x=-a_2/2$ and $y=L/2$ at all z values. The variation of the axial stress components along the z-axis is similar for both loads at L4 and L9. When the L4 beam theory is used, the distribution character and values of the shear stress component are different from L9. As the order of the polynomial of the Lagrange expansion functions changes from L4 to L9, the distribution character changes along the z-axis.

Figure 4 shows the distributions of axial and shear stress components in L and Nl regimes for different ${\mathrm{P}}^{*}$ values in L9. The distributions of the axial stress components for ${\mathrm{P}}_6^{*}$ are close for the linear (L) and the NL regimes. For ${\mathrm{P}}_{18}^{*}$ and ${\mathrm{P}}_{34}^{*}$, the distribution of the axial stress component is the same along the z-axis in the L regime. In the NL regime, it changes linearly along the z-axis of the beam. The shear stress component has a similar distribution for all ${\mathrm{P}}^{*}$ values in the L regime. For ${\mathrm{P}}_{34}^{*}$ only, the stress component acts at smaller levels in the lower parts of the beam. In the NL regime, the distribution for all ${\mathrm{P}}^{*}$ values is similar.

3.2. Thin-Walled Single-Cell Composite Box Beam

In this section, the geometric properties of the thin-walled composite box beam are taken into account as in Figure 1, while the material distribution of that beam is shown in Figure 5. Orthotropic materials with alignment angle as [0${}^{°}$/90${}^{°}$/0${}^{°}$] cross-ply are used here. The properties of these materials are detailed above.

Large Deflection of Cantilever Thin-Walled Single-Cell Composite Box Beam for Post-Buckling

In this section, the behaviour of the cantilever thin-walled single-cell composite box beam subjected to large deflection due to axial loading is examined. Figure 6 shows the equilibrium curves in NL regimes in beam models with polynomial degrees L4 and L9 subjected to axial loading effect. Accordingly, the displacement distribution of the box beam at L4 and L9 is similar up to $u_z/L=0.5$. Afterwards, L4 and L9 create a similar curve, but the L4 model is effective with slightly larger levels than L9.

Figure 7 shows the distribution of axial and shear stress components for two different ${\mathrm{P}}^{*}$ values of L4 and L9 of the composite box beam subjected to axial loading for the NL regime. Here, the distribution along z-axis is analyzed for $x=-a_2/2$ and $y=L/2$. The variation of the axial stress component is similar for both loads at L4 and L9. The distribution of the shear stress component along the z-axis at different load effects is similar for L4, but varies for L9.

Figure 8 shows the distributions of axial and shear stress components in L and Nl regimes for different ${\mathrm{P}}^{*}$ values in L9. The distributions of the axial stress components for ${\mathrm{P}}_6^{*}$ are close to each other for L and NL regimes. For ${\mathrm{P}}_{18}^{*}$ and ${\mathrm{P}}_{30}^{*}$, the distribution of the axial stress component is similar along the z-axis in the L regime. In the NL regime, the axial stress along the z-axis changes the behaviour of the beam towards the upper parts at ${\mathrm{P}}_{18}^{*}$ and ${\mathrm{P}}_{30}^{*}$ loads. As seen in Figure 8, when the shear stress is examined, for all load values in linear behaviour, the distribution along the z-axis occurs similarly. In the NL behaviour, as the load increases, the intensity of the stress increases at the lower parts of the beam and decreases at the upper parts.

3.3. Thin-Walled Single-Cell FG Box Beam

In this section, the mechanical behaviour of a cantilever thin-walled single-cell FG beam subjected to large deflection under an axial load is investigated. The geometric properties of the FG beam are as in Figure 1, and the cross section of the FG beam consisting of 5 layers is shown in Figure 9. The outer layer of the beam is ceramic, while the inner layer is metal. The properties of c and m materials are detailed above. The volume fraction of the c material of the beam follows the power law as [55]

$$B{V}_c(x,z)={\left(1-\frac{x,z}{t}\right)}^m$$

(21)

for the volume fraction of the m material of the beam,

$$V_m(x,z)=1-V_c(x,z)$$

(22)

where m is the compositional gradient exponent in the x- and z-directions, and t is the beam thickness. The Mori–Tanaka scheme [56] is used for grading the material properties according to this volume fraction. Details are not included in the name of brevity.

Large Deflection of Cantilever Thin-Walled Single-Cell FG Box Beam for Post-Buckling

In this section, two different compositional gradient exponent values of the FG beam are taken, $\mathrm{m}=0.5$ and 5, representing the ceramic-rich and metal-rich composition, respectively. For these compositional gradient exponents, the displacement and stress distributions of the FG beam subjected to the axial load effect are investigated in the NL regime. Figure 10 shows the equilibrium curves of the FG box beam with polynomial degrees L4 and L9 subjected to axial loading effect for different compositional gradient exponents. The behaviour characteristics of the displacement component at L4 and L9 are similar for both compositions. However, ${\mathrm{P}}^{*}$ levels increase approximately 2.5 times in ceramic-rich composition compared to metal-rich composition.

Figure 11 shows the distribution of axial and shear stress components for two different ${\mathrm{P}}^{*}$ values of L4 and L9 of ceramic-rich FG box beam subjected to axial loading for the NL regime. The variation of the axial stress component is similar for both loads at L4 and L9. The distribution of the shear stress component along the z-axis at different load effects is similar for L4, but for L9, as the load increases, a sudden rise occurs in the lower parts of the beam. Figure 12 shows the distribution of axial and shear stress components for two different ${\mathrm{P}}^{*}$ values of L4 and L9 of the metal-rich FG box beam subjected to axial loading for the NL regime. The variation of the axial stress component along the z-axis is similar for L4 and L9, but as the load increases, the level of the stress component decreases at the upper parts of the beam, and the distribution changes. The distribution of the shear stress component is similar for both loads in L4 and L9.

Figure 13 shows the distribution of the axial and shear stress components in L and NL regimes for different ${\mathrm{P}}^{*}$ values of the ceramic-rich FG box beam in L9. The distributions of the axial stress components for ${\mathrm{P}}_6^{*}$ are close to each other for L and NL regimes. For ${\mathrm{P}}_{18}^{*}$ and ${\mathrm{P}}_{34}^{*}$, the distribution of the axial stress component is similar along the z-axis in the L regime. In the NL regime, axial stress along the z-axis changes the behaviour of the beam towards the top at ${\mathrm{P}}_{18}^{*}$ and ${\mathrm{P}}_{34}^{*}$ loads. As seen in Figure 13, when the shear stress is examined, the distribution along the z-axis is similar for all load values in L behaviour. In the NL regime, while ${\mathrm{P}}_6^{*}$ and ${\mathrm{P}}_{18}^{*}$ are similar, at ${\mathrm{P}}_{34}^{*}$ the levels decrease at the bottom of the beam. The distribution of axial and shear stress components in L and NL regimes for different ${\mathrm{P}}^{*}$ values of the metal-rich FG box beam is shown in Figure 14 for L9. The distribution of the axial stress component in L and NL regimes is similar at ${\mathrm{P}}_6^{*}$. However, for ${\mathrm{P}}_{18}^{*}$ and ${\mathrm{P}}_{34}^{*}$, the distribution of the axial stress component decreases towards the upper part of the beam. In Figure 14, when the shear stress is examined, the linear distribution and nonlinear distribution have similar distribution characters for all loads, but the levels of the stress component increase in the NL distribution.

3.4. Thin-Walled Two-Cell Composite Box Beam

The geometric properties of the thin-walled two-cell box beam are given in Figure 15, where $a_1=0.2$ m, $a_2=0.1$ m, and $t=0.005$ m, while the material distribution of that beam is shown in Figure 16. As can be seen here, orthotropic material with different alignment angles is used. The composite material of that beam is arranged as [0${}^{°}$/90${}^{°}$/0${}^{°}$] cross-ply, and orthotropic properties of these materials are detailed above.

Large Deflection of Thin-Walled Two-Cell Composite Box Beam for Post-Buckling

In this section, the buckling behaviour of the cantilever twin-cell composite box beam, where large deflections occur under axial load, is examined. The equilibrium curve of a two-cell box beam with polynomial degree L9 is shown in Figure 17 for both NL and L regimes.

Figure 18 shows the distributions of axial and shear stress components at two different ${\mathrm{P}}^{*}$ values of the composite box beam subjected to axial loading for the NL regime at L9. The variation of the axial stress component is similar for both loads. However, the stress level increases at the upper parts of the beam as the load increases. The distribution of the shear stress component along the z-axis at different load effects is similar.

Figure 19 shows the distributions of axial and shear stress components of beam in L and Nl regimes for different ${\mathrm{P}}^{*}$ values in L9. The distributions of the axial stress components are similar for all loads in the L regime. However, the ${\mathrm{P}}_6^{*}$ load affects the top of the beam at lower levels. The distribution of the same stress component in the NL regime is similar for all loads. In the shear stress distribution, the distribution does not change for different load values both in the L regime and in the NL regime.

3.5. Thin-Walled Two-Cell FG Box Beam

The geometric properties of the two-cell beam are as in Figure 15, and the cross section of the FG beam consisting of three layers is shown in Figure 20. Here, the outer layer of the beam represents the ceramic, the inner layer represents the metal and the middle layer represents the graded region. In this study, the buckling behaviour of the FG beam subjected to axial loading for two compositional gradient exponents is investigated. Equations (21) and (22) are used in the FG region grading and the ceramic and metal material properties are as given above.

Large Deflection of Thin-Walled Two-Cell FG Box Beam for Post-Buckling

In this section, the buckling behaviour of the cantilever twin-cell FG beam with large deflections subjected to axial loading is investigated for different compositional gradient exponents. The equilibrium curve of two-cell box FG beam with polynomial degree L9 is shown in Figure 21 for both NL and L regimes. The curves are plotted along the beam’s y-axis at the position $x=0$ and $z=0$. Here, the displacement distribution for both compositional gradient exponents is similar.

Figure 22 shows the distribution of axial and shear stress components for two different ${\mathrm{P}}^{*}$ values in L9 of the ceramic-rich FG box beam subjected to axial loading for the NL regime. As seen in Figure 22a,b, although the distribution character is similar, when the load increases, the levels increase at the upper parts of the beam. The distribution of the shear stress component along the z-axis at different load effects is similar. Figure 23 shows the distribution of axial and shear stress components for two different ${\mathrm{P}}^{*}$ values at L9 of the metal-rich FG box beam subjected to axial loading. The distribution character of the axial stress component along the z-axis is similar for different loads, but the levels of the stress component increase in the upper parts of the beam as the load increases. The distribution of the shear stress component is similar for both loads.

Figure 24 shows the distributions of the axial and shear stress components of the ceramic-rich FG box beam in L and NL regimes for different ${\mathrm{P}}^{*}$ values in L9. The distributions of the axial stress components are similar for all loads in the L regime. The distribution of the axial stress component in the NL regime is similar for all loads. However, the load ${\mathrm{P}}_6^{*}$ effects at lower levels. The shear stress distribution does not change for different load values both in the L regime and in the NL regime.

Figure 25 shows the distributions of the axial and shear stress components of the metal-rich FG box beam in L and NL regimes for different ${\mathrm{P}}^{*}$ values in L9. The distributions of axial stress components are similar for all loads in both L and NL regimes. The shear stress distribution does not change for different load values both in the L regime and in the NL regime. The levels of the axial and shear stress components change depending on the change of the composition (Figure 24 and Figure 25).

4. Conclusions

In this study, the NL buckling behaviour of thin-walled single-cell isotropic box beam, single-cell composite box beam, single-cell FG box beam, two-cell composite box beam, and two-cell FG box beam structures is investigated. Here, the focus is particularly on displacement and stress distributions. In the numerical model, using the CUF combined formulation of geometrically NL theories, the kinematics of the generic 1D model is expressed as an arbitrary expansion of the primary displacement unknowns. In addition, correct stress/strain distributions are obtained by using Lagrange extensions to formulate layerwise models in the CUF area. As a result, the proposed model in the light of CUF in this study has been demonstrated to provide accurate and effective results in determining the geometrically NL behaviour of isotropic, composite, and FG beam structures with different geometric properties subjected to axial loading.