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

**Milan Žmindák**^{1,}, **Daniel Riecky**^{2}, **Zoran Pelagić**^{1}, **Martin Dudinsky**^{1}

^{1}Department of Applied Mechanics, University of Žilina, Žilina, Slovak Republicy

^{2}Plastic Omnium Auto Exteriors, Ltd., Slovak Republic

### Abstract

This paper deals with the application of meshless methods for the analysis of composite plates. The main attention is focused on the implementation of the Meshless Local Petrov Galerkin (MLPG) formulation for multilayered orthotropic plates. At first for this purpose the implementation of homogenization theory was needed and analyzes were made to obtain homogenized material properties of composite plates. The software for homogenization of material properties uses direct homogenization method that is based on volume average of stresses on the representative volume element (RVE). Homogenization is performed by a multi software approach, by linking MATLAB and ANSYS software. Then data obtained are used in analyzes performed in user own software, which is based on the MLPG method. Strain, stress and displacement fields were analyzed. Results obtained by MLPG were compared with those obtained by FEM programs, ANSYS and ABAQUS.

### At a glance: Figures

**Keywords:** composite plates, homogenization, meshless, local Petrov-Galerkin method

*American Journal of Mechanical Engineering*, 2013 1 (7),
pp 461-469.

DOI: 10.12691/ajme-1-7-62

Received October 18, 2013; Revised October 31, 2013; Accepted November 11, 2013

**Copyright**© 2014 Science and Education Publishing. All Rights Reserved.

### Cite this article:

- Žmindák, Milan, et al. "Meshless Local Petrov-Galerkin Formulation for Static Analysis of Composite Plates Reinforced by Unidirectional Fibers."
*American Journal of Mechanical Engineering*1.7 (2013): 461-469.

- Žmindák, M. , Riecky, D. , Pelagić, Z. , & Dudinsky, M. (2013). Meshless Local Petrov-Galerkin Formulation for Static Analysis of Composite Plates Reinforced by Unidirectional Fibers.
*American Journal of Mechanical Engineering*,*1*(7), 461-469.

- Žmindák, Milan, Daniel Riecky, Zoran Pelagić, and Martin Dudinsky. "Meshless Local Petrov-Galerkin Formulation for Static Analysis of Composite Plates Reinforced by Unidirectional Fibers."
*American Journal of Mechanical Engineering*1, no. 7 (2013): 461-469.

Import into BibTeX | Import into EndNote | Import into RefMan | Import into RefWorks |

### 1. Introduction

Plate and shell structures are the most widely used structural members in mechanical and civil engineering thanks to their good weight to load carrying capability ratio. Pipes, aeroplanes, hulls of ships ^{[1]}, bridge reinforcements as well as reinforced roofs and membranes are good examples.

Thanks to good material properties, especially to good weight to strength ratio fiber-reinforced materials came into practice. Laminates are composed of multiple layers where fiber orientation in layers can be different ^{[2]}. By this way it is possible to obtain required material properties in required direction.

The finite element method (FEM) is one of the most widely used and most popular numerical methods for analyzing plate structures ^{[3]}. Although the method is stable, well developed and has reached extensive development during last decades, it also has some limitations. One of the assumptions in FEM is an existence of the mesh created from the finite elements. Mesh quality affects the accuracy of obtained solution and stability of convergence. Creating well defined meshes can be very time consuming in some cases. Generally it is recommended to use finite elements with shapes which are the most similar to ideal shapes ^{[4, 5]}. Degenerated geometry or deformed shape can have a negative effect on the solution accuracy. This is very important in analyzing nonlinear problems with large deflections or large strains such as metal forming, fragmentation after impact, etc. In these cases an improper mesh causes serious decrease of accuracy or failure of computation.

In last years an increase of interest in new type of numerical methods known as meshless methods was observed ^{[6, 7, 8, 9, 10]}. These methods are interesting due to their flexibility and ability of solving boundary value problems without predefined mesh. Computational model in these methods is represented by set of nodes distributed within global domain and its boundary. These nodes do not have to be connected into explicitly defined elements. Information about relations between the nodes is theoretically not necessary before analysis. By this way some problems with the mesh in FEM can be overpassed. Problems with remeshing can be overpassed by adding or removing some nodes as necessary. In some methods the resultant stress fields are globally continuous and this simplifies consequent analyses.

One of the areas where meshless methods are convenient to use is analysis of plate and shell structures. These methods are useful due to flexibility of meshless algorithms and ability of meshless approximation functions to obtain interpolating field with high order of continuity by simple way. In some cases it is possible to overcome some “locking” effects in more simple way than in FEM. Meshless methods are relatively new concept in computational mechanics. Compared to FEM formulations there are less meshless formulations available for plate and shell structures. Additional research and development of general meshless methods able to successfully solve various problems in plate structural analyzes are therefore necessary. Meshless approaches for problems of continuum mechanics have attracted much attention during the past decade especially owing to their high adaptivity and low costs to prepare input data for numerical analyses.

