Nonlinear bending analysis of matrix cracked functionally graded graphene origami-enabled auxetic metamaterial plates using the C 0 -HSDT meshless method

Abstract

The auxetic metamaterials characterized by negative Poisson’s ratio (NPR) are attracting significant interest due to their extraordinary and captivating mechanical properties. This paper presents a meshless model based on the Improved Moving Kriging (IMK) interpolation to investigate the nonlinear bending of functionally graded (FG) graphene origami (GOri)-enabled auxetic metamaterials (GOEAM) plates with matrix cracks. The model is based on the C⁰-higher order shear deformation theory (C⁰-HSDT) framework with seven variables and incorporates the von-Karman nonlinearity. The material properties of FG-GOEAM plate are determined by a genetic programming (GP)-assisted micromechanics model, while the degraded stiffness of cracked layers is incorporated via a self-consistent micromechanics model. The convergence and effectiveness of the proposed method is examined by comparing it with existing literature results. Firstly, the convergence and validity of the method presented in this paper are verified through comparison with existing literature. Then, A comprehensive study is conducted to examine key parametric effects on the tunability of nonlinear bending deflection and normal stress in matrix-cracked FG-GOEAM plates. The results show that higher GOri content enhances the stiffness of the FG-GOEAM plate, leading to lower deflection, whereas increased GOri folding degree results in greater deflection. Additionally, the deflection and in-plane stress of the matrix-cracked FG-GOEAM plate are higher than those of the intact plate, while shear stress is lower, with these effects becoming more pronounced as crack density increases.

Introduction

Since the graphene was discovered by Novoselov et al. ¹ in 2004, it has attracted significant attention from the research and engineering communities due to its superior thermal, mechanical, and electrical properties compared to traditional carbon fiber materials. ² Graphene and its derivatives have found applications as fillers in polymer composite materials in various fields, including flexible electronics, biosensors, aerospace, and civil engineering.^3,4 Graphene nanoplatelets (GPLs), ⁵ composed of single-layered atoms in a two-dimensional carbon nanostructure, possess exceptional features such as extremely high surface area, high aspect ratio, and strong matrix adhesion. Experimental studies conducted by Rafiee et al. ⁶ have demonstrated that the addition of a small weight fraction (0.1 ± 0.002%) of GPLs as fillers in a pure epoxy matrix can lead to significant enhancements in the Young’s modulus by approximately 31% and fracture toughness by over 50% compared to the pure epoxy matrix alone. However, it should be noted that excessive incorporation of graphene into the matrix can result in aggregation issues, impeding the full exploitation of its reinforcing effects. To address this issue, Yang et al.^7,8 introduced the concept of functional graded into graphene platelet-reinforced composite materials in 2017, resulting in the development of functionally graded graphene platelet-reinforced composites (FG-GPLRC). Subsequently, within a short span of a few years, a significant number of research papers have emerged focusing on the mechanical properties of FG-GPLRC beams and plates,^9–13 including both linear and nonlinear bending, free and nonlinear vibration, buckling and post-buckling behavior. Recently, Ni et al. ¹⁴ conducted Taylor series expansion and differential quadrature (DQ) method to analyze the nonlinear vibration of FG-GPLRC circular membrane with dielectric properties, where the effective medium theory and rule of mixture are used to evaluate the material properties of the membrane. Using the same method, Ni et al.^15,16 also studied the nonlinear vibration problem of FG-GPLRC dielectric membranes containing internal pores, employing a two-step hybrid mechanical model to determine the material properties of the multiphase composites. Additionally, based on first-order shear deformation theory (FSDT) and von Karman nonlinear theory, the nonlinear dynamics of cracked FG-GPLRC dielectric beams subjected to electrical field and mechanical excitation are investigated by Ni et al. ¹⁷

In recent years, metamaterials exhibiting a negative Poisson’s ratio (NPR) have emerged as novel materials with unconventional physical and mechanical properties, making them a prominent research area in the engineering community.^18,19 The graphene has proven to have highly tunable negative Poisson’s ratio (NPR) characteristics. ²⁰ Through the utilization of graphene origami (GOri) techniques, ²¹ Zhao et al. ²² developed a novel GOri-enabled auxetic metal metamaterials (GOEAM) with tunable NPR and enhanced mechanical properties. The material properties of the GOEAM, including Poisson’s ratio, Young’s modulus, and mass density can be accurately and efficiently estimated by genetic programming (GP)-assisted micromechanical models. ²³ Zhao et al.^24–27 investigated the tunable nonlinear bending, nonlinear free vibration, buckling and post-buckling behavior of FG-GOEAM beams by using the FSDT and DQ method. Lv et al. ²⁸ adopted the FSDT and two-step perturbation method to study the thermal postbuckling of the FG-GOEAM plate under Pasternak foundation. An et al. ²⁹ applied the hyperbolic shear deformation theory (HPSDT) and generalized finite difference method to examine the linear bending and buckling analysis of FG-GOEMA irregular plates. Murari et al.^30–33 combined FSDT and DQ method to study the linear bending, free vibration, nonlinear vibration, and postbuckling of FG-GOEMA beams immersed in fluid. Besides, Murari et al. ³⁴ investigated the vortex-induced vibration of FG-GOEMA plates attached to a circular cylinder in laminar flow by using the finite element method (FEM). From the above, it can be observed that there is no literature currently available reporting the nonlinear bending response of FG-GOEAM plates, which is one of the main objectives of this paper.

Matrix cracks are commonly observed in engineering structures such as shipbuilding, civil engineering, marine, and naval industries, typically resulting from the manufacturing process, fatigue, or excessive local stress concentrations during the structure’s service life. Matrix cracks are the primary form of initial damage in composite laminates, leading to a reduction in the stiffness and strength of the laminated composite plate. Laws et al.^35,36 proposed the self-consistent model to evaluate the properties of matrix-cracked composites. Based on the differential self-consistent method, Hoiseth ³⁷ developed an analytical expression for the effective stiffness of cross-ply laminates with transverse cracks, which shows a good agreement with FEM calculations. Subsequently, the self-consistent model was successfully employed to model laminated composites with matrix cracks. Fan and Wang^38–41 conducted research on the nonlinear bending, vibration, and buckling of the matrix cracked hybrid laminated plates containing CNTRC layers in the thermal environment or the plates on viscoelastic foundations. They utilized the third-order shear deformation theory (TSDT) and von Karman’s nonlinear theory, and employed the two-step perturbation technique to solve the governing equations. Lei et al.^42,43 used FSDT and element-free kp-Ritz method to investigate the linear bending, free vibration and buckling of matrix cracked hybrid laminated plates containing CNTRC layers. Based on the FSDT and element-free IMLS-Ritz method, Guo et al.^44,45 studied the flutter characteristics of matrix cracked FG-GPLRC plates and the vibration characteristics of rotating matrix cracked FG-GPLRC cylindrical shells.

The classical plate theory (CPT) relies mainly on the Kirchhoff-Love thin plate assumption, neglecting the transverse shear deformation effects. The FSDT considers the influence of shear deformation to accommodate medium-thick plates. However, it yields inaccurate shear strain/stress distributions and fails to satisfy the zero traction boundary conditions on the top and bottom surfaces of the plate, requiring the introduction of shear correction factors. Consequently, researchers have developed various higher-order shear deformation theories (HSDT), including TSDT by Reddy, ⁴⁶ sinusoidal shear deformation theory (SSDT) by Touratier, ⁴⁷ inverse hyperbolic shear deformation theory (IHSDT) by Grover et al. ⁴⁸ , and exponential shear deformation theory (ESDT) by Karama et al. ⁴⁹ . These HSDTs eliminate the need for shear correction factors and can accurately predict transverse shear stresses, making them suitable for thick plates. However, they require the generalized displacements to have C ¹ continuity, whereas traditional Lagrangian interpolation-based finite element shape functions only possess C⁰ continuity. To overcome this limitation, Shankara et al. [50] introduced two additional artificial variables in the displacement field, proposing the C⁰-TSDT with seven variables. They combined this theory with FEM to verify its effectiveness in analyzing the free vibration problems of laminated plates. Afterward, the C⁰-HSDT theory, combined with meshless method, FEM and isogeometric analysis (IGA), has been extensively applied in the mechanical performance analysis of laminated plates,^51,52 FG-CNTRC plates, ⁵³ FG-GPLRC plate^54,55 and magneto-electro-elastic functionally graded porous plates. ⁵⁶

The Element-Free Galerkin (EFG) method,⁵⁸⁷ a significant branch of meshless methods, typically constructs shape functions using the Moving Least Squares (MLS) approximation. However, overcoming the Kronecker delta condition in MLS necessitates the use of penalty functions, Lagrangian interpolation, and full transformation methods, ⁵⁸ which can adversely affect computational efficiency and accuracy. Gu et al. ⁵⁹ introduced the MK interpolation function to solve the Kronecker delta issue, allowing for boundary conditions to be applied similarly to the FEM. Nonetheless, the Gaussian function used in MK interpolation depends on the parameter θ, whose optimal value is notoriously difficult to determine. To address this, Thai and workers^60–62 proposed an IMK meshless method using quartic spline functions to define the correlation function. Combined with refined plate theory, this method has been demonstrated to provide stability and enhanced computational accuracy in the static bending, free vibration, and buckling analysis of isotropic and sandwich plates.

So far, no previous work has been done on the nonlinear bending analysis of FG-GOEAM plates with matrix cracks. Hence, this paper aims to explore the nonlinear bending response of matrix cracked FG-GOEAM plates using the C⁰-HSDT in conjunction with an improving moving kriging meshless (IMK) method. The GP-assisted and self-consistent micromechanical models are used to model the material properties of FG-GOEAM plate and the properties of matrix cracked layers, respectively. The nonlinear behavior is considered by employing the von-Karman strain assumptions, and the meshless nonlinear governing equations are solved numerically by the Newton-Raphson iterative method. The numerical illustrations show the effects of various parameters on the nonlinear bending and stress of matrix cracked FG-GOEAM plates.

Problem formulations

