Fracture as a material sink

Abstract

Cracks are created by massive breakage of molecular or atomic bonds. The latter, in its turn, leads to the highly localized loss of material, which is the reason why even closed cracks are visible by a naked eye. Thus, fracture can be interpreted as the local material sink. Mass conservation is violated locally in the area of material failure. We consider a theoretical formulation of the coupled mass and momenta balance equations for a description of fracture. Our focus is on brittle fracture and we propose a finite strain hyperelastic thermodynamic framework for the coupled mass-flow-elastic boundary value problem. The attractiveness of the proposed framework as compared to the traditional continuum damage theories is that no internal parameters (like damage variables, phase fields, etc.) are used while the regularization of the failure localization is provided by the physically sound law of mass balance.

Introduction

Within the framework of continuum mechanics there are surface and bulk material failure models.

Surface failure models are known by name of cohesive zone models (CZMs). In the latter case, continuum is enriched with discontinuities along surfaces—cohesive zones—with additional traction-displacement-separation constitutive laws. These laws are built qualitatively as follows: traction increases up to a maximum and then goes down to zero via increasing separation (Barenblatt 1959; Camacho and Ortiz 1996; de Borst 2001; Gong et al. 2012; Moes et al. 1999; Needleman 1987; Park et al. 2009; Rice and Wang 1989; Tvergaard and Hutchinson 1992; Xu and Needleman 1994). If the location of the separation surface is known in advance (e.g., fracture along weak interfaces) then the use of CZM is natural. Otherwise, the insertion of cracks in the bulk in the form of the separation surfaces remains an open problem, which includes definition of the criteria for crack nucleation, orientation, branching and arrest. Besides, the CZM approach presumes the simultaneous use of two different constitutive models: one for the cohesive zone and another for the bulk, for the same real material. Certainly, a correspondence between these two constitutive theories is desirable yet not promptly accessible. The issues concerning the CZM approach have been discussed by Needleman (2014), the pioneer of the field.

Bulk failure models are known by name of Continuum Damage Mechanics (CDM). In the latter case, material failure or damage is described by constitutive laws including softening in the form of the falling stress-strain curves (Benzerga et al. 2016; Dorfmann and Ogden 2004; Gurson 1977; Kachanov 1958; Klein and Gao 1998; Lemaitre and Desmorat 2005; Menzel and Steinmann 2001; Simo 1987; Volokh 2004; 2007; Voyiadjis and Kattan 1992). Remarkably, damage nucleation, propagation, branching and arrest naturally come out of the constitutive laws. Unfortunately, numerical simulations based on the the bulk failure laws show the so-called pathological mesh sensitivity, which means that the finer meshes lead to the narrower damage localization areas. In the limit case, the energy dissipation in failure tends to zero with the diminishing size of the computational mesh. This physically unacceptable mesh sensitivity is caused by the lack of a characteristic length in the traditional formulation of continuum mechanics. To surmount the latter pitfall gradient- or integral- type nonlocal continuum, formulations are used where a characteristic length is incorporated to limit the size of the spatial failure localization (Borden et al. 2012; de Borst and van der Giessen 1998; Francfort and Marigo 1998; Hofacker and Miehe 2012; Lasry and Belytschko 1988; Peerlings et al. 1996; Pijaudier-Cabot and Bazant 1987; Silling 2000). The regularization strategy rooted in the nonlocal continua formulations is attractive because it is lucid mathematically.

Unluckily, the generalized nonlocal continua theories are based (often tacitly) on the physical assumption of long-range particle interactions while the actual particle interactions are short-range—on nanometer or angstrom scale. Therefore, the physical basis for the nonlocal models appears disputable. A more physically based treatment of the pathological mesh sensitivity of the bulk failure simulations should likely include multi-physics coupling. Such an attempt to couple mass flow (sink) and finite elastic deformation within the framework of brittle fracture is considered in the present work.

Basic idea

Cracks are often thought of as material discontinuities of zero thickness. Such idealized point of view is probably applicable to nano-structures with perfect crystal organization. In the latter case fracture appears as a result of a separation—unzipping—of two adjacent atomic or molecular layers—Fig. 1 (left).

