# Simulation numérique d'impact en dynamique rapide de matériaux hyper-élastiques par la méthode X-FEM

Patrick Rozycki, Nicolas Moes, Steven Le Corre, Malek Zarroug

► To cite this version:

Patrick Rozycki, Nicolas Moes, Steven Le Corre, Malek Zarroug. Simulation numérique d'impact en dynamique rapide de matériaux hyper-élastiques par la méthode X-FEM. 10e colloque national en calcul des structures, May 2011, Giens, France. hal-00592731

HAL Id: hal-00592731

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

Submitted on 3 May 2011

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
# Simulation numérique d'impact en dynamique rapide de matériaux hyper-élastiques par la méthode X-FEM

P. Rozycki<fnref n="1" />, C. Dubois<fnref n="2" />, N. Moës<fnref n="1" />, S. Le Corre<fnref n="1" />, M. Zarroug<fnref n="3" />

**Résumé** — La simulation numérique en dynamique rapide de structures composées de matériaux hyper-élastiques peut s'avérer fastidieuse surtout si l'on souhaite obtenir une bonne qualité de résultats. La distorsion excessive des éléments finis du maillage au cours du temps, qui conduit à une dégradation du pas de temps critique, impose par exemple l'utilisation de techniques de type ALE. Nous proposons dans cet article une méthode alternative consistant à utiliser un maillage régulier unique via la méthode X-FEM et dont le principal atout est de simplifier les procédures de remaillage.

**Mots clefs** — X-FEM, Dynamique explicite, Elasticité non-linéaire.

## 1 Introduction

Bien que les intérêts et avantages que possède la simulation numérique par éléments finis dans de nombreux domaines de la mécanique ne soit plus à démontrer, il n'en demeure pas moins que dans certains cas particuliers, des limitations peuvent brider son application. Citons par exemple le comportement de structures qui incluent des matériaux de type mousse lors d'impact ou de choc. En effet, lors des simulations impliquant de grandes transformations, les éléments finis composant la structure peuvent subir des distorsions importantes, voir leur jacobien devenir négatif imposant ainsi un arrêt brutal du calcul [1]... A cela peut s'ajouter le fait que, lorsque ces simulations sont réalisées à l'aide de schémas d'intégration explicites, le pas de temps chute éventuellement de manière drastique, accroissant ainsi considérablement les temps de calcul [2]. Des méthodes de type ALE (Arbitrary Lagrangian Eulerian) peuvent cependant apporter une solution à ces problèmes induits par l'utilisation d'un maillage lagrangien unique au cours du calcul explicite. L'avantage principal qui consiste à juxtaposer une description eulérienne en plus de la description lagrangienne est d'éliminer la dépendance entre l'évolution du maillage en regard du matériau offrant ainsi une meilleure qualité de maille. Cependant, comme il a été dit, dans le cadre de calculs explicites, cela a pour désavantage de diminuer fortement le pas de temps critique et il peut aussi s'avérer difficile de choisir la bonne évolution du maillage [1]. Par ailleurs, la mise en œuvre peut être parfois aussi, dans le cas de problèmes industriels, très fastidieuse.

Depuis quelques années la méthode des éléments finis étendues (X-FEM) a su s'imposer comme une solution très efficace pour tout problème incluant des discontinuités telles que des fissures, des trous, des inhomogénéités matérielles [3], [4]... Dans ce cadre, le maillage n'est plus tenu d'être conforme aux discontinuités à traiter et notamment dans le cas qui nous intéresse : les surfaces libres. L'étape de construction de nouveaux maillages si besoin est, devient alors triviale puisque l'on peut choisir rapidement des maillages structurés réguliers. Bien qu'il existe des couplages qui ont aussi été proposés entre X-FEM et des méthodes de type ALE dans le cadre des bi-m matériaux [5], [6], nous proposons dans cet article une application différente. Les travaux présentés [7], [8] se basent sur l'idée de remaillage à des instants particuliers pour lesquels le maillage courant devient de mauvaise qualité. Grâce à un algorithme spécifique de transfert de données, de nouveaux maillages viennent se substituer aux anciens. On choisit de travailler sur des maillages réguliers pour préserver les pas de temps critiques et donc les temps de calcul. Du fait de l'utilisation de la méthode X-FEM, les procédures peuvent être automatisées et ne requierent donc plus d'intervention « utilisateur ».