The research described in this paper focuses on a multilayer plate with dimensions: length a in the x direction, width b in the y direction, and thickness h in the z direction. This plate is composed of multilayer GOEAMs, which are formed by the Cu metal matrix and Gori. The plate’s individual layers contain cracks modeled as elliptic cylindrical cavities. Each layer has the same thickness, and is assumed to be isotropic and homogeneous, but the content of GOri varies in a stepwise manner along the thickness direction of the plate. Based on the weight fraction (W_Gr) variation of GOri from layer to layer, four distribution patterns were designed: U-W_Gr pattern, X-W_Gr pattern, O-W_Gr pattern and V-W_Gr pattern. The distribution patterns are illustrated in Figure 1, where N_L represents the total number of layers. In the illustration, a darker color indicates a higher concentration of dispersed GOri in the corresponding layer, implying that it has a higher auxetic characteristic.

Figure 1.

Schematic configuration of matrix cracked FG-GOEAM plate.

Functionally graded metamaterial plates model

The volume fractions $V_{\mathrm{Gr}}^{(k)}$ of GOri in the kth (k = 1, 2, 3, …, N_L) layer of the four distribution patterns are determined by the following formula ²⁶ :

$$V_{\mathrm{Gr}}^{(k)}=\{\begin{array}{cc}\multicolumn{1}{c}{V_{\mathrm{Gr}}^{*}}& \mathrm{U}-W_{\mathrm{Gr}}\\ \multicolumn{1}{c}{2V_{\mathrm{Gr}}^{*}|2k-N_L-1|/N_L}& \mathrm{X}-W_{\mathrm{Gr}}\\ \multicolumn{1}{c}{2V_{\mathrm{Gr}}^{*}(1-|2k-N_L-1|/N_L)}& \mathrm{O}-W_{\mathrm{Gr}}\\ \multicolumn{1}{c}{V_{\mathrm{Gr}}^{*}(2k-1)/N_L}& \mathrm{V}-W_{\mathrm{Gr}}\end{array}$$

(1)

Where $V_{\mathrm{Gr}}^{*}$ represents the total volume fraction of the GOri, and is calculated by

$$V_{\mathrm{Gr}}^{*}=\frac{{\rho }_{\mathrm{Cu}}{W}_{\mathrm{Gr}}}{{\rho }_{\mathrm{Cu}}{W}_{\mathrm{Gr}}+{\rho }_{\mathrm{Gr}}(1-W_{\mathrm{Gr}})}$$

(2)

The Young’s modulus (E_c) and Poisson’s ratio (ν_c) can be determined by a combination of mixing rules and a GP-assisted micromechanical model ²³ :

$$\begin{array}{c}\multicolumn{1}{c}{E_{\mathrm{c}}=\frac{1+\xi \eta V_{\mathrm{Gr}}^{(k)}}{1-\eta V_{\mathrm{Gr}}^{(k)}}{E}_{\mathrm{Cu}}\times f_E(H_{\mathrm{Gr}},V_{\mathrm{Gr}}^{(k)},T)}\\ \multicolumn{1}{c}{v_{\mathrm{c}}=(v_{\mathrm{Gr}}{V}_{\mathrm{Gr}}^{(k)}+v_{\mathrm{Cu}}{V}_{\mathrm{Gu}}^{(k)})\times f_v(H_{\mathrm{Gr}},V_{\mathrm{Gr}}^{(k)},T)}\end{array}$$

(3)

$$\eta =\frac{(E_{\mathrm{Gr}}/E_{\mathrm{Cu}})-1}{(E_{\mathrm{Gr}}/E_{\mathrm{Cu}})+\xi },\xi =2(l_{\mathrm{Gr}}/t_{\mathrm{Gr}})$$

(4)

where $V_{\mathrm{Gr}}^{(k)}$ and $V_{\mathrm{Gu}}^{(k)}$ are the volume fraction of GOri nano-filler and Cu matrix, respectively, and $V_{\mathrm{Gu}}^{(k)}+V_{\mathrm{Gr}}^{(k)}=1$ . l_Gr and t_Gr are the length and thickness of graphene, respectively. The GP algorithm is employed to determine the expression of the modification coefficient.

$$\begin{array}{c}\multicolumn{1}{c}{f_E(H_{\mathrm{Gr}},V_{\mathrm{Gr}}^{(k)},T)=1.11-1.22V_{\mathrm{Gr}}^{(k)}-0.134(T/T_0)}\\ \multicolumn{1}{c}{+0.559V_{\mathrm{Gr}}^{(k)}(T/T_0)-5.5H_{\mathrm{Gr}}{V}_{\mathrm{Gr}}^{(k)}}\\ \multicolumn{1}{c}{+38H_{\mathrm{Gr}}(V_{\mathrm{Gr}}^{(k)}{)}^2-20.6H_{\mathrm{Gr}}^2(V_{\mathrm{Gr}}^{(k)}{)}^2}\end{array}$$

(5)

$$\begin{array}{c}\multicolumn{1}{c}{f_v(H_{\mathrm{Gr}},V_{\mathrm{Gr}}^{(k)},T)=1.01-1.43V_{\mathrm{Gr}}^{(k)}+0.165(T/T_0)}\\ \multicolumn{1}{c}{-16.8H_{\mathrm{Gr}}{V}_{\mathrm{Gr}}^{(k)}-1.1H_{\mathrm{Gr}}{V}_{\mathrm{Gr}}^{(k)}(T/T_0)}\\ \multicolumn{1}{c}{+16H_{\mathrm{Gr}}^2(V_{\mathrm{Gr}}^{(k)}{)}^2}\end{array}$$

(6)

where T is the ambient temperature, and T₀ = 300 K is the reference temperature. The H_Gr is the H atom coverage that is used to quantify the degree of folding in Gori.

Figure 2 illustrates the distribution of Poisson’s ratio in each layer of FG-GOEAM plates with different folding degrees. It can be observed that as the folding degree increases, more layers in the FG plates exhibit NPR characteristics. Additionally, it is evident that the surface layers of the X-W_Gr pattern exhibit higher NPR values, while the inner layers show positive values. In contrast, the O-W_Gr pattern demonstrates the opposite behavior. On the other hand, the U-W_Gr pattern has the same Poisson’s ratio in each layer. For the V-W_Gr pattern, the Poisson’s ratio gradually increases from the top layer (negative) to the bottom layer (positive). Therefore, the definition of a metamaterial beam is based on the presence of NPR layers in the FG plate.

Figure 2.

The Poisson’s ratio distributions along the direction of thickness for the FG-OEAM plates with various GOri folding degrees and distribution patterns (W_Gr = 2.5 wt%, T = 300 K).

Degraded stiffness of a cracked composite lamina

Matrix cracks in composite structures can result in a decrease in stiffness. A micromechanical model known as the Self-Consistent Model (SCM)^36,37 is utilized to assess the degraded stiffness. By employing a self-consistent estimation, the overall compliance matrices L^(k) of the kth layer in a matrix-cracked composite can be given as

$${\mathbf{L}}^{(k)}={\overline{\mathbf{L}}}^{(k)}+\frac{1}{4}\mathrm{\pi }\beta {\mathbf{\Lambda }}^{(k)}$$

(7)

where the compliance matrix ${\overline{\mathbf{L}}}^{(k)}$ represents the perfect composite lamina. The crack density parameter β is determined by the ratio of the layer thickness h_i to the distance between adjacent cracks d. Consequently, β provides an estimate of the average number of cracks within the length of h_i. When β = 1, the crack spacing is equal to the layer thickness. When β = 0, there are no cracks present. The matrix Λ^(k) contains only three non-zero components: ${\Lambda }_{22}^{(k)}$ , ${\Lambda }_{44}^{(k)}$ and ${\Lambda }_{66}^{(k)}$ . This indicates that the introduction of cracks affects only three compliance coefficients: $L_{22}^{(k)}$ , $L_{44}^{(k)}$ and $L_{66}^{(k)}$ . The remaining six terms in the compliance matrix L^(k) are unchanged, and their expressions can be written as follows:

$$\begin{array}{l}{L}_{11}^{(k)}={\overline{L}}_{11}^{(k)}=\frac{1}{E_c^{(k)}},L_{12}^{(k)}={\overline{L}}_{12}^{(k)}=-\frac{v_c^{(k)}}{E_c^{(k)}},L_{13}^{(k)}={\overline{L}}_{13}^{(k)}=-\frac{v_c^{(k)}}{E_c^{(k)}},\\ L_{23}^{(k)}={\overline{L}}_{23}^{(k)}=-\frac{v_c^{(k)}}{E_c^{(k)}},L_{33}^{(k)}={\overline{L}}_{33}^{(k)}=\frac{1}{E_c^{(k)}},L_{55}^{(k)}={\overline{L}}_{55}^{(k)}=\frac{2\left(1+v_c^{(k)}\right)}{E_c^{(k)}}.\end{array}$$

(8)

The other non-zero elements in the compliance matrix ${\overline{\mathbf{L}}}^{(k)}$ are

$${\overline{L}}_{22}^{(k)}=\frac{1}{E_c^{(k)}},{\overline{L}}_{44}^{(k)}=\frac{2(1+v_c^{(k)})}{E_c^{(k)}},{\overline{L}}_{66}^{(k)}=\frac{2(1+v_c^{(k)})}{E_c^{(k)}}$$

(9)

The $L_{44}^{(k)}$ can be given by the following equation:

$$L_{44}^{(k)}={\overline{L}}_{44}^{(k)}+\frac{1}{2}\beta (\beta L_{55}^{(k)}+\sqrt{{\beta }^2{\left(L_{55}^{(k)}\right)}^2+4{\overline{L}}_{44}^{(k)}{L}_{55}^{(k)}})$$

(10)

In order to obtain $L_{22}^{(k)}$ and $L_{66}^{(k)}$ , the following equation is introduced:

$$\begin{array}{c}\multicolumn{1}{c}{X_1={\overline{L}}_{22}^{(k)}{L}_{33}^{(k)}-{\left(L_{23}^{(k)}\right)}^2,X_2=L_{33}^{(k)}{\overline{L}}_{66}^{(k)},}\\ \multicolumn{1}{c}{X_3=L_{11}^{(k)}{L}_{33}^{(k)}-{\left(L_{13}^{(k)}\right)}^2,X_4=2(L_{12}^{(k)}{L}_{33}^{(k)}-L_{13}^{(k)}{L}_{23}^{(k)}).}\end{array}$$

(11)

By solving the equation below, Y can be obtained:

$$Y-X_1-\beta \sqrt{3Y\sqrt{X_{3}Y}+(X_2+X_4)Y-X_1\sqrt{X_{3}Y}}=0$$

(12)

Then, $L_{22}^{(k)}$ and $L_{66}^{(k)}$ can be given by the following equation:

$$L_{22}^{(k)}=\frac{Y+{\left(L_{23}^{(k)}\right)}^2}{L_{33}^{(k)}},L_{66}^{(k)}={\overline{L}}_{66}^{(k)}+(L_{22}^{(k)}-{\overline{L}}_{22}^{(k)})\sqrt{\frac{X_3}{Y}}$$

(13)

The relationship between the stiffness coefficients $Q_{\mathrm{ij}}^{(k)}$ and compliance coefficients $L_{\mathrm{ij}}^{(k)}$ are expressed as:

$$\begin{array}{l}{Q}_{11}^{(k)}=\frac{L_{22}^{(k)}}{L_{11}^{(k)}{L}_{22}^{(k)}-{\left(L_{12}^{(k)}\right)}^2},Q_{12}^{(k)}=\frac{-L_{12}^{(k)}}{L_{11}^{(k)}{L}_{22}^{(k)}-{\left(L_{12}^{(k)}\right)}^2},Q_{22}^{(k)}=\frac{L_{11}^{(k)}}{L_{11}^{(k)}{L}_{22}^{(k)}-{\left(L_{12}^{(k)}\right)}^2},\\ Q_{66}^{(k)}=\frac{1}{L_{66}^{(k)}},Q_{44}^{(k)}=\frac{1}{L_{44}^{(k)}},Q_{55}^{(k)}=\frac{1}{L_{55}^{(k)}}.\end{array}$$

(14)

Theoretical formulations

The C⁰-HSDT formulation

In accordance with the properties of the equivalent higher-order shear deformation theory with five field variables, the assumed displacement field of the plate can be expressed in a general form as

Follows:

$$\{\begin{array}{c}\multicolumn{1}{c}{u(x,y,z)=u_0(x,y)-z\frac{\partial w_0}{\partial x}+f(z){\beta }_x(x,y)}\\ \multicolumn{1}{c}{v(x,y,z)=v_0(x,y)-z\frac{\partial w_0}{\partial y}+f(z){\beta }_y(x,y)}\\ \multicolumn{1}{c}{w(x,y,z)=w_0(x,y)}\end{array}$$

(15)

where u, v and w represent the displacements at an arbitrary point in the x, y, and directions, respectively; u₀, v₀ and w₀ correspond to the displacements at midplane; θ_x and θ_y denote the rotations about the y-axis and x-axis of the midplane, respectively; f(z) is the transverse shear shape function, which can satisfy the traction-free boundary condition at the top and bottom plate surfaces. Up to now, a considerable number of HSDTs that contain different shear shape functions have been proposed. Table 1 lists several shear functions. There are studies in the literature^{48, 63} indicate that the inverse hyperbolic shear function offers higher computational accuracy compared to the exact 3D solutions. Therefore, this study selects it for analysis.

Table 1.

Several shear functions and its derivatives.

Source	$f(z)$	$f\prime (z)$
TSDT by Reddy 46	$z-\frac{4}{3}{z}^3/h^2$	$1-4z^2/h^2$
ESDT by Karama 72	$z{e}^{-2{(z/h)}^2}$	$(1-\frac{4}{h^2}{z}^2)e^{-2{(z/h)}^2}$
SSDT by Touratier 47	$\frac{h}{\pi }\sin \left(\frac{\pi }{h}z\right)$	$\cos \left(\frac{\pi }{h}z\right)$
IHSDT by Grover 48	${\sinh }^{-1}\left(\frac{3z}{h}\right)-\frac{6}{\sqrt{13}h}z$	$\frac{3}{\sqrt{9z^2+h^2}}-\frac{6}{\sqrt{13}h}$

Although meshless methods facilitate the construction of higher-order shape functions, there are still difficulties in applying the partial derivative of the displacement w_0,x, w_0,y with respect to the boundary conditions. For example, Auricchio et al. ⁶⁴ proposed the use of stream function formulations, which require incorporating corresponding neighboring nodes at the boundary nodes to enforce the boundary conditions. This approach leads to a rapid increase in the number of nodes in the problem domain, increasing computational complexity and sacrificing computational accuracy. Thus, to reduce the continuity from C ¹ to C⁰, an artificial constraint is imposed by setting −∂w₀/∂x = θ_x and −∂w₀/∂y = θ_y. This constraint increases the number of field variables from five to seven, namely u = [u₀, v₀, w₀, θ_x, θ_y, β_x, β_y]^T, the new displacement field model is given as ⁵⁰

$$\{\begin{array}{c}\multicolumn{1}{c}{u(x,y,z)=u_0(x,y)+z{\theta }_x+f(z){\beta }_x(x,y)}\\ \multicolumn{1}{c}{v(x,y,z)=v_0(x,y)+z{\theta }_y+f(z){\beta }_y(x,y)}\\ \multicolumn{1}{c}{w(x,y,z)=w_0(x,y)}\end{array}$$

(16)

Strain-displacement relations

In the ESL theory, by neglecting the transverse strain in the z-direction and assuming small strains, moderate rotations and large displacements, the nonlinear strain-displacement relations can be described in the von-Karman form as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\epsilon }}_{\mathrm{p}}=\left(\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_y}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}}\end{array}\right)=\left(\begin{array}{c}\partial u/\partial x\hfill \\ \multicolumn{1}{c}{\partial v/\partial y}\\ \multicolumn{1}{c}{\partial u/\partial y+\partial v/\partial x}\\ \multicolumn{1}{c}{\partial u/\partial z+\partial w/\partial x}\\ \multicolumn{1}{c}{\partial v/\partial z+\partial w/\partial y}\end{array}\right)}\\ \multicolumn{1}{c}{+\frac{1}{2}\left(\begin{array}{c}{(\partial w/\partial x)}^2\hfill \\ \multicolumn{1}{c}{{(\partial w/\partial y)}^2}\\ \multicolumn{1}{c}{2(\partial w/\partial x)(\partial w/\partial y)}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right)={\mathbf{\epsilon }}^{0\mathrm{L}}+{\mathbf{\epsilon }}^{0\mathrm{NL}}}\end{array}$$

