## Abstract

In glass research, the effect and influence of pre-deformation by thermal or chemical treatment is of great importance when configuring different mechanical properties or scratch resistance on the surface of glasses. In particular, such pre-deformation affects dynamic fracture or damage evolution when glass structures are under impact or collision conditions. Peridynamics provides a seamless approach for the simulation of dynamic damage evolution of the system under aggressive environments. Revising the pair interaction of each material point, the effect of pre-deformation is implemented, and the corresponding damage evolution can be simulated conveniently. Our approach is composed of two steps: first, a static solution is found via energy minimization with thermal boundary conditions in the peridynamic platform. Second, comparing the initial and the pre-deformed structures from the energy minimization, the effect of residual deformation, strengthening and reactive behaviour of brittle structures are seamlessly simulated. The developed methods are applied to the Prince Rupert’s drop and Bologna vial, which are classic examples of strengthened glasses. This study reports the first complete and successful simulation of dynamic behaviour of strengthened glasses, and a significant contribution in simulating residual stress behaviour in any material.

## 1. Introduction

Implanting pre-stress or residual deformation in glasses by chemical or thermal treatment is a common method to strengthen glasses [1]. Compared to plain unstrengthened glasses, strengthened glasses have a much higher surface resistance against scratches and may show drastic differences when they are damaged or fractured, yielding dicing (many cubic pieces) or pulverization (fine fragmentation). Therefore, it is of great importance to understand and model the effects of such residual deformation for the design and manufacture of strengthened glasses or structures. To simulate such pre- or residual deformation in dynamic cases like impact conditions, there have been various ways such as switching between implicit and explicit solvers [2] or coupling them [3] in modern computational solid mechanics simulation tools.

Peridynamics is a state-of-the-art simulation method for solid mechanics and structures, inspired by the method of molecular dynamics or the interaction of pairs of material points [4–7]. Constitutive relations are managed by the interactions among material points, and dynamic damage evolution can be simulated seamlessly. Briefly, peridynamics replaces the divergence of stress in Cauchy’s equation of motion with a volume integral of a pairwise force function acting among a neighbourhood of material points enclosed within a spherical radius called the horizon. The displacements and internal forces are permitted to have discontinuities and other singularities [8].

The method has been welcomed by the scientific communities in brittle fracture, damage and solid mechanics, and extensive research has been done, producing impressive results in various problems, such as crack path prediction [9–11], cracking in composite materials [12,13], bi-material strip delamination [14], J-integral calculation [15] and impact damage [16–18].

In this study, it is found that by comparing the initial/deformed pair interactions of peridynamics, the effect of pre- or residual deformation can be simulated seamlessly and efficiently. Our approach is twofold; first, energy minimization has been explored to determine the steady solution with thermal deformation. This method finds a static solution in the form of peridynamics. Secondly, using the thermally deformed structure, which is found from energy minimization, dynamic simulations are performed. Comparing the initial/deformed pair interactions, the state of pre-tension or pre-compression is clarified. The corresponding implementation of pre-tension or pre-compression will be described below.

Using the developed methods, a couple of glass structures, which have been well known for their nonlinear or interesting behaviour, have been studied and tested. In particular, we have confirmed that the strengthening of glass is a result of pre-compression, and structural reactivity, which is observed as small explosions or quick fragmentation, arises due to pre-tension.

## 2. Theory of peridynamics

Peridynamics employs a continuum model using a set of material points. Each material point represents a fractional volume of continuous media; the behaviour of each point is dependent upon nearby material points. Simply, in the case of a bond-based force field, we may presume that those material points are connected to each other using constant-stiffness (or linear) springs. The range of coupling is determined using a horizon, which corresponds to a cut-off radius, like that used in short-range potentials in molecular dynamics. Briefly, without body force, the peridynamic equation of motion for a configuration of positions **x** and time *t* will be given as
*R* is the domain of integration, **u** is the displacement vector, *ρ* is density and **f** is the response function which is described as