Fig. 1

Schematic drawing of cracks with zero or finite thickness

In the case of the bulk material with a sophisticated heterogeneous organization, the crack appears as a result of the development of multiple micro-cracks triggered by the massive breakage of molecular or atomic bonds—Fig. 1 (right). The bond breakage is not confined to two adjacent molecular layers, and the process involves thousands layers within an area or volume with the representative characteristic size l.

It is interesting to note that material failure does not require the breakage of all molecular or atomic bonds within a representative volume. Only fraction of these bonds should be broken for the material disintegration. For example, in the case of natural rubber, roughly speaking, every third bond should be broken within a representative volume to create a crack (Volokh 2013a).

The local bond failure leads to the highly localized loss of material. The latter, in our opinion, is the reason why even closed cracks are visible by a naked eye. Thus, material flows out of the system during the fracture process. The system becomes open from the thermodynamic standpoint. However, cracks usually have very small thickness and the amount of the lost material is negligible as compared to the whole bulk. The latter observation allows ignoring the additional supply of momenta and energy in the formulation of the initial boundary value problem described in the next sections.

Thermodynamics

Following the approach of continuum mechanics, we replace the discrete molecular structure of materials by a continuously distributed set of material points which undergo mappings from the initial (reference), Ω₀, to current, Ω, configuration: x↦y(x). The deformation in the vicinity of the material points is described by the deformation gradient F=Grady(x).

In what follows we use the Lagrangean description with respect to the initial or reference configuration and define the local mass balance in the form

$$ \frac{d\rho}{dt}=\text{Div}\mathbf{s}+\xi, $$

(1)

where ρ is the referential (Lagrangean) mass density; s is the referential mass flux; ξ is the referential mass source (sink); and Divs=∂ s _i /∂ x _i in Cartesian coordinates.

We further assume that failure and, consequently, mass flow are highly localized and the momenta and energy balance equations can be written in the standard form without adding momenta and energy due to the mass alterations.

In view of the assumption above, we write momenta and energy balance equations in the following forms accordingly

$$ \frac{d(\rho\mathbf{v})}{dt}=\text{Div}\mathbf{P}+\rho\mathbf{b},\quad\mathbf{P}\mathbf{F}^{\mathrm{T}}=\mathbf{F}\mathbf{P}^{\mathrm{T}}, $$

(2)

and

$$ \frac{d(\rho e)}{dt}=\mathbf{P}:\dot{\mathbf{F}}+\rho r-\text{Div}\mathbf{q}, $$

(3)

where $\mathbf {v}=\dot {\mathbf {y}}$ is the velocity of a material point; b is the body force per unit mass; P is the first Piola-Kirchhoff stress and (DivP) _i =∂ P _ij /∂ x _j ; e is the specific internal energy per unit mass; r is the specific heat source per unit mass; and q is the referential heat flux.

Entropy inequality reads

$$ \frac{d(\rho\eta)}{dt}\geq\frac{1}{T}(\rho r-\text{Div}\mathbf{q})+\frac{1}{T^{2}}\mathbf{q}\cdot\text{Grad}T, $$

(4)

where T is the absolute temperature.

Substitution of (ρ r−Divq) from (3) to (4) yields

$$ \rho\dot{\eta}+\dot{\rho}\eta\geq\frac{1}{T}(\rho\dot{e}+\dot{\rho}e-\mathbf{P}:\dot{\mathbf{F}})+\frac{1}{T^{2}}\mathbf{q}\cdot\text{Grad}T, $$

(5)

or, written in terms of the internal dissipation,

$$ D_{\text{int}}=\mathbf{P}:\dot{\mathbf{F}}-\rho(\dot{e}-T\dot{\eta})-\dot{\rho}(e-T\eta)-\frac{1}{T}\mathbf{q}\cdot\text{Grad}T\geq0. $$

(6)

We introduce the specific Helmholtz free energy per unit mass

$$ w=e-T\eta, $$

(7)