(17)

where ε^0L and ε^0NL are the linear and nonlinear strain terms, respectively.

By substituting equation (16) into equation (17), the in-plane ε can be obtained as follows:

$$\mathbf{\epsilon }={\left\{\begin{array}{ccc}{\epsilon }_{\mathrm{xx}}\hfill & {\epsilon }_{\mathrm{yy}}& {\gamma }_{\mathrm{xy}}\end{array}\right\}}^{\mathrm{T}}={\mathbf{\epsilon }}_{0\mathrm{L}}+z{\mathbf{\kappa }}_1+f(z){\mathbf{\kappa }}_2+{\mathbf{\epsilon }}_{0\mathrm{NL}}$$

(18)

where the in-plane forces, moments and shear strain components are expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\epsilon }}_{0\mathrm{L}}=\left\{\begin{array}{c}{u_0}_{,x}\hfill \\ \multicolumn{1}{c}{{v_0}_{,y}}\\ \multicolumn{1}{c}{{u_0}_{,y}+{v_0}_{,x}}\end{array}\right\},{\mathbf{\kappa }}_1=\left\{\begin{array}{c}{\theta }_{x,x}\hfill \\ \multicolumn{1}{c}{{\theta }_{y,y}}\\ \multicolumn{1}{c}{{\theta }_{x,y}+{\theta }_{y,x}}\end{array}\right\},}\\ \multicolumn{1}{c}{{\mathbf{\kappa }}_2=\left\{\begin{array}{c}{\beta }_{x,x}\hfill \\ \multicolumn{1}{c}{{\beta }_{y,y}}\\ \multicolumn{1}{c}{{\beta }_{x,y}+{\beta }_{y,x}}\end{array}\right\}}\end{array}$$

(19)

The nonlinear strain terms in equation (18) can be expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\epsilon }}_{0\mathrm{NL}}=\left\{\begin{array}{c}\frac{1}{2}{(w_0,x)}^2\hfill \\ \multicolumn{1}{c}{\frac{1}{2}{(w_0,y)}^2}\\ \multicolumn{1}{c}{w_0,y·w_0,x}\end{array}\right\}=\frac{1}{2}\left[\begin{array}{cc}{w}_0,x\hfill & 0\\ \multicolumn{1}{c}{0}& w_0,y\\ \multicolumn{1}{c}{w_0,y}& w_0,x\end{array}\right]\left[\begin{array}{cc}{w}_0,\hfill & x\\ \multicolumn{1}{c}{w_0,}& y\end{array}\right]=\frac{1}{2}{\mathbf{A}}_{\mathbf{\theta }}\mathbf{\theta }}\end{array}$$

(20)

Similarly, substituting equation (16) into equation (17), the transverse shear strain γ can be obtained as follows:

$$\mathbf{\gamma }={\left\{\begin{array}{c}{\gamma }_{\mathrm{xz}}{\gamma }_{\mathrm{yz}}\hfill \end{array}\right\}}^{\mathrm{T}}={\mathbf{\epsilon }}_s+f\prime (z){\mathbf{\kappa }}_s$$

(21)

where

$${\mathbf{\epsilon }}_s=\left\{\begin{array}{c}{w_0}_{,x}+{\theta }_x\hfill \\ \multicolumn{1}{c}{{w_0}_{,y}+{\theta }_y}\end{array}\right\},{\mathbf{\kappa }}_s=\left\{\begin{array}{c}{\beta }_x\hfill \\ \multicolumn{1}{c}{{\beta }_y}\end{array}\right\}$$

(22)

Constitutive relations

According to the elastic stress-strain constitutive relationship, the stress components of the kth layer can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{\left\{\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xy}}}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{yz}}}\end{array}\right\}}^{(k)}=\left[\begin{array}{ccccc}{Q}_{11}^{(k)}\hfill & Q_{12}^{(k)}& 0& 0& 0\\ \multicolumn{1}{c}{Q_{21}^{(k)}}& Q_{22}^{(k)}& 0& 0& 0\\ \multicolumn{1}{c}{0}& 0& Q_{66}^{(k)}& 0& 0\\ \multicolumn{1}{c}{0}& 0& 0& Q_{55}^{(k)}& 0\\ \multicolumn{1}{c}{0}& 0& 0& 0& Q_{44}^{(k)}\end{array}\right]{\left\{\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_y}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}}\end{array}\right\}}^{(k)}}\end{array}$$

