A ductile fracture model based on continuum thermodynamics and damage

A ductile fracture model based on continuum thermodynamics and damage

Journal Pre-proof

A ductile fracture model based on continuum thermodynamics and damage Razanica S., Larsson R., Josefson B.L. PII: DOI: Reference:

S0167-6636(19)30256-X https://doi.org/10.1016/j.mechmat.2019.103197 MECMAT 103197

To appear in:

Mechanics of Materials

Received date: Revised date: Accepted date:

31 March 2019 4 August 2019 3 October 2019

Please cite this article as: Razanica S., Larsson R., Josefson B.L., A ductile fracture model based on continuum thermodynamics and damage, Mechanics of Materials (2019), doi: https://doi.org/10.1016/j.mechmat.2019.103197

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights • A ductile failure model derived based continuum thermodynamics and damage is presented. • The model response draws from a continuum damage degraded effective material model. • As main prototype for the effective material, the Johnson-Cook continuum model is used. • A novel thermo-mechanically motivated damage driving energy is proposed. • A damage threshold controlling the onset of inelastic damage driving energy is used. • The model is computationally efficient and shows well controlled damage progression.

1

A ductile fracture model based on continuum thermodynamics and damage Razanica, S.a,∗, Larsson, R.a , Josefson, B. L.a a

Division of Material and Computational Mechanics Department of Industrial and Materials Science Chalmers University of Technology, SE-412 96 G¨ oteborg, Sweden

Abstract The paper presents an approach to ductile failure modeling derived based on continuum thermodynamics and damage. A continuum damage enhanced formulation of the effective material is used to describe the degradation of the response. From the thermo-mechanically motivated dissipation rate, a novel damage driving energy that involves both stored energy and dissipative contributions due to inelasticity is presented. This damage driving energy is combined with a damage threshold that controls the onset of inelastic damage driving dissipation. The damage evolution law is formulated as a balance of produced dissipation and produced fracture area induced dissipation, involving damage progression velocity and length-scale parameters without any non-local gradient term. As main prototypes for the effective material and damage threshold, the Johnson-Cook continuum and failure models are exploited. From the verification examples, satisfactory convergence properties of the model are obtained. The model capabilities to represent real thermodynamic ductile failure processes are demonstrated for a specimen made of a ductile Weldox steel subjected to tensile split-Hopkinson Tension Bar tests. The model is computationally efficient and shows well controlled damage progression in the FE-application. Keywords: Continuum damage; Ductile fracture; Damage driving energy; Johnson-Cook

1. Introduction Machining is a complex process but also one of the most common manufacturing processes, setting the specific geometrical dimensions and surface finish for the final product. To optimize the machining process, a FE-based computational failure mechanics approach is adopted, cf. refs. [1], [2], which qualitatively models large inelastic deformations at high strain rates and elevated temperatures occurring in the cutting region. The FE-modeling of the machining process is strongly related to the ability of the constitutive and failure models to represent the complex phenomena including localized damage that occur in the vicinity of the cutting tool [3], [4]. These complex phenomena along with the localized failure patterns ∗

Corresponding author: [email protected] (Professor Ragnar Larsson)

Preprint submitted to Elsevier

October 10, 2019

can be modelled using thermo-viscoplasticity combined with a damage evolution model. In this context, state of the art in many commercial platforms is that ”damage” is modelled as an instant drop of the stress (or element deletion) whenever the failure threshold signals, cf. e.g. [5], [6] for different approaches to define the failure threshold. However, this procedure is problematic since it does not respect the fracture mechanics of the ductile failure and it leads to pathological FE-mesh dependence, cf. e.g. [7]. Rooted in the early work by Gurson [8], ductile damage is often attributed to the growth of micro voids. The evolution of ductile damage is typically formulated in (local) Gurson type plasticity models coupled to void volume growth, assumed to coalesce into fracture surface in the FE-application. Various formulations of void growth have been considered, for existing and nucleated voids and the morphology of the voids, e.g. Needleman and Tvergaard [9], Pardoen and Hutchinson [10]. As is often the case of such models, localization of deformation corresponding to void coalescence is highly mesh dependent in the FE-application. Extensions from the local Gurson model concept has therefore been made by, e.g. Tvergaard and Needleman [11], Reusch et al. [12], to non-local formulations; however, without completely convincing capability to represent truly mesh convergent behavior. Another way to represent ductile fracture is to generalize recent developments of continuum phase field models for brittle fracture processes, e.g. cf. [13], [14], [15], [16], into the corresponding one for ductile fracture, e.g. [17], [18]. In Miehe et al. [18], the proposed phase field modeling is introduced into a thermo-plasticity framework in conjunction with a non-local damage model. Furthermore, it has also been argued that viscous regularization of the continuum material model coupled to damage via visco-plasticity may remove the pathological mesh dependence, [19], [20]. However, in [20] it was found that the effective visco-plastic regularization length decreases with increasing damage or deformation. Furthermore, Borvik et al. [21] proposed a coupled constitutive model with respect to viscoplasticity and ductile damage for the representation of penetration and impact problems. To include damage in the constitutive equations, a modified version of the Johnson-Cook (JC) constitutive model [22] is used in [21] to represent the effective material along with the assumption of strain equivalence [23]. The modeling response shows good agreement with the experimental results; however, there is no convergence study with respect to mesh refinements. In this contribution, a damage enhanced continuum of Lemaitre type [23] is proposed to capture ductile fracture processes, where the effective (or undamaged) material response is described by the visco-plastic JC-model, cf. [7],[24]. To obtain the proper thermodynamic coupling, the damage modeling framework is formulated based on continuum thermodynamics, where the dissipated energy due to the adopted damage formulation for ductile failure respects the basic postulates of positive entropy flux and energy balance, cf. also ref. [18]. Here, the resulting mechanical dissipation is subdivided into a degraded effective contribution and a damage driving dissipation, involving elastic and inelastic contributions in terms of the damage threshold, cf. [7]. The inelastic contribution is controlled by the damage threshold that is based on the JC failure model [5]. This threshold yields a more physically based ductile fracture response as compared to the developments in [18], where ”barrier functions” for the damage driving energy are defined based on e.g. critical plastic strain or 3

stress levels for fracture. Another approach is adopted in [17], where the plasticity coupling to the damage driving energy is formulated directly in the damage degradation function. As to the adopted damage evolution law herein, a major novel contribution is the formulation of input dissipation rate to the fracture area development based on a thermomechanically motivated damage driving energy, accounting for elastic and inelastic effects. This is a clear development from brittle (elastic) failure in [25] to ductile failure modeling, where the damage-plasticity coupling is described by the damage driving energy. To maintain computational efficiency of the model, the gradient damage dependence is neglected in the fracture area production as done in [25]. In this way the fracture area production is completely local, consisting of a damage induced convective part and another damage rate induced contribution. As opposed to the rate independent case corresponding to infinite fracture surface progression, we find in the FE-application that the damage field exhibits a damage pattern that generally localizes in a consistent way within a fixed damage width, depending on the fracture area propagation velocity and the internal length parameters. Indeed the local character of the damage evolution law contributes to computational efficiency and robustness. To highlight the significance of the novel damage driving energy formulation, the response and consequent energy dissipation in a 1D tension bar is considered as a start-up example. To verify our proposed ductile damage model, the convergence behavior (including mesh convergence) is studied for a plane strain sheet loaded in tension. The significance of the proposed damage evolution law as compared to the corresponding ”static” one for a hollow plane strain sheet is also considered. In this example, the effect of the strain rate and temperature dependencies in the continuum model are studied as well. In addition, model validation is carried out against tensile split-Hopkinson Tension Bar tests for preheated specimens of a ductile Weldox 460 E material in [26]. 2. Scalar damage enhanced material model representing ductile fracture In this section a continuum damage enhanced thermo-dynamical framework for ductile fracture is outlined based on a scalar damage enhancement of the effective (undamaged) material response, cf. also [27] and [28] for more details. In order to define the visco-plastic response of the effective material, the ”effective” energy dissipation rate is developed to handle high deformation rates, hardening and temperature effects. 2.1. Continuum damage enhanced thermoplasticity framework We consider particles, with position X in the reference configuration B0 , in motion with the non-linear deformation ϕ[X, t] from X ∈ B0 to positions x ∈ B in the current configuration. According to standard notation, F = ∂ϕ/∂X is the deformation gradient and J = det[F] = ρ0 / ρ is the volume deformation measure, ρ0 and ρ are the densities in reference and current configurations, respectively. Consider also the deformation gradient subdivided into an elastic component Fe , a purely volumetric thermal component Fθ and the irreversible inelastic part Fp , forming the total deformation gradient multiplicatively as F = F · Fp

