A 3D computational meshfree model for the mechanical and thermal buckling analysis of rectangular composite laminated plates with embedded delaminations

2.1 Point interpolation using RBFs

Considering that a scalar functions u(x) that is defined in problem domain Ω represented by a set of scattered nodes, the point interpolation augmented with polynomials can be written as follows:

with the constraint condition:

where R_i (r) is the RBF, k is the number of points in the neighborhood of the point of interest x, Ψ_j (x) is monomial in the space coordinates x^T=[xy], and m is the number of polynomial basis functions. When m=0, pure RBFs are used. Otherwise, the RBF is augmented with m polynomial basis functions. Coefficients a_i and b_j are interpolation constants to be determined. In the RBF R_i (x), the variable is only the distance between the point of interest x and a node at x_i , i.e. r=(x-xi)2+(y-yi)2.

Coefficients a_i and b_j can be determined by enforcing Equation (1) to be satisfied at these n nodes surrounding the point of interest x, which leads to n linear equations. Equation (1) can be rewritten in matrix form as

where

The following expression can be obtained:

In the following section, Ψ^T(x)=[1 xy] is adopted. Then, the shape function Φ(x) can be defined by

and

2.2 LRPIM formula of the Hamilton canonical equation

For isotropic, orthotropic, or anisotropic elasticity solids, a local weak form equivalent to the modified H-R variational principle in 3D Cartesian coordinate system can be expressed as follows [43, 46, 47]:

where Q=[uvw]^T is the displacement vector and u, v, and w are the displacement components along x, y, and z coordinates, respectively. P=[σ_xzσ_yzσ_zz ]^T is the stress vector, σ_xz , σ_yz , and σ_zz are the transverse stresses. The superscript T signifies a matrix transposition. H is Hamiltonian. v is referred to the volume considered. S_A denotes the surface over the volume V. B_pq̅=[p_x (u-u̅) p_y (v-v̅) p_z (w-w̅)]^T, where p_i (i=x, y, z) are the stress boundary conditions in three coordinate directions. u̅, v̅, and w̅ indicate the prescribed displacement boundary conditions along x, y, and z coordinates, respectively. B_p̅_q =[p̅_xup̅_yvp̅_zw]^T, where p̅_i (i=x, y, z) are the prescribed stress boundary conditions in x, y, and z directions, respectively. The characteristic coefficients λ₁=[λ_x -1 λ_y -1 λ_z -1]^T and λ₀=[λ_xλ_yλ_z ]^T are introduced especially. The value of λ_i (i=x, y, z) are 1 or 0. λ_i (i=x, y, z)=1 denotes the stress boundary cases in three coordinates directions. λ_i (i=x, y, z)=0 denotes the displacement boundary cases in three coordinates directions.

For the mechanical buckling analysis of a laminate subjected to uniform in-plane pressures p_bx and p_by , respectively, in x and y directions, and for the thermal buckling analysis of a laminate subjected to uniform temperature rise ΔT, the Hamiltonian function H in Equation (10) can be stated as follows:

where

and W⌢=Diag(W⌢), where W⌢ is weight function.

For the isotropic and orthotropic materials,

where s₁=1/c₅₅, s₂=1/c₄₄, s₃=1/c₃₃, s₄=-c₁₃/c₃₃, s₅=-c₂₃/c₃₃, s6=c11-c132/c33,s₇=c₁₂-c₁₃c₂₃/c₃₃, s8=c22-c232/c33,s₉=c₆₆, and c_ij (i, j=1, 2, 3, 4, 5, 6) are known as the stiffness coefficients, and

where λ_m and λ_t are the mechanical and thermal buckling factors, respectively; NTx=∑μ=1N((s6)μ(αx)μ+(s7)μ(αy)μ)hμΔT and NTy=∑μ=1N((s7)μ(αx)μ+(s8)μ(αy)μ)hμΔT denote the thermal loading in x and y directions, respectively; α_x =α₁cos²θ+α₂sin²θ and α_y =α₁sin²θ+α₂cos²θ are the thermal expansion coefficients in the x and y directions, respectively; α₁ and α₂ are the longitudinal and transverse thermal expansion coefficients, respectively; θ denotes the angle between the positive x-axis and the fiber direction, measured in counterclockwise direction; (s₆)_μ , (s₇)_μ , (s₈)_μ , (α_x )_μ , (α_y )_μ , and h_μ are the stiffness coefficients, thermal expansion coefficients, and thickness of the μ-th sublaminate; N is the total layer numbers of the laminate; and F=-[f_xf_yf_z ]^T denote the vectors of body forces or external loadings.