### 2. Composite Homogenization

Plates can be reinforced by long glass, carbon or kevlar fibres. The first step is to obtain homogenized lamina material properties. These can be obtained by using homogenization techniques.

**2.1. Homogenization Techniques**

There are various homogenization methods ^{[11, 12]}. Direct homogenization is based on the volume average of field variables, as stress, strain and energy density. Effective properties can be calculated from effective properties definitions. The diameter and the calculation of field variables can be calculated numerically, for example by FEM, BEM and the geometry and microstructural properties can be generalized.

Indirect homogenization is based on the Eshelby solution of self-deformation for one inclusion in an infinite matrix – equivalent inclusion method ^{[13]}. This method does not use averaging field variables and the effective properties can be obtained by deducing the volume fractions and the inclusion geometry also as the component properties. According to this pattern these methods were developed: self-consistent scheme ^{[14]}, generalized self-consistent scheme Christensen and Lo ^{[15]}, differential methods ^{[16]}, Mori-Tanaka method ^{[17]}. They are often used to calculate various properties of composites. But the generalized microstructural morphology, which is often present in real materials, cannot be handled out deterministically in these models. Constitutive responses of the component phases are limited and the estimated results with large disagreements are not reliable enough. These models cannot for they insufficient representation of microscopic stresses and strains catch the effects of local non-homogeneities.

An alternate approach to direct and indirect homogenization is the variational method, which can determine the upper and lower borders of the elasticity modulus ^{[18]}.

A relatively new approach or homogenization of microstructures consisting of mathematical homogenization based on a two-scale extension of the displacement field. This comes from the analysis of physical systems containing two or more scales ^{[19]}. This approach is good for multiphase materials, in which are natural scales of microscopic scales, characteristic heterogeneity or local discontinuity spacing. Macroscopic scales are characterized by the body dimensions. This method can be called mathematical homogenization.

**2.2. Homogenization Results**

This part describes the procedure of homogenization of material properties of composites of using the method of representative volume element (RVE). For the analysis of the material properties an own software was programmed in MATLAB language and a part of the solution was carried out in ANSYS software. The RVE consists of volume elements and is then loaded by unit strains in various directions. The effective lamina properties are obtained from the volume fractions of stress values obtained by loading the RVE.

**Fig**

**ure**

**1**. Fiber arrangement in matrix: a)hexagonal array, b) square array

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*_{1} is the *x-*direction, in this case the fiber direction, *a*_{2} is the y direction, orthogonal to the fiber direction, *a*_{3} in direction z, transverse vertical to the fiber direction.

### 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*_{1}*, x*_{2}). 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*_{3}.

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*_{1}, *x*_{2}]^{T} is the position vector, *w*_{α }( *x*_{1}, *x*_{2}*,* *t*) are plane rotations about the axes for indices *α = 1,2* and *w*_{3}(*x*_{1}, *x*_{2}, *t*) is the deflection of the *z *plane (*x*_{1}, *x*_{2}).

**Fig**

**ure**

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

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.

**Fig**

**ure**

**3**. Layer coordinate notation

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*^{3}*/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) |

**Fig**

**ure**

**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

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*_{3}(**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*_{1}^{h}*, w*_{2}^{h}*, w*_{3}^{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*_{1}*, x*_{2} 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_{1}(x) ; Q _{2}(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**^{i} and also for the node on the boundary .

The point **x**^{i} located on the global boundary Γ the boundary subdomain is split into *L*^{i}_{Q} and *Γ*^{i}_{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.

**Fig**

**ure**

**5.**Flowchart representing process of calculation

### 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)_{2} / 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.

**Fig**

**ure**

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

**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.

**Fig**

**ure**

**7.**Deformation course through the plates thickness. a) deformation

*ε*

_{11,}

_{ }b) deformation

*ε*

_{22}

*.*

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.

**Fig**

**ure**

**8.**Stress distribution in the middle of each ply M40J vf04 a) σ

_{11, }b) σ

_{22}.

**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.

**Fig**

**ure**

**9.**Deformation course through the plates thickness. a) deformation ε

_{11, }b) deformation ε

_{22}

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.

**Fig**

**ure**

**10.**Stress distribution in the middle of each ply S2Glass vf04 a) σ

_{11. }b) σ

_{22}

**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.

**Fig**

**ure**

**11.**Deformation course through the plates thickness. a) deformation ε

_{11, }b) deformation ε

_{22}

**Fig**

**ure**

**12.**Stress distribution in the middle of each ply K49vf04 a) σ

_{11, }b) σ

_{22}

### 6. Conclusion

This paper was prepared in order to implement and apply "new" numerical methods which are used to solve problems in mechanics, specifically analyzing the stress conditions of composite plate structures. The presented methods are very attractive but the theoretical understanding of these methods and their applications are not as extensive as the knowledge of the finite element method. The meshless method is a new tool that after sufficient debugging of the method parameters (like integration, the integration of size, size of weight function) provides accurate results.

