# Comment gérer en dynamique non régulière des contacts l'indétermination et l' interpénétration dans les milieux granulaires ?

Pierre Alart

► To cite this version:

Pierre Alart. Comment gérer en dynamique non régulière des contacts l'indétermination et l' interpénétration dans les milieux granulaires ?. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722060

HAL Id: hal-01722060

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

Submitted on 2 Mar 2018

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 CC0 1.0 - Universal - International License
# Comment gérer en dynamique non régulière des contacts l'indétermination et l'interpénétration dans les milieux granulaires?

Pierre Alart<fnref n="1" />

**Résumé** — On présente une étude exploratoire visant à maîtriser la convergence du schéma d'intégration de l'approche NSCD (NonSmooth Contact Dynamics) notamment lors de la simulation de granulats denses fortement confinés. Deux thèmes sont traités : l'indétermination du réseau d'efforts de contact et l'interpénétration résiduelle due à la formulation en vitesse-impulsion du contact unilatéral. Au-delà de l'optimisation numérique visée, l'approche proposée peut motiver une stratégie de représentation multi-modèle séparant les comportements global-statique-continu et local-dynamique-discret.

**Mots clés** — Mécanique non régulière, Dynamique des contacts, Milieux granulaires, Indétermination

## 1 Introduction

L'algorithme de Gauss-Seidel Non Linéaire (Non Linear Gauss-Seidel - NLGS) est le solveur générique associé à la Dynamique Non Régulière des Contacts (Non Smooth Contact Dynamics - NSCD) [1]. L'intérêt principal de la formulation NSCD est de permettre de simuler le comportement d'une collection de corps avec des régimes variés pouvant coexister au sein d'un même système : statique (solide), dynamique lente, dynamique rapide (type fluide). En effet il est possible de résoudre, sur un pas de temps, de nombreux contacts simultanés. Le pas de temps peut être beaucoup plus grand que pour la Dynamique Moléculaire qui requiert un pas de temps d'autant plus petit que le système est de grande taille et dense.

Pourtant la méthode présente deux inconvénients. Le premier est lié à la formulation en vitesse de la condition de non-interpénétration. Cette formulation est exacte en temps continu, mais induit des interpénétrations parasites dans un schéma à pas de temps fixe (time-stepping integration). Le second inconvénient réside dans la convergence lente d'un algorithme de type Gauss-Seidel. Pour un granulat dense en régime quasi-statique la méthode peut même diverger, ce qui permet de relier cette défaillance à l'indétermination caractérisée par de multiples réseaux de forces auto-équilibrés. Les deux inconvénients se cumulent quand une granulat est fortement confiné par des conditions aux limites en vitesses imposées.

Nous proposons ici une étude exploratoire d'une stratégie numérique visant à accélérer la convergence en tenant compte de l'indétermination et en contrôlant l'interpénétration dans un algorithme de prédiction-correction (prédiction linéaire, correction non régulière modifiée). Une telle approche constitue une piste intéressante pour séparer les contributions globale-statique et locale-dynamique dans les contraintes internes et motiver par exemple une stratégie multi-modèle (discret-local vs continu-global) dans une méthode de décomposition de domaine [2, 3]. Les quelques développements numériques ont été effectués sur des échantillons granulaires de taille réduite (de 5 à 500 grains) sous les logiciels MAPLE et MATLAB.

## 2 Dynamique non régulière

En suivant un schéma d'intégration pas à pas (time-stepping) les variables sur un pas de temps sont présentées dans le Tableau 1. Partant des paramètres de configuration x, on dérive les vitesses généralisées V. Dans les granulats confinés il est utile de distinguer les vitesses libres $V_f$, inconnues, et les vitesses données $V_g$. La vitesse relative au contact $\alpha$ est fournie par une relation linéaire des vitesses généralisées $v_\alpha = H_\alpha^T V$. Par dualité l'impulsion généralisée $R_\alpha$ associée à l'impulsion locale $r_\alpha$ satisfait la relation $R_\alpha = H_\alpha r_\alpha$. Le schéma d'intégration de nature implicite est présenté dans Algorithme 1.

