Nonlocal theoretical inquiry into pultruded GFRP plate dynamics: Integrating experimental and numerical analyses

Abstract

Few studies combining experimental and theoretical analyses of pultruded fiber-reinforced polymer (FRP) composites are available in the literature. In this study, experimental methods and a nonlocal theoretical approach were used to investigate the flexural performance of pultruded glass-FRP (GFRP) plates. In the theoretical method used, the Eringen’s nonlocal elasticity theory was adopted so that small-scale material effects on pultruded FRP composite plates were taken into account. Based on the plate theory of high-order shear deformation, the kinematic relationships of the composite plates were established. The stiffness matrix coefficients and tangential stiffness matrix were shown along with the nonlocal force vector. Also, finite element modeling was performed with the help of ABAQUS for the accuracy of experimental and theoretical results. For the analysis, a variety of laminated plate types were taken into consideration. The flexural performance was investigated utilizing flexural tests under point and distributed loads. To comprehend the behavior of the pultruded GFRP plates more fully, macro mechanical damage analyses were carried out. Numerical values were compared with experimental findings, theoretical solution results, and finite element simulation results which were found to be very similar.

Introduction

Pultrusion-produced composite materials are used extensively all over the world. The use of profiles created by the pultrusion technique has gained more attention in recent decades from a variety of sectors.^1–3 Alternative materials such as glass fiber-reinforced polymer (GFRP) are frequently employed in manufacturing or construction for high-performance applications.^4–7 Fiber-reinforced polymers (FRPs) stand out as ideal materials for serving as structural load-carrying elements due to their exceptional advantages.^8–11 Their superior properties, including high strength and durability, make them well-suited for bearing structural loads effectively.^12–20 As a result, FRPs are widely utilized in various industries where robust and reliable structural components are essential.^21–25 Pultrusion is one of the most efficient of various composite manufacturing processes.^26–29 The benefits of pultrusion over other composite production techniques include its high production rate, ³⁰ better efficiency,^31,32 lower production costs, ³³ the ability to produce composite materials at a much faster pace, and the capacity to create profiles with almost infinite lengths. ³⁴ Matrix-based composites come in thermoplastic and thermoset varieties. ³⁵ Unidirectional fibers or fiber reinforcement in the form of fabric are drawn through an impregnation bath for thermoset pultrusion. ³⁶ The impregnated section is then passed through collimating plates to remove surplus resin before being placed into a heated forming die. After cooling, the component is ready for use after being impregnated with hot melt thermoplastic polymer with fiber reinforcing. The cutting saw is utilized to cut the formed profile to the right length once the polymerized composite has been pulled from the forming die by the pulling mechanism. ³⁷ Generally, pultruded components are formed of an embedded fiber discontinuous phase known as reinforcement and a continuous polymeric phase known as the matrix. Fillers can occasionally be used to boost the interaction between the matrix and reinforcement, reduce costs, or enhance the matrix’s properties.^38–42 Many researchers have documented or assessed the mechanical properties of the component elements to examine the mechanical behavior of pultruded pieces. ⁴³ Glass, carbon, or aramid fibers are typically used to make the pultruded FRP composites. ⁴⁴

In order to evaluate multilayered composite structures, many models were created as the usage of composite materials. The literature has evolved the Kirchhoff plate theory, first-order and higher-order shear deformation theories, also known as equivalent single layer theories, to assess isotropic and composite plates.^45–51 Various mathematical models have been extensively presented by Tornabene.^52–56 The Kirchhoff plate theory, sometimes referred to as the classical plate theory, is appropriate for thin plates since it disregards interlaminar shear deformation and rotational inertia. The Kirchhoff plate theory produces false findings for composites by ignoring transverse shear effects. Based on Reissner and Mindlin assumptions,^57,58 the first-order shear deformation plate theory is presented for the thick plates and includes the transverse shear effects with a linear change of transverse shear strain through the plate thickness. ⁵⁹ This theory is necessary as a shear correction to account for the variability of transverse shear stress and shear strain through thickness. Since it considers a linear transverse shear distribution across the plate thickness, the first-order shear deformation plate theory does not satisfy the requirement of zero transverse shear stress at the top and bottom surfaces of plates. Other factors that affect the shear correction factor are plate border circumstances, loading conditions, layer orientation, and geometrical features. ⁶⁰ More accurate models, including higher-order shear deformation theories, assume quadratic, cubic, and more variations surface parallel displacements throughout the whole length of the laminates to forecast how the structure would behave for thick to thin portions.^61–64

In contrast to experimental testing, modifications of conventional classical continuum theories, which mathematically estimate the size effects, may be a more affordable and trustworthy alternative technique of analysis. Classical theory in continuum mechanics assumes that stress at an arbitrary point is a function of its own stress, thus, it cannot capture the scale effect of advanced composites in microstructures. ⁶⁵ Recently, the accessible non-classical continuum models of nonlocal elasticity theory, modified couple stress theory, and strain gradient theory have been applied in the mechanical analysis of small-scale structures. ⁶⁶ To account for the scale effect in elasticity, the nonlocal elasticity theory has been proposed by Eringen. ⁶⁷ By include extra degrees of freedom in the kinematics assumptions for each material point, nonlocal elasticity models may be created.^68–74 For the axisymmetric bending analysis of circular plates under general stress, Duan and Wang ⁷⁵ deduced precise solutions of the nonlocal Kirchhoff plate theory model. The nonlocal parameter’s impacts on the circular graphene sheets deflection, radial moment, circumferential moment, and shear force were examined under uniform loads, with either clamped or just supported boundary conditions. Ma’en and Al-Kouz ⁷⁶ have focused into the free vibration analysis of nonlocal orthotropic Kirchhoff plates. Kirchhoff plates are described at the micro/nanoscale by the nonlocal elasticity theory of Eringen, which takes into account the small-scale effect. Both the first-order shear deformation plate theory and the nonlocal Kirchhoff plate theory models for graphene sheets were introduced by Kananipour, ⁷⁷ however, they were used in the static bending analysis with differing boundary conditions. For functionally graded circular plates, Hosseini-Hashemi et al. ⁷⁸ created a nonlocal first order shear deformation plate theory. Also, closed-form solutions were found for the natural frequencies of plates under various boundary conditions. Using von Karman nonlinear stresses, Golmakani and Rezatalab ⁷⁹ provided a nonlinear nonlocal first-order shear deformation plate theory model for the nonlinear bending analysis of orthotropic embedded plates. Based on the revised plate theory, Narendar ⁸⁰ suggested a nonlocal higher-order shear deformation theory model for the buckling analysis of isotropic nanoplates. Belkorissat et al. ⁸¹ developed a simple nonlocal higher-order shear deformation theory model for functionally graded nano-plates. Based on the quasi-3D sinusoidal theory with five unknowns, Bessaim et al. ⁸² developed a nonlocal quasi-3D model for the dynamic analysis of plates.

Several scholars utilize the usual concept that the highest deflection occurs in the mid-plan of the thickness. When the deflection field’s fluctuation across the thickness is taken into account, this assumption appears to be almost incorrect. The transverse deflections and thickness changes are confirmed to cause the highest deflections to occur at the plate’s top surface.

