A method for coupling atoms to continuum mechanics for capturing dynamic crack propagation

. Conventionally, dynamic crack propagation is modelled using fracture mechanics (either linear elastic, or with an extension to confined plasticity). Herein, we propose a different view, based on a coupling between an atomic description at the crack tip and a classical continuum description away from it. The paper presents the theoretical background and some first numerical results. RÉSUMÉ. La modélisation de la propagation de fissure repose principalement sur une bonne utilisation de la mécanique de la rupture (dans le cadre élastique linéaire ou bien dans des extensions à la plasticité confinée). Nous proposons ici une approche différente, tant sur la modélisation numérique que physique, qui consiste à coupler une vision atomique en pointe de fissure à une description classique continue. Cet article présente le cadre théorique du couplage ainsi que quelques résultats numériques en dimension 1.


Introduction
In most approaches to date, fracture is modelled using a continuum approach, either using linear-elastic fracture mechanics, or with an extension to take care of confined inelastic yielding, or crack bridging effects, such as in cohesive-zone models.In this contribution we shall describe a novel approach to fracture.The basicideais that the crack tip is described by an atomistic model while the surrounding will be described as a continuum, discretized via a finite element method that exploits the partition-of-unity property of finite element shape functions.Indeed, byu s i n ge nriched interpolation functions to capture cracks that may be incompatible with the underlying mesh structure no remeshing is needed.At the same time, there is no need to use asymptotic enriched functions at the continuum level, since the crack tip is described by an atomistic model.At the atomistic scale, either a classical molecular dynamics description or quantum mechanics can be employed.The latter approach relies on a very detailed description: The Schrödinger equation is solved as a function of the electronic configuration (Marx et al., 2000).The major inconveniency of such a description is its cost.A more pragmatic approach, that is used here (Zhou et al., 1997;Sutmann, 2002), employs a classical molecular approach with interatomic interactions being captured by a potential.
The coupling of the two descriptions is the major challenge of this contribution.Indeed, away from the crack tip we have a classical continuum, but in its vicinity a discrete atomistic model is used.A crucial difficulty of this approach is that characteristic parameters in space and time are much smaller at the atomic level than those at the continuum level, which are used in the finite element simulation.
Many coupling methods have been developed (Liu et al., 2006;Xiao et al., 2004).Herein, we propose a weak coupling between the two models, namely an energy coupling.Its major advantage is that energy has the same physical meaning in both domains.To be more precise, the continuum equations are cast in a weak form andare coupled to the atomic description via a partition of the energy.

Changing scales in the crack propagation problem
Dynamick crack propagation can be described as follows: A discontinuity, the crack Γ, is included in a structure described as the closure of the connex open set Ω.The displacements are fixed on a boundary ∂ 1 Ω and external loads are given on ∂ 2 Ω.Finally, a volumic forces field, g d , is applied in Ω (Figure 1).
We consider a subdomain Ω m of Ω that includes the crack tip.The goal is to link an atomic description in the subdomain Ω m to a continuum approach in Ω\Ω m ,i .e .outside the zone of the crack tip.In order to separate the two zones, Ω m can be defined spatially as the crak tip domain to which the plastic zone is confined, so that the macro subdomain outside this plastic zone can be captured by an elastic stress-strain relation.

Formulations
Denoting by u and d the displacements in the continuum and in the discrete zone, respectively, we have the initial-value problem in the continuum subdomain: For x ∈ Ω M (t) and t ∈ [0; T ] , knowing u (x, 0) , u (x, 0) , find (u, σ) ∈ U ad × S ad such that: ρü = div σ + g d [1] with: and in the discrete subdomain Ω m , where we build a discrete grid of N a atoms, the second initial-value problem reads as follows: with: The interatomic forces are given by differentation of a potential.A good description relies on the proper choice of this potential.To test the coupling method a classical potential, such as that of Lennard-Jones or of Morse, suffices: The parameters ε, σ, D and a are material properties and r is the interatomic distance, with r e that at equilibrium.We note that other potentials have been developed, allowing for a better description of the metallic bond, for example -EAM potential (Daw, 1989).The atomic problem relies on the resolution of non-linear equations with non-convex potentials.The study of existence and uniqueness properties for this kind of problem is not the aim of this contribution, see e.g., (Lions et al., 1965) for more information.
In order to get an energetic framework and allow for a discretization, the continuum problem is rewritten in a weak format, as follows: or, written in a more concise way: The atomic problem becomes:

Coupling models
In order to achieve an efficient coupling between the two problems introduced before and to specify the coupling conditions, we assume a weak coupling on the common zone.We consider the mechanical energy as a primordial quantity.Dualizing formulations and writing them in a weak format allows us to obtain a global description that preserves the descriptive properties of each model, and focuses on the quantity of interest, the energy, which must not depend on the model.
In the common, or handshaking, zone Ω c = Ω M ∩ Ω m , a velocity coupling is adopted using a weak format.Moreover, the energy is distributed between the two models inside the coupling zone in a partition of the unity sense (Dhia et al., 2005).More precisely, on the whole domain Ω we split the energy between the two models in the following sense: and we subsequently obtain the weak formulation for the distribution of the energy:

. Partition of unity for the energy distribution
As stated before, the displacement and velocity fields in the domains Ω M and Ω m have a different nature.One field is a continuous in its definition space, while the other has a discrete character and is only defined at the geometrical points that correspond to the atoms.In order to construct a velocity coupling, we assume an equality between the velocities as we go from one model to the other.This condition has to be written in a weak form, which means in a "global", or "integral" manner.
Accordingly, a new space is constructed, denoted by M and called the "mediator space", on which we will project the fields u and ḋ in order to compare them.The nature of M is determined by the discrete property of the atomic fields.Indeed, no extrapolation outside the atomic positions is possible if we wish to keep a physical sense at the fine scale.Then, M has to be a sub-space of the physical atomic space.
More precisely, through the operator Π we can project the velocities onto a discrete subset Ω c of the atomic positions included in Ω c .Considering that M has been constructed as a Hilbert space, we introduce a scalar product c from M × M into R, and we write the velocity coupling as: This formulation allows us to finally write the global equations coupled with Lagrange multipliers:

Discretization
We now introduce two main discretizations.First, the macro-problem in Ω M has to be discretized with a finite element interpolation, and subsequently, a time discretization is needed to solve the dynamic global problem.Moreover, we will use a Heaviside enrichment at the crack in order to avoid remeshing by exploiting the partition-of-unity property of finite element shape functions (Moës et al., 1999;Remmers et al., 2003).

Spatial discretization for the continuum problem
The weak formulation allows us to adopt a finite element form for the macroproblem.With the shape functions N i and nodal unknown vectors u i , using the Heaviside step function to take into account the discontinuity, one obtains: With the latter notation, the bilinear form a M (., .)and the linear form, l M (.), become: Moreover, with the classical scalar product, the coupling term in the continuum becomes: which leads to the typical matrix formulation: Remark 1: The added term F L M is a fictitious force due to the coupling with Lagrange multipliers.This force has non-zero value only in the coupling zone Ω c .
Remark 2: The matrices M, K and the vector F M contain information about the energy distribution, i.e. in the domain Ω c , the repartition function α is used to build their elementary terms.
Remark 3: The vector Λ stands for the Lagrange unknowns and its size is equal to the Ω c subset cardinal times the dimension of the considered space.

