Reliable computation of contact force in FRP composite laminates under transverse impact

A simple and computationally efficient adaptive finite element analysis strategy has been adopted for accurate and reliable evaluation of contact force under transverse impact. Contact of a spherical and cylindrical impactor on FRP composite laminates are considered. Adaptive mesh refinement enables the finite element mesh to be obtained iteratively and automatically for the solution to have the desired level of accuracy. The refined meshes influence the calculated contact force and contact period appreciably. Influences of impactor velocity, mass of impactor, mass of the plate on discretization error as well as on contact force history have also been studied and are found to be sensitive.


Introduction
FRP composite laminates, in-spite of having many inherent advantages are susceptible to transverse impact. Especially in case of low-velocity impact, damages such as matrix cracking, fiber breakage and delamination occur at certain plies. Such defects generally remain embedded in nature and are apparently not visible. Once such a defect occurs at a particular location in a laminate, it becomes the weakest point of the laminate and failure initiates from that location leading to the final fracture of the laminate. Therefore, accurate predictions of such defects are of significant importance from the design point of view of structural components made of FRP composites. Since the state of stress at a point dictates the occurrence of such defects, accurate determination of contact force is of prime importance. Many works have already been reported in this direction studying the effects of important parameters such as impactor mass, impactor velocity, plate mass, plate thickness on the magnitude of contact force and the contact force history. Caprino et al. [3] developed a simple model based on energy consideration for studying the response of composites under low-velocity impact and compared the results with experiments. Sun and Chen [2] reported the effects of impact velocity, initial stress, mass and size of the impactor on contact force history, deflection and strain in an initially stressed graphite epoxy plate subjected to impact. Hsi-Yung and Springer [4] developed a method for prediction of location and size of delamination in a laminated plate subjected to low velocity impact using 3D finite element and dimensional analysis. His-Yung and Chang [5] performed a 3D finite element analysis for determination of contact force history and went on to predict delamination initiation in the laminate due to transverse impact. Cantwell and Morton [22] conducted tests to compare the damage initiation from low and high velocity impact on CFRP laminates. Choi and Chang [6] performed experiments to study impact damages in composite laminates due to low velocity impact. Choi and Hong [7] proposed a simple method for predicting impact force history using finite element analysis with modified Hertzian contact law. They also determined the frequency characteristics of the impact force history from modal analysis and compared with natural frequency of the system. Gao and Kim [15] performed a 3D penalty finite element analysis for dynamic contact problem, which can account for span or thickness and stacking sequence of the laminate. They further evaluated the impact behaviour of curved composite laminates. Lam and Sathiyamoorthy [13] presented a theoretical model for analysis of low velocity impact of multiple masses with different initial velocities on laminated beam specimen. However, without assessment of the reliability of the results it is hardly reasonable to use numerical methods like the finite element method in contact impact problems of plate etc. Designing of an appropriate finite element mesh for contact impact problems is not an easy task and it could be very erroneous if made on the basis of intuition alone. Hence, a mesh obtained automatically with densities varying according to the requirement of the problem is always a welcome choice. A simple, reliable and computationally efficient adaptive finite element strategy becomes very useful to satisfy such requirement.
At the present state of knowledge, the work recorded on adaptive finite element analysis under dynamic loads is limited. Wiberg and Li [16] proposed a postprocessed type of a posteriori estimates in space and also in time while using direct integration for dynamic response evaluation. It updated the spatial mesh and time step so that the discretization errors were controlled. This process is not only time consuming but also gives rise to new errors as nodal values need to be interpolated from the previous mesh to the newly generated mesh whenever a mesh is changed. Combe et al. [9] proposed an error estimate based on the concept of error in the constitutive relation for transient dynamic loading. Wilson and Joo [10] arrived at the final mesh using Ritz vector as basis of transformation. In their investigation, the authors made use of modal participation and amplification factor and obtained error estimates based on Babuska's criterion using amplified Ritz modes. Cook and Avrashi [20] discussed the procedure for estimating the discretization error of the finite element modeling as applied to the calculation of natural frequency of vibrations. Meshes were obtained corresponding to each mode. Dutta and Ramakrishnan [1] proposed a measure for discretization, which is a logical extension of Zienkiewicz and Zhu [18] error criterion by involving time integration to consider the variation of response with time. Mahajan and Dutta [19] utilized the adaptive strategy suggested by Dutta and Ramakrishnan [1] for determination of accurate impact force on isotropic plate under low velocity impact. The adaptive strategy as proposed by Dutta and Ramakrishnan [1] is simple to implement, compu-tationally efficient and quite reliable. In the direction of adaptive finite element study on FRP composites, Filho and Dow [8] developed a procedure to estimate discretization error in laminated composite plate finite element model and a refinement strategy to achieve an acceptble level of accuracy. Thus, literature review reveals that many works have already been reported in the areas of adaptive finite element study of structures under general loading. Further, a few literatures are also available on adaptive study of contact force determination under transverse impact on isotropic plates. Since the choice of appropriate finite element mesh greatly influences the contact force, stress distribution in the vicinity of point of contact, the use of adaptive finite element method will greatly help in eliminating the uncertainty of accuracy of those quantities. This is especially important where one would be interested in predicting the damages due to impact. Unless the stresses responsible for damages are evaluated accurately, damage prediction may become highly erroneous. Thus, the major objective of the paper is the reliable evaluation of contact-impact force and stress distribution. In this paper, the adaptive strategy given by Dutta and Ramakrishnan [1] for space discretization has been adopted for the accurate evaluation of contact force on FRP composite laminates under low velocity impact. However, time discretization generally does not pose any problem for contact impact analysis as the time step size adopted for simulation of contact history is fine enough than what is demanded for the control of error due to time discretization.