(23)

The Eq. (23) can be rewritten as in a compact form:

$$\left\{\begin{array}{c}\mathbf{\sigma }\hfill \\ \multicolumn{1}{c}{\mathbf{\tau }}\end{array}\right\}=\left[\begin{array}{cc}{\mathbf{Q}}^{\mathrm{b}}\hfill & 0\\ \multicolumn{1}{c}{0}& {\mathbf{Q}}^{\mathrm{s}}\end{array}\right]\left\{\begin{array}{c}\mathbf{\epsilon }\hfill \\ \multicolumn{1}{c}{\mathbf{\gamma }}\end{array}\right\}$$

(24)

The in-plane forces, moments and shear forces can be obtained by

$$\mathbf{N}={\left\{\begin{array}{ccc}{N}_x\hfill & N_y& N_{\mathrm{xy}}\end{array}\right\}}^{\mathrm{T}}=\sum _{k=1}^{N_L}{{\int }_{z_k}^{z_{k+1}}\left\{\begin{array}{ccc}{\sigma }_x^{(k)}\hfill & {\sigma }_x^{(k)}& {\tau }_{\mathrm{xy}}^{(k)}\end{array}\right\}}^{\mathrm{T}}\mathrm{d}z$$

(25)

$$\mathbf{M}={\left\{\begin{array}{ccc}{M}_x& M_y& M_{xy}\end{array}\right\}}^{\mathrm{T}}=\sum _{k=1}^{N_L}{{\int }_{z_k}^{z_{k+1}}\left\{\begin{array}{ccc}{\sigma }_x^{(k)}& {\sigma }_x^{(k)}& {\tau }_{xy}^{(k)}\end{array}\right\}}^{\mathrm{T}}z\mathrm{d}z$$

(26)

$$\mathbf{P}={\left\{\begin{array}{ccc}{P}_x\hfill & P_y& P_x\end{array}\right\}}^{\mathrm{T}}=\sum _{k=1}^{N_L}{{\int }_{z_k}^{z_{k+1}}\left\{\begin{array}{ccc}{\sigma }_x^{(k)}\hfill & {\sigma }_x^{(k)}& {\tau }_{\mathrm{xy}}^{(k)}\end{array}\right\}}^{\mathrm{T}}f(z)\mathrm{d}z$$

(27)

$${\mathbf{Q}}^{\mathrm{s}}={\left\{\begin{array}{cc}{Q}_x^{\mathrm{s}}\hfill & Q_y^{\mathrm{s}}\end{array}\right\}}^{\mathrm{T}}=\sum _{k=1}^{N_L}{{\int }_{z_k}^{z_{k+1}}\left\{\begin{array}{cc}{\tau }_{\mathrm{xz}}^{(k)}\hfill & {\tau }_{\mathrm{yz}}^{(k)}\end{array}\right\}}^{\mathrm{T}}\mathrm{d}z$$

(28)

$${\mathbf{R}}^{\mathrm{s}}={\left\{\begin{array}{cc}{R}_x^{\mathrm{s}}\hfill & R_y^{\mathrm{s}}\end{array}\right\}}^{\mathrm{T}}=\sum _{k=1}^{N_L}{{\int }_{z_k}^{z_{k+1}}\left\{\begin{array}{cc}{\tau }_{\mathrm{xz}}^{(k)}\hfill & {\tau }_{\mathrm{yz}}^{(k)}\end{array}\right\}}^{\mathrm{T}}f\prime (z)\mathrm{d}z$$

(29)

By inserting Eq. (23) into Eqs. (25)-(29), the stress resultants are obtained as

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}\mathbf{N}\hfill \\ \multicolumn{1}{c}{\mathbf{M}}\\ \multicolumn{1}{c}{\mathbf{P}}\\ \multicolumn{1}{c}{{\mathbf{Q}}^{\mathrm{s}}}\\ \multicolumn{1}{c}{{\mathbf{R}}^{\mathrm{s}}}\end{array}\right\}=\left[\begin{array}{ccccc}\mathbf{A}\hfill & \mathbf{B}& \mathbf{C}& 0& 0\\ \multicolumn{1}{c}{\mathbf{B}}& \mathbf{D}& \mathbf{E}& 0& 0\\ \multicolumn{1}{c}{\mathbf{C}}& \mathbf{E}& \mathbf{F}& 0& 0\\ \multicolumn{1}{c}{0}& 0& 0& {\mathbf{A}}^{\mathrm{s}}& {\mathbf{B}}^{\mathrm{s}}\\ \multicolumn{1}{c}{0}& 0& 0& {\mathbf{B}}^{\mathrm{s}}& {\mathbf{D}}^{\mathrm{s}}\end{array}\right]\left\{\begin{array}{c}{\mathbf{\epsilon }}_{0\mathrm{L}}+{\mathbf{\epsilon }}_{0\mathrm{NL}}\hfill \\ \multicolumn{1}{c}{{\mathbf{\kappa }}_1}\\ \multicolumn{1}{c}{{\mathbf{\kappa }}_2}\\ \multicolumn{1}{c}{{\mathbf{\epsilon }}_s}\\ \multicolumn{1}{c}{{\mathbf{\kappa }}_s}\end{array}\right\}}\\ \multicolumn{1}{c}{=\left[\begin{array}{cc}{\mathbf{D}}_{\mathrm{b}}\hfill & 0\\ \multicolumn{1}{c}{0}& {\mathbf{D}}_{\mathrm{s}}\end{array}\right]\hat{\mathbf{\epsilon }}}\end{array}$$

(30)

where

$$\begin{array}{c}\multicolumn{1}{c}{(A_{\mathrm{ij}},B_{\mathrm{ij}},C_{\mathrm{ij}},D_{\mathrm{ij}},E_{\mathrm{ij}},F_{\mathrm{ij}})=}\\ \multicolumn{1}{c}{\sum _{k=1}^{N_L}{\int }_{z_k}^{z_{k+1}}{Q}_{\mathrm{ij}}^{\mathrm{b}}\{1,z,f(z),z^2,\mathrm{zf}(z),f{(z)}^2\}\mathrm{d}z,i,j=1,2,6}\end{array}$$

(31)

$$(A_{\mathrm{ij}}^s,B_{\mathrm{ij}}^s,D_{\mathrm{ij}}^s)=\sum _{k=1}^{N_L}{\int }_{z_k}^{z_{k+1}}{Q}_{\mathrm{ij}}^{\mathrm{s}}\{1,f\prime (z),{(f\prime (z))}^2\}\mathrm{d}z,i,j=4,5$$

(32)

The total potential energy principle

The strain energy of the matrix cracked FG-GOEAM plate is expressed as

$$U_{\mathrm{p}}=\frac{1}{2}{\int }_{\Omega }({\mathbf{\epsilon }}^{\mathrm{T}}{\mathbf{D}}_{\mathrm{b}}\mathbf{\epsilon }+{\mathbf{\gamma }}^{\mathrm{T}}{\mathbf{D}}_s\mathbf{\gamma })\mathrm{d}\Omega$$

(33)

The strain energy contributed by the Winkler elastic medium is expressed as:

$$U_{\mathrm{e}}={\int }_{\Omega }{K}_w·w^2\mathrm{d}\Omega$$

(34)

The virtual work done by transverse loading is given as

$$W_{\mathrm{e}}=-{\int }_{\Omega }q(x,y)w\mathrm{d}\Omega$$

(35)

By summing equations (33)–(35), we can obtain the total energy functional for the matrix cracked FG-GOEAM plate as follows:

$$\Pi =U_{\mathrm{p}}+U_{\mathrm{e}}+W_{\mathrm{e}}$$

(36)

An improving moving kriging meshless method

Let’s consider a domain Ω in ${\mathbb{R}}^2$ with the boundary Γ, consisting of N nodes. Within this domain, there exist n arbitrary x nodes located in the support domain Ω _x , as illustrated in Figure 3.

Figure 3.

Circular support domain of nodes in 2D model.

The Moving Kriging interpolation function u ^h (x) can be defined as follows [59]:

$${\mathbf{u}}^h(\mathbf{x})=[{\mathbf{p}}^{\mathrm{T}}(\mathbf{x})\mathbf{A}+{\mathbf{r}}^{\mathrm{T}}(\mathbf{x})\mathbf{B}]\mathbf{u}(\mathbf{x})$$

(37)

where

$$\begin{array}{c}\multicolumn{1}{c}{\{\begin{array}{c}\mathbf{p}(\mathbf{x})={\left[\begin{array}{c}{p}_1(\mathbf{x})p_2(\mathbf{x})p_3(\mathbf{x})\cdots p_m(\mathbf{x})\hfill \end{array}\right]}^{\mathrm{T}},\mathbf{A}={\left({\mathbf{P}}^{\mathrm{T}}{\mathbf{R}}^{-1}\mathbf{P}\right)}^{-1}{\mathbf{P}}^{\mathrm{T}}{\mathbf{R}}^{-1}\hfill \\ \multicolumn{1}{c}{\mathbf{r}(\mathbf{x})={\left\{\begin{array}{c}R({\mathbf{x}}_1,\mathbf{x})R({\mathbf{x}}_2,\mathbf{x})\dots R({\mathbf{x}}_n,\mathbf{x})\hfill \end{array}\right\}}^{\mathrm{T}},\mathbf{B}={\mathbf{R}}^{-1}(\mathbf{I}-\mathbf{PA})}\end{array}}\end{array}$$

(38)

in which I is n -dimensional identity matrix. P is an n×m matrix that stores the values of the polynomial basis functions. The correlation matrix R is an n×n square matrix. These two matrixes can be expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{P}=\left[\begin{array}{ccc}{p}_1\left({\mathbf{x}}_1\right)\hfill & \cdots & p_m\left({\mathbf{x}}_1\right)\\ \multicolumn{1}{c}{\cdots }& \ddots & \cdots \\ \multicolumn{1}{c}{p_1\left({\mathbf{x}}_n\right)}& \cdots & p_m\left({\mathbf{x}}_n\right)\end{array}\right],}\\ \multicolumn{1}{c}{\mathbf{R}=\left[\begin{array}{ccc}R({\mathbf{x}}_1,{\mathbf{x}}_1)\hfill & \cdots & R({\mathbf{x}}_1,{\mathbf{x}}_n)\\ \multicolumn{1}{c}{\cdots }& \ddots & \cdots \\ \multicolumn{1}{c}{R({\mathbf{x}}_n,{\mathbf{x}}_1)}& \cdots & R({\mathbf{x}}_n,{\mathbf{x}}_n)\end{array}\right]}\end{array}$$

