# A NonSmooth Contact Dynamics-based multi-domain solver

Frédéric Perales, Frédéric Dubois, Yann Monerie, Bruno Piar, Laurent Stainier

► To cite this version:

Frédéric Perales, Frédéric Dubois, Yann Monerie, Bruno Piar, Laurent Stainier. A NonSmooth Contact Dynamics-based multi-domain solver: Code coupling (Xper) and application to fracture. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2010, 21 (3-6), pp.242-253. 10.3166/ejcm.19.389-417. hal-00806626

HAL Id: hal-00806626

https://hal.science/hal-00806626

Submitted on 11 Oct 2017

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 Attribution 4.0 International License
# A Non Smooth Contact Dynamics-based multi-domain solver : code coupling (Xper) and application to fracture.

Frederic Perales<sup>a,c</sup> — Frederic Dubois<sup>b,c</sup> — Yann Monerie<sup>a,c</sup> — Bruno Piar<sup>a,c</sup> — Laurent Stainier<sup>d</sup>

<sup>a</sup> Institut de Radioprotection et de Sûreté Nucléaire, DPAM, BP3, 13115 Saint Paul Lez Durance Cedex, France, frederic.perales,yann.monerie,bruno.piar@irsn.fr

<sup>b</sup> Laboratoire de Mécanique et Génie Civil, CC048, Place E. Bataillon, 34095 Montpellier Cedex 5, dubois@lmgc.univ-montp2.fr

<sup>c</sup> Laboratoire de Micromécanique et d'Intégrité des Structures, IRSN DPAM-CNRS - Université de Montpellier 2

<sup>d</sup> Institut de Recherche en Génie Civil & Mécanique (UMR 6183 CNRS), 1 rue de la Noë, BP 44321 Nantes Cedex 3, laurent.stainier@ec-nantes.fr

<footnote>RÉSUMÉ. Cet article présente d'une part une stratégie de modélisation dédiée à la simulation micromécanique des interactions entre corps, et d'autre part, sa mise en œuvre numérique. Cette stratégie repose sur une formulation de type décomposition de domaines d'une méthode multicorps périodique dans le cadre de l'approche Non Smooth Contact Dynamics de Moreau (Moreau, 1988). Les potentialités de cette méthode sont illustrées par la complexité des interactions possibles: interactions entre éléments d'une discrétisation, entre corps discrétisés ou rigides, en compression (contact) lente ou sous impact, en glissement (frottement) ou en traction (fissuration-rupture), etc. La plateforme numérique associée, Xper, repose sur une architecture orientée objet composée de bibliothèques indépendantes spécifiquement développées pour: (1) la résolution des équations aux dérivées partielles, (2) la gestion des interactions superficiques complexes et l'intégration en temps associée, (3) la modélisation des comportements volumiques. La pertinence numérique de l'approche est illustrée sur des exemples de fissuration de matériaux hétérogènes.</footnote>
**ABSTRACT.** This paper presents a micromechanical modeling strategy for complex multibody interactions and the associated numerical framework. The strategy rests on a periodic multi-body method in the framework of the Non Smooth Contact Dynamics approach of Moreau (Moreau, 1988) extended to classical domain decomposition problems. Many complex interactions can be taken into account : interactions between discrete elements, between discrete or rigid bodies, (quasistatic) contact or impact, friction or adhesion, decohesion (cracking), etc. The associated numerical platform, Xper, is composed of three independent libraries with Object Oriented Programming. The libraries are specifically developed for : (1) the solution of systems of partial differential equations (PDEs), (2) the modeling of complex interaction problems and the time discretization, (3) the integration of complex non linear constitutive models. The ability of this computational approach is illustrated by two examples of fracture of heterogeneous materials.

**MOTS-CLÉS** : NonSmooth Contact Dynamics, Modèle de Zone Cohésive Frottante, Méthode de Décomposition de Domaine, Structures Périodiques, Fissuration Dynamique, Matériaux Hétérogènes

**KEYWORDS** : NonSmooth Contact Dynamics, Frictional Cohesive Zone Model, Domain Decomposition Method, Periodic Media, Dynamic Fracture, Heterogeneous Media
# 1. Introduction

Dynamic crack propagation in heterogeneous materials is a complex problem and has become a challenge in many engineering domains. Since the heterogeneities (microcavities, hard inclusions, brittle precipitates, ...) have major effects on the overall dynamic fracture, a relevant way to investigate the problem is to consider models at the scale of the heterogeneities. At this scale, the fracture can be described using micromechanical concepts, and recent advances in computer technologies make possible the simulation of the overall nonlinear response of a given microstructure.

In this context, a new computational micromechanical approach is developed to analyse the effects of the microstructure heterogeneity on the overall material behavior submitted to static or transient loadings. This approach is based both on the concept of Frictional Cohesive Zone Model and on a multibody method in the context of the NonSmooth Contact Dynamics (NSCD). In particular, the NSCD approach aims to solve dynamic frictional contact problems without regularization nor penalization techniques (Moreau, 1988; Jean, 1999). Since periodic formulations are well adapted to micromechanical studies, a two field Finite Element is written and the framework is extended to this formulation (Perales *et al.*, 2006; Perales *et al.*, 2008).

The scope of the NSCD framework can simply be extended to classical domain decomposition problems. The domain decomposition methods (DDM) solve a standard boundary value problem by splitting it in smaller problems on subdomains and managing “continuity” conditions between subdomains. Many techniques exist to enforce continuity of the solution in quasi-static (Dodds *et al.*, 1980; Magoules *et al.*, 2006; Glowinski *et al.*, 1990) or dynamic (Herry *et al.*, 2002; Gravouil *et al.*, 2001). In this paper, the NSCD framework is written as a dual Schur approach where the continuity condition is replaced by any interaction law. The applications of this method concern parallelism, multiphysics or multibody interactions.

Since the implementation of the micromechanical framework would involve a high programming cost, the development strategy is to reuse and to extend existing specialized libraries. The software developed here, called Xper, is based on the coupling of three libraries (using Fortran90/C++ mixed programming) and thus takes advantage of the ability of each of them. Each library has a clear meaning from the mechanical point of view :

- LMGC90 is dedicated to the surfacic interaction part : model, solving method. It relies on the NonSmooth Contact Dynamics approach (LMGC90, 2009),

- PELICANS is dedicated to the bulk part : periodic Finite Element modeling (PELICANS, 2009),

- MatLib is dedicated to complex constitutive models (Stainier *et al.*, 2003). It is embeded in PELICANS.

The ability of the proposed strategy and of the platform is illustrated on the dynamic fracture of metal matrix composite with brittle inclusions under transient loading and on intergranular fracture in a periodic medium.
The paper is organized as follows. In section 2, the main elementary classes of multibody interactions that can summarize most of the complex interaction problems are recalled. In section 3, the interactions modeling (multibody, NSCD, periodic and DDM) is described. The aforementioned case study are presented in section 4.

## 2. The modeling strategy

The main purpose of this paper is to model and simulate divided and fractured media as a multibody system with interactions. From a mechanical point of view, these interactions correspond to any compression-sliding-traction situation. In particular, we have here in mind to model contact, impacts, dry or lubricated friction, adhesion, cohesion, multiple cracking, failure or any combinations as frictional contact with cohesion (Raous *et al.*, 1999; Perales *et al.*, 2006). Moreover these interactions can take place between rigid or deformable media. From a geometrical point of view, any complex interaction belongs to one of the three following situations : (1) the body-to-body interaction, (2) the cluster-to-cluster interaction, (3) the body-to-cluster interaction. In order to classify these three situations, the two first are detailed below.

– *Body-to-body interaction.* For example, this situation occurs when any surface of a Finite Element mesh is susceptible to surface interaction with its surrounding (see Figure 1). This class of interactions includes the failure of bulk materials or of heterogeneous materials, or transgranular failure. In this kind of applications, the framework is known as the cohesive/volumetric finite element approach (Xu *et al.*, 1994; Cachardo *et al.*, 1996; Raous *et al.*, 1999; Jean *et al.*, 2001; Perales *et al.*, 2006). This concept is extended to any interaction of rigid bodies. Therefore, the body-to-body interaction stands for mesh-to-mesh interaction, rigid-to-rigid interaction or mesh-to-rigid interaction,

<figure><img src="image_2.png" /><figcaption>**Figure 1.** *body-to-body interaction : each finite element is a body.*</figcaption></figure>
– Cluster-to-cluster interaction. The bodies can be cluster-meshed domains or riygids. A surface behavior may be introduced between the bodies (for example, see Figure 2 for meshed bodies). This class of interactions includes intergranular failure (Zavattieri *et al.*, 2001; Vincent *et al.*, 2007; Vincent, 2007), inclusion/matrix interface or granular media (Azema *et al.*, 2006; Renouf *et al.*, 2005; Chetouane *et al.*, 2005).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** *cluster-to-cluster : surface behavior between meshed bodies.*</figcaption></figure>