and, consequently, we have

$$ e=w+T\eta,\quad\dot{e}=\dot{w}+\dot{T}\eta+T\dot{\eta}. $$

(8)

Substituting (8) in (6) we get

$$ D_{\text{int}}=\mathbf{P}:\dot{\mathbf{F}}-\rho(\dot{w}+\dot{T}\eta)-\dot{\rho}w-\frac{1}{T}\mathbf{q}\cdot\text{Grad}T\geq0. $$

(9)

Then, we calculate the Helmholtz free energy increment

$$ \dot{w}=\frac{\partial w}{\partial\mathbf{F}}:\dot{\mathbf{F}}+\frac{\partial w}{\partial T}\dot{T}, $$

(10)

and substitute it in (9) as follows

$$ D_{\text{int}}=\left(\mathbf{P}-\rho\frac{\partial w}{\partial\mathbf{F}}\right):\dot{\mathbf{F}}-\rho\left(\frac{\partial w}{\partial T}+\eta\right)\dot{T}-\dot{\rho}w-\frac{1}{T}\mathbf{q}\cdot\text{Grad}T\geq0. $$

(11)

The Coleman-Noll procedure suggests the following choice of the constitutive laws

$$ \mathbf{P}=\rho\frac{\partial w}{\partial\mathbf{F}},\quad\eta=-\frac{\partial w}{\partial T}. $$

(12)

and, consequently, the dissipation inequality reduces to

$$ D_{\text{int}}=-\dot{\rho}w-\frac{1}{T}\mathbf{q}\cdot\text{Grad}T\geq0. $$

(13)

We further note that the process of the bond breakage is very fast as compared to the dynamic deformation process and the mass density changes in time as a step function. So, strictly speaking, the density rate should be presented by the Dirac delta in time. We will not consider the super fast transition to failure, which is of no interest on its own, and assume that the densities before and after failure are constants and, consequently,

$$\dot{\rho}=\text{Div}\mathbf{s}+\xi=0, $$

or

$$ \text{Div}\mathbf{s}+\xi=0. $$

(14)

Then, the dissipation inequality reduces to

$$ D_{\text{int}}=-\frac{1}{T}\mathbf{q}\cdot\text{Grad}T\geq0, $$

(15)

which is obeyed because the heat flows in the direction of the lower temperature.

It remains to settle the boundary and initial conditions.

Natural boundary conditions for zero mass flux represent the mass balance on the boundary ∂Ω₀

$$ \mathbf{s}\cdot\mathbf{n}=0, $$

(16)

where n is the unit outward normal to the boundary in the reference configuration.

Natural boundary conditions for given traction $\bar {\mathbf {t}}$ represent the linear momentum balance on the boundary ∂Ω₀

$$ \mathbf{P}\mathbf{n}=\bar{\mathbf{t}}, $$

(17)

or, alternatively, the essential boundary conditions for placements can be prescribed on ∂Ω₀

$$ \mathbf{y}=\bar{\mathbf{y}}. $$

(18)

Initial conditions in Ω₀ complete the formulation of the coupled mass-flow-elastic initial boundary value problem

$$ \mathbf{y}(t=0)=\mathbf{y}_{0},\quad\mathbf{v}(t=0)=\mathbf{v}_{0}. $$

(19)

Remark The fact that we ignore the process of the transition to failure and use (14) instead of (1) might be difficult to comprehend at first glance. To ease the comprehension the reader might find it useful to consider the analogy between fracture and the buckling process in thin-walled structure. The pre-buckled and post-buckled states of a structure are usually analyzed by using a time-independent approach. The very process of the fast dynamic transition to the buckled state is of no interest and it is normally ignored in analysis by dropping the inertia terms from the momentum balance equation: $\frac {d(\rho \mathbf {v})}{dt}=\text {Div}\mathbf {P}+\rho \mathbf {b}=\mathbf {0}$ or DivP+ρ b=0. By analogy with the buckling analysis we are only interested in the pre-cracked and post-cracked states while the transition (bond rupture) process can be ignored. The latter is the reason why the mass balance can be written in the simplified form: $\frac {d\rho }{dt}=\text {Div}\mathbf {s}+\xi =0$ or Divs+ξ=0. It is also important to emphasize that the proposed simplification does not affect the natural boundary condition (16). This boundary condition is the expression of the mass balance on the boundary, which is obtained by using the standard Cauchy tetrahedron argument.