<footnote n="1">¹ LMGC – UMR 5508, Université Montpellier 2 / CNRS, Pierre.Alart@univ-montp2.fr</footnote>
<table><caption>Tableau 1 – Notations générales</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>[t<sup>-</sup>, t<sup>+</sup>]</td><td>Intervalle de temps (dt = t<sup>+</sup> - t<sup>-</sup>, pas de temps)</td></tr><tr><td>x<sup>-</sup> (resp. x for x<sup>+</sup>)</td><td>Coordonnées généralisées à t<sup>-</sup> (resp. t<sup>+</sup>)</td></tr><tr><td>V<sup>-</sup> (resp. V for V<sup>+</sup>)</td><td>Vitesses généralisées à t<sup>-</sup> (resp. t<sup>+</sup>)</td></tr><tr><td>V<sub>f</sub><sup>-</sup> (resp. V<sub>g</sub><sup>-</sup>)</td><td>Vitesses libres (resp. données) à t<sup>-</sup></td></tr><tr><td>R = H<sub>f</sub> r</td><td>Torseurs d'impulsions des contacts par grain</td></tr><tr><td>g (for g(x))</td><td>Interstices de contact</td></tr><tr><td>v (for v<sup>+</sup>) : v = H<sup>T</sup> V = H<sub>f</sub><sup>T</sup> V<sub>f</sub> + H<sub>g</sub><sup>T</sup> V<sub>g</sub></td><td>Vitesses relatives aux contacts</td></tr><tr><td>r</td><td>Impulsions de contact sur [t<sup>-</sup>, t<sup>+</sup>] par contact</td></tr><tr><td>M (for M<sub>f</sub>(x))</td><td>Matrice d'inertie</td></tr></tbody></table>

Le système d'équations et d'inéquations (1) est résolu par un algorithme de type Gauss-Seidel Non Linéaire détaillé dans Algorithme 2. D'inconnues principales $v$ et $r$, le système comporte deux lignes : la dynamique réduite aux contacts via l'opérateur de Delassus ($W = H_f^T M^{-1} H_f r$) et une équation représentant formellement l'ensemble des relations de contact, de frottement, voire de comportement aux contacts plus complexes (adhésion, décohésion, pont capillaire,...). Dans un granulat la matrice $W$ est singulière avec un noyau d'autant plus grand que le confinement est fort.