Analytical techniques as well as numerical approaches can be employed to solve the governing equations resulting from dimension dependent models. Nevertheless, owing to straightforward geometry, stress, and boundary conditions, the use of analytical techniques is restricted to a certain advanced composite. For instance, the Levy approach only works with rectangular plates that have two opposing sides that are simply supported, and the remaining two sides can have any arbitrary boundary condition, whereas the Navier method only works with rectangular plates that have simple supports on all four sides. It is impracticable to search for analytical solutions to real-world problems such as generic geometry, loading, and boundary conditions since size dependent models are more mathematically difficult than classical models. In order to tackle such problems, numerical approaches like the finite element method (FEM), Ritz and Galerkin methods, differential quadrature method, meshless method, etc., are known as some of the most efficient methods.^83,84 FEM is the most potent instrument among the many numerical approaches, and it is frequently employed for the study of structures.^85,86 Reddy ⁸⁷ created the nonlinear nonlocal Kirchhoff plate theory and first-order shear deformation plate theory models for the von Karman-based nonlinear bending analysis of nano-plates. For the creation of finite element solutions, the variational statement of these models was also supplied. Reddy and Chao ⁸⁸ developed a displacement FEM for the analysis of thick plates, based on first-order plate theory, describing the transverse shear. Fantuzzi et al. ⁸⁹ studied free vibration analysis of laminated composite plates using FEM and differential quadrature techniques based on Fourier broadening.

It is evident from the literature analysis that surface effects and small-scale effects are crucial in micro structural plates exposed to compressive in-plane loads, and as a result, the bending problem of laminated composite plates has to be examined in light of these effects. Throughout the past decade, substantial efforts have been given to the creation of size-dependent models based on higher-order continuum mechanics method. Currently, there is no combined research on pultruded GFRP plates in terms of experimental application, nonlocal theories, and numerical modeling and analysis. In this study, bending analyses of GFRP plates produced by pultrusion method were investigated experimentally, numerically, and theoretically. Since the pultrusion consists of GFRP plate veil layer, reinforcement phase, and matrix layer with various sands, finite element solution based on nonlocal elasticity theory was performed. It covers several layer orientations for rectangular plates with varying sizes. The Eringen’s nonlocal elasticity theory is utilized here. The weak forms are derived using the virtual displacement concept, and the displacement finite element models are created using the weak forms. In the polymer matrix of the modeled system, shear and transverse forces coexist. A rectangle conforming element with four nodes and eight degrees of freedom are employed. Together with the nonlocal force vector, the stiffness matrix coefficients and tangent stiffness matrix are demonstrated. Experimental, numerical, and theoretical results are presented comparatively. In addition, macro damage analysis is conducted.

Experimental study

This section uses test procedures on specifically manufactured pultruded GFRP composite plate specimens to examine the accuracy of theoretical and numerical solution techniques as well as the behavior of the material (Figure 1). The fibers used in the pultruded GFRP plates are E-glass. 4800 tex fiber bundles were utilized in the glass roving and 450 g/m² nonwoven fabric was employed for the glass mat. In compliance with standards, GFRP composite plates were used for the testing.

Figure 1.

Pultruded GFRP plate.

The elastic modulus and shear modulus of the square pultruded GFRP sections were determined previously by Madenci et al.^38,39 and are listed in Table 1.

Table 1.

Mechanical properties of pultruded GFRP.

Property	Mean value (MPa)
Longitudinal tensile modulus of elasticity	23,000
Transverse tensile modulus of elasticity	7000
Longitudinal tensile strength	240
Transverse tensile strength	50
Longitudinal compressive strength	150
Transverse compressive strength	70
Shear strength	25

Specimen design

The pultruded GFRP strip is divided into four groups, resulting in a total of 24 examples. Polyester resin and GFRPs were used to create these specimens. The design properties of the pultruded GFRP specimens are given in Table 2. The thickness of all specimens was 6 mm. The plate section was 130 mm × 130 mm.

Table 2.

Design of specimens.

Specimen	Fiber orientation	Load type	Boundary conditions
SFSF_D_0°	0°	Distributed
SFSF_D_90°	90°	Distributed
CFCF_D_0°	0°	Distributed
CFCF_D_90°	90°	Distributed
SFSF_P_0°	0°	Point
SFSF_P_90°	90°	Point
CFCF_P_0°	0°	Point
CFCF_P_90°	90°	Point

Burnout test

The ratios of matrix and fiber of the pultruded GFRP plates were determined using a burnout test on the specimens produced from the strip. The combustion test was carried out in accordance with the ASTM D-2584 standard. The pultruded GFRP specimens, which were made by slicing the specimens into pieces 5 cm wide, were cooked in the oven at 600 °C for an hour. The matrix fully burnt off at this period, leaving only the fibers. Once the specimens were taken out of the furnace, fiber weights were calculated. Following the burnout test, the glass and veil layers that were removed in accordance with the direction of the fibers are depicted in detail in Figure 2.

Figure 2.

Burnout test: (a) FRPs of plate, (b) weighing after burning test, and (c) oven.

Tensile test

The pultruded GFRP exhibits varying resistance to pressures coming from various directions because of its fibrous and orthotropic structure. To ascertain the tensile characteristics of the pultruded GFRP composite beams, tensile tests were carried out both parallel and perpendicular to the direction of the fiber. On the traction machine with the hydraulic loader, tensile force was delivered to the ready specimens at a steady pace. By gathering the essential information, it was possible to compute the mechanical properties of the pultruded GFRP. A universal testing equipment was utilized to test the samples with nominal length, breadth, and thickness measurements of 250 mm × 25 mm × 3 mm under tensile pressures in order to assess their tensile strength. The technique outlined in ASTM D3039 was used.^90,91 One mm was applied each minute by the forces (Figure 3).

Figure 3.

(a) Tensile and compressive tests setup, (b) tensile test specimen, and (c) compressive test specimen.

Compressive tests

The samples, which had dimensions of 150 mm × 25 mm × 3 mm, were put under axial compression to find out how strong they were (length-width-thickness). The distance between the noses was 40 mm. The compressive properties of stiff polymers were evaluated using the ASTM D695 standard test procedure. ⁹² One mm was applied each minute by the forces (Figure 3).

Flexural test setup

A comprehensive experimental study was performed to comprehend the behavior of the pultruded GFRP plates under different loading conditions. Three different parameters were selected. These are boundary conditions, type of loading, and fiber direction. For boundary conditions simple-free-simple-free (SFSF) and clamped-free-clamped-free (CFCF) were selected. In order to provide clamped boundary conditions, the pultruded GFRP plates were squeezed via steel plates using clamps. For loading, point loading (P) and distributed loading (D) were adopted. Two different directions were considered as 0° and 90°. The specimens were tested under servo controlled hydraulic piston with capacity of 50 t. The load and displacement values were recorded. The speed of experiments was 2 mm/min. The loading was terminated when excessive deformation or reduction in loading was detected. For each configuration, three repetitions were done to ensure accurate findings. Figure 4 illustrates the setup of the flexural test.

Figure 4.

Flexural test setup and configurations: (a) specimen SFSF_D, (b) specimen SFSF_P, (c) specimen CFCF_D, and (d) specimen CFCF_P.

Theoretical formulation

Nonlocal elasticity model

Eringen maintains that the stress at a point of comparison in a continuum system requires not only the strain at that point but also the strains at all the nodes. This presumption requires that the constitutive description has a length scale parameter. The strain tensor at every other place on an elastic body is what determines the stress tensor at point x, in accordance with nonlocal elasticity theory.^5–7 For homogeneous elastic solids, the nonlocal stress-tensor components ${\sigma }_{ij}$ at each point x in the solid can be expressed as:

$${\sigma }_{ij}(x)=\underset{\Omega }{\int }\alpha \left(|x^{\prime }-x|,\tau \right)t_{ij}(x^{\prime })d\Omega (x^{\prime })$$

(1)

where $t_{ij}(x^{\prime })$ are the parts of the local stress tensor that are accessible at position x which are associated to the strain tensor components ${\epsilon }_{ij}$ as:

$$t_{ij}=C_{ijkl}{\epsilon }_{kl}$$

(2)