where *s* is the stretch between material points *x* and *x*′, *α* is a coefficient of thermal expansion of a material and *T* is the mean value of the temperature between points. Detailed description of peridynamic thermoelasticity can be found in the work by Kilic & Madenci [19]. Also the dilatational component can be included in the force field, coupling the effect of volume change. A detailed theoretical background of peridynamics can be found in the work by Silling and others [4–7]. In addition to the background theory, influence of material response, mesh sensitivity and other characteristics of peridynamics can be found in recent references [11,20].

Regarding implementation, a common approach in peridynamics is to define linked lists (or *family*) for each material point—checking the nearby material points, pairs of interactions are defined using the reference configurations, and their interactions are updated every time step. When material points exist near free surfaces, their interactions are truncated and a smaller number of neighbouring material points are found. To ensure consistent elastic energy density, appropriate scaling must be performed; this work follows the surface correction method of Kilic & Madenci [21]. When a structure deforms, pair distances change, being extended or compressed. As an extended stretch reaches a criterion (=*s*_{0}, critical stretch), we may assume that the bond pair breaks, yielding damage of the structure. Damaged material points behave like free particles when they interact with nearby material points, bearing only repulsive forces [22].

As the force field is described in terms of the pair distances, various force fields have been developed as functions of these distances. Using a linear peridynamic solid [23] as the basic force field, the effect of thermal residual deformation can be tracked by the comparison between the deformed and the un-deformed structure, as done in this work. In other words, the pre-tension or pre-compression by thermal deformation can be calculated, and dynamic failure of a structure with residual deformation can be seamlessly simulated.

### (a) Static solution via energy minimization in peridynamics

Dynamic simulation methods employ time integration over the problem domain, which bears the stress wave of the structure. As hyperbolic partial differential equations, a time step is determined from the elastic modulus and the resolution of discretization [24]. Therefore, the time length of a simulation is limited by the time-step size, and finding equilibrated structures might not be suitable in the time frame of dynamic simulations. Mediating such constraints, dynamic relaxation [25] would be adopted to find steady solutions using dynamic solvers, and a few peridynamic implementations have been reported [14,26]. However, dynamic relaxation requires damping constants or viscous forces, and this may produce ambiguity in the calculation of damping magnitude.

Motivated from the *geometry optimization* in computational chemistry [27], energy minimization can be implemented using peridynamic thermoelasticity. A common approach is to solve a Hessian matrix, using Newton–Raphson or Secant methods. But the Hessian matrix will be *N*×*N*, as *N* is the number of degrees of freedom (d.f.) in the system. As common peridynamic applications employ more than several tens of thousands of material points, the Hessian matrix will be on the order of 10^{4}×10^{4} or larger, which is computationally expensive. Therefore, we use the conjugate gradient method, which solves the nonlinear system using iteration. A detailed description and formulation is summarized in the work by Andrei [28]. Also note that our approach is equivalent to the energy minimization method by Le *et al.* [26] while our implementation is oriented for the peridynamic thermoelasticity [19].

Energy minimization employs no damping for static or steady solutions so the ambiguity of damping size can be avoided. However, energy minimization might be computationally more expensive than other methods; an appropriate method would be determined by the characteristics of the problem.

As energy minimization finds a fully equilibrated (or steady) structure, pre-tension or pre-compression naturally arises in the given structure. Such pre-deformation can be found by comparing the initial/deformed structures, and the effect of residual deformation can be employed in the dynamic simulations as shown below.

### (b) Force field with the effect of thermal residual deformation