Problem definition
In this investigation, FRP composite plates subjected to both point-loading and line-loading impacts have been studied. FRP laminated composite plates with different boundary conditions have been considered. Spherical and cylindrical impactor strike the laminated composite plate as shown in Fig. 1(a and b). Appropriate contact laws have been assumed for the forcedeformation relationship. The computations of contact force have been carried out for different impactor mass, impactor velocities and plate mass, while ensuring that the numerical errors are kept within the desired limit.

Finite element analysis
A Cartesian co-ordinate system is chosen in such a way that the x-y plane coincides with the midplane (b) Line contact (a) Point contact of the plate. The components of displacements are assumed to be as follows: where t is the time, u, v, w are the displacements in the x,y,z directions respectively, w o denotes the transverse displacements of the midplane, ϕ x and ϕ y are the middle-surface slopes in the xz and yz planes respectively. The parameters w * o , ϕ * x , ϕ * y are the higher order terms in the Taylor's series expansion [14]. The stresses and strains in the laminates are calculated by transient finite element method by considering a nine noded element for the assumed displacement model.
The element stiffness matrix is given by where, |J| is the determinant of standard Jacobian matrix, B is the strain displacement matrix and D is the constitutive matrix. The element mass matrix is given by where, shape function matrix, in which I 1 and I 2 are the usual inertias corresponding to w o and θ x (also θ y ) degrees of freedom, while I 3 and I 4 are the higher order inertia terms corresponding to the higher order degrees of freedom w * o and θ * x (also θ * y ) respectively and these are defined explicitly as follows: where ρ L is the material density of the Lth layer.

Analysis of contact impact
When a plate is impacted by a mass, the magnitude of contact force, which results because of the impact, is not known a priori. This contact force needs to be calculated before the plate motion is analysed. The evaluation of contact force depends on a contact law, which relates the contact force with the indentation.

