Dynamic simulation of composite layered plates reinforced by unidirectional fibers subjected low velocity impact

In the recent years a big focus is subjected to the response of structures subjected to out-of-plane loading such as blasts, impact, etc. not only to homogenous materials, but also to heterogeneous materials, such as composites. In the present study, low-velocity impact response of composite laminates was studied using ANSYS Workbench finite element code (FEM). The material of composite structure is orthotropic material reinforced by carbon fibers and have been taken from the Workbench material database. Plate consists of layers which are reinforced with unidirectional fibers in hexagonal and square array. Layer is considered as homogeneous transversely isotropic and layer stacking sequence is symmetrical or unsymmetrical.


Introduction
The impact response of structural components is a long-standing issues in many technological fields. In the last few years, the modeling of fiber/polymer composites using finite element method (FEM) has gained a great interest of researchers [1][2][3]. Nowadays, carbon fiber reinforced composites, CFRs, are widely used in various engineering fields where high specific strength and stiffness are mainly required. These characteristics are crucial for applications in aeronautics and military industries [4].
Impact is typically a dynamic process in which the kinetic energy of impactor at the beginning of the process is transformed into other types of energy. However, low resistance of composite materials in transverse impact necessitates analyzing of their damage mechanics to ensure safe implementation during service life [5]. Composites, unlike metals, are damaged in several ways including fiber breakage, matrix cracking, fiber/matrix deboning and delamination. Damage reduces stiffness and load carrying capability of composites, and it propagates by service load and eventually could lead to catastrophic failure. Gaining insight on damage in composites is a difficult task due to the various damage mechanisms involved and difficulties in conducting experimental tests. For these reasons, simulation of damage and damage modeling is an important aim [6,7]. Significant work is done to model composite impact. Wang et al. [8] conducted an experimental and numerical effort to study damage and failure of a laminated composite under low velocity impact. An elastic-plastic damage model for three dimensional FRC composite under low velocity impact using ABAQUS/Explicit code was developed [9]. The global -local techniques with connection with FEM was used for FEM simulation for impact induced damage in composites and cohesive elements. Cohesive elements and Hashin's criteria were employed in their research [10]. In [11] employed a 3D failure model based on Continuum Damage Mechanics (CDM) for low velocity impact test using LS-DYNA explicit finite element code. Damage modeling usually achieved using CDM [12]. Shi et al. [13] employed the stress based criteria for damage initiation and fracture mechanics technique to model damage evolution. The FEM model was implemented in ABAQUS/Explicit with a used-defined subroutine (VUMAT) and NDT X-ray radiography is used to observe the damage and failure after low velocity impact [14].
At the present time mainly FEM is used for numerical modeling of composite structures. In general it is well-known method [15,16]. There are other numerical methods that can be more efficient than the FEM. The boundary element method (BEM) is also suitable for solving some of problems with boundary and initial conditions. Disadvantage of classical BEM is that must be known fundamental solutions for general continuous nonhomogeneous continuum [17,18]. Other methods are meshless methods that do not require mesh generation. This fact decreases computational time for data transfer from one calculation step to next one. Using meshless method based on a Petrov -Galerkin formulation remove the well-known disadvantages of linear FEM. Nowadays a local Petrov -Galerkin formulation is used also for analyzing of FGM composite shells [19][20].
Low-velocity impact can be considered a quasi-static case, the upper limit can vary from 1 to 10 m/s, depending on the stiffness, material properties and impactor weight and stiffness. In low velocity impact, dynamic response is important because the contact duration is long enough for the entire structure to respond to the impact and as a result more energy is absorbed. For low-velocity shock damage, the mode of failure and energy absorption depends greatly on the sample size, stiffness and initial conditions [21].
In this paper low velocity impact response of composite laminates reinforced by carbon fibers was investigated using ANSYS Workbench finite element code and the various FEM models was used to verify the validity of the results.

