Electric vehicles can have only a minor role in reducing transport's energy and environmental challenges

1.   Introduction

Fracture mechanics is a fundamental field of study in mechanical sciences. Over the years, extensive research has focused on studying fracture processes under static loading conditions, providing valuable insights into crack initiation, propagation, and structural failure. However, in many real-world scenarios, materials are subjected to dynamic loading conditions, where cracks can propagate at high speeds and exhibit unique fracture behaviors ^{[7,39,43,45,64,73,77,80]}. Understanding the mechanics of dynamic fracture is hence crucial for enhancing the safety and reliability of structural components and advancing the design of resilient materials. The mechanics underlying dynamic fractures presents many more intricacies than underlying quasi-static fracture ^{[58,59,60,61]}. Experimental studies, however, show its practical relevance^{[6,17,20,23,25,68,79,81]} and also elucidate its loading rate-dependent behavior ^{[18,24,26,27,30,40,44,70,72,78]}. More specifically, the study of dynamic fracture finds application in engineering ^[8,18,54], automotive safety ^[33], medicine ^[69], material science ^[28], and biology ^[74]. Examples of dynamic fracture events include rapid crack propagation in structural components under time-dependent loading due to, e.g., impacts, explosions, or earthquakes ^[24], fragmentation of brittle materials under high-velocity impacts, and crack propagation in geological formations due to seismic activity. The study of dynamic fracture mechanics aims to unravel the underlying mechanisms governing the initiation, propagation, and arrest of cracks under these extreme loading conditions.

This paper presents a physico-mathematical model that addresses dynamic fracture propagation in brittle materials by applying a purely hemi-variational continuum mechanics-based strain gradient theory. This paper is based on experimental analyses conducted in ^[41], which delved into the dependence of the fracture path on the loading rate. In a theoretical paper aimed at predicting mathematically the fracture path and the related branching phenomena, the same authors utilized a microplane model ^[42] in a 3D numerical simulation. In works ^[5,22,55,71], in order to describe the experimental results on dynamic fracture propagation reported in ^[41], different approaches were used like peridynamic, microplane, and phase-field damage modeling. The recent work of Nguyen et al. ^[38] gives a reformulation, in terms of sprain and spress, of a subset of micromorphic theories discussed more than half a century ago by Germain ^[21], Mindlin ^[34], and Eringen ^[16], whose micromechanical foundations for granular systems have been elaborated in a variational approach for finite deformation ^[36], geometrically nonlinear strain gradient models ^[4], and n-th order micromorphic materials ^[37]. Following this granular micromechanical approach and utilizing hemivariational methods, a group of researchers in a series of papers ^{[10,35,47,49,50,52,67]} have clearly established the physical meaning of boundary conditions associated with strain gradient theories ^{[2,3,11,12,13,14,15,62,63,66,76]}. Moreover, they have shown that for heterogeneous systems, strain gradient theories with appropriate elastic and dissipative energies can predict a multitude of phenomena, such as evolution of chirality and size of localization/fracture zone, that are oftentimes ignored or overlooked in classical approaches ^[9,19,29,31]. The strain gradient with the inertia gradient effects approach was used in ^[46,65], the treatment being purely analytical, without application of numerical methods such as the finite element method (FEM). In the papers of Rao et al. ^[56,57], strain gradient effects in the modeling of fracture mechanics have been considered with a two-scale asymptotic analysis. In this study, we employ an isotropic model to investigate the behavior of experiments presented in ^[41]. Our objective is to tune the evolution of damage based on discerning the local strain/stress state, namely, tension/compression. Specifically, we aim at mitigating, avoiding actually, its occurrence in compression. The utilization of an isotropic model accounts for uniformity of the microstructure with respect to directions.

The paper is organized as follows. In Section 2, we introduce the strain-gradient continuum model–elastic and damage components-employed to investigate numerically the propagation of dynamic fracture. In Section 3, we apply the presented model to an experimental case study that is then used for comparison with experimental data and other theoretical approaches available in the literature by solving through the weak form package of the software COMSOL Multiphysics and the weak form of its governing equations. Finally, conclusions are drawn in Section 4.

2.   Formulation of the problem

2.1.   Preliminary definitions