As discussed above, a thermally deformed peridynamic structure can be found from energy minimization. The solution from energy minimization corresponds to a meta-stable structure, having pre-compression or pre-tension due to thermal deformation. As a structure deforms by thermal loading, the pair distances change accordingly. The potential curves, as shown in figure 1, show that the states of pairs move from the ground state into higher energy states due to pre-compression or pre-tension. By comparing each bond distance between the initial bond state (=*ξ*_{0}) and the thermally deformed (or meta-stable) bond state (=*ξ*_{m}) from the results of energy minimization, the amount of pre-tension/compression can be determined,
*s*_{0}) than the ground state (=*ξ*_{0}). Pre-compression is the contrary; each pair interaction has to be extended longer than any tensile state to reach *s*_{0}, increasing the resistance to failure. In terms of potential energy, pre-tension gives an additional energy of Δ*E* to the pair interaction, as shown figure 2*a*, and this energy needs to be released or discharged when the bond pair fails. In the case of pre-compression, the energy given by thermal deformation would be consumed when each pair is extended to reach *s*_{0}, yielding no stored energy released.

In particular, the release of the stored residual energy (=Δ*E*) is believed to be responsible for the dynamic damage evolution of strengthened glasses. As discussed above, commonly strengthened glasses show higher strength against external loading but when failure initiates, the damage rate increases faster than plain glasses, yielding more fragmentation, or even pulverization like an explosion.

When the stored energy is released by the failure of pre-tensile bond pairs, we may assume that the stored energy will be converted into the kinetic energy of the material points of the broken bonds, instantaneously. The strain energy is divided evenly to both material points,
*m*_{i}*v*_{i}=*f*_{i}*δt* then

## 3. Applications

Using the developed methods, the following models have been tested. By thermal treatment, pre-tension/compression are configured in the structure and are well known to show highly nonlinear behaviour when damaged. Therefore, the following applications would be good candidates to test the developed methods qualitatively.

As peridynamic simulations are performed, an appropriate critical stretch (=*s*_{0}) must be determined, using the sizes of grid, horizon and material properties. Prior to this work, we have performed many sets of tension simulations using dog-bone style specimens. Using the grid size of 0.2–0.5 mm, each elastic strain limit has been compared with available experimental data [30]. Soda lime glass shows a wide range of strength for different compositions and heat treatments; we have selected a maximum value for an extreme case like Rupert’s drop and an intermediate value for the Bologna vial, as shown below. Note that the main motivation for using the strength-based critical stretch as opposed to a toughness-based value is to allow the large tensile stretch caused by thermal deformation.

### (a) Prince Rupert’s drop

This odd glass bead has been reported in the sixteenth century at the Royal Society of Britain [31]. It has a tadpole shape—a semi-spherical head and a long-thin tail. It is formed by dripping a melted glass rod into water, which will cool down the surface of the melted glass bead quickly, forming high compressive stress in the surface layer. The inside of the glass bead remains relatively hot and expanded. As the surface is cooled and contracted, the inside gains tensile pre-stress. Finally, the Rupert’s drop configures opposing stress states on the inside and surface, showing highly nonlinear behaviour. The head has extremely high resistance to loading, while the tail is highly fragile; any slight load even a pinch by fingers breaks the tail, disintegrating the entire structure like an explosion [32]. In particular, this *reactive* behaviour, showing explosive disintegration, is one kind of dynamic fracture phenomena and has attracted scientific research [33,34]. This historically famous glass structure has been recognized as an ancestor of modern tempered or strengthened glasses [35].

The peridynamic model has been produced using uniform grids of 0.2 mm grid distances, yielding 106 553 material points. The actual shape of the Rupert’s drop looks like a tadpole, as it is made from the molten glass bead. In this modelling, we simplified the shape as a combination of a hemisphere and an elongated cone. The head has a shape of a hemisphere with a radius of 5 mm, and the tail has a radius of 1.7 mm and is extended up to 25 mm, yielding the length of the glass bead as around 30 mm. Between the head and tail parts, a second-order polynomial has been used to make a smooth geometry. Note that even the narrow region has more than 15 material points to represent the elasticity of the structure. In order to configure pre-tension and pre-compression in the glass bead, the skin of the bead has been set at Δ*T*=−400 K and the core at Δ*T*=400 K. The figure 3*a* shows the configured temperature distribution. Even though an actual Rupert’s drop is produced from a molten glass drop, which will have a temperature around 1500°C, residual deformation would decrease due to annealing. We configure these boundary conditions as a numerical experiment, in order to simulate the behaviour of Rupert’s drop qualitatively.