Modeling approaches
The behaviour of impacting bodies is a topic in rigid-body dynamics that does not arise in formulating other problems in mechanics [22]. The classical theory of rigid-body impact is primarily based on impulse-momentum balance relations and does not involve mathematical difficulties as is the transient stress wave propagation. For the perfectly elastic two-body collision between smooth surfaces, the impulse-momentum and conservation of mechanical energy laws are sufficient to determine the post-impact velocities [23]. However, when impacts produce some indentation or permanent deformation, the law of conservation of energy is replaced by the restitution relationship. The Hertz theory predicts the stress distribution in the contact zone between two solid bodies of revolution. In addition, absolute values of the contact size, compression, and maximum pressure were obtained. The Hertz theory also allows one to calculate the normal and shear stress distribution inside the solid.
In general there are four types of analysis for low speed collisions, associated with particle impact, rigid body impact, transverse impact on flexible bodies (i.e. transverse wave propagation or vibrations) and axial impact on flexible bodies [24].
Transverse impact on flexible bodies 4.
Axial impact on flexible bodies For bodies that are hard (i.e. with small compliance), only very small deformations are required to generate very large contact pressures. Although the contact area remains small in comparison with cross-sectional dimensions of either body, the contact pressure is large, and it gives a large stress resultant, or contact force. From an analytical point of view, the most important consequence of the small compliance of hard bodies is that very little movement occurs during the very brief period of contact; i.e., despite large contact forces, there is insufficient time for the bodies to displace significantly during impact. This observation forms a fundamental hypothesis of rigid body impact theory, namely, that for hard bodies, analyses of impact can consider the period of contact to be vanishingly small.

Impact models in rigid-body dynamics
Kinetic energy E k (t), which will occur when an object with the mass m (impactor + weight) hits a composite plate at a velocity of v 0 , which could be represented as follows The velocity values in relation to the impactor can be represented through the following formula by using the F force value, which is experimentally obtained [25] v Thus impact energy, which is created in the composite plate, is obtained as the total of elastic energy E e and absorbed energy E a values, so The elastic energy region represents the quantity of the rebound energy from the composite specimen.

Finite element models
With the rapid development of computational mechanics, great progress has been made in numerical studies of the problem. By use of the FEM, many contact problems with complicated body shapes can be solved with useful accuracy [26]. Depending on the purpose of the analysis, different modeling techniques for composites can be used:  Microscopic modelling -the matrix and reinforcement material are both modeled separately as deformable continua.  Macroscopic modelling -the composite is modeled as a single orthotropic material or a single fully anisotropic material.  Mixed modelling -the composite is modeled by a number of discrete, macroscopically modeled reinforced layers.  Discrete reinforcement modelling-reinforcement modeled with discrete elements or other modeling tools (e.g., rebar).
 Submodeling -useful for studying stress concentrations around the tips of reinforcing fibers Fig.1 shows the impact of the rigid body (impactor) on the laminate composite plate. At the t moment, the equilibrium equation can be deduced as: Fig. 1. Impact of the impactor on the composite laminate plate.