with F = Fe · Fθ 4

1/3

and Fθ = Jθ 1

(1)

where Jθ = det [Fθ ] = 1 + αθ (θ − θ0 ), αθ is the coefficient of thermal expansion, θ > 0 is the absolute temperature, θ0 is a reference temperature and 1 is the second order unit tensor. In Eq.(1), the “bar” indicates that F is a representative of the reversible deformation, including elastic and thermal components. Let us also introduce the spatial velocity gradient l = F˙ · F−1 , which in view of the multiplicative split induces an additive decomposition of l in recoverable and inelastic portions l and lp as l = ¯l + lp

−1 with l = F˙ · F

and lp = F · Lp · F

−1

(2)

where Lp =: F˙ p · F−1 p is the plastic velocity gradient. Related to the non-symmetric spatial velocity gradient l, let us also introduce the following symmetric total and plastic rate of deformation tensors as d = lsym ,

  −1 ¯ ·b ¯ = −Lv b dp = lsym p

  ¯ = F · F−1 · b˙ · F−t · Ft with Lv b

(3)

where the •sym denotes the symmetric portion of the second order tensor • and Lv [•] is the ¯ = F · Ft was also introduced. Lie derivative. Here, the reversible Finger tensor b As a guide to the model developments, we consider the basic postulates of energy balance and positive entropy flux S˙ − Qθ , written in the reference configuration B0 as    Z Z Z  h · ∇ dB ≥ 0 (4) ρ0 edB ˙ = (τ : l − Jh · ∇)dB, S˙ − Qθ = ρ0 s˙ + J θ B0 B0 B0 where ρ0 is the density, e is the internal energy, s > 0 is the entropy density and h is the spatial heat flux vector. Moreover, τ is the Kirchhoff stress tensor and ∇ is the spatial gradient operator. Upon combining the local formats of the relations in Eq.(4) one obtains ρ0 e˙ = τ : l − Jh · ∇,

τ : l − ρ0 (e˙ − θs) ˙ − Jθ−1 h · ∇θ ≥ 0

(5)

which are sometimes also named the energy equation and the Clausius Duhem inequality, respectively. Let us also introduce the Helmholtz free energy ψ and the internal energy e with the argumentation     ¯ k, α, s , ψ = ψ b, ¯ k, α, θ e = e b, (6) where k is an internal variable accounting for isotropic material hardening and 0 ≤ α ≤ 1 is the scalar internal damage variable that represents material degradation in terms of the damage degradation function 0 ≤ f [α] ≤ 1 defined as f [α] = (1 − α)2

(7)

The relation between e and ψ is determined via the Legendre transformation     ¯ k, α, s = inf ψ b, ¯ k, α, θ + sθ e b, θ 5

(8)

where

 ∂ψ ∂e ∂e ˙ ∂ψ ˙ ¯ k, α , θ= and :A= : A with A = b, ∂θ ∂s ∂A ∂A Based on this relationship, elimination of the internal energy in Eq.(5) yields   ˙ ˙ τ : l − ρ0 ψ + sθ −Jθ−1 h · ∇θ ≥ 0 | {z } s=−

(9)

(10)

D≥0

It is assumed (as usual) that the mechanical dissipation rate D and the purely thermal portion are independent dissipative processes. The thermal portion −h·∇θ ≥ 0 is automatically satisfied by Fourier’s law of heat conduction, i.e. h = −kθ ∇θ where kθ is the coefficient of heat conduction. As to the mechanical contribution D ≥ 0 to the dissipation rate we have that   ˙ ˙ D = τ : l − ρ0 ψ + sθ = τ : dp + κk˙ + Aα˙ ≥ 0 (11) corresponding to the constitutive state equations of the Kirchhoff stress τ , the micro hardening stress κ and the elastic damage driving energy A defined as ¯ · ∂ψ , τ = 2ρ0 b ¯ ∂b

κ = −ρ0

∂ψ , ∂k

A = −ρ0

∂ψ ∂α

(12)

Upon combining Eq.(11) with Eq.(12) the energy equation Eq.(5a) is then reformulated as 2 2 ¯ · ∂ ψ : d − 2ρ0 θb ¯ · ∂ ψ : dp + θ ∂κ k˙ − Jh · ∇ cv θ˙ = D + ρ0 2θb ¯ ¯ ∂θ ∂θ∂ b ∂θ∂ b

(13)

2

where cv = −ρ0 θ ∂∂θψ2 is the specific heat capacity. 2.2. Scalar damage enhanced viscoplastic effective material model A scalar damage enhanced viscoplastic effective material model is formulated to promote damage evolution by shear deformation in the effective material. The model is obtained by considering the isochoric part of the effective material degraded with the degradation function f [α] in Eq.(7) as       ¯ + g[θ]ψˆmic [k] + ρ0 ψˆvol J, ¯ θ + ρ0 ψˆth [θ] ψ = f [α]ρ0 ψˆiso b (14)

where ψˆiso is the isochoric part of the free energy, ψˆmic represents stored energy due to internal (isotropic) hardening processes in the material, ψˆvol is the stored free energy due to volumetric deformation and ψˆth is the thermal part describing stored free energy induced by the temperature. Moreover, J¯ = Je Jθ = JJθ denotes the reversible volumetric deformation (here assumed to be without any inelastic component) and the notation •ˆ denotes quantities of the effective material. To more specific, the visco-plastic JC constitutive model [22] is 6

chosen to be the prototype for the effective material response. Pertinent to the JC model a temperature degradation function g[θ] degrades the micro hardening in Eq.(14) defined as   θ ≤ θtrans  1   g[θ] =

1−

  0

θ−θtrans θmelt −θtrans

m

θtrans < θ < θmelt

(15)

θ ≥ θmelt

where m is an exponent of the temperature degradation, θtrans and θmelt are room and melting temperatures, respectively. Moreover, the stored energy functions in Eq.(14) are explicityly defined as
1 ¯ iso − 3 , ρ0 ψˆvol = 1 K log[J] ¯ 2 − 3Kαθ (θ − θ0 ) log[J] ¯ ρ0 ψˆiso = G 1 : b 2 2    (16) B θ mic 1+n th ˆ ˆ (−k) , ρ0 ψ = ρ0 cv θ 1 − log ∗ ρ0 ψ = 1+n θ ¯ is the isochoric portion of b. ¯ In (16) G is the shear modulus, K is the ¯ iso = J¯− 23 b where b bulk modulus and θ0 is the reference temperature. Note that θ0 6= θtrans generally. In-line with the JC-model, B is the isotropic hardening parameter, k represents the plastic strain (to be explicitly defined later) and n is the hardening exponent. As to the thermal part, cv is the specific heat capacity of the material and θ∗ is an arbitrary temperature. It may be noted that ψˆvol does not involve any inelastic component and ψˆiso is temperature independent, i.e ¯ = J − 32 be . ¯ iso = J¯− 23 b b In view of Eq.(12) and Eq.(14), the continuum stress, the internal hardening stress and the elastic damage driving energy are obtained as   iso iso mic 0 ˆ ˆ ¯ τ = f [α]ˆ τ − Jp1, κ = f [α]ˆ κ, A = −f [α]ρ0 ψ + g[θ]ψ (17) where τˆ iso is the isochoric Kirchhoff stress tensor, κ ˆ is the micro hardening stress and p is the pressure of the effective material. These state equations correspond to the dissipation rate written as ˆ + Aα˙ ≥ 0 D = f [α]D

with

ˆ = τˆ iso : dp + κ D ˆ k˙ ≥ 0

(18)

ˆ is the portion of D that corresponds to the dissipation in the effective material, cf. where D also the developments in refs [7] and [24]. In view of Eq.(16), the explicit expressions of the effective stress response (=response of the undamaged material) for our prototype JC-model become ∂ ψˆiso ¯ ∂ ψˆmic iso dev ¯ ˆ τ = 2ρ0 ¯ · b = Gb , κ ˆ = −ρ0 g[θ] = g[θ]Bk n ∂k ∂b (19) vol ˆ
∂ψ −1 ¯ − 3αθ (θ − θ0 ) p = −ρ0 ¯ = −J¯ K log[J] ∂J where it may be noted that the pressure response is ”undamaged” during the entire damage process, i.e. p = pˆ. 7