<footnote n="1">¹GeM, Ecole Centrale de Nantes, France, {patrick.rozycki,nicolas.moes,steven.le-corre}@ec-nantes.fr</footnote>

<footnote n="3">³AREVA, Paris la Défense, France, celine.dubois@areva.com</footnote>

<footnote n="2">²DRIA/SARA/STEO/MAST, PSA-Peugeot Citroën, France, malek.zarroug@mpsa.com</footnote>
# 2 X-FEM en dynamique rapide

## 2.1 X-FEM dans le cas de surfaces complexes non maillées

Comme il a été dit précédemment la méthode X-FEM permet de gérer de nombreux types de discontinuités. Dans notre cadre d'étude, nous nous restreignons à un type en particulier pour lequel il n'est pas nécessaire d'avoir un enrichissement à proprement parler : la discontinuité engendrée par des surfaces complexes non maillées, comme par exemple des trous de « vide » dans la matière ou des surfaces extérieures libres de structure. Les éléments traversés par la frontière de ces surfaces (décrites par des Level Set) sont donc constitués d'une partie de « matière » et de « vide » (Figure 1).

<figure><img src="image_2.png" /><figcaption>FIG. 1 – Exemple de surface non maillée et d'éléments finis X-FEM</figcaption></figure>

Les fonctions d'interpolation qui sont utilisées, sont les fonctions d'interpolation classiques $\varphi_i(x)$ d'éléments finis multipliées par une fonction d'enrichissement $V(x)$ qui prend la valeur de 1 si le nœud est dans le « vide » et 0 dans le cas contraire [9], [10] :

$$u(x) = \sum_{i} \varphi_{i}(x)V(x)u_{i} = \sum_{i} \hat{\varphi}_{i}(x)u_{i}$$

Ce choix revient à intégrer les fonctions sur le domaine « matière » uniquement ; dans le cas d'éléments finis linéaires, on retrouve que la matrice raideur de l'élément X-FEM est proportionnelle à celle d'un élément « plein » d'un rapport équivalent à la fraction de matière. Par contre, pour la matrice masse consistante et la matrice masse diagonalisée (utilisée en dynamique rapide), il n'en va pas de même [11].

## 2.2 Technique de diagonalisation de matrices de masse et pas de temps critique

Dans les codes de dynamique explicites, pour des raisons principalement de gain de temps de calcul et de mise en œuvre, des techniques de « lumping » sont utilisées pour diagonaliser la matrice des masses ; son inversion par la suite en devient triviale. Pour l'élément fini X-FEM développé nous avons dû repenser à la technique de diagonalisation employée : cette dernière ne permettait pas de découpler le pas de temps critique induit (pour respecter la condition CFL) de la fraction matière $\varepsilon$. Dans le cadre d'un élément fini barre X-FEM, nous avons :

$$K_{\text{élémentaire}} = \frac{\varepsilon E S}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}; M_{\text{élémentaire}} = \frac{\varepsilon \rho S L}{2} \begin{bmatrix} 2 - \varepsilon & 0 \\ 0 & \varepsilon \end{bmatrix} \Rightarrow \Delta t_c \le \frac{2}{\omega_{\max}} = \frac{L}{c} \sqrt{\varepsilon(2 - \varepsilon)}$$

L'idée a donc été de répartir la masse de façon équitable sur chacun des nœuds [11] ; dans ce cas de diagonalisation, la matrice masse retrouve la propriété de proportionnalité qu'avait la matrice raideur. De ce fait, la pulsation propre du système devient indépendante de la fraction matière :
$$K_{\text{élémentaire}} = \frac{\varepsilon E S}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}; \quad M_{\text{élémentaire}}^{\text{lumped}} = \frac{\varepsilon \rho S L}{2} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \Rightarrow \Delta t_c \le \frac{2}{\omega_{\text{max}}} = \frac{L}{c}$$

L'utilisation d'éléments X-FEM avec la technique de diagonalisation particulière permet, en dynamique rapide, de conserver un pas de temps identique à celui d'un élément fini « plein ». La méthode devient, de ce fait, avantageuse lors de l'emploi de schémas explicites, et donc compétitive.