For the particular case of periodic media, the clusters can be subject to periodic conditions (see Figure 3). This class includes the failure of periodic media or of Representative Volume Elements (Perales *et al.*, 2008; Pelissou *et al.*, 2009).

The overall behavior is then obtained by coupling the standard *volumetric* behavior inside the meshes (or rigid bodies) and the *surface* properties between the bodies, taking complex interactions into account. In this study, this coupling is based on the NonSmooth Contact Dynamics (NSCD) framework (Moreau, 1988; Jean, 1999).

In the framework of the NSCD method, a two level resolution is carried out : the standard or periodic *volumetric* problem is solved at global level - by PELICANS (PELICANS, 2009) and MatLib (Stainier *et al.*, 2003) - and the non smooth contact problem is treated at a local level - by LMGC90 (LMGC90, 2009) -.

The corresponding two field modeling framework is detailed next.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** *Periodic media : the clusters are reproduced by periodicity*</figcaption></figure>

# 3. A two field modeling framework

## 3.1. Global level and periodic problem

A specific bulk model is derived to model a dynamical system under periodic conditions. Since the classical dynamic problem can be seen as a subproblem of the periodic one, we focus on the formulation of the periodic problem. In what follows, the main points of the formulation are described (see (Perales et al., 2008) for more details).

