Meshless Local Petrov-Galerkin Formulation for Static Analysis of Composite Plates Reinforced by Unidirectional Fibers

Veiw figure View Figure

Homogenized lamina RVE consists of a fiber and epoxy matrix. The fibers are from three material types, carbon, glass, polyaramid. We assumed cylindrical fiber shapes and an ideal cohesion between the fiber and the matrix. Used carbon fibers are industrially labeled as T300 and M40J. The glass fiber label is EGlass and S2Glass. Polyaramide fibers have the label K49. Fiber material properties are listed in Table 1. and the matrix properties are listed in Table 1 ^[20].

The RVE dimensions are calculated for the hexagonal fiber configuration in Figure 1a from the relations (1) and for a square configuration the RVE dimensions are in Figure 1b, calculated from the relations in (2).

(1)

(2)

Where a₁ is the x-direction, in this case the fiber direction, a₂ is the y direction, orthogonal to the fiber direction, a₃ in direction z, transverse vertical to the fiber direction.

Table 1. Fiber material properties

Download as

Veiw figure View Table

Table 2. Matrix material properties

Download as

Veiw figure View Table

Table 3. RVE Material properties for volume fraction Vf=0.4

Download as

Veiw figure View Table

3. Governing Equations

Classical laminated plate theory was formulated by deriving classical plate theory for composite plates. In this theory the plate composes of N orthotropic layers with total thickness h. The midsurface of layered plate is located in the region Ω, in plane (x₁, x₂). Axis x3 ≡ z is perpendicular to the midsurface Figure 2. k-th layer is located between coordinates from z = z_k to z = z_k+1 in thickness direction x₃.

Deformation of the plate is described by the “Reissner-Mindlin” plate theory ^[21]. In this theory the shear strains in thickness direction are constant and correction coefficients are necessary for calculation of transverse shear forces. Spatial displacement field caused by transverse loading can be expressed, in terms of for displacement components in the form of

(3)

Where x =[ x₁, x₂]^T is the position vector, w_α( x₁, x₂, t) are plane rotations about the axes for indices α = 1,2 and w₃(x₁, x₂, t) is the deflection of the z plane (x₁, x₂).

Figure 2. Composite plate a) geometry and displacements, b) moments and shear forces

Download as

Veiw figure View Figure

Linear deformation can be expressed

(4)

If k-th is lamina made of an orthotropic material, then the relation between stresses σ_ij and strains ε_mn is expressed by the constitutive equation for stress tensor

(5)

With assumption of homogeneous coefficients of constitutive tensor for k-th lamina.

From the equation (4) it can be seen that strains are continuous through the plate thickness. Discontinuous material coefficients cause the stiffness change in the interfaces and hence the stresses in lamina interfaces are discontinuous.

Equation (5) for plate problems is usually written as a tensor of elastic constants of the second order. Constitutive equations for orthotropic material and plane stress have then the form of

(6)

Where

Where , is Young's elasticity modulus in xα (α = 1, 2) -axis direction, , and are shear elasticity moduli and are Poisson's ratios of k-th lamina.

Figure 3. Layer coordinate notation

Download as

Veiw figure View Figure

Bending moments and can be expressed in integral form as

(7)

(8)

Where κ = 5/6 in the Reisner-Mindlin plate theory, z coordinate is considered as is shown in Figure 2. By substituting (4) and (6) into resultant moments and forces in (7), (8) it is possible to express bending moments M_αβ and shear forces Q_α, α, β=1,2 for orthotropic plate in terms of displacements and rotations. In the case of layer-wise continuous material properties the following relations are obtained

(9)

For indices α a β in (9) Einstein summation rule does not apply and material parameters D_αβ and C_αβ are given by relations

(10)

For homogeneous plate relations (10) are simplified to the form of

(11)

Plate is subjected to transverse distributed loading q(x, t). If each lamina has homogeneous density in thickness direction, equations of motion for Reissner linear theory for thick plates can be written in the form of

(12)

(13)