Il est à noter qu'il existe d'autres approches, plus adaptées à la propagation de fissures, concernant la méthode X-FEM et les schémas d'intégration en temps [12].

## 2.3 Pas de temps critique et problème dynamique non linéaire

Notre propos se réfère à l'impact de structures en mousse, matériaux cellulaires de faible densité. Bien qu'il existe des descriptions appropriées du comportement de ce type de matériau, dans un souci de simplification, nous utilisons une loi hyper élastique classique Néo-Hookénne compressible. Pour les problèmes non linéaires, l'évaluation du pas de temps critique est légèrement différente du cas linéaire : il est donné par le système linéarisé (1) que l'on vient ensuite pondérer par un coefficient compris généralement entre 0.8 et 0.9 [13].

$$K_{\text{élémentaire}}^{\text{tangent}} - \omega^2 M_{\text{élémentaire}}^{\text{lumped}} \quad (1)$$

La matrice raideur tangente est obtenue par dérivation du vecteur des efforts intérieurs par rapport au champ de déplacement. Ils sont donnés pour le noeud $i$ par (2) où $P$ est le premier tenseur de Piola-Kirchhoff et $\hat{\phi}_i(x)$ les fonctions de forme X-FEM.

$$F_i^{\text{intérieur}} = \int_{A_{XFEM}} \frac{\partial \hat{\phi}_i}{\partial x} P dA \quad (2)$$

On constate que l'évaluation de ces efforts intérieurs n'implique pas uniquement que la dérivée des fonctions d'interpolation ; comme pour le cas linéaire, il existe une relation de proportionnalité à la fraction de matière pour les efforts intérieurs et donc pour la matrice raideur tangente. Cela nous permet aussi en utilisant la technique de diagonalisation précédente de conserver un pas de temps critique identique à celui qui est obtenu pour des éléments finis classiques, remplis de matière.

## 3 Remaillage et projection

La Figure 2 présente le problème de remaillage dans le contexte X-FEM. Le maillage bleu (M1) correspond à un premier calcul X-FEM lagrangien présenté dans sa configuration déformée à l'instant choisi pour le remaillage. Comme il s'agit d'un calcul X-FEM, il est à noter que les bords de la structure déformée ne sont pas conformes avec les éléments du maillage. On souhaite poursuivre le calcul précédent sur un nouveau maillage (M2), régulier et représenté en rouge sur la Figure 2.

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Exemple de problème de remaillage dans le contexte X-FEM</figcaption></figure>
## 3.1 Méthode de projection des quantités

Parmi les classes de méthodes qui existent et puisque nous nous intéressons aux matériaux hyper élastiques pour lesquels il n'y a pas de variables internes à transporter, le choix s'est porté sur une méthode par collocation (méthode la plus courante). Elle se décompose en trois étapes principales : l'identification de l'emplacement de chacun des nouveaux nœuds (M2) dans l'ancien maillage (M1), la détermination des nouvelles coordonnées locales de ce nouveau nœud (M2) dans l'ancien élément fini (M1) et enfin l'utilisation de l'interpolation sur l'ancien maillage (M1) pour calculer les nouvelles valeurs nodales du nœud (M2).

L'une des problématiques en dynamique explicite est que le transfert devrait inclure non seulement les champs de déplacement mais aussi de vitesse et d'accélération. Il est donc difficile d'assurer la conservation des quantités de mouvement sur le nouveau maillage (M2). Suite aux travaux dans [14], nous avons choisi de projeter uniquement les prédicteurs des champs de déplacement et de vitesse de l'ancien maillage (M1) vers le nouveau (M2). Une étape de rééquilibrage est ensuite effectuée en résolvant les équations du mouvement sur le nouveau maillage et se termine en mettant à jour partiellement le champ de vitesse. Cette procédure permet d'assurer la stabilité de la projection et de garantir l'équilibre énergétique entre les instants. Il subsiste cependant des apports ou pertes énergétiques inévitables lors de l'étape de projection qui ne sont pas estimables *a priori*. Il est nécessaire de quantifier leur importance lors des simulations.

## 3.2 Scan conversion