Because the component problem’s integral form solution is complicated, the following simplified differential form equation can be utilized instead.

$$\begin{array}{l}\Im {\sigma }^{nl}(\mathbf{x})=\left(1-\mu {\nabla }^2\right){\sigma }^{nl}=T\\ \mu ={\tau }^2{l}^2\end{array}$$

(3)

where ${\nabla }^2$ is the Laplacian operator and $\Im =\left(1-\mu {\nabla }^2\right)$ indicates nonlocal differential.

As opposed to the linear algebraic equations between the stress resultants and strains in a local theory, the nonlocal theory results in differential relations involving the stress resultants and the strains. In equation (4), we present these relations for homogeneous isotropic beams under the assumption that the nonlocal behavior is negligible in the thickness direction. The nonlocal constitutive relationship can be written as:

$${\sigma }_{xx}-\mu \frac{{\partial }^2{\sigma }_{xx}}{\partial x^2}={\overline{Q}}_{ij}{\epsilon }_{xx}$$

(4)

Kinematic relations

Transverse shear stresses in laminated composite plates cannot be ignored when assessing structures since they have lower shear modulus than in-plane modulus. Transverse shear strains and shear stresses can vary in a parabolic manner and are not required to use a shear correction. According to Reddy, ⁹³ the displacement field has been enlarged to the third degree of the thickness coordinate.

$$\left\{\begin{array}{c}u(x,y,z)\\ v(x,y,z)\\ w(x,y,z)\end{array}\right\}=\left\{\begin{array}{c}{u}_0(x,y)\\ v_0(x,y)\\ w_0(x,y)\end{array}\right\}+z\left\{\begin{array}{c}{\varphi }_x\\ {\varphi }_y\\ 0\end{array}\right\}-\frac{4z^3}{3h^2}\left\{\begin{array}{c}{\varphi }_x+\frac{\partial w_0}{\partial x}\\ {\varphi }_y+\frac{\partial w_0}{\partial y}\\ 0\end{array}\right\}$$

(5)

where $u_0$ , $v_0$ , and $w_0$ represent the in-plane displacements of a point on the mid-plane of the plate in x, y, and z directions (Figure 5).

Figure 5.

Geometry of a laminated composite plate.

The linear transverse and in-plane shear stresses are provided by:

$$\begin{array}{l}\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_{xx}^{(0)}\\ {\epsilon }_{yy}^{(0)}\\ {\gamma }_{xy}^{(0)}\end{array}\right\}+z\left\{\begin{array}{c}{\epsilon }_{xx}^{(1)}\\ {\epsilon }_{yy}^{(1)}\\ {\gamma }_{xy}^{(1)}\end{array}\right\}+z^3\left\{\begin{array}{c}{\epsilon }_{xx}^{(3)}\\ {\epsilon }_{yy}^{(3)}\\ {\gamma }_{xy}^{(3)}\end{array}\right\}\\ \left\{\begin{array}{c}{\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right\}=\left\{\begin{array}{c}{\gamma }_{yz}^{(0)}\\ {\gamma }_{xz}^{(0)}\end{array}\right\}+z^2\left\{\begin{array}{c}{\gamma }_{yz}^{(2)}\\ {\gamma }_{xz}^{(2)}\end{array}\right\}\end{array}$$

(6)

where

$$\left\{\begin{array}{c}{\epsilon }_{xx}^{(0)}\\ {\epsilon }_{yy}^{(0)}\\ {\gamma }_{xy}^{(0)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u_0}{\partial x}\\ \frac{\partial v_0}{\partial y}\\ \frac{\partial u_0}{\partial y}+\frac{\partial v_0}{\partial x}\end{array}\right\},\left\{\begin{array}{c}{\epsilon }_{xx}^{(1)}\\ {\epsilon }_{yy}^{(1)}\\ {\gamma }_{xy}^{(1)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial {\varphi }_x}{\partial x}\\ \frac{\partial {\varphi }_y}{\partial y}\\ \frac{\partial {\varphi }_x}{\partial y}+\frac{\partial {\varphi }_y}{\partial x}\end{array}\right\}$$

(7)

$$\left\{\begin{array}{c}{\epsilon }_{xx}^{(3)}\\ {\epsilon }_{yy}^{(3)}\\ {\gamma }_{xy}^{(3)}\end{array}\right\}=-\frac{4}{3h^2}\left\{\begin{array}{c}\frac{\partial {\varphi }_x}{\partial x}+\frac{{\partial }^2{w}_0}{\partial x^2}\\ \frac{\partial {\varphi }_y}{\partial y}++\frac{{\partial }^2{w}_0}{\partial y^2}\\ \frac{\partial {\varphi }_x}{\partial y}+\frac{\partial {\varphi }_y}{\partial x}+2\frac{{\partial }^2{w}_0}{\partial x\partial y}\end{array}\right\}$$

(8)

$$\left\{\begin{array}{c}{\gamma }_{yz}^{(0)}\\ {\gamma }_{xz}^{(0)}\end{array}\right\}=\left\{\begin{array}{c}{\varphi }_y+\frac{\partial w_0}{\partial y}\\ {\varphi }_x+\frac{\partial w_0}{\partial x}\end{array}\right\},\left\{\begin{array}{c}{\gamma }_{yz}^{(2)}\\ {\gamma }_{xz}^{(2)}\end{array}\right\}=-\frac{4}{h^2}\left\{\begin{array}{c}{\varphi }_y+\frac{\partial w_0}{\partial y}\\ {\varphi }_x+\frac{\partial w_0}{\partial x}\end{array}\right\}$$

(9)

Stress-strain relations

The layer’s material behavior is described via the Eringen’s nonlocal elasticity model, which relates the local stresses ${\sigma }_{ij}$ to the nonlocal stresses ${\sigma }_{ij}^{nl}$ via the following differential relationship.

$$\Im {\sigma }_{ij}^{nl}={\sigma }_{ij}$$

(10)

Equation (10) may be used to obtain the relationships between the local and nonlocal stress resultants as:

$$\begin{array}{l}\Im N_{ij}^{nl}=N_{ij};\Im M_{ij}^{nl}=M_{ij};\Im P_{ij}^{nl}=P_{ij};\Im N_{ij}^{nl}=N_{ij}\\ \Im Q_i^{nl}=Q_i;\Im R_i^{nl}=R_i\end{array}$$

(11)

$$\left\{\begin{array}{c}{N}_{ij}\\ M_{ij}\\ P_{ij}\end{array}\right\}=\underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}{\sigma }_{ij}\left\{\begin{array}{c}1\\ z\\ z^3\end{array}\right\}dz;\left\{\begin{array}{c}{Q}_i\\ R_i\end{array}\right\}=\underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}{\sigma }_{iz}\left\{\begin{array}{c}1\\ z^2\end{array}\right\}dz$$

(12)

For the pultruded GFRP plate, the constitutive relations of nonlocal theory may be expressed as:

$$\left(1-\mu {\nabla }^2\right)\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\tau }_{xy}\\ {\tau }_{yz}\\ {\tau }_{xz}\end{array}\right\}=\left[\begin{array}{ccccc}{\overline{Q}}_{11}& {\overline{Q}}_{12}& {\overline{Q}}_{16}& 0& 0\\ {\overline{Q}}_{21}& {\overline{Q}}_{22}& {\overline{Q}}_{26}& 0& 0\\ {\overline{Q}}_{16}& {\overline{Q}}_{26}& {\overline{Q}}_{66}& 0& 0\\ 0& 0& 0& {\overline{Q}}_{44}& {\overline{Q}}_{45}\\ 0& 0& 0& {\overline{Q}}_{45}& {\overline{Q}}_{55}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\\ {\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right\}$$

(13)

In which $\mu =(e_0{a}^2)$ and transformed stiffness coefficients can be mentioned as:

$$\begin{array}{l}{\overline{Q}}_{11}^{(k)}=Q_{11}^{(k)}{\mathrm{c}}^4+2(Q_{12}^{(k)}+2Q_{66}^{(k)}){\mathrm{c}}^2{\mathrm{s}}^2+Q_{22}^{(k)}{\mathrm{s}}^4\\ {\overline{Q}}_{12}^{(k)}=(Q_{11}^{(k)}+Q_{22}^{(k)}-4Q_{66}^{(k)}){\mathrm{c}}^2{\mathrm{s}}^2+Q_{12}^{(k)}({\mathrm{c}}^4+{\mathrm{s}}^4)\\ {\overline{Q}}_{22}^{(k)}=Q_{11}^{(k)}{\mathrm{s}}^4+2(Q_{12}^{(k)}+2Q_{66}^{(k)}){\mathrm{c}}^2{\mathrm{s}}^2+Q_{22}^{(k)}{\mathrm{c}}^4\\ {\overline{Q}}_{16}^{(k)}=(Q_{11}^{(k)}-Q_{12}^{(k)}-2Q_{66}^{(k)}){\mathrm{c}}^3\mathrm{s}+(Q_{11}^{(k)}\\ -Q_{22}^{(k)}+2Q_{66}^{(k)})\mathrm{c}{\mathrm{s}}^3\\ {\overline{Q}}_{26}^{(k)}=(Q_{11}^{(k)}-Q_{12}^{(k)}-2Q_{66}^{(k)})\mathrm{c}{\mathrm{s}}^3+(Q_{11}^{(k)}\\ -Q_{22}^{(k)}+2Q_{66}^{(k)}){\mathrm{c}}^3\mathrm{s}\\ {\overline{Q}}_{66}^{(k)}=(Q_{11}^{(k)}+Q_{22}^{(k)}-2Q_{12}^{(k)}){\mathrm{c}}^2{\mathrm{s}}^2+Q_{66}^{(k)}{({\mathrm{c}}^2-{\mathrm{s}}^2)}^2\\ {\overline{Q}}_{44}^{(k)}=Q_{44}^{(k)}{\mathrm{c}}^2+Q_{55}^{(k)}{\mathrm{s}}^2\\ {\overline{Q}}_{45}^{(k)}=(Q_{44}^{(k)}-Q_{55}^{(k)})\mathrm{c}\mathrm{s}\\ {\overline{Q}}_{55}^{(k)}=Q_{44}^{(k)}{\mathrm{s}}^2+Q_{55}^{(k)}{\mathrm{c}}^2\end{array}$$

(14)

where c is $\cos \theta$ , s is $\sin \theta$ , and $\theta$ is the angle between the x-axis and the material orientation. The k-th orthotropic lamina’s terms of engineering constants are:

$$\begin{array}{l}{Q}_{11}^{(k)}=\frac{E_1^{(k)}}{1-{\upsilon }_{12}{\upsilon }_{21}},Q_{12}^{(k)}=\frac{{\upsilon }_{12}{E}_2^{(k)}}{1-{\upsilon }_{12}{\upsilon }_{21}},\\ Q_{22}^{(k)}=\frac{E_2^{(k)}}{1-{\upsilon }_{12}{\upsilon }_{21}},\\ Q_{44}^{(k)}=G_{13},Q_{55}^{(k)}=G_{23},Q_{66}^{(k)}=G_{12}\end{array}$$

(15)

Governing equations

The concept of virtual displacements, whose variational statement for the k-th layer specialized as the source of the governing equations for a layer inside the laminate, is:

$$0=\underset{T}{\int }(\delta U+\delta V-\delta K)dt$$

(16)

The governing equations of the nonlocal theory are derived in terms of local stress findings by applying the L operator to both sides of the equations:

$$\begin{array}{l}•\frac{\partial N_{xx}}{\partial x}+\frac{\partial N_{xy}}{\partial y}=0\\ •\frac{\partial N_{yy}}{\partial y}+\frac{\partial N_{xy}}{\partial x}=0\\ •\frac{\partial {\overline{Q}}_x}{\partial x}+\frac{\partial {\overline{Q}}_y}{\partial y}+\frac{\partial }{\partial x}\left[N_{xx}\frac{\partial w_0}{\partial x}+N_{xy}\frac{\partial w_0}{\partial y}\right]\\ +\frac{\partial }{\partial y}\left[N_{yy}\frac{\partial w_0}{\partial y}+N_{xy}\frac{\partial w_0}{\partial x}\right]\\ +\frac{4}{3h^2}\left[\frac{{\partial }^2{P}_{xx}}{\partial x^2}+2\frac{{\partial }^2{P}_{xy}}{\partial x\partial y}+\frac{{\partial }^2{P}_{yy}}{\partial y^2}\right]\\ =-q\left[1-\mu {\nabla }^2\right]\\ •\frac{\partial {\overline{M}}_{xx}}{\partial x}+\frac{\partial {\overline{M}}_{xy}}{\partial y}-{\overline{Q}}_x=0\\ •\frac{\partial {\overline{M}}_{yy}}{\partial y}+\frac{\partial {\overline{M}}_{xy}}{\partial x}-{\overline{Q}}_y=0\end{array}$$

(17)

The following constitutive equations are created by adding the stress-strain relations to equations (12)–(17).

$$\begin{array}{l}\left\{\begin{array}{c}{N}_{xx}\\ N_{yy}\\ N_{xy}\end{array}\right\}=\left\{\begin{array}{ccc}{A}_{11}& A_{12}& A_{16}\\ A_{12}& A_{22}& A_{26}\\ A_{16}& A_{26}& A_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(0)}\\ {\epsilon }_{yy}^{(0)}\\ {\gamma }_{xy}^{(0)}\end{array}\right\}\\ +\left\{\begin{array}{ccc}{B}_{11}& B_{12}& B_{16}\\ B_{12}& B_{22}& B_{26}\\ B_{16}& B_{26}& B_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(1)}\\ {\epsilon }_{yy}^{(1)}\\ {\gamma }_{xy}^{(1)}\end{array}\right\}\\ +\left\{\begin{array}{ccc}{E}_{11}& E_{12}& E_{16}\\ E_{12}& E_{22}& E_{26}\\ E_{16}& E_{26}& E_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(3)}\\ {\epsilon }_{yy}^{(3)}\\ {\gamma }_{xy}^{(3)}\end{array}\right\}\end{array}$$

(18a)

$$\left\{\begin{array}{c}{M}_{xx}\\ M_{yy}\\ M_{xy}\end{array}\right\}=\left\{\begin{array}{ccc}{B}_{11}& B_{12}& B_{16}\\ B_{12}& B_{22}& B_{26}\\ B_{16}& B_{26}& B_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(0)}\\ {\epsilon }_{yy}^{(0)}\\ {\gamma }_{xy}^{(0)}\end{array}\right\}+\left\{\begin{array}{ccc}{D}_{11}& D_{12}& D_{16}\\ D_{12}& D_{22}& D_{26}\\ D_{16}& D_{26}& D_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(1)}\\ {\epsilon }_{yy}^{(1)}\\ {\gamma }_{xy}^{(1)}\end{array}\right\}+\left\{\begin{array}{ccc}{F}_{11}& F_{12}& F_{16}\\ F_{12}& F_{22}& F_{26}\\ F_{16}& F_{26}& F_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(3)}\\ {\epsilon }_{yy}^{(3)}\\ {\gamma }_{xy}^{(3)}\end{array}\right\}$$

(18b)