Contact of a sphere with a plate
The projectile is modeled as an elastic body with a spherical nose as shown in Fig. 1(a). During loading and unloading, the contact force distribution is determined using the Hertizian contact law. Thus, the contact force F can be related to the indentation depth α, the distance between the center of the projectile's nose and mid surface of the plate [4] as where K is the modified constant of the Hertz contact theory and is given by where r, V i and E i are the local radius, the poisons ratio and Young's modulus of the material of the impactor respectively. E yy is the transverse modulus normal to the fiber direction in the upper most composite layer. Upon unloading, the contact force is simulated by the following relation developed by Yang and Sun [21].
where f m is the maximum contact force just before unloading, α m is the maximum indentation corresponding to f m and α 0 is the permanent indentation during this loading/ unloading process.
The permanent indentation can be determined from the following expression [21] α 0 = 0 when α m < α cr where α ar is the critical indentation and is approximately 0.1016 mm for glass/epoxy and 0.08026 mm for graphite/epoxy [21]. By neglecting damping, the governing equation can be written as At time t + ∆t, Eq. (8) is written as where M and K are the mass and stiffness matrices and q, d,d are the force, displacement and acceleration vectors respectively. Newmark-β method is employed to obtain the solution to this equation. Accordingly, the velocity and acceleration vectors at t + ∆t are written asḋ The parameters β, γ are constants whose values depend on the finite difference scheme used in the calculations. Here, we have used the constant average acceleration method, which is implicit and unconditionally stable. For this method β is 0.25 and γ is 0.5.   Substituting Eq. (11) into Eq. (9), we obtain whereK is the effective stiffness matrix andF(t + ∆t) is the effective force vector. The parameters are defined asK where h(t) is given as µ Fig. 9. The impactor velocity and the plate center velocity corresponding to Mesh 3.
The unknown in Eq. (12) are the displacement vector d (t+∆t) and the force vector q(t + ∆t), since the displacement, velocity and acceleration at time "t" are known at every point inside the plate.
In the absence of body forces, for a plate initially at rest, the only concentrated load is contact force caused by the impactor. We define a scalar force "f", which is a point force acting perpendicular to the plate at the contact point and has a magnitude equal to contact force. The contact force q is then written as where u is a unit vector, which is equal to −1 in the direction of the contact force at the point of contact and zero elsewhere. The displacement vector d is expressed as the sum of the displacements due to the force h and the contact force q.
Thus we can writê andKd q t+∆ = q(t + ∆t) The solution of Eqs (19) and (20) now proceeds in two steps. First, the forces h and the displacement d h are calculated from Eqs (15) and (19) respectively. Second, the contact force vector q and the displacement d q are calculated as follows. At time t + ∆t, Eq. (16) is written as Equations (20) and (21) yield For a unit contact force (f (t + ∆t) = 1), Eq. (22) becomeŝ where d u t+∆t is the displacement caused by the unit contact force. From Eqs (22) and (23) Equations (17) and (24) give In Eq. (25), the unknowns are displacement vector d and the scalar contact force "f" at time t + ∆t. In order to find these two unknowns, another expression for contact force is required. The expression for contact force during the loading and unloading process has been given in Eqs (4) and (6)   where "f" is the time varying constant force and "m" is the impactor mass. By combining Eqs (26-28) and the contact laws as given in Eqs (4) and (6), we obtain the following expression for the contact force: During loading During unloading  The contact force "f" is calculated either by Eq. (26) during loading or by Eq. (27) during unloading by employing Newton-Raphson method.

Contact of a cylinder with a plate
A Hertzian distribution of pressure is assured to act on the plate, which is uniform along the length of the barreled cylinder. In order to reduce the concentration of the stress at the ends, the axial profile of the cylinder should be slightly barreled. The displacement at the centre of the contact [12] is where 'f is the total contact force, '2l' is the length of the cylinder in contact, and a = 4f r ΠE * and 'r' is the radius of the cylinder, E i and E are Young's modulus of impactor and target respectively, v i and v are Poisson's ratio of impactor and target respectively. Equation (31) is considered as the relationship between indentation depth and contact force for both loading and unloading. Equation (31) has been originally derived from the contact between two elastic, isotropic bodies. However, in the absence of any relationship specifically for FRP composite laminates, Eq. (31) has been utilized in the present analysis. The same solution procedure which is used for spherical contact is followed here as well and leads to the expression for contact force during loading and unloading as f (t + ∆t) = π1E * 1.886 + ln(1/a)