Constitutive equations

Constitutive law for the mass source is the very heart of the successful formulation of the theory and the reader is welcome to make a proposal.

We choose, for example, the following constitutive law, whose motivation is clarified below,

$$ \xi(\rho,\rho_{0},w,\phi)=\beta(\rho_{0}H(\zeta)\exp[-(w/\phi)^{m}]-\rho), $$

(20)

where ρ ₀=ρ(t=0) is a constant initial density; β>0 is a material constant; ϕ is the specific energy limiter per unit mass, which is calibrated in macroscopic experiments; m is a dimensionless material parameter, which controls the sharpness of the transition to material failure on the stress-strain curve; and H(ζ) is a unit step function, i.e. H(ζ)=0 if ζ<0 and H(ζ)=1 otherwise.

The switch parameter ζ, which is necessary to prevent from material healing, will be explained below.

Constitutive law for the Lagrangean mass flux can be written by analogy with the Fourier law for heat conduction

$$ \mathbf{s}=\kappa H(\zeta)\exp[-(w/\phi)^{m}]J(\mathbf{F}^{\mathrm{T}}\mathbf{F})^{-1}\text{Grad}\rho, $$

(21)

where κ>0 is a mass conductivity for the isotropic case, which might depend on the deformation process.

The exponential factor is necessary in (21) to suppress diffusion in the failed material.

Substitution of (21) and (20) in (14) yields

$$ \text{Div} \left(l^{2}H(\zeta)\exp[-(w/\phi)^{m}]J(\mathbf{F}^{\mathrm{T}}\mathbf{F})^{-1}\text{Grad}\frac{\rho}{\rho_{0}}\right) +H(\zeta)\exp[-(w/\phi)^{m}]-\frac{\rho}{\rho_{0}}=0, $$

(22)

where

$$ {l=\sqrt{\kappa/\beta}} $$

(23)

is the characteristic length, which might depend on the deformation process.

It is remarkable that we, actually, do not need to know κ and β separately and the knowledge of the characteristic length is enough. For example, the estimate of the characteristic length for rubber is l=0.2 mm (Volokh 2011) and for concrete it is l=2.6 cm (Volokh 2013b).

To justify the choice of the constitutive Eq. (20) for the mass source/sink we note that in the case of the homogeneous deformation and mass flow the first term on the left hand side of (22) vanishes and we obtain

$$ \rho=\rho_{0}H(\zeta)\exp[-(w/\phi)^{m}]. $$

(24)

Substituting this mass density in the hyperelastic constitutive law we have

$$ \mathbf{P}=\rho_{0}H(\zeta)\exp[-(w/\phi)^{m}]\frac{\partial w}{\partial\mathbf{F}}=H(\zeta)\exp[-(W/\varPhi)^{m}]\frac{\partial W}{\partial\mathbf{F}}, $$

(25)

where

$$ W=\rho_{0}w,\quad\varPhi=\rho_{0}\phi $$

(26)

are the Helmholtz free energy and energy limiter per unit referential volume accordingly.

Constitutive law (25) presents the hyperelasticity with the energy limiters - see Volokh (2007, 2013a, 2016) for the general background. Integrating (25) with respect to the deformation gradient, we introduce the following form of the strain energy function

$$ \varPsi(\mathbf{F},\zeta)=\varPsi_{\mathrm{f}}-H(\zeta)\varPsi_{\mathrm{e}}(\mathbf{F}), $$

(27)

where

$$ {\varPsi}_{\mathrm{e}}(\mathbf{F})=\frac{\varPhi}{m}{\varGamma}\left(\frac{1}{m},\frac{W(\mathbf{F})^{m}}{{\varPhi}^{m}}\right), \quad{\varPsi}_{\mathrm{f}}={\varPsi}_{\mathrm{e}}(\mathbf{1}). $$

