Numerical techniques to simulate crack propagation can roughly be divided into sharp and diffuse interface methods. Two prominent approaches to quantitative dynamic fracture analysis are compared here. Specifically, an adaptive cohesive element technique and a phase-field fracture approach are applied to simulate Hopkinson bar experiments on the fracture toughness of high-performance concrete. The experimental results are validated numerically in the sense of an inverse analysis. Both methods allow predictive numerical simulations of crack growth with an a priori unknown path and determine the related material parameter in a quantitative manner. Reliability, precision, and numerical costs differ however.
- Split-Hopkinson bar experiment
- cohesive elements
- phase-field fracture
- inverse analysis
- dynamic fracture
- crack propagation
- crack tracking algorithms
One of the main challenges in computational mechanics is the prediction of cracks and fragmentation in dynamic fracture. There are high demands on the modeling side, but mainly the complicated structure and the nonregular behavior of the cracks turn numerical simulations into a difficult task. Every crack in a solid forms a new surface of a priori unknown position, which needs to be identified. Different discretization techniques have been developed to solve such problems, for example the cohesive element technique [1, 2, 3], the extended finite element method [4, 5], eroded finite elements or eigenfracture strategies [6, 7], and phase-field approaches [8, 9, 10, 11, 12, 13].
The numerical techniques to treat the moving boundary problem of crack propagation can roughly be divided into two different strategies: sharp interface and diffuse interface modeling. The sharp interface approach describes a crack as a new boundary in a solid of domain undergoing a deformation in a time . For a known crack path this is the natural way to capture the mechanics of fracture. In dynamic fracture, however, with a priori unknown ways of crack propagation, kinks, and branching, sophisticated tracking methods need to be employed to localize the boundaries by the position and to enforce (free) boundary conditions along the moving crack. Unfortunately, such surface tracking methods require a high numerical effort and tend to fail for great changes in the topology of the solid, such as the fragmentation into small particles and their further movement.
An alternative way to describe moving boundaries are diffuse interface models where the cracks are smeared over a small but finite length
. Here an additional field
characterizes the state of the material and marks the intact or broken state. The set of evolving crack surfaces is replaced by a crack-surface density
, which is typically a function of the marker field
being not only a function of field
but also its gradient
, it regularizes (or
Here we compare a sharp interface method with crack tracking algorithm and a diffuse interface method for its usability in material identification. Background for our comparison are our experimental investigations on the fracture toughness of ultra-high performance concrete (UHPC). Specifically, we use the cohesive element technique and the phase-field fracture approach to simulate spalling experiments performed with concrete specimen in a Hopkinson-Bar (HB) setup.
UHPC is a class of advanced cementitious-based composites whose mechanical strength and durability surpass classical concrete. Typically, UHPC composites are fine grained, almost homogeneous mixtures of small aggregates of cement, a certain amount of silica, other supplements, and a low water content—and so they are more similar to brittle ceramics than to construction concrete. UHPCs are still under development and in order to optimize their composition mechanical tests have to provide material data. Hereby classical experiments determine the concrete’s elasticity as well as its compressive and flexural strength under static loading conditions. For the dynamic properties, however, such as dynamic tensile resistance and fracture energy, it is more complicated to ensure reproducible test conditions. Here numerical simulations in the sense of an inverse analysis are helpful to evaluate the reliability of the obtained material data.
HB spalling experiments are test arrangements to determine the failure strength of brittle materials, see [15, 16, 17, 18, 19]. In these tests the experimental setup of a classical HB is modified in such a way, that the induced pressure impulse is transmitted via an incident bar into the specimen, see Figure 1. Within the specimen a superposition of transmitted and reflected waves determines the stress state. For details of the experimental work we refer to another work , here we just use the experimental setup to compare two numerical techniques employed for
The remaining paper is organized as follows. In the next section we provide shortly the governing equations of elasto-dynamics and fracture mechanics. Then we introduce the cohesive element technique in Section 3 and the phase-field fracture method in Section 4. Both sections conclude with a short study on the influence of the relevant model parameters. In Section 5 the simulations of the HB spalling experiment are described in detail and a range of values for the fracture parameters is derived. The inverse analysis is presented in Section 6. Here we provide several numerical simulations and evaluate both methods. Such a quantitative comparison is new and has not yet been presented before. In particular, predictive applications of the phase-field approach to fracture are not common by now. A summary of the pros and cons of both methods in Section 7 concludes the paper.
2. Governing equations
We consider a body of domain with external boundary . The body’s displacement field at point and time is denoted by ; its velocity and acceleration fields are and . A crack splits the body into subbodies and induces a displacement jump on as . The displacements satisfy the Dirichlet boundary conditions at . The body is loaded with traction at boundary ; it holds .
Linear-elastic material is presumed to follow Hooke’s law with elastic strain energy density,
where the Lamé material parameters and are formulated with Young’s modulus and Poisson number . Equivalently we use the Hookean material tensor with components . The strain tensor describes small deformations ; the mechanical stress tensor follows as . It holds the balance of linear momentum,
where is the mass density and a prescribed body force density. Boundary conditions are prescribed as
with normal vector and traction ; denotes a jump. Eq. (4) 3 corresponds to traction-free crack boundaries. Additional initial conditions may apply.
2.2 Fracture mechanics
Let the evolving internal cracks be represented by a set of boundaries . According to the linear-elastic fracture theory of Griffith and Irwin [21, 22], a material fails upon attainment of a critical surface-energy density. The crack growth corresponds to the creation of new surfaces and hence the internal work of the body is composed of
where is commonly known as Griffith’s critical energy release rate (Griffith energy). An optimum of Eq. (5) corresponds to crack growth.
The specific energy
corresponds to the energy
Another fracture criterion is the crack tip opening displacement with critical value . The underlying theory of Dugdale  and Barenblatt  explains crack growth as loss of cohesion in a cohesive zone. If the decohesion is modeled in a nonlinear way, that is, for a given relation between the vector of cohesive traction and crack opening , a standard application of the -integral will establish a link between the critical energy release rate and the critical crack opening displacement. Starting with the relation and choosing a contour for the evaluation of the -integral that surrounds the cohesive zone gives for an effective crack opening ,
2.3 Weak form of the problem and finite element discretization
The motions of a solid can be characterized by recourse to Hamilton’s principle of stationary action. The action of a motion within a closed time interval defines a functional with the Lagrangian function as the difference of kinetic and potential energy, . Stationarity demands the first variation to vanish,
for all admissible test functions with on ; denotes the Hilbert space of weak square-integrable functions and its first derivatives. Eq. (8) must hold for all times , which leads to the weak momentum balance. For a body with stationary crack it is equivalent to
For discretization the domain is subdivided into a finite set of nonoverlapping elements. We make use of a conforming ansatz and approximate the displacement field and its variation with
where are the piecewise ansatz functions collected in a matrix and the vectors , contain all unknown nodal displacements of the nodes and its nodal variations , respectively. Plugging Eq. (10) into Eq. (8) and a straightforward calculation gives the global system of equations
where and denote the mass and stiffness matrix of the finite element discretization, and is the vector of external forces. The Hookean tensor is reformulated in a matrix , and we specify
Discretization in time is performed by an implicit Euler method, that is, with time step and, thus, velocity and acceleration are approximated by
3. The cohesive element technique
The nucleation and the propagation of cracks are efficiently modeled through the cohesive zone model where fracture is assumed to happen along an extended crack tip triggered by tractions on the crack flanks, [23, 24]. A particularly appealing aspect of the cohesive zone model is that it fits naturally in the framework of finite element analysis and leads directly to the cohesive element technique introduced by Needleman, Ortiz, and co-workers [25, 26, 27]. The main idea of this approach is to add cohesive interfaces between the continuum elements that are able to model crack growth, see Figure 2. We employ this classical cohesive element approach combined with an automatic fragmentation and cohesive surface insertion procedure. The method has proven to be reliable and efficient for numerous applications, see among others [2, 28, 29, 30].
In cohesive theories, the displacement jump across a cohesive surface
plays the role of a deformation measure while the tractions furnish the work-conjugate stress measure. They follow a traction separation law—the cohesive law—which models locally the loss of material resistance during cracking. If the cohesive element has attained a critical opening displacement , no tractions can be transfered and the adjacent continuum elements are de-facto disconnected.
Typically, isotropic and anisotropic materials behave differently in crack opening (mode-I separation) and sliding (mode-II and mode-III separation) and, therefore, normal and tangential components of the displacement jump across the surface have to be treated differently. Given a vector field over , its normal and tangential components are and , and, the corresponding jump components and follow from Eq. (13). To further simplify the formulation of mixed-mode cohesive laws, we follow  and introduce an effective opening displacement
Here the parameter assigns different weights to the sliding and normal opening components. This allows us to formulate the traction as a function of the effective opening displacement only.
3.2 Cohesive laws
A cohesive law defines the relation between crack opening displacements d and tractions on the crack flanks . Generally, a cohesive law can be stated in the general form , where is the specific fracture energy describing the dissipation in the cohesive zone. It is subject to the restrictions imposed by material frame indifference, material symmetry, and the isotropy of sliding. The most general dependence of has the form and with Eq. (14) follows . Thus, a key benefit of the potential structure of the cohesive law is that it reduces the identification of the cohesive law from the three components of to a single scalar function. Then an effective traction follows,
An appropriate choice of cohesive variable is the maximum attained (effective) crack opening displacement . Loading of the cohesive surface is then characterized by the conditions and ; all other states correspond to unloading. Figure 3 shows different common loading envelopes, whereas unloading is commonly assumed to be linear in to the origin,
The simplest cohesive law for brittle materials has a linear loading envelope
The two parameter cohesive strength and critical opening displacement determine via Eq. (6) the specific fracture energy
There are several modifications of Eq. (17), for example bilinear laws for concrete  or convex cohesive laws for ductile materials . A cohesive law that can be adapted to brittle and ductile behavior is the universal binding law of Smith and Ferrante .
where is the position of the maximal traction. Note that here the cohesive stress does not vanish at the critical separation . The corresponding fracture energy is . Unloading at follows a linear relation of the form Eq. (16). Upon closure, the cohesive surfaces are subject to unilateral contact constraints, including friction. Therefore, suitable contact conditions have to be applied under compression. Friction in the cohesive zone is regarded as an independent phenomena, which is not modeled here.
3.3 Finite element implementation
The class of cohesive elements considered here consists of two surface elements, which coincide in the reference configuration of the solid. Each surface element has nodes; the total number of nodes of the cohesive element is . The particular triangulation depicted in Figure 2 is compatible with edge or line elements of four nodes.
Basis of the finite element implementation is Hamilton’s principle given in Eq. (7). Inserting the balance of linear momentum and the static boundary conditions gives the deformation power with variation
where are the total cohesive surfaces and is the variation of the separation vector given in Eq. (13). The first term of Eq. (20) 1 as well as the remaining energy contributions correspond to standard finite element forms and will not be repeated here. The variation of the cohesive energy leads with ansatz Eq. (10) in Eq. (20) 2 for one cohesive element to
The kinetic energy does not have any support in the cohesive element and only the external virtual work has to be determined. For one cohesive element it is
The tangent stiffness matrix follows by its consistent linearization, with the result . To avoid singularities at in the derivatives of Eq. (17) some numerical modifications are required, cf. . Finally, the equations added by the cohesive elements to the system Eq. (11) can be formulated as with , where symbolizes the incremental solution procedure, which is required for the nonlinear crack-opening problem.
3.4 Adaptive meshing
Since in most problems the expected crack path is not known the decision where the cohesive elements should be inserted has to be made during the simulation. The analysis proceeds incrementally in time. Our decision criterion is based on the effective tensile stress given in Eq. (15), which has to exceed a threshold. This means, in every time step of the calculations, this condition is checked for each internal face. The faces that met the criterion are flagged for subsequent processing. A cohesive element will be inserted at the flagged face and in this manner, the shape and location of a successive crack front is itself an outcome of the calculations.
Within the finite element mesh the insertion of cohesive elements requires topological changes. The local sequential numbering of the corner-nodes defines the orientation; the mid-side node is subsequently duplicated. Owing to the variable environment of the edges in the triangulation, the data structure has to be adapted as illustrated in Figure 4 with case 1: the marked edge with nodes
3.5 Effect of the crack initiation criteria in the cohesive model
Here we illustrate the influence of the cohesive strength and critical opening displacement on the cohesive element simulations. Exemplarily we investigate a plane mode-I tension test under quasistatic conditions; however, the observations have been confirmed also in other models.
The test specimen with material data has a unit domain ; all lengths are given in meter. On the vertical boundaries the displacements in -direction are constrained and on the lower and upper side and , a traction boundary condition is prescribed, . A crack is predefined in the area , , by including cohesive elements that are completely open. We employ a cohesive law with linear envelope given in Eq. (17), set the cohesive strength and the critical crack opening displacement . Following Eq. (18) this corresponds to a critical energy release rate of .
During the simulation a cohesive element will be added if the effective traction Eq. (15), here , exceeds a critical value, . We identify the insertion criterium with the cohesive strength . Because the specimen is pulled with constant load, the computation is stopped at a crack length of 0.7 m.
Figure 5 demonstrates the crack evolution in a mesh of 25 × 25 squares, each divided into eight triangular finite elements with linear shape functions. The computed stresses in Figure 5 correspond very well to the analytical solution. They disappear at the crack flanks and are maximal at the crack tip with the typical butterfly shape.
At next we apply the traction MPa as a linearly increasing load within 200 load steps. A cohesive element will be added if the effective tensile stress exceeds . The crack length is determined via . Figure 6 shows the crack growth normalized to the crack growth at MPa. Obviously, the insertion criteria has an influence on the appearance of cohesive elements and so—indirectly—on the crack propagation. In consequence, the insertion stress should correspond to the cohesive strength to avoid artificial numerical stiffness.
4. The phase-field fracture approach
The evolving crack in a solid with potential energy of Eq. (5) is represented in the phase-field fracture approach by an additional continuous field
. It has a value of
in the intact material and indicates for
The parameter has the unit of a length and is a measure for the width of the diffuse interface zone, see Figure 7. Moreover, this parameter weights the influence of the linear and the gradient term whereby the gradient enforces a regularization of the sharp interface. Another ansatz for the crack-surface density function is the fourth-order form
which has been used, for example, in [38, 39]. However, in a finite element discretization ansatz given in Eq. (24) requires -continuous basis functions. In  we investigated this formulation in more detail, applied quadratic NURBS-ansatz functions, and found that a smoother crack and a better convergence rate can be achieved.
Here we use the ansatz of Eq. (23) and the corresponding total potential energy reads
in the elastic strain energy density is derived from Eq. (2) by means of a substitute-material approach,
is a degradation function, which is here of the form
being a small parameter introduced for numerical reasons only. With
At this point the evolution equation for the phase-field parameter
where the nonnegative function
is the mobility with unit [
denotes the variational or total derivative of
4.1 Elastic strain energy split
The quadratic form of the elastic strain energy density does not distinguish between tensile and pressure states in the material. A direct use of the formulations Eq. (25) or Eq. (27) would allow a crack to grow also in a compressive regime, which clearly contradicts the physics of the underlying problem. For that reason a split of the elastic strain energy into a tensile and a compressive part is necessary. We use the degradation function and state . This split is based on the spectral decomposition of the strain tensor , where denote the principal strains and the corresponding principal directions, . Based on this representations and using the Macaulay brackets
we define the positive and the negative parts of the strain tensor as . The positive parts contain contributions due to positive dilatation and contributions due to positive principal strains. Only this part of the strain energy is responsible for crack growth. A similar decomposition can be deduced for the stress tensor . This finally leads to an elastic energy density function, which only accounts for tension, .
Furthermore, the irreversibility of the crack growth has to be considered. This can be done by Dirichlet constraints on the phase-field parameter, that is if . The same effect has a product of with the mobility parameter , where we additionally formulate the evolution Eq. (26) in a dimensionless form, .
The numerical solution of the phase-field fracture model within the finite element framework leads to a coupled-field problem. The weak formulation of the mechanical field is derived in the usual way with the result given in Eq. (9). The weak formulation of the phase-field equation is set up analogously
with test functions and .
with ansatz functions for nodes. Plugging this ansatz into Eq. (28) leads after a straightforward calculation to a finite element system of equations. Additionally we specify the phase-field driving force .
Here we assumed the same ansatz functions for
After the discretization in time by using Eq. (12) the following system of global equations is obtained:
Note that the coupled problem of Eq. (11) and Eq. (30) is nonlinear. The solution of the implicit problem is obtained with recourse to a Newton-Raphson method. The necessary linearization (tangent stiffness matrix) can be calculated monolithically or by recourse to a staggered scheme. For the HB experiment we employ an explicit time discretization, which simplifies the solution.
4.3 Effect of the coefficients
With the problem formulation at hand we now illustrate the influence of the parameter mobility
The kinematic mobility parameter
is responsible for the rate of phase-field evolution. For higher values of
The influence of the length-scale parameter
on the crack is now demonstrated. For definiteness, the domain where the phase-field parameter
whereby the Lebesgue-measure of a set is denoted with and the diffuse interface zone is included in the relevant set.
Please note, that the effect of the length-scale parameter is two-fold. On one hand, it enters the material because, in the sense of Griffith’s criteria for crack growth, the fracture energy density competes with the elastic energy density determined by or combinations thereof. Here has the effect of a material parameter. On the other hand, is determined by the mesh size because it has to be large enough to enable the approximation of a diffuse interface. This clearly requires an adaption on the mesh size, .
In Figure 9 the crack width
5. Simulation of the HB-spalling experiment
A classical Split-Hopkinson-Pressure Bar consists of a steel projectile (striker), an incident bar and a transmission bar. The specimen is placed between the bars and an analysis of the propagating waves allows to deduce its Young’s modulus. For our UHPC mixture the result is
GPa; details of the experiments can be found in [20, 44]. For the
The aim of spalling experiments is to determine a material’s resistance to fracture, specifically its fracture energy or specific energy and/or its tensile resistance . The spalling test presumes brittle materials which can sustain compression but fails under tension. Plastic deformations do not matter, neither do temperature effects; all tests are conducted under ambient conditions.
5.1 Finite element discretization
The UHPC specimen has a length of 200 mm and a diameter of 20 mm. Because of the cylindrical symmetry of the problem we can use an axialsymmetric finite element model. This model maps a fully three-dimensional material behavior with the reduced effort of a plane mesh, which allows us to do extensive parametric studies.
A first challenge was the correct reproduction of the incident and reflected stress pulses in the specimen. From the strain pulse measured in the incident bar, the difference in impedance, and a low-amplitude pulse measured in the specimen we conclude on the shape of the transmitted wave. It is applied on the (left) boundary as a pressure impulse of trapezoidal form with , , and, , where and and and . Within the of simulation time the stress wave travels through the specimen, is reflected and reaches the left end again. Further propagation (and further cracks) are not considered. The average velocity of the specimen before spallation is m/s. For time discretization we use a special central difference scheme with a weighted displacement field, which results in stress pulses that largely correspond to the measured data, see Figure 10; for numerical details we refer to [34, 45].
5.2 Fracture parameter
The dynamic tensile resistance of a brittle material
is usually defined as the maximum tension a material can sustain. A higher stress results in fracture—in our experiments in spallation. There are several attempts but no established way to deduce the tensile resistance directly from the data measured in the experiments, see  for a discussion. One way is to measure the incident and reflected waves, “shift” them to the position of fracture
5.2.1 Cohesive element technique
Our specimen is meshed uniformly with triangular finite elements (2560 elements) and a mixed-mode cohesive law is employed. We start with the linear envelope given in Eq. (17) and an effective opening displacement m for , N/m, . As outlined in Section 3.5 a cohesive element will be added if the effective tensile stress given in Eq. (15) exceeds the value of . In tension, the elements can subsequently open, in compression contact conditions apply. Once the opening has exceeded , a crack has formed. The simulation stops after a spallation plane has built or after .
Figure 11 shows one symmetry half of the specimen at the end of the simulation for different values of . At first, with , the cohesive stress is obviously too high, no elements are inserted and no crack can grow. For a small, localized crack zone develops. Lowering MPa gives a wider zone and, moreover, the cracked zone moves toward the free end. This follows from the fact that the superposed pulse reaches the value of earlier. However, this position does not correspond to the measured crack position. In all our experiments the crack appeared at mm, with coordinate starting at the free end.
Obviously, the cohesive stress plays a significant role in the simulation of crack growth with cohesive elements. If it is too low, an unrealistic scattered crack zone will be computed, which does not give fully opened cracks (assuming a fixed fracture energy). However, some energy is dissipated here and the crack plane cannot appear at the right position anymore.
We further studied the influence of the critical crack opening displacement . Here we observe for higher values of a longer time of crack opening whereas has no influence on the position of the crack. Studies with the exponential cohesive law given in Eq. (19) show a very similar behavior in variations of and , cf. Figure 12.
5.2.2 Phase-field fracture
In phase-field simulations the fracture energy is the essential parameter for crack growth. Other parameter, like the mobility, are of numerical nature and can be calibrated. Further, the length-scale parameter is responsible for the width of the diffuse interface. Because of the fact, that depends on the mesh size, the mesh has to be very fine in the domain of a potential crack. To avoid numerical artifacts we chose here a uniform mesh of 5760 triangular elements; the mesh size is mm. We set and and define regions with phase-field parameter as cracked.
For a large critical energy release rate, , no crack is computed. Obviously, here such values are too high. Reducing the value leads to the appearance of a zone with ; however, it depends on the defining threshold if this is already considered to be a crack. For a sharp crack zone appears at the expected position. Smaller values, for example in Figure 13(left), show a wider cracked zone but the position of the crack is nearly the same, see Figure 13(right). Specifically we define the diffuse zone as
and the cracked zone is given by .
Figure 14 shows the effect of the specific fracture energy on the time of crack evolution. In opposite to the cohesive model the time difference between crack formation and final state grows linearly with . Similar is the situation for an decrease of parameter . A value of increases the time of crack formation at constant .
Further studies have been performed to simulate the experiments, cf. , and with the collected knowledge we conclude, that for UHPC the tensile strength , the critical opening displacement , and the specific fracture energy are in the ranges and .
6. Inverse analysis: determination of the
The aim of our study was to evaluate the two fracture simulation methods for the use in an inverse analysis where the measured data of the experiment are used to deduce fracture parameters and the obtained results are used to simulate the experiment. The deduced parameters are considered “correct” when the difference between experiment and simulation is small.
6.1 Derivation of the specific fracture energy from spallation
After spallation two fragments result with the crack located at the position where the stress exceeds the tensile resistance first. Depending on the energy of the incoming wave the same process may continue in both fragments with the results of additional cracks. The total fracture energy
corresponds to the amount of work necessary to form such a new surface. In order to determine
we balance the energy before and after crack initiation, that is, at time
where refer to the velocities of the fragments 1 and 2. Because we presume the UHPC to behave linear elastically, the loss of kinetic energy between initial and spalled state will completely be related to the fracturing process. Consequently, we state to be the fracture energy of the specimen.
In order to deduce the fracture energy Eq. (33) needs to be referred to the fractured surface of the specimen. Ideally, the crack is smooth and perpendicular to the specimen’s axis. Then the specific fracture energy simply follows as with . In practice, there may be spalled and rough crack surfaces and so we replaced by the measured surface size.
For the inverse analysis we proceed as follows: we use the data obtained in our previous simulations to define a range of input data , simulate the spallation experiment, and from the computed velocities and masses we deduce via Eq. (33) the “measured” value . To this end we determine the velocities of the fragments by , , and, , with the kinetic energy of the finite elements, .
6.2 Cohesive element technique
We mesh the specimen uniformly with 2560 elements and employ the linear cohesive law from Eq. (17) with parameter . The cohesive stress is MPa to ensure crack formation at the expected position and we vary the critical opening displacement to obtain values of between 20 N/m and 150 N/m. All computations result in similar fragments with masses and . Also the fractured surfaces show little variations, cm2.
In Figure 16 the computed specific fracture energy is shown. In the optimal case is , this is marked with a dotted line. The computed values show a clear proportionality between and its simulated counterpart . This generally validates the chosen method of energy balance. However, the value of deduced from the simulation is higher than the real one. This corresponds to a higher difference in the kinetic energy, that is, for the given initial velocity the velocity of the fragments is underestimated, that is there is too much energy dissipated. This likely results from the fact that some of the adaptively inserted cohesive elements are “useless,” they are inserted in a wrong place and do not fully open. This partial opening costs energy, which is dissipated but does not contribute to spallation. Therefore we observe with discrepancies between 30 and 60%, see Figure 16.
6.3 Phase-field fracture
For phase-field simulations we use a finer mesh of 5760 elements, and set and . At first we accelerate the specimen to , then the actual simulation starts. Again, during spallation the velocity of the fragments develops differently, the first fragment moves faster then the second one. The time of fracture depends on the input value , the lower it is the easier both fragments split, cf. Figure 14.
Since the phase-field model is a diffuse interface approach, which does give a discrete distinction between the both states, it is necessary to specify tolerances for the states and . Here we consider and set for symmetry . A lower “crack initiation” point clearly needs less time and energy to open the crack. On the one hand, if crack initiation equals crack opening, , that is the crack appears immediately after initiation and the dissipated energy will be close to zero. The input value would work only as a threshold for the fragments to loose contact and (almost) no dissipation will take place. On the other hand, if we raise , the crack may not fully open, which gives inaccurate fragment velocities. A stable state we found for and so we set and .
In Figure 17 the computed specific fracture energy over the input value is shown. Again, the optimal case is marked with a dotted line and the computed values show a clear proportionality between and . The value of is now lower than the input value with relative difference less than 10%. Reason for this is basically the definition of crack initiation with , which attributes the beginning of the crack still to the full specimen. Please note that for very small values of N/m the specimen cracks basically without dissipation whereas too high values of prevent any crack. In the range of interest, N/m both values and correlate very well with an R-coefficient of 0.9969. Specifically we have . Therefore, it is possible and useful to employ a phase-field fracture simulation for a quantitative analysis of the experiments.
In the previous we compared the possibilities of a sharp interface method and a diffuse interface method for crack nucleation and quantitative dynamic fracture analysis. Exemplarily, we validated investigations on the fracture toughness of high-performance concrete in a Hopkinson bar spallation experiment whereby, in particular, the fracture energy values have been determined. Both methods, the cohesive element technique and the phase-field fracture approach, allow numerical simulations of crack growth with an a priori unknown path, and both methods allow to determine the related material parameter in a quantitative manner. Reliability, precision, and numerical costs differ however. Pros and cons of both methods are summarized in the following.
7.1 Model parameters
The core of the cohesive zone model is a cohesive law, , which describes the forces between the crack flanks as a function of separation. Such cohesive laws allow for pure mode-I cracks in the sense of Griffith as well as for mixed-mode cracks, for example by using the effective traction and separation. Essential cohesive parameters are the critical cohesive stress , the critical separation , and the weight , which relates shear and tension. These parameters depend on the specific material and can be determined experimentally whereby is implicitly given via the specific crack energy . Further specifications of the cohesive law may require additional material parameters, for example, in the classical exponential Rose-Ferrante law an additional parameter needs to be set. All these parameters have a clear physical meaning.
Sensitive for the cohesive element technique is the critical traction for adaptive insertion of the cohesive elements, which has no direct physical background but strongly influences energy dissipation and numerical efficiency. Wrongly inserted elements may dissipate energy but do not contribute to fracture and skew the simulation results.
The phase-field approach to fracture is based on an evolution equation that essentially refers to the elastic strain energy density of the material. The remaining relevant parameters are the mobility , the specific crack energy , and the length-scale parameter . The mobility has only numerical character and controls the rate of phase-field decrease during crack formation. The critical length-scale parameter is a measure for the width of the diffuse interface, which relates to both, the finite element mesh size and the material properties. The latter enters the crack growth criteria via the term in the crack resistance and therefore needs to be calibrated carefully. The only material parameter in the phase-field model with a clear physical background is the Griffith energy .
7.2 Numerical implementation and computation
In the cohesive zone model cohesive elements are adaptively inserted between the continuum elements to describe the crack opening. The continuum elements themselves are not directly affected and the crack can only propagate along the element boundaries, which results in a certain mesh dependence. The adaptive insertion of cohesive elements require a continuous update of the data structure, which leads to a significant programming effort and also increases the costs of computation. Additionally, the cohesive zone has to be equipped with contact constraints in order to prevent penetration in case of unloading. In total, the numerical implementation of an adaptive cohesive zone model becomes very complex.
In contrast, the structure of the finite element mesh in the phase-field approach remains constant during the simulation. The phase-field parameter can decrease to zero at each node and the crack is able to propagate theoretically everywhere in the whole domain. Essential requirement is a very fine mesh with .
For numerical computation of the coupled fields
there exist two different types of solution. On the one hand it is possible to determine the displacements
and the phase-field parameter
7.3 Constraints and driving forces
The major advantage of the cohesive zone model is that the crack properties can be mapped exactly. The local opening is known, the crack width is the separation of the crack flanks, and the corresponding normal and tangential forces follow from the cohesive law. Since the loading and unloading processes are distinguished, compressive forces do not contribute to crack growth. The irreversibility of crack propagation is guaranteed and although the separation can decrease during the simulation the crack will not “heal.” Regarding mixed-mode problems it is positive that the normal and tangential traction components are weighted by a parameter so that different crack opening modes can be realized for each specific material under consideration.
In phase-field fracture by definition a continuous function
Summarizing we state that both methods are mechanically consistent and have a clear variational structure. The cohesive element technique is difficult to implement but provides a strong physical background. For static computations with expected way of crack propagation it is definitely preferred because it allows cohesive laws, which may consider anisotropy, friction, and other material specific properties. In general dynamic applications of unknown crack path, however, its numerical drawbacks, together with the fact that a suboptimal insertion may lead to wrong predictions, dominate. Here the phase-field approach to fracture is clearly the better choice. Unknown crack paths can simply be followed—as long as the mesh resolution is fine enough. The major drawback of phase-field fracture is its parameter sensitivity. Also, extensions to more complex fracture models, which account, for example, for sliding, anisotropy, and interlocking, contradict the original variational derivation and are still an open problem.
Pandolfi A, Ortiz M. Solid modeling aspects of three-dimensional fragmentation. Engineering with Computers. 1998; 14:287-308
Scheider I, Brocks W. Cohesive elements for thin-walled structures. Computational Materials Science. 2006; 37(1–2):101-109
Xu XP, Needleman A. Numerical simulations of dynamic interfacial crack growth allowing for crack growth away from the bond line. International Journal of Fracture. 1995; 74:253-275
Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics. 2002; 69(7):813-833
Stazi FL, Budyn E, Chessa J, Belytschko T. An extended finite element method with higher-order elements for curved cracks. Computational Mechanics. 2003; 31:38-48
Pandolfi A, Ortiz M. An eigenerosion approach to brittle fracture. International Journal for Numerical Methods in Engineering. 2012; 92:694-714
Schmidt B, Fraternali F, Ortiz M. Eigenfracture: An eigen deformation approach to variational fracture. SIAM Journal on Multiscale Modeling and Simulation. 2009; 7:1237-1266
Borden MJ, Verhoosel CV, Scott MA, Hughes TJR, Landis CM. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering. 2012; 217–220:77-95
Bourdin B. The variational formulation of brittle fracture: Numerical implementation and extensions. In: Volume 5 of IUTAM Symposium on Discretization Methods for Evolving Discontinuities, IUTAM Bookseries, Chapter 22; Springer, Netherlands. 2007. pp. 381-393
Henry H, Levine H. Dynamic instabilities of fracture under biaxial strain using a phase field model. Physics Review Letters. 2004; 93:105505
Karma A, Kessler DA, Levine H. Phase-field model of mode III dynamic fracture. Physical Review Letters. 2001; 81:045501
Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering. 2010; 199:2765-2778
Verhoosel CV, de Borst R. A phase-field model for cohesive fracture. International Journal for Numerical Methods in Engineering. 2013. In press
Ambrosio L, Tortorelli VM. Approximation of functional depending on jumps by elliptic functional via γ-convergence. Communications on Pure and Applied Mathematics. 1990; 348–349:13-16
Diaz-Rubio FG, Perez JR, Galvez VS. The spalling of long bars as a reliable method of measuring the dynamic tensile strength of ceramics. International Journal of Impact Engineering. 2002; 27:161-177
Klepaczko JR, Brara A. An experimental method for dynamic tensile testing of concrete by spalling. International Journal of Impact Engineering. 2001; 25(4):387-409
Schuler H, Mayrhofer C, Thoma K. Spall experiments for the measurement of the tensile strength and fracture energy of concrete at high strain rates. International Journal of Impact Engineering. 2006; 32(10):1635-1650
Weerheijma J, Van Doormaal JCAM. Tensile failure of concrete at high loading rates: New test data on strength and fracture energy from instrumented spalling tests. International Journal of Impact Engineering. 2007; 34:609
Zhang L, Hu S-S, Chen D-X, Yu Z-Q, Liu F. An experimental technique for spalling of concrete. Experimental Mechanics. 2009; 49(4):523-532
Khosravani MR, Wagner P, Fröhlich D, Weinberg K. Dynamic fracture investigations of ultra-high performance concrete by spalling tests. Engineering Structures. 2019; 201:109844
Griffith AA. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of London. 1921; 221:163-198
Irwin GR. Elasticity and plasticity: Fracture. In: Függe S, editor. Encyclopedia of Physics. 1958
Dugdale DS. Yielding of steel sheets containing clits. Journal of the Mechanics and Physics of Solids. 1960; 8:100-104
Barenblatt GI. The mathematical theory of equilibrium of cracks in brittle fracture. Advances in Applied Mechanics. 1962; 7:55-129
Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. International Journal of Solids and Structures. 1996; 33:2899-2938
Ortiz M, Pandolfi A. A class of cohesive elements for the simulation of three-dimensional crack propagation. International Journal for Numerical Methods in Engineering. 1999; 44:1267-1282
Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids. 1994; 42:1397-1434
Ferrara A, Pandolfi A. Numerical modelling of fracture of human arteries. Computer Methods in Biomechanics and Biomedical Engineering. 2008; 11:553-567
Pandolfi A, Weinberg K. A numerical approach to the analysis of failure modes in anisotropic plates. Engineering Fracture Mechanics. 2011; 78:2052-2069
Yu C, Pandolfi A, Ortiz M, Coker D, Rosakis AJ. Three-dimensional modeling of intersonic shear-crack growth in asymmetrically loaded unidirectional composite plates. International Journal of Solids and Structures. 2002; 39:6135-6157
Bazant ZP. Concrete fracture models: Testing and practice. Engineering Fracture Mechanics. 2002; 69:165-205
Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. Journal of the Mechanics and Physics of Solids. 1992; 40:1377-1397
Rose JH, Ferrante J, Smith JR. Universal binding energy curves for metals and bimetallic interfaces. Physical Review Letters. 1981; 47:675-678
Dally T. Vergleich von Kohäsivelement-Methode und Phasenfeld-Methode anhand ausgewählter Probleme der Bruchmechanik. Technical Report at Universität Siegen; 2017
Pandolfi A, Ortiz M. An efficient adaptive procedure for three-dimensional fragmentation simulations. Engineering with Computers. 2002; 18:48-159
Bourdin B, Francfort GA, Marigo JJ. The variational approach to fracture. Journal of Elasticity. 2008; 9:5-148
Francfort GA, Marigo J-J. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids. 1998; 46:1319-1342
Borden MJ, Hughes TJR, Landis CM, Verhoosel CV. A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering. 2014; 273(0):100-118
Hesch C, Stefan S, Dittmann M, Franke M, Weinberg K. Isogeometric analysis and hierarchical refinement for higher-order phase-field models. Computer Methods in Applied Mechanics and Engineering. 2016; 303:185-207
Weinberg K, Hesch C. A high-order finite-deformation phase-field approach to fracture. Continuum Mechanics and Thermodynamics. 2017; 29:935-945
Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering. 2010; 83:1273-1311
Bilgen C, Kopanicakova A, Krause R, Weinberg K. A phase-field approach to conchoidal fracture. Meccanica. 2018; 53:1203-1219
Weinberg K, Dally T, Schuß S, Werner M, Bilgen C. Modeling and numerical simulation of crack growth and damage with a phase field approach. GAMM-Mitteilungen. 2016; 39(1):55-77
Weinberg K, Khosravani MR. On the tensile resistance of UHPC at impact. The European Physical Journal Special Topics. 2018; 227:167-177
Park KC, Lim SJ, Huh H. A method for computation of discontinuous wave propagation in heterogeneous solids: Basic algorithm description and application to one-dimensional problems. International Journal for Numerical Methods in Engineering. 2012; 91:622-643
Erzar B, Forquin P. An experimental method to determine the tensile strength of concrete at high rates of strain. Experimental Mechanics. 2010; 50(7):941-955
Bilgen C, Weinberg K. On the crack-driving force of phase-field models in linearized and finite elasticity. Computer Methods in Applied Mechanics and Engineering. 2019; 353(15):348-372