By means of LRPIM shape function, the stress vector P and the displacement vector Q at any quadrature point can be written as follows:

where P_κ =[σ_xzκ (z) σ_yzκ (z) σ_zzκ (z)]^T, Q_κ =[u_κ (z) v_κ (z) w_κ (z)]^T, N=Diag(Φ)_3×3, and notation z is the thickness of the problem domain. The subscript κ denotes the variables of a point.

Substituting Equation (15) into Equation (10) and applying the tool of variational calculus and integrating, the expression of the first term of Equation (10) can be obtained as follows:

where

where

where Ω_q is local quadrature domain, Γ is the global boundary of surface S_A , Γ_qi is the internal boundary of the local quadrature domain, which does not intersect with Γ, Γ_qu is the part of the essential boundary that intersects with the local quadrature domain, Γ_qt is the part of the natural boundary that intersects with the local quadrature domain, V⌢I1 and V⌢I2 are the matrices that collects the derivatives of the weight functions, W⌢I is a matrix of weight functions, and n is a matrix of vectors of the unit outward normal on the boundary.

The boundary term can be reduced to the following form:

where X11=∮Γ(G3N)TNdΓ,X12=X12T=∮Γ(NT(G4N)+(G4N)TN)dΓ,Y11=−∮ΓNTΛ0dΓ,Y12=−∮Γ(G4NT)dΓ, and Y22=∮Γ(G3NT)dΓ, in which

and Λ₀=Diag(λ_x , λ_y , λ_z ).

Adding Equation (19) to the right-hand side of Equation (16), an LRPIM formula of the Hamilton canonical equation for the local quadrature domain can be expressed as

where E11=K11T+X11T,E₁₂=K₁₂+X₁₂, E₂₁=K₂₁, and E₂₂=-K₁₁-X₁₁ are the equivalent stiffness matrices for local quadrature domain and F_κ1(z)=Y₁₁P̅_κ (z)+Y₁₂Q̅_κ (z)+ς_κ and F_κ2(z)=-Y₁₁Q̅(z) are the equivalent external load vectors for local quadrature domain.

2.3 Numerical implementation for LRPIM

Gauss quadrature is used to evaluate the integrations in Equations (16) and (19). The integrations are performed based on the local quadrature domains centered by field nodes. Obviously, more reasonable results can be obtained using the weight function with the local property, which decreases in magnitude as the distance from a quadrature point x_Q to a field node x_i increase. The following quartic spline weight functions only depend on the distance between two points is used in the present work:

where d_i =|x_Q -x_i | is the distance between node x_i to the quadrature point x_Q and r_w is the size of the support for the weight function.

For each Gauss quadrature point x_Q , the LRPIM shape functions are constructed to obtain the integrand. For field node x_i , there exist three local domains, which are the local quadrature domain Ω_q (size r_q ), the local weight function domain Ω_w , where w_i ≠0 (size r_w ), and the local support domain Ω_s (size r_s ). These domains are arbitrary as long as the condition r_q ≤r_w is satisfied. As the quartic spline weight function will be zero along the internal boundary Γ_qi when the sizes of local integration domain and weight domain are the same, r_q =r_w is adopted in the present work. r_q and r_s are defined as r_q =α_qd_ci and r_s =α_sd_ci , respectively, where α_q and α_s are dimensionless sizes chosen to control the actual domain sizes and d_ci is the shortest distance between the node x_i and neighbor nodes.

Note that the LRPIM shape functions have a complex feature and different forms in each small integration region; therefore, the derivatives of shape functions might have an oscillation. In addition, the sum of all local quadrature domains should cover the whole domain. The overlapping of interpolation domains makes the integrand in the overlapping domain very complicated. For the reasons above, the local quadrature domain Ω_q should be divided into some regular smaller partitions to guarantee the accuracy of the numerical integration [52].

3.1 Governing equations of upper and lower sublaminates