(28)

Here Ψ_f and Ψ_e(F) designate the constant bulk failure energy and the elastic energy respectively; $\varGamma (s,x)=\int _{x}^{\infty }t^{s-1}e^{-t}dt$ is the upper incomplete gamma function.

The switch parameter ζ∈(−∞,0] is defined by the evolution equation

$$ \dot{\zeta}=-H\left(\epsilon-\frac{{\varPsi}_{\mathrm{e}}}{{\varPsi}_{\mathrm{f}}}\right),\quad\zeta(t=0)=0, $$

(29)

where 0<ε≪1 is a dimensionless precision constant.

The physical interpretation of (27) is straightforward: material is hyperelastic for the strain energy below the failure limit - Ψ_f. When the failure limit is reached, then the strain energy becomes constant for the rest of the deformation process precluding the material healing. Parameter ζ≤0 is not an internal variable. It is a switch: ζ=0 for the reversible process; and ζ<0 for the irreversibly failed material and dissipated strain energy.

For illustration, we present the following specialization of the intact strain energy for a filled Natural Rubber (NR) (Volokh 2010)

$$ {W=\rho_{0}w=\sum\limits_{k=1}^{3}c_{k}(I_{1}-3)^{k},\quad J}=\det\mathbf{F}=1, $$

(30)

where c ₁=0.298 MPa, c ₂=0.014 MPa, c ₃=0.00016 MPa and the failure parameters are m=10, and Φ =82.0 MPa.

The Cauchy stress, defined by σ=J ⁻¹ P F ^T, versus stretch curve for the uniaxial tension is shown in Fig. 2 for both cases with and without the energy limiter. Material failure takes place at the critical limit point in correspondence with tests conducted by Hamdi et al. (2006).

Fig. 2

Uniaxial tension of natural rubber: Cauchy stress [MPa] versus stretch. Dashed line specifies the intact model; solid line specifies the model with energy limiter

For the implications and experimental comparisons of the elasticity with energy limiters, the reader is advised to look through Volokh (2013a; 2016), for example. We completely skip this part for the sake of brevity.

Thus, the proposed constitutive law for the mass source is motivated by the limit case of the coupled formulations in which the deformation is homogeneous.

Conclusion

Crack in a bulk material is not an ideal unzipping of two adjacent atomic layers. It is rather a massive breakage of atomic bonds diffused in a volume of characteristic size. The massive bond breakage is accompanied by the localized loss of material. Thus, material sinks in the vicinity of the crack. Evidently, the law of mass conservation should be replaced by the law of mass balance, accounting for the mass flow in the vicinity of the crack. The coupled mass-flow-elasticity problem should be set for analysis of crack propagation.

In the present work, we formulated the coupled problem based on the thermodynamic reasoning. We assumed that the mass loss related to the crack development was small as compared to the mass of the whole body. In addition, we assumed that the process of the bond breakage was very fast and the mass density jumped from the intact to failed material abruptly allowing to ignore the transient process of the failure development. These physically reasonable assumptions helped us to formulate a simple coupled initial boundary value problem. In the absence of failure localization into cracks the theory is essentially the hyperelasticity with the energy limiters. However, when the failure starts localizing into cracks the diffusive material sink activates via the mass balance equation and it provides the regularization of numerical simulations. The latter regularization is due to the mass diffusion—first term on the left hand side of (22).

The attractiveness of the proposed framework as compared to the traditional continuum damage theories is that no internal parameters (like damage variables and phase fields) are used while the regularization of the failure localization is provided by the physically sound law of mass balance.

A numerical integration procedure for the formulated coupled initial boundary value problem is required and it will be considered elsewhere.

Finally, it should be noted that in the present work we focused on brittle fracture of elastomers, concrete, etc. In the case of ductile fracture of polycrystalline metals, for example, the dislocation-triggered plasticity should be taken into account. The marriage of the mass sink approach and theories of ductile failure would be interesting.