La méthode de projection nécessite une étape d'identification des nouveaux nœuds (M2) sur l'ancien maillage (M1). Comme nous utilisons la méthode X-FEM et que nous souhaitons l'appliquer sur des grilles de maillage régulières, l'algorithme de recherche qui nous est apparu le plus efficient au regard de ces contraintes, est un algorithme issu de l'informatique « géométrique » : le scan conversion [15]. Le but d'un scan conversion est de déterminer efficacement les pixels situés l'intérieur d'un polygone, c'est-à-dire sans avoir testé l'ensemble de la grille de pixels (Figure 3).

<figure><img src="image_4.png" /><figcaption>FIG. 3 – Déroulement de l'algorithme de scan conversion d'un polygone en 2D</figcaption></figure>

Pour des maillages uniformes 2D la technique peut donc être utilisée directement. Dans le cadre de maillages uniformes 3D utilisant des tétraèdres, il est nécessaire de définir avant l'intersection d'un plan « mobile » avec le tétraèdre considéré : une fois cette étape effectuée, on se retrouve dans le cas précédent de polygone (Figure 4).
<figure><img src="image_5.png" /><figcaption>FIG. 4 – Technique de scan d'un tétraèdre</figcaption></figure>

Dans notre cas d'étude du remaillage sur maillage régulier, le nouveau maillage (M2) est assimilé à la grille de pixels. L'opération de scan conversion est réalisée sur tous les éléments tétraédriques du maillage (M1) en configuration déformée. En ce qui concerne les éléments de l'ancien maillage (M1) traversés par la Level Set, un scan conversion de leur partie « matière » et de leur partie « vide » est fait séparément : il sera possible ainsi de distinguer les nœuds du nouveau maillage (M2) contenus dans la matière de ceux contenus dans le vide.

## 3.2 Construction de la Level Set sur le nouveau maillage

La Level Set, contrairement aux champs de déplacement, n'est pas un champ lagrangien mais eulérien. Lorsque le maillage se déplace, la Level Set stockée à chaque nœud représente la distance du nœud à l'interface par rapport à sa configuration initiale et non par rapport à la configuration déformée. Par conséquent, la Level Set ne devient plus une fonction distance lors de la déformation du maillage initial. Le transfert par collocation est donc impossible car cela induirait une position erronée de l'interface (Figure 5).

<figure><img src="image_6.png" /><figcaption>FIG. 5 – Exemple d'un transport par collocation de la Level Set</figcaption></figure>

Afin de limiter les opérations et de bénéficier pleinement des informations obtenues lors du scan conversion, il est apparu plus judicieux de discrétiser l'interface et de calculer ensuite la distance afin de reconstruire la Level Set sur le nouveau maillage (M2). Un algorithme a été proposé afin de sélectionner les éléments proches de l'interface de l'ancien maillage (M1). Via les informations obtenues par le scan conversion, on sélectionne les nœuds du nouveau maillage (M2) proches de l'interface. Via les informations obtenues par le découpage du maillage initial (M1) par la Level Set, on limite le nombre de facettes. Le scan conversion réalisé séparément pour la Level Set en fonction des domaines « matière » et de « vide » permettra de définir enfin le signe à affecter à la Level Set après calcul de la distance. L'application de l'algorithme sur des cas « classiques » et des cas
« pathologiques » a permis de mettre en avant une robustesse des développements menés [8].

# 4 Implémentation d'une gestion de contact par pénalités

La gestion des contacts se fait dans notre cas par la méthode des pénalités. En effet, c'est le choix qui s'impose puisque les interfaces ne sont pas maillées (représentation par Level Set). L'algorithme de la méthode de contact va nous imposer une évaluation d'efforts de contact si une pénétration est détectée. Suite à cela, une réévaluation des champs de déplacement et de vitesse est menée. La pénalisation consiste à augmenter la raideur des éléments de contact. Dans le cadre de notre élément fini XFEM, cela impose donc, pour satisfaire les conditions de stabilité, un pas de temps critique adéquat. En raison de la pénalisation, ce dernier devient fonction de la fraction « matière » et tendra vers zéro si la fraction « matière » tend vers 0. Pour contourner ce problème, nous avons adopté une technique de mass scaling [8] : en ajoutant une quantité de masse équivalente à ce qui est perdu par la fraction de « vide », on montre que le pas de temps nécessaire tend vers une limite finie, cette fois-ci, quand la fraction « matière » tend vers zéro. En comparaison avec des éléments finis classiques, sur une étude monodimensionnelle on constate que le pas de temps critique dans notre cas X-FEM est toujours supérieur au pas de temps critique « standard » (Figure 6). Par conséquent, c'est un grand atout, et ce, même si nous utilisons une technique de mass scaling (la masse ajoutée est bornée et ne changera pas de manière importante la masse totale du système étudié).

