Inelastic damage using continuum damage mechanics in composite plate reinforced by unidirectional fibers

It is well that a finite element method is very popular simulation method to predict the physical behavior of systems and structures. In the last years an increase of interest in a new type of numerical methods known as meshless methods was observed. The paper deals with application of radial basis functions on modelling of inelastic damage using continuum damage mechanics of layered plate composite structures reinforced with long unidirectional fibers. For numerical simulations of elastic-plastic damage of layered composite plates own computational programs were implemented in MATLAB programming language. We will use the Newton-Raphson method to solve nonlinear systems of equations. Evaluation damage during plasticity has been solved using return mapping algorithm. The results of elastic-plastic damage analysis of composite plate with unsymmetrical laminate stacking sequence are presented.


Introduction
Carbon fiber-reinforced composite materials are finding increased application in many areas [1].The finite element method (FEM) is one of the most widely used and most popular numerical methods for damage of plate structures [2,3].In the last years an increase of interest in a new type of numerical methods known as meshless methods or meshfree was observed [4].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.One of the areas where meshless methods are convenient to use is the analysis of plate and shell structures [5].
This paper deals with the application of radial basis function (RBF) on modelling of inelastic damage using continuum damage mechanics (CDM) of layered plate composite structures reinforced with long unidirectional fibers.For numerical simulations of elastoplastic damage of layered composite plates own computational algorithms were implemented into the MATLAB programming language.For solving a system of nonlinear system equations the Newton-Raphson method was used.The evaluation of damage during plasticity has been solved using return mapping algorithm described in [6].The results of elasto-plastic damage of a composite plate with unsymmetrical laminate stacking sequence are also presented.

Point interpolation methods
One of the common characteristics of all mesh-free methods is their ability to construct functional approximation or interpolation entirely from information at a set of scattered nodes or particles, among which there is no pre-specified connectivity or relationships.The point interpolation method (PIM) was proposed by Liu [7] to replace moving least square (MLS) approximation for creating shape functions in meshfree settings.In contrast to mesh-based methods, algorithms based on collocation methods are frequently called meshless or meshfree methods.The terminology "meshless" is not always clear, and it is sometimes also used for Galerkin or Petrov-Galerkin methods [8].In both cases, a variational formulation of a physical problem is solved in an average sense around node points.However, these types of methods still require geometrical information about the connectivity of the nodes, i.e. the volume surrounding a node, and are therefore not true meshless methods.
The goal of this section is to describe mathematical formulation of computational models based on (RBF) for linear analysis of layered composite structures.

Radial point interpolation method
The RBF was originally introduced in the 1970s to multivariate scattered data approximation and function interpolation.Many types of meshless RBF methods have been proposed in the literature [9,10].Radial point interpolation method (RPIM) can be included among the methods that represent the functions of unknown field variables by series.This method, in addition to the least squares (MLS) method, belongs to the most commonly used interpolation / approximation techniques in meshless methods.
The function defined on the local sub-area is interpolated at point [ ] as follows where denotes the number of nodes supporting interpolation at a given point.These points can be selected using the concept of a local support domain or concept of the domain of influence [8].Vector is is vector of unknown functions and is radial basis function with distance defined as The most widely used RBF in the RPIM are the functions listed in Table 1 [7], where is the average distance between nodes in the local sub/area.In these functions, there are several shape parameters that must be determined prior to performing the analysis.One of the advantages of the shape functions assembled using the RPIM is the property of the Kronecker delta function, which implies that no special techniques are required to introduce the essential boundary conditions.The RPIM which uses only RBF is not polynomial consistent because cannot reconstruct the polynomial fields accurately.For this reason, polynomial members are added to the base functions.The function will be interpolated as follows where is coefficient for radial base and is coefficient for polynomial base .The number of RBF n is determined by the number of nodes in the support region and the number of polynomial bases m can be selected based on the reproducibility requirement [10].
The coefficients and can be determined from the condition that interpolation passes through all nodes, which can be written in the matrix form as , where is vector of function values of the field variable in n nodes.The moment matrix corresponding to RBF is given as where is the distance independent of the direction, so it must be meet that ( ) .Then it is valid that matrix and also matrix (defined below) are symmetric.Note that matrix is invertible for arbitrary scattered nodes [11].The moment matrix is given as: Polynomial terms must meet extra requirements to guarantee a unique approximation of function.The following restrictions are often used [4] ∑ .
Combining Eq. ( 4) and Eq. ( 8) yields the following set of equations in the matrix form or , -, -.
By solving the above system of equations we get the following relation for interpolation in point where matrices and are defined as follows From Eq. 11 the matrix of shape functions is Shape function for -th node is given as The first and the second derivative of shape function we get differentiation (15) In Fig. 1 to Fig. 3 is described the shape function as well as its first and second derivations.These functions are assembled for a central node in a square area represented by 9x9 nodes, The domain of influence concept was used with a circular shape for construction of the shape function.A multi-quadratic RBF was used to which linear polynomial members were added, i.e. m = 3.The average node spacing was calculated as the average of the distances to the four closest nodes.The following parameters were used to construct shape functions: =4,0, =4.0, =1.03.