where t  and t  are the density and dynamic friction coefficient of laminate at t moment, respectively where, S  represents the stress boundary.
The equivalent integration of the equilibrium equation and the load boundary condition can be expressed, as follows: As there will be geometric nonlinearity during the deformation of composite under the lowvelocity impact load, strain tensor at the t moment can be expressed, as follows Decompose the above equation in to linear and nonlinear terms There is the following relationship of t ij  and t kl  at the time where tt ijkl Q  represents material elasticity matrix at tt  moment, and it can be obtained by coordinate transformation Using equations (8) and (9), the stress equilibrium equation at each moment and obmiting higher-order items can finally be deduced as [27] ∫( where t i T and n i T are the surface force at t moment and n th step in numerical analysis respectively; ij  represents the strain at t moment; and, ij   is the nonlinear term of strain increment. Equation (11) can be solved by FEM and after discretization we get the following equilibrium equations [28] ̈+ ̇+ = (12) with initial conditions given by ( = 0) = 0 , ̇( = 0) =̇0 and M, C and K are the global mass, damping and stiffness matrices, F is the vector of external loads. Explicit dynamics solves the above equation using the following algorithm: Explicit dynamics uses the principle of conservation of energy to monitor the solution accuracy. It calculates overall energy at each cycle. If the energy error reaches a threshold, the solution is regarded as unstable and stops [29]

Finite element model of composite plate
In the next simulations were considered composite plate with dimensions 100 x 100 x 2 mm composed from eight layers. The shell geometry of the composite structure was developed using ANSYS/Workbench. Due to the symmetry, only quarter of the geometry was modelled to save the computational cost ( Fig. 2 and Fig. 3). This module involves the formation of conventional and volumetric shells. This module defines the individual layers of the composite structure, the type of integration rule, symmetry, material properties, thickness, orientation and the number of integration points of the layer.   Thick mesh models and quarter shell model consisted of 934 linear shell elements with 1001 nodes (Fig. 4). The incident body was modeled as an ideal rigid body with 5442 elements and 1245 nodes. The whole model consisted of 6376 elements and 2246 nodes. Mesh of collision surface is refined while the other parts are coarsely meshed to reduce the computational time. The composite is stayed on the base to simulate real situation during impact. The base and drop mass is modeled as rigid bodies. For the 3D dynamic analyses of the composite plate model is used as the simulation platform. Therefore, to perform the analyses of the composite plate, the domains are divided into smaller subdomains (consisting of geometric primitives, such as hexahedra/tetrahedral in 3D and quadrilaterals/triangles in 2D). The governing equations are then discretized and resolved in each of these subdomains. The design and simulations are performed on 16 JAOIGR, 2.4 GHz CPU and 256 GB of RAM.

Material properties
The material of present composite is quasi-isotropic laminate [31]. Material parameters of the laminate plate are listed in Table 1.

Analysis settings
To model the impact, ANSYS Workbench/Autodyne Explicit is used. ANSYS Autodyn explicit solver is used for the solution of discretized FE equations. ANSYS Autodyn simulates the response of materials to short duration severe loadings from impact, high pressure or explosions. It is best suited for simulating large material deformation or failure. Explicit time integration is well suited to dynamic contact-impact problems for the following reason [30] 1. The time steps are small because of stability requirements, so the discontinuities due to contact-impact wreak less havoc. 2. Neither linearization nor a Newton solver is needed, so the deleterious effects of discontinuities on Newton solvers are avoided.
The Lagrangian FE solvers enable fast and efficient solutions when looking at structural components subjected to shock loading and large deformations. A problem encountered when dealing with numerical time integration in the presence in the solution of unphysical high-frequency oscillations resulting from a spatial discretization, which adds to the loss of accuracy in the solution. Indeed, even for stable integration schemes and in the absence of hourglass in the element formulations, most from eigen modes they have no physical meaning, but numerical. To control these spurious frequencies, numerical dissipation can be introduced in time integration scheme. Although numerical dissipation has been widely used in combination with implicit schemes to enhance the convergence of the iterations, some recent works have shown interest in numerical dissipation when dealing with explicit schemes [32].

Numerical results
In the calculations, a fall of a body weighing 364g + 30.7g = 394.7g is taken from the following heights: 100 mm, 150 mm, 200 mm, 250 mm, 300 mm, 350 mm and 400 mm. The calculations were performed using shell and volume models. The shell models were designed as symmetric quarter models and full models. Calculations show that full shell model has better convergence than a quarter model -the generation of high hourglass energy. We note, that hourglass modes are nonphysical, zero-energy modes of deformation that produce zero strain and no stress. Hourglass modes occur only in under-integrated (single integration point (solid, shell, and thick shell elements. The disadvantage of full shell models is 4 times more number of elements and longer calculation times. The volume models show a different character of the solution than the shell models -a later increase in force and also a different (flatter) character of the course of contact force. The abbreviations DC, D, DF, DV, etc. are used to indicate the curves in the graph legend. The first letter of the abbreviation denotes a variant of the arrangement of the fibers in the composite. The second abbreviation C, F and V denotes coarse shell model, fine shell model and quarter volume model. Fig. 5 shows the values of contact forces for variant D of the arrangement of fibers when the contact forces values reach the highest values. In the following figures, we compare the contact forces only for impact height 400 mm (impact velocity is 2800.9 mm/s). In Fig. 6 we  can see comparison of course of contact force for full and a quarter shell model for height H = 400 mm. Fig. 7 compares of contact force for whole shell model and half volume model for H = 400 mm. The contact force for the arrangement of the fibers C has not been evaluated because the hourglass energy is too high and the convergence of the solution is very bad. Further test are need to obtain acceptable results.
The interesting behaviour of the contact force in variant B of the fiber arrangement was verified on the modified variant B (working name B-IMPERF). Variant B-IMPERF has a fiber arrangement slightly modified from the original variant B as follows: A = [1/−1/0/ 0] . We wanted to find out if the contact force would drop rapidly after reaching its initial maximum. Fig. 8 shows that the courses of variant B and variant B-IMPERF are very similar (they are identical in the critical area) and variant B-IMPERF began to diverge later than variant B. This indicates a physically realistic behaviour. Fig. 9 presents course of total displacements for variant B in time 9.00e-05 s, when the contact force dropped. The maximum value of the displacement is 0.25117 mm.

Conclusion
We have analysed by the FEM transient elastic deformations of a fiber reinforced laminate impacted at normal incidence by a rigid sphere moving at a slow speed. The shell models were designed as symmetric quarter models and complete models. It turns out that a