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

Pascal Aubertin, René de Borst, Julien Réthoré

► To cite this version:

Pascal Aubertin, René de Borst, Julien Réthoré. A method for coupling atoms to continuum mechanics for capturing dynamic crack propagation. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2012, 17 (5-7), pp.651-662. 10.3166/remn.17.651-662. hal-00381486

HAL Id: hal-00381486

https://hal.science/hal-00381486v1

Submitted on 27 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A method for coupling atoms to continuum mechanics for capturing dynamic crack propagation

Pascal Aubertin — René de Borst — Julien Réthoré

Laboratoire de Mécanique des Contacts et des Structures

INSA de Lyon - CNRS UMR 5259

18-20, rue des Sciences, F-69621 Villeurbanne cedex

pascal.aubertin@insa-lyon.fr ; r.d.borst@tue.nl ; julien.rethore@insa-lyon.fr

**ABSTRACT.** 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.

**KEYWORDS:** molecular dynamics, finite elements, energy coupling, crack propagation.

**MOTS-CLÉS :** dynamique moléculaire, éléments finis, couplage faible, propagation de fissure.
# 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 basic idea is 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, by using enriched 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 and are coupled to the atomic description via a partition of the energy.

# 2. Problem statement

## 2.1. *Changing scales in the crack propagation problem*

Dynamick crack propagation can be described as follows: A discontinuity, the crack $\Gamma$, is included in a structure described as the closure of the connex open set $\Omega$. The displacements are fixed on a boundary $\partial_1\Omega$ and external loads are given on $\partial_2\Omega$. Finally, a volumic forces field, $g_d$, is applied in $\Omega$ (Figure 1).

We consider a subdomain $\Omega_m$ of $\Omega$ that includes the crack tip. The goal is to link an atomic description in the subdomain $\Omega_m$ to a continuum approach in $\Omega \setminus \Omega_m$, i.e. outside the zone of the crack tip. In order to separate the two zones, $\Omega_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.
<figure><img src="image_2.png" /><figcaption>Figure 1. Global problem</figcaption></figure>

## 2.2. 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:

$$ \text{For } x \in \Omega_M(t) \text{ and } t \in [0; T], \text{ knowing } (u(x, 0), \dot{u}(x, 0)), \\ \text{find } (u, \sigma) \in \mathcal{U}^{ad} \times \mathcal{S}^{ad} \text{ such that:} \tag{1} $$

$$ \rho \ddot{u} = \mathbf{div} \, \sigma + g_d $$

with:

$$ \mathcal{U}^{ad} = \{ u = u(x,t) \in [\mathcal{H}^1(\Omega_M)]^3 ; u = u_d \text{ on } \partial_1 \Omega, \forall t \in [0, T] \} $$

$$ \mathcal{S}^{ad} = \{ \sigma = \mathbb{K}\epsilon(u(x,t)) \in [\mathcal{L}^2(\Omega_M)]^6 ; \sigma.n = F_d \text{ on } \partial_2 \Omega, \forall t \in [0, T] \} $$

and in the discrete subdomain $\Omega_m$, where we build a discrete grid of $N_a$ atoms, the second initial-value problem reads as follows:

$$ \text{For } 1 \le i \le N_a(t) \text{ and } t \in [0; T], \text{ knowing } (d(0), \dot{d}(0)), \\ \text{find } (d, f) \in \mathcal{D}^{ad} \times \mathcal{F}^{ad} \text{ such that:} \tag{2} $$

$$ m_i \ddot{d}_i = f_i $$
with:

$$ \mathcal{D}^{ad} = \left\{ d = (d_i(t))_{1 \le i \le N_a}, \forall t \in [0, T] \right\} $$

$$ \mathcal{F}^{ad} = \left\{ f = \left( f_i(t) = -\nabla_i W(d(t)) \right)_{1 \le i \le N_a}, \forall t \in [0, T] \right\} $$

The interatomic forces are given by differentiation 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:

$$ W_{L-J} = 4\epsilon \left( \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right), \quad W_{Morse} = D \left(1 - e^{-a(r-r_e)}\right)^2 \quad [3] $$

