Elastic plate under low velocity impact: Classical continuum mechanics vs peridynamics analysis

1.   Introduction

Thin-walled structural glass panels are usually subjected to various quasi-static (snow weight, wind pressure), dynamic (hailstorm, wind gust) and non-stationary thermal loadings. In order to investigate deformation and stress states in monolithic and laminated glass components under dynamic loadings several approaches are developed. Glass plies and glass laminates with shear-stiff interlayers can be analyzed with first order shear deformation theories of beams and plates ^[1,2]. Layer-wise and zig-zag theories are proposed for glass laminates with soft interlayers such as beams ^[3,4], plates ^[5] and shells ^[6]. User-defined finite elements for laminates with soft interlayers are developed and utilized in ^[7,8,9]. To account for transverse normal deformation, stress concentrations and contact behavior as a result of ball drop, three dimensional continuum mechanics theory and solid finite elements are applied in ^[4,10,11,12] among others. Critical zones of failure initiation in glass plies can be predicted applying the continuum damage mechanics. Various scalar and tensor-valued damage state variables for brittle materials are proposed in ^[13,14] among others. To evaluate onset of failure propagation of pre-cracked structures, fracture mechanics criteria, for example the the energy criterion ^[15] and cohesive zone models ^[16,17] are developed. For monolithic and laminated glass plates under impact loading, the analysis of fracture process including the formation of damage patterns, as well as post-fracture regimes after the failure of a glass ply are of practical importance. In the classical continuum mechanics (CCM) the damage evolution equation is formulated in terms of the local characteristics of the stress and strain state. With this approach damage initiation zones in a component can be well predicted, e.g., ^[18,19]. However, the consideration of fracture process leads in general to a loss of well-posedness of the initial-boundary value problem of CCM ^[20], and as a result to a spurious mesh dependence of finite element solutions. To regularize the finite element formulations, phase field and gradient-enhanced damage theories were developed ^[21,22]. Examples of crack propagation and crack branching simulations with the phase field models are presented in ^[23,24], among others. The phase field theory is applied in ^[25] to analyze fracture in monolithic and laminated glass plates. The results illustrate that zones of damage initiation are well predicted. However, the standard phase field formulations fail to capture complex damage patterns, usually observed in float glass specimens. In addition, non-local long-range interaction forces are of key importance in the analysis of fracture phenomena ^[26,27].

Peridynamics (PD) is a non-local theory which operates with long-range force interactions ^[28]. The deformation gradient, its higher gradients and gradients of internal state variables are not required. In addition to contact surface forces, long range bond force interactions are considered. Many studies of complex fracture processes such as crack initiation ^[29,30], crack branching ^[31], crack kinking ^[32] and crack interaction with initial heterogeneities ^[33,34] are recently performed with PD. For the analysis of glass panels PD plate theories can be applied ^[35,36].

Within the framework of PD, new parameters, such as the horizon are introduced. Furthermore, the treatment of constraints, applied external forces and contact conditions requires special considerations due to the non-locality. Several available analytical solutions to PD equations ^[37] indicate that fictitious domains are required to apply constraints and loads. Therefore an important step is to calibrate the PD model and to assess the adequacy of results.

The aim of this paper is to compare the classical continuum mechanics and the PD models as they predict time-dependent deformation states in the monolithic glass plate after the ball drop. In our study damage processes in glass will be ignored and the material behavior will be assumed to be linear elastic. The paper is organized as follows. In Section 2 governing equations of classical continuum mechanics and PD are presented in order to highlight the basic features and differences between the approaches. In Secttion 3 mechanical models for the ball drop simulation are discussed in detail. An emphasis is placed on the discretization including the finite element mesh, the PD node lattice and the time stepping, as well as constraints and contact conditions in both finite element and PD models. Results of analysis are presented in Section 4. Time-dependent deflections of the plate after the ball drop will be analyzed in detail, and the results based on the finite element and PD analysis will be compared.

2.   Governing equations