By assembling the equivalent stiffness matrixes and equivalent external load vectors of local quadrature domain, the LRPIM formula of Hamilton canonical equation for a layer with n nodes can be expressed by

where P(z)=[σxz1 σxz2 ⋯ σxzn σyz1 σyz2 ⋯ σyzn σzz1 σzz2 ⋯ σzzn]T denotes the transverse stress vector, Q(z)=[u1 u2 ⋯ un v1 v2 ⋯ vn w1 w2 ⋯ wn]T denotes the displacement vector, [P(z) Q(z)]^T is the state variables, A12=A12T and A21=A21T are the equivalent stiffness matrices for the whole domain, and Fp(z)=[fp1 fp2 ⋯ fpn]T and Fq(z)=[fq1 fq2 ⋯ fqn]T are the equivalent external load vector for the whole domain.

The exact solution to Equation (23) of the μ-th sublaminate can be expressed by

where

where P_μ (h_μ ) and Q_μ (h_μ ) are the stress and displacement vectors, respectively, relative to the bottom surface of the μ-th sublaminate. P_μ (0) and Q_μ (0) are the stress and displacement vectors, respectively, relative to the top surface of the μ-th sublaminate. h_μ is the thickness of the μ-th sublaminate.

Based on the transfer matrix technique, which uses interlaminar displacement and stress compatibility conditions, i.e. the state variables of the two surfaces are regarded as the same, the recurrence formulation of a laminate with l-layer is

where h_k (k=1, 2, …, l) is the thickness of each layer.

Equation (26) is the governing equation of a l-layer composite laminates. It can be recast into a matrix form as follows:

where [P(h_l ) Q(h_l )]^T and [P(0) Q(0)]^T are the state variables in the bottom and top surface of the laminate with l layer. [T11T12T21T22] is known as the equivalent stiffness matrix, and [Γp Γq]T=(∏k=2lTk)F1(h1)+(∏k=3lTk)F2(h2)+⋯+Fl(hl) is the equivalent external load vector.

Applying the procedure above to the upper and lower sublaminates, respectively, two control equations that are similar to Equation (27) can be obtained. For convenience, the following expressions with superscripts are used to express the control equations of upper and lower sublaminates. For the upper sublaminates,

and for the lower sublaminates,

where superscripts bu and tu denote the bottom and top surfaces of upper sublaminates, respectively, and superscripts bl and tl denote the bottom and top surfaces of lower sublaminates, respectively.

3.2 Spring layer model of interfaces with delaminations

The spring layer model, which is similar to that used in [44–47], is employed to describe the bonding conditions between upper and lower sublaminates:

where σxztl, σyztl, σzztl, σxzbu, σyzbu, and σzzbu denote the out-of-plane stresses of top surface of lower sublaminate and bottom surface of upper sublaminate. u^tl, v^tl, w^tl, u^bu, v^bu, and w^bu denote the displacements of top surface of lower sublaminate and the bottom surface of upper sublaminate. K_x , K_y , and K_z are three bonding stiffness constants in the x, y, and z directions, respectively. Obviously, the displacements will be continuous across the interface when K_i (i=x, y, z)→∞, which implies a perfect bonding, whereas K_i (i=x, y, z)=0 indicates that the upper sublaminates and the lower sublaminates are completely delaminated from each other.

To prevent interpenetration phenomenon between delaminated sublaminates in the delaminated region, a unilateral frictionless contact interface characterized by a zero stiffness for opening relative displacements (Δw≥0) and a positive stiffness for closing relative displacements (Δw<0) can be introduced, which can be expressed as follows:

where sign is the signum function.

Assuming that the FEM background cells scheme used in spring layer with n field nodes is the same as that used in each layer of sublaminates, the matrix form of Equation (30) can be expressed as