Adaptive analysis
In order to develop reliability on the obtained contact history, an estimate of error involved in the finite element solution is evaluated. A solution of desired accuracy can be achieved iteratively based on the information obtained from the error analysis. The error is calculated based on the discontinuous stress field. Finite element stresses are calculated where D, B and z are usual elasticity, strain displacement matrices and displacement vector respectively. The approximate solution σ containing discretization errors differs from the accurate solution σ * and the difference is the error e σ Thus, e σ = σ * − σ A recent and elegant technique for the evaluation of accurate stress is done using superconvergent patch recovery method by Zienkiewicz and Zhu [17]. In the recovery process, it is assumed that the accurate nodal values belong to a polynomial expansion σ * e of the same complete order "p" as that present in the basis function "N" and which is valid over an element patch surrounding the particular assembly node considered. Such a 'patch' represents a union of elements containing this vertex node. This polynomial expansion will be used for each component of σ * e and we get σ * e = P a (36) where "P" contains the appropriate polynomial terms and "a" is a set of unknown parameters. For a general plate and shell problem, the polynomial can be expressed as a function of x, y and z. However, since only plate problems are solved in this paper, x and y are considered only as the parameters in the polynomial expression. Thus it can be written as and a = [a 1 , a 2 , a 3 , a 4 , . . . . . . . . .] T (38) The determination of the unknown parameters "a" of the expression given in Eq. (36) is best made by ensuring a least square fit of this to the set of superconvergent points existing in the patch considered. To do this we minimize where (x i , y i ) are the co-ordinates of a group of sampling points, n=mk is the total number of sampling points and k is the number of the sampling points on each element m j (m j = 1, 2, . . . . . . m) of the element patch. The minimization condition of F(a) implies that "a" satisfies This can be solved in matrix form as Once the parameter "a" are determined the recovered nodal vales of σ * are simply calculated by insertion of appropriate co-ordinate into the expression for σ * e . The stresses in the nodes inside the patch are evaluated using "a" from Eq. (41). The stresses at boundary nodes can be determined using interior patches.
The mathematical basis and broad details of the algorithm for discretization error measure are described in reference [1]. However, the complete procedure of discretization error measure is given below in an algorithmic form.