Governing equations and discretization
Reddy's theory of shear deformations is one of many higher order theories.Higher-order theory can better describe deformation kinematics, do not require shear correction coefficients, and more accurately describe the distribution of transverse shear stresses.These theories, however, include internal higher order power quantities that are difficult to interpret physically and require higher calculation costs.In principle, it is possible to extend the displacement field expressed by thickness dimensions to any degree.Due to algebraic complexity and higher computational costs, theories higher than the third order are not used.The reason for which the relations for displacement were extended to the cubic member is to achieve the quadratic distribution of the transverse shear deformations and the stresses across the thickness of each layer.Therefore, shear correction coefficients are not required in these theories.In Reddy´s theory, the displacement field is expressed by the unknown , , , and as follows . (20) The strain components is obtained by the equations of linear elasticity where Integrating the stress over the thickness of the plate we can determine the stress resultants where the symbols and can take and , N denotes the normal and in-plane shear forces, M means the bending and twisting moments, Q is the shear force, P and S are higher-order stress resultants, respectively.If we use interpolation, then we can discretize equations (26) By placing them in the Reddy governing equations [12] we obtain for all nodes a system of linear equations The shape functions (14) assembled using RBF possess the delta function property, so the imposition of essential boundary conditions can be applied similarly as in FEM.
The material model was verified for T300/914 Carbon/Epoxy [13].The elastic properties of T300/914 lamina arte given in Table 2.The strength and critical integrity values are shown in Table 3.The indices t and c at critical values of the damage variable denote tensile stress, and compression stress.The damaged in-plane shear moduli at imminent failure is = 3.410 GPa.The plastic strain and damage thresholds are = = 17.8 MPa.The parameters required to define the damage surface and the plasticity surface as well as the parameters needed to define the hardening functions were calculated according to the procedure described in [6].Only two terms were taken in Prony's series in calculation of the parameters.These parameters are listed in Tab. 4 -6.Subsequently, the material model with these input parameters was used to determine the material response under a shear load in plane 1-2.This response was tested at one point with the prescribed values.The results obtained using the programmed material model were compared with the results obtained from experimental data published in the literature [13].

Damage and material model for elasto-plastic damage
For polymer matrix composites reinforced by strong and stiff fibers, damage and its conjugate thermodynamic force was used.The formulation of the damage and material model is based on linking CDM with the classical plasticity theory by using thermodynamic principles of irreversible processes.This model is suitable for damage prediction and plastic deformations in composite materials reinforced by long uniaxial fibers with a polymer matrix.To represent the damage, the diagonal damage tensor D was used.Experimental data for T300/914 Carbon Epoxy are used to identify model.The development of damage together with the development of plastic deformation is solved incrementally and is evaluated at each collocation point.The input and output data for the material model are described in Fig. 4, where input variables are: ̅ is the constitutive matrix in the effective configuration, , , , a are parameters for damage and plasticity surface, , , , are parameters for hardening functions, , are damage threshold and yield threshold, are damage components, are components of plastic strain, p -is yield hardening variable and is damage hardening variable.The output variables are: D i are damage variables, stress components and tangent constitutive matrix at the end of the last iteration.

Numerical results and discussion
As a verification example, a composite plate with the dimensions 100x100mm with a thickness of 2.5mm and ply stack orientation [0/45/-45/90] is selected.The plate is clamped on the edges and subjected to pressure load q* = 1MPa.The material data for each ply are: = 142 GPa, = 10.3GPa, = 6.42 GPa, = 3.71 GPa, = 0.21.The implemented computational models were verified by comparing the results with results obtained from the FEM program ANSYS with a very fine mesh (40 000 quadratic shell91type shell elements).
Unlike metal structures, when evaluating the results of the analysis of the composite structure it is necessary to evaluate quantities such as stresses or deformations in the coordinate system of the lamina in the individual locations of structure.This is because a single lamina of fiber reinforced composite behaves as an orthotropic material.The degree of stiffness reduction can be calculated from the values of the individual variables of damage and integrities.From the evaluation of the values of the damage variable the degree of stiffness reduction in the direction of the fibers can be deduced (Fig. 5).The degree of stiffness reduction in the direction perpendicular to the fibers in plane 1-2 can be deduced from the values of the damage variable (Fig. 6).To evaluate the degree of stiffness reduction caused by the shear load in plane 1-2, the value of the integrity component must be evaluated (Fig. 7).In the case of the analysed composite plate with an asymmetrical layout, the most critical places from the point of view of the macro cracks caused by reduction of stiffness in the fiber direction are in the vicinity of the center of the edges are parallel to the x-axis, especially in the [45 °] layer.The most critical macro-crack spots due to the reduction of stiffness in the direction perpendicular to the fibers in the laminate plane are locations in the center of the edges parallel to the x-axis in the [0 °] layer.

Conclusion
In the presented paper a developed algorithm for the analysis of elasto-plastic damage by using RPIM based on strong form of governing equations based on Reddy's plate theories is presented.The shape functions constructed by RPIM have the Kronecker delta function property, thus there is no need to apply special techniques to apply essential boundary conditions.From the RPIM results we can see that implemented computational model be able to calculate the nodal displacements with very good accuracy using a relatively small number of nodes.The stresses , and are also calculated with good accuracy.However, when the transverse shear stress is compared with FEM calculations the errors are larger.In the case of the analysed composite plate with non-symmetrical layer arrangement, the maximum damage is in the fiber direction around the middle point parallel to the x-axis, especially in the layer [45 °].In the direction perpendicular to the fiber, the maximum damage is in the area of the middle point of the x-axis parallel to the layer [0 °].The maximum damage in shear plane 1-2 is in the central area of the edges parallel to the xaxis in the layer [0 °].For a given load of the analysed composite plate only very small plastic deformations were observed.

Fig. 4 .
Fig. 4. Description of input and output data for elasto-plastic material model

Table 1 .
Typical RBF with dimensionless shape parameter used in RPIM

Table 3 .
Strength and critical integrity properties for a T300/914 unidirectional lamina

Table 4 .
Damage parameters for defining surface damage for a T300/914 unidirectional lamina

Table 5 .
Plasticity parameters for defining plasticity surface for a T300/914 unidirectional lamina