Using the energy minimization mentioned above, the glass bead deforms accordingly. As the inner/surface of the bead is configured as a higher/lower temperature, each region bears tension/compression after the energy minimization. Figure 3*b* shows the resulting pressure contour along the bottom half cross-section.

Using the thermally deformed structure, several impact tests have been performed. They are dynamic tests and no fixed boundary condition has been used. The employed material properties are summarized in table 1. Note that material failure properties of the impactor is not configured, making it have an infinite elasticity. Also the critical stretch of the glass has been extended artificially, because the value based on toughness does not allow the extension made by the thermal deformation. In the Rupert’s drop tests, 3.0×10^{−3} has been used as the critical stretch. From tension simulation tests using dog-bone specimens, this critical stretch has yielded 0.002 maximum strain which corresponds to 145 MPa strength, as an extreme condition [30].

Modelling an impactor as a 3 mm radius sphere with 13 997 material points, the impactor hits the glass bead at a velocity of 15 m s^{−1}. For comparison, we have tested a plain glass bead, which has no pre-deformation, and the difference in behaviour is discussed.

Firstly, the impact sphere is applied to the tail of the Rupert’s drop. As shown in figure 4*a*, we can observe the explosive or reactive behaviour of the Rupert’s drop. Initial damage is found near the impact site, and the damage propagates into the whole body. Figure 4*b* shows how the damage evolves inside of the Rupert’s drop. The impact initiates the damage in the core of the bead, where it bears high pre-tension facilitating failure. Then the damaged domains influence nearby pre-tension domains, and the damage cascades. As the inner core is under pre-tension, damage evolves quickly and finally the damage influences the surface domain as well, which is under pre-compression. As discussed above, this *structural reactivity* is produced when the stored strain energy due to thermal deformation is released, and the implementation of released energy describes such behaviour successfully.

For further analysis, the speed of damage propagation has been estimated. Using figure 4*b*, the speed of the damage front is calculated to be around 4.5 km s^{−1}, which is almost three times larger than the experimentally reported value (≈1.5 km s^{−1}) [38]. This is closer to the elastic sound wave velocity of soda lime glass (5.3 km s^{−1}), because no damping and no friction is included in the interactions between damaged material points. Bobaru & Hu [11] reported the effect of horizon size over the cracking, and different grid sizes and horizon sizes as 80–150% of the employed configuration have been tested over the same model but any significant change of propagation speed has not been observed yet. For better description of post-fracture or post-damage behaviour, we expect that micro-damage evolution or internal frictions would be required. They are beyond the current scope and will be studied in the future.

Results from the plain glass bead is shown in figure 4*c*, showing a brittle breakage of the tail part. As no pre-deformation is configured, no reactive behaviour is observed, yielding a passive failure state. Compared with the results of the pre-deformed Rupert’s drop, the difference in behaviour is clear, showing the importance of structural reactivity.

Secondly, now the impactor hits the head of the glass bead at 25 m s^{−1}. The damage comparison along a cross-section is shown in figure 5. Note that the impact velocity in this test is higher than the tail test in order to produce more damage. Compared to the tail, the head is physically larger and resists the loading. From figure 5*a*, no damage is found at the impact site of the Rupert’s drop. When a plain glass bead is tested, localized indent damage is found near the impact site, As expected, the Rupert’s drop shows higher strength than a plain glass bead in the head impact, because it has a deep pre-compression layer in the head.