Consider a periodic multibody medium $\Omega_0 = \bigcup_{e} \Omega_0^e$. At any boundary of a body $\Omega_0^e$, mixed boundary conditions (given by the interactions depending on the displacement jump) are introduced. In this framework, the deformation gradient field $\mathbf{F} = \nabla \mathbf{u} + \mathbf{I}$ and the first Piola-Kirchhoff stress field $\mathbf{\Pi}$ are assumed to be periodic with the same period as the medium. The average fields over the periodic medium are denoted by $\bar{\mathbf{F}}$, $\bar{\nabla}\mathbf{u}$ and $\bar{\mathbf{\Pi}}$. The fields $\mathbf{F}$, $\nabla\mathbf{u}$ and $\mathbf{\Pi}$ fluctuate around their average values. The local deformation gradient field can be thus split into an overall field (the field if the medium were homogeneous) and a fluctuation denoted $\nabla\mathbf{u}^{\#}$, which takes the presence of heterogeneities into account. Since fracture is expected, the heterogeneities are not only due to the inclusions but also the cracks in the structure. The global displacement field $\mathbf{u}$ admits the following decomposition : $\mathbf{u} = (\bar{\mathbf{F}} - \mathbf{I}) \cdot \mathbf{X} + \mathbf{u}^{\#}$ where $\mathbf{X}$ is the initial position vector and $\mathbf{I}$ is the second-order identity tensor. The Finite Element formulation becomes a two field Finite Element formulation. The local periodic dynamic problem can be written :

find the periodic displacement field $\mathbf{u}^{\#}$, the deformation gradient field $\mathbf{F}$ and the stress field $\mathbf{\Pi}$ verifying :
– relations for each body $\Omega_0^e$:

$$ \left\{ \begin{array}{ll} \nabla \cdot \Pi = \rho \frac{d^2}{dt^2} \mathbf{u}^{\#} & \text{in } \Omega_0^e, \\ \mathbf{F}'(\mathbf{F})' & \text{(Constitutive law)}, \\ \mathbf{F} = \bar{\mathbf{F}} + \nabla \mathbf{u}^{\#} & \text{in } \Omega_0^e, \\ \Pi \cdot N = \mathbf{T}(\mathbf{u}^{\#}) & \text{on } \partial\Omega_0^e, \end{array} \right. \qquad [1] $$

– average relations in the medium $\Omega_0$:

$$ \left\{ \begin{array}{ll} \bar{F}_{ij} & = \bar{F}_{ij}^{\text{imp}}, \\ \bar{\Pi}_{kl} & = \bar{\Pi}_{kl}^{\text{imp}}, \end{array} \right. \qquad [2] $$

where $\mathbf{T}$ is a boundary force given by the considered interaction law, $\mathbf{N}$ is the unit outward normal vector of the body, $\rho$ is the density, the jump symbol $[\mathbf{f}] = \mathbf{f}^+ - \mathbf{f}^-$ is defined as the difference of a field $\mathbf{f}$ on the two facing surfaces (here the superscripts + and - denote the two opposite surfaces), $\bar{F}_{ij}^{\text{imp}}$ are the components of the prescribed macroscopic transformation gradient and $\bar{\Pi}_{kl}^{\text{imp}}$ are the components of the prescribed macroscopic stress, with $ij \neq kl$ and $\{i, j, k, l\} \in \{1, 2, 3\}$.

Note that the classical dynamic problem can be obtained by replacing the periodic displacement field $\mathbf{u}^{\#}$ by the standard displacement field $\mathbf{u}$ in the first equation of [1].

The framework is dedicated to the study of periodic problems embedding non regular interactions. Considering the standard NSCD algorithm, the non smooth contact problem is treated at the local level (Jean, 1999; Jean et al., 2001). At this level, two main points have to be underlined in the writing of the periodic problem:

– dynamics are dedicated to the treatment of the non regular conditions and has to be only introduced at the local level,

– the presence of the heterogeneities are taken into account only by the fluctuation field $\mathbf{u}^{\#}$.

After choosing the admissible spaces, the broken Sobolev space $U^{\#} = \{\mathbf{v} \in [L^2(\Omega_0)]^m, \mathbf{v}|_{\Omega_0^e} \in [H^1(\Omega_0^e)]^m \quad \forall \Omega_0^e, \mathbf{v} \text{ periodic}\}$ and the space of linear transformation $V = \mathcal{L}(\mathbb{R}^m)$ ($m$ is the space dimension) for the velocity field and its periodic part respectively, and considering the kinematically admissible virtual velocity field $\mathbf{v}^* = (\mathbf{v}^*)^{\#} + \dot{\bar{\mathbf{F}}}^* \cdot \mathbf{X}$, the weak unit cell value problem is then obtained:
find $\mathbf{u}^\# \in U^\#$ and $\bar{\mathbf{F}} \in V$ such that :

$$
\left\{
\begin{array}{ll}
\displaystyle \sum_{e=1}^{\tau} \int_{\partial \Omega_0^e} \mathbf{T}([\mathbf{u}^#]) \cdot (\mathbf{v}^*)^\# dS & \\
\displaystyle - \sum_{e=1}^{\tau} \int_{\Omega_0^e} \mathbf{\Pi}(\mathbf{F}) : \nabla(\mathbf{v}^*)^\# dx & \forall (\mathbf{v}^*)^\# \in U^\#; \\
\displaystyle = \sum_{e=1}^{\tau} \int_{\Omega_0^e} \rho \frac{d^2}{dt^2} \mathbf{u}^\# \cdot (\mathbf{v}^*)^\# dx & \\
\displaystyle \sum_{e=1}^{\tau} \int_{\Omega_0^e} \mathbf{\Pi}(\mathbf{F}) : \dot{\mathbf{F}}^* dx = |\Omega_0| \bar{\mathbf{\Pi}} : \dot{\mathbf{F}}^* & \forall \dot{\mathbf{F}}^* \in V.
\end{array}
\right.
\quad [3]
$$

3.2. Non Smooth Contact Dynamics strategy and its extension to periodic problem

3.2.1. Standard dynamic equation and the extension to periodic

Considering that some discontinuities may appear in velocity time evolution, the
standard dynamic problems are written in a semi-discrete form for each body $\Omega_0^e$ :

$$
M \cdot d\dot{\mathbf{q}} = F(\mathbf{q}, \dot{\mathbf{q}}, t)dt + d\mathbf{p}, \quad [4]
$$

$$
\mathbf{q}(t) = \mathbf{q}(t_0) + \int_{t_0}^{t} \dot{\mathbf{q}} dt, \qquad [5]
$$

where *M* is the mass matrix, **q**, **dot_q** and *ddot_q* are respectively the discrete displacement,
velocity and differential measure of velocity, *d*p represents the differential measure of
interaction impulse and *F*(**q**, **dot_q**, *t*) represents the internal and external forces without
the contribution of interaction *d*p (Jean, 1999; Jean *et al.*, 2001). The measure *d*p may
contain both smooth and non-smooth contributions. In this framework, the derivatives
are written in a distribution sense.

Practically the differential measure equation [4] needs to be integrated over an ar-
bitrary time step (even reduced to a shock instant), which gives a momentum balance :

$$
M \cdot (\dot{\mathbf{q}}_{i+1} - \dot{\mathbf{q}}_i) = \int_{t_i}^{t_{i+1}} \mathbf{F}(\mathbf{q}, \dot{\mathbf{q}}, t) dt + \mathbf{p}, \quad [6]
$$

$$
\mathbf{q}(t_{i+1}) = \mathbf{q}(t_i) + \int_{t_i}^{t_{i+1}} \dot{\mathbf{q}} dt. \quad [7]
$$

Using an implicit time integration scheme of [6] and [7] (θ-method) one obtains a
discrete form :
$$ M \cdot (\dot{\mathbf{q}}_{i+1} - \dot{\mathbf{q}}_i) = h(\theta F(\mathbf{q}_{i+1}, \dot{\mathbf{q}}_{i+1}, t) + (1-\theta)F(\mathbf{q}_i, \dot{\mathbf{q}}_i, t)) + \mathbf{p}, \quad [8] $$

$$ \mathbf{q}_{i+1} = \mathbf{q}_i + h(\theta \dot{\mathbf{q}}_{i+1} + (1-\theta)\dot{\mathbf{q}}_i), \quad [9] $$

where $h = t_{i+1} - t_i$.

The linearized formulation of [8] is obtained through a Newton-Raphson method (superscript k stands for iterations) (Jean, 1999; Jean et al., 2001) :

$$ \left\{ \begin{aligned} \tilde{M}^k(\dot{\mathbf{q}}_{i+1}^{k+1} - \dot{\mathbf{q}}_{i+1}^k) &= \mathbf{p}_{\text{free}}^{k+1} + \mathbf{p}_{i+1}^{k+1} \\ \tilde{M}^k &= M + h\theta \frac{\partial F(\mathbf{q}_{i+1}^k, \dot{\mathbf{q}}_{i+1}^k, t)}{\partial \dot{\mathbf{q}}} + h^2\theta^2 \frac{\partial F(\mathbf{q}_{i+1}^k, \dot{\mathbf{q}}_{i+1}^k, t)}{\partial \mathbf{q}} \\ \mathbf{p}_{\text{free}}^{k+1} &= -M(\dot{\mathbf{q}}_{i+1}^k - \dot{\mathbf{q}}_i) + h[(1-\theta)\mathbf{F}_i + \theta \mathbf{F}_{i+1}^k] \end{aligned} \right. \quad [10] $$

which leads to :

$$ \left\{ \begin{aligned} \dot{\mathbf{q}}_{i+1}^{k+1} &= \dot{\mathbf{q}}_{\text{free}}^{k+1} + \mathbf{w}_{i+1}^k \mathbf{p}_{i+1}^{k+1} \\ \dot{\mathbf{q}}_{\text{free}}^{k+1} &= \dot{\mathbf{q}}_{i+1}^k + \mathbf{w}_{i+1}^k \mathbf{p}_{\text{free}}^{k+1} \end{aligned} \right. \quad [11] $$

where $\mathbf{w}_{i+1}^k$ denotes the inverse of the iteration matrix, $h$ is the length of time subinterval $]t_i, t_i + 1]$, subscript $i$ the quantities at a time $t_i$ and $i + 1$ the quantities at time $t_{i+1}$, $\theta \in [0.5, 1]$ and $\mathbf{p}_{i+1}^{k+1}$ is the impulse value.

In a periodic case (see section 3.1 for the formulation), the semi-discrete form [4]-[5] becomes for each body $\Omega_0^e$ (Perales et al., 2008) :

$$ \left\{ \begin{aligned} M \cdot d\dot{\mathbf{q}}^{\#} &= F(\mathbf{q}^{\#}, \dot{\mathbf{q}}^{\#+}, \bar{\mathbf{d}}, \dot{\bar{\mathbf{d}}}, t)dt + d\mathbf{p}, \\ 0 &= G(\mathbf{q}^{\#}, \dot{\mathbf{q}}^{\#+}, \bar{\mathbf{d}}, \dot{\bar{\mathbf{d}}}, t) + \mathbf{K}(t), \end{aligned} \right. \quad [12] $$

$$ \left\{ \begin{aligned} \mathbf{q}(t) &= \mathbf{q}(t_0) + \int_{t_0}^{t} \dot{\mathbf{q}} dt, \\ \bar{\mathbf{d}}(t) &= \bar{\mathbf{d}}(t_0) + \int_{t_0}^{t} \dot{\bar{\mathbf{d}}} dt \end{aligned} \right. \quad [13] $$

where $G(\mathbf{q}^{\#}, \dot{\mathbf{q}}^{\#}, \bar{\mathbf{d}}, \dot{\bar{\mathbf{d}}}, t)$ and $\mathbf{K}(t)$ represents respectively the macroscopic stress and the macroscopic prescribed stress, and $\bar{\mathbf{d}}$ and $\dot{\bar{\mathbf{d}}}$ are respectively the discrete average deformation gradient and the first time derivative.

The linearized formulation [10] is then rewritten as :
$$
\left\{
\begin{aligned}
\tilde{\mathbf{M}}^k \cdot \left\{ \frac{((\dot{\mathbf{q}}_{i+1})^\#)^{k+1} - ((\dot{\mathbf{q}}_{i+1})^\#)^k}{\dot{\mathbf{d}}_{i+1}^{k+1} - \dot{\mathbf{d}}_{i+1}^k} \right\} &= \mathcal{P}_{\text{free}}^{k+1} + \left\{ \begin{matrix} \mathbf{p}_{i+1}^{k+1} \\ 0 \end{matrix} \right\}, \\
\tilde{\mathbf{M}}^k &= \begin{pmatrix}
\mathcal{M} - h\theta \frac{\partial \mathbf{F}}{\partial \dot{\mathbf{q}}^\#} - h^2\theta^2 \frac{\partial \mathbf{F}}{\partial \mathbf{q}^\#} & -h\theta \frac{\partial \mathbf{F}}{\partial \dot{\mathbf{d}}} - h^2\theta^2 \frac{\partial \mathbf{F}}{\partial \bar{\mathbf{d}}} \\
-h\theta \frac{\partial \mathbf{G}}{\partial \dot{\mathbf{q}}^\#} - h^2\theta^2 \frac{\partial \mathbf{G}}{\partial \mathbf{q}^\#} & -h\theta \frac{\partial \mathbf{G}}{\partial \dot{\mathbf{d}}} + h^2\theta^2 \frac{\partial \mathbf{G}}{\partial \bar{\mathbf{d}}}
\end{pmatrix} \\
\mathcal{P}_{\text{free}}^{k+1} &= \left\{ \begin{matrix} -\mathcal{M} \cdot ((\dot{\mathbf{q}}_{i+1})^\# - \dot{\mathbf{q}}_i^\#) + h[(1-\theta)\mathbf{F}_i + \theta\mathbf{F}_{i+1}^k] \\ h\{(1-\theta)(\mathbf{G}_i + \mathbf{K}_i) - \theta(\mathbf{G}_{i+1}^k + \mathbf{K}_{i+1}^k)\} \end{matrix} \right\}
\end{aligned}
\right.
\quad [14]
$$

The unknowns of the periodic problem are $\dot{\mathbf{q}}^\#$ and $\dot{\bar{\mathbf{d}}}$. A mapping $\mathcal{P}$ permits to recover the standard discrete velocity $\dot{\mathbf{q}}$ from $\dot{\mathbf{q}}^\#$ and $\dot{\bar{\mathbf{d}}}:$

$$ \dot{\mathbf{q}} = \mathcal{P}^{-1} \left\{ \begin{matrix} \dot{\mathbf{q}}^\# \\ \dot{\bar{\mathbf{d}}} \end{matrix} \right\}, \quad [15] $$

with $\mathcal{P}$ such that :

$$
\left\{
\begin{array}{l}
\mathcal{P} : \mathbb{R}^{N_{\mathbf{u}}^{\text{dof}}} \longrightarrow \mathbb{R}^{N_{\mathbf{u}}^{\text{dof}} + N_{\bar{\mathbf{F}}}^{\text{dof}}} \\
\dot{\mathbf{q}} \longmapsto \left\{ \begin{matrix} \dot{\mathbf{q}}^\# \\ \dot{\bar{\mathbf{d}}} \end{matrix} \right\} = \left\{ \begin{matrix} \dot{\mathbf{q}} - \text{Rel} \langle \nabla \dot{\mathbf{q}} \rangle_d \cdot \mathbf{X}_d \\ \langle \nabla \dot{\mathbf{q}} \rangle_d \end{matrix} \right\},
\end{array}
\right.
\quad [16]
$$

where $\text{Rel} \in \mathbb{R}^{N_{\mathbf{u}}^{\text{dof}}} \times \mathbb{R}^m$ is a discrete mapping such that the periodic velocity is given by $\dot{\mathbf{q}}^\# = \dot{\mathbf{q}} - \text{Rel} \langle \nabla \dot{\mathbf{q}} \rangle_d \cdot \mathbf{X}_d$, $\mathbf{X}_d$ is the discrete position vector, $\langle \cdot \rangle_d$: $\mathbb{R}^{N_{\mathbf{u}}^{\text{dof}}} \times \mathbb{R}^m \longrightarrow \mathbb{R}^{N_{\bar{\mathbf{F}}}^{\text{dof}}}$ is the discrete average mapping over $\Omega_0$, $m$ is the space dimension and $N_{\mathbf{u}}^{\text{dof}}$ and $N_{\bar{\mathbf{F}}}^{\text{dof}}$ are the number of degrees of freedom of the discrete velocity $\dot{\mathbf{q}}$ and of the first derivative average deformation gradient $\dot{\bar{\mathbf{d}}}$ respectively.

At the iteration $k+1$, the discrete velocity and the discrete free velocity write :

$$ \dot{\mathbf{q}}_{i+1}^{k+1} = \mathcal{P}^{-1} \left\{ \begin{matrix} \dot{\mathbf{q}}_{i+1}^{k+1} \\ \dot{\bar{\mathbf{d}}}_{i+1}^{k+1} \end{matrix} \right\}, \quad \dot{\mathbf{q}}_{\text{free}}^{k+1} = \mathcal{P}^{-1} \left\{ \begin{matrix} ((\dot{\mathbf{q}}_{\text{free}}^\#)^{k+1}) \\ \dot{\bar{\mathbf{d}}}_{\text{free}}^{k+1} \end{matrix} \right\}. \quad [17] $$

The local problem at the contact level is thus solved using the standard NSCD algorithm without any modification. The two field periodic formulation can thus be seen as a simple extension at the global level of the NSCD algorithm, that is to say the resolution of the periodic Finite Element problem.
### 3.2.2. Kinematic relations for frictional contact

Writing interaction behavior in a discrete form requires to define mappings $H$ between bulk unknowns $(\dot{\mathbf{q}}, \mathbf{p})$ and interaction unknowns (see Figure 4) : relative velocity $(\dot{\mathbf{U}})$ and impulse $(\mathbf{R})$. Considering one interaction (index $\alpha$) between two bodies ($c$ and $a$) and using classical kinematic relations, the relative velocity writes :

$$ \dot{\mathbf{U}}^{\alpha} = H_{c}^{T\alpha} \dot{\mathbf{q}}_{c} - H_{a}^{T\alpha} \dot{\mathbf{q}}_{a} = [H_{c}^{T\alpha}, -H_{a}^{T\alpha}] \cdot \begin{pmatrix} \dot{\mathbf{q}}_{c} \\ \dot{\mathbf{q}}_{a} \end{pmatrix} \quad [18] $$

Due to duality consideration, one may write the contribution of the interaction $\alpha$ to the global impulse as :

$$ \begin{pmatrix} \mathbf{p}_c \\ \mathbf{p}_a \end{pmatrix} = \begin{pmatrix} H_c^\alpha \\ -H_a^\alpha \end{pmatrix} \mathbf{R}^\alpha \quad [19] $$

Using the linear mappings $H$, the dynamic system of equations is rewritten as :

$$ \left\{ \begin{aligned} \dot{\mathbf{U}}^{\alpha} &= \dot{\mathbf{U}}_{\text{free}}^{\alpha} + \mathbf{W}^{\alpha\alpha} \mathbf{R}^{\alpha} \\ \mathbf{W}^{\alpha\alpha} &= H_{c}^{T\alpha} \mathbf{w} H_{c}^{\alpha} - H_{a}^{T\alpha} \mathbf{w} H_{a}^{\alpha} \end{aligned} \right. \quad [20] $$

One can generalize the previous form in case of a multi-interaction situation, simply by including an additional contribution :

$$ \left\{ \begin{aligned} \dot{\mathbf{U}}^{\alpha} &= \dot{\mathbf{U}}_{\text{free}}^{\alpha} + \mathbf{W}^{\alpha\alpha} \mathbf{R}^{\alpha} + \sum_{\beta \neq \alpha} \mathbf{W}^{\alpha\beta} \mathbf{R}^{\beta} \\ \mathbf{W}^{\alpha\beta} &= H_{c}^{T\alpha} \mathbf{w} H_{c}^{\beta} - H_{a}^{T\alpha} \mathbf{w} H_{a}^{\beta} \end{aligned} \right. \quad [21] $$

To close the problem, one needs an interaction law relating relative velocity and impulse. Numerous choices are possible depending on the phenomenology of the interaction, for example, unilateral condition (velocity Signorini condition (Moreau, 1988)), friction (Coulomb's law) or cohesion (Jean et al., 2001; Perales et al., 2006).

As detailed further, the general traction-sliding-compression interactions presented in the 'modeling strategy' (section 2) have to include sliding-compression relations which are causing impulse conditions. For sake of clarity, unilateral frictional contact laws of Signorini-Coulomb type are presented but other choices can be made (e.g. Tresca friction law). These relations write :

$$ -R_N \in \partial I_{\mathbb{R}_+}(U_N), \quad [22] $$

$$ \mathbf{R}_T \in \partial_{\dot{\mathbf{U}}_T} \left( \mu |R_N| \|\dot{\mathbf{U}}_T\| \right), \quad [23] $$

where $I_K$ is the indicator function of the set $K$, $\mu$ is the Coulomb friction coefficient and $\mathbf{U}$ is decomposed such as $\mathbf{U} = U_N \mathbf{n} + \mathbf{U}_T$ with $\mathbf{n}$ the unit normal vector of the frictional contact zone.
The principle of the global to local mapping and the associated variables are summarized in a standard case and in a periodic case respectively in Figure 4 and in Figure 5. In particular, Figure 5 shows that the periodic mapping concerns only the global level of the algorithm and that the kinematic relations and local NSCD solver are not affected.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** *The standard NSCD algorithm*</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** *The NSCD algorithm extended to the periodic formulation*</figcaption></figure>
### 3.3. NSCD as a Domain Decomposition Method

In classical domain decomposition methods (DDM) one assumes contiguous subdomains and tries to enforce continuity of primal (displacement) and/or dual (force) interface unknowns. Considering quasi-static problems, there exist various techniques to ensure this continuity which differ in the treatment of the interface between subdomains. In the following, we consider techniques without overlapping between subdomains. Among them one can consider primal Schur-type approaches (where continuity of the displacement field is enforced at the interfaces) (Dodds et al., 1980), dual Schur-type approaches (where equilibrium of interface forces is enforced through Lagrange multiplier) (Magoules et al., 2006) and mixed approaches where one tries to achieve both conditions (Glowinski et al., 1990). Considering dynamical problems one can potentially use the three previous approaches. But as mentioned by Gravouil et al. (Herry et al., 2002; Gravouil et al., 2001), in dynamics, the choice of which kinematic quantity is continuous at the interface (displacement, velocity or acceleration) is difficult.

In the case of contact problems various possibilities exist. Dureissex et al. (Dureisseix et al., 2001) propose an extension of a dual approach (FETI). Alart et al. (Alart et al., 2003) propose an approach where contact is embedded in a continuous domain. More recently, in Nineb (Nineb et al., 2007) and Icéta (Iceta et al., 2009), domain decomposition approaches are proposed for non smooth problems, where the continuity condition is enforced through the bulk model.

Our purpose in this part is to show that the NSCD method may be written as a dual-Schur approach where the basic continuity condition may be replaced by any interaction law. In that case the continuity condition between subdomains is not written on the bulk part but through an interaction law.

For the sake of simplicity we will only consider in the following two domains in interaction. In the context of the dual Schur formulation one may rewrite the standard discrete problem ([10] and [18]) involving two subdomains (a and c) :

$$ \left( \begin{array}{ccc} \tilde{\mathbf{M}}_{a}^{k+1} & 0 & H_{a} \\ 0 & \tilde{\mathbf{M}}_{c}^{k+1} & -H_{c} \\ H_{a}^{T} & -H_{c}^{T} & 0 \end{array} \right) \left( \begin{array}{c} \dot{\mathbf{q}}_{a}^{k+1} \\ \dot{\mathbf{q}}_{c}^{k+1} \\ \mathbf{R}^{k+1} \end{array} \right) = \left( \begin{array}{c} \tilde{\mathbf{M}}_{a}^{k+1} \dot{\mathbf{q}}_{a}^{k} + \mathbf{p}_{a, \text{free}}^{k+1} \\ \tilde{\mathbf{M}}_{c}^{k+1} \dot{\mathbf{q}}_{c}^{k} + \mathbf{p}_{c, \text{free}}^{k+1} \\ -\dot{\mathbf{U}}^{k+1} \end{array} \right) \quad [24] $$

where the basic continuity condition of DDM written in velocity ($\dot{\mathbf{U}}^{k+1} = H_a^T \dot{\mathbf{q}}_a^{k+1} - H_c^T \dot{\mathbf{q}}_c^{k+1}$) is replaced by a more general implicit interaction law ($\text{law}(\mathbf{U}^{k+1}, \mathbf{R}^{k+1}) = \text{true}$).

One can compute the solution splitting the problem in two stages :

1) free motion

$$ \left( \begin{array}{ccc} \tilde{\mathbf{M}}_{a}^{k+1} & 0 & 0 \\ 0 & \tilde{\mathbf{M}}_{c}^{k+1} & 0 \\ 0 & 0 & 0 \end{array} \right) \left( \begin{array}{c} \dot{\mathbf{q}}_{a, \text{free}}^{k+1} \\ \dot{\mathbf{q}}_{c, \text{free}}^{k+1} \\ \mathbf{R}^{k+1} \end{array} \right) = \left( \begin{array}{c} \tilde{\mathbf{M}}_{a}^{k+1} \dot{\mathbf{q}}_{a}^{k} + \mathbf{p}_{a, \text{free}}^{k+1} \\ \tilde{\mathbf{M}}_{c}^{k+1} \dot{\mathbf{q}}_{c}^{k} + \mathbf{p}_{c, \text{free}}^{k+1} \\ 0 \end{array} \right) \quad [25] $$
The free motion is computed at the global level and is independent of interaction conditions. It can be computed by any classical FEM library such as PELICANS.

2) interaction motion

$$ \left( \begin{array}{ccc} \tilde{\mathbf{M}}_a^{k+1} & 0 & H_a \\ 0 & \tilde{\mathbf{M}}_c^{k+1} & -H_c \\ H_a^T & -H_c^T & 0 \end{array} \right) \left( \begin{array}{c} \dot{\mathbf{q}}_{a,\text{inter}}^{k+1} \\ \dot{\mathbf{q}}_{c,\text{inter}}^{k+1} \\ \mathbf{R}^{k+1} \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \\ (H_c^T \dot{\mathbf{q}}_{c,\text{free}}^{k+1} - H_a^T \dot{\mathbf{q}}_{a,\text{free}}^{k+1}) - \dot{\mathbf{U}}^{k+1} \end{array} \right) \quad [26] $$

The interaction motion can be computed at the local level solving two subproblems (omitting $^{k+1}$):

$$ \left\{ \begin{array}{l} \dot{\mathbf{U}} = \dot{\mathbf{U}}_{\text{free}} + \mathbf{W}\mathbf{R} \\ \text{law}(\mathbf{U}, \mathbf{R}) = \text{true} \end{array} \right. \qquad [27] $$

where $\dot{\mathbf{U}}_{\text{free}} = H_c^T \dot{\mathbf{q}}_{c,\text{free}} - H_a^T \dot{\mathbf{q}}_{a,\text{free}}$ and $\mathbf{W} = H_c^T \tilde{\mathbf{M}}_c^{-1} H_c - H_a^T \tilde{\mathbf{M}}_a^{-1} H_a$, and

$$ \dot{\mathbf{q}}_{a,\text{inter}} = \tilde{\mathbf{M}}_a^{-1} H_a \mathbf{R} \qquad [28] $$

$$ \dot{\mathbf{q}}_{c,\text{inter}} = \tilde{\mathbf{M}}_c^{-1} H_c \mathbf{R} \qquad [29] $$

The local solver needs the projection of the free velocity on the boundary of the subdomains and the projection of the inverse of the pseudo mass matrix $\tilde{\mathbf{M}}_a$ or $\tilde{\mathbf{M}}_c$.

The resulting velocity is:

$$ \dot{\mathbf{q}}_a^{k+1} = \dot{\mathbf{q}}_{a,\text{free}}^{k+1} + \dot{\mathbf{q}}_{a,\text{inter}}^{k+1} \qquad [30] $$

$$ \dot{\mathbf{q}}_c^{k+1} = \dot{\mathbf{q}}_{c,\text{free}}^{k+1} + \dot{\mathbf{q}}_{c,\text{inter}}^{k+1} \qquad [31] $$

Note that the DDM problem can simply be extended to the periodic formulation using the relations [14] and the mapping [16].

### 3.4. Multi domain solver

The classical NSCD approach, and its periodic extension, rely on a Non Linear Gauss Seidel (NLGS) algorithm to solve the problem.

The spirit of the method is the following. Considering, one by one, the local systems to solve for each contact $\alpha$:

$$ \left\{ \begin{array}{l} \dot{\mathbf{U}}^{\alpha} = \dot{\mathbf{U}}_{\text{free}}^{\alpha} + \mathbf{W}^{\alpha\alpha}\mathbf{R}^{\alpha} + \sum_{\beta \neq \alpha} \mathbf{W}^{\alpha\beta}\mathbf{R}^{\beta} \\ \text{Law}(\mathbf{g}^{\alpha}, \mathbf{U}^{\alpha,n}, \mathbf{U}^{\alpha,t}, \mathbf{R}^{\alpha,n}, \mathbf{R}^{\alpha,t}) = \text{true} \end{array} \right. \qquad [32] $$
the contributions due to other contact ($\beta \neq \alpha$) are frozen taking updated values if $\beta < \alpha$ or old values if $\beta > \alpha$.

The $\alpha$ block of equations can be solved using :

## 2D context

- An explicit uncoupled resolution if $W^{\alpha\alpha}$ is diagonal

- A coupled $(n, t)$ graph intersection (Jean, 1999)

- A pseudo-potential approach (bi-potential) (De Saxcé *et al.*, 1991)

- Linear Complementarity Problem (LCP) as local solver

## 3D context

- An explicit resolution if $W^{\alpha\alpha}$ is diagonal (Renouf, 2004)

- A pseudo-potential approach (De Saxcé *et al.*, 1991)

- A Generalized Newton algorithm (Alart *et al.*, 1991)

- Linear Complementarity Problem (LCP) as local solver

## 3.5. Global-to-local strategy and dedicated platform

The corresponding global-to-local strategy is summarized for one iteration of the Newton-Raphson algorithm in Figure 6.

The strategy rests on the local/global levels :

- at the global level, the Finite Element method (including the periodic extension) is taken into account,

- at the local level, the standard non smooth contact framework is used.

For lower programming cost, the software strategy development retained here was not to develop the entire software *from scratch* but to reuse and extend existing libraries. It permits to take advantages of each library update while developing the whole software.

The software platform, called Xper ('eXtended cohesive zone model and PERiodic homogenization'), respects the two levels strategy (Figure 6) :

- at the global level, the extended Finite Element method is managed by the PELICANS library (PELICANS, 2009). This software, developed by the French 'Institut de Radioprotection et de Sûreté Nucléaire' (IRSN) is a toolbox for the implementation of various numerical methods dedicated to the solution of systems of partial differential equations (PDEs). To take into account complex non linear constitutive models, the library is coupled with the MatLib library, developed by Stainier (Stainier *et al.*, 2003),

- at the local level, the complex interactions between bodies are taken into account by the LMGC90 library (LMGC90, 2009), developed by Dubois and Jean (Dubois *et al.*, 2009).
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Global-to-Local coupling : one iteration of Newton-Raphson. In dashed boxes : periodic case.</figcaption></figure>
al., 2003). LMGC90 is a platform for the modeling interaction problems including multi-physics.

For more details concerning the implementation of the software, see the appendix A.

In the following, both the abilities of the numerical strategy and the Xper software are illustrated on fracture of heterogeneous media.

# 4. Application to the fracture of heterogeneous materials

## 4.1. Cohesive zone interaction law

Examples presented deal with fracture of heterogeneous materials. The fracture model is based on a cohesive-volumetric micromechanical approach involving Frictional Cohesive Zone Models (FCZM). These models rest on coupling traction-separation interaction law to some frictional contact model. The cohesive-friction coupling is a key concept of fracture of heterogeneous media: whatever the overall loading of such a media (including pure mode I loadings), interfaces between phases with different Poisson ratio are locally subjected to combination of shear and tensile loading.

The interaction law relating the displacement jump $\mathbf{U}$ to the stress $\mathbf{R}$, used to close the NSCD framework, is obtained introducing a cohesive stress, called $\mathbf{R}^{coh}$, in the Signorini-Coulomb problem [22]-[23] (Perales et al., 2006):

$$ -(R_N + R_N^{coh}) \in \partial I_{\mathbb{R}_+}(U_N), \qquad [33] $$

$$ (\mathbf{R}_T + \mathbf{R}_T^{coh}) \in \partial_{\dot{\mathbf{U}}_T} (\mu |R_N + R_N^{coh}| |\dot{\mathbf{U}}_T|). \qquad [34] $$

$$ \mathbf{R}^{coh} = \beta \left( C_N \mathbf{n} \otimes \mathbf{n} + C_T \frac{\mathbf{U}_T \otimes \mathbf{U}_T}{\|\mathbf{U}_T\|^2} \right) \cdot \mathbf{U} \qquad [35] $$

where $\mathbf{n}$ is the unit normal vector of the cohesive zone, $I_K$ is the indicator function of the set $K$, $\mu$ is the Coulomb friction coefficient, $C_N$ and $C_T$ denote respectively the initial normal and tangential stiffness of the perfect interface (MPa/m).

The surface variable $\beta$, initially introduced by Fremond (Fremond, 1987), plays the role of a surface damage variable. The evolution law of this variable is governed by eqs. [36] and [37], where the function $g$ describes the weakening process leading
from perfect interface to crack ($\beta = 1$: the interface is undamaged, $0 < \beta < 1$: the interface is partially damaged and $\beta = 0$: the interface is fully damaged) :

$$ \beta = \min(g(\|\mathbf{U}\|), g(\|\mathbf{U}\|_{\max})), \qquad [36] $$

$$ g(x) = \begin{cases} \beta_0 & \text{if } x \le \delta_0, \\ \beta_0 \frac{\delta_0}{x} \left( 1 - \left( \frac{x - \delta_0}{\delta_c - \delta_0} \right)^2 \right) & \text{if } \delta_0 < x < \delta_c, \\ 0 & \text{if } x \ge \delta_c, \end{cases} \qquad [37] $$

where $\delta_0 = \frac{R_{\max}}{2} \left( \frac{1}{C_N} + \frac{1}{C_T} \right)$, $\delta_c = \frac{3}{2} \left( \frac{w}{R_{\max}} + \frac{\delta_0}{6} \right)$, $0 \le \beta_0 \le 1$ is an initial surface damage, $w$ is a reference fracture energy ($J/m^2$), $R_{\max}$ is the maximum value of the cohesive stress (MPa), $\|\mathbf{U}\|_{\max}$ is the maximum value reached by $\|\mathbf{U}\|$ during the fracture process. In a 2D case, Figure 7 shows respectively (a) the normal behavior (with $\|\mathbf{U}_T\| = 0$) and (b) the tangential behavior (with $R_N$ constant) associated with the FCZM [33]-[37].

<figure><img src="image_8.png" /><figcaption>**Figure 7.** The 2D FCZM: (a) normal behavior ($U_T = 0$) and (b) tangent behavior ($U_N = 0$, $R_N$ constant).</figcaption></figure>

This Frictional Cohesive Zone Model takes into account the progressive damage between two bodies and the post-fracture frictional contact on the created crack lips. Different models can be used according to the material (brittle, ductile) by specifying the different evolution laws [36]-[37] (Michel et al., 1994; Tvergaard, 1990; Alfano et al., 2001; Perales et al., 2006).
## 4.2. Examples

In what follows, the finite element discretization is based on linear displacement triangular elements that are arranged in a "crossed-triangle" quadrilateral pattern. The analysis considers 2D plane-strain conditions.

The considered microstructure is composed of a metal matrix (Zircaloy-4) and inclusions ($\delta$-hydrides). The Zircaloy-4 behavior is assumed to be elastoplastic (J2 plasticity, Young Modulus $E = 99\text{GPa}$, Poisson's ratio $\nu = 0.325$, Yield in tension $\sigma_0 = 450\text{MPa}$, Hardening Modulus $H_Y = 850\text{MPa}$) (Balourdet et al., 1999; Cazalis et al., 2005) and hydrides to be neo-Hookean ($E = 135\text{GPa}$, $\nu = 0.32$) (Yamanaka et al., 1999; Yamanaka et al., 2001). We assume that the Zircaloy-4 and the hydrides inclusions have the same density $\rho = 7800\text{kg/m}^3$. The FCZM coefficients of Zircaloy-4, zirconium hydrides and Zircaloy-hydride interface are respectively: $C_N^{Zr} = 2 \times 10^{18}\text{Pa/m}$, $w^{Zr} = 0.5\text{J/m}^2$, $R_{\max}^{Zr} = 1\text{GPa}$, $C_N^{ZrH} = 2C_N^{Zr}$, $w^{ZrH} = 0.8w^{Zr}$, $R_{\max}^{ZrH} = 1.25R_{\max}^{Zr}$, $C_N^{Zr-ZrH} = 2C_N^{Zr}$, $w^{Zr-ZrH} = 0.001w^{Zr}$, $R_{\max}^{Zr-ZrH} = 0.045R_{\max}^{Zr}$, if the interface is considered as 'weak' or $w^{Zr-ZrH} = 1000w^{Zr}$, $R_{\max}^{Zr-ZrH} = 45R_{\max}^{Zr}$ if the interface is considered as 'strong'.

Moreover, we assume a low friction coefficient $\mu = 0.05$ and same compliance for the normal and tangential behaviors $C_N = C_T$.

The PELICANS library manages the finite element part of the problem : the structure geometry, the inclusions distribution and morphology, the boundary conditions, the periodicity, the discretization and the global solver. The LMGC90 library takes into account the FCZM between each body. The material properties of the matrix and the inclusions are managed by the MatLib library.

### 4.2.1. *Body-to-body interaction : influence of interface type on the fracture of a brittle heterogeneous material*

This example deals with the fracture of a bimaterial. Each mesh is considered as an independent body (body-to-body interaction, Figure 1 case). In particular, the influence of the interface behavior on the fracture of metal matrix composites is investigated. The considered composite is representative of hydried Zircaloy-based alloys at high burnup which compose cladding of nuclear fuel rods after many years in Pressurized Water Reactor.

The structure is composed by a metal matrix (Zircaloy-4) and rectangular aligned inclusions ($\delta$-hydrides). The structure geometry is a square with length $L = 50\mu\text{m}$. Horizontal velocity boundary conditions on the left and right vertical faces are prescribed (respectively $V = -1\text{m/s}$ and $V = 1\text{m/s}$). A precrack is introduced at the bottom perpendicular to the loading (Figure 8).

Two cases of bonding strength value between the two phases, strong and weak, are considered.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Boundary condition (arrows), matrix (light gray), inclusions (dark gray) and precrack (white line)</figcaption></figure>

Figure 9 shows the fracture features for the weak and strong matrix/inclusion cases. The crack path is significantly influenced by the interface matrix/inclusions bonding strength. For the weak interfaces, the cracks propagate inside the matrix and along the inclusions boundaries due to the formation of microcracks along weak interfaces. In case of strong interfaces, the cracks propagate through the inclusions due to the high cohesive strength with the matrix. These results are consistent with the following criterion (Raous *et al.*, 2002; Raous *et al.*, 2001; He *et al.*, 1989; Siegmund *et al.*, 1997; Martin *et al.*, 1998; Xu *et al.*, 1998) : the transition between deflection and penetration occurs for a ratio of interface fracture energy $w^{int}$ to inclusion fracture energy $w^i$ in the range [0.013, 0.25]. In other words, when $w^{int} < 0.013w^i$ the matrix crack is delfected at the interface (case of weak interface); when $w^{int} > 0.25w^i$ the matrix crack propagates through the inclusion (case of strong interface).

Figure 10 shows the evolution of the stress during cracking process. During the loading, the stress increases linearly. When the precrack propagates into the matrix to the first inclusion, the stress becomes to decrease. In case of strong interfaces, it continues to decrease rapidly due to the propagation of the cracks through the inclusions while in case of weak interfaces, it increases due to the formation of microcracks and their coalescence (precrack is arrested at the matrix/inclusion interface).

The energy release rate decreases with increasing interfacial bonding strength. The energy release for strong interfaces is lower than for weak interfaces. These results show the importance of the interface bonding strength and, in conclusion, a strong interface is more deleterious than a weak interface.

### 4.2.2. *Body-to-body with periodic boundary interaction : fracture of a ductile heterogeneous material*

Since the previous example can be considered as a calculus of a structure behavior, we consider now the same type of heterogeneities but from the point of view of the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Rupture features of the two cases : weak (left) and strong (right)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Stress (MPa) vs strain (%) for the two cases : weak (solid curve) and strong (dashed curve)</figcaption></figure>

material behavior. Moreover, a plastic deformation inside the matrix is considered in order to underline the abilities of the MatLib library.

This example deals with the fracture of a bimaterial with periodic conditions and without precrack (body-to-body interaction, Figure 1 case, with periodic conditions). The example is the same as the previous one (section 4.2.1), except that the boundary conditions are periodic and there is no precrack in the structure. In addition, the bonding strength value between the two phases is assumed to be strong.
A macroscopic strain gradient is prescribed along the horizontal direction.

Without any initial precrack, a multifracture initiates in the inclusions, especially at the locii of high concentration of inclusions (Figure 11). Then, the cracks propagation occur through growth and coalescence and lead to the failure of the structure.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Initiation of cracks (white line) in the inclusions</figcaption></figure>

Figure 12 shows the overall stress-strain curve. The overall behavior is ductile over an uniaxial stress about 530MPa (the yield stress of the matrix being 450MPa): as expected, the presence of elastic inclusions with higher strength than the matrix increases the overall yield stress.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Stress (MPa) vs strain (%) for the two cases: weak (solid curve) and strong (dashed curve)</figcaption></figure>
### 4.2.3. Cluster-to-cluster with periodic boundary interactions : intergranular fracture

This example deals with intergranular fracture in order to show the ability of the method in a cluster-to-cluster strategy. Each grain is a meshed domain related to each other with FCZM (cluster-to-cluster interaction, Figure 2 case). The considered structure is composed by metal grains (Zircaloy) and by two multi-grains inclusions ($\delta$-hydride). The structure geometry is a periodic square with length $L = 100\mu m$ (periodic boundary interaction, Figure 3 case), see Figure 13. A periodic Voronoi tessellation is used. A macroscopic strain gradient rate is prescribed along the horizontal direction.

The bonding strength value between the hard grains and the soft grains is assumed to be weak.

Figure 13 (right) shows the fractured periodic structure. Again as expected, the cracks initiate at the grain boundaries of the hard grains. Two main cracks propagate along the grain boundaries of the soft grains and they join together to form a single one which leads to the failure of the material.

Figure 14 shows the overall stress-strain curve. The overall behavior is brittle due to the weak bonding strength of the grain boundaries.

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Periodic meshed bodies : shades of gray represents the grains and the inclusions are delimited by the white borders (left) and fractured cell (right).</figcaption></figure>

## 5. Conclusion and perspective

This paper presented a numerical framework for the modeling of dynamic crack propagation in heterogeneous materials. The underlying model rests on a multibody approach and complex interactions behavior in a periodic extension of the standard NonSmooth Contact Dynamics framework. This framework has been rewritten as a dual Schur approach with complex interaction law between subdomains instead of “basic” continuity condition. It can be seen as an extension of the classical nonoverlapping domain decomposition methods. In particular, it allows to model complex in-
<figure><img src="image_15.png" /><figcaption>**Figure 14.** *Stress (MPa) vs strain (%) for the periodic domain.*</figcaption></figure>

teractions, discrete media, periodic media, parallelism or multiphysics. The associated
developed software has the following capabilities : dynamic, finite strain, non-linear
behaviors, periodic Finite Element, crack initiation and propagation, nonsmooth post-
fracture behavior. Two examples dealing with the influence of the matrix/inclusion
interfaces on the crack propagation and the intergranular fracture in a heterogeneous
material have been presented to illustrate the ability of Xper and of the modeling stra-
tegy.

This strategy and Xper code can be applied to nuclear safety, for example, to de-
terminate the fuel behavior at high burnup in a nuclear reactor or the ageing effect on
the behavior of nuclear power plant equipment.

A natural extension of the proposed modeling framework is the parallelism. Each
body can be considered as an independent meshed subdomain. In that case, the sub-
problems are solved on parallel computers. Another application is the multiphysics
in which each mathematical problem is posed on a different domain. Furthermore, one
can introduce complex interactions law between the subdomains.
# Appendix

## A. Overview of the software implementation

The architecture design of the software is based on object oriented techniques. Object oriented programming (OOP) provides a clear modular structure for programs and makes it easy to maintain and modify existing code. The fundamental concepts are (Martin, 1996; Martin, 2003; Meyer, 1997) : inheritance, encapsulation, abstraction and polymorphism. One of the principal advantages of OOP techniques over procedural programming techniques is that they enable programmers to create modules that do not need to be changed when a new type of object is added. A programmer can simply create a new object that inherits many of its features from existing objects. This makes object-oriented programs easier to modify.

For low programming cost, the strategy retained here is not to develop the software 'from scratch' but to reuse and extend existing libraries. It permits to take advantages of the libraries updates while developing the whole software. The framework can be splited in two distinct levels (see section 3) : a global Finite Element resolution (including periodic extension) and a local non smooth contact resolution. The Xper architecture respects this local/global levels by the coupling of three libraries : LMGC90, PELICANS and MatLib. The libraries are described in the following.

### A.1. Existing libraries

#### A.1.1. LMGC90

The LMGC90 library (LMGC90, 2009) is developed by Dubois and Jean (Dubois et al., 2003). The software is an open platform under the CECILL License (CECILL, 2005) for modeling interaction problems between elements including multi-physics. It allows to model : granular material made of rigid or deformable bodies and with complex interactions (contact, friction, cohesion, wear, etc.), discrete media, masonry, fracture, etc. The modeling approach is based on an hybrid or extended Finite Element Method (FEM) - Discrete Element Method (DEM), using various numerical strategies such as Molecular Dynamics (MD) or NSCD. In particular, the NSCD algorithm allows to take into account cohesive zone models with contact and friction (see section 3 and 4.1).

The programming language used is Fortran90. Although Fortran90 is a procedural programming language, LMGC90 is implemented in the form of modules using Object Oriented Programming. In particular, the code is open for extension and closed for modification, it respects the 'Open-Closed' principle (Meyer, 1997). LMGC90 can be used as a library or as a standalone software through a macro language.
### A.1.2. PELICANS

The PELICANS library (PELICANS, 2009), developed by the French 'Institut de Radioprotection et de Sûreté Nucléaire' (IRSN) under CECILL-C License (CECILL, 2005), is a toolbox for the implementation of various numerical methods dedicated to the solution of systems of partial differential equations (PDEs).

The PELICANS platform is written in standard C++ and is object-oriented. The constitutive classes are classified in two groups : 'plug-points' and 'service provider'. The 'plug-points' classes play the role of mother classes and permit user classes to plug into the platform. The 'service provider' classes make functionalities available to users. The platform respects the Component-Based Development principles, such as the Design by Contract, command-query separation, Liskov substitution principle, inheritance, naming and self-documentation issues (Martin, 1996; Martin, 2003; Meyer, 1997).

### A.1.3. MatLib

The MatLib library is developed by Stainier (Stainier et al., 2003). MatLib is a material constitutive models library. It is based on a variational formalism of thermo-mechanical constitutive updates (Ortiz et al., 1999; Yang et al., 2006). This formalism rests on potential combination naturally involving an object oriented structure. The architecture permits the creation of new models from existing models and their implementation in the library. The models are accessible to users through a common interface.

The library is written in C++.

## A.2. Software architecture

The Xper rests on the coupling of the three libraries (see Figure 15) using OO and Mixed Programming.

In the current version of the code, the LMGC90 library plays the role of the master program in the coupling. It manages, in particular, the time discretization and the Newton loops. In the future, the master program will be an interpreted program written with high-level and OO programming language like Python and independent of the three libraries.

The choice of the master program being made, each library has a clear mechanical meaning (Figure 15) :

- PELICANS manages the periodic finite element at the global level,

- MatLib provides the non linear constitutive models (global level),

- LMGC90 manages the surface behavior between the bodies at the local level.
The volumetric behavior is obtained by the coupling of the finite element library developed from PELICANS and the constitutive models library MatLib. Since the two libraries are written in C++ language, the coupling is strong. For example, a MatLib object can be instantiated from a PELICANS class. A PELICANS-MatLib interface permits to manage the MatLib objects in the library developed from PELICANS, preserving the independence of each library.

The surface behavior is managed by LMGC90. The coupling between LMGC90 and PELICANS needs Fortran90/C++ Mixed Programming techniques. Two interfaces are developed for passing parameters from C++ to Fortran90 and back, one in PELICANS and the other in LMGC90.

This strategy permits to compile all the libraries independently, and then to link them.

Note that minor changes have been made in the LMGC90 and MatLib libraries.

<figure><img src="image_16.png" /><figcaption>**Figure 15.** *Role of each library*</figcaption></figure>

### A.3. Xper software

The main features of the Xper software are summarized as follows :

- dynamics,

- finite strain,

- non-linear behaviors,

- heterogeneous materials,
- periodic Finite Element,

- complex surface behavior (i.e. contact and friction),

- mixing FEM-DEM.

Note that for the preprocessing, the postprocessing and the linear algebra, the internal capabilities of the PELICANS and LMGC90 libraries are completed by the coupling to external software. These functionalities are thus available in the software.

To reduce the high cost of software development, existing codes have been reused instead of developing them 'from scratch'. According to this strategy, the implementation and the validation of the Xper project, including all the features presented, represents 1 man-year. This is a relatively low development costs compared to developing entirely the features of the libraries in a new platform.

## Bibliographie

Alart P., Barboteu M., Renouf M., « Parallel computational strategies for multicontact problems : applications to cellular and granular media. », *International Journal for Multiscale Computational Engineering*, vol. 1, p. 419-430, 2003.

Alart P., Curnier A., « A mixed formulation for frictional contact problems prone to newton like solution methods », *Computer Methods in Applied Mechanics and Engineering*, vol. 92, p. 353-375, 1991.

Alfano G., Crisfield M., « Finite element interface models for the delamination analysis of laminated composites : mechanical and computational issues », *International Journal for Numerical Methods in Engineering*, vol. 50, p. 1701-1736, 2001.

Azema E., Radjaï F., Peyroux R., Dubois F., Saussine G., « Vibrational dynamics of confined granular materials », *Physical Review E*, vol. 74, n° 3, p. 031302, 2006.

Balourdet M., Bernaudat C., Basini V., Hourdequin N., « The PROMETRA program : Assessment of mechanical properties of Zircaloy-4 fuel properties during a RIA. », *Transactions of the 15th International Conference on Structural Mechanics in Reactor Technology*, Seoul, Korea, p. II-485-492, August 15-20, 1999.

Camacho G., Ortiz M., « Computational modelling of impact damage in brittle materials », *International Journal of Solids and Structures*, vol. 33, p. 2899-2938, 1996.

Cazalis B., Bernaudat C., Yvon P., Desquines J., Poussard C., Averty X., « The PROMETRA program : a reliable material database for highly irradiated Zircaloy-4, Zirlo™ and M5™ fuel claddings. », *18th International Conference on Structural Mechanics in Reactor Technology*, Beijing, China, p. 383-393, August 7-12, 2005. SMiRT18 - C02_1.

CECILL, More information on the web site : http://www.cecill.info/index.en.html, 2005.

Chetouane B., Dubois F., Vinches M., Bohatier C., « NSCD discrete element method for modeling masonry structures », *International Journal for Numerical Methods in Engineering*, vol. 64, p. 65-94, 2005.

De Saxcé G., Feng Z.-Q., « New inequality and functional for contact with friction : the implicit standard material approach », *Mechanics of Structures and Machines*, vol. 19, p. 301-325, 1991.
Dodds J., Lopez L., « Substructuring in linear and nonlinear analysis », *International Journal of Numerical Methods in Engineering*, vol. 15, p. 583-597, 1980.

Dubois F., Jean M., « Une plateforme de développement dédiée à la modélisation des problèmes d'interaction », in M. Potier-Ferry, M. Bonnet, A. Bignonnet (eds), *6ème colloque national en calcul de structures*, vol. 1, Giens, p. 111-118, 2003.

Dureisseix D., Farhat C., « A numerically scalable domain decomposition method for the solution of frictionless contact problems », *International Journal for Numerical Methods in Engineering*, vol. 50, p. 2643-2666, 2001.

Fremond M., « Adhérence des solides. », *Journal de Mécanique Théorique et Appliquée*, vol. 6, n° 3, p. 383-407, 1987.

Glowinski R., Le Tallec P., « Augmented Lagrangian interpretation of the non-overlapping Schartz alternating method », *Domain Decomposition method SIAM*, Philadelphia, USA, p. 224-231, 1990.

Gravouil A., Combescure A., « A multi-time-step explicit-implicit method for non-linear structural dynamics », *International Journal for Numerical Methods in Engineering*, vol. 50, p. 199-225, 2001.

He M. Y., Hutchinson J. W., « Crack deflection at an interface between dissimilar elastic materials », *International Journal of Solids and Structures*, vol. 25, p. 1053-1067, 1989.

Herry B., Valentin L. D., Combescure A., « An approach to the connection between subdomains with non-matching meshes for transient mechanical analysis », *International Journal for Numerical Methods in Engineering*, vol. 55, p. 973-1003, 2002.

Iceta D., Alart P., Dureisseix D., « A multilevel domain decomposition solver suited to nonsmooth mechanical problems », in M. Bercovier, M. Gander, R. Kornhuber, O. Widlund (eds), *Lecture notes in Computational Science and Engineering*, vol. 70, Springer Verlag, p. 113-120, 2009.

Jean M., « The non-smooth contact dynamics method », *Computer Methods in Applied Mechanics and Engineering*, vol. 177, p. 235-257, 1999.

Jean M., Acary V., Monerie Y., « Non-smooth contact dynamics approach of cohesive materials », *Phil. Trans. R. Soc. Lond.*, vol. A359, p. 2497-2518, 2001.

LMGC90, Web site : http://www.lmgc.univ-montp2.fr/~dubois/LMGC90/, 2009.

Magoules F., Roux F., « Lagrangian formulation of domain decomposition methods : A unified theory », *Applied Mathematical Modelling*, vol. 30, p. 593-615, 2006.

Martin E., Peters P., Leguillon D., Quenisset J., « Conditions for matrix crack deflection at an interface in ceramic matrix composites », *Materials Science and Engineering*, vol. A250, p. 291-302, 1998.

Martin R. C., « The Liskov Substitution Principle », *C++ Report*, 1996.

Martin R. C., *Agile Software Development. Principles, Patterns and Practices*, Alan Apt Series, Prentice Hall, Pearson Education, Inc. Upper Saddle Reaver, New Jersey 07458., 2003.

Meyer B., *Object-Oriented Software Construction*, second edn, Prentice Hall PTR, 1997.

Michel J.-C., Suquet P., « An analytical and numerical study of the overall behaviour of metal-marix composites », *Modelling Simulation in Materials Science and Engineering*, vol. 2, p. 637-658, 1994.
Moreau J.-J., « Unilateral contact and dry friction in finite freedom dynamics », in J.-J. Moreau, P. Panagiotopoulos (eds), *Non Smooth Mechanics and Applications*, vol. 302 of CISM - *Courses and Lectures*, Springer, Vienna, p. 1-82, 1988.

Nineb S., Alart P., Dureisseix D., « Domain decomposition approach for non-smooth discrete problems, example of a tensegrity structure », *Computers and Structures*, vol. 85, p. 499-511, 2007.

Ortiz M., Stainier L., « The variational formulation of viscoplastic constitutive updates », *Computer methods in applied mechanics and engineering*, vol. 171, p. 419-444, 1999.

PELICANS, Web site : https://gforge.irsn.fr/gf/project/pelicans/, 2009.

Pelissou C., Baccou J., Monerie Y., Perales F., « Determination of the size of the representative volume element for random quasi-bittle composites », *International Journal of Solids and Structures*, 2009. In Press.

Perales F., Bourgeois S., Chrysochoos A., Monerie Y., « Two field multibody method for periodic homogenization in fracture mechanics of nonlinear heterogeneous materials », *Engineering Fracture Mechanics*, vol. 75, p. 3378-3398, 2008.

Perales F., Monerie Y., Chrysochoos A., « Nonsmooth fracture dynamics of functionally graded materials », in J. Cirne, R. Dormeval, et al. (eds), *8th International Conference on Mechanical and Physical Behaviour of Materials under Dynamic Loading*, vol. 134, EURODYMAT, Journal de Physique IV, Dijon, France, p. 367-372, sep, 2006.

Raous M., Cangémi L., Cocu M., « A consistent model coupling adhesion, friction and unilateral contact », *Computer Methods in Applied Mechanics and Engineering*, vol. 177, n° 3-4, p. 383-399, 1999.

Raous M., Monerie Y., « Unilateral contact, friction and adhesion : 3D crack in composite material », in J. Martins, M. Monteiro Marques (eds), *3d Contact Mechanics International Symposium*, Collection Solid Mechanics and its Applications, Kluwer, Peniche (Portugal), p. 333-346, 17-21 June, 2001.

Raous M., Monerie Y., « Unilateral contact, friction and adhesion : 3D cracks in composite materials », in J. Martins, M. M. Marques (eds), *Collection Solid Mechanics and its Applications, Contact Mechanics*, Kluwer, vol. 302 of CISM - *Courses and Lectures*, p. 333-346, 2002.

Renouf M., Optimisation numérique et calcul parallèle pour l'étude des milieux divisés bi- et tridimensionnels, PhD thesis, Université Montpellier II, 2004.

Renouf M., Bonamy D., Dubois F., Alart P., « Steady surface flows in rotating drum : Numerical Simulation », *Physics of Fluids*, vol. 17, p. 103303, 2005.

Siegmund T., Fleck N. A., Needleman A., « Dynamic crack growth across an interface », *International Journal of Fracture*, vol. 85, p. 381-402, 1997.

Stainier L., Dubois F., Peyroux R., « MatLib, une bibliothèque portable de modèles constitutifs pour la mécanique non-linéaire des solides : concepts et implémentation », in M. Potier-Ferry, M. Bonnet, A. Bignonnet (eds), *6ème Colloque National en Calcul des Structures*, vol. 3, Giens, p. 25-32, 2003.

Tvergaard V., « Effect of fibre debonding in a whisker-reinforced metal. », *Materials Science and Engineering*, vol. A125, p. 203-213, 1990.

Vincent P.-G., Modélisation micromécanique de la croissance et de la percolation de pores sous pression dans une matrice céramique à haute température, PhD thesis, Université d'Aix Marseille 1 - Provence, 2007.
Vincent P.-G., Monerie Y., Suquet P., « A ductile damage model for bi-porous saturated materials at high strain rate », *Multi-phase and multi-components materials under dynamic loading*, EMMC10, Poland, 2007.

Xu G., Bower A. F., Ortiz M., « The influence of crack trapping on the toughness of fiber reinforced composites. », *Journal of the Mechanics and Physics of Solids*, vol. 46, n° 10, p. 1815-1833, 1998.

Xu X.-P., Needleman A., « Numerical simulations of fast crack growth in brittle solids », *Journal of Mechanics and Physics of Solids*, vol. 42, n° 9, p. 1397-1434, 1994.

Yamanaka S., Yamada K., Kurosaki K., Uno M., Takeda K., Anada H., Matsuda T., Kobayashi S., « Thermal properties of zirconium hydride. », *Journal of Nuclear Materials*, vol. 294, p. 94-98, 2001.

Yamanaka S., Yoshioka K., Uno M., Katsura M., Anada H., Matsuda T., Kobayashi S., « Thermal and mechanical properties of zirconium hydride. », *Journal of Alloys and Compounds*, vol. 293-295, p. 23-29, 1999.

Yang Q., Stainier L., Ortiz M., « A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids », *Journal of Mechanics and Physics of Solids*, vol. 54, p. 401-424, 2006.

Zavattieri P. D., Raghuram P. V., Espinosa H. D., « A computational model of ceramic microstructures subjected to multi-axial dynamic loading. », *Journal of the Mechanics and Physics of Solids*, vol. 49, p. 27-68, 2001.