(39)

For the two-dimensional elasticity problem, the basis functions p^T(x) are expressed using the following quadratic expressions:

$${\mathbf{p}}^{\mathrm{T}}(\mathbf{x})=\left\{\begin{array}{cccccc}1\hfill & x& y& x^2& \mathrm{xy}& y^2\end{array}\right\},(m=6)$$

(40)

Currently, the Gaussian model is the most widely used correlation function R(x _I , x _J ), and its expression is given by:

$$R({\mathbf{x}}_I,{\mathbf{x}}_J)=\exp (-\theta r_{\mathrm{IJ}}^2)$$

(41)

where $r_{\mathrm{IJ}}=\Vert {\mathbf{x}}_I-{\mathbf{x}}_J\Vert ,\theta >0$ is parameter which is strongly dependent on the coordinates in the support domain. However, the quality of the MK method is heavily influenced by the correlation parameter θ, which can sometimes lead to instability in the numerical model. Determining the optimal value for θ is still a matter of debate and remains uncertain. To overcome this drawback, Thai et al.^61,65 introduced the use of fourth-order spline functions to represent the correlation function, which can be expressed as follows:

$$R({\mathbf{x}}_I,{\mathbf{x}}_J)=1-6{\left(\frac{\theta r_{\mathrm{IJ}}}{a_0}\right)}^2+8{\left(\frac{\theta r_{\mathrm{IJ}}}{a_0}\right)}^3-3{\left(\frac{\theta r_{\mathrm{IJ}}}{a_0}\right)}^4,0\le \frac{\theta r_{\mathrm{IJ}}}{a_0}\le 1$$

(42)

where θ=1 and a₀ denotes the maximum distance between a pair of nodes within the support domain Ω _x .

Rewrite equation (37) into another form as

$${\mathbf{u}}^h(\mathbf{x})=\sum _{I=1}^n{N}_I(\mathbf{x})u_I$$

(43)

where the MK shape function N_I is given as

$$N_I(\mathbf{x})=\sum _{j=1}^m{p}_j(\mathbf{x})A_{\mathrm{jI}}+\sum _{k=1}^n{r}_k(\mathbf{x})B_{\mathrm{kI}}$$

(44)

and the first-order derivatives of the MK shape functions are computed as follows:

$$\{\begin{array}{c}{N}_{I,x}(\mathbf{x})=\sum _{j=1}^m{p}_{j,x}(\mathbf{x})A_{\mathrm{jI}}+\sum _{k=1}^n{r}_{k,x}(\mathbf{x})B_{\mathrm{kI}}\hfill \\ \multicolumn{1}{c}{N_{I,y}(\mathbf{x})=\sum _{j=1}^m{p}_{j,y}(\mathbf{x})A_{\mathrm{jI}}+\sum _{k=1}^n{r}_{k,y}(\mathbf{x})B_{\mathrm{kI}}}\end{array}$$

(45)

The MK interpolation function satisfies the Kronecker delta property, allowing for the direct imposition of essential boundary conditions, similar to standard FEM.

$${\varphi }_I\left({\mathbf{x}}_J\right)={\delta }_{\mathrm{IJ}}=\{\begin{array}{cc}1\hfill & (I=J,I=1,2,\dots ,n)\\ \multicolumn{1}{c}{0}& (I\ne J,J=1,2,\dots ,n)\end{array}$$

(46)

In this study, the circular influence domain Ω _x is defined as follows:

$$d_m=\alpha d_c$$

(47)

where a is the scale factor. The value of 2.4 is adopted for parameter a in the numerical simulations, as suggested in the literature.^61,65d_c is a characteristic length that is taken to be equal to the nodal spacing with regular distribution nodes.

Discrete system equations

From Eq. (43), the approximations of displacements for the matrix cracked FG-GOEAM plates is obtained as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{u}}_0^h=\left(\begin{array}{c}{u}_0^h\hfill \\ \multicolumn{1}{c}{v_0^h}\\ \multicolumn{1}{c}{w_0^h}\\ \multicolumn{1}{c}{{\theta }_x^h}\\ \multicolumn{1}{c}{{\theta }_y^h}\\ \multicolumn{1}{c}{{\beta }_x^h}\\ \multicolumn{1}{c}{{\beta }_y^h}\end{array}\right)=\sum _{l=1}^n{N}_I(x)\left(\begin{array}{c}{u}_l\hfill \\ \multicolumn{1}{c}{v_I}\\ \multicolumn{1}{c}{w_l}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{xl}}}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{yI}}}\\ \multicolumn{1}{c}{{\beta }_{\mathrm{xI}}}\\ \multicolumn{1}{c}{{\beta }_{\mathrm{yl}}}\end{array}\right)=\sum _{I=1}^n{N}_I(x){\mathbf{u}}_I}\end{array}$$

(48)

Substituting equation (48) into equations (19)–(20), the linear and nonlinear strain components can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\epsilon }}_{0\mathrm{L}}=\sum _{I=1}^n{\mathbf{B}}_I^{\mathrm{m}}{\mathbf{u}}_I,{\mathbf{\kappa }}_1=\sum _{I=1}^n{\mathbf{B}}_{1I}^{\mathrm{b}}{\mathbf{u}}_I,{\mathbf{\kappa }}_2=\sum _{i=1}^n{\mathbf{B}}_{2I}^{\mathrm{b}}{\mathbf{u}}_I,}\\ \multicolumn{1}{c}{{\mathbf{\epsilon }}_{0\mathrm{NL}}=\frac{1}{2}\sum _{I=1}^n{\mathbf{A}}_{\mathbf{\theta }}(\mathbf{u})\overline{\mathbf{Y}}{\mathbf{u}}_I,}\end{array}$$

(49)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{B}}_I^{\mathrm{m}}=\left[\begin{array}{ccccccc}{N}_{I,x}\hfill & 0& 0& 0& 0& 0& 0\\ \multicolumn{1}{c}{0}& N_{I,y}& 0& 0& 0& 0& 0\\ \multicolumn{1}{c}{N_{I,y}}& N_{I,x}& 0& 0& 0& 0& 0\end{array}\right],}\\ \multicolumn{1}{c}{{\mathbf{B}}_{1I}^{\mathrm{b}}=\left[\begin{array}{ccccccc}0\hfill & 0& 0& N_{I,x}& 0& 0& 0\\ \multicolumn{1}{c}{0}& 0& 0& 0& N_{I,y}& 0& 0\\ \multicolumn{1}{c}{0}& 0& 0& N_{I,y}& N_{I,x}& 0& 0\end{array}\right]}\end{array}$$

(50)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{B}}_{2I}^{\mathrm{b}}=\left[\begin{array}{ccccccc}0\hfill & 0& 0& 0& 0& N_{I,x}& 0\\ \multicolumn{1}{c}{0}& 0& 0& 0& 0& 0& N_{I,y}\\ \multicolumn{1}{c}{0}& 0& 0& 0& 0& N_{I,y}& N_{I,x}\end{array}\right],}\\ \multicolumn{1}{c}{{\mathbf{A}}_{\mathbf{\theta }}(\mathbf{u})={\left[\begin{array}{ccccccc}{w_0}_{,x}\hfill & 0& {w_0}_{,y}& 0& 0& 0& 0\\ \multicolumn{1}{c}{0}& {w_0}_{,y}& {w_0}_{,x}& 0& 0& 0& 0\end{array}\right]}^{\mathrm{T}}}\end{array}$$

(51)

$$\overline{\mathbf{Y}}=\left[\begin{array}{ccccccc}0\hfill & 0& N_{I,x}& 0& 0& 0& 0\\ \multicolumn{1}{c}{0}& 0& N_{I,y}& 0& 0& 0& 0\end{array}\right]$$

(52)

Substituting equation (48) into equation (21), the shear strain components can be rewritten as

$${\mathbf{\epsilon }}_{\mathrm{s}}=\sum _{I=1}^n{\mathbf{B}}_{0I}^{\mathrm{s}}{\mathbf{u}}_I,{\mathbf{\kappa }}_{\mathrm{s}}=\sum _{I=1}^n{\mathbf{B}}_{1I}^{\mathrm{s}}{\mathbf{u}}_I$$

(53)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{B}}_{0I}^{\mathrm{s}}=\left[\begin{array}{ccccccc}0\hfill & 0& N_{I,x}& N_I& 0& 0& 0\\ \multicolumn{1}{c}{0}& 0& N_{I,y}& 0& N_I& 0& 0\end{array}\right],}\\ \multicolumn{1}{c}{{\mathbf{B}}_{1I}^{\mathrm{s}}=\left[\begin{array}{ccccccc}0\hfill & 0& 0& 0& 0& N_I& 0\\ \multicolumn{1}{c}{0}& 0& 0& 0& 0& 0& N_I\end{array}\right]}\end{array}$$

(54)

Substituting equation (48) into equation (36) and taking its variation, we can obtain the following nonlinear bending equilibrium equation

$$({\mathbf{K}}_L+{\mathbf{K}}_{\mathrm{NL}}+{\mathbf{K}}_e)\mathbf{u}=\mathbf{F}$$

(55)

where K_L and K_NL are the linear and nonlinear stiffness matrix, respectively. K_e is Winkler foundation stiffness matrix, respectively. F is load vector. The expressions for these matrixes and the load vector are calculated by:

$${\mathbf{K}}_{\mathrm{L}}={\int }_{\Omega }{\mathbf{B}}^{\mathrm{T}}{\mathbf{D}}_{\mathrm{b}}\mathbf{B}\mathrm{d}\Omega +{\int }_{\Omega }{\mathbf{S}}^{\mathrm{T}}{\mathbf{D}}_{\mathrm{s}}\mathbf{S}\mathrm{d}\Omega$$

(56)