Time integration scheme
In order to solve the coupled system we use a standard central difference scheme.This scheme is widely used in classical molecular dynamics and is also knowna s Verlet algorithm.Yet, as it is explicit, its stability relies on an appropriate choice of the time step size.For the continuum problem other time schemes could have been used, but for simplicity the central difference scheme has also been adopted for this subdomain.
The global coupling equations are: Similar to the continuum model, the equation for the atomic scale includes the repartition function β, and a fictitious coupling force f L m is introduced.The matrix m is diagonal with elementary term βm i , and the coupling force is: The system can be rewritten with the three main vectors (U, d, Λ): The time scheme relies on discretization with a time step Δt and has four stages: -given the quantities at step n, compute the displacements U n+1 and d n+1 , -compute the accelerations Ün+1 and dn+1 , neglecting the Lagrange forces, -compute the predictive velocities U * n+1 and ḋ * n+1 , -adjust these velocities to give the final velocities Un+1 and ḋn+1 by taking into account the coupling terms and Lagrange multipliers Λ n+1 .
The three first steps are simple and do not need further explanation.The lasts t e p computes the coupling terms and enforces condition (9).
From ( 9) and ( 18) with the matrix M standing for M lumped , the coupling condition becomes: where Λn+1 = 1 2 (Λ n+1 + Λ n ).The new Lagrange multipliers Λn+1 are subsequently computed by solving: [20] with: and b n+1 stands for the weak coupling condition on the predictive velocities.Thus, this term is a measure of the error compared to the solution that satisfies thes y s t e m (18).

A multiscale time decomposition
As we have very different orders of magnitude in the models, it is useful to construct a multiscale time integration scheme.Indeed, at the atomic scale, we need a very fine time scale to achieve sufficient accuracy and to satisfy the critical time step condition.The concept is simple: Let Δt be the fine time step, i.e. for the atomistic problem, and Δ T be the coarse time step, such that Δ T = kΔt with k ∈ N * .The entire atomistic resolution is done at the fine scale, but the macro-problem is solved only at coarse time steps.For consistency between the models, the coupling condition has been enforced at each fine time step, via an interpolation of the macro acceleration.Indeed, there is no need to compute all the macro quantities on the fine grid, but we must approximate the acceleration in order to update the velocities for the coupling condition.

One-dimensional results
As an example we consider a longitudinal bar, which is described with finite elements and with molecular dynamics, with a coupling region in between.The bar is submitted to a traction wave, the initial configuration being displaced on the leftmost ten elements.The right-hand side is fixed.The whole domain is 88.1563nm long and we put 176 atoms in the molecular domain.The interatomic distance is r e = 0.3253nm, and the finite elements size is h = 4r e .In the molecular model, we use a Lennard Jones potential, with ε = 43.306zJand a mass m = 0.00448yg.T h e s e parameters come from FCC Al properties, and the proper material propertiesforthe finite element model are derived from the atomic properties.In Figure 3, we see the displacements in the bar as the wave propagates.The zoom 3(e) shows how the atoms fit the travelling wave in the coupling zone.
We now focus on the mechanical energy transfer when the wave passes through the coupling zone from the finite elements to the molecular domain.In the first case (Figure 4(a)), the coupling length L c covers 4 elements.The wave passes from one domain to the other, but there is a non-negligible energy loss.When L c covers 10 elements, i.e. the second case (Figure 4(b)), the mechanical energy transfer has improved.A good transfer is reached when γ = L c L i → 1, and then we have negligible energy losses (< 1%).On the contrary, when the coupling length is too small compared to the wave length, there is some reflection in the finite element domain.In Figure 4, we have plotted the energy transfer ratio -when the wave passes from the finite element to the molecular domain -as a function of the number of atoms involved in the coupling domain, for different multscale ratios (from h = r e to h = 10r e ).We observe that it does not depend on the multiscale aspect.Only the number of atoms has an influence.To achieve a good energy balance, a sufficient number of atoms in the coupling zone is needed: With approximately 20 atoms we obtain less than 2% of energy loss.

Conclusion
We have given the theoretical concepts of a coupling method between two models that are physically different.The formulation has been written in a weak format in order to preserve the accuracy of each model and the mechanical energy as the fundamental quantity.We have studied this method in a one-dimensional example and have observed the energy transfer as information passes from one domain to the other.The coupling length and the number of atoms involved in the transition zone have a major influence on the energy balance.

Figure 3 .
Figure 3. Wave propagation in a one-dimensional beam

Figure 4 .
Figure 4. Mechanical energy transfer from FE to MD

Figure 5 .
Figure 5. Energy transfer ratio depending on the multiscale ratio