Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2022.34.0098 Acta Polytechnica CTU Proceedings 34:98–104, 2022 © 2022 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague PHASE-FIELD DAMAGE OF LAMINATED GLASS IN EXPLICIT DYNAMICS Jaroslav Schmidt∗, Tomáš Janda Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Praha 6, Czech Republic ∗ corresponding author: jaroslav.schmidt@fsv.cvut.cz Abstract. Laminated glass is aesthetic but also functional modern material, which has spread from the automotive industry to the civil engineering in the form of design as well as structural elements. Proper design of this material requires understanding not only its intact phase, but also its post-breakage behavior. To describe crack initiation and propagation the phase field damage model seems to be suitable tool despite the fact that this approach requires a fine mesh to sufficiently interpolate sharp discontinuity. In this contribution we investigate possibilities of explicit phase field dynamic model with specific use on glass plates, especially applied for dynamic impact load. The dependence of the result on residual stiffness term is examined and incorrect crack speed is observed. Usability of phase-field model is still uncertain. Keywords: Phase field damage model, laminated glass, dynamic damage, glass fracture. 1. Introduction The so called phase-field damage model firstly pre- sented in [1] is nowadays very popular approach and it is widely studied. The disadvantage of such a model is the need for a fine mesh to simulate brittle failure, but the main advantage remains the prediction of crack initiation and branching without additional ad hoc criteria. These phenomena are predicted within a variational principle. Despite the fact that it was academicly investigated in dynamic regime [2], there is still a few articles focused on application on real structures. This paper focuses on such an analysis by modeling the cracking of glass under dynamic loading. The dynamic model including impact is briefly pre- sented in Section 2. The phase-field model formulation follows partially from [3] and [4]. For comprehen- sive overview see [5]. The implementation details are summarized in Section 3 followed by FEM model de- scription in Section 4. The first numerical results are benchmarked in Section 5 and finally the main results are presented in Section 6 and concluded in Section 7. 2. Phase-field model We suppose the elastic damageable material occupying space Ω is loaded by external traction forces t on part Ωt and the current body deformation is described by displacement field u(x) and damage field d(x). The latter one represents the relative measure of damage as a transition of the intact material d = 0 to a fully developed crack d = 1. For dynamical systems, the common principle of minimum potential is replaced by Hamilton’s principle of stationary action and fields u(t) and d(t), t ∈ 〈0,T〉 are obtained by minimization of action A(u,u̇,d) = ∫ T 0 L(u,u̇,d) dt, (1) where L is so called Lagrangian which form is postu- lated as L(u,u̇,d) = K(u̇) − Ψe(u,d) − Ψs,l(d) + P(u), (2) where K is kinetic energy defined as quadratic form of displacement K(u̇) = ∫ Ω 1 2 ρu̇ · u̇ dx. (3) This formulation assumes that kinetic energy is not affected by damage, therefore the functional is d- independent. The second term in equation (2) is elastic potential. It can be additively decomposed on part Ψ+ which is degraded by damage through function g(d) and intact part Ψ−, i.e. Ψe(u,d) = ∫ Ω g(d)ψ+(u) + ψ−(u) dx, (4) where ψ• represents elastic energy density of Ψ•. In what follows we always assume that g(d) = (1 −d)2 + kres, where kres is small (kres << 1) auxiliary residual relative stiffness which stabilizes the calculations. The third part of (2) is dissipated energy Ψs,l which is approximated by following regularization Ψs,l = Gf cα ∫ Ω α(d) l + l|∇d|2 dx, (5) where Gf is the critical energy release rate and l is length of the process zone. For brittle material, l can be treated as numerical parameter with value as 98 https://doi.org/10.14311/APP.2022.34.0098 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 34/2022 Phase-field damage of laminated glass in explicit dynamics small as possible, but not smaller than element size. The function α and associated coefficient cα drives the phase-field model type. Two common options are commonly used in literature: α = d with cα = 8/3 and α = d2 with cα = 2 however we restrict to the first one. Finally, the last ingredient is external work of trac- tion forces t defined as P(u) = ∫ Ωt t ·u dx. (6) Minimization of the action integral (1) with addi- tional per partes integration leads to common conti- nous governing equations for displacement field div ( g(d) ∂ψ+ ∂u + ∂ψ− ∂u ) = ρü (7) supplemented by damage-evolution inequality Gf cα ( α′(d) l − 2l∆d ) = −g′(d)ψ+, ḋ > 0 (8) Gf cα ( α′(d) l − 2l∆d ) > −g′(d)ψ+, ḋ = 0 (9) and approriate boundary conditions. Strong-form equations (7)-(9) have a generally unknown solution, therefore some numerical method must be used. 2.1. Plate models The model and its properties are determined by the choice of dissipation function α(d), degradation func- tion g(d) and by active/passive elastic energy densities ψ±. Four common choices for the latter one are pre- sented in this subsection. Despite the fact that we restrict our attention to model P-VD, others are pre- sented. This is because we intend to investigate them in the future research. 3D spectral decomposition It is assumed that material behaves linearly elastically in pre-breakage stage, therefore it can be described by two independent Lame’s coefficients λ and ν. Since the damage variable starts evolving, the only positive part of energy must be degraded, therefore ψ±(u) = 1 2 λ〈±tr(ε)〉2 + µε± : ε±, (10) where ε is conventional strain tensor defined as sym- metric part of displacement field gradient ε = ∇su. Positive and negative component of strain tensor fol- lows from the principal strains εi and the dyadic product of eigenvectors pi of ε, thus ε± = ∑ i ±〈±εi〉pi ⊗pi. (11) We denote this model as 3D − SD 3D volumetric decomposition The second com- mon choice of decomposition is volumetric-deviatoric split, which assumes that only compression part of volumetric strain component remains intact, thus ψ+(u) = 1 2 K〈tr(ε)〉2 + µεD : εD (12) ψ−(u) = 1 2 K〈−tr(ε)〉2, (13) where K stands for bulk modolus and εD is deviatoric part of strain ε. This model is denoted as 3D − VD. Mindlin-Reissner plate Besides full 3D model, we investigate the spatialy reduced Mindlin-Reissner plate model. This model assumes linearly distibuted strain and stress field across thickness and kinematics is driven by middle-surface longitudinal displacement field u, out-of-plane scalar deflection w and vectorial tensor of cross-section rotations ϕ. Based on the linearity assumption the longitudinal strain ε and shear strain γ can be established as ε(xs,z) = ∇su(xs) + ∇sSϕ(xs)z, (14) γ(xs) = 1 2 (∇w(xs) + Sϕ(xs)) , (15) where z is out-of-plane coordinate meanwhile xs = {x,y} is in-plane coordinate system of middle-surface Ωs. The S is auxiliary matrix S = [ 0 1 −1 1 ] . (16) Mindlin-Reissner plate model ignores normal stress across thickness, therefore in each infinitesimal layer plane-stress situation occurs. Therefore stress fields are defined as σ(ε) = DPS : ε, (17) τ(γ) = 2νI : γ, (18) where DPS is plane-stress elastic tensor defined again through Lame’s coefficients λ and ν and I is identity tensor. The last ingredient is thickness integration for spatial reduction. It is impossible to perform it in closed form therefore numerical integration is employed and active/pasive elastic energy become Ψ+ = 1 2 ∫ Ωs g(d) (∑ i ∆ziεi : σ+i + ξhγ ·τ ) dxs, (19) Ψ− = 1 2 ∫ Ωs ∑ i ∆ziεi : σ−i + (1 − ξ)hγ ·τ dxs. (20) Summation represents dividing thickness into layers i, each represented by in-plane stress σi and strain εi with sub-thickness ∆zi. Parameter ξ decides what amount of shear is degraded; ξ = 1 means shear is fully degraded meanwhile ξ = 0 is shear damage-free 99 Jaroslav Schmidt, Tomáš Janda Acta Polytechnica CTU Proceedings state. Still the stress tensor σ±i must be somehow decomposed. Based on this decomposition this model is denoted as P − SD if spectral decomposition is used and is denoted P − VD for volumetric-deviatoric split. The degradation function still remains g(d) = (1 − d)2 +kres, but we found that it is appropriate to select different parameter kres for shear, which is denotes by ksres. 2.2. Discrete dynamic model First of all we suppose that the solution is evaluated only in discrete time instants 0 = t0 < t1 < t2 < ... < tN = T selected from the original time interval 〈0,T〉. The action integral (1) is approximated based on this time discretization as N−1∑ i=0 ∆tiLd(ui,ui+1,di,di+1), (21) where common notation ui = u(ti) and di = d(ti) is used and also ∆ti = ti+1 − ti = ∆t is supposed to be constant. Function Ld is discretized Lagrandian, which approximate the action integral on given sub- interval. For purpose of this article, relatively simple expres- sion leading to central difference scheme is used, i.e. Ld(ui,ui+1,di,di+1) = L(ui, ui+1 −ui ∆t ,di). (22) As a result the original minimization in continuous time domain fell apart to individual minimizations with respect to each displacement ui and damage di, where i = 0, . . . ,N. The stationarity condition for displacement field ui leads to discrete Euler-Lagrange equation in form ∂Ld(ui−1,ui,di−1,di) ∂ui + ∂Ld(ui,ui+1,di,di+1) ∂ui = 0, (23) which can be expressed implicitly in weak form δK(δu) − δΨe(δu) + δP(δu) = 0,∀δu (24) where δK = 1 ∆t2 ∫ Ω ρ(ui+1 − 2ui + ui−1) · δu dx, (25) δΨe = ∫ Ω ( g(di) ∂ψ+(ui) ∂ui + ∂ψ−(ui) ∂ui ) · δu dx, (26) δP = ∫ Ω t · δui dx. (27) The presented weak form can be further modified to obtained governing equations in strong form, but we employ finite element method for solving spatial distribution of ui which requires weak form as an input. The governing weak form for damage is identical with quasi-static version. The stationary condition δA(u,u̇,d,δd) = 0, ∀δd (28) becomes∫ Ω Gf cα ( α′(d) l + l∇d∇δd ) + g′(d)ψ+ dx = 0, ∀δd (29) with damage irreversibility condition di(x) ≤ di+1(x), ∀x, i = 0, . . . ,N − 1 (30) This approach therefore induces the following ex- plicit non-incremental dynamic procedure: (i) solve explicitly displacement in time ti by linear problem (24), (ii) solve damage field di by variational ineqaulity problem (29)+(30), (iii) increment time. 2.3. Nonlinear Hertz force The impact of steel impactor with weight mimp is in- cluded in the model through Hertz-law. The impactor is characterized by scalar displacement uimp, which is binded to plate by nonlinear contact force F(u,uimp) = k〈uz(ximp) −uimp〉3/2, (31) where ximp is impact position on body Ω, uz is out-of- plane component of u = (ux,uy,uz) and k is contact stiffness. Although the force is nonlinear the pseudo potential of this force can be still founded. This potential PHertz = 2 5 k〈uz(ximp) −uimp〉5/2 (32) is added to external work (6) and it induces force (31) applied to impact point ximp and negative force (31) applied to impactor. Additionaly the kinetic energy must be enhanced by Kimp = 1 2 mimpu̇impu̇imp (33) with impactor displacement rate u̇imp. This part of kinetic energy is again approximated by discretized Lagrangian and discrete Euler-Lagrange equation in- duces following additional member of weak form (26), i.e. δKimp = mimp uimp,i+1 − 2uimp,i + uimp,i−1 ∆t2 δuimp (34) 3. Implementation For solving spatial distribution of ui, the common finite element method (FEM) is used. The displace- ment field is approximated by nodal values ri through selected basis functions. Since the FEM is well estab- lished, the details of method is not presented here. To handle simultaneously deflection of elastic body and 100 vol. 34/2022 Phase-field damage of laminated glass in explicit dynamics L1 = 1.5 m L 2 = 1. 5 m h = 0.02 m E = 72 GPa ν = 0.22 ρ = 2500 kg/m3 x y Figure 1. Setup of investigated task. impactor position, the all matrices and vectors are ex- tended by one additionaly degree of freedom and final object is marked by tilde. For example displacement vector r̃i or mass matrix M̃ mean r̃i = {ri,uimp}, M̃ = [ M 0 0 mimp ] . (35) With this in hand, the displacement problem (24) after FEM assemble process and localization becomes 1 ∆t2 M̃(r̃i+1 − 2r̃i + r̃i−1) + K̃r̃i + F̃ (r̃i) = 0̃, (36) which is linear in parameter ri+1. The matrices K̃ and M̃ are stiffness and mass matrices. Presentation of their precise form goes beyond the scope of this paper. The assembling and localization of damage problem (29) leads to linear form Hdi + fd = 0, (37) but let us remark several comments. First of all, the discretized damage field di still undergoes irreversibil- ity condition di � di+1 i = 0, . . . ,N − 1, (38) where � operator stands for element-wise compari- son. From that reason the linear problem (37) must be solved by nonlinear solver based on variational in- equalities to satisfy condition (38). The second remark is about matrices K̃, H and damage-driving force fd. These are ui resp. di nonlinearly dependent and must be assembled in each time steps. It significantly slows down the calculations. 4. Setup The all simulations are performed on the same task – unsupported square glass plate hit by steal impactor. This setup follows from real experimental arrangement, see [6]. The plate has dimensions 1.5 × 1.5 m and thickness 0.02 m, but only one quarter is modeled due to symmetry, see Figure 1, where material parameters are also presented. The contact stiffness is defined as k = 4 3 √ R 1−ν2 E + 1−ν2imp Eimp , (39) where R = 0.05 m is impactor head radius and Eimp = 210 GPa and νimp = 0.3. Further weight of impactor is mimp = 48.7 kg and initial velocity-free height of impactor is 0.3 m. Finally the phase-field parameter l is choosen as element length and Gf = (8lf2t )/(3E), where ft = 55 MPa is tensile strength. 5. Benchmarks Phase-field damage model requires fine mesh with lower-order elements rather than higher-order ele- ments with lower resolution. It is contrary to plate theory, where lower-order elements exhibit shear lock- ing and must be somehow enhanced to overcome this issue. The plate element enrichment leads to addi- tional degrees of freedom which reduces computability and the ability to do repeated calculations. The finer resolution is more important in this analysis than exact shear implementation, therefore suitability of "classic" low-order plate elements is investigated in this benchmark. The first benchmark is behavior of impacted plate in elastic regime. The task setup presented in Sec- tion 4 is used and the response is calculated in time interval 〈0; 0.01〉 s with time step ∆t = 4 · 10−7 s. Resolutions of mesh 50 × 50 and 100 × 100 elements are investigated with two different element typolo- gies: triangles and squares. The reduced integration is employed for shear terms, which enhanced rectangle elements, but triangles still locking. The reason for not using square elements is because they exhibits un- desirable zero-energy modes. The results of each case is presented in Figure 2, where evolution of impactor acceleration is plotted. This is one of the interesting investigated quantity, which can be experimentaly val- idated in future research. Also the different behavior of two cases stands out more than in impactor dis- placement function. It is evident from Figure 2 that behavior is slighly different for triangles in case of mesh 50 × 50, but both line are comparably the same for case 100 × 100 and shear locking does not cause significant deviation. This statement holds for given time resolution and fraction h/L of course. Never- theless this setup follows from real experiment and we assume that results are adequate. The deflections behaves similarly and do not suffer for larger deviation for mesh 100 × 100. However, the usage of triangle elements exhibit par- asitic behavior in the damaged state as is presented in second benchmark. The task is the same as in previous benchmark, but initial crack is incorporated 101 Jaroslav Schmidt, Tomáš Janda Acta Polytechnica CTU Proceedings 0 0.2 0.4 0.6 0.8 1 ·10−2 −600 −400 −200 0 Time t [s] Im p. ac ce le ra ti on [m /s 2 ] 50x50 triangles rectangles 0 0.2 0.4 0.6 0.8 1 ·10−2 −600 −400 −200 0 Time t [s] 100x100 triangles rectangles Figure 2. Evolution of impactor acceleration for triangles (solid) or squares (dashed) mesh plotted for resolution 50 × 50 (red) and 100 × 100 (blue). 0 0.2 0.4 0.6 0.8 0 2 4 6 ·10−3 x coordinate [m] D efl ec ti n [m ] rectangles triangles triangles-full Figure 3. Deflection distribution on line {(x, y)|y = 0.75}. on plate in this case. The initial crack is implemented through Dirichlet boundary condition d = 0.99 on the line {(x,y)|ζ > |x−L1/4|}, where ζ parameter is cho- sen such that the damage d = 0.99 is enforced at least over whole width of one element. The simplified but still sufficient analysis is performed in this benchmark using so called hybrid approach [7] for phase-field. Furthermore, the damage d and updated matrix K̃ from (36) is calculated in each third step to speed up analysis. This approach allows some minor wave prop- agation through crack, but wave dispersion through crack is still significant. The result also depends on the parameter ξ from (19)-(20). First of all, let us investigate case ξ = 0, therefore shear is not degraded. The results are plotted by solid lines in Figure 3, where plate deflection along line {(x,y)|y=0.75} is plotted in time instant t = 0.003 s. The square-topology mesh M1 M2 M3 Figure 4. Three investigated meshes. behaves reasonable, but triangle-topology mesh does not stop the wave on crack interface. It is caused probably due to element locking. On the other hand, triangles with ξ = 1 behaves almost like squares with ξ = 0 as a result of shear degradation by damage. It significantly reduces the shear locking artifacts. Based on these conclusions, the triangle-based mesh with shear reduced by degradation function is used in the next section. 6. Results The calculations are time-consuming and therefore only preliminary results of model P-VD are presented in this contribution. The model is still under inves- tigation and there are still many unresolved issues. The same setup, defined in Section 4, is used, i.e. the steel impactor hit the thin glass plate. Three different meshes M1, M2 and M3 are tested to detect possible mesh-dependency, see Figure 4. The model with fully degraded kinematic variables, i.e. ξ = 1.0,kres = ksres = 0.0, does not behave cor- rectly. The impactor hit the plate causing damage lo- calization under impactor. The plate has zero stiffness and impactor flies through fully damaged elements and no longer influences the rest of plate. As a result, the cracks are not evolving further. This behavior is illustrated in Figure 5 by solid line, where again the distribution of plate deflection along {(x,y)|y = 0.75} 102 vol. 34/2022 Phase-field damage of laminated glass in explicit dynamics 0 0.2 0.4 0.6 0.8 0 1 2 3 ·10−3 x coordinate [m] D efl ec ti n [m ] ksres = 0.0 ksres = 0.05 Figure 5. Deflection distribution on line {(x, y)|y = 0.75}. (a) (b) (c) (d) Figure 6. Damage distribution on meshes M1 (a, c) and M2 (b, d) for resolution 50 × 50 (a, b) and 100 × 100 (c, d). The blue represents intact material whereas red represents d close to 1. and in time t = 1.4 ms is plotted. The example is eval- uated on mesh M3. The sort of way-out can be shear coupling in damage state, i.e. kres = 0.0,ksres 6= 0. It prevents large slope of deflection curve, see dashed line in Figure 5, where ksres = 0.05. The observed behavior is more pleasant and it leads to crack evolving as pre- sented below. Moreover the crack is not thickening after additional load. In contrast with that, the case kres 6= 0 causes thickening. Mesh dependency The meshes M1, M2 and M3 in two resolutions 50 × 50 and 100 × 100 was used to investigate mesh dependency. The setup is again the same with constant time step ∆t = 2·10−7 s and final time instant tN = 0.01 s. The results of damage fields are plotted in Figure 6, where (a) is mesh M1 50×50, (b) is mesh M2 50×50, (c) is mesh M1 100×100 and finally (d) is mesh M2 100 × 100 all evaluated in end time. The mesh M3 behaves similarly as M1 and is not presented here. There is evident that crack topologies are different on coarse mesh, however this discrepancy is removed Figure 7. Visualisation of plate deflection. after refinement. The slight mesh dependency can still remain in phase-field models for unstable solutions, but it still remains more or less objective approach in contrast with other techniques. Time refinement The damage solution (d) in Fig- ure 6 and its visualisation of deflection in Figure 7 is therefore appropriate and objective without any further improvement during time refinement as seen in the graph 8. The red and blue line are almost iden- tical. An interesting parasitic behavior appears when the residual stiffness is prescribed also to the bending members (kres = 0.01), see black line in graph 8. The residual term causes unwanted harmonic oscillations between impactor and plate. 0 2 4 6 8 ·10−3 −400 −200 0 Time t [s] Im p. ac c. [m /s 2 ] ∆t = 2 · 10−7, kres 6= 0 ∆t = 1 · 10−7, kres = 0 ∆t = 2 · 10−7, kres = 0 0 0.2 0.4 0.6 0.8 1 ·10−3 −350 −300 −250 −200 Time t [s] Im p. ac c. [m /s 2 ] DETAIL Figure 8. Evolution of impactor acceleration for two different time steps (red and blue) and for case kres 6= 0 (black) with detail below. Crack speed The main drawback of presented ap- proach is its inability to correctly predict crack ve- locity. It is evident from near-the-failure evolution of 103 Jaroslav Schmidt, Tomáš Janda Acta Polytechnica CTU Proceedings t = 0.00492 s t = 0.00522 s 0.5135 m 0.7210 m Figure 9. Crack development in two different time instants for mesh M1 100 × 100. crack in Figure 9. The crack length increment divided by time increment gives us the numerically predicted crack speed vPF = 692 m/s. However the crack speed of glass is about vreal ≈ 1500 m/s [[8]], which is still twice as much. 7. Conclusion The explicit variationally consistent dynamic phase- field damage model was presented in this contribu- tion with additional implementation details and pre- liminary results. The plate model with volumetric- deviatoric decomposition was successfully imple- mented, but it still exhibits some issues – especially wrong crack velocity and no crack branching. This two facts may be related, because crack branching is caused by large crack velocity. It is under further in- vestigation as well as implementation of other models (particularly P-SD). Acknowledgements This publication was supported by the Czech Science Foundation, the grant No. 19-15326S and by the Grant Agency of the Czech Technical University in Prague, grant No. SGS21/038/OHK1/1T/11. References [1] B. Bourdin, G. A. Francfort, J.-J. Marigo. The variational approach to fracture. Journal of elasticity 91(1):5–148, 2008. [2] M. Hofacker, C. Miehe. Continuum phase field modeling of dynamic fracture: variational principles and staggered fe implementation. International journal of fracture 178(1):113–129, 2012. [3] Y. Shen, M. Mollaali, Y. Li, et al. Implementation details for the phase field approaches to fracture. Journal of Shanghai Jiaotong University (Science) 23(1):166–174, 2018. [4] K. Pham, H. Amor, J.-J. Marigo, C. Maurini. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics 20(4):618–652, 2011. [5] J.-Y. Wu, V. P. Nguyen, C. T. Nguyen, et al. Phase-field modeling of fracture. Advances in Applied Mechanics 53:1–183, 2020. [6] A. Zemanová, P. Konrád, P. Hála, et al. Experimental study on the gradual fracture of layers in multi-layer laminated glass plates under low-velocity impact. arXiv preprint arXiv:210712151 2021. [7] M. Ambati, T. Gerasimov, L. De Lorenzis. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Computational Mechanics 55(2):383–405, 2015. [8] F. Barstow, H. Edgerton. Glass-fracture velocity. Journal of the American Ceramic Society 22(1-12):302–307, 1939. 104 Acta Polytechnica CTU Proceedings 34:98–104, 2022 1 Introduction 2 Phase-field model 2.1 Plate models 2.2 Discrete dynamic model 2.3 Nonlinear Hertz force 3 Implementation 4 Setup 5 Benchmarks 6 Results 7 Conclusion Acknowledgements References