As discussed above, the release of the stored energy is found to contribute to the pulverization or quick disintegration of the structure. This *explosion* is not from a chemical reaction, but due to the discharge of the *structurally* stored energy. Therefore, this *structural reactivity* is of great importance in the dynamic damage evolution of strengthened glasses. Rupert’s drop is an example of an extremely strengthened glass, and these methods were successful to represent the nonlinear and reactive behaviour quantitatively. As a second application, the Bologna vial has been tested as a moderately strengthened glass structure.

### (b) Bologna vial

This vial has a shape of a bulb or flask, and thermal treatment is done differently on the inside and outside of the vial, similar to the Rupert’s drop. While the outside of the bottle is strengthened by rapid cooling, the inside experiences slower cooling, forming a tensile state. Therefore, the outside of the bottle shows high strength against impact or scratch while the inside is extremely weak, shattering the bottle with a small amount of force. [39]. Due to this nonlinear behaviour, this bottle has been used in magic shows or effects, called the *Devil’s Flask*.

There could be various shapes of Bologna vials but an egg shape bottle has been adopted in this study, which is commonly shown in the Corning Museum of Glass. Simplifying the shapes using ellipsoidal functions, a peridynamic model has been built. Between surfaces of *x*^{2}/10^{2}+*y*^{2}/10+*z*^{2}/15=1 and *x*^{2}/13^{2}+*y*^{2}/13+*z*^{2}/25=164 000 material points are allocated and an end tab is configured as truncation at *z*=10. This end tab allows an impactor to enter the bottle to access the inner side for impact. The configured shape yields an outer ellipsoid with 13–25 mm radii and an inner ellipsoid of 10–15 mm radii. The material points are produced between these two ellipsoids. A uniform grid with material points every 0.5 mm has been used and temperature conditions were configured as Δ*T*=−400 to 250, on the outside and the inside of the bottle. Figure 6 shows the distribution of temperature and the produced pressure, along the bottom half cross-section of the vial.

The employed material properties are same as the values used in the Rupert’s drop tests, but the critical stretch has been adjusted to 2.3×10^{−3}, yielding approximately 0.0017 maximum tensile strain. For impact tests, an impactor was modelled using a long-rod, with a 100 mm length and 3 mm radius on the hemispherical tip. This impactor hits the apex of the glass bottle at 20 m s^{−1} through the open hole and from the outside.

Impact test results through the open hole are summarized in figure 7. From overall shapes, completely different damage evolution is observed. In particular, along the half cross-sections, the source of the different damage morphology is clarified. For a plain glass bottle, damage is found near the impact site, leaving dent damage. However, the Bologna vial, which has thermal pre-deformation, shows a complete shattering of the structure with pulverization or complete damage near the impact site. Similar to the Rupert’s drop, the pre-tension inside of the bottle promotes more damage, and much greater damage evolves, cracking the whole structure.

The same configuration and material properties are used for the impact tests from the outside. As shown in figure 8, the impactor hits the apex of the glass bottle from the outside, and each damage pattern is compared. Even though both of the bottles show similar resistance against the impact, different local damage is observed. The plain glass bottle leaves a slight dent-like damage near the impact site, similar to the previous result, but the impact on the Bologna vial shows no visible damage, as the outside is strengthened by pre-compression. The damage magnitude of the Bologna vial shows one order lower damage rate than a plain glass bottle when the apex is impacted from the outside.

Similar to the Rupert’s drop, the Bologna vial shows the effect of thermal residual deformation. With pre-tension, structures are more vulnerable and failure is promoted. As damage initiates, the stored energy is released and *structural reactivity* appears. Even though the Bologna vial shows less reactivity than the Rupert’s drop, it still shows some pulverization near the impact site and fragmentation of the structure. With pre-compression, the structure is more able to resist against the external loading and shows higher strength against impact. The developed methods simulate these nonlinear or reactive behaviours successfully, revealing the role of residual deformation in dynamic damage evolution.

## 4. Discussion