Where the upper dot denotes the time derivative by t, with the order equal to the dot number, the indexes α, β = 1, 2

(14)

are the global inertial characteristics of the laminate plate. For a plate with a constant density through the whole thickness are these characteristics expressed as I_M=ρh³/12, I_Q= ρh.

In MLPG method the local weak form is assembled on local subdomain , which is a small domain for each node within the global domain ^[6]. These local subdomains overlap each other and cover the whole global domain Ω, Figure 4.

(15)

(16)

Figure 4. Local boundaries for weak formulation, the support domain Ω for MLS approximation of the trial function, and support domain of weight function around node X_i, Ω_Q is integration domain around given node

Download as

Veiw figure View Figure

where w^*_αγ and w^*is the weight or testing function. By applying to the equations (13) and (14) the Gauss divergence theorem we obtain

(17)

(18)

is boundary of local subdomain, M_α(x)=M_αβ(x)n_β(x) is the normal bending moment and n_α is the outward normal on the boundary . Local weak forms (17) and (18) are starting points for deriving local boundary integral equations on the basis of proper test function. Unit step function is chosen in each subdomain for test functions w^*_αγ(x) and w^*(x)

(19)

Then the local weak form (15) and (16) transforms on the following local integral equations

(20)

(21)

where w_α(x;t) is the trial function corresponding to rotations α = 1, 2, the trial function w₃(x; t) corresponds to the transverse displacement. These trial functions are assembled in MLS approximation over nodes in local subdomain around the point x.

4. Numerical Implementation of MLPG

In general, a meshless method uses a local interpolation to represent the trial function with the values (or the fictitious values) of the unknown variable at some randomly located nodes ^[21]. To approximate the distribution of the generalized displacements (rotations and deflection) in over a number of randomly located nodes X_j (j = 1, 2, ... , n), the MLS approximant of is defined by

(22)

Where vector w^h = [w₁^h, w₂^h, w₃^h]^T, ϕ^αis the MLS approximated shape function for the corresponding node.

(23)

Where ϕ^α_,_k(x) is the partial shape function derivative ϕ^α(x) for the k = x₁, x₂ direction.

Substituting equations (22) and (23) into the bending moment definition in equation (9) by using the expression M_α(x; t) = M_αβ (x; t) n_β (x), becomes

(24)

Where vector w^*αis defined as a column vector w^*α=[;]^T, the matrices are related to the normal vector n(x) at the area boundary defined by

(25)

matrix which defines the gradients of shape functions which is written as

(26)

The influence of material orientation for the composite is included in and , which are defined in equations (10), the shear force equations are

(27)

where Q(x) = [Q₁(x) ; Q ₂(x)]^T and

(28)

Substituting the MLS discretized moment and force field (24) and (27) into the local integral equations (20) and (21) local discrete integral equations will be formed

(29)

(30)

where is the prescribed moment on the boundary and

(31)

The equations (29) and (30) are for the subdomains in the neighborhood of the inner node xⁱ and also for the node on the boundary .

The point xⁱ located on the global boundary Γ the boundary subdomain is split into Lⁱ_Q and Γⁱ_QM (this part of the global boundary on which are the bending moments prescribed), Figure 2

In the local weak form (15) and (16) are not associated with penalization parameters or Lagrange multiplier, because the boundary conditions on (part of the global boundary, on which displacements and rotations are prescribed) can be directly formulated, by using the interpolation formulation (22).

(32)

where is the generalized displacement vector described on the boundary. For a clamped plate the first three elements of the vector (displacement and rotations) are zero for the clamped edge and equation (32) is used in all boundary nodes in this case.

For a simply supported plate only the third member of the displacement vector is prescribed the rotations are unknown.

Then equation (29) and the third member of the equation (32) are applied on the nodes laying on the global boundary. On the nodes on the global boundary, on which no boundary condition are prescribed, both local integration equations (29) and (30) are applied.

4.1. Simplification of Formulation for Static Analysis