Errors in the values of strain and stress may have several sources. The accuracy of the MLPG method is influenced by several factors like rounding errors when compiling approximation and rounding errors caused by numerical integration. Mentioned errors are reflected in the differences between the reference values and calculated values. Both FEM programs were calculating deflection in the middle layer as the arithmetic average of the deformation, or strain, the lower and upper interface of the plate.

### Acknowledgement

The authors gratefully acknowledge for support the Slovak and Technology Assistance Agency registered under number APVV-0169-07, Slovak Grant Agency VEGA 1/1226/12.

### References

[1] | Davies, P., Rajapakse, Yapa, D.S., Durability of Composites in Marine Environment, Springer Science + Business Media Dordrecht, 2014. | ||

In article | |||

[2] | Piovar, S., Kormanikova, E., Statical and Dynamical Analysis of Composite Sandwich Plates, Bulletin of the Transilvania University of Brasov, Series I: Engineering Sciences, 4 (53), No.1, 2011. | ||

In article | |||

[3] | Zmindak, M., Dudinsky, M., Finite Element Implementation of Failure and Damage Simulation in Composite Plates, In: HU, N., eds: Composites and Their Properties, InTech, Rijeka, 2012, 131-152. | ||

In article | |||

[4] | Liu, G.R., Quek,S.S., The Finite Element Method: Practical Course, Elsevier Science, Ltd., 2003. | ||

In article | |||

[5] | Saga, M., Kopas, P., Vasko, M., “Some Computational Aspects of Vehicle Shell Frames Optimization Subjected to Fatigue Life” Communications, 12 (4), 2010, p. 73-79. | ||

In article | |||

[6] | Atluri, S. N. , The Meshless Method, (MLPG) For Domain & BIE Discretizations, Tech Science Press, 2004. | ||

In article | |||

[7] | Sladek, J., Sladek, V., Krivacek, J., Wen, P.; Zhang, Ch., “Meshless Local Petrov-Galerkin (MLPG) Method for Reissner-Mindlin Plates under Dynamic Load” Computer Meth. Appl. Mech. Engn., 196, 2007, 2681-2691. | ||

In article | CrossRef | ||

[8] | Sladek, J., Sladek, V., Atluri, S. N., “Application of the Local Boundary Integral Equation Method to Boundary-value Poblems.,Int. Appl. Mech., 38 (2), 2002, 1025-1047. | ||

In article | CrossRef | ||

[9] | Soares, D. Jr., Sladek, J., Sladek, V., Zmindak, M., Medvecky, S : “Porous Media Analysis by Modified MLPG Formulations”, CMC: Computers, Materials,& Continua, 27(2), 2012, pp. 101-127. | ||

In article | |||

[10] | Zmindak, M., Riecky, D., Soukup, J.,”Failure of Composites with Short Fibres”. Communications, 12 ( 4), 2010, 33-39. | ||

In article | |||

[11] | Quin, Q.H., Yang., Q,S., Macro-Micro Theory on Multifield Coupling Behavior of Heterogenous Materials, Springer, 2008. | ||

In article | |||

[12] | M. Sejnoha, J. Zeman, Micromechanics in Practice, WIT Press, Southampton, UK, 2013. | ||

In article | |||

[13] | Eshelby, J. D., The Determination of Elastic Field of an Ellipsoidal Inclusion and Related Problems, In:Proc. R. Soc. Londýn, 1957 276-396. | ||

In article | |||

[14] | Hill, R., “A self-consistent Mechanics of Composite Materials”, J. Mech. Phys. Solids, 13, 1965, 213-222. | ||

In article | CrossRef | ||

[15] | Christensen, R. M., Lo. K. H. (1979), “Solutions for Effective Shear Properties in Three Phase Sphere and Cylinder Models”, J. Mech. Phys. Solids, 27, 315-330. | ||

In article | CrossRef | ||

[16] | Norris, A. N., “A Differential Scheme for the Effective Moduli of Composites”, Mechanics of Materials, 1985, 1-16. | ||

In article | CrossRef | ||

[17] | Mori, T.- Tanaka, K., “Average Stress in Matrix and Average Elastic Energy of Materials with Misfitting Inclusion”. Acta Mat., 21, 1973, 571-574. | ||

In article | CrossRef | ||

[18] | Hashin, Z., Shtrikman, S., “On Some Variational Principles in Anisotropic and Nonhomogeneous Elasticity”, J. Mech. Phys. Solids, vol. 10, 1962, 335-342. | ||

In article | CrossRef | ||

[19] | Bensoussan, A., Lion, J. L., Papanicolaou, G., Asymptotic Analysis for Periodic Structures, Amsterdam Holland, 1978 | ||

In article | |||

[20] | Wallenberger, F. T., Bingham, P. A., Fiberglass and Glass Technology, Springer, 2000. | ||

In article | |||

[21] | Reddy J.N., Mechanics of Laminated Composite Plates: Theory and Analysis, CRC Press, Boca Raton, 1997. | ||

In article | |||