where Ptl=[σxz1tl σxz2tl ⋯ σxzntl σyz1tl σyz2tl ⋯ σyzntl σzz1tl σzz2tl ⋯ σzzntl]T denotes the transverse stress vector of the top surface of lower sublaminate, Pbu=[σxz1bu σxz2bu ⋯ σxznbu σyz1bu σyz2bu ⋯ σyznbu σzz1bu σzz2bu ⋯ σzznbu]T denotes the transverse stress vector of the bottom surface of upper sublaminate, Qtl=[u1tl u2tl ⋯ untl ν1tl ν2tl ⋯ νntl w1tl w2tl ⋯ wntl]T denotes the displacement vector of the top surface of lower sublaminate, Qbu=[u1bu u2bu ⋯ unbu ν1bu ν2bu ⋯ νnbu w1bu w2bu ⋯ wnbu]T denotes the displacement vector of the bottom surface of upper sublaminate, I=Diag(1)_3n×3_n , R=Diag(Rx1Rx2⋯RxnRy1Ry2⋯RynRz1Rz2⋯Rzn),Rik=c44(1)/[Kikh] (i=x, y, z; k=1, 2, ⋯, n) are the dimensionless compliance coefficients of the interfaces, where c44(1)=G23(1)=E2(1)/E2(1)2(1+ν23(1)) is an elastic coefficient of the first lamina, G23(1), E2(1), and ν23(1) are the shear modulus, elastic modulus, and Poisson ratio of the first lamina, respectively. In the present computing program, the following value of dimensionless compliance coefficients can avoid the numerical instabilities and acquire reliable numerical results. In the delaminated region, R_x =R_y =R_z =10⁸, whereas in the undelaminated region, R_x =R_y =R_z =0.

3.3 Global governing equation for mechanical and thermal buckling analysis of the laminate with delamination

Substituting Equation (32) into Equation (29), the following expression can be obtained:

Eliminating the vector [P^buQ^bu]^T by Equation (28), the final governing equation can be obtained as follows:

where

is the equivalent stiffness matrix of the laminate with delaminations.

where Θ_p and Θ_q are the equivalent body force vectors.

For the mechanical and thermal buckling analysis of delaminated laminated plates, the equivalent body force vectors vanish from Equation (34), and then the final governing equation can be expressed as follows:

Considering that the top surface and bottom surface are stress free, i.e. the stress vector P^bl=P^tu=0. The governing equation for buckling analysis can be deduced from Equation (37), which can be expressed as follows:

The determinant condition of Q^tu=0 existing nonzero solution is that B₁₂ is a singular matrix. Thus, the eigenpolynomial must be zero, namely

The mechanical and thermal buckling factors λ_m and λ_t can be solved from Equation (39). The critical mechanical and thermal buckling factors λ_mcr and λ_tcr , respectively, coincides with the smallest values of λ_m and λ_t . Note that the elements in matrix B₁₂ are the transcendental functions of the mechanical and thermal buckling factors. This leads to a nonlinear eigenvalue problem, and iteration methods should be used to solve this problem. In this present work, the Newton-Raphson method is used to obtain the required eigenvalues.

Moreover, it is noted that the present mathematical model are implemented by means of programming based on the Mathematica 8.0.1 (Stephen Wolfram, Champaign, Illinois, USA), and it is a original program written for analyzing the mechanical and thermal buckling analysis of delaminated laminated plates.

4.1 Mechanical buckling analysis of an intact laminated plate

Simply supported square symmetric cross-ply (0/90/0…) plates subjected to an in-plane uniform uniaxial compressive loading p_bx are considered for numerical illustration. The plane geometric sketch of the plate is shown in Figure 2. Each layer is a unidirectional fiber reinforced composite and has the same thickness.

The aspect ratio is a/h=b/h=10 and h is the thickness of the plate. The material engineering elastic constants of each layer are

Figure 2:

Plane geometric sketch of delaminated plate under uniaxial load.

The critical buckling load is dimensionless and presented as

where σ_kr is the critical buckling load, which is the product of the critical buckling factor λ_mcr and the in-plane uniform pressure p_bx .

In this example, the whole plate in the xy plane is discretized by 289 field nodes (17×17) distributed regularly and uniformly for convenient in programming and calculation. For Gauss quadrature, the 8×8 quadrature scheme (i.e. 64 Gauss points in each regular smaller partition of local quadrature domain) is employed to evaluate local domain integrals. The dimensionless critical buckling loads of the intact simply supported square laminates with different layer number, and E₁/E₂ ratios are compared in Table 1. The present results are in good agreement with the existing results based on 3D elasticity theory [1], FEM [3], CLT [8], HSDT [8], and LWT [38]. It can be observed that the critical buckling loads increase with the increase in E₁/E₂ ratios and layer numbers, which means the plate stiffness and buckling strength increase with the increase in E₁/E₂ ratios and layer numbers.