Numerical results and discussion
In order to validate the code, an isotropic steel square plate 0.2 m long, 0.2 m wide and 0.008 m thick rigidly fixed on its four edges is considered [5]. The plate is subjected to an impact induced by a steel ball of 0.02 m diameter to its center with a velocity of 1.0 m/s. A quadrant of the square plate is considered for the analysis and in order to have reliable results, the adaptive analysis is carried out for a prescribed domain discretization error limit η = 5%. Three analyses are needed to bring the overall domain discretization error η near to 5%. Meshes have been generated using codes developed based on advancing front technique. Corresponding to each analysis, a new mesh is generated with changes taking place globally. However, the mesh generation is not automatically adapted. It is achieved interactively after each run based on the error information of previous run. The initial mesh (Fig. 2a) has 9 elements and the overall domain discretization error η = 46.12%. The intermediates mesh (Fig. 2b) is having 266 elements and η =11.87 %. The final adaptively refined mesh (Fig. 2c) has 569 elements with η = 5.51%. Figure 3 shows the plot of contact force vs. time corresponding to all the meshes considered during the adaptive process and it is clear from the Fig. 3 that the finite element meshes influence the value of the contact force so calculated. Maximum value of the contact force has reduced from 1407.041 N corresponding to initial mesh (Fig. 2a) to 1283.905 N corresponding to final mesh (Fig. 2c). Figures 4 and 5 show the plot of displacement and velocity of the center of the plate and impactor corresponding to the refined mesh (Fig. 2c). Contact period is also observed to increase by a few microseconds, as the mesh is refined. The solution compares well with series solution of Karas [11] and the finite element solution as presented in ref. [5]. The study on isotropic plate demonstrates that control on discretization errors is required for reliable results.
A Composite plate 0.0762 m long, 0.0762 m wide made of T300/934 graphite/epoxy prepregs with stacking sequence of [0/−45/45/90] 2S clamped along all four edges is considered for study. The FRP laminate is impacted at the center by an aluminum sphere of diameter 0.0127 m with a velocity of 25.4 m/s. The ply thickness is 0.00015875 m. The properties of graphiteepoxy are given in Table 1. The analysis has been carried out by considering a quadrant of the symmetric plate and a finite element mesh having 9 elements similar to as shown in Fig. 2a. Overall domain discretization error η = 34.98%. Three iterations are required to attain the desired accuracy of near to 5%. The intermediates mesh (Fig. 6a) is having 208 elements and η =13.67%. The final mesh (Fig. 6b) has 565 elements and the overall domain discretization error η = 5.43%. Figure 7 shows the plot of contact force vs. time corresponding to all the meshes considered during the adaptive process. The impactor loses contact with the target and regains contact after a period of time depending upon the flexibility of the plate. The relative displacements between the impactor and the target and the contact force gradient have been utilized as indices to trace the loading, unloading and reloading. The separated parts represent initial loading /unloading and reloading / unloading stages. It is again evident from Fig. 7 that the finite element meshes influence the value of the contact force so calculated. Maximum value of the contact force has reduced from 2297.44 N corresponding to initial mesh (with nine elements) to 2000.4 N corresponding to final mesh (Fig. 6b). Figures 8 and 9 show the plot of displacement and velocity of the center of the plate and impactor corresponding to the refined mesh (Fig. 6b). The variation of six stress components with time near the point of contact is plotted as shown in Fig. 10 corresponding to two meshes having nine elements and 569 elements (Fig. 6b). It is quite evident from the stress plot that there is a substantial change in the stress values as the refinement in finite element meshes are carried out and a reliable stress distribution can be ensured by providing a control on discretization error through adaptive analysis. The solution compares well with the results presented by Wu and Chang [5].
Further, the same composite plate, which is used in the case of a point contact, is considered again with fixed condition on all four edges. The plate in this case is subjected to an impact induced by a steel cylinder of 0.0127 m diameter and 0.0254 m length to its center with a velocity of 25.4 m/s. The analysis is carried out by considering a quadrant again and the initial finite element mesh is also having 9 elements similar to as shown in Fig. 2a with overall domain discretization error η = 34.05 %. Adaptively refined mesh (Fig. 11b) corresponding to a target accuracy of 5 % is obtained after three iterations. The intermediates mesh (Fig. 11a) is having 168 elements and the overall domain discretization error η = 8.92%. The final mesh (Fig. 11b) has 496 elements and η = 5.12%. Figure 12 shows the plot of contact force vs. time corresponding to all the meshes considered during the adaptive process. Maximum value of the contact force has again reduced from 4703.35 N corresponding to initial mesh (with nine elements) to 4603.17 N corresponding to final mesh (Fig. 11b).
Parameters like impactor velocity, impactor mass and plate mass have been varied to study their influence on discretization error. The impact of an aluminum ball on FRP composite laminate is considered. The geometrical and material data are the same as the point contact problem for FRP composite laminate described earlier. A finite element mesh similar to as shown in Fig. 2a is considered for the study. When the velocity of the impactor is varied keeping the mass of the impactor unchanged, it is observed that the domain discretization error increases ( Table 2). Hence mesh refinement requirement also increases when the velocity of the impactor is increased. Further, by considering a mesh refined for a target level of 5% accuracy (Fig. 6b), which ensures a reliable response, velocity of the impactor is varied again. Figure 13 shows the variation of impact force with respect to time for different velocities of impactor for constant impactor mass and it is observed that the impact force increases, the contact period decreases as the velocity of the impactor increases. The impact force increased from 1154 N (v = 15 m/s) to 3204.76 N (v = 40 m/s) and the contact period decreased from 87 µsec (v = 15 m/s) to 80 µsec (v = 40 m/s). The displacements and velocities of impactor and plate also vary appreciably, when subjected to different impactor velocities. Table 3 shows the values of overall domain discretization error η for cases with different Mi/Mp ratios with fixed velocity of impactor, where Mi is the mass of impactor and Mp is the mass of the target plate. The discretization error decreases with the increase in the ratio of Mi/Mp and hence coarser and coarser meshes would be required. The analysis has been carried out considering a nine-element mesh similar to as shown in Fig. 2a. A refined mesh (Fig. 6b.) is considered again for the analysis for different Mi/Mp ratio's with fixed velocity of impactor. It is observed (Fig. 14) that the impact force and contact period increase with the increase in Mi/Mp ratio. The impact force increased from 926.01 N (for Mi/Mp = 0.016) to 2297.78 N (for Mi/Mp =0.221) and contact period increased from 29 µsec (Mi/Mp = 0.016) to 103 µsec (Mi/Mp = 0.221).

Conclusion
A discretization error measure based on Z-Z type of error estimate is effectively used. For a prescribed limit of η, the optimal mesh is obtained iteratively. The computed contact forces change with the change in finite element mesh, where mesh is refined nearer the point of application of the contact force as the stress gradient will be high there. Discretization error as well as contact forces are affected by the parameter like impactor velocity and ratio of mass of impactor to mass of plate.