Let us consider a 2D continuum body $\mathfrak{V}\subset E^{2}$, where the domain $E^{2}$ is the Euclidean two-dimensional space. The position of every point in the body is established through the utilization of $X$ coordinates in a frame of reference. The set of kinematic descriptors for such a model contains both the displacement field $u = u(X, t)$ and the damage field $\omega = \omega(X, t)$. The damage state of a material point $X$ is therefore characterized, at time $t$, by a scalar internal variable $\omega$, which is assumed to take values within the range [0, 1]. The case $\omega = 0$ corresponds to the undamaged state, while $\omega = 1$ corresponds to complete failure. We assume the material to be not self-healing and, hence, $\omega$ is assumed to be a nondecreasing function of time. If we consider the material as a microstructured one, it becomes significant to allow the deformation energy density to vary based on the second gradient of the displacement. Thus, because of objectivity, the potential energy $U$ is assumed to be a function of the strain tensor $G$, of its gradient $\nabla G$, and of the damage $\omega$, i.e., $U = U(G, \nabla G, \omega)$. Here, $G = \frac{1}{2}\left(F^{T}F-I\right)$, where $F = \nabla\chi$ is the deformation gradient tensor, $\chi = X+u$ is the placement function, and $F^{T}$ is the transpose of $F$. For the small displacement approximation that is dealt with in this work, the strain tensor $G$ is assumed to be equal to $E = Sym(\nabla u)$, the symmetric part of the displacement gradient.

2.2.   The Energy functional

The energy functional $\mathcal{E}(u, \omega)$ depends on the displacement $u$ and the damage $\omega$ and is given as follows ^[53]:

where $\partial\mathfrak{B}$ are the regular parts of the boundary of $\mathfrak{B}$, the symbol $[\partial\partial\mathfrak{B}]$ denotes the regular parts of the boundary of $\partial\mathfrak{B}$, i.e., its wedges, and $n$ is the unit external normal. The symbols (·) and (:) indicate the scalar product between vectors or tensors, respectively. The quantities $b^{ext}$ and $m^{ext}$ are, respectively, the external body force and double force (per unit area), while the quantities $t^{ext}$ and $\tau^{ext}$ are, respectively, the external force and double force (per unit length). Finally, the quantity $f^{ext}$ is the external concentrated force which is applied on the vertices $[\partial\partial\mathfrak{B}]$.

2.3.   Variational principle

In order to get governing equations for this model, we resort to the variational principle thoroughly presented in ^{[32,48,49,51,53,75]}, where the case of one-dimensional bodies was considered. We assume that the descriptors $u(X, t)$ and $\omega(X, t)$ verify the condition

where $\nu$ and $\beta$ are compatible virtual velocities and dots represent derivation with respect to time.

2.4.   Reduction from the 3D Mindlin elastic strain gradient energy to the 2D explicit form

A general form of the deformation energy density of an elastic 3D isotropic strain gradient elastic material is provided in ^[34]. In our case, we have

where $\mathbb{H}$ is the Heaviside function, $TrG$ is the trace of $G$, $r$ affects on the 1st gradient residual inside the fracture, and $m$ affects on microstructure due to damage. Parameters $r$ and $m$ should be as close as possible to the value 1. $\gamma$ is a constitutive parameter representing the threshold of the trace of $G$ controlling the bi-modular behavior of damage, which stipulates that damage evolution occurs only when the trace of the strain tensor $G$ is positive and exceeds a predefined threshold. By implementing this condition, we ensure that damage progression is assessed solely in regions experiencing tensile stresses. The first gradient part $U^I_e$ of the deformation energy density above reads as

and the strain gradient part $U^{II}_e$ reads as

where $\lambda$, $\mu$, $c_{3}$, $c_{4}$, $c_{5}$, $c_{6}$, $c_{7}$ $\in\mathbb{R}$ are independent parameters. Our model accounts for an isotropic material, because the key focus is on the examination of damage evolution exclusively under tension while avoiding its manifestation in compression, aimed at capturing the behavior of cementitious materials. The quantities $\lambda$ and $\mu$ in Eq (2.4) are referred to as Lamé parameters and it is a classical result in mechanics that they can be easily related to the Young modulus $Y$ and to the Poisson's ratio $\nu$. For the 2D case, we have