Let $\mathbf {X}$ be the position vector for a point $P$ in a reference configuration, and $X^i, i = 1, 2, 3$ are the corresponding coordinates, as shown in Figure 1. The position vector of this point in the actual configuration (designated by $P'$) is specified by $\mathbf x$.

The motion motion of the continuum is defined by the following mapping

The basic problem of continuum mechanics is to compute the function ${ {\boldsymbol{\varPhi}}}$ for all vectors $\mathbf {X}$ within the body in the reference configuration, for the given time interval $t_0\leq t \leq t_n$ as well as for given external loads and temperature. The displacement vector $\mathbf {u}$ is defined as follows Eq 2 (Figure 1)

In this section the governing equations of the CCM and PD are presented and compared.

2.1.   Classical continuum mechanics

Classical continuum mechanics operates with local measures of deformation related to the differential line, area and volume elements. The governing equations of continuum mechanics and finite element solution procedures are discussed in ^[38,39], among others. The directed line element in a differential neighborhood of $P'$ in the actual configuration is computed as follows Eq 3

where $\mathbf F$ is the deformation gradient $\mathbf G_i$, $i = 1, 2, 3$ is the local basis, $\mathbf G^i$ is the corresponding dual basis and $\delta_i^k$ is the Kronecker symbol. For many applications in the structural mechanics local strains and rotations can be assumed small. In this case the deformation gradient takes the following form Eq 4

where $\pmb\varepsilon$ is the tensor of infinitesimal strains, $\pmb \varphi$ is the vector of infinitesimal rotations and $\mathbf I$ is the second rank unit tensor. From Eqs 4 the following compatibility condition for the strain tensor can be derived form Eq 5:

Within the CCM only contact surface forces of interaction between any part of the solid and the remainder are introduced. Applying the balance of linear momentum for a part of the solid and the localization, the following equation can be established Eq 6:

where $\rho$ is the density, $\pmb\sigma$ is the Cauchy stress tensor and $\mathbf b$ is the body force vector. The balance of angular momentum results in Eq 7,

The constitutive equation for the stress tensor can be assumed in the form of the generalized Hooke's law. For isotropic materials it takes the following form Eq 8,

where $\lambda$ and $G$ are the Lamé's parameters, Eq 9

$E$ is the Young's modulus and $\nu$ is the Poisson's ratio.

The differential Eqs 4, 6 must be supplemented by the boundary conditions Eq 10:

where $\bar{\pmb u}$ is the given displacement vector, $\bar{\pmb p}$ is the vector of given surface forces and $\pmb\nu$ is the outward unit normal to $A_p$. The vectors $\bar{\pmb p}$ and $\bar{\pmb u}$ can, in general, be functions of coordinates and time. Furthermore the initial conditions for the displacement and velocity vectors must be specified.

2.2.   Peridynamics

PD is a nonlocal continuum theory, in which the material points are assumed to interact through long-range forces ^[28,40,41]. For a recent review on applications of peridynamics we refer to ^[42]. Unlike the CCM, differential line elements are not introduced. Instead, the bond vector $\pmb \xi$ with the finite length connecting two material points $\mathbf X$ and $\mathbf X'$ in the reference configuration is considered as follows Eq 11:

The PD deformation state $\underline {\mathbf Y}\langle \pmb \xi \rangle$ for the bond vector in the actual configuration is defined by Eq 12:

The balance of linear momentum is formulated in the following integral form ^[28,41]

where ${\underline {\mathbf T}}[\mathbf X, t]\langle\mathbf X'-\mathbf X\rangle$ is the PD force density state for the bond $\pmb \xi$. ${\mathcal H_{\mathbf X}}$ is a finite size spherical neighborhood centered at $\mathbf X$, with the radius $\delta$ which is called the horizon. The PD equation of motion Eq13 can also be formulated as follows Eq 14 ^[28]

The PD force density state is related to the deformation state by a constitutive equation. In our study we apply the linear PD solid material model developed in Eqs 15, 16^[43]

with

where ${\underline \omega}\langle \pmb \xi\rangle$ is the influence function. In Eq 16 the bulk modulus $K$ is defined by Eq 17,

The nonlocal dilatation $\theta$, the weighted volume $m$ and the deviatoric part of extension ${\underline e}^d \langle \pmb \xi\rangle$ are defined as follows

For a numerical solution, the domain is descritized into a finite number of nodes with finite volumes corresponding to each node. The integral in the right-hand side of Eq 13 is replaced by a finite sum over the neighbors of the point $\mathbf x$ inside the horizon, while a time integration scheme is applied to solve the initial value problem ^[44]. In our study we use the open-source code Peridigm. For a more detailed overview on the numerical methods we refer to ^[44,45,46].

3.   Models for the ball drop simulation

In this study a drop of a solid steel ball with the mass of 1 kg with initial velocity of 8.85 m/s on a glass plate is simulated. The square glass plate with the length $l = 500$ mm and the thickness $h = 9.9$ mm is considered, as shown in Figure 2.

The corresponding material properties of steel and glass are presented in Table 1.

3.1.   Classical continuum mechanics model

Based on the CCM the impact of the ball on the glass plate is analyzed using the solid hexagonal finite elements (FE) with 20 nodes and 3 degrees of freedom per node. The interaction of the ball and the plate is considered as a one-way contact by the surface to surface algorithm. The resistance caused by air during impact was neglected. The initial time step size for an explicit time integration is chosen as 5$\times 10{}^{-6}$ s.

For the studied plate several computations with different FE meshes were performed. For each FE size maximum deflections were calculated to determine the computational error. The mesh with the element size of 4 mm was chosen such that the total number of elements in the model is 67256. A representative FE mesh is illustrated in Figure 3.

To model the support frame, the area with the width $l_\mathrm{s} = 10$ mm on the bottom surface, see Figure 4 was created. The set of nodes on the frame was selected and deflections were set to zeros.

3.2.   Peridynamic model

To generate the node lattice for the PD model, the plate and the ball are discretized first with a FE mesh generator. From the parent FE mesh the node lattice is created such that the nodes correspond to the centroids of the finite elements. Figure 5 illustrates a parent FE mesh and the corresponding PD discretization.

The total number of nodes is 557496 for the plate, 6084 for the ball and 6364 for the frame, respectively.

Following the recommendation in ^[41] the horizon size is specified as $\delta = 3.015\Delta x$, where $\Delta x$ is the grid spacing. The parent finite element mesh is non-regular leading to the non-regularity of nodes distribution in the PD lattice, see Figure 5. The grid spacing $\Delta x$ for the node set of the plate is specified as an averaged value based on the radii of the inscribed circle and the circumscribed circle of an equilateral triangle, where the edge length is equal to a given plane element size, and a given element size in plate thickness direction. The horizon for the plate and for the frame is 6.96 mm while for the ball it takes the value of 4.47 mm.

The contact between the ball and the plate is modeled applying the short-range force approach, as proposed in ^[44,46] and implemented in Peridigm software. Spring-like repulsive forces of interaction are assumed for nodes of two bodies that are close to each another.

For the solution the explicit time-step method based on the central difference scheme to approximate the second time derivative of the displacement vector, as implemented in Peridigm is used. Details of the method and the critical time step estimation are presented in ^[44]. The time step size in PD analysis was 10$^{-7}$ s.

4.   Results

4.1.   CCM solution

Figure 6 illustrates the displacement of the ball bottom node and the maximum deflection of the plate bottom central node as a function of time.

The result for the contact pressure between the ball and the plate vs time is illustrated in Figure 7. During the introduced time interval several peaks of the contact pressure are observable, Figure 7. In addition, regimes with zero contact pressure can be recognized.

These regimes are explained by the vibration of the plate followed the first impact. Figure 8 illustrates the maximum principal stress vs time at the central point of the bottom surface of the plate. Experimental data for the strength of float glass are presented in^[29].It is worth to note that differences between the surfaces of the float glass should be taken into account, since one surface during manufacturing is in contact with tin (called tin side) and the other with air (called air side). For the tin side the fracture stress value estimated from ring bending tests is in the range 185–198.7 MPa while for the air side significantly higher stress values and much bigger scatter within the range of 438.0$-$713.6 MPa are obtained in ^[29]. If the bottom side of the considered plate is the tin side, then crack initiation in the midpoint of the bottom side is expected after the first stress peak, as shown in Figure 8. If the bottom side of the plate is the air side than the first impact peak would not lead to the failure and cracks would be observable after the second peak if the lowest value of the fracture stress in the scatter band is considered.

4.2.   PD vs CCM solutions

To verify that the loading on the plate was the same by using two approaches, Figure 9 illustrates the displacement of the ball as a function of time. The results based on the finite element and the PD analysis agree well. Figures 10, 11 illustrate the numerical results for the deflection in two points of the bottom surface of the plate vs time. Within the time interval $0\leq t\le 0.002$ s a very good agreement between the results of classical and PD solutions are observed. For the remaining time interval $0.002\leq t\leq 0.005$ s a slight shift of the PD solution towards the time axis is visible.

This can be explained by the different time step procedures implemented inside ANSYS and Peridigm codes. Despite this slight discrepancy, details of plate deflections in the course of vibrations after the ball drop agree well by two approaches. In addition, Figures 12–14 illustrate the deflection fields at three time points. Good agreement between the deflection distributions obtained by FEM and PD can be observed. This indicates that model settings for the non-local analysis including the horizon size, the short-range force contact and the support conditions are well suited.

5.   Conclusions

The aim of this paper was to compare the CCM and the non-local PD approaches in modeling and simulation of the ball drop on the glass plate. The obtained results for the deflection in selected points of the plate vs time as well as deflection fields of the plate at several time steps show very good agreement. This indicates that the non-local theory is well applicable to reproduce the displacement fields. Let us note, that unlike the classical continuum mechanics, both damage and fracture patterns in glass can be reproduced by PD. To this end the constitutive equation with a damage parameter should be calibrated based on experimental data on failure initiation and crack propagation in float glass. PD models of failure initiation and early stages of crack growth under quasi-static loading are discussed in ^[29]. Future studies should be related to the analysis of damage patterns in float glass under low velocity impact.

Conflict of interest

All authors declare no conflicts of interest in this paper.