<figure><img src="image_7.png" /><figcaption>FIG. 6 – Comparaison du rapport des pas de temps FEM et X-FEM pour différentes pénalités</figcaption></figure>

# 5 Exemple d'application

Nous considérons un exemple 3D impliquant un remaillage au cours de l'impact [8] ; cet exemple concerne l'impact d'une sphère sur un plan rigide. La sphère, de rayon $0,3 \text{ m}$, a une vitesse initiale de $10 \text{ m/s}$. Un maillage tétraédrique régulier isotrope dont la taille de maille est de $0,05 \text{ m}$ est utilisé. La géométrie de la sphère est décrite par une fonction Level Set et le plan rigide n'est pas discrétisé mais simplement représenté par une fonction analytique. L'écart initial entre la sphère et le plan est de $0,005 \text{ m}$. Durant l'impact, la sphère déformée est remaillée par notre méthode de projection avec un maillage identique au premier (Figure 7). Les simulations sont effectuées à la fois sur un code de calcul explicite commercial et sur notre code X-FEM.

Les points principaux que l'on peut observer sont, d'une part que, les choix que nous avons effectués tant sur les plans éléments finis, techniques de diagonalisation, remaillage et gestion du contact permettent de représenter correctement le comportement de la sphère en impact au regard d'un calcul éléments finis classiques (c.f. calcul des énergies de la Figure 8). D'autre part, il est à noter que le remaillage proposé permet de mieux contrôler le pas de temps critique de la simulation numérique. En effet, plutôt que de continuer à décroître sur un maillage unique pour respecter la condition de stabilité, le remaillage permet de conserver un pas de temps critique plus stable. Nous ne le présentons
pas dans cet article, mais il en va de même pour des études portant sur les aspects (aspect ratio) des éléments finis.

<figure><img src="image_8.png" /><figcaption>FIG. 6 – Maillages utilisés au cours du calcul</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 7 – Evolution des énergies et du pas de temps critique pour l’exemple de la sphère</figcaption></figure>

# 6 Conclusions et perspectives

Dans ces travaux, nous avons montré que nous pouvions aborder sous un autre angle les problèmes d'impact sur structures de type mousse. L'idée d'utiliser les avantages de la méthode X-FEM nous permette de simplifier grandement les problèmes liés au remaillage et ce, notamment en utilisant des maillages réguliers. De ce fait, les procédures de recherche d'appartenance de nœuds entre maillages pour la méthode de collocation sont optimisées. Par ailleurs, la méthode que nous avons développée pour le transport de la Level Set s'est montrée aussi efficace. Par ailleurs, nous avons aussi implémenté une méthode de gestion de contact par pénalités dans le cadre X-FEM qui a su démontrer ses atouts. L'ensemble des développements ont été testé sur plusieurs exemples et ont montré des résultats plus que satisfaisants tant sur la qualité des calculs que sur les temps de simulation. Il est possible avec la stratégie de remaillage adoptée de conserver un aspect géométrique d'éléments et un pas de temps raisonnable.

L'objectif des travaux étaient à but industriel. De ce fait, ils se poursuivent maintenant vers l'optimisation de cette méthode pour des cas non « académiques ». Pour ce faire, une des voies qui est explorées maintenant est l'utilisation de maillage Octree ; l'avantage du classement des entités sur ce type de maillage permettrait de connaitre encore plus rapidement les informations au moment du transfert. La gestion des étapes de remaillage et de raffinement autour de certains points spécifiques dans les maillages s'effectue aussi avec un temps amoindri. L'introduction et l'utilisation de ce genre de maillage nous autoriseraient à penser que les étapes relatives aux remaillages et au transport des
champs pourraient être plus rapides et automatiques.