The coefficients $c_{3}$–$c_{7}$ can instead be identified in terms of Young modulus and Lamé parameters, and an internal characteristic length–also referred to as the intergranular distance–within the granular micro-mechanics-based framework presented in ^[4] is as follows:

We now introduce the dissipated energy $U_d$ as

where $k_{\omega}$ represents the resistance to damage. The kinetic energy $U_{k}$ is instead assumed to having the following expression:

where $\rho$ is material density per unit area and $u_{1}$ and $u_{2}$ are the displacement components in the horizontal and vertical directions, respectively. It is easy to prove ^[1] that having a nonzero kinetic energy in the form above is equivalent to have $b^{ext} = -\rho\ddot{u}$.

3.   Numerical results

The weak form package of the commercial software COMSOL Multiphysics has been employed to perform the numerical simulation presented in this section, based on the modeling presented so far.

3.1.   Description of the numerical experiments

In the current section, we present numerical simulations to show the capabilities of the presented model to describe initiation and growth of damage localization zones. Following ^[41], we consider a 2D square specimen with a notch (pre-crack). The followig type of numerical experiment is considered: the left edge of the notch (the blue line) is fixed, i.e., cannot move, along the direction $\hat{e_{x}}$, while the right edge of the notch (the red line) is subjected to an applied velocity in the $\hat{e_{x}}$ direction (from left to right). Such an applied velocity will be also referred to as crack-opening velocity.

Figure 1 illustrates this loading condition. In the original paper ^[41], there are twelve cases of different applied velocities, the loading-rates being in the range from 0.045 m/s to 4.298 m/s. In this section, only four cases of tested velocities are considered for the study, which are Case 1: 0.491 m/s; Case 2: 1.375 m/s; Case 3: 3.268 m/s, and Case 4: 4.298 m/s. In the current study, the applied velocity is assumed to depend upon time–the data of the actual loading history used in the experiments was not reported in the original paper, as follows:

with $t_{0}$ = 100 $\mu s$ and $t$ = 500 $\mu s$. The prescribed displacement velocity as a function of time for $v_{0}$ = 3.268 m/s is reported in Figure 2.

The constitutive parameters used in the numerical simulations are reported in Table 1, where $t_{f}$ is the time horizon of the simulations and $\triangle t$ the time step. It should be noted that since simulations were performed for the 2D case, to take into account the (uniform) thickness of the specimen–25 mm–some constitutive parameters were adapted by multiplying them by this value.

3.2.   Convergence analysis

For convergence analysis, we employed two kinds of numerical simulations: elastic and with damage. In the former kind of simulations, to avoid damage and have only an elastic behavior, the threshold on the trace of $G$ was set to 1, an extremely large value, that was never reached in the simulations. In the elastic simulations, it was shown that strain gradient regularization gives better convergence with respect to maximum finite element size (the numbers of degrees of freedom (DOF), and finite elements grow with respect to maximum finite element size). Meanwhile without it, the mesh should be more refined to get the same accuracy.

From Table 2 and Figure 5, we can see that, without implementation of the strain gradient regularization, there is no convergence for the considered values of the maximum finite element size, however the implementation of this term leads us to convergence with refinement of mesh. From Table 3 and Figure 8, we can see that enhanced convergence was achieved making use of the strain gradient regularization, convergence is indeed obtained for maximum finite element size equal to 0.002 m. The same approach to mesh and convergence analysis utilized for the purely elastic case was employed also in the damage case, with applied displacement velocity $v_{0}$ = 3.268 m/s.

Figure 9 shows the different utilized meshes, each characterized by a different maximum element size. Figure 10 shows synoptically reaction force versus time diagrams (in the first row) and damage contour plots at complete failure (i.e., the fracture pattern, in the second row) for different maximum finite element sizes. It is seen that the fracture pattern is similar in all cases. Figure 11 compares the fracture thickness for different meshes. For the maximum finite element size equal to 0.004 m, we can see only one damaged element in the crack-transverse direction. However, as the element size decreases, the number of damaged elements in the crack-transverse direction increases. Most importantly, the thickness of the crack is almost the same for all simulations, which leads us to conclude the mesh independence of the results. From Figure 12 we can see that the evolution of the reaction force in time is almost the same for the different considered maximum element sizes, except that it attains its maximum value at different time instants. Additionally, fracture reaches the specimen's boundary at different times. A maximum mesh size of 0.001 m. was chosen for the simulations that follow in the next subsection.