Equations (29) and (30) were simplified for static loading, the time has no physical meaning, but enables to define the stepwise increasing loading.

(33)

(34)

Then resulting system of equations (29) and (30) has the form

(35)

(36)

which is the local integral equation form, which is solved by the software for plate bending problems. The homogenized data obtained are used in analyzes performed in user own software, which is based on the MLPG method. A simple flowchart representing a process of calculation by MLPG is given in Figure 5.

Figure 5. Flowchart representing process of calculation

Download as

Veiw figure View Figure

5. Numerical Results

In the presented analyzes we considered a composite plate with dimensions L_x and L_y, where L_x = 0,24m and L_y = 0.2m, Figure 6. It consists from six plies, every ply has a thickness Δz = 0.00025m, the total thickness of the plate is h = 0.00115m. The material orientation of the is [45 / 90 / (0)₂ / 90 / -45]. The material properties of analyzed materials are presented in Table 3.

This part presents the results for the deformations and stresses. Both are presented graphically and in tables. The average percentage error (APE) was calculated according to equation

Where N is total number of nodes in given domain, u^ref(x^a) is the reference value in node x^a, u^MLPG(x^a), uMLPG(xa) is value calculated by means of MLPG in node x^a.

A rectangular plate with two types of boundary conditions was considered: a) clamped plate and b) simply supported plate.

Figure 6. Plate dimensions and point distribution, in which the results were compared

Download as

Veiw figure View Figure

5.1. Analysis 1: M40J Fiber

Figure 7 shows the deformation course through the plates thickness in the point of maximum deflection of the plate made from the M40J material with vf04, the results are presented in Table 4.

Figure 7. Deformation course through the plates thickness. a) deformation ε_11,b) deformation ε₂₂.

Download as

Veiw figure View Figure

Table 4. Comparison of deformations ε₂₂and ε₂₂ from MLPG and FEM, M40J vf04

Download as

Veiw figure View Table

Figure 8 shows the stress course in the middle of each ply in the point of maximum deflection of the plate. The results are presented in Table 5.

Figure 8. Stress distribution in the middle of each ply M40J vf04 a) σ_11,b) σ₂₂.

Download as

Veiw figure View Figure

Table 5. Comparison of stresses σ₂₂and σ₂₂ from MLPG and FEM, M40J vf04

Download as

Veiw figure View Table

5.2. Analysis 2: S2Glass Fiber

Figure 9 shows the deformation course through the plates thickness in the point of maximum deflection of the plate made from the S2Glass material with vf04, the results are presented in Table 6.

Figure 9. Deformation course through the plates thickness. a) deformation ε_11,b) deformation ε₂₂

Download as

Veiw figure View Figure

Table 6. Comparison of deformations ε₂₂and ε₂₂ from MLPG and FEM, S2Glass vf04

Download as

Veiw figure View Table

Figure 10 shows the stress course in the middle of each ply in the point of maximum deflection of the plate. The results are presented in Table 7.

Figure 10. Stress distribution in the middle of each ply S2Glass vf04 a) σ_11.b) σ₂₂

Download as

Veiw figure View Figure

Table 7. Comparison of stresses σ₂₂and σ₂₂ from MLPG and FEM, S2Glass vf04

Download as

Veiw figure View Table

5.3. Analysis 3: K49 Fiber

Figure 12 shows the stress course in the middle of each ply in the point of maximum deflection of the plate. The results are presented in Table 9.

Figure 11. Deformation course through the plates thickness. a) deformation ε_11,b) deformation ε₂₂

Download as

Veiw figure View Figure

Table 8. Comparison of deformations ε₂₂and ε₂₂ from MLPG and FEM, K49 vf04

Download as

Veiw figure View Table

Figure 12. Stress distribution in the middle of each ply K49vf04 a) σ_11,b) σ₂₂

Download as

Veiw figure View Figure

Table 9. Comparison of stresses σ₂₂and σ₂₂ from MLPG and FEM, K49 vf04

Download as

Veiw figure View Table