$${\mathbf{K}}_{\mathrm{NL}}={\int }_{\Omega }(\frac{1}{2}{\mathbf{B}}_{\mathrm{L}}^{\mathrm{T}}\mathbf{D}{\mathbf{B}}_{\mathrm{L}}+{\mathbf{B}}_{\mathrm{NL}}^{\mathrm{T}}\mathbf{D}{\mathbf{B}}_{\mathrm{L}}+\frac{1}{2}{\mathbf{B}}_{\mathrm{NL}}^{\mathrm{T}}\mathbf{D}{\mathbf{B}}_{\mathrm{NL}})\mathrm{d}\Omega$$

(57)

$${\mathbf{K}}_{\mathrm{e}}={\int }_{\Omega }{\mathbf{B}}_{\mathrm{e}}^{\mathrm{T}}{K}_w{\mathbf{B}}_{\mathrm{e}}d\Omega$$

(58)

$$\mathbf{F}={\int }_{\Omega }q(x,y){\left[\begin{array}{ccccccc}0\hfill & 0& N_I& 0& 0& 0& 0\end{array}\right]}^{\mathrm{T}}\mathrm{d}\Omega$$

(59)

where

$$\mathbf{B}=\left[\begin{array}{c}{\left({\mathbf{B}}^{\mathrm{m}}\right)}^{\mathrm{T}}{\left({\mathbf{B}}_1^{\mathrm{b}}\right)}^{\mathrm{T}}{\left({\mathbf{B}}_2^{\mathrm{b}}\right)}^{\mathrm{T}}\hfill \end{array}\right],\mathbf{S}=\left[\begin{array}{c}{\left({\mathbf{B}}_0^{\mathrm{s}}\right)}^{\mathrm{T}}{\left({\mathbf{B}}_1^{\mathrm{s}}\right)}^{\mathrm{T}}\hfill \end{array}\right],$$

(60)

$${\mathbf{D}}_{\mathrm{b}}=\left[\begin{array}{ccc}\mathbf{A}\hfill & \mathbf{B}& \mathbf{C}\\ \multicolumn{1}{c}{\mathbf{B}}& \mathbf{D}& \mathbf{E}\\ \multicolumn{1}{c}{\mathbf{C}}& \mathbf{E}& \mathbf{F}\end{array}\right],{\mathbf{D}}_{\mathrm{s}}=\left[\begin{array}{cc}{\mathbf{A}}^{\mathrm{s}}\hfill & {\mathbf{B}}^{\mathrm{s}}\\ \multicolumn{1}{c}{{\mathbf{B}}^{\mathrm{s}}}& {\mathbf{D}}^{\mathrm{s}}\end{array}\right],\mathbf{D}=\left[\begin{array}{cc}{\mathbf{D}}_{\mathrm{b}}\hfill & 0\\ \multicolumn{1}{c}{0}& {\mathbf{D}}_{\mathrm{s}}\end{array}\right]$$

(61)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{B}}_L={\left[\begin{array}{cc}\mathbf{B}\hfill & \mathbf{S}\end{array}\right]}^{\mathrm{T}},{\mathbf{B}}_{\mathrm{NL}}={\left[\begin{array}{cc}{\mathbf{B}}^{\mathrm{mNL}}\hfill & 0\end{array}\right]}^{\mathrm{T}},}\\ \multicolumn{1}{c}{{\mathbf{B}}_{\mathrm{e}}=\left[\begin{array}{ccccccc}0\hfill & 0& N_I& 0& 0& 0& 0\end{array}\right]}\end{array}$$

(62)

Solution procedure

The present study focuses on plate structures subjected to transverse loads, which typically do not exhibit snap-back and snap-through behavior. Therefore, for nonlinear static analysis, the Newton-Raphson iterative method is employed to solve the nonlinear equilibrium equations

$$({\mathbf{K}}_{\mathrm{L}}+{\mathbf{K}}_{\mathrm{e}}+{\mathbf{K}}_{\mathrm{NL}}(\mathbf{u}))\mathbf{u}=\mathbf{F}$$

(63)

The discretized Newton–Raphson scheme is formulated as:

$$\mathbf{R}(\mathbf{u})=({\mathbf{K}}_{\mathrm{L}}+{\mathbf{K}}_{\mathrm{e}}+{\mathbf{K}}_{\mathbf{NL}}(\mathbf{u}))\mathbf{u}-\mathbf{F}$$

(64)

$${}^{i+1}\mathbf{u}={}^i\mathbf{u}+\mathbf{\Delta }\mathbf{u}$$

(65)

where R(u) represents the residual force and ⁱ u represents the displacement at the ith iteration. The increment in displacement, denoted as Δu, is expressed as follows:

$$\mathbf{\Delta }\mathbf{u}=\frac{\mathbf{F}-{({\mathbf{K}}_{\mathrm{L}}+{\mathbf{K}}_{\mathrm{e}}+{\mathbf{K}}_{\mathrm{NL}}\left({}^i\mathbf{u}\right))}^i\mathbf{u}}{{\mathbf{K}}_{\mathrm{T}}}$$

(66)

where K_T represents the tangent stiffness matrix, which is mathematically defined as:

$${\mathbf{K}}_{\mathrm{T}}=\frac{\partial \mathbf{R}\left({}^i\mathbf{u}\right)}{{\partial }^i\mathbf{u}}={\mathbf{K}}_{\mathrm{L}}+{\mathbf{K}}_{\mathrm{e}}+{\mathbf{K}}_{\mathrm{NL}}+{\mathbf{K}}_{\mathrm{g}}$$

(67)

where

$${\mathbf{K}}_{\mathrm{NL}}={\int }_{\Omega }({\mathbf{B}}_{\mathrm{L}}^{\mathrm{T}}\mathbf{D}{\mathbf{B}}_{\mathrm{L}}+{\mathbf{B}}_{\mathrm{NL}}^{\mathrm{T}}\mathbf{D}{\mathbf{B}}_{\mathrm{L}}+{\mathbf{B}}_{\mathrm{NL}}^{\mathrm{T}}\mathbf{D}{\mathbf{B}}_{\mathrm{NL}})\mathrm{d}\Omega$$

(68)

In Eq. (67), K_g represents the geometric stiffness matrix, which can be computed once the displacement is determined at the previous iteration. The expression for K_g is as follows:

$${\mathbf{K}}_{\mathrm{g}}={\int }_{\Omega }{\overline{\mathbf{Y}}}^{\mathrm{T}}\left[\begin{array}{cc}{N}_x\hfill & N_{\mathrm{xy}}\\ \multicolumn{1}{c}{N_{\mathrm{xy}}}& N_y\end{array}\right]\overline{\mathbf{Y}}\mathrm{d}\Omega$$

(69)

The process is repeated until the convergence criterion is met. hence, a relative displacement norm criterion is given in the following form:

$$\frac{\Vert {}^{i+1}\mathbf{u}-{}^i\mathbf{u}\Vert }{\Vert {}^i\mathbf{u}\Vert }\le \chi$$

(50)

where χ represents the predefined tolerance which is set to 10^-2 in the present study.

Numerical results and discussion

Firstly, this section presents the nonlinear bending analysis of isotropic, orthotropic, laminated and FG-GPLRC plates for the convergence and effectiveness. Then, detailed parametric studies are carried out to present the effect of cracked density, foundation coefficients, distribution patterns, GOri contents, folding degree, thickness-to-width ratios, aspect ratios and boundary conditions on the non-dimensional deflection and stress of the FG-GOEMA plates. The numerical integration is performed by the background mesh, with a Gauss quadrature scheme applied using a 4 × 4 point configuration in each mesh element.

Three kinds of boundary conditions are considered as follows:

SSSS-1: Movable simply supported