$$ \left\{ \begin{array}{l} W r - v = -v^d = -H_f^T V_f^- - H_g^T V_g^+ - H_f^T M^{-1} R^d \\ \mathcal{R}(v,r) = 0 \end{array} \right. \qquad (1) $$

Parmi les interactions, le contact unilatéral doit être réécrit en vitesse, alors qu'il est naturellement exprimé par des relations de complémentarité entre interstice, noté $g_\alpha$, et force de contact. Le lemme de viabilité de Moreau [1] justifie la formulation d'une loi de Signorini vitesse-impulsion sous forme d'un algorithme sur la composante normale de la vitesse relative $v_{\alpha,n}$ et de l'impulsion de contact $r_{\alpha,n}$ pour chaque contact $\alpha$ :

$$ \begin{array}{l} \text{Si } \tilde{g}_{\alpha} \ge 0, \text{ alors } r_{\alpha,n} = 0 \\ \text{Si } \tilde{g}_{\alpha} < 0, \text{ alors } 0 \le v_{\alpha,n} \perp r_{\alpha,n} \ge 0 \end{array} \qquad (2) $$

En temps continu le lemme de viabilité garantit l'équivalence de la formulation en vitesse et de la formulation en interstice. En temps discret la notation $\tilde{g}_\alpha$ constitue une prédiction de l'interstice et non une évaluation exacte. Plusieurs stratégies sont envisageables sans conséquence sur les développements ultérieurs ; ici, la valeur en fin du pas précédent est retenue, $\tilde{g}_\alpha = g_\alpha$. En temps discret les interpénétrations peuvent s'accumuler, et ce d'autant plus que l'Algorithme 2 peine à converger. C'est notamment le cas quand la matrice $W$ est singulière, avec un grand espace de réseaux de forces de contact auto-équilibrés. **L'indétermination du réseau d'efforts de contact est donc un facteur essentiel dans l'accumulation des interpénétrations comme dans le ralentissement ou de mise en échec du solveur.**

On peut illustrer les difficultés algorithmiques sur de petits cas très simples simulés sous MAPLE en se focalisant sur le nombre d'itérations à convergence, limité à 150 dans les exemples présentés ci-après si l'algorithme de converge pas en deçà, et sur la violation moyenne, c'est à dire la moyenne de l'interpénétration sur l'ensemble des contacts. Les effets dus au coincement peuvent être mis en évidence sur le cas de 2 billes de taille identique serrées entre 2 parois. Mais sans aller jusqu'à des situations aussi pathologiques, considérons un cas similaire avec 5 billes (Figure 1). Une simulation sur 40 pas de temps montre la saturation fréquente du nombre d'itérations à 150 (Figure 2a) et l'accumulation de violation de l'interpénétration (Figure 2b), d'autant plus marquée quand l'algorithme peine à converger.

## 3 Lever l'indétermination

Pour lever l'indétermination, la première idée est de régulariser $W$ par son noyau quitte à introduire un multiplicateur de Lagrange pour contraindre le second membre à être orthogonal au noyau. On peut montrer cependant que l'on peut introduire une prédiction du réseau d'impulsions autorisant à ne pas itérer sur le multiplicateur conduisant ainsi à un schéma prédicteur-correcteur (Algorithme 3).
**Algorithm 1** Schéma d'intégration de l'approche NSCD

for *istep* = 1,...,*nstep* **do**

Détecter les contacts and

Calculer l'interstice prédit ($\tilde{g} = g^-$)

Définir les contacts potentiels $\alpha$ ($\tilde{g}_{\alpha} \le 0$) and

Définir les opérateurs $H, H_f, H_g$

Initialiser la configuration ($x^-, V^-, r^-$) au temps $t^-$

Calculer l'opérateur de Delassus $W$ et le membre droit $v^d$

Résoudre le système suivant par l'algorithme NLGS (Algorithm 2)

$$ egin{cases} W r - v = -v^d \\ \mathcal{R}(v, r) = 0 \end{cases} $$

Actualiser la configuration ($r^+ = r$)

$$ egin{cases} V_f^+ = V_f^- + M^{-1} (H_f r^+ + R^d) \\ x^+ = x^- + dt V^+ \end{cases} $$

end for

**Algorithm 2** Solveur NLGS

for *k* = 1,...,convergence **do**

for $\alpha = 1,...,n_{contact}$ **do** {boucle sur les contacts}

Déterminer ($r_{\alpha}^k, v_{\alpha}^k$) tel que

$$ egin{cases} W_{\alpha\alpha} r_{\alpha}^k - v_{\alpha}^k = -v_{\alpha}^d - \sum_{\beta<\alpha} W_{\alpha\beta} r_{\beta}^k - \sum_{\alpha<\beta} W_{\alpha\beta} r_{\beta}^{k-1} \\ \mathcal{R}(r_{\alpha}, v_{\alpha}) = 0 \end{cases} $$

end for

end for

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Confinement de 5 billes : configurations initiale et finale.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2 – Evolution du nombre d'itérations et de la violation sur 40 pas</figcaption></figure>

**Algorithm 3** Un algorithme de prediction-correction utilisant le noyau de W (NLKerGS-P)

Prédiction : $\tilde{r} = 0$ or $r^{\text{lin}}$ (Algorithm 4) ou $r^{\text{unil}}$ ou $r^{\text{elas}}$ (Algorithm 5)

Correction : trouver $(r, v)$, via le solveur NLGS, tel que

$$ \begin{cases} (W + NN^T) r - v = -v^d + NN^T \tilde{r} \\ \mathcal{R}(r, v) = 0 \end{cases} $$

Dans un premier temps le prédicteur linéaire est basé sur une régularisation de la seule loi de contact unilatéral, le frottement n'étant pas pris en compte dans la prédiction. Une fois le contact prédit, les relations de complémentarité (2) sont approchées par une relation linéaire (3) où $\eta$ est un petit paramètre.

$$ \begin{aligned} \text{si } g_{\alpha}^{-} \geq 0, \text{ alors } \tilde{r}_{\alpha,n} &= 0 \\ \text{si } g_{\alpha}^{-} < 0, \text{ alors } \tilde{r}_{\alpha,n} &= -\frac{dt}{\eta} \tilde{v}_{\alpha,n} \quad (\eta \to 0) \end{aligned} \quad (3) $$

L'implémentation de cette procédure passe par la sélection des contacts actifs via un opérateur diagonal $D^{-}$ permettant alors de résoudre un système linéaire d'inconnues les vitesses généralisées libres et d'en déduire une prédiction des impulsions de contact (Algorithm 4).

**Algorithm 4** Etape de prédiction linéaire

$$ \text{Définir la matrice } D^{-} : D^{-}_{\alpha\alpha,n} = \begin{cases} 0 & \text{if } g_{\alpha}^{-} \geq 0 \\ 1 & \text{if } g_{\alpha}^{-} < 0 \end{cases} $$
$$ \text{Résoudre le système linéaire : } (M + \frac{dt}{\eta} H_f D^{-} H_f^T) \tilde{V}_f = R^d + M V_f^{-} - \frac{dt}{\eta} H_f D^{-} H_g^T V_g^{+} $$
$$ \text{Définir l'impulsion prédite : } r^{\text{lin}} = -\frac{dt}{\eta} D^{-} H_f^T \tilde{V} $$

La réduction du nombre d'itérations (Figure 3a) est alors spectaculaire pour un résultat très comparable à celui obtenu avec le solveur générique. Le coût à payer est alors une "dérive de l'interpénétration"(Figure 3b). Cette dérive apparaît également avec le solveur générique mais de manière plus localisée et l'interpénétration moyenne a tendance à saturer, au moins tant que la sollicitation reste régulière (vitesse constante de la paroi droite). Avec notre approche cette dérive est plus diffuse à travers l'échantillon autorisant une approche globale pour la contrôler, mais l'accumulation est sensible et régulière. Un prédicteur "unilatéral" résolvant exactement les relations de contact unilatéral (2), mais sans gérer le frottement, a été testé (les résultats ne sont pas donnés ici). Il s'agit d'un solveur non linéaire de type NLGS (éventuellement gradient conjugué projeté), plus lourd à mettre en oeuvre, même s'il converge
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Evolution du nombre d’itérations et de la violation (5 disques) : algorithme générique, prédiction linéaire, prédiction élastique pour 2 valeurs du paramètre $\epsilon$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Configurations finales pour l’algorithme générique et 2 prédicitions (lineaire, élastique avec $\eta = 10^{-5}$)</figcaption></figure>

rapidement (sans frottement), sans paramètre de régularisation, mais n’améliorant pas significativement la dérive d’interpénétration due au prédicteur linéaire. Un prédicteur plus "riche", linéarisant le frottement ou imposant un contact adhérant par exemple, s’avère trop contraignant et amplifie encore la dérive d’interpénétration. Résultat au demeurant logique puisque c’est la composante contact unilatéral des lois d’interactions qui est approchée en vitesse en temps discret.

# 4 Contrôler l’interpénétration

Pour éviter une accumulation de l’interpénétration il est naturel d’introduire une dose d’élasticité dans les contacts (à défaut d’une solution plus élégante). L’étape de prédiction est alors modifiée (Algorithm 5).

<caption>**Algorithm 5** Etape de prédiction élastique</caption>

Définir la matrice $D^{-} : D_{\alpha\alpha,n}^{-} = \begin{cases} 0 & \text{if } g_{\alpha}^{-} \ge 0 \\ 1 & \text{if } g_{\alpha}^{-} < 0 \end{cases}$

Résoudre le système linéaire : $(M + \frac{dt^2}{\epsilon} H_f D^{-} H_f^T) \tilde{V}_f = R^d + M V_f^{-} - \frac{dt}{\epsilon} H_f D^{-} (g^{-} + dt H_g^T V_g^+)$

Définir l’impulsion prédite : $r^{\text{elas}} = -\frac{dt}{\epsilon} D^{-} (g^{-} + dt H^T \tilde{V})$
Le gain en nombre d'itérations est préservé (Figure 3a) et l'interpénétration est corrigée au cours du processus d'autant plus rapidement que le paramètre $\epsilon$ est petit, c'est à dire la raideur de contact plus grande (Figure 3b). On peut même obtenir sur la violation moyenne une amélioration de l'indicateur par rapport à la méthode générique, mais ceci est trompeur car si la violation moyenne diminue, des vibrations parasites constituées de micro-rebonds entre les grains engendrent des interpénétrations locales importantes et une éventuelle déstabilisation globale du système sans interprétation physique.

Un prédicteur élastique peut donc remédier à l'accumulation d'interpénétrations, mais au prix d'artefacts dynamiques associés à une élasticité de contact d'origine pour l'instant purement numérique. Alors même que les milieux granulaires sont sujets au phénomène dit d'effondrement inélastique [4] où la dissipation du système croît très rapidement avec le nombre de particules. L'élasticité n'est plus alors que résiduelle. On peut bien entendu envisager de combiner prédicteurs linéaire et élastique, mais on retrouvera les difficultés de calibrage discutées dans [5]. Il reste sur ce chapitre des questions ouvertes afin de choisir ce paramètre en fonction de la dynamique globale du processus simulé.

# 5 Optimiser la mise en œuvre numérique

La prise en compte du noyau $N$ de $W$ exige son calcul, ce qui est coûteux. De plus la matrice $W$ initialement creuse devient pleine par régularisation, $\tilde{W} = W + \omega NN^T$. Pour une exploitation efficace dans un logiciel ayant vocation à traiter des centaines de milliers de contacts la mise en œuvre de la méthode doit être revue. En premier lieu l'algorithme de Gauss-Seidel non linéaire doit être modifié en effectuant sur la partie pleine de $\tilde{W}$ une itération de type Jacobi plutôt que Gauss-Seidel. En effet les matrices $N$ et $Z = NN^T$ ne doivent pas être stockées ; seul le calcul de produits matrice-vecteur peut être effectué. La méthode de Gauss-Seidel non linéaire est ainsi modifiée et détaillée dans Algorithme 6.

**Algorithm 6 :** Solveur NLGS-Jac sans stocker la matrice $Z$

Mettre en oeuvre une procédure de calcul de $Nx$ et $N^T x$

Calculer $N^T e_\alpha$, $\alpha = 1, \dots, n_{contact}$

Calculer la partie (bloc-)diagonale de $Z = NN^T : Z_{\alpha\alpha} = e_\alpha^T NN^T e_\alpha$, $\alpha = 1, \dots, n_{contact}$

Calculer les termes supplémentaires au membre droit : $\tilde{z} = Z\tilde{r}$ ($y = N^T\tilde{r}$, puis $\tilde{z} = Ny$)

for $k = 1, \dots, convergence$ do

Calculer $z^{k-1} = Zr^{k-1}$ ($y = N^T r^{k-1}$, puis $z^{k-1} = Ny$)

for $\alpha = 1, \dots, n_{contact}$ do {Boucle sur les contacts}

Déterminer $(r_\alpha^k, v_\alpha^k)$ tel que

$$ \left\{ \begin{array}{l} \tilde{W}_{\alpha\alpha} r_{\alpha}^{k} - v_{\alpha}^{k} = -v_{\alpha}^{d} + \tilde{z}_{\alpha} - \sum_{\beta<\alpha} W_{\alpha\beta} r_{\beta}^{k} - \sum_{\alpha<\beta} W_{\alpha\beta} r_{\beta}^{k-1} - z_{\alpha}^{k-1} + Z_{\alpha\alpha} r_{\alpha}^{k-1} \\ \mathcal{R}(r_{\alpha}, v_{\alpha}) = 0 \end{array} \right. $$

end for
end for

La première phase de l'étude concerne le choix du facteur $\omega$. Une première batterie de tests semble montrer qu'un $\omega$ de l'ordre d'un vingtième de la norme euclidienne de $W$ permet d'obtenir une vitesse de convergence de NLGS-Jac comparable à celle de NLGS générique. Une étude plus approfondie est à mener sur ce sujet.

En second lieu une représentation creuse du noyau doit être effectuée en utilisant les travaux pointus de Le Borne [6]. Dans [6] la motivation est de déterminer une base du noyau de l'opérateur divergence discret en contexte éléments finis. Nous proposons ici une extension de ces travaux à l'opérateur de connectivité $H$ entre variables globales $(V, R)$ et locales $(v, r)$. En effet on peut démontrer que le noyau de $W$ est celui de l'opérateur $H_f$, noté $H$ par souci de simplicité dans la suite, défini dans le tableau des notations Tableau 1. Sans entrer dans les détails la méthode de Le Borne consiste à représenter (et stocker) une base du noyau $N$ par 3 matrices creuses, $L$ triangulaire inférieure, $U$ triangulaire supérieure et surtout $R_1$ factorisée de Choleski de la matrice $HH^T$. Une procédure permet ensuite de calculer le produit $Nx$ ou $N^T x$ pour tout vecteur $x$. Des outils de renumérotation optimale sont également requis et
<figure><img src="image_6.png" /><figcaption>Fig. 5 – Echantillon à 500 grains - grains flottants : un seul contact (gris), aucun contact (noir)</figcaption></figure>

disponibles sous MATLAB. Nous avons donc testé cette approche sur des échantillons de 30 à 500 grains, préalablement élaborés sur la plateforme LMGC90 (http://www.lmgc.univ-montp2.fr/LMGC90/).

La méthode suppose que $H^T$ soit de rang maximal afin que la factorisation de Choleski existe. Ce n'est pas toujours le cas dans les assemblages granulaires, mais l'on peut adapter la méthode à ces situations. En effet prenons l'exemple à 500 billes représenté sur la Figure 5. La matrice $H$ a 1500 lignes correspondant aux degrés de libertés des 500 billes ($3 \times 500$) et 2044 colonnes associées aux deux composantes des forces de 1022 contacts. Le rang de $H^T$ n'est que de 1487 ; les 13 degrés de liberté d'écart proviennent des billes "flottantes" : 7 billes n'ayant qu'un seul contact et 1 degré des liberté, et 2 billes sans aucun contact, soit 3 degrés de liberté, soit $7 \times 1 + 2 \times 3 = 13$ degrés de liberté. Les billes et les contacts isolés ainsi identifiés peuvent être éliminés de $H$, fournissant une matrice $H$ réduite de rang maximal et disposant du même noyau. La méthode de Le Borne peut alors être mise en oeuvre.

<table><caption>Tableau 2 – Nombre de coefficients non nulls dans la représentation creuse d'une base du noyau de $H$ pour l'échantillon à 500 billes. (a) $H$ a été préalablement réduite.</caption><thead><tr><th>Matrices</th><th>H</th><th>W</th><th>W&#771;</th><th>R<sub>1</sub></th><th>L</th><th>U</th><th>Total (R<sub>1</sub>, L, U)</th></tr></thead><tbody><tr><td>Non zeros</td><td>9 758 (a)</td><td>27 858</td><td>4 177 936</td><td>32 583</td><td>12 487</td><td>43 254</td><td>88 324</td></tr></tbody></table>

La Tableau 2 fournit l'encombrement mémoire engendré par la représentation creuse du noyau de $H$ pour l'échantillon à 500 billes ; près de 90 000 coefficients non nuls, soit 3 fois l'encombrement de $W$ lui-même à comparer également aux 4 millions de termes qu'exigerait un stockage plein de $W$.

Le Tableau 3 détaille le nombre d'itérations nécessaire à la convergence de trois méthodes pour ré-soudre le même problème non régulier : l'algorithme NLGS générique (Algorithme 2) mais appliqué à la matrice $W$ et non $W$ (et par conséquent utilisant un stockage plein), l'algorithme NLGS-Jac (Algorithme 6) utilisant un stockage plein de $N = KerH$, enfin le même algorithme avec une représentation creuse de $N$. Le surcoût en itérations de NLGS-Jac reste limité. Le gain en temps CPU est sensible, les calculs
<table><caption>Tableau 3 – Nombre d'itérations et temps CPU - échantillon à 500 billes.</caption><thead><tr><th>Method or Phase</th><th>NLGS (full)</th><th>NLGS-Jac (full)</th><th>Compute Z (sparse)</th><th>Compute DZ (sparse)</th><th>NLGS Z (sparse)</th><th>Total (sparse)</th></tr></thead><tbody><tr><td>iterations</td><td>20</td><td>23</td><td></td><td></td><td>23</td><td></td></tr><tr><td>CPU time (s)</td><td>25.7</td><td>33.2</td><td>0.01</td><td>1.</td><td>7.</td><td>8.01</td></tr></tbody></table>

préparatoires, calcul de la représentation creuse de Z et de la partie diagonale de de Z, représentant un coût marginal. On ne peut comparer aisément ce coût avec celui de la résolution standard (Algorithme 2 appliqué à W), car, d'un pas à un autre le nombre d'itérations à convergence peut être erratique (de quelques dizaines à plusieurs milliers) sans que la solution soit comparable du fait de la non unicité de la solution. Bien entendu ces premières tendances demandent à être confirmé sur des échantillon de grande taille dans une plateforme logicielle adaptée.

# 6 Conclusions

L'étude proposée ici est essentiellement exploratoire, même si elle a exigé des développements numériques importants sur des logiciels standards (MAPLE et MATLAB) et une multitude de petites simulations. La stratégie mise en avant répond à quelques questions posées par l'approche NSCD, en essayant d'isoler prioritairement le noyau de l'opérateur de Delassus dans les causes des convergences erratiques de l'algorithme NLGS. Cependant d'autres questions apparaissent qui nécessiteront des études approfondies. La tâche la plus ardue dans un premier temps reste l'implémentation des outils algébriques de la représentation creuse du noyau de W dans une plateforme dédiée aux éléments discrets.

# Références

[1] Moreau J.J., *Numerical aspects of sweeping process*, Computer Methods in Applied Mechanics and Engineering, Vol. 177, Issues 3-4, 329-349, 1999.

[2] Alart, P., Iceta, D. and Dureisseix, *A nonlinear domain decomposition formulation with application to granular dynamics*, Computer Methods in Applied Mechanics and Engineering, Vol 205-208, 59-67, 2012.

[3] Visseq, V., Martin, A., Iceta, D., Azema, E., Dureisseix, D. and Alart, P., *Dense granular dynamics analysis by a domain decomposition approach*, Computational Mechanics, 709-723, 2012.

[4] McNamara S. and Young W.R. *Inelastic collapse and clumping on a one dimensional granular medium*, Phys. Fluids A 4, 496-504, 1991

[5] Jean M., *Numerical simulation of granular materials*, ISTE, 2010, pp. 149-315.

[6] Le Borne S., *Block computation and representation of a sparse nullspace basis of a rectangular matrix*, Linear Algebra and its Applications, Volume 428, Issues 11-12, 2455-2467, 2008