The choice of critical stretch and grid-size in this work was chosen based on uniaxial strength calibration for the critical stretch and the grid size determined for convenience due to the large model dimensions. This causes a discrepancy between the fracture toughness and the strength of the material. Consider the equation for the critical stretch, which can be found in most published literature on peridynamics,
*s*_{0} to be the critical strain measured from experiments and use the material fracture energy G_{IC} then the horizon, *δ*, is constrained. In this work, the Prince Rupert’s model has *s*_{0}=3×10^{−3} and *G*_{IC}=7.7 J m^{−2}. This would require *δ*=9.85 μm with grid size of a third of that (=2.4 μm), yet for convenience, due to the large size of the model, and a grid size of 200 μm was used. This made the toughness 469 J m^{−2}. Thus the examples shown in this paper are for instructive purposes only and are not a true representation of the material. In practice, fitting the critical stretch to the material strength is good for studying fracture initiation when no obvious flaw exists in the material, but for modelling crack propagation, an accurate value of toughness is required, and for those models an initial flaw is highly suggested.

## 5. Conclusion

In this study, we have employed energy minimization for thermal deformation and revised the force field to implement the effect of thermal residual deformation in glasses or brittle structures. Pre-compression increases the resistance to tensile failure, requiring larger extension of each bond pair, and pre-tension promotes tensile failure, releasing stored energy.

Using the developed methods, dynamic failure of thermally strengthened glasses has been simulated, confirming the strengthened effect due to pre-compression and structural reactivity by pre-tension. Within a single dynamic solver, the effect of thermal residual deformation could be included in the dynamic damage evolution, enabling seamless simulations of brittle structures under extreme conditions.

As applications, the well-known Prince Rupert’s drop has been modelled and tested. Thermally deforming a glass bead, pre-tension/compression has been configured in the inside and the skin of the structure, respectively. The tail part has a thin skin of pre-compression, and it is more fragile than the head, which has a thick skin of pre-compression. As the tail of the bead is impacted, the internal damage grows quickly, cascading damage into most of the structure. The thermal deformation energy is stored structurally, and the stored energy is released as the pre-tension material points break, yielding pulverization as strong structural reactivity. As the head is impacted, it resists the failure, showing higher strength than the tail or even a plain glass bead. Even though well-known nonlinear properties, such as strengthening in the head and reactivity in the tail, are simulated successfully, further effort on the modelling of interactions between damaged material points would be required. Modelling of micro-damage or introducing friction would be done in the future.

The second application, a Bologna vial has been configured. As thermal residual deformation is implanted, the behaviour of dynamic damage evolution has been tested. With a long-rod impactor, the inside/outside of the vial was hit. As the pre-tension exists on the inside of the vial, tensile failure is facilitated and the glass is shattered. The impact test on the outside shows the contrary: higher strength than a plain glass bottle. Compared to the Rupert’s drop, the density of the stored thermal energy is relatively lower and strong reactivity is not observed. But the stored energy is still released and influences the dynamic failure of the structure, cracking the bottle as structural reactivity.

Through this study, the method and implementation to simulate the effect of residual deformation have been developed. Enabling seamless simulations of dynamic damage evolution with peridynamics could be a breakthrough in the research of glass and other brittle solids.

## Data accessibility

This work does not have any experimental data.

## Authors' contributions

B.J. developed the methods, the peridynamics code, and models, analysing the results. R.S. and I.A. revised the mechanism of the strengthening and reactivity, consulting the results. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This work was supported by Corning Incorporated.

## Acknowledgements

We thank Sam Zoubi at Corning Incorporated for giving us this unique opportunity to study peridynamics. Also, we thank Prof. Erdogan Madenci from University of Arizona for helpful discussion on Peridynamics. ByoungSeon Jeon appreciates the help and support from Chris Wilson, Charles E Forsyth III, James E Ross, Matthew McKenzie and Sergey Buduchin at Corning Incorporated.

- Received April 8, 2015.
- Accepted October 27, 2015.

- © 2015 The Author(s)