$$\{\begin{array}{c}{v}_0=w_0={\theta }_y={\beta }_y=0\mathrm{at}x=0,a\hfill \\ \multicolumn{1}{c}{u_0=w_0={\theta }_x={\beta }_x=0\mathrm{at}y=0,b}\end{array}$$

(75)

SSSS-2: immovable simply supported

$$\{\begin{array}{c}{u}_0=v_0=w_0={\theta }_y={\beta }_y=0\mathrm{at}x=0,a\hfill \\ \multicolumn{1}{c}{u_0=v_0=w_0={\theta }_x={\beta }_x=0\mathrm{at}y=0,b}\end{array}$$

(76)

CCCC: clamped

$$\begin{array}{c}\multicolumn{1}{c}{u_0=v_0=w_0={\theta }_x={\theta }_y={\beta }_x={\beta }_y=0\mathrm{at}x=0,a\mathrm{and}y=0,b}\end{array}$$

(77)

In this paper, unless otherwise stated, the non-dimensional deflection, thickness, load factor, stress and Winkler parameter are defined as follows:

$$\{\begin{array}{c}\multicolumn{1}{c}{\overline{w}=w(\frac{a}{2},\frac{b}{2})/h,\overline{z}=\frac{z}{h},\overline{q}=\frac{q_0{a}^4}{E_{*}{h}^4},{\overline{\sigma }}_x=\frac{h^2}{q_0{a}^2}{\sigma }_x(\frac{a}{2},\frac{b}{2}),}\\ \multicolumn{1}{c}{{\overline{\tau }}_{\mathrm{yz}}=\frac{h}{q_{0}b}{\tau }_{\mathrm{yz}}(\frac{a}{2},0),k_w=\frac{K_w{b}^4}{D},D=\frac{E_{\mathrm{Cu}}{h}^3}{12(1-{\nu }_{\mathrm{Cu}}^2)}}\end{array}$$

(78)

where the letter E_* can represent E for isotropic plates, E₂ for laminated plates, and E_Cu for FG-GOEAM plates.

Nonlinear bending analysis of isotropic and orthotropic plates

Firstly, a convergence study is conducted for the nonlinear bending analysis of a clamped (CCCC) thin plate subjected to uniform distributed transverse load. The material and geometric parameters of the plate are: E = 206.84 GPa, ν = 0.316, a/b = 1, b/h = 100. Five node discretization schemes of 15 × 15, 17 × 17, 19 × 19, 21 × 21, and 23 × 23 are considered. The non-dimensional load-central deflections are illustrated in Figure 4. It can be observed that the results obtained using the proposed method with 21 × 21 nodes agree well with the analysis results of Tran et al. ⁶⁶ based on the five-variable IGA-TSDT. Therefore, a discretization with 21 × 21 nodes is employed for all further analysis. Figure 5 shows the relationship between normal in-plane stress ${\overline{\sigma }}_x={\sigma }_x(h^2/q_0{a}^2)(a/2,b/2,h/2)$ and load parameters $\overline{q}=q_0{a}^4/E{h}^4$ , and is also observed that the present solution matches well with IGA-TSDT results of Tran et al. ⁶⁶ Moreover, Figure 6 illustrates the computational efficiency assessed through computational time, comparing the time consumed by the same IMK used with conventional FSDT with five variables. The programs are executed in the same software environment using MATLAB 2021, with hardware consisting of an Intel Core i7-14700F CPU at 2.10 GHz and 32GB RAM. It can be observed that the computational time for IMK based on C⁰-HSDT is higher than that for IMK based on FSDT, and the difference in computational time between the two theoretical frameworks increases with the number of nodes.

Figure 4.

Convergence analysis of load-central deflection-curve for the clamped isotropic thin plate.

Figure 5.

Comparison of the normal in-plane stress for the clamped isotropic thin plate.

Figure 6.

Comparison of the computational efficiency of the present method with the IMK based on FSDT.

Additionally, we carried out a nonlinear bending analysis of a simply supported (SSSS1) orthotropic thin plates under uniformly distributed transverse load. The plate length a = b = 12in, thickness h = 0.138in, and material properties E₁ = 3 × 106psi, E₂ = 1.28 ×106psi, G₁₂ = G₁₃ = G₂₃ = 0.37 × 106psi, ν₁₂ = 0.32. The central deflection-load curve obtained from the proposed method, together with results obtained by way of experiment, ⁶⁷ Navier-CPT, ⁶⁷ FEM-FSDT, ⁶⁸ and IGA-TSDT, ⁶⁶ are shown in Figure 7. From Figure 7, it can be observed that the present result is agreement with to the results obtained using IGA-TSDT and IGA-FSDT. Additionally, compared to the result of Navier-CPT, the numerical solutions obtained from the IGA-FSDT and IGA-TSDT are closer to the experimental results. This can be attributed to the fact that FSDT and TSDT account for shear effects.

Figure 7.

Comparison of the central deflection for the simply supported orthotropic thin plates.

Linear and nonlinear bending analysis of laminated plates

Further, the method will be applied to analyze the linear and nonlinear bending analysis of a simply supported (SSSS1) square laminated (0°/90°/90°/0°) plate under sinusoidal distributed load (SDL) q = q₀sin(πx/a)sin(πx/b). The material properties are: E₁/E₂ = 25, G₁₂/E₂ = G₁₃/E₂ = 0.5, G₂₃/E2 = 0.2, ν₁₂ =0.25. The non-dimensional central linear deflection and stresses for the various span-to-thickness ratio a/h are calculated by the present method are tabulated in Table 2. Specifically, when a/h = 4, this paper selected the shear functions TSDT, SSDT, and ESDT from Table 1, and combined them with IMK interpolation to provide numerical results, which are also listed in Table 2. The non-dimensional quantities used in Table 2 are as follow:

$$\{\begin{array}{c}\multicolumn{1}{c}{\overline{w}*=\frac{100E_2{h}^3}{q_0{b}^4}w(\frac{a}{2},\frac{b}{2},0),{\overline{\sigma }}_x=\frac{h^2}{q_0{b}^2}{\sigma }_x(\frac{a}{2},\frac{b}{2},\frac{h}{2}),}\\ \multicolumn{1}{c}{{\overline{\sigma }}_y=\frac{h^2}{q_0{b}^2}{\sigma }_y(\frac{a}{2},\frac{b}{2},\frac{h}{4}),{\overline{\sigma }}_{\mathrm{xy}}=\frac{h^2}{q_0{b}^2}{\sigma }_{\mathrm{xy}}(0,0,\frac{h}{2}),}\\ \multicolumn{1}{c}{{\overline{\tau }}_{\mathrm{yz}}=\frac{h}{q_{0}b}{\tau }_{\mathrm{yz}}(\frac{a}{2},0,0),{\overline{\tau }}_{\mathrm{xz}}=\frac{h}{q_{0}b}{\tau }_{\mathrm{xz}}(0,\frac{b}{2},0)}\end{array}$$

(79)

Table 2.

Non-dimensional central linear deflection ${\overline{w}}^{*}$ and stresses for the simply supported square laminated (0°/90°/90°/0°) plate under sinusoidal distributed load (SDL).

a/h	Methods	${\overline{w}}^{*}$	$x$	$y$	$\mathrm{xy}$	$\mathrm{xz}$	$\mathrm{yz}$
4	Exact 3D 73	1.9540	0.7200	0.6660	0.0467	0.2700	0.2910
	Navier-TSDT 46	1.8937	0.6651	0.6322	0.0440	0.2064	0.2390
	IMK-TSDT	1.9021	0.7054	0.6306	0.0459	0.2137	0.2459
	Navier-SSDT 47	1.9088	0.6830	0.6349	0.0450	0.2162	0.2462
	IMK-SSDT	1.9126	0.7069	0.6335	0.0467	0.2212	0.2519
	Navier-ESDT 49	1.9193	0.7004	0.6367	0.0459	0.2264	0.2532
	IMK-ESDT	1.9203	0.7063	0.6357	0.0471	0.2281	0.2573
	Navier-IHSDT 48	1.9257	0.7255	0.6390	0.0473	0.2500	0.2700
	Present	1.9288	0.6992	0.6398	0.0466	0.2444	0.2705
10	Exact 3D 73	0.7430	0.5590	0.4010	0.0275	0.3010	0.1960
	Navier-TSDT 46	0.7147	0.5456	0.3888	0.0268	0.2640	0.1530
	Navier-SSDT 47	0.7198	0.5490	0.3910	0.0270	0.2788	0.1582
	Navier-ESDT 49	0.7240	0.5525	0.3929	0.0272	0.2941	0.1634
	Navier-IHSDT 48	0.7284	0.5578	0.3947	0.0275	0.3287	0.1760
	Present	0.7286	0.5602	0.3946	0.0274	0.3331	0.1779
20	Exact 3D 73	0.5170	0.5430	0.3080	0.0230	0.3280	0.1560
	Navier-TSDT 46	0.5060	0.5393	0.3043	0.0228	0.2825	0.1230
	Navier-ESDT 49	0.5088	0.5411	0.3057	0.0229	0.3160	0.1311
	Navier-IHSDT 48	0.5102	0.5425	0.3064	0.0230	0.3542	0.1412
	Present	0.5101	0.5432	0.3062	0.0229	0.3613	0.1422
100	Exact 3D 73	0.4390	0.5390	0.2760	0.0220	0.3370	0.1410
	Navier-TSDT 46	0.4343	0.5387	0.2708	0.0213	0.2897	0.1390
	Navier-ESDT 49	0.4350	0.5380	0.2700	0.0210	0.3240	0.1180
	Navier-IHSDT 48	0.4345	0.5388	0.2710	0.0214	0.3643	0.1271
	Present	0.4340	0.5374	0.2698	0.0212	0.3711	0.1269

It can be observed from the Table 2 that the present results are closer to their respective Navier-IHSDT solutions. Compared to the solution provided by the Navier-TSDT, ⁴⁶ Navier-SSDT, ⁴⁷ Navier-ESDT ⁴⁹ and Navier-IHSDT, ⁴⁸ the relative change occurs in the in-plane stress ${\overline{\sigma }}_x$ are 6.05%, 3.5%, 0.8% and 3.6%, respectively. Additionally, it can be observed that the methods utilizing the inverse hyperbolic shear function exhibit a higher level of agreement with the Exact 3D solution, particularly for the transverse shear stress.

Table 3 lists the comparison between linear and nonlinear bending solutions of the above laminated plate with various span-to-thickness ratios. It is evident that the present results show good agreement with the results based on the IGA-TSDT ⁵⁸ method. Figure 8 shows the stress distribution through the thickness of the plate (a/h = 10) via the increases of load parameter. As expected, the nonlinear effects caused by strain-displacement become more dominant under higher loads. This means that for the same normalized external load, lower normalized stresses are developed. It is worth noting that as the load increases, asymmetry in the through-thickness distribution predominantly affects the in-plane stress ${\overline{\sigma }}_x$ , while the transverse shear stress ${\overline{\tau }}_{\mathrm{yz}}$ maintains their symmetry. This observation is consistent with the findings reported by Tran et al. ⁶⁶

Table 3.

Non-dimensional central deflection $\overline{w}$ of the linear and nonlinear bending for the simply supported square laminated (0°/90°/90°/0°) plate under sinusoidal distributed load (SDL).

a/h	$\overline{q}$	Linear				Nonlinear
		FEM- FSDT 74	IGA- TSDT 66	Exact- 3D 73	Present	FEM- FSDT 74	IGA- TSDT 66	Present
4	50	0.856	0.947	0.977	0.964	0.679	0.720	0.727
	100	1.712	1.894	/	1.929	1.079	1.121	1.133
	200	3.423	3.787	/	3.858	1.611	1.656	1.677
	300	5.135	5.681	/	5.787	1.988	2.045	2.065
10	50	0.331	0.357	0.372	0.364	0.324	0.347	0.354
	100	0.662	0.715	/	0.729	0.612	0.650	0.660
	200	1.324	1.430	/	1.457	1.067	1.115	1.127
	300	1.986	2.144	/	2.186	1.410	1.461	1.473
20	50	0.245	0.253	0.259	0.255	0.243	0.250	0.253
	100	0.490	0.506	/	0.510	0.473	0.487	0.491
	200	0.980	1.012	/	1.020	0.876	0.896	0.901
	300	1.470	1.518	/	1.531	1.202	1.226	1.231
100	50	0.216	0.217	0.217	0.217	0.215	0.216	0.216
	100	0.432	0.434	/	0.434	0.423	0.424	0.424
	200	0.865	0.868	/	0.868	0.797	0.799	0.798
	300	1.297	1.303	/	1.302	1.112	1.115	1.114

Figure 8.

Influence of load parameters $\overline{q}$ on the through thickness distribution of stresses for the simply supported square laminated (0°/90°/90°/0°) plate ( $\overline{w}$ , a/h = 10).

Nonlinear bending analysis of FG-GPLRC plates

In this example, a clamped (CCCC) FG-GPLRC square plate subjected to uniformly distributed load (UDL) is considered. The geometry of the plate and GPLs are: a×b×h = 0.45 × 0.45 × 0.045 m, l_GPL ×w_GPL × h_GPL = 2.5 ×1.5 × 1.5 nm. The material properties are given as: E_GPL = 1.01TPa, E_m =3.0 GPa, ν_GPL = 0.186, ν_m = 0.34. The total number of layers of FG-GPLRC plate is set to be N_L = 10 and the non-dimensional UDL is defined as $\overline{q}=q_0{a}^2(1-{\nu }_{\mathrm{m}}^2)/E\mathrm{m}{h}^2$ . Figures 9–10 present the comparison of the non-dimensional linear deflections of FG-GPLRC plates under various GPL distribution patterns and weight fractions obtained in this study with previous results from available Gholami and Ansari ⁶⁹ , respectively. It is observed that the present results are in good agreement with that reported by Gholami and Ansari ⁶⁹ using the SSDT and generalized differential quadrature (GDQ).

Figure 9.

Comparison of the central deflection for the clamped FG-GPLRC plates (g_GPL = 0.8%).

Figure 10.

Comparison of the central deflection for the clamped X-GPLRC plates.

Nonlinear bending analysis of FG-GOEAM plates

In this section, the tunable nonlinear bending response of matrix cracked FG-GOEAM plates is investigated with the focus on the effects of crack density, distribution patterns, GOri contents, folding degree and foundation coefficients on the nonlinear bending deflection and normal stress for the FG metamaterial plates with various boundary condition. The Cu is utilized as the metal matrix for the metamaterials, with Young’s modulus E_Cu = 65.79GPa, Poisson’s ratio ν_Cu = 0.387 and mass density ρ_Cu = 8.8 g/cm³ at 300 K. The GOri filler with Young’s modulus E_Gr = 929.57 GPa, Poisson’s ratio ν_Gr = 0.22 and mass density ρ_Gr = 1.8 g/cm³ at 300 K. Additionally, the pristine graphene’s dimensions are specified as a length l_Gr = 83.76 Å and a thickness of t_Gr = 3.4 Å. Unless otherwise specified, the FG-GOEAM plate is subjected to uniformly distributed loads (UDL) under movable simply supported (SSSS1) conditions. The plate aspect ratio is 1, the width to thickness ratio is 10, the GOri contents is 2.5%, and the GOri folding degree is 80%.

Figure 11(a) and (b) show the effect of GOri distribution patterns on the nonlinear bending and in-plane through thickness of FG-GOEAM plates, respectively. For comparison, the results of geometric linearity are also plotted in the Figure 12. Four types of GOri distribution patterns, that is, UD-W_Gr, X-W_Gr, O-W_Gr, and V-W_Gr are considered. It can be seen that the X-W_Gr plate has the lowest, while the V-W_Gr plate has the highest deflections among the four cases in Figures 11(a) and 12(a). It can be also seen from the Figures 11(a) and 12(a) that the load-central deflection for the plate of UD-W_Gr and O-W_Gr types are very close to each other when the load factor is sufficiently smaller. Additionally, Figures 11(b) and 12(b) reveals that in the case of linear bending, the in-plane stress distribution through the thickness of the UD-W_Gr, X-W_Gr, and O-W_Gr plates, which exhibits a symmetric distribution of GOri content, is symmetric about the z = 0. However, for the nonlinear bending, the in-plane stress deviates from central symmetry about z = 0 due to the mid-plane tensile effects.

Figure 11.

Effect of Gori distribution patterns on the nonlinear bending response of FG-GOEAM plates: (a) load-central deflection; (b) normal in-plane stress versus thickness.

Figure 12.

Effect of Gori distribution patterns on the linear bending response of FG-GOEAM plates: (a) load-central deflection; (b) normal in-plane stress versus thickness.

Figures 13(a) and (b) present the effect of crack density on the nonlinear bending and in-plane through thickness of X-W_Gr plate, respectively. The results indicate that X-W_Gr plate containing matrix cracks has higher deflection and in-plane stress than the intact one, and increase with the increase of crack density.

Figure 13.

Effect of crack density on the nonlinear bending response of X-W_Gr plates: (a) load-central deflection; (b) normal in-plane stress versus thickness.

Figure 14 illustrates the nonlinear bending responses of matrix cracked FG-GOEAM plates tuned by GOri content. The significant enhancement effect of GOri content can be observed in Figure 14(a), where the deflections decrease significantly with increasing GOri content. This can be attributed to the fact that incorporating more GOri into the metal matrix increases the plate stiffness, resulting in lower deflections. In Figure 14(b), an intriguing observation is found that the deflection of FG-GOEAM plates do not exhibit significant variations when the GOri content is below 1.25 wt%. Conversely, the deflection experiences a dramatic decline as the GOri content increases beyond this threshold. This phenomenon can be attributed to the transformation of the GOri/Cu nanocomposite into the GOEAM when the GOri weight fraction reaches 1.25 wt%, as reported by Zhao et al.^24,26 This transformation results in an increased NPR and Young’s modulus, thus enhancing the plate stiffness.

Figure 14.

Nonlinear bending responses of matrix cracked FG-GOEAM plates tuned by GOri content: (a) load-central deflection; (b) central deflection versus GOri content; (c) normal in-plane stress versus GOri content; (d) normal shear stress versus GOri content.

The normal in-plane stress ${\overline{\sigma }}_x=\frac{h^2}{q_0{b}^2}{\sigma }_x$ (a/2, b/2, h/2) and shear stress ${\overline{\tau }}_{\mathrm{yz}}=\frac{h}{q_{0}b}{\tau }_{\mathrm{yz}}$ (a/2, 0, 0) of the matrix cracked FG-GOEAM plate can be also tuned by GOri content, as shown in Figure 14(c) and (d), respectively. The in-plane stress variation in FG-GOEAM plates with respect to the GOri content is observed to be dependent on the density parameters β. For β = 0, the in-plane stress decreases as the GOri content increases. Conversely, for β = 1.0, the in-plane stress increases with increasing GOri content. Notably, when the β = 0.5, the in-plane stress remains relatively stable as the GOri content varies. The shear stress of FG-GOEAM plate is observed to decrease as the density parameter increases. Moreover, the shear stress of UD-W_Gr plate shows a slight increase with increasing GOri content, whereas the shear stress of X-W_Gr plate decreases with increasing GOri content.

Figure 15 depicts the nonlinear bending responses of matrix cracked FG-GOEAM plates tuned by GOri folding degree. According to the reported by Zhao et al.^24,26, when the folding degree of GOri/Cu nanocomposites reaches 50%, they undergo a transition to GOEAM due to the appearance of NPR. It can be observed from the Figure 15(a) and (b) that the deflection of FG-GOEAM plate increases with increase of the GOri folding degree. It can be seen from the Figure 15 and d that in-plane stress ${\overline{\sigma }}_x=\frac{h^2}{q_0{b}^2}{\sigma }_x$ (a/2, b/2, h/2) and shear stress ${\overline{\tau }}_{\mathrm{yz}}=\frac{h}{q_{0}b}{\tau }_{\mathrm{yz}}$ (a/2, 0, 0) of the matrix cracked FG-GOEAM plate overall decrease with the increase of the GOri folding degree, except for a few points.

Figure 15.

Nonlinear bending responses of matrix cracked FG-GOEAM plates tuned by GOri folding degree: (a) load-central deflection; (b) central deflection versus GOri folding degree; (c) normal in-plane stress versus GOri folding degree; (d) normal shear stress versus GOri folding degree.

Figure 16(a), (b), (c), and (d) present various boundary conditions, foundation coefficients, load form and width-to-thickness ratios on the nonlinear bending response of FG-GOEAM plates, respectively. As expected, the deflection of the FG-GOEAM plate decreases with the increase of boundary constraints, foundation coefficient, and width-to-thickness ratio. The load-deflection curves of the FG-GOEAM plate subject to a sinusoidal distributed load (SDL) are lower than those of the plate subjected to uniformly distributed load (UDL). The above phenomenon has also been validated in the nonlinear bending analysis of functionally graded graphene-reinforced composite (FG-GRC) laminated plates and functionally graded carbon nanotube-reinforced (FG-CNTRC) composites plates in Shen et al.’s^70,71 literature.

Figure 16.

Effect of other parameters on the nonlinear bending response of FG-GOEAM plates: (a) boundary condition; (b) load form; (c) foundation coefficients; (d) width-to-thickness ratios.

Conclusions

To assess the effect of crack density, distribution pattern, GOri content, folding degree, foundation coefficients and boundary conditions on the tunability of the nonlinear bending response of matrix cracked FG-GOEAM plates, a meshless model in the basis of Improved Moving Kriging (IMK) interpolation and theoretical framework of C⁰-HSDT with seven variables is proposed for it. The meshless model is suitable for analyzing thin, moderately thick, and thick plates, and it circumvents the challenges related to constructing high-order shape functions and applying clamped boundary conditions. The GP-assisted and self-consistent micromechanical models are used to estimate the material properties of FG-GOEAM plate and the properties of matrix cracked layers, respectively. The von-Karman strain assumptions are used, and the meshless nonlinear governing equations are numerically solved by the Newton-Raphson iterative method. The numerical results indicate that:

• (1) The present method exhibits good convergence and demonstrates high computational accuracy for thin plates, moderately thick plates, and thick plates. However, compared to IMK based on traditional FSDT, the computational efficiency is reduced.

• (2) The rise of GOri content leads to an increase in the stiffness of the FG-GOEAM plate, resulting in a decrease in deflection, while the growth of GOri folding degree shows the opposite trend.

• (3) The deflection and in-plane stress of the FG-GOEAM plate with matrix cracks are higher than those of the intact plate, while the shear stress is exactly the opposite, and this phenomenon becomes more pronounced as the crack density increases.

• (4) The relationship between crack density and the curves of normal stress versus GOri folding degree and GOri content displays a consistent pattern, further demonstrating the tunable behavior of the FO-GOEAM plate.

• (5) The increase of the boundary conditions, foundation coefficients, and width-to-thickness ratio leads to the decrease of the deflection of the FG-GOEAMs plates.