The parameters $\epsilon$, $\sigma$, $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:

$$ \begin{gathered} \forall v^* \in \dot{\mathcal{U}}^{ad,0}, \text{ knowing } (u(x,0), \dot{u}(x,0)), \\ \text{find } u \in \mathcal{U}^{ad} \text{ such that:} \end{gathered} $$

$$ \int_{\Omega_M} \rho \ddot{u} \cdot v^* d\Omega + \int_{\Omega_M} \varepsilon(u) : \mathbb{K} : \varepsilon(v^*) d\Omega = \int_{\partial_2 \Omega} F_d \cdot v^* dS + \int_{\Omega_M} g_d \cdot v^* d\Omega \quad [4] $$

or, written in a more concise way:

$$ a_M(u, v^*) = l_M(v^*) \quad [5] $$

The atomic problem becomes:

$$ \begin{gathered} \forall w^* \in \dot{\mathcal{D}}^{ad,0}, \text{ knowing } (d(0), \dot{d}(0)), \\ \text{find } d \in \mathcal{D}^{ad} \text{ such that:} \end{gathered} $$

$$ \sum_{i=1}^{N_a} m_i \ddot{d}_i \cdot w_i^* + \sum_{i=1}^{N_a} \nabla_i W(d) \cdot w_i^* = a_m(d, w^*) = 0 \quad [6] $$
## 2.3. 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 $\Omega_c = \Omega_M \cap \Omega_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 $\Omega$ we split the energy between the two models in the following sense:

$$ \left\{ \begin{array}{ll} \alpha(x) = 1 & \text{for } x \in \Omega_M \setminus \Omega_c \\ \beta(x) = 1 & \text{for } x \in \Omega_m \setminus \Omega_c \\ \alpha(x) + \beta(x) = 1 & \text{for } x \in \Omega_c \end{array} \right. \qquad [7] $$

and we subsequently obtain the weak formulation for the distribution of the energy:

$$ \begin{gathered} \forall (v^*, w^*) \in \dot{\mathcal{U}}^{ad,0} \times \dot{\mathcal{D}}^{ad,0}, \text{ knowing } (u(x,0), \dot{u}(x,0), d(0), \dot{d}(0)), \\ \text{find } (u,d) \in \mathcal{U}^{ad} \times \mathcal{D}^{ad} \text{ such that:} \end{gathered} $$

$$ \alpha a_M(u, v^*) + \beta a_m(d, w^*) = \alpha l_M(v^*) \qquad [8] $$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Partition of unity for the energy distribution</figcaption></figure>

As stated before, the displacement and velocity fields in the domains $\Omega_M$ and $\Omega_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 $\mathcal{M}$ and called the “mediator space”, on which we will project the fields $\dot{u}$ and $\dot{d}$ in order to compare them. The
nature of $\mathcal{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, $\mathcal{M}$ has to be a sub-space of the physical atomic space.

More precisely, through the operator $\Pi$ we can project the velocities onto a discrete subset $\overline{\Omega}_c$ of the atomic positions included in $\Omega_c$. Considering that $\mathcal{M}$ has been constructed as a Hilbert space, we introduce a scalar product $c$ from $\mathcal{M} \times \mathcal{M}$ into $\mathbb{R}$, and we write the velocity coupling as:

$$ \forall \mu^* \in \mathcal{M}, \quad c(\mu^*, \Pi\dot{u} - \Pi\dot{d}) = \langle \mu^*, \Pi\dot{u} - \Pi\dot{d} \rangle_{\mathcal{M}} = 0 \quad [9] $$

This formulation allows us to finally write the global equations coupled with Lagrange multipliers:

$$ \forall (v^*, w^*, \mu^*) \in \dot{\mathcal{U}}^{ad,0} \times \dot{\mathcal{D}}^{ad,0} \times \mathcal{M}, \text{ knowing } (u(x,0), \dot{u}(x,0), d(0), \dot{d}(0)), \\ \text{find } (u, d, \lambda) \in \mathcal{U}^{ad} \times \mathcal{D}^{ad} \times \mathcal{M} \text{ such that:} $$

$$ \alpha a_M(u, v^*) + \beta a_m(d, w^*) + c(\lambda, \Pi v^* - \Pi w^*) = \alpha l_M(v^*) \\ c(\mu^*, \Pi \dot{u} - \Pi \dot{d}) = 0 \quad [10] $$

### 3. Discretization

We now introduce two main discretizations. First, the macro-problem in $\Omega_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).

#### 3.1. Spatial discretization for the continuum problem

The weak formulation allows us to adopt a finite element form for the macro-problem. 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:

$$ \forall x \in \Omega_M, \quad u_h(x) = \sum_{i \in \mathcal{N}_i} U_i N_i(x) + \sum_{e \in \mathcal{N}_e} U_e^e N_i(x) \mathcal{H}(x) = N^T U \quad [11] $$

With the latter notation, the bilinear form $a_M(.,.)$ and the linear form, $l_M(.)$, become:

$$ \alpha a_M(u_h, v_h^*) = V^{*T} M \ddot{U} + V^{*T} K U \quad [12] $$

$$ \alpha l_M(v_h^*) = V^{*T} F_M \quad [13] $$
Moreover, with the classical scalar product, the coupling term in the continuum becomes:

$$c(\lambda, \Pi v_h^*) = V^{*T} C_M \Lambda = V^{*T} F_M^L \quad [14]$$

which leads to the typical matrix formulation:

$$M \ddot{U} + KU = F_M - F_M^L = F_M - C_M \Lambda \quad [15]$$

**Remark 1:** The added term $F_M^L$ is a fictitious force due to the coupling with Lagrange multipliers. This force has non-zero value only in the coupling zone $\Omega_c$.

**Remark 2:** The matrices $M$, $K$ and the vector $F_M$ contain information about the energy distribution, i.e. in the domain $\Omega_c$, the repartition function $\alpha$ is used to build their elementary terms.

**Remark 3:** The vector $\Lambda$ stands for the Lagrange unknowns and its size is equal to the $\bar{\Omega}_c$ subset cardinal times the dimension of the considered space.

## 3.2. 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 known as 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:

$$\left\{
\begin{array}{l}
M \ddot{U} + KU = F_M - F_M^L \\
md\ddot{d} = f + f_m^L \\
\forall \mu^* \in \mathcal{M} \quad , \quad c(\mu^*, \Pi \dot{u}_h - \Pi \dot{d}) = 0
\end{array}
\right.
\quad [16]$$

Similar to the continuum model, the equation for the atomic scale includes the repartition function $\beta$, and a fictitious coupling force $f_m^L$ is introduced. The matrix $m$ is diagonal with elementary term $\beta m_i$, and the coupling force is:

$$f_m^L = C_m \Lambda \quad [17]$$

The system can be rewritten with the three main vectors $(U, d, \Lambda)$:

$$\left\{
\begin{array}{l}
M \ddot{U} + KU = F_M - C_M \Lambda \\
md\ddot{d} = f + C_M \Lambda \\
C_M^T \dot{U} = C_M^T \dot{d}
\end{array}
\right.
\quad [18]$$
The time scheme relies on discretization with a time step $\Delta t$ and has four stages:

- given the quantities at step $n$, compute the displacements $U_{n+1}$ and $d_{n+1}$,

- compute the accelerations $\ddot{U}_{n+1}$ and $\ddot{d}_{n+1}$, neglecting the Lagrange forces,

- compute the predictive velocities $\dot{U}_{n+1}^*$ and $\dot{d}_{n+1}^*$,

- adjust these velocities to give the final velocities $\dot{U}_{n+1}$ and $\dot{d}_{n+1}$ by taking into account the coupling terms and Lagrange multipliers $\Lambda_{n+1}$.

The three first steps are simple and do not need further explanation. The last step computes the coupling terms and enforces condition (9).

From (9) and (18) with the matrix $\tilde{M}$ standing for $M_{\text{lumped}}$, the coupling condition becomes:

$$C_M^T (\dot{U}_{n+1}^* - \tilde{M}^{-1} C_M \tilde{\Lambda}_{n+1} \Delta t) = C_m^T (\dot{d}_{n+1}^* + m^{-1} C_m \tilde{\Lambda}_{n+1} \Delta t) \quad [19]$$

where $\tilde{\Lambda}_{n+1} = \frac{1}{2} (\Lambda_{n+1} + \Lambda_n)$. The new Lagrange multipliers $\tilde{\Lambda}_{n+1}$ are subsequently computed by solving:

$$A \tilde{\Lambda}_{n+1} = b_{n+1} \quad [20]$$

with:

$$\left\{ \begin{aligned} A &= (C_M^T \tilde{M}^{-1} C_M + C_m^T m^{-1} C_m) \Delta t \\ b_{n+1} &= C_M^T \dot{U}_{n+1}^* - C_m^T \dot{d}_{n+1}^* \end{aligned} \right. \quad [21]$$

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 the system (18).

### 3.3. 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 $\Delta t$ be the fine time step, i.e. for the atomistic problem, and $\Delta T$ be the coarse time step, such that $\Delta T = k\Delta t$ with $k \in \mathbb{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.
<figure><img src="image_4.png" /><figcaption>(a) Displacements at $t = 0\,\text{ps}$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Displacements at $t = 0.080\,\text{ps}$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Displacements at $t = 0.160\,\text{ps}$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) Displacements at $t = 0.240\,\text{ps}$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) Zoom at $t = 0.080\,\text{ps}$</figcaption></figure>

<caption>**Figure 3.** *Wave propagation in a one-dimensional beam*</caption>
## 4. 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.3253$nm, and the finite elements size is $h = 4r_e$. In the molecular model, we use a Lennard Jones potential, with $\epsilon = 43.306$zJ and a mass $m = 0.00448$yg. These parameters come from FCC Al properties, and the proper material properties for the 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 $\gamma = \frac{L_c}{L_i} \to 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.

<figure><img src="image_9.png" /><figcaption>**Figure 4.** *Mechanical energy transfer from FE to MD*</figcaption></figure>

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 multiscale 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.
<figure><img src="image_10.png" /><figcaption>**Figure 5.** *Energy transfer ratio depending on the multiscale ratio*</figcaption></figure>

**5. 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.

**Acknowledgements**

We thank the DGA (Délégation Générale pour l'Armement) for financial support.
6. References

Daw M. S., “Model of metallic cohesion : the embedded-atom method”, *Physical Review B*, vol. 39, n° 11, p. 7441-7452, 1989.

Dhia H. B., Rateau G., “The Arlequin Method as a flexible engineering design tool”, *International journal for numerical methods in engineering*, vol. 62, p. 1442-1462, 2005.

Lions J. L., Strauss W. A., “Some non linear evolution equations”, *Bulletin de la S.M.F.*, vol. 93, p. 43-96, 1965.

Liu W. K., Karpov E. G., Park H. S., *Nano Mechanics and Materials*, Wiley, 2006.

Marx D., Hutter J., “Ab initio molecular dynamics: Theory and Implementation”, *Modern methods and algorithms of quantum chemistry*, NIC Series, vol. 1, p. 301-449, 2000.

Moës N., Dolbow J., Belytschko T., “A finite element method for crack growth without remeshing”, *International Journal for Numerical Methods in Engineering*, vol. 46, p. 131-150, 1999.

Remmers J. J. C., de Borst R., Needleman A., “A cohesive segments method for the simulation of crack growth”, *Computational Mechanics*, vol. 32, n° 1-2, p. 69-77, 2003.

Sutmann G., “Classical Molecular Dynamics”, *Quantum simulations of complex many-body systems: from theory to algorithms - Lecture Notes, NIC Series*, vol. 10, p. 211-254, 2002.

Xiao S. P., Belytschko T., “A bridging method for coupling continua with molecular dynamics”, *Comput. Methods Appl. Mech. Engrg.*, vol. 193, p. 1645-1669, 2004.

Zhou S. J., Beazley D. M., Lomdahl P. S., Holian B. L., “Large-scale Molecular Dynamics simulations of three-dimensional ductile failure”, *Physical Review Letters*, vol. 78, n° 3, p. 479-482, 1997.