2.3. Viscoplastic evolution rules for the effective material ˆ ≥ 0 and To ensure positive dissipation rate D ≥ 0 in Eq.(18), it suffices to consider D α˙ ≥ 0 (since A > 0 always). Following the JC model , the static “von Mises” yield function is introduced as r 3 iso |ˆ τ | (20) φs = τˆe − (g[θ]A + κ ˆ ) , τˆe = 2 where τˆe is the von Mises stress of the effective stress response and A is the initial yield stress. Here it is assumed that g[θ] degrades the yield stress A in the same way as the effective hardening in Eq.(17). To account for strain rate hardening effects according to the JC model , a rate dependent yield function is formulated in terms of the quasi-static one as   ˙ φ = φs − (g[θ]A + κ ˆ ) C log (21) ˙0 where C is a viscosity parameter, ˙ ≥ ˙0 is the effective strain rate (cf. Eq.(25) below) and ˙0 is the reference strain rate parameter. A constrained maximum dissipation principle in the effective material is formulated as    ˙ max iso ˆ = iso ≥0 (22) D τˆ : dp + κ ˆ k˙ − λφ = λ g[θ]A + (g[θ]A + κ ˆ ) C log τˆ , κ ˆ, λ ˙0 where the last equality corresponds to the optimality conditions in plastic loading ˆ ∂D ∂φ = dp − λ iso = 0, iso ∂ τˆ ∂ τˆ

ˆ ∂D = k˙ + λ = 0, ∂ˆ κ

∂D =φ=0 ∂λ

(23)

Hence, the evolution equations are obtained as dp = λ

∂φ 3 τˆ iso , k˙ = −λ with φ ≤ 0, λ ≥ 0, φλ = 0 = λ 2 |ˆ ∂ τˆ iso τ iso |

(24)

Here, the evolution rules for dp and k˙ are complemented by the Karush-Kuhn-Tucker conditions to distinguish elastic and plastic loading. We also choose ˙ = −k˙ = λ ≥ ˙0 . In view of Eq.(21), these conditions are satisfied by the JC-viscoplastic evolution rule   1 < φs > ≥ ˙0 ≥ 0 (25) ˙ := λ = ˙0 exp C g[θ]A + κ ˆ where < • > is the positive part function (or Macaulay brackets). Here it may be noted that the reference strain parameter is generally a quite small number; however, to handle the situations at very slow loading rates λ < ˙0 , we set C = 0 in Eq.(21) corresponding to rate independent response as done in [24]. In view of Eq.(18), the reduced effective dissipation rate is obtained as ˆ = λ (g[θ]A+ < φs >) ≥ 0 D (26) 8

2.4. Energy equation revisited We are now in position to establish the explicit expression for the energy equation Eq.(13), where it is concluded in view of Eq.(16) that 2 2 ¯ · ∂ ψ : dp = 0 ¯ · ∂ ψ = −3Kαθ θ1 ⇒ 2ρ0 θb ρ0 2θb ¯ ¯ ∂θ∂ b ∂θ∂ b

(27)

since dp is deviatoric in view of Eq.(24). With the total dissipation rate in Eq.(11) rewritten in terms of the plastic work wp = f [α]ˆ τe λ, i.e. D = wp + κk˙ + Aα, ˙ one obtains   ∂κ p k˙ − Jh · ∇ (28) cv θ˙ = w + Aα˙ − 3Kαθ θ1 : d + κ + θ ∂θ Furthermore, it is assumed that adiabatic conditions prevail (h · ∇ = 0), whereby a simplified formulation for the energy equation is obtained as   ∂ˆ κ cv θ˙ = w + Aα˙ − 3Kαθ θ1 : d − f [α] κ ˆ+θ λ ∂θ p

(29)

where it may be remarked that elastic damage also contributes to the heat generation, whereas the thermal stress Kαθ θ and the effective hardening κ ˆ require a portion of the energy to develop. 2.5. Continuum damage enhancement of the effective material This sub-section describes the continuum damage degradation of the effective material in terms of a progressive damage evolution model. A key issue in this model is the damage input energy to the fracture area progression, in terms of both elastically stored and dissipated energies. In order to maintain the proper control of the input energy to the damage process, a damage threshold is introduced. The model is derived from the balance between input damage driving energy from the deformation process and energy dissipation due to produced fracture area. 2.5.1. Damage driving energy with damage threshold In order to allow for a thermo-mechanically motivated damage driving energy that accounts for inelastic deformations, the dissipation R t rate D ≥ 0 in Eq.(18) is reformulated, cf. ˆ ˆ dt as also [7], in the total effective dissipation DT = 0 D .

ˆ T + (A + B) α˙ ≥ 0 with B = −f 0 [α]D ˆT ˆ + Aα˙ = f [α]D D = f [α]D

(30)

where the damage driving energy includes the elastic and inelastic contributions A and B, respectively. The induced elastic contribution A is initiated instantaneously as dictated by Eq.(18), whereas the inelastic contribution B is subdivided into an initial part B i and the 9

damage driving part B f . Here the initial part B i is considered as the portion of the effective dissipation that is required to generate onset of inelastic induced damage production. In ˆ T is subdivided additively written order to clarify this point, the total effective dissipation D as   ˆ Ti + D ˆ f = Bi + Bf B = −f 0 [α] D (31) T where initial B i -portion and failure portion B f of B are separated by the damage threshold time t = tf defined as ˆ f [t] , D ˆ Ti = ˆT = D ˆ Ti [tf ] + D D T

Z

0

tf

ˆ dt , D ˆ f [t] = D T

Z

t

tf

ˆ dt D

(32) .

.

ˆ T = f [α]D ˆf − Upon combining the relations Eq.(30), Eq.(31) and Eq.(32), we find that f [α]D T B i α˙ for t ≥ tf . In view of Eq.(30) the dissipation rate is finally obtained in the piece-wise form  ˆ+ f [α]D A α˙ 0 ≤ t < tf   |{z}   el. dam. driv. .
D[t] = (33) ˆ f + A + B f α˙ t ≥ tf  f [α]D  T  | {z }  el.-pl. dam. driv.

where it is emphasized that different expressions (elastic or elastic-plastic) for the damage driving energy are obtained before and after the threshold time has been reached. Different damage processes can be monitored from Eq.(33): 1) elastic damage corresponding to tf → ∞ and 2) elastic-plastic damage for tf → 0. In order to allow for a flexible modeling of ductile fracture, a finite threshold time is introduced for the onset of inelastic damage. In practice, the damage threshold is based on e.g. JC failure criterion [5] or the accumulated plastic work using the Cockcroft-Latham model [29] to mention two possibilities. The implications of these cases are discussed in the first numerical example. 2.5.2. A progressive damage evolution model A progressive damage evolution model in line with ref. [25] is advocated. In this model it is first noted that the damage field α ∈ B0 may exhibit portions in B0 with both diffuse and localized characters. Let us also assume that the damage field is localized within a finite localization band of width lc , as shown in Figure 1. To describe the diffuse fracture zone, the fracture area functional Af is introduced in terms of the fracture area density function γ[α, α] ˙ representing the progression and localization of the damage field α ∈ B0 defined as  Z Z α 0 α˙ 0 α + ∗ dα0 (34) Af = γdB with γ := lc v B0 0 where the internal length parameter lc describes the diffuse character of the fracture area and v ∗ is the fracture area progression speed parameter, controlling the damage evolution 10

α. ˙ In the damage zone we have from Eq.(34) the area progression γ˙ =

1 1 αα˙ + ∗ α˙ 2 ≥ 0 lc v

(35)