$$\left\{\begin{array}{c}{P}_{xx}\\ P_{yy}\\ P_{xy}\end{array}\right\}=\left\{\begin{array}{ccc}{E}_{11}& E_{12}& E_{16}\\ E_{12}& E_{22}& E_{26}\\ E_{16}& E_{26}& E_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(0)}\\ {\epsilon }_{yy}^{(0)}\\ {\gamma }_{xy}^{(0)}\end{array}\right\}+\left\{\begin{array}{ccc}{F}_{11}& F_{12}& F_{16}\\ F_{12}& F_{22}& F_{26}\\ F_{16}& F_{26}& F_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(1)}\\ {\epsilon }_{yy}^{(1)}\\ {\gamma }_{xy}^{(1)}\end{array}\right\}+\left\{\begin{array}{ccc}{H}_{11}& H_{12}& H_{16}\\ H_{12}& H_{22}& H_{26}\\ H_{16}& H_{26}& H_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_{xx}^{(3)}\\ {\epsilon }_{yy}^{(3)}\\ {\gamma }_{xy}^{(3)}\end{array}\right\}$$

(18c)

$$\left\{\begin{array}{c}{Q}_{yz}\\ Q_{xz}\end{array}\right\}=\left\{\begin{array}{cc}{A}_{44}& A_{45}\\ A_{45}& A_{55}\end{array}\right\}\left\{\begin{array}{c}{\gamma }_{yz}^{(0)}\\ {\gamma }_{xz}^{(0)}\end{array}\right\}+\left\{\begin{array}{cc}{D}_{44}& D_{45}\\ D_{45}& D_{55}\end{array}\right\}\left\{\begin{array}{c}{\gamma }_{yz}^{(2)}\\ {\gamma }_{xz}^{(2)}\end{array}\right\}$$

(18d)

$$\left\{\begin{array}{c}{R}_{yz}\\ R_{xz}\end{array}\right\}=\left\{\begin{array}{cc}{D}_{44}& D_{45}\\ D_{45}& D_{55}\end{array}\right\}\left\{\begin{array}{c}{\gamma }_{yz}^{(0)}\\ {\gamma }_{xz}^{(0)}\end{array}\right\}+\left\{\begin{array}{cc}{F}_{44}& F_{45}\\ F_{45}& F_{55}\end{array}\right\}\left\{\begin{array}{c}{\gamma }_{yz}^{(2)}\\ {\gamma }_{xz}^{(2)}\end{array}\right\}$$

(18e)

$$\left\{A_{ij};B_{ij};D_{ij};E_{ij};F_{ij};H_{ij}\right\}=\sum _{k=1}^N\underset{z_k}{\overset{z_{k+1}}{\int }}{\overline{Q}}_{ij}^{(k)}(1;z;z^2;z^3;z^4;z^6)dz(i,j=1,2,6)$$

(19a)

$$\left\{A_{ij};D_{ij};F_{ij}\right\}=\sum _{k=1}^N\underset{z_k}{\overset{z_{k+1}}{\int }}{\overline{Q}}_{ij}^{(k)}(1;z^2;z^4)dz(i,j=4,5)$$

(19b)

Finite element analysis

The generalized displacements are approximated over an element ${\Omega }^e$ by the expressions below:

$$\begin{array}{l}{u}_0(x,y,t)=\sum _{i=1}^m{u}_i^e(t){\psi }_i^e(x,y)\\ v_0(x,y,t)=\sum _{i=1}^m{v}_i^e(t){\psi }_i^e(x,y)\\ w_0(x,y,t)=\sum _{i=1}^m{\overline{\Delta }}_i^e(t){\phi }_i^e(x,y)\\ {\varphi }_x(x,y,t)=\sum _{i=1}^m{X}_i^e(t){\psi }_i^e(x,y)\\ {\varphi }_y(x,y,t)=\sum _{i=1}^m{Y}_i^e(t){\psi }_i^e(x,y)\end{array}$$

(20)

Substitution of equation (20) into the weak form equation (17) yields the finite element equations for bending.

$$\left[\begin{array}{ccccc}\left[K^{11}\right]& \left[K^{12}\right]& \left[K^{13}\right]& \left[K^{14}\right]& \left[K^{15}\right]\\ {\left[K^{12}\right]}^T& \left[K^{22}\right]& \left[K^{23}\right]& \left[K^{24}\right]& \left[K^{25}\right]\\ {\left[K^{13}\right]}^T& {\left[K^{23}\right]}^T& \left[K^{33}\right]& \left[K^{34}\right]& \left[K^{35}\right]\\ {\left[K^{14}\right]}^T& {\left[K^{24}\right]}^T& {\left[K^{34}\right]}^T& \left[K^{44}\right]& \left[K^{45}\right]\\ {\left[K^{15}\right]}^T& {\left[K^{25}\right]}^T& {\left[K^{35}\right]}^T& {\left[K^{45}\right]}^T& \left[K^{55}\right]\end{array}\right]\left\{\begin{array}{c}\left\{u^e\right\}\\ \left\{v^e\right\}\\ \left\{{\Delta }^e\right\}\\ \left\{X^e\right\}\\ \left\{Y^e\right\}\end{array}\right\}=\left\{\begin{array}{c}\left\{F^1\right\}\\ \left\{F^2\right\}\\ \left\{F^3\right\}\\ \left\{F^4\right\}\\ \left\{F^5\right\}\end{array}\right\}$$

(21)

The linear stiffness coefficients $K_{ij}^{\alpha \beta }$ are defined by:

$$\begin{array}{l}{K}_{ij}^{1\alpha }=\underset{{\Omega }^e}{\int }\left(\frac{\partial {\psi }_i}{\partial x}{N}_{1j}^{\alpha }+\frac{\partial {\psi }_i}{\partial y}{N}_{6j}^{\alpha }\right)dxdy\\ K_{ij}^{2\alpha }=\underset{{\Omega }^e}{\int }\left(\frac{\partial {\psi }_i}{\partial x}{N}_{6j}^{\alpha }+\frac{\partial {\psi }_i}{\partial y}{N}_{2j}^{\alpha }\right)dxdy\\ K_{ij}^{3\alpha }=\underset{{\Omega }^e}{\int }\left(\begin{array}{l}\frac{\partial {\phi }_i}{\partial x}{\widehat{Q}}_{1j}^{\alpha }+\frac{\partial {\phi }_i}{\partial y}{\widehat{Q}}_{2j}^{\alpha }\\ -\frac{4}{3h}\left(\frac{{\partial }^2{\phi }_i}{\partial x^2}{P}_{1j}^{\alpha }+2\frac{{\partial }^2{\phi }_i}{\partial x\partial y}{P}_{6j}^{\alpha }+\frac{{\partial }^2{\phi }_i}{\partial y^2}{P}_{2j}^{\alpha }\right)\end{array}\right)\\ dxdy\\ K_{ij}^{4\alpha }=\underset{{\Omega }^e}{\int }\left(\frac{\partial {\psi }_i}{\partial x}{\widehat{M}}_{1j}^{\alpha }+\frac{\partial {\psi }_i}{\partial y}{\widehat{M}}_{6j}^{\alpha }+{\psi }_i{\widehat{Q}}_{1j}^{\alpha }\right)dxdy\\ K_{ij}^{5\alpha }=\underset{{\Omega }^e}{\int }\left(\frac{\partial {\psi }_i}{\partial x}{\widehat{M}}_{6j}^{\alpha }+\frac{\partial {\psi }_i}{\partial y}{\widehat{M}}_{2j}^{\alpha }+{\psi }_i{\widehat{Q}}_{2j}^{\alpha }\right)dxdy\end{array}$$

(22)

The above coefficients are defined using the notation indicated below:

$$\begin{array}{l}{N}_{1j}^1=A_{11}\frac{\partial {\psi }_j}{\partial x}+A_{16}\frac{\partial {\psi }_j}{\partial y};N_{1j}^2=A_{16}\frac{\partial {\psi }_j}{\partial x}+A_{12}\frac{\partial {\psi }_j}{\partial y};\\ N_{1j}^3=-\frac{4}{3h}\left(E_{11}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2E_{16}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+E_{12}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ N_{1j}^4={\widehat{B}}_{11}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{16}\frac{\partial {\psi }_j}{\partial y};N_{1j}^5={\widehat{B}}_{12}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{16}\frac{\partial {\psi }_j}{\partial y};\\ N_{6j}^1=A_{16}\frac{\partial {\psi }_j}{\partial x}+A_{66}\frac{\partial {\psi }_j}{\partial y};N_{6j}^2=A_{66}\frac{\partial {\psi }_j}{\partial x}+A_{26}\frac{\partial {\psi }_j}{\partial y};\\ N_{6j}^3=-\frac{4}{3h}\left(E_{16}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2E_{66}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+E_{26}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ N_{6j}^4={\widehat{B}}_{16}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{66}\frac{\partial {\psi }_j}{\partial y};N_{6j}^5={\widehat{B}}_{26}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{66}\frac{\partial {\psi }_j}{\partial y};\\ N_{2j}^1=A_{12}\frac{\partial {\psi }_j}{\partial x}+A_{26}\frac{\partial {\psi }_j}{\partial y};N_{2j}^2=A_{26}\frac{\partial {\psi }_j}{\partial x}+A_{22}\frac{\partial {\psi }_j}{\partial y};\\ N_{2j}^3=-\frac{4}{3h}\left(E_{12}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2E_{26}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+E_{22}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ N_{2j}^4={\widehat{B}}_{12}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{26}\frac{\partial {\psi }_j}{\partial y};N_{2j}^5={\widehat{B}}_{22}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{26}\frac{\partial {\psi }_j}{\partial y};\end{array}$$

(23a)

$$\begin{array}{l}{\widehat{M}}_{1j}^1={\widehat{B}}_{11}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{16}\frac{\partial {\psi }_j}{\partial y};{\widehat{M}}_{1j}^2={\widehat{B}}_{16}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{12}\frac{\partial {\psi }_j}{\partial y};\\ {\widehat{M}}_{1j}^3=-\frac{4}{3h}\left({\widehat{F}}_{11}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2{\widehat{F}}_{16}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+{\widehat{F}}_{12}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ {\widehat{M}}_{1j}^4={\widehat{D}}_{11}\frac{\partial {\psi }_j}{\partial x}+{\widehat{D}}_{16}\frac{\partial {\psi }_j}{\partial y};{\widehat{M}}_{1j}^5={\widehat{D}}_{12}\frac{\partial {\psi }_j}{\partial x}+{\widehat{D}}_{16}\frac{\partial {\psi }_j}{\partial y};\\ {\widehat{M}}_{6j}^1={\widehat{B}}_{16}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{66}\frac{\partial {\psi }_j}{\partial y};{\widehat{M}}_{6j}^2={\widehat{B}}_{66}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{26}\frac{\partial {\psi }_j}{\partial y};\\ {\widehat{M}}_{6j}^3=-\frac{4}{3h}\left({\widehat{F}}_{16}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2{\widehat{F}}_{66}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+{\widehat{F}}_{26}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ {\widehat{M}}_{6j}^4={\widehat{D}}_{16}\frac{\partial {\psi }_j}{\partial x}+{\widehat{D}}_{66}\frac{\partial {\psi }_j}{\partial y};{\widehat{M}}_{6j}^5={\widehat{D}}_{26}\frac{\partial {\psi }_j}{\partial x}+{\widehat{D}}_{66}\frac{\partial {\psi }_j}{\partial y};\\ {\widehat{M}}_{2j}^1={\widehat{B}}_{12}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{26}\frac{\partial {\psi }_j}{\partial y};{\widehat{M}}_{2j}^2={\widehat{B}}_{26}\frac{\partial {\psi }_j}{\partial x}+{\widehat{B}}_{22}\frac{\partial {\psi }_j}{\partial y};\\ {\widehat{M}}_{2j}^3=-\frac{4}{3h}\left({\widehat{F}}_{12}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2{\widehat{F}}_{26}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+{\widehat{F}}_{22}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ {\widehat{M}}_{2j}^4={\widehat{D}}_{12}\frac{\partial {\psi }_j}{\partial x}+{\widehat{D}}_{26}\frac{\partial {\psi }_j}{\partial y};{\widehat{M}}_{2j}^5={\widehat{D}}_{22}\frac{\partial {\psi }_j}{\partial x}+{\widehat{D}}_{26}\frac{\partial {\psi }_j}{\partial y};\end{array}$$

(23b)

$$\begin{array}{l}{P}_{1j}^1=E_{11}\frac{\partial {\psi }_j}{\partial x}+E_{16}\frac{\partial {\psi }_j}{\partial y};P_{1j}^2=E_{12}\frac{\partial {\psi }_j}{\partial x}+E_{16}\frac{\partial {\psi }_j}{\partial y};\\ P_{1j}^3=-\frac{4}{3h}\left(H_{11}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2H_{16}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+H_{12}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ P_{1j}^4={\widehat{F}}_{11}\frac{\partial {\psi }_j}{\partial x}+{\widehat{F}}_{16}\frac{\partial {\psi }_j}{\partial y};P_{1j}^5={\widehat{F}}_{12}\frac{\partial {\psi }_j}{\partial x}+{\widehat{F}}_{16}\frac{\partial {\psi }_j}{\partial y};\\ P_{2j}^1=E_{12}\frac{\partial {\psi }_j}{\partial x}+E_{26}\frac{\partial {\psi }_j}{\partial y};P_{2j}^2=E_{22}\frac{\partial {\psi }_j}{\partial x}+E_{26}\frac{\partial {\psi }_j}{\partial y};\\ P_{2j}^3=-\frac{4}{3h}\left(H_{21}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2H_{26}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+H_{22}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ P_{2j}^4={\widehat{F}}_{12}\frac{\partial {\psi }_j}{\partial x}+{\widehat{F}}_{26}\frac{\partial {\psi }_j}{\partial y};P_{2j}^5={\widehat{F}}_{22}\frac{\partial {\psi }_j}{\partial x}+{\widehat{F}}_{26}\frac{\partial {\psi }_j}{\partial y};\\ P_{6j}^1=E_{16}\frac{\partial {\psi }_j}{\partial x}+E_{66}\frac{\partial {\psi }_j}{\partial y};P_{6j}^2=E_{26}\frac{\partial {\psi }_j}{\partial x}+E_{66}\frac{\partial {\psi }_j}{\partial y};\\ P_{6j}^3=-\frac{4}{3h}\left(H_{16}\frac{{\partial }^2{\phi }_j}{\partial x^2}+2H_{66}\frac{{\partial }^2{\phi }_j}{\partial x\partial y}+H_{26}\frac{{\partial }^2{\phi }_j}{\partial y^2}\right);\\ P_{6j}^4={\widehat{F}}_{16}\frac{\partial {\psi }_j}{\partial x}+{\widehat{F}}_{66}\frac{\partial {\psi }_j}{\partial y};P_{6j}^5={\widehat{F}}_{26}\frac{\partial {\psi }_j}{\partial x}+{\widehat{F}}_{66}\frac{\partial {\psi }_j}{\partial y};\end{array}$$

(23c)

$$\begin{array}{l}{\widehat{Q}}_{1j}^3={\widehat{A}}_{45}\frac{\partial {\phi }_j}{\partial x}+{\widehat{A}}_{55}\frac{\partial {\phi }_j}{\partial y};\\ {\widehat{Q}}_{2j}^3={\widehat{A}}_{44}\frac{\partial {\phi }_j}{\partial x}+{\widehat{A}}_{45}\frac{\partial {\phi }_j}{\partial y};\\ {\widehat{Q}}_{1j}^4={\widehat{A}}_{55}{\psi }_j;{\widehat{Q}}_{1j}^5={\widehat{Q}}_{2j}^4={\widehat{A}}_{45}{\psi }_j\end{array}$$

(23d)

The elements of the force vector are given by:

$$\begin{array}{l}{F}_i^1=\underset{{\Gamma }^e}{\oint }\left(N_{xx}{\widehat{n}}_x+N_{xy}{\widehat{n}}_y\right)ds\\ F_i^2=\underset{{\Gamma }^e}{\oint }\left(N_{xy}{\widehat{n}}_x+N_{xy}{\widehat{n}}_y\right)ds\\ F_i^3=\underset{{\Omega }^e}{\int }(1-\mu {\nabla }^2)q{\phi }_{i}dxdy+\underset{{\Gamma }^e}{\oint }\left\{\begin{array}{l}({\overline{Q}}_x{\widehat{n}}_x+{\overline{Q}}_y{\widehat{n}}_y)+\left(N_{xx}\frac{\partial w_0}{\partial x}+N_{xy}\frac{\partial w_0}{\partial y}\right){\widehat{n}}_x\\ +\left(N_{yy}\frac{\partial w_0}{\partial y}+N_{xy}\frac{\partial w_0}{\partial x}\right){\widehat{n}}_y\\ +\frac{4}{3h^2}\left[\left(\frac{\partial P_{xx}}{\partial x}+\frac{\partial P_{xy}}{\partial y}\right){\widehat{n}}_x+\left(\frac{\partial P_{yy}}{\partial y}+\frac{\partial P_{xy}}{\partial x}\right){\widehat{n}}_y\right]\\ -\frac{4}{3h^2}\left[\left(P_{xx}+P_{xy}\right){\widehat{n}}_x+\left(P_{yy}+P_{xy}\right){\widehat{n}}_y\right]\end{array}\right\}ds\\ F_i^4=\underset{{\Gamma }^e}{\oint }\left(M_{xx}{\widehat{n}}_x+M_{xy}{\widehat{n}}_y\right){\psi }_{i}ds\\ F_i^5=\underset{{\Gamma }^e}{\oint }\left(M_{xy}{\widehat{n}}_x+M_{yy}{\widehat{n}}_y\right){\psi }_{i}ds\end{array}$$

(24)

Numerical study

The aim of the finite element analysis is to verify the obtained results. For this purpose, ABAQUS was utilized to construct numerical models. Shell elements with S4R were used for modeling. The Hashin failure criteria were adopted to achieve failure of the plates. Tensile fiber mode, fiber compressive mode, tensile matrix mode, and compressive matrix mode are the failure modes defined in the Hashin failure criteria. These failure modes can be computed from the equations listed in Table 3. The pultruded GFRP was modeled employing the Hashin failure criteria, and the parameters used in the Hashin failure are presented in Table 4. ⁹⁴ The supports and loading noses were made from solid elements with C3D8R. Steel materials properties with elastic modulus of 200 GPa and yielding of 420 MPa were utilized for the supports. All degrees of freedom were restrained under the supports. Frictional surfaces were considered between the pultruded plate and the supports and loading noses. Friction coefficients of 0.5 was considered. ³⁸ Figure 6 shows the typical numerical modeling.

Table 3.

Hashin failure criteria.

Mode description	Equation
Tensile fiber mode	$F_f^t={\left(\frac{{\sigma }_{11}}{X^T}\right)}^2+{\left(\frac{{\tau }_{12}}{S^L}\right)}^2$
Fiber compressive mode	$F_f^c={\left(\frac{{\sigma }_{11}}{X^C}\right)}^2$
Tensile matrix mode	$F_m^t={\left(\frac{{\sigma }_{22}}{Y^T}\right)}^2+{\left(\frac{{\tau }_{12}}{S^L}\right)}^2$
Compressive matrix mode	$F_m^c={\left(\frac{{\sigma }_{22}}{2S^T}\right)}^2+\left[{\left(\frac{Y^C}{2S^T}\right)}^2-1\right]\frac{{\sigma }_{22}}{Y^C}+{\left(\frac{{\tau }_{12}}{S^L}\right)}^2$
where σ11 and σ22 are components of the effective stress tensor, XT and XC are th e tensile and compressive strengths in the direction of fiber, YT and YC are the tensile and compressive strengths in the direction of matrix, SL and ST are the shear strength of longitudinal and transverse directions, respectively.

Table 4.

Parameters used in Hashin failure criteria.

Property of CFRP composites		Used value
Longitudinal elasticity modulus (GPa)	E 1	23,000
Transverse elasticity modulus (GPa)	E 2	7000
Poisson ratio	v 12	0.21
Shear modulus (GPa)	G 12	3400
Shear modulus (GPa)	G 13	3000
Shear modulus (GPa)	G 23	2000
Tensile strength in fiber direction (MPa)	X T	290
Compressive strength in fiber direction (MPa)	X C	240
Tensile strength in matrix direction (MPa)	Y T	34
Compressive strength in matrix direction (MPa)	Y C	70
Longitudinal shear strength (MPa)	S L	25
Transverse shear strength (MPa)	S T	25
Energy required to fully damage a ply using fiber tension only (mJ/mm2)	G ft	70
Energy required to fully damage a ply using fiber compression only (mJ/mm2)	G fc	50
Energy required to fully damage a ply using matrix tension only (mJ/mm2)	G mt	2
Energy required to fully damage a ply using matrix compression only (mJ/mm2)	G mc	10

Figure 6.

Examples of numerical model and mesh.

The load–displacement graphs obtained from the numerical simulations are exhibited in Figures 7, 9, 11, 13, 15, 17, 19, and 21. The results display that the numerical models captured the obtained results. The accuracy of the numerical models for specimens under SFSF conditions is higher than the specimens under CFCF. The numerical results exhibit higher capacity than those of the tested specimens. The reason for these can be attributed to the setup of the tests. Obtaining perfect clamped boundary conditions in the tests is difficult. Slight rotation of boundary conditions may result in decrease in the capacity.

Figure 7.

Load–displacement graphs of specimen SFSF_D_0°.

Results and discussion

Material properties

A pultruded GFRP plate that had been employed in the trials served as the source of the specimen for the burnout test. The specimen was weighed, and volume/density calculations were done before the burnout test. The matrix was totally burnt off as a result of this procedure, leaving just the glass roving, glass mat, and CaCO₃ behind. These items were weighed once more. The pultruded GFRP was evaluated and validated in terms of fiber configuration and stacking order after the burnout test. Fiber and matrix volume fractions are found to be 54% and 46%, respectively. E₁ and E₂ are 20 GPa and 6 GPa, respectively. Tensile and compressive strengths in the longitudinal direction is 320 MPa while these are 80 MPa in the transverse direction. Detailed information about the material properties was reported in theprevious studies.^14,15

Macro damage analysis

Load–displacement graphs of the tested pultruded GFRP plates are illustrated in Figures 7, 9, 11, 13, 15, 17, 19, and 21. The tests for each configuration were repeated 3 times. Damage analyses were done based on macro damage photos in Figures 8, 10, 12, 14, 16, 18, 20, and 22. Figures 7–14 display the outcomes for the specimens with distributed loading.

Figure 8.

Damages of specimen SFSF_D_0°.

Figure 9.

Load–displacement graphs of specimen SFSF_D_90°.

Figure 10.

Damages of specimen SFSF_D_90°.

Figure 11.

Load–displacement graphs of specimen CFCF_D_0°.

Figure 12.

Damages of specimen CFCF_D_0°.

Figure 13.

Load–displacement graphs of specimen CFCF_D_90°.

Figure 14.

Damages of specimen CFCF_D_90°.

Figures 7 and 8 demonstrate the specimens with 0° under SFSF boundary conditions. The maximum loading was measured as 2.3 kN on average. The first damage was formed in surfacing veil layer. As loading increased, matrix cracks in the fiber directions of 90°, in other words, intralaminar cracks were progressed toward the center of the plate. On the other hand, Figures 9 and 10 depict the specimens with 90° under the same loading and boundary conditions. This time, the maximum loading modified to 10.3 kN. The progressive damage was totally different due to the fiber directions. Initial damage was occurred at veil surfaces. Then, large delamination damage (interlaminar cracks) formed. After this damage, fibers were ruptured, and delamination damage was seen in different locations.

In Figures 11 and 12, clamped boundary conditions were utilized. The damage pattern was similar to the SFSF case. The main failure mode was intralaminar cracks. The maximum loading was measured as 5.4 kN on average which was 134% higher that the case with SFSF. On the other hand, Figures 13 and 14 indicate the specimens with 90° under same loading and boundary conditions. Initially, intralaminar cracks were formed. Then, surface veil damage was progressed. As the loading increased, fibers were triggered to rupture. The maximum loading was measured as 16.5 kN on average.

Figures 15–18 present the findings for the specimens with point loading. Figures 15 and 16 depict the specimens with 0° under SFSF boundary conditions. When the specimen reached its capacity, intralaminar and interlaminar cracks were seen at the center line simultaneously. As the loading increased, intralaminar cracks were seen at different locations toward the supports. The maximum loading was recorded as 1.7 kN on average. On the other hand, Figures 17 and 18 show the specimen with 90° and the same conditions. In this case, a local damage was witnessed at the center of the specimen. Local indication damage was initially detected underneath the loading nose. As the loading increased, the fibers started to break. The damage was not able to propagate toward the ends. The maximum loading was 6.9 kN which was approximately 3 times higher than the specimens with 0° and the same conditions.

Figure 15.

Load–displacement graphs of specimen SFSF_P_0°.

Figure 16.

Damages of specimen SFSF_P_0°.

Figure 17.

Load–displacement graphs of specimen SFSF_P_90°.

Figure 18.

Damages of specimen SFSF_P_90°.

Figures 19–22 illustrate the findings for the specimens with clamped boundary conditions. Figures 19 and 20 display the specimens with 0°. The maximum loading was recorded as 3.8 kN. At this loading, intralaminar and interlaminar cracks were observed at the center line simultaneously. Interlaminar cracks were progressed as the loading increased. Figures 21 and 22 indicate the specimens with 90° and the same conditions. The local indication damage was initially detected. Then, intralaminar and interlaminar cracks were witnessed at the center line at the same time. As the loading increased, the fibers underneath the loading nose started to break and intralaminar and interlaminar cracks progressed to grow. The maximum loading was around 7.6 kN which was double of the capacity of the specimens with 0° and the same conditions.

Figure 19.

Load–displacement graphs of specimen CFCF_P_0°.

Figure 20.

Damages of specimen CFCF_P_0°.

Figure 21.

Load–displacement graphs of specimen CFCF_P_90°.

Figure 22.

Damages of specimen CFCF_P_90°.

Figure 23 demonstrates the failure modes obtained from numerical study. It is seen that the predicted failure modes are very similar to those observed in the experimental results. HSNFTCRT = 1 represents the fiber tension damage and HSNMTCRT = 1 points out the matrix tension damage. As similar to the obtained experimental results, the specimens with 0° suffered from matrix related damages while the specimens with 90° failed from fiber breakage.

Figure 23.

Failure modes of FEM modeling.

Static analysis

For composite plates under bending load, the theoretical nonlocal high-order FEM and FEM ABAQUS findings were converged. The displacement results were obtained at the halfway point and compared to the experimental displacement findings. The maximum displacement values of the pultruded GFRP plate are presented in Tables 5 and 6. The nonlocal parameter μ has been given a value of 0 to 5. It has been shown that the nonlocality effect has a considerable impact on the displacement value. The value of the displacement increases as the nonlocal parameter rises. Displacement increased as the layer angle changed from 0° to 90°. All the analysis methods captured similar numerical values. As the nonlocal parameter decreased, the displacement value also decreased. If we take the numerical results obtained experimentally as reference, FEM ABAQUS and nonlocal FEM results are very close to the experimental results.

Table 5.

Center displacements (mm) of a pultruded GFRP plate exposed to distributed loads as a result of lower integration and nonlocality.

Boundary conditions	Fiber orientation	Experimental	Numerical	Theoretical nonlocal
				µ = 0	µ = 1	µ = 3	µ = 5
SFSF	0°	4.89	4.93	4.94	4.98	5.06	5.14
	90°	9.14	9.16	9.18	9.22	9.3	9.37
CFCF	0°	4.14	4.18	4.19	4.23	4.31	4.38
	90°	5.94	5.96	5.98	6.02	6.10	6.18

Table 6.

Center displacements (mm) of a pultruded GFRP plate exposed to point loads as a result of lower integration and nonlocality.

Boundary conditions	Fiber orientation	Experimental	Numerical	Theoretical nonlocal
				µ = 0	µ = 1	µ = 3	µ = 5
SFSF	0°	4.28	4.30	4.32	4.36	4.43	4.49
	90°	6.97	7.0	7.12	7.16	7.23	7.30
CFCF	0°	7.99	8.02	8.24	8.28	8.35	8.41
	90°	7.54	7.56	7.59	7.63	7.7	7.69

Conclusions

In the present study, by taking into account the Eringen’s nonlocal differential constitutive equations, the governing equations of higher-order deformation plate theory are reformulated using the virtual work principle for the pultruded GFRP composite plates. The nonlocal finite element model is created and the weak forms are obtained. The burnout test, tensile test, and compressive test were conducted to determine the material properties. For the accuracy of the experimental and theoretical results, FEM modeling and analysis were performed in the ABAQUS program. The impacts of anisotropy and nonlocality on the bending behavior of laminated composite plates were illustrated using examples with various stacking sequences and boundary conditions.

• • It has been noted that the Eringen’s nonlocal model causes laminated composite plates to lose some of their structural stiffness; as a result, the displacement rises as the nonlocal parameter is increased.

• • The through-the-thickness displacement is constant for 2D elasticity calculations, but the employed quasi-3D technique assesses any potential change in the kinematic response for both negative and positive layers. Most of the times, the top surface of the plate experiences the largest displacement.

• • As expected, the specimens with CFCF boundary conditions exhibited 3.5 times higher capacity than that of the specimens with SFSF boundary conditions. The specimens exposed to distributed loading depicted 61% higher capacity than the specimens exposed to point loading

• • The main failure mode of the specimens with 0° was intralaminar cracks regardless of loading type and boundary conditions. On the other hand, the main failure modes of the specimens with 90° were interlaminar cracks (delimination) and fiber breakage.

• • When the specimens were subjected to the point loading, local indication damages were detected in the specimens with 90°. For the specimens with 0°, the damage patterns were almost similar regardless of loading type.

• • The numerical models with the Hashin failure criteria were able to capture the behavior of the pultruded GFRP plates. The accuracy of the numerical results was higher for the specimens with SFSF boundary conditions. This can be attributed for the difficulties to ensure clamped boundary conditions experimentally.