3.3.   Numerical results with varying loading rate.

According to the convergence analysis presented above, a set of numerical simulations was performed making use of the same set of constitutive parameters in Table 1 and a mesh with maximum element size equal to 0.001 m. The difference between the different performed simulations lies in the value of the applied displacement velocity. The obtained results were compared with the experimental data.

Figures 13–15 show, respectively, a comparison of numerical simulation results with experimental ones in terms of reaction force versus imposed displacement diagram and fracture pattern, contour plots of damage evolution, and contour plots of kinetic energy density evolution for applied crack-opening velocity $v_{0}$ = 0.491 m/s.

Figures 16–18 are analogous but refer to an applied crack-opening velocity $v_{0}$ = 1.375 m/s, while Figures 19–21 and Figures 22–24 do refer to an applied crack-opening velocity $v_{0}$ = 3.268 m/s and $v_{0}$ = 4.298 m/s, respectively. It is possible to observe that the path and character of the fracture match the experimental data, except for the higher loading rate. The reaction force versus imposed displacement diagram properties match relatively well, too. Furthermore, we observe that, for low loading rate, the propagation is almost vertical, leading to standard mode-Ⅰ crack propagation, while for increasing crack-opening velocity the fracture deviates from the vertical direction and branching phenomena are observed. For the applied crack-opening velocity $v_{0}$ = 1.375 m/s, the computed kink angle is the same as in the experimental data. However, at the end of the fracture process, the computed fracture touches a different point of the boundary. For the applied crack-opening velocity $v_{0}$ = 3.375 m/s, one can see simple branching, while for the applied crack-opening velocity $v_{0}$ = 4.298 m/s, we can see multiple branching. With this difference, we can explain that it does not correspond that all twelve experimental specimens were done from the same concrete and have exactly the same properties. As well, it is not enough that the statistics and some errors could occur during the experiments. Anyway numerical simulations are the way to predict the behavior and it could be improved. Table 4 reports synoptically, for comparison, the reaction force measured experimentally and computed with various approaches–including the one presented in this paper–for different applied velocities.

4.   Conclusions

In this paper, we have focused on the loading rate effect on fracture propagation in brittle materials with strain gradient effects. We have demonstrated that, the slower the applied loading velocity, the less the maximum reaction force; and the faster the applied loading velocity, the higher must be the imposed displacement to reach the maximum value of the reaction force; cfr. Figure 25. It is also observed that the nonsymmetric increasing loading rate leads to a deviation from standard mode-Ⅰ crack propagation consisting of crack deviation and simple and multiple branching that cannot be revealed in the quasi-static case. It is worth it to note that, as expected, kinetic energy concentrates in the specimen's region, directly putting in movement by the loading and by the occurrence of fracture, which allows crack opening and, hence, movement of the body particles in that region. Increasing the applied loading velocity, one can also see the growth of some orders of magnitude of the kinetic energy, which suggests that this is the key quantity underlying the loading rate dependent behavior of fracture propagation.

Author contributions

L.P, A.M. and E.B. conceived the presented idea and verified the analytical methods. V.M., L.P, A.M., D.T. and E.B developed the theory. V.M., L.P, F.J.L.T., A.M., D.T., F.F. and E.B discussed the theoretical model. V.M., D.T. and E.B performed the computations. V.M., L.P, F.J.L.T., C.D.S., A.M., D.T. and E.B discussed the numerical results. All authors contributed to the final manuscript.

Use of AI tools declaration

The Authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Emilio Barchiesi is an editorial boardmember for [Networks and Heterogeneous Media] and was not involved inthe editorial review or the decision to publish this article.

The research reported in the present contribution was carried out as part of the project "A Fluid-Structure Interaction tool for the protection of Clean Energy Production sites (FSI-CEP)" funded by the MUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2022 PNRR-grant P20227CSJ5.

Conflict of interest

The Authors declare that there is no conflict of interest.