4.2 Mechanical buckling analysis of laminated plate with a delamination

4.2.1 Mechanical buckling analyses of delaminated plates under uniaxial compressive load

The simply supported five-layer square plates studied in Section 4.1 with central embedded square delaminations subjected to an in-plane uniform uniaxial compressive load are considered for a comparison purpose. The aspect ratio and material engineering elastic constants of each layer are the same as that described in Section 4.1. The effect of delamination locations and delaminated areas on the buckling behavior of laminated plates is investigated. The dimensionless critical buckling load is defined the same as that given in Section 4.1.

The results of dimensionless critical buckling loads of delaminated plates with different delamination locations and delaminated areas are presented in Figure 3. Interface 1 represents (0/90/0//90/0) delamination (where // denotes the delamination location) and interface 2 represents (0/90/0/90//0) delamination. Figure 3 shows that the present results are in good agreement with the existing results based on LWT [38]. It can easily be observed that the dimensionless critical buckling loads decrease with the increase in the delaminated areas, especially when delaminated area exceeds 6%. In addition, the downtrend of dimensionless critical buckling loads is most obvious when the delamination lies in midplane, which means that the buckling strength of delaminated plate decreases even more when the delamination is closer to midplane.

Figure 3:

Dimensionless critical buckling load for different delamination locations and areas of plates.

4.3 Thermal buckling analysis of an intact laminated plate

Simply supported antisymmetric angle-ply (±45)₃ square plates subjected to uniform temperature rise ΔT are considered for numerical illustration. The material engineering elastic constants of each layer and thermal expansion coefficients are

The following dimensionless buckling temperature has been used in this study:

In this section, the discrete method and quadrature scheme are the same as that used in the last two sections. The results of the critical buckling temperature T_cr of (±45)₃ simply supported square laminates with different side-to-thickness ratios are presented in Table 3. Good agreement can be observed between the present results and that obtained based on CLT [7] (except when a/h≤20), FSDT [5], HSDT [2], and LWT [6]. The exception in results obtained by CLT for thicker plates occurs because the CLT cannot account for the effect of interlaminar shear deformation. The critical buckling temperatures decrease with the increase in side-to-thickness ratios, and the discrepancy of the results is decreased as the side to thickness ratio is increased.

4.4 Thermal buckling analysis of laminated plate with a delamination

Symmetric angle-ply (0/±45/0₂/±45/90₂/±45)_s rectangular plates with central elliptical and rectangular delaminations subjected to uniform temperature rise ΔT are considered for a comparison purpose. The dimensions of the delaminated laminated plates are shown in Figure 4, which can also refer to Wang et al. [27].

Figure 4:

Dimensions of the laminated plate with elliptical or rectangular delaminations (mm).

The major axis of elliptical delamination and the length of rectangular delaminations are a₁; the minor axis of elliptical delamination and the width of rectangular delaminations are b₁; the delaminations lies between the second layer and the third layer in the laminated plates. The two edges (x=0, a) perpendicular to x axis are clamped support, and the other two edges (y=0, b) are free. The typical material properties of the plates and thermal expansion coefficients are E₁=131 Gpa, E₂=13 Gpa, G₁₂=6.4 Gpa, ν₁₂=0.34, α₁=1×10^-6/°C, α₂=22×10^-6/°C.

In this section, the discrete method and quadrature scheme are the same as that used in the above sections. The results of critical buckling temperature ΔT_cr of the symmetric angle-ply rectangular laminated plates with elliptical or rectangular delaminations are presented in Table 4.

Good agreement can be observed between the present results and that obtained by a Rayleigh-Ritz analysis based on the Whitcomb delaminated buckling model [27]. As some uncontrollable external factors can never be avoided, a maximum error of <20% between the experimental results and the theoretical solutions is acceptable. Although the area of elliptical delamination is smaller than that of rectangular delamination, the critical buckling temperatures of plates with elliptical delaminations are lower than that of plates with rectangular delaminations. Therefore, it can be inferred that the thermal buckling behavior of laminated plates with delaminations is sensitive to the shape of delaminations.