The interpretation is shown in the close-ups in Figure 1, considering isotropic fracture area progression in the microstructure with the representative length l < lc . From Eq.(35) the fracture area production is represented by the term αα/l ˙ c due to growth (or convection) of existing damage induced voids, whereas the term α˙ 2 /v ∗ represents void nucleation. Clearly, the conventional rate independent fracture area evolution is obtained when the damage propagation speed parameter v ∗ → ∞, cf. also refs. [7], [25]. To complete the damage evolution model, a global balance relationship is established between produced energy dissipation rate due to the fracture surface area production and input damage dissipation rate from the damage driving energy in Eq.(33). To this end, the fracture energy dissipation functional due to fracture area development is G = Gc Af , where Gc is the crack surface energy release rate. As to input of disspated energy, the inelasticity driven damage driving energy dissipation rate developed in sub-section 2.5.1 is considered. As suggested, it is the damage driving part of Eq.(33) that is the input energy dissipation rate to the damage process. From Eq.(34) and Eq.(35), the global balance relationship is thus formulated as Z Z Gc ˙ = AT αdB ˙ Gc γ˙ = ∗ (v ∗ α + lc α) ˙ α˙ = AT α˙ (36) G = Gc γdB ˙ v lc B0 B0 {z } | {z } | input

produced