Enfin, nous nous tournons aussi vers une autre réécriture des schémas numériques d’intégration explicites afin d’explorer l’optimisation en temps de calcul. Pour des structures de taille importante, une méthode classique parmi d’autres serait d’utiliser des sous-domaines, par exemple pour diminuer les temps de calcul. Cependant, dans ce genre de méthode, le recollage peut utiliser des multiplicateurs de Lagrange, ce qui peut être pénalisant. La méthode proposée dans [16] permet de partitionner la structure et, grâce à la définition de certains filtres, de n’effectuer que certaines opérations en prédiction, d’autres en correction et ceci sans introduction de données supplémentaires telles les multiplicateurs de Lagrange.

Références

[1] D.J. Benson, N. Moës, T. Belytschko. *Computational methods in lagrangian and eulerian hydrocodes*, Computer Methods in Applied Mechanics and Engineering, Elsevier, page235- page394, 1992.

[2] T. Belytschko, *W. Liu and B. Moran*. Nonlinear finite elements for continua and structures, Wiley, 2000.

[3] C. Daux, N. Moës, J. Dolbow, N. Sukumar, T. Belytschko. *Arbitrary branched and intersecting cracks with the extented finite element method*, International Journal of Numerical Methods in Engineering, Wiley, page1741- page1760, 2000.

[4] N. Moës, J. Dolbow, T. Belytschko. *A finite element method for crack growth without remeshing*, International Journal of Numerical Methods in Engineering, Wiley, page131- page150, 1999.

[5] AR. Khoei, M. Anahid, K. Shahim. *An extended arbitrary Lagrangian-Eulerian finite element method for large deformation of solids mechanics*, Finite Element Analysis Design, Elsevier, page401- page416, 2008.

[6] E. Vitali, D.J. Benson, T. Belytschko. *An extended finite element formulation for contact in multi-material arbitrary Lagrangian-Eulerian calculations*. , International Journal of Numerical Methods in Engineering, Wiley, page1420- page1444, 2006.

[7] C. Dubois, S. Le Corre, M. Zarroug, P. Rozycki, N. Moës. *Impact on highly compressible media in explicit dynamics using the X-FEM*, Computational Mechanics, Springer, page 329- page 348, 2006.

[8] C. Dubois, *Application de la méthode des éléments finis étendus à la simulation des grandes transformations en dynamique rapide*, thèse de doctorat, Ecole Centrale de Nantes, 2010.

[9] N. Sukumar, N. Moës, T. Belytschko. *Modelling holes and inclusions by level sets in the extended finite element method*, Computer Methods in Applied Mechanics and Engineering, Elsevier, page6183-page6200, 2001.

[10]N. Moës, M. Cloirec, P. Cartraud, J.F. Remacle. *A computational approach to handle complex microstructure geometries*, Computer Methods in Applied Mechanics and Engineering, Elsevier, page3163-page3177, 2003.

[11]P. Rozycki, N. Moës, E. Bechet, C. Dubois. *X-FEM explicit dynamics for constant strain elements to alleviate mesh constraints on internal or external boundaries*, Computer Methods in Applied Mechanics and Engineering, Elsevier, page349-page363, 2008.

[12]T. Menouillard, J. Rethore, N. Moës, A. Combescure, H. Bung. *Mass lumping strategies for X-FEM explicit dynamics: application to crack propagation*, International Journal of Numerical Methods in Engineering, Wiley, page447-page474, 2008.

[13]T. Belytschko, *W. Liu and B. Moran*. Nonlinear finite elements for continua and structures, Wiley, 2000.

[14]J. Rethore, A. Gravouil, A. Combescure. *A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing*, Computer Methods in Applied Mechanics and Engineering, Elsevier, page4493-page4510, 2004.

[15]D Breen, S. Mauch, R. Whitaker. *3D scan conversion of csg models into distance, closest-point and color volume*, Proceedings of Volume Graphics, Springer, page135-page158, 2000.

[16]F Casadei, JP. Halleux. *Binary spatial partitioning of the central-difference time integration scheme for explicit fast transient dynamics*, International Journal of Numerical Methods in Engineering, Elsevier, page1436-page1476, 2008.