where AT is the damage driving energy that becomes different depending on whether the threshold t = tf for the onset of inelastic damage has been reached. In view of Eq.(33), the damage driving energy is defined as ( A 0 ≤ t ≤ tf AT = (37) f A + B t > tf

In view of the balance relation (36) between produced and input dissipation rates to the damage zone, a maximum dissipation rate formulation is considered. Like in Eq.(21), the dynamic damage loading function is φα = AT −

Gc ∗ (v α + lc α) ˙ v ∗ lc

(38)

Whenever damage loading occurs, a constrained maximum of the damage driving portion of the total dissipation rate is formulated as Dα := AT α˙ =

Gc max AT α˙ − µφα = ∗ (v ∗ α + lc α) ˙ α˙ AT , µ v lc 11

(39)

γ[α] .

α=0 .

α=0

n

lc

.

.2

γ = αα + α* lc v .

.

B0

α>0 −

lc 2

lc 2

n

.

v*

αα lc

v* lh < lc

v* . α2 v*

v*

lh Figure 1: Localized damage progression zone in B0 ; the first close-up shows the fracture area progression within damage localization zone of internal width lc . The second close-up shows the interpretation of damage area progression in terms of void growth for existing and nucleated voids in the localization zone. Note that the representative length l of the material microstructure represents a portion of the localized damage width lc .

corresponding to

∂Dα ∂φα = α˙ − µ = 0, ∂AT ∂AT 7

∂Dα = φα = 0 ∂µ

(40)

To describe the damage loading situations, the Karush-Kuhn-Tucker conditions are specified φα ≤ 0, α˙ ≥ 0, φα α˙ = 0

(41)

Respecting these damage loading conditions and Eq.(38), the damage evolution is obtained in a simple format as lc α˙ = v ∗ < αs [α] − α >

with αs =

AT [α]lc Gc

(42)

where αs represents source of damage energy from elastic-plastic deformation. Please note once again that the damage model corresponds to a local balance law for the damage production, where (in the case of damage loading) the local damage production lc α˙ and the 12

convection v ∗ α balance input damage production v ∗ αs . In particular, static damage development corresponds to v ∗ → ∞ leading to α = αs [α]. In the case of damage unloading α˙ = 0. We also note that an additional spatial (or non-local) production term may be added to the balance relation Eq.(42) in terms of gradient damage dependence, cf. e.g. Larsson et al. [25]. 3. Representative numerical examples 3.1. Preliminaires To study the capabilities of the proposed damage enhanced JC-model, a set of representative numerical examples highlighting the model behavior, the model convergence properties and model validation are considered. To assess the model behavior, a uniaxial tensile test of ductile damage development at the material point level is studied. In particular, we demonstrate the implications of the novel concept for elastic-plastic damage driving energy on the evolution of damage. A plane strain sheet with a central imperfection at tensile extension is used for the model verification. Model sensitivities regarding mesh refinements and choices of the fracture area propagation velocity are studied. In addition, the damage evolution law is investigated for a hollow plane strain sheet. Here, the significance of the void nucleation term (α˙ 2 /v ∗ ) in the damage evolution law is compared to the static damage evolution (α = αs ). For the model validation, experimental results for a tensile split-Hopkinson test of circular test specimens from [26] with varying initial temperature are modeled. The JC-model along with the damage enhancement is implemented for dynamic loading in the commercial finite element software ABAQUS/ExplicitTM via a VUMAT user subroutine. For computational efficiency of the explicit dynamic implementation, the hypo-inelastic concept is used based the Green-Naghdi stress rate. In the present context of isotropic elasticity/inelasticity and “small” elastic strains, this concept may be shown to be equivalent to the outlined hyper-inelastic formulation, cf. ref. [24]. As to the temperature generation, the adiabatic condition is assumed whereby Eq.(29) is used to determine the temperature. We consider the elastic damage driving energy as A = −f 0 [α]ρ0 ψˆiso . The assumption is thus that the stored inelastic free energy in Eq.(17) does not contribute to the damage evolution. Another key issue for the damage drive in ductile failure is the damage threshold t = tf defining the onset of inelastic damage driving energy B f . For this, the JC failure criterion [5] is considered, stating that t = tf is obtained whenever the equivalent plastic strain k approaches the plastic failure strain pf = −k = h[r, λ, θ, t]. The threshold time is thus defined via k + h [r, λ, θ, tf ] = 0

     θ − θtrans p λ 1 + d5 with r = − h = (d1 + d2 exp [−d3 r]) 1 + d4 log ˙0 θmelt − θtrans τˆe

(43)

The parameters d2 and d3 account for pressure dependent initiation via the “triaxiality” ratio r in (43); moreover, d4 accounts for the strain rate effect, whereas the d1 -parameter relates to the amount of plastic deformation that is required to initiate onset of damage evolution 13

and d5 relate to temperature dependence. In practice, this criterion is never exactly fulfilled. It has to be checked at the end of each time step whether t > tf is signaled. 3.2. Uni-axial tensile test of ductile damage at material point level In this example, a uniaxial tensile test of ductile damage at material point level is carried out to study the relevance and implications of the damage threshold. The example considers a pearlitic steel with the JC material and failure parameters given in Table 1, with an applied 1D loading rate ˙ = 1000 s−1 . In turn, these parameters were calibrated against Table 1: Johnson-Cook material and failure parameters for the pearlitic plane strain sheet from [30].

A [MPa] B [MPa] 550 500

C 0.0804

n 0.4

m 1.68

θmelt [◦ C] 1400

d1 0.25994

d2 0.61368

d3 2.5569

d4 -0.027652

d5 0.6

˙0 [s−1 ] 0.001

E [GPa] 190

ν 0.3

ρ [kg/m3 ] cv [J/Kg◦ C] θtrans [◦ C] αθ [mm/mm◦ C] 7850 452 20 12*10−6 uniaxial tensile test

Out[115]=

σ [MPa]

1500

1000

elastic-plastic damage, no threshold, t f =0

500

threshold defined by JC fracture criterion elastic damage, t f →∞ ϵ =ϵ f

0 0

50

100

150

200

ϵ [%]

Figure 2: Uniaxial Cauchy stress σ response versus Lagrange strain  for tensile failure representation of pearlite. Different response curves are obtained depending on the choice of damage threshold time tf .

experimental data for pearlite in ref [24]. As to the fracture energy release rate Gc and the internal length lc , it appears that these parameters control the damage energy dissipation Dα , whereas the fracture area progression speed v ∗ stabilizes the FE-solution sufficiently. 14

For computational efficiency v ∗ should be chosen as “high” as possible without seriously influencing the response, cf. [25]. In this example v ∗ = 1250000 mm/s is used, which is in the order of the crack propagation speed. In order to choose a reasonable value of the fracture energy for the pearlitic steel, we 2 consider the relation Gc = KIc /E 0 where E 0 = E/(1 − ν 2 ) for plane strain and KIc is the fracture toughness. Here, the fracture toughness of pearlitic steel (AISI 1080) from [31] is employed, whereby the fracture energy release rate is Gc = 12 N/mm. Central to the damage dissipation is the internal length parameter lc via the ratio lc /Gc in Eq.(42), which for a given fracture energy Gc increases the ductility as lc is reduced. In particular, before the damage threshold has been reached (or for elastic damage when tf → ∞), it appears that lc controls the amount of damage induced by the elastic damage driving energy. In the present paper, we take on the view that the influence of elastic damage driving energy should be kept as small as possible. To maintain this effect the lc -value cannot be too large. On the other hand, in order to obtain fast stress degradation pertinent to pearlite, the lc -value cannot be too small. In the present example, a representative value is lc = 0.05 mm. Another way to avoid the issue with the elastic damage driving energy, and to allow for a flexible choice of lc , is simply to set AT = B f , thus neglecting A as is done in the last example. However, in this case thermodynamic consistency is lost, as discussed in conjunction with the dissipation energy developed during a ductile failure process in the uniaxial tensile test below. uniaxial tensile test

0.8

0.6

Out[119]=

α[t]

elastic-plastic damage, no threshold, t f =0 threshold defined by JC fracture criterion elastic damage, t f →∞

0.4

0.2

0.0 0.0

0.5

1.0

1.5

2.0

t [ms]

Figure 3: Damage development versus time during tensile failure of pearlite. The different damage curves are obtained because of the different choices of the damage threshold time tf . Note the limited damage generated by the elastic stored energy obtained for tf → ∞.

Figure 2 shows the uniaxial Cauchy stress response versus the Lagrange strain for three different (extreme) cases of damage evolution. In the first extreme case tf is zero, corresponding to complete elastic-plastic damage and completely brittle response. Note the obtained 15

uniaxial tensile test

400

elastic-plastic damage, no threshold, t f =0 threshold defined by JC fracture criterion elastic damage, t f →∞

Out[123]=

θ[t]

300

200

100

0 0.0

0.5

1.0

1.5

2.0

t [ms]

Figure 4: Adiabatic temperature increase versus time during ductile failure of pearlite.

“fast” decay of the stress that resembles “instantaneous stress drop” obtained when using the element removal technique. The corresponding damage evolution and adiabatic temperature increase are shown in Figures 3 and 4. A very small temperature increase is obtained because of the rapid stress drop. The other extreme case is that of pure elastic damage driving energy AT = A corresponding to tf → ∞. Note also the limited effect of the elastic damage driving energy, cf. Figure 3, for the present choice lc = 0.05 mm. As a consequence the uniaxial response is governed by the JC-model solely, where a slight thermal softening effect due to the temperature increase 500◦ C is obtained. The intermediate case is defined by the JC failure criterion Eq.(43) predicting f = 45 % for the tensile test. The stress response up to the damage threshold point is followed by the damage induced steep stress decay shown in Figure 2. The drop of the stress is reflected by the damage evolution and temperature increase in Figures 3–4. Note the desired limited elastic induced damage in the beginning of the deformation process. In this case the adiabatic temperature increase is not as pronounced as for the case of pure elastic damage. Let us also highlight the damage driving energy concept outlined in sub-section 2.5.1. To illustrate this concept, consider the total mechanical energy dissipation from the integrated dissipation rate in Eq.(33) as (R t R Z t f ˆ dt + tf Aα˙ dt f [α] D 0 ≤ t ≤ tf 0 0
Rt DT [t] = D dt = (44) f f ˆ [t] + DT [tf ] + f [α]D A + B α ˙ dt t ≥ t f 0 T tf Associated with Eq.(44), Figure 5 shows the integrated dissipation DT when the damage threshold is determined by the JC failure criterion. Prior to the onset of inelastic damage, it may be noted that some damage is initiated due to the induced elastic damage driving energy A, whereas B f = 0 in this range. After the onset of inelastic damage, the total 16

ˆ f and the damage dissipation is characterized by the degraded effective dissipation f [α]D T driving part, consisting of both elastic and inelastic contributions. . In view of the thermodynamic consistency. and Figure 5, it is noted that DT = D ≥ 0
ˆ f is generally undetermined. Here, it may and A + B f α˙ ≥ 0, whereas the sign of f [α]D T .

ˆ f < 0 corresponding to an under predicted be remarked that it is inconsistent that f [α]D T damage driving part in Eq.(33). However, as α → 1 the degraded effective dissipation ˆ f → 0. Hence, from Eq.(44) the entire dissipation goes into the damage driving part f [α]D T of the dissipation, i.e. Z t
A + B f α˙ dt (45) DT = DT [tf ] + tf

This behavior is clearly obtained in the uniaxial case shown in Figure 5. However, the observed convergence between the total and the damage driving dissipation hinges on that both the elastic and the inelastic damage driving contributions are accounted for. 3.3. Tensile loaded plane strain sheet with a central imperfection

The tensile loaded plane strain sheet in Figure 6 is used to study the sensitivity of the finite element solution for the parameters of the damage evolution law in Eq.(42). Besides this study, mesh convergence of the response for varying element size le is also considered. As shown in Figure 6, a slight in-homogeneity is provoked by an imperfection introduced in a 2x2 mm square in the center of the sheet with a degraded yield stress. In this example, Dissipation energies 600

t=t f 500

DT [t f ]

Out[43]=

DT [N/mm]

400

300

DT [t] from Eq.(44) t  DT [t f ]+∫t (+ℬ f )αⅆt

200

f

f f[α]DT 100

0 0.0

0.5

1.0

1.5

2.0

t [ms]

Figure 5: Dissipative energy versus time during ductile failure process in the uniaxial tensile test. Note the decay of the dissipation rate D = D˙ T induced by increasing damage.

17

v, R

0.7σ 0

50mm

2 × 2mm

50mm

Figure 6: Plane strain sheet with a central imperfection and boundary conditions.

4-node plain strain elements with 1-point integration and hourglass control are used. The tensile loaded sheet is free to move horizontally at the bottom surface except for the fully constrained lower left corner. The upper boundary is subjected to a vertical deformation velocity of v¯ = 2500 mm/s corresponding to high strain rates in the sheet. The same pearlite parameters and internal length lc = 0.05 mm as in the previous example are used. In view of the boundary conditions of this example, an almost instantaneous drop in the forcedisplacement response is obtained. The steep decay of the reaction force is manifested by damage initiation in the center of the sheet followed by a damage progression into a localized shear mode. Please note that the lc –parameter determines the width of the localization zone in a qualitative fashion. It turns out that the actual width of the localization also depends on the fracture energy release Gc and the boundary conditions. 3.3.1. Fracture area progression speed and damage mode To understand the influence of the fracture area velocity, a parametric study is conducted 2 for the v ∗ -parameter when keeping the FE-discretization (le = 0.256 mm) and the internal length (lc = 0.05 mm) fixed. The results are shown in Figure 7 for a range of propagation velocities covering rate dependent (v ∗ = 10000 mm/s) and almost rate independent cases (v ∗ = 1250000 mm/s) of the damage evolution. As expected, the smallest propagation velocity yields a slight increase of the ductility. For the higher values of v ∗ , it is clearly observed that a converged force-displacement response is obtained for a finite v ∗ . The converged response is achieved when the fracture area propagation velocity is between 250000 mm/s - 1250000 mm/s. These values are in line with the experimentally evaluated crack propagation speed in steel, e.g the crack propagation velocity for 4340 steel is approximately 1000000 mm/s [32]. 18

7

x 10

4

6 pí Éé W=bñí Éå ë á ç =c ê ~ã ÉW=OUT qç í ~ä =qá ã ÉW=MKMMQRVO

Force [N]

5

4

3

2

1

0 0

v $ = 10000 mm/s v $ = 250000 mm/s v $ = 500000 mm/s v $ = 1000000 mm/s v $ = 1250000 mm/s 2

4

6

Displacement [mm]

8

10

12

Figure 7: Reaction force versus prescribed displacement for tensile loaded plane strain sheet. The response curves correspond to different fracture area propagation velocities v ∗ . The FE-discretization considers the element length le = 0.256 mm and the length-scale parameter is lc = 0.05 mm. Note the identical damage pattern at the final load step obtained for all choices of v ∗ , where blue means α ≈ 0 and red means α ≈ 1.

The damage mode associated with the converged response is shown in Figure 7. Here, the width of the localized damage mode is significantly larger than the lc – and le –parameters. 3.3.2. Influence of fracture energy and the length scale parameter It may be noted that for ”high” values of v ∗ the damage dissipation rate in Eq.(39) is governed by the ratio Gc /lc , corresponding to constant global damage dissipation rate for all choices βGc and βlc where β is a scaling factor. The assumption is that the converged damage mode is the same in the FE application. To see this, consider the plane strain sheet with the constant ratio Gc /lc = 240 N/mm2 (for the motivated values Gc = 12 N/mm and lc = 0.05 mm) with β = 1, 2, 4, 8, 16, 100. Figure 8 shows the FE element response for the different cases. As suggested, the very similar force vs prescribed displacement responses and localized damage modes (as in Figure 7) for all values of β are obtained. It is thus concluded that for a fixed ratio Gc /lc there is a range of Gc and lc -values where the global response and localized damage modes are the same. 3.3.3. Influence of varying element size le and mesh convergence Let us consider next the mesh convergence properties of the proposed damage enhanced JC-model for uniform mesh refinements. The chosen element sizes le of the sheet are shown in Figure 9, corresponding to a total model where the crudest mesh contains ca NEL=625 and the finest mesh yields ca NEL=65000. Moreover, based on the findings of the previous 19

7

x 10

4

6

Force [N]

5

4

3

2

1

0 0

Gc Gc Gc Gc Gc Gc

= 12 N/mm, = 24 N/mm, = 48 N/mm, = 96 N/mm, = 192 N/mm, = 1200 N/mm, 2

lc lc lc lc lc lc 4

= 0.05 = 0.10 = 0.20 = 0.40 = 0.80 = 5.00

mm (REF) mm mm mm mm mm 6

Displacement [mm]

8

10

12

Figure 8: Reaction force versus prescribed displacement for plane strain sheet for a fixed ratio Gc /lc = 240 N/mm2 . The FE-discretization considers le = 0.256 mm and the fracture area propagation velocity was set to v ∗ = 1 250 000 mm/s.

sub-section, the fracture area velocity is v ∗ =1 250 000 mm/s. From the global load vs prescribed displacement response of the plate shown in Figure 9, we find that mesh convergence is obtained as the mesh is refined in the range 0.196 mm ≤ le ≤ 2 mm, corresponding to the damage mode of finite width in Figure 7. The observed mesh convergence is also in line with the results in [25], where mesh convergent behaviour is obtained for brittle failure using the same fracture area evolution model. For all the considered FE-discretizations, damage initiates in the middle of the sheet, where the imperfection is located, and evolve in the shear mode through the sheet as illustrated in Figure 7. 3.4. Tensile loaded hollow plane strain sheet In this example the tensile loaded hollow sheet shown in Figure 10 is used to study the significance of strain rate and temperature dependencies. The same geometry and boundary conditions were considered in Miehe et al. [18]. A fixed non-structured mesh with bilinear plane strain elements is used for the analyses. The element size is le ≈ 0.025 mm corresponding to NEL=5319. In order to integrate the FE-nodal forces, reduced 1 × 1integration combined with hourglass control is adopted. For comparison, we also consider the significance of the void nucleation term α˙ 2 /v ∗ in the damage evolution law related to the corresponding static damage evolution rule, when α = αs is used in Eq. (42). The same pearlitic steel parameters as in the previous examples are considered, cf. Table 1. Like in the other examples, the internal length is lc = 0.05 mm. Three cases for the continuum model of the effective material is considered: 20

7

x 10

4

6

Force [N]

5

4

3

2

1

0 0

le le le le le

= 2.000 = 1.000 = 0.435 = 0.256 = 0.196

mm mm mm mm mm 5

Displacement [mm]

10

15

Figure 9: Reaction force versus prescribed displacement for plane strain sheet in tension when v ∗ = 1 250 000 mm/s. The different curves represent uniform mesh refinement with the element size le .

1. a plasticity model with pure isotropic hardening at constant temperature θ = θ0 , named ”plasticity”. We thus set C = 0 corresponding to the rate independent case. 2. a visco-plastic model at constant temperature, named ”visco-plasticity”. 3. the model contains all dependencies, named ”thermo-visco-plasticity”. The results are shown in Figure 11 for the reaction force versus prescribed displacement responses. Here, the solid black curve represents the ”Plasticity” case, corresponding to the response at very slow loading rate. The dotted black curve shows the corresponding response using the static damage evolution α = αs . Similar response curves of the sheet are obtained for the Visco-Plasticity and the Thermo-Visco-Plasticity cases; however, with a much higher peak load. From Figure 11 it is observed that the strain rate has a significant influence on the response and the thermal effect is quite small in this example. It is also concluded that the proposed damage evolution law (including the nucleation term α˙ 2 /v ∗ ) yields a much more stable response as compared to the conventional ”static” damage law when α = αs [α]. In order to demonstrate the ductile failure process of the sheet, we consider a sequence of deformed configurations for the thermo-visco-plastic model during the loading in Figure 12. The first configuration shows the sheet at onset of damage (as dictated by the JC-failure criterion). The contours represent the effective plastic strain k, showing that onset occurs when k ≈ 20%. Figure 12b shows the sheet in the post peak failure region, where the ductile failure is presented by removing the elements that have reached α ≈ 1. The last configuration in Figure 12c shows the deformed configuration at the complete failure state, arrived at in the final load step. Note the residual strain state after the ductile failure.

21

v, R

r

4r

4r Figure 10: Tensile loaded hollow plane strain sheet with r=0.5mm. The base is completely clamped and the top is clamped horisontally and subjected to a prescribed velocity v¯ = 2500 mm/s vertically. 1800

Plasticity Plasticity ( α = α s ) Visco-Plasticity Visco-Plasticity ( α = α s ) Thermo-Visco-Plasticity Thermo-Visco-Plasticity ( α = α s )

1600 1400

Force [N]

1200 1000 800 600 400 200 0 0

0.05

0.1

0.15

0.2

Displacement 2 [mm]

0.25

0.3

0.35

Figure 11: Reaction force versus prescribed displacement curves of the hollow plane strain sheet in Figure 10. The curves represent different assumptions on the continuum model and the adopted damage law. Plasticity/visco-plasticity means rate independent/dependent behavior in the JC-model at isothermal conditions with θ0 = θtrans =20◦ C. For comparison, the corresponding response curves using the static damage evolution rule (α = αs ) are included, dotted lines.

3.5. Tensile split Hopkinson test of preheated cylindrical specimen In this sub-section the capability of the ductile damage failure model to predict results from the tensile split Hopkinson test of [26] is investigated. In the split Hopkinson test cylin22

0.82 0.75 0.68 0.62 0.55 0.48 0.41 0.34 0.27 0.21 0.14 0.07 0.00

(a) At onset of damage.

(b) In post peak regime.

(c) At failure state

Figure 12: Modes of ductile failure at different stages in the deformation process of the hollow sheet. The contours represent the effective plastic strain −k. Elements that have reached failure α ≈ 1 are removed to display the ductile failure process for the thermo-visco-plastic model.

drical specimens with different notch geometries were subjected to tensile loading at high strain rates at the initial temperatures 100◦ C, 300◦ C and 500◦ C. In the current validation the specimen shown in Figure 13 is considered. The material of the specimen is the high 14

v

axisymmetry 3 15

5

5

15 [mm]

10

Figure 13: Geometry, dimensions and boundary conditions of tensile loaded split-Hopkinson cylindrical specimen from [26]. The prescribed velocity is v¯ = 3250 mm/s at the right end.

strength and ductility Weldox steel 460 E, whose JC-constitutive and failure parameters are given in Table 2. All the parameters in Table 2 are from [26], apart from the strain rate parameters C and d4 , which are modified as compared to those in [26]. In order to fit the model data in [26] to the current work, a least square fit of the strain-rate dependence term was conducted. As a result, the parameters C and d4 in Table 2 were obtained. Furthermore, the fracture energy release rate for Weldox 460 E is Gc = 24 N/mm using the fracture toughness based on Domex 450 in [33], which has similar material characteristics as Weldox 460 E. In order to fit the model with the experiment, a representative value of the length scale is lc = 0.125 mm. In this way the damage mode can be captured with a reasonable mesh size. The specimen is considered axisymmetric, analyzed using 4-node elements with reduced 1-point integration and hourglass control. Following [26] and [34], a prescribed elongation 23

9

rate at the right end of the specimen is applied, corresponding to an overall strain rate between 450-550 s-1 . The progression of the damage variable is followed until the degradation f [α] = 0.05 is reached, equivalent to α ≈ 0.78 in the integration point of the elements. Then, the element is considered as fully damaged and it is decoupled from the FE-analysis (named ”element deletion”), by neglecting the element contribution to the FE-nodal force vector. Table 2: Johnson-Cook material and failure parameters for the Weldox steel of the preheated specimen from [26].

A [MPa] B [MPa] 490 383 d1 0.0705

d2 1.732

E [GPa] 200

ν 0.33

C 0.012843

n 0.45

m 0.94

θmelt [◦ C] 1527

d3 0.54

d4 -0.012912

d5 0

˙0 [s−1 ] 0.0005

ρ [kg/m3 ] cv [J/Kg◦ C] θtrans [◦ C] αθ [mm/mm◦ C] 7850 452 20 11*10−6

3.5.1. Model validation In order to validate the model, the influence of the choice of fracture area progression velocity v ∗ on the force-elongation response is considered first. This behavior is shown in Figure 14 at an initial temperature of 300 ◦ C in the specimen. It is evident from Figure 14 that convergence is rapidly obtained. We also note that no influence on the fracture area propagation velocity is noticed if the initial temperature is decreased to 100 ◦ C or increased to 500 ◦ C. Therefore, v ∗ =250 000 mm/s is used in the simulations. Figure 15 shows the simulated and experimental reaction force vs prescribed displacement relations for the specimen. A fairly good agreement is obtained for the initial temperatures 100◦ C and 300◦ C. Both the experiment and the model exhibit a temperature induced softening response of the specimen. This softening effect is due geometry changes, possibly combined with thermal softening in the adiabatic stress-strain response. The model prediction for the initial temperature of 500◦ C does not corroborate as nicely as for the lower initial temperatures, which is most likely due to the almost linear decrease in yield strength for increased temperature as discussed in [26]. We also investigate the effect of neglecting elastic damage driving energy in Eq.(37), where the elastic contribution A is generally included in AT from the start t > 0. If only the inelastic contribution B f is considered in AT , damage evolution initiates when the onset criterion is fulfilled (t ≥ tf ). These effects are also included in Figure 15. It is clearly seen that the effect of the elastic contribution is not so substantial in this example. In fact, the difference between the representations (AT = A + B f or AT = B f ) is reduced with increased initial temperature.

24

4000 3500 3000

Force [N]

2500 2000 1500 1000 500 0 0

v$ v$ v$ v$ v$ v$

= 10000 mm/s = 250000 mm/s = 500000 mm/s = 1000000 mm/s = 1250000 mm/s = 1500000 mm/s 0.5

1

1.5

Displacement [mm]

2

2.5

3

Figure 14: Reaction force versus prescribed displacement for tensile split-Hopkinson bar in Figure 13. The convergence for increasing propagation velocities v ∗ at an initial temperature of 300 ◦ C is displayed. The FE-discretization is set to le =0.075 mm and the internal length is lc = 0.125 mm. Only the inelastic contribution B f to the damage driving energy is considered. 4500 4000 3500

Force [N]

3000 2500 2000 1500 1000 500 0 0

θ0 θ0 θ0 θ0 θ0 θ0 θ0 θ0 θ0

= 100 °C = 100 °C = 100 °C = 300 °C = 300 °C = 300 °C = 500 °C = 500 °C = 500 °C

(AT = A + B f ) (A T = B f ) ( EXP ) (AT = A + B f ) (AT = B f ) ( EXP ) (AT = A + B f ) (AT = B f ) ( EXP

0.5

1

1.5

Displacement [mm]

2

2.5

3

Figure 15: The influence of damage driving energy (AT ) for the smooth geometry having the initial specimen temperatures 100 ◦ C, 300 ◦ C and 500 ◦ C. The numerical results are compared to experimental data obtained from [26].

3.5.2. Mesh sensitivity of the finite element solution In order to investigate the mesh sensitivity of the finite element solution, a set of FE- discretizations with the element length le ranging from 0.125 mm - 0.025 mm are considered. Like in the previous analyses, the length scale parameter is lc =0.125 mm. The resulting mesh sensitivity of the ductile damage process is shown in terms of the global reaction force vs prescribed displacement relations in Figure 16 when the initial temperature is 300◦ C. 25

In Figure 16 mesh convergence is observed after the onset of elastic-plastic damage driving energy in the structural softening range until α ≈ 0.5 is approached. When α is further increased, a failure stage follows where the post peak force vs prescribed displacement responses converge towards steeper softening slopes as le decreases. For the two finest meshes, with le ≤ 0.05 mm, branching of the damage mode into two paths is observed. For the smallest element length le = 0.025 mm, one of these branches is dominant in the very end of the deformation process, as demonstrated by the non-symmetric failure mode of the specimen shown in Figure 17. In the final stage of the failure analysis when using the second finest mesh, with le = 0.05 mm, a symmetric failure mode is obtained, which explains the slight difference in post peak response for the two finest meshes in Figure 16. In addition, when element deletion is applied for f [α] = 0.05 in the fracture zone, a vibration of the specimen is observed, as manifested by the oscillating response in Figure 16 for le = 0.05 mm. For the coarsest mesh with le = 0.125 mm a uniform damage zone is observed as shown in Figure 17, where the red color indicates a fully damaged zone (f [α]= 0.05). Hence, it is concluded that the present formulation has good mesh convergence behavior up to high values of the damage (α ≈ 0.5 or f [α] ≈ 0.25) in this example. In the very final failure state, the mesh convergence of the FE solution is less conclusive; partly because of the branching of the damage mode for the two finest meshes, and partly because of element deletion, which is known to destroy the mesh convergence property as discussed in [7]. 4000

le le le le le

3500 3000

Force [N]

2500

= 0.125 = 0.100 = 0.075 = 0.050 = 0.025

mm mm mm mm mm

Damage onset

2000 1500

Fracture

1000 500 0 0

0.5

1

1.5

Displacement [mm]

2

2.5

3

Figure 16: Reaction force vs prescribed displacement for tensile split-Hopkinson bar with the initial temperature 300 ◦ C. The figure shows the mesh sensitivity for increasing FE-resolution, where the element sizes range between le = 0.125 mm to le = 0.025 mm. The length scale parameter was set to lc = 0.125 mm.

4. Concluding remarks In the present paper we have presented a thermodynamically motivated damage enhanced model for ductile fracture processes. To obtain and to explain the damage/plasticity coupling, the model development was placed in context of continuum thermodynamics. A 26

le = 0 .125mm

le = 0 .025mm

Figure 17: Details of the ductile split-Hopkinson tensile bar fracture indicated in Figure 16 for the finest FE discretization le = 0.025 mm. The initial temperature is 300◦ C. The close-ups show the captured damage pattern for the coarsest and finest FE discretizations. Here, ”red” means fully developed damage.

novel feature of the ductile damage model is the (input) damage driving dissipation rate, covering elastic and inelastic contributions separated by the damage threshold. The final danage model is arrived at from the balance between input and produced dissipation rate, in terms of the diffuse fracture area parameters v ∗ and lc without the gradient damage term. In this fashion, a computationally efficient model amenable for implementation in commercial platforms like ABAQUS/ExplicitTM , cf. also [25], is obtained. As the prototype for the effective material the visco-plastic JC-constitutive model was considered, whereas the JC failure model was used to determine the damage threshold. There are several enhancements of the JC-continuum/failure models available in the literature. However, in the present work the original JC-models [22], [5] were considered to suffice. Also values for material parameters were available for these models. One may note, though, that the proposed framework opens up for (more refined) alternative modeling of the effective material and the damage threshold depending on the application. The adopted damage threshold for the elastic-plastic damage driving energy makes the formulation flexible for ductile failure modeling. For example, if the threshold time tf → ∞, a pure elastic damage driving process is obtained, whereas tf → 0 yields a completely elastic-plastic driven failure process. Connected to the first example, it was shown that the approach is reasonable in the sense that the total dissipation generated (including elastic and inelastic potions) enters into the damage driving process. For thresholds 0 < tf < ∞, the strategy is to limit the influence of elastic damage driving portion, whereby the internal length lc should be kept low. To allow for more flexibility in the lc -value, it was proposed to simply neglect elastic damage driving portion A. However, this choice makes the model less thermodynamically consistent, whereby Gc might not properly reflect the fracture energy. As to the practical choice of the internal length scale parameter lc we make the following general conclusions: 1. it is chosen so that the observed ductility and localization mode are properly captured 27

by the model in calibration situations. Note that for a given fracture energy release rate Gc , a decreased lc -value yields a more ductile behavior. 2. it cannot be too small since it affects the width of the localization zone, cf. also [25], which in turn must be resolved by an reasonable element size. 3. it has a value that is not too large in order to reduce the influence of elastic damage driving energy. From our experience satisfactory convergence properties are obtained for ductile as well as brittle failure processes, cf. [25]. A common problem in mesh convergence studies of dynamic ductile fracture is localized branches of the damage mode developing with enhanced FE-resolution. In the plane strain sheet example with a central imperfection, where no branching of the damage pattern could be observed, we obtained mesh convergent behavior towards a localized shear mode. In the hollow plane stain sheet example, it was also observed that the nucleation term α˙ 2 /v ∗ indeed stabilizes the response as compared to standard static damage evolution. From the numerical examples, the fracture area propagation velocity v ∗ has a significant effect on the response. An interesting model feature is that when the propagation velocity is increased, convergence in force–displacement response and damage mode is obtained for a finite v ∗ ; this conclusion was also made in [25] for brittle damage processes. In some applications, we have to increase v ∗ so that it approaches the crack propagation speed to obtain this convergence. We also considered the behavior for fixed ratios Gc /lc (with βGc and βlc ). Here, it was concluded that the force-displacement response is fairly constant as long as the same damage mode is captured for the considered range of β-values. Hence, the internal length lc does not determine the width of the localization zone alone. The width also depends on the fracture energy release and boundary conditions. In the validation example our model was compared to the split-Hopkinson tension bar experiments in [26]. The model was able to well predict the response up to the initial temperature 300◦ C. For the higher temperature 500◦ C the agreement with experiments was not as good, which may be due to a less known temperature dependence of the yield stress. It was also shown that mesh convergence was obtained up to α ≈ 0.5. For higher damage values, an enhanced mesh refinement resolved a localized branching of the damage mode, associated with less conclusive mesh convergent behavior in the very final part of the ductile failure analysis. 5. Acknowledgement This research was financially supported by The Ekmanska Foundation (Ekmanska Stiftelsen) and the Swedish national research program FFI (Strategic Vehicle Research and Innovation). The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Chalmers Centre for Computational Science and Engineering (C3SE).

28

References [1] A. Tay, M. Stevenson, G. de Vahl Davis, Using the finite element method to determine temperature distributions in orthogonal machining, Proceedings of the Institution of Mechanical Engineers 188 (1) (1974) 627–638. [2] R. Ambati, H. Yuan, FEM mesh-dependence in cutting process simulations, The International Journal of Advanced Manufacturing Technology 53 (2011) 313–323. [3] P. Arrazola, T. Ozel, D. Umbrello, M. Davies, I. Jawahir, Recent advances in modelling of metal machining processes, CIRP Annals - Manufacturing Technology 62 (2013) 695–718. [4] G. Ljustina, R. Larsson, M. Fagerstr¨om, A FE based machining simulation methodology accounting for cast iron microstructure, Finite Elements in Analysis and Design 80 (2014) 1–10. [5] G. Johnson, W. Cook, Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures, Engineering Fracture Mechanics 21 (1985) 31–48. [6] T. Mabrouki, J. Rigal, A contribution to a qualitative understanding of thermo-mechanical effects during chip formation in hard turning, Journal of Materials Processing Technology 176 (2006) 214–221. [7] R. Larsson, S. Razanica, B. Josefson, Mesh objective continuum damage models for ductile fracture, International Journal for Numerical Methods in Engineering 106 (2015) 840–860. [8] A. L. Gurson, Continuum theory of ductile rupture by void nucleation and growth: Part I-yield criteria and flow rules for porous ductile media, Journal of Engineering Materials and Technology 99 (1977) 2–15. [9] A. Needleman, V. Tvergaard, An analysis of ductile rupture in notched bars, Journal of the Mechanics and Physics of Solids 32 (1984) 461–490. [10] T. Pardoen, J. Hutchinson, An extended model for void growth and coalescence, Journal of the Mechanics and Physics of Solids 48 (2000) 2467–2512. [11] V. Tvergaard, A. Needleman, Nonlocal effects on localizaion in a void-sheet, International Journal of Solids and Structures 34 (18) (1997) 2221–2238. [12] F. Reusch, B. Svendsen, D. Klingbeil, Local and non-local Gurson-based ductile damage and failure modelling at large deformation, European Journal of Mechanics A/Solids 22 (2003) 779–792. [13] O. Abiri, L. Lindgren, Non-local damage models in manufacturing simulations, European Journal of Mechanics - A/Solids 49 (2015) 548–560. [14] 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) (2014) 383–405. [15] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations, International Journal for Numerical Methods in Engineering 83 (2010) 1273–1311. [16] C. Miehe, L. M. Schanzel, H. Ulmer, Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids, Computer Methods in Applied Mechanics and Engineering 294 (2015) 449–485. [17] M. Ambati, T. Gerasimov, L. De Lorenzis, Phase-field modeling of ductile fracture, Computational Mechanics 55 (5) (2015) 1017–1040. [18] C. Miehe, M. Hofacker, L. Schanzel, F. Aldakheel, Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo–elasticplastic solids, Computer Methods in Applied Mechanics and Engineering 294 (2015) 486–522. [19] V. Dias da Silva, A simple model for viscous regularization of elasto-plastic constitutive laws with softening, Communications in Numerical Methods in Engineering 29 (2004) 547–568. [20] M. Niazi, H. Wisselink, T. Meinders, Viscoplastic regularization of local damage models: revisited, Computational Mechanics 51 (2) (2013) 203–216. [21] T. Borvik, O. S. Hopperstad, T. Berstad, M. Langseth, A computational model of viscoplasticity and ductile damage for impact and penetration, European Journal of Mechanics - A/Solids 20 (2001) 685– 712. [22] G. Johnson, W. Cook, A constitutive model and data for metals subjected to large strains, high

29

[23] [24] [25] [26]

[27] [28] [29] [30] [31] [32] [33] [34]

strain rates and high temperatures, Proceedings 7th International Symposium on Ballistics, Hague, Netherlands, April 19–21, 1983. J. Lemaitre, A short course in damage mechanics, Springer-Verlag, 1992. G. Ljustina, M. Fagerstr¨om, R. Larsson, Hypo- and hyper-inelasticity applied to modeling of compact graphite iron machining simulations, European Journal of Mechanics - A/Solids 37 (2012) 57–68. R. Larsson, R. Gutkin, M. Rouhi, Damage growth and strain localization in compressive loaded fiber reinforced composites, Mechanics of Materials 127 (2018) 77–90. T. Borvik, O. S. Hopperstad, S. Dey, E. V. Pizzinato, M. Langseth, C. Albertini, Strength and ductility of weldox 460 e steel at high strain rates, elevated temperatures and various stress triaxialities, Engineering Fracture Mechanics 72 (2005) 1071–1087. C. Miehe, J. Simo, Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation, Computer Methods in Applied Mechanics and Engineering 98 (1992) 41–104. M. Ristinmaa, M. Wallin, N. Ottosen, Thermodynamic format and heat generation of isotropic hardening plasticity, Acta Mechanica 194 (2007) 103–121. M. G. Cockcroft, D. J. Latham, Ductility and the workability of metals, Journal of The Institute of Metals 96 (1968) 33–39. G. Ljustina, M. Fagerstr¨om, R. Larsson, Rate sensitive continuum damage models and mesh dependence in finite element analyses, The Scientific World Journal (2014) doi:10.1155/2014/260571. H. Roy, N. Parida, S. Sivaprasad, S. Tarafder, K. K. Ray, Acoustic emissions during fracture toughness tests of steels exhibiting varying ductility, Materials Science and Engineering A 486 (2008) 562–571. T. A. Zehnder, J. A. Rosakis, Dynamic fracture initiation and propagation in 4340 steel under impact loading, International Journal of Fracture 43 (1990) 271–285. B. Sundstr¨ om, Handbok och formelsamling i h˚ allfasthetsl¨ara, Instant Book AB, Stockholm, 2007. T. Borvik, O. S. Hopperstad, T. Berstad, On the influence of stress triaxiality and strain rate on the behaviour of structural steel. Part II. Numerical study, European Journal of Mechanics - A/Solids 22 (2003) 15–32.

30