définition géométrique des surfaces en synthèse d'images. Les nœuds sont au nombre de 4 par élément, coïncidant avec les sommets des tétraèdres. Ainsi, puisque des tétraèdres voisins partagent les mêmes nœuds, la méthode garantit la continuité de la solution entre deux éléments. Enfin, la fonction d'interpolation se réduit à sa plus simple expression, c'est-à-dire une fonction linéaire. Ce modèle particulier d'élément est appelé élément de type P1.

Si l'on considère l'équation (3.7) qui régit le mouvement d'un milieu continu déformable en utilisant comme inconnue le déplacement $u(t) = r(t) - r^0$ subi par la particule plutôt que sa position instantanée $r(t)$, la méthode des éléments finis de type P1 donne le système différentiel suivant :

$$M\ddot{u} + D\dot{u} + Ku = f \quad (3.14)$$

Le système différentiel (3.14) est très similaire au système (3.11) obtenu par l'utilisation des différences finies dans l'approche de Terzopoulos. Pour des raisons de performance, on suppose dans la plupart des cas que la matrice de raideur est constante. Cette hypothèse, qui signifie que les efforts élastiques sont proportionnels aux déplacements, me place dans le cadre de la théorie de l'élasticité linéaire. L'équation (3.14) n'est alors valide que pour des petits déplacements restreignant ainsi le champ d'application de la méthode des éléments finis.

Parmi les premiers travaux utilisant la méthode des éléments finis, on peut citer celui de Gourret et Thalmann [Gourret 1989] dans lequel ils exploitent cette méthode pour simuler les déformations de la peau d'une main serrant une balle. Même en négligeant la dynamique pour ne résoudre que des problèmes statiques successifs, les temps de calculs demeurent très importants. Collier et al [Collier 1991] ainsi que Eischen et al [Eischen 1996] s'inspirent des méthodes de simulation des membranes pour mettre au point un modeleur de textiles tissés pour en étudier le drapé. Dans ce cas de figure aussi, la méthode des éléments finis se révèle lourde et donc peu adaptée à des applications temps réel.

Pour pallier aux problèmes de temps de calcul, Bro-Nielsen et Cotin [Bro-Nielsen 1996] proposent une méthode dite de condensation qui consiste à traiter les nœuds de la surface et ceux de l'intérieur séparément. En considérant que seul le déplacement des nœuds de la surface est visuellement pertinent, ils diminuent drastiquement le nombre de nœuds à animer. Pour les besoins de sa thèse, Cotin [Cotin 1997] propose une autre démarche d'optimisation qui repose sur la propriété de superposition des efforts en élasticité linéaire. Partant de l'idée que pour une déformation double l'intensité des efforts sera double, il pré-calcule les efforts d'élasticité interne pour un élément de référence unitaire. Les efforts à appliquer à un élément de la configuration instantanée sont obtenus par homothétie, autorisant des applications temps réel en simulation chirurgicale. Mais tous les modèles précédents reposent sur l'invariance de la matrice des raideurs $K$, restreignant le champ d'application de ces méthodes aux petits déplacements. Plus récemment, Picinbono [Picinbono 2001] propose une utilisation de la méthode des éléments finis dans un contexte temps réel en utilisant une énergie élastique non linéaire en s'appuyant sur le modèle de St Venant-Kirchhoff. La non-linéarité de la déformation vis-à-vis du déplacement engendre une baisse significative de la performance, incompatible avec leur optique temps-réel. Picinbono utilise en fait un modèle adaptatif lui permettant de basculer entre une représentation linéaire ou non-linéaire de la déformation en fonction de
l'ordre de grandeur du déplacement. Dans la mesure où le contexte applicatif consiste en des manipulations ponctuelles d'un organe virtuel inerte, le modèle adaptatif se révèle tout à fait efficace.

### 2.2.3 Modèle « dynamique modale »

La démarche de Pentland et Williams [Pentland 1989] qualifiée de dynamique modale est en fait une stratégie de simplification des équations (3.14) issues de la méthode des éléments finis. Ils supposent que la matrice d'amortissement est une combinaison linéaire de la matrice des masses et de la matrice des raideurs conformément au modèle de Raleigh (i.e. $\mathbf{D} = \alpha \mathbf{M} + \beta \mathbf{K}$). Ils déterminent alors une base de vecteurs propres $(\mathbf{\Phi}^i)_{1\le i \le n}$ M-orthonormale de la matrice $\mathbf{M}^{-1}\mathbf{K}$.
On a alors les relations suivantes :

$$ 
\begin{gathered}
\mathbf{K}\mathbf{\Phi}^i = \omega_i^2 \mathbf{M}\mathbf{\Phi}^i \quad \text{avec } \omega_i^2 \text{ valeur propre de } \mathbf{M}^{-1}\mathbf{K} \\
\mathbf{\Phi}^T \mathbf{M}\mathbf{\Phi} = \mathbf{I} \quad \text{où } \mathbf{\Phi} = (\mathbf{\Phi}^1 \quad \dots \quad \mathbf{\Phi}^n)
\end{gathered}
\qquad (3.15) $$

De l'équation (3.15), on déduit l'égalité suivante :

$$ \mathbf{K}\mathbf{\Phi} = \mathbf{M}\mathbf{\Phi}\mathbf{\Omega}^2 \quad \text{avec } \mathbf{\Omega}^2 = \begin{pmatrix} \omega_1^2 & & & \\ & \ddots & & \\ & & \omega_n^2 & \\ & & & \end{pmatrix} \Rightarrow \mathbf{\Phi}^T \mathbf{K}\mathbf{\Phi} = \mathbf{I} \qquad (3.16) $$

Sous ces conditions, le système différentiel (3.14) admet une nouvelle expression.

$$ \ddot{\mathbf{u}} + \bar{\mathbf{D}}\dot{\mathbf{u}} + \mathbf{\Omega}^2 \bar{\mathbf{u}} = \mathbf{\Phi}^T \mathbf{f} \quad \text{avec} \quad \mathbf{u} = \mathbf{\Phi} \bar{\mathbf{u}} \qquad (3.17) $$

où la matrice $\bar{\mathbf{D}}$ est diagonale. Ainsi si l'on pose $\bar{\mathbf{f}} = \mathbf{\Phi}^T \mathbf{f}$, le système diagonal (3.17) est équivalent aux $n$ équations linéaires indépendantes suivantes :

$$ \ddot{\bar{u}}_i + \bar{D}_{ii} \dot{\bar{u}}_i + \omega_i^2 \bar{u}_i = \bar{f}_i \quad 1 \le i \le n \qquad (3.18) $$

En changeant de base de représentation, le problème initial formulé en termes de déplacement s'interprète grâce à l'équation (3.18) comme une superposition de modes vibratoires indépendants. Un premier avantage de ce type d'approche réside dans la nature fortement parallèle du nouveau problème. Les équations étant indépendantes ; elles peuvent être résolues simultanément au sein d'une architecture parallèle. L'autre avantage réside dans la possibilité de supprimer certains modes vibratoires sans altérer le comportement global du corps déformable. Puisque la fréquence du mode vibratoire est proportionnelle aux valeurs propres, les équations associées à de grandes valeurs propres peuvent être alors négligées.

### 2.2.4 Modèle « dynamique globale »

Cette approche proposée par Witkin et Welch [Witkin 1990] s'appuie sur les équations de Lagrange pour décrire le mouvement d'un corps déformable.
$$ \frac{d}{dt} \left( \frac{\partial K}{\partial \dot{x}_i} \right) - \frac{\partial K}{\partial \dot{x}_i} = Q_i \qquad (3.19) $$

Ils choisissent délibérément de ne considérer que des déformations équivalentes à des transformations linéaires. Ainsi, la position instantanée d'une particule de l'objet déformable est déduite de sa position **p** dans la configuration de référence et de la matrice de transformation **R**

$$ x_i = \sum_{j=1}^{3} R_{ij} p_j + R_{i4} \qquad (3.20) $$

On considère alors les coefficients de la matrice **R** comme étant les nouvelles inconnues de Lagrange du système. Ce changement d'inconnues permet de déboucher sur l'expression d'un nouveau système différentiel constitué de 16 équations seulement.

$$ \ddot{\mathbf{R}}\mathbf{M} = \mathbf{\tilde{Q}} \quad \text{avec} \quad \tilde{Q}_{ij} = \sum_{k=1}^{n} f_k \frac{\partial x_k}{\partial R_{ij}} \qquad (3.21) $$

En plus de diminuer considérablement le nombre des inconnues, cette méthode en animant le corps dans sa globalité fait disparaître tous les phénomènes oscillatoires permettant ainsi d'utiliser des pas de temps relativement élevés lors de l'intégration numérique. Malheureusement, cette méthode se restreint aux transformations linéaires. Deux ans plus tard, Baraff et Witkin [Baraff 1992] proposent une extension de la démarche aux transformations polynomiales mais se limitent au cas des polynômes du second degré.

## 2.3 Synthèse

<table><caption>Tableau 1 Approches classiques pour l'animation dynamique de corps déformables</caption><thead><tr><th rowspan="2">Modèle</th><th rowspan="2">Structuré</th><th>Mécanique</th><th>Préservation</th><th rowspan="2">Multi Echelle</th></tr><tr><th>des Milieux<br>Continus</th><th>Continuité<br>Matérielle</th></tr></thead><tbody><tr><td>Masses-Ressorts</td><td>OUI</td><td>NON</td><td>NON</td><td>OUI</td></tr><tr><td>Particules</td><td>NON</td><td>NON</td><td>NON</td><td>OUI</td></tr><tr><td>SPH</td><td>NON</td><td>OUI</td><td>OUI</td><td>OUI</td></tr><tr><td>Modèle Terzopoulos</td><td>OUI</td><td>OUI</td><td>NON</td><td>OUI</td></tr><tr><td>Eléments finis</td><td>OUI</td><td>OUI</td><td>NON</td><td>NON</td></tr><tr><td>Dynamique modale</td><td>OUI</td><td>OUI</td><td>NON</td><td>OUI</td></tr><tr><td>Dynamique globale</td><td>OUI</td><td>OUI</td><td>OUI</td><td>NON</td></tr><tr><td>Variétés Splines Dynamiques</td><td>OUI</td><td>OUI</td><td>OUI</td><td>OUI</td></tr></tbody></table>

En guise de synthèse, je rassemble les approches précédemment citées au sein du tableau 1. Le modèle qualifié de « variétés splines dynamiques » correspond au modèle alternatif et innovant
que je propose. Les propriétés de ce modèle, mentionnées dans les colonnes du tableau, sont développées dans la section suivante.

## 3 Modèle « Variétés Splines Dynamiques »

Tout comme les modèles présentés précédemment, le modèle que je propose se plie à la contrainte que constitue la réduction du nombre de paramètres caractéristiques. J'ai ainsi opté pour une modélisation qualifiée de **modélisation spline** des objets déformables continus. Pour préciser cette notion de modélisation spline, j'emprunterai cette remarque formulée par Carole Blanc dans sa thèse [Blanc 1994] :

> « Les objets splines sont des modèles mathématiques qui permettent d'associer une représentation paramétrique continue à un ensemble fini de points de l'espace »

La configuration d'un corps déformable est entièrement déterminée à partir d'un réseau fini de points caractéristiques de l'espace $\mathbb{R}^3$, appelés **points de contrôle**. La position de toute particule constitutive de l'objet déformable est alors obtenue par combinaison linéaire des positions des points de contrôle, pondérée par des **fonctions de mélange**. Ce choix de modélisation a de nombreux avantages dont :

* La garantie de la finitude des degrés de liberté correspondant aux coordonnées des points de contrôle. En effet, les positions de toutes les particules de l'objet spline sont entièrement conditionnées par les positions des points de contrôle.

* La souplesse et la richesse de formes qu'offre le modèle spline. En effet, l'utilisateur responsable de la modélisation dispose d'un outil mathématique puissant lui permettant de réaliser des objets de synthèse aux formes très variées.

* La popularité du modèle spline n'est plus à prouver. Il s'est largement imposé avec l'épanouissement des techniques de conception et de fabrication assistées par ordinateur. Le modèle spline est donc un élément incontournable dans la culture de l'informatique graphique au vu du nombre important de logiciels de modélisation qui l'utilisent. Cet état de fait me donne la garantie de pouvoir concevoir aisément les objets que je souhaite par la suite animer.

De par l'identification paramétrique des particules, ce modèle peut être comparé à celui proposé par Terzopoulos. En effet, le modèle que je propose englobe à la fois les définitions de courbes, surfaces et volumes splines. Mais contrairement au modèle de Terzopoulos, j'attache une grande importance à la nature continue de l'objet déformable. Je considère une répartition continue de la matière en accord avec la continuité du modèle géométrique. Dans ce modèle, la masse n'est donc pas répartie parmi un ensemble discret de points matériels caractéristiques, le distinguant des modèles masses-ressorts ou à base de particules.

La particularité de la définition géométrique d'une variété spline (combinaison linéaire de points de contrôle pondérée par les valeurs des fonctions de mélange) me permet de procéder à de nombreuses simplifications formelles des équations du mouvement. Ainsi, le système différentiel issu des simplifications, demeure équivalent aux équations de départ tout en exhibant des propriétés structurelles intéressantes (matrice des masses composée de trois blocs identiques constants par exemple). Il est important de mentionner que je n'ai recours à aucune
étape de discrétisation supplémentaire. Cette démarche de conservation des équations continues du mouvement est similaire à l'approche utilisée pour le modèle SPH : des équations d'état continues sont plus proches de la nature du milieu qu'elles simulent et permettent donc de calculer une solution plus proche de la solution exacte.

Enfin la généricité de modèle spline permet d'envisager des formes alternatives pour la génération du mouvement. Au lieu d'animer un objet complexe par l'intermédiaire de ses paramètres caractéristiques souvent trop nombreux, on peut considérer son évolution de manière globale en le supposant plongé dans un volume paramétrique défini à partir d'un réseau régulier de points de contrôle. Cette démarche est assimilable à la méthode des « Free Form Deformations » [Sederberg 1986] [Coquillart 1990] [Coquillart 1991] où chaque point de l'objet à déformer est considéré à paramètres constants au sein d'un volume englobant. Contrairement aux approches déjà proposées, l'innovation de mon approche repose sur le fait que je considère les points de contrôle du volume englobant comme les nouvelles inconnues des équations du mouvement. Cette démarche peut être perçue comme une généralisation du modèle de « dynamique globale » proposé par Witkin et Welch dans la mesure où il ne se restreint pas à des transformations linéaires ou quadratiques. Cette démarche est détaillée dans le chapitre consacré à la réduction de paramétrage.
<figure><img src="image_34.png" /></figure>
<figure><img src="image_35.png" /></figure>

# VARIETES SPLINES DE DIMENSION D

Pour les besoins de sa thèse, Jean-Michel Nourrit a proposé un modèle de fils définis à partir de courbes splines matérielles. Ce choix, couplé avec une description du mouvement à l'aide des équations de Lagrange, révèle de nombreux avantages dont la garantie de la continuité de la définition géométrique du système mais aussi de sa répartition massique. Ces premiers résultats satisfaisants nous ont encouragé à poursuivre dans cette voie en adoptant une démarche dictée par la volonté du groupe d'animation d'aboutir à la réalisation d'une plate-forme d'animation générique.

Au sein de cette partie, je commence par proposer une expression générale d'une courbe spline, indépendante de la nature de ses fonctions de mélange, généralisant ainsi le travail effectué par Jean-Michel Nourrit.

Dans un deuxième temps, au vu de la similarité des formulations splines associées à des courbes, des surfaces et des volumes, je propose une deuxième abstraction dans l'étape de modélisation en donnant une nouvelle formulation indépendante à la fois des fonctions de mélange et désormais de la dimension intrinsèque de l'objet. Ainsi, les courbes, surfaces et volumes splines se

trouvent être rassemblés dans un cadre théorique unifié que je choisis de qualifier de variété spline de dimension d.

L'implantation informatique du concept de variétés splines se veut tout aussi générique. Je détaille comment, en exploitant les notions de polymorphisme et de généricité de la conception objet, j'aboutis à une structure informatique des variétés splines fidèle à leur définition mathématique.

Le point fort de cette modélisation réside dans notre capacité à conserver une description continue du milieu sans avoir recours à une étape de discrétisation supplémentaire. Une des lacunes principales du modèle initial proposé par Jean-Michel résidait dans l'expression de efforts élastiques internes qui auparavant étaient modélisés par des ressorts disposés dans un réseau discret de points caractéristiques de la variété. Je propose alors une première formalisation des efforts issues d'une énergie continue de déformation se limitant pour le moment au cas des courbes.

Enfin, j'exploite la richesse du modèle de variétés splines pour proposer une méthode de réduction de paramétrage que j'applique à la simulation de textiles tricotés.
<figure><img src="image_36.png" /></figure>
# GEOMETRIE DES VARIETES SPLINES

## 1 Introduction

Ce chapitre a pour objectif de détailler la démarche que j’ai entreprise pour bâtir le modèle géométrique de variété spline. Je commence par donner une expression générique d’une courbe spline, indépendante de la nature de ses fonctions de mélange. Cette expression garantit dans tous les cas de figures la propriété de localité par l’introduction du concept de segments splines. Puis, au vu des similarités existantes dans les expressions unidimensionnelle, bidimensionnelle et tridimensionnelle des objets splines, je propose un modèle novateur, considéré comme l’unification des représentations linéiques, surfaciques et volumiques. Je termine en présentant le modèle informatique codant une variété spline qui lui aussi reflète cette volonté d’unification de la représentation. Le modèle de variété spline dans sa dimension géométrique, cinématique et dynamique a fait l’objet de publications dont une dans la revue internationale The Visual Computer [Remion 2001].

## 2 Courbe spline

Afin de préciser ce que j’entend par modélisation spline, j’emprunterai à Carole Blanc la définition d’une courbe spline qu’elle proposa dans sa thèse de doctorat [Blanc 1994] :

Soient un ensemble de points $\mathbf{p}^i \in \mathbb{R}^3$ et des fonctions $b_i : [0,1] \to \mathbb{R}$ avec $1 \le i \le n$. On appelle **courbe spline** l’ensemble continu C des points dont la position est donnée par la formule suivante :

$$ \mathbf{p}(\omega) = \sum_{i=1}^{n} b_i(\omega) \mathbf{p}^i \quad \omega \in [0,1] \qquad (1.1) $$

Les points $\mathbf{p}^i \in \mathbb{R}^3$ sont appelés les **points de contrôle** et l’ensemble de ces points, noté $\mathbb{R}$, est appelé **réseau de contrôle** de la courbe spline. Les fonctions $b_i$ correspondent aux **fonctions de mélange** (i.e. **d’influence**) associées respectivement aux points $\mathbf{p}^i \in \mathbb{R}^3$ (Figure 25).
<figure><img src="image_37.png" /><figcaption>Figure 25 Courbe spline définie à partir de n points de contrôle.</figcaption></figure>

Il me semble important de signaler que la définition d'une courbe spline (1.1) est entièrement indépendante de la nature de ses fonctions de mélange. Cette définition fournit une expression mathématique unique, englobant toutes les familles de courbes splines. Ainsi les courbes interpolantes telles que les courbes d'Hermite [Barnhill 1974], de Catmull-Rom [Catmull 1974], TCB [Kochanek 1984] ou bien les courbes approximantes telles que les courbes de Bézier [Bézier 1966, 1977] [Farin 1990], les B-splines [Riesenfeld 1973] [Foley 1989], les β-splines [Barsky 1981, 1983], les NURBS [Piegl 1991] sont réunies sous une même formulation et autorisent une certaine abstraction dans le processus de modélisation spline vis-à-vis de la nature des fonctions de mélange.

Cependant, la formulation proposée en (1.1) met en évidence une forte dépendance entre les différents points de contrôle $\mathbf{p}^i$, suggérant que les $n$ points de contrôle interviennent dans le calcul de la position $\mathbf{p}(\omega)$. Je choisis alors d'adopter une définition alternative pour une courbe spline. En fait, à partir d'un réseau de contrôle commun $\mathbb{R}$, je construis plusieurs courbes splines contiguës que je nommerai par la suite **segments splines**. Le terme de courbe spline désigne la réunion de ces segments contigus. Cette idée (illustrée par la Figure 26) peut être mathématiquement formalisée de la manière suivante :

Soient un ensemble de points $\mathbf{p}^i \in \mathbb{R}^3$ avec $1 \le i \le n$, un nombre $ns$ de segments et des fonctions $b_i^j : [0,1] \to \mathbb{R}$ avec $1 \le i \le n$, $1 \le j \le ns$. On appelle **segments splines** les $ns$ ensembles continus $C^j$ des points dont la position est donnée par l'expression suivante :

$$ \mathbf{p}^j(\omega) = \sum_{i=1}^{n} b_i^j(\omega) \mathbf{p}^i \quad 1 \le j \le ns, \omega \in [0,1] \quad (1.2) $$

L'expression de la courbe spline C est déduite de la relation (1.3) :

$$ C = \bigcup_{j=1}^{ns} C^j \quad (1.3) $$
<figure><img src="image_38.png" /><figcaption>Figure 26 Courbe spline composée de ns segments.</figcaption></figure>

Ainsi, un point de la courbe spline C est identifié par son numéro de segment $j$ et son indice paramétrique $\omega$ dans ce segment. On définit alors une **fonction de position globale** pour l'ensemble des points de la courbe spline.

$$ \mathbf{p}: \begin{cases} \{1, \dots, ns\} \times [0,1] \to \mathbb{R}^3 \\ (j, \omega) \mapsto \mathbf{p}(j, \omega) = \mathbf{p}^j(\omega) \end{cases} \quad (1.4) $$

L'introduction de cette notion de segments splines a plusieurs conséquences importantes. En l'occurrence la définition (1.2) permet d'imposer la propriété de localité à une courbe spline et ce, quelque soit la nature des fonctions de mélange. Cette propriété de localité garantit que l'évaluation de la position $\mathbf{p}(j, \omega)$ ne dépend que d'un petit nombre de points de contrôle. Géométriquement, le déplacement d'un point de contrôle n'affecte qu'une partie de la courbe. De façon pratique, cette propriété est réalisée en imposant que, pour un segment $j$ dépendant d'un sous-ensemble $\mathbb{R}^j$ du réseau de contrôle $\mathbb{R}$, les fonctions de mélange $b_i^j$ soient identiquement nulles en dehors du sous ensemble $\mathbb{R}^j$. La Figure 27 illustre cette propriété de localité en montrant son effet sur l'allure générale d'une courbe spline.

<figure><img src="image_39.png" /><figcaption>Figure 27 Illustration de la localité.</figcaption></figure>

Afin d'illustrer le modèle de « segments splines » que je propose, je donne quelques exemples de courbes splines classiques exprimées dans mon formalisme. Mais auparavant, je rappelle plusieurs définitions relatives aux propriétés des courbes splines.
Une courbe spline vérifie la **propriété de normalité** si ses fonctions de mélange $b_i^j$ vérifient l'égalité suivante :

$$ \forall j, 1 \le j \le ns: \sum_{i=1}^{n} b_{i}^{j}(\omega) = 1 \quad \forall \omega \in [0,1] \qquad (1.5) $$

Les courbes splines vérifiant l'égalité (1.5) sont alors qualifiées de **courbes normales**.

Une courbe spline est une **courbe d'interpolation** des points $\mathbf{p}^k$ ($1 < k < n$) si elle passe par l'ensemble des points de contrôle $\mathbf{p}^k$ ($1 < k < n$).

Les courbes splines qui ne sont pas des courbes d'interpolation sont des **courbes d'approximation**.

## 2.1 Courbe d'Hermite

<figure><img src="image_40.png" /><figcaption>Figure 28 Courbe d'Hermite composée de 3 segments, 6 points de contrôle.</figcaption></figure>

Les fonctions de mélange s'expriment de la façon suivante :

$$ b_{i}^{j}(\omega) = \begin{cases} H_{i-j}(\omega) & \text{si } j \le i \le j+3 \\ 0 & \text{sinon} \end{cases} \qquad (1.6) $$

où les fonctions $H_i$ sont les polynômes d'Hermite définis comme :

$$ \begin{aligned} H_1(\omega) &= -\omega + 2\omega^2 - \omega^3 \\ H_2(\omega) &= 1 - 2\omega^2 + \omega^3 \\ H_3(\omega) &= \omega + \omega^2 - \omega^3 \\ H_4(\omega) &= -\omega^2 + \omega^3 \end{aligned} \qquad (1.7) $$

Les courbes d'Hermite sont des courbes normales, d'interpolation des points $\mathbf{p}^k$ ($1 < k < n$). Et les courbes d'Hermite sont des courbes de continuité $C^1$. En effet, deux segments splines d'Hermite $j$ et $j+1$ coïncident au point $\mathbf{p}(j,1) = \mathbf{p}(j+1,0) = \mathbf{p}^{j+2}$. De plus, la pente à droite et à gauche au point $\mathbf{p}^{j+2}$ est égale à la moitié du vecteur d'extrémités $\mathbf{p}^{j+1}$ et $\mathbf{p}^{j+3}$ (Figure 28).
## 2.2 Courbe de Catmull-Rom

<figure><img src="image_41.png" /><figcaption>Figure 29 Courbes de Catmull-Rom composées de 3 segments, 6 points de contrôle pour des tensions égales à 0, 1 et 2.</figcaption></figure>

Les fonctions de mélange d'une courbe de Catmull-Rom s'expriment de la façon suivante :

$$b_i^j(\omega) = \begin{cases} C_{i-j}(\omega) & \text{si } j \le i \le j+3 \\ 0 & \text{sinon} \end{cases} \qquad (1.8)$$

où les fonctions $C_i$ sont les polynômes définis comme:

$$\begin{align*}
C_1(\omega) &= \rho(-\omega + 2\omega^2 - \omega^3) \\
C_2(\omega) &= 1 + (\rho - 3)\omega^2 + (2 - \rho)\omega^3 \\
C_3(\omega) &= \rho\omega + (3 - 2\rho)\omega^2 + (\rho - 2)\omega^3 \\
C_4(\omega) &= \rho(-\omega^2 + \omega^3)
\end{align*} \qquad (1.9)$$

Le paramètre $\rho$ est un paramètre de tension qui influe sur la norme de la pente autour de chaque point de contrôle. Les courbes de Catmull-Rom sont des courbes normales, d'interpolation des points $\mathbf{p}^k$ ($1 < k < n$). Et les courbes de Catmull-Rom sont aussi des courbes de continuité $C^1$, pour les mêmes raisons que les courbes d'Hermite. D'ailleurs, les courbes d'Hermite sont un cas particulier de courbes Catmull-Rom pour une tension unitaire (Figure 29).

## 2.3 Courbe de Kochanek (i.e. courbe TCB)

Les fonctions de mélange s'expriment de la façon suivante :

$$b_i^j(\omega) = \begin{cases} K_{i-j}(\omega) & \text{si } j \le i \le j+3 \\ 0 & \text{sinon} \end{cases} \qquad (1.10)$$

où les fonctions $K_i$ sont les polynômes définis comme:
$$
\begin{align*}
K_1(\omega) &= -k_1\omega^3 + 2k_1\omega^2 - k_1\omega \\
K_2(\omega) &= (2 + k_1 - k_3 - k_4)\omega^3 + (k_3 + 2(k_4 - k_1) - 3)\omega^2 + (k_4 - k_1)\omega + 1 \tag{1.11} \\
K_3(\omega) &= (k_4 + k_3 - k_2 - 2)\omega^3 + (k_2 - k_3 - 2k_4 + 3)\omega^2 + k_4\omega \\
K_4(\omega) &= k_2\omega^3 - k_2\omega^2
\end{align*}
$$

$$
\begin{align*}
\text{avec } k_1 &= \frac{1}{2} (1 - \text{tension}) (1 + \text{continuité}) (1 + \text{bias}) \\
k_2 &= \frac{1}{2} (1 - \text{tension}) (1 + \text{continuité}) (1 - \text{bias}) \\
k_3 &= \frac{1}{2} (1 - \text{tension}) (1 - \text{continuité}) (1 + \text{bias}) \\
k_4 &= \frac{1}{2} (1 - \text{tension}) (1 - \text{continuité}) (1 - \text{bias})
\end{align*}
$$

Les courbes de Kochanek sont des courbes normales, d'interpolation des points **p**<sup>*k*</sup> (1&lt;*k* &lt; *n*)
de continuité C<sup>1</sup>.

2.4 Courbe de Bézier

Les fonctions de mélange d'une courbe de Bézier de degré 5 s'expriment de la façon suivante :

$$
b_i^j(\omega) = \begin{cases} B_{i-4(j-1)+1}^5(\omega) & \text{si } 4(j-1)+1 \le i \le 4(j-1)+5 \\ 0 & \text{sinon} \end{cases} \quad (1.12)
$$

où les fonctions $B_i^5$ sont les polynômes de Bernstein définis comme:

$$
B_i^5(\omega) = C_5^i (1-\omega)^{5-i} \omega^i \quad \text{où} \quad C_5^i = \frac{5!}{(5-i)!i!} \qquad (1.13)
$$

Les courbes de Bézier de degré 5 sont des courbes normales, d'approximation mais
interpolantes aux points **p**<sup>5k</sup>. Les courbes de Bézier sont de continuité C<sup>0</sup> seulement. On peut
garantir une continuité C<sup>1</sup> autour d'un point **p**<sup>*k*</sup>, extrémité commune de deux segments splines
de Bezier, s'il est milieu du segment [**p**<sup>*k*−1</sup>, **p**<sup>*k*+1</sup>].

$$
\mathbf{p}^k = \frac{\mathbf{p}^{k-1} + \mathbf{p}^{k+1}}{2} \qquad (1.14)
$$

3 Extension aux dimensions supérieures

3.1 Surface et volume spline

Le modèle spline est un formalisme mathématique intrinsèquement riche et puissant. Hormis sa
capacité d'abstraction vis-à-vis des fonctions de mélange, la définition d'une courbe spline (1.2)
est naturellement transposable aux cas des surfaces et des volumes splines.
Ainsi, si l'on considère un réseau bidimensionnel de points de contrôle $\mathbf{p}^{i_1, i_2}$ associés à des fonctions de mélange de deux variables réelles $\omega_1$ et $\omega_2$, on définit des sous-surfaces splines $S^{i_1, i_2}$ communément appelées « patches » dont la position des points est donnée par l'équation (1.15) :

$$ 
\begin{gathered}
\mathbf{p}^{i_1, i_2} (\omega_1, \omega_2) = \sum_{i_1=1}^{n_1} \sum_{i_2=1}^{n_2} b_{i_1, i_2}^{i_1, i_2} (\omega_1, \omega_2) \mathbf{p}^{i_1, i_2} \\
1 \le j_1 \le ns_1, 1 \le j_2 \le ns_2, (\omega_1, \omega_2) \in [0,1]^2
\end{gathered}
\qquad (1.15)
$$

<figure><img src="image_42.png" /><figcaption>Figure 30 Surface spline composée de (2,2) patches.</figcaption></figure>

Les patches $S^{i_1, i_2}$ sont, en dimension 2, les éléments analogues aux segments $C^j$ d'une courbe spline. De plus, du fait de l'accroissement de la dimension, la notion de succession de segments splines se transpose en une notion de juxtaposition de patches (Figure 30).

<figure><img src="image_43.png" /><figcaption>Figure 31 Volume spline composé de (2,2,2) patches.</figcaption></figure>

L'expression mathématique d'un volume spline est déduite de la même démarche utilisée pour expliciter une surface spline et permet d'aboutir à la définition synthétisée par l'équation (1.16).

$$ 
\begin{gathered}
\mathbf{p}^{i_1, i_2, i_3} (\omega_1, \omega_2, \omega_3) = \sum_{i_1=1}^{n_1} \sum_{i_2=1}^{n_2} \sum_{i_3=1}^{n_3} b_{i_1, i_2, i_3}^{i_1, i_2, i_3} (\omega_1, \omega_2, \omega_3) \mathbf{p}^{i_1, i_2, i_3} \\
1 \le j_1 \le ns_1, 1 \le j_2 \le ns_2, 1 \le j_3 \le ns_3, (\omega_1, \omega_2, \omega_3) \in [0,1]^3
\end{gathered}
\qquad (1.16)
$$

On remarque une forte similarité entre les équations (1.2), (1.15) et (1.16) qui ne diffèrent que par la dimension intrinsèque de l'objet spline (courbe, surface ou volume) qui conditionne le système d'indexation et le nombre des sommes. Il est donc apparu pertinent de mettre en place une nouvelle expression de ces objets splines d'une part indépendante de la nature de leurs
fonctions de mélange mais aussi indépendante de leur dimension intrinsèque. De cette manière, je parlerai indifféremment de courbes, de surfaces et de volumes splines en utilisant une formulation mathématique unique. Pour rendre compte de cette unité d'expression, je nommerai désormais ces courbes, surfaces et volumes splines à l'aide du terme de **variété spline**. Le terme de variété, emprunté au vocabulaire de la géométrie différentielle [Berger 1987], désigne un ensemble de dimension *d* ayant certaines propriétés topologiques que je ne détaillerai pas dans ce document.

Afin de pouvoir expliciter l'expression mathématique unifiée d'une variété spline de dimension *d*, il me faut introduire la notion de **multi-indice**.

## 3.2 Multi-indice de dimension *d*

Un **multi-indice** est un *d*-uplet que l'on note de la façon suivante : $\mathbf{i}=(i_1,i_2,...,i_d) \in \mathbb{N}^d$. Parmi ces multi-indices, on en distingue deux particuliers dont je ferai grand usage : le multi-indice **unité** noté $\mathbf{1}=(1,1,...,1) \in \mathbb{N}^d$ et le multi-indice **élément neutre** $\mathbf{0}=(0,0,...,0) \in \mathbb{N}^d$. Je définis ici le **cardinal** d'un multi-indice de cette manière :

$$ \forall \mathbf{i} \in \mathbb{N}^d : \mathrm{Card}(\mathbf{i}) = i_1 \times i_2 \times \dots \times i_d \quad (1.17) $$

On munit $\mathbb{N}^d$, l'ensemble des multi-indices, d'une relation d'ordre partiel $\le$ :

$$ \forall \mathbf{i}, \mathbf{j} \in \mathbb{N}^d : \mathbf{i} \le \mathbf{j} \Leftrightarrow \forall k \in \{1,...,d\} : i_k \le j_k \quad (1.18) $$

## 3.3 Variété spline de dimension *d*

En s'appuyant sur les notations introduites pour les multi-indices, on donne la définition d'une variété spline de dimension *d*, unification des expressions mathématiques (1.2), (1.15) et (1.16), indépendante de la nature des fonctions de mélange et de la dimension intrinsèque de l'objet spline.

Soient un ensemble de points $\mathbf{p}^i \in \mathbb{R}^3$ avec $1 \le \mathbf{i} \le \mathbf{n}$, $\mathbf{n} \in \mathbb{N}^d$ le nombre de sous variétés et des fonctions de *d* variables réelles $b_i^j : [0,1]^d \to \mathbb{R}$ avec $1 \le \mathbf{i} \le \mathbf{n}, 1 \le \mathbf{j} \le \mathbf{n}$. On appelle **sous variétés splines**, les variétés $V^j$ dont la position des points est donnée par l'expression suivante :

$$ \mathbf{p}^j(\omega) = \sum_{1 \le i \le \mathbf{n}} b_i^j(\omega) \mathbf{p}^i \quad 1 \le \mathbf{j} \le \mathbf{n}, \omega \in [0,1]^d \quad (1.19) $$

La variété spline V correspond à la réunion des sous variétés $V^j$ :

$$ V = \bigcup_{1 \le j \le \mathbf{n}} V^j \quad (1.20) $$
Comme pour le cas des courbes splines, un point de la variété spline est identifié par son indice de sous-variété **j** et son indice paramétrique **ω**. On définit alors une fonction de position globale de la variété spline V dont l’expression est donnée par l’équation (1.21) :

$$
\mathbf{p} : \left\{ \begin{array}{l} \mathbb{N}_{\mathrm{ns}}^{d} \times [0,1]^{d} \rightarrow \mathbb{R}^{3} \\[1em] (\mathbf{j}, \omega) \mapsto \mathbf{p}(\mathbf{j}, \omega) = \mathbf{p}^{\mathrm{j}}(\omega) \end{array} \right. \quad \text{où} \quad \mathbb{N}_{\mathrm{ns}}^{d} = \{ \mathbf{j} \in \mathbb{N}^{d} : \mathbf{1} \le \mathbf{j} \le \mathbf{n}\} \qquad (1.21)
$$

<figure><img src="image_44.png" /><figcaption>Figure 32 Variété spline de dimension 1.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Figure 33 Variété spline de dimension 2.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Figure 34 Variété spline de dimension 3.</figcaption></figure>
Ce modèle géométrique de variété spline, générique d'une part vis-à-vis des fonctions de mélange mais aussi vis-à-vis de la dimension intrinsèque de l'objet, correspond au modèle mathématique que j'ai choisi pour représenter des corps déformables continus à degrés de liberté finis.

## 3.4 Coordonnées généralisées

Conformément au formalisme Lagrangien, une configuration géométrique cohérente d'une variété spline est caractérisée par un nombre fini de paramètres réels suffisants et indépendants. Ce *n*-uplet de réels, noté $\eta$ est appelé vecteur des **coordonnées généralisées** de la variété spline. Les coordonnées du vecteur $\eta$ correspondent aux degrés de liberté du système. Pratiquement, les coordonnées généralisées d'une variété spline sont les coordonnées de ses **n** points de contrôle. Dans le cadre d'une configuration étendue incluant le temps, le vecteur $\eta$ évolue donc dans un espace de configuration $\mathcal{U}$ de dimension $3n+1$ où $n = \text{Card}(\mathbf{n})$.

Dans la mesure où les points de contrôle de la variété spline sont indexés par des multi-indices de $\mathbb{N}^d$, il me faut mettre en place une fonction de transcodage $\alpha_n$ permettant d'établir la relation entre un multi-indice $i$ et l'indice entier $i$ correspondant aux coordonnées du point de contrôle $i$ dans le vecteur des coordonnées généralisées $\eta$. La fonction $\alpha_n$ est une bijection de l'ensemble $\mathbb{N}_n^d = \{i \in \mathbb{N}^d : 1 \le i \le n\}$ vers l'ensemble $\{1, \dots, n\}$ où $n = \text{Card}(\mathbf{n})$, $\mathbf{n}$ correspondant au « nombre » de points de contrôle de la variété spline. L'expression de la fonction de transcodage $\alpha_n$ pour les dimensions 1,2 et 3 est donnée par l'équation (1.22).

$$ \alpha_n(i) = \begin{cases} i_1 & \text{si } d=1 \\ i_1 + (i_2-1)n_1 & \text{si } d=2 \\ i_1 + (i_2-1)n_1 + (i_3-1)n_1 n_2 & \text{si } d=3 \end{cases} \quad (1.22) $$

Ainsi, les coordonnées des points de contrôle de la variété spline et les coordonnées du vecteur $\eta$ sont liées par l'expression (1.23) :

$$ \forall 1 \le i \le n: \eta_{\alpha_n(i+(k-1)n)} = p_k^i, \quad 1 \le k \le 3 \quad (1.23) $$

Le vecteur des coordonnées généralisées $\eta$ est composé du temps puis successivement des coordonnées $x$, $y$ et $z$ des points de contrôle. De manière plus synthétique, on écrit :

$$ \eta = \left( \underbrace{\eta_0}_{\text{temps}}, \underbrace{\eta_1, \dots, \eta_n}_{x}, \underbrace{\eta_{n+1}, \dots, \eta_{2n}}_{y}, \underbrace{\eta_{2n+1}, \dots, \eta_{3n}}_{z} \right)^T \quad (1.24) $$

## 3.5 Fonction de positionnement

A partir de la relation (1.23) je suis en mesure d'exprimer aisément la **fonction de positionnement** qui, à une configuration géométrique donnée par le vecteur des coordonnées généralisées $\eta$ et à un point de la variété repéré par son numéro de sous-variété $j$ et son indice paramétrique $\omega$, donne la position de ce point dans l'espace $\mathbb{R}^3$.
$$
\Phi : \begin{cases} \mathbb{U} \times \mathbb{N}_{\mathrm{ns}}^d \times [0.1]^d \to \mathbb{R}^3 \\ (\boldsymbol{\eta}, \mathbf{j}, \boldsymbol{\omega}) \mapsto \Phi(\boldsymbol{\eta}, \mathbf{j}, \boldsymbol{\omega}) \end{cases} \tag{1.25}
$$
$$
\text{où } \Phi_k(\boldsymbol{\eta}, \mathbf{j}, \boldsymbol{\omega}) = \sum_{1 \le i \le n} b_i^j(\boldsymbol{\omega}) \eta_{\alpha_n(i)+(k-1)n} \quad 1 \le k \le 3
$$

Au vu de l'équation (1.25), il apparaît que la fonction de positionnement $\Phi$ est une application linéaire du vecteur des coordonnées généralisées $\eta$. On peut donc associer à la fonction $\Phi$ une représentation matricielle $\mathbf{B}(\mathbf{j}, \boldsymbol{\omega})$ dont les propriétés structurelles sont détaillées dans l'équation

$$
\mathbf{B}(\mathbf{j}, \omega) = \begin{pmatrix} 0 & \mathbf{b}^j(\omega) & \cdots & \mathbf{0} \\ 0 & \vdots & \mathbf{b}^j(\omega) & \vdots \\ 0 & \mathbf{0} & \cdots & \mathbf{b}^j(\omega) \end{pmatrix} \tag{1.26}
$$

Le terme $\mathbf{b}^j(\omega)$ correspond à une matrice ligne composée de l'ensemble des fonctions de mélange :

$$
\mathbf{b}^j(\omega) = (b_{\alpha_n^{-1}(1)}^j(\omega) \cdots b_{\alpha_n^{-1}(n)}^j(\omega)) \quad (1.27)
$$

4 Codage informatique

Lorsque j'ai entrepris l'implantation informatique du modèle de variété spline, je voulais que le modèle informatique reflète cette idée d'unification présente dans le modèle mathématique, aussi bien vis-à-vis des fonctions de mélange que de la dimension intrinsèque de la variété. Ainsi, je fus naturellement amené à privilégier une méthode de programmation orientée objet afin de bénéficier des concepts de polymorphisme et de généricité.

En premier lieu, j’ai commencé par définir trois classes **Dimension** (Figure 35) présentant pour les dimensions 1, 2 et 3 des particularités fonctionnelles semblables :

□ une classe **E** codant un multi-indice *d*-dimensionnel. La classe **E** intervient dans la numérotation des sous-variétés $V^j$ et des points de contrôle $\mathbf{p}^i$ de la variété spline V.

□ une classe nichée **R** codant un vecteur réel *d*-dimensionnel. La classe **R** est utilisée pour représenter l'indice paramétrique **ω** identifiant une particule au sein d'une sous-variété $V^j$.

□ une classe nichée **ParcoursE** codant un « parcoureur » de valeurs multi-indicielles. La classe **ParcoursE**, munie, entre autres, de méthodes d'incrémentation et de décrémentation, permet de balayer un ensemble de multi-indices apparaissant dans une somme pour l'évaluation d'une position **p**(j,ω) par exemple.
<figure><img src="image_47.png" /><figcaption>Figure 35 Représentation schématique des trois classes **Dimension**.</figcaption></figure>

J'ai alors implanté une classe abstraite **FonctionMelange** paramétrée par l'un de ces types **Dimension** codant une fonction de mélange générique. La définition de la classe **FonctionMelange** s'appuyant sur un patron de classe paramétré par un type **Dimension** garantit la généricité de l'expression au vu de la dimension. Tandis que la nature abstraite de la classe permet d'unifier la prise en compte des différentes fonctions de mélange, déduites par héritage et polymorphisme. Ainsi, une classe concrète correspondant à une fonction de mélange donnée est obtenue via l'instantiation pour un type **Dimension** concret d'une classe dérivée de **FonctionMelange** (Figure 36).

<figure><img src="image_48.png" /><figcaption>Figure 36 Représentation schématique de l'architecture des classes **FonctionMelange**.</figcaption></figure>

Cette architecture de classe permet en l'occurrence, grâce au mécanisme d'héritage et de méthodes virtuelles, de procéder au sein des classes dérivées à des optimisations structurelles : chaque classe a la charge d'allouer l'espace mémoire qui lui est nécessaire pour stocker ses fonctions de mélange (4 polynômes de degré 3 pour une fonction de mélange de Kochanek, *m* polynômes de degré *m* pour une fonction de mélange de Bézier). D'autre part, la surcharge de méthodes virtuelles permet de générer du code spécifique exploitant les propriétés mathématiques des fonctions mélange. Ainsi, pour des fonctions de mélange définies à partir de polynômes, je suis en mesure de pré-calculer leurs dérivées de façon exacte. Ces dérivées interviennent dans le calcul de la dérivée de la fonction de positionnement, utilisée pour
déterminer des vecteurs tangents à la variété spline ou pour calculer l'abscisse curviligne en tout point de la variété, nécessaire à l'évaluation des longueurs.

Enfin, la classe **Spline** est l'objet informatique qui modélise une variété spline de dimension *d*. Pour les mêmes raisons énoncées pour les fonctions de mélange, la classe **Spline** est un patron de classe abstrait, paramétré par un type **Dimension**. Ainsi, une classe concrète correspondant à une variété spline de nature et de dimension données est obtenue via l'instantiation pour l'un des types **Dimension** concret d'une classe dérivée de **Spline** (Figure 37). Le choix d'une classe dérivée conditionne la nature des fonctions de mélange $b_i^j(\omega)$ utilisées ainsi que la relation entre le nombre de points de contrôle **n** et le nombre de sous-variétés splines **ns**.

<figure><img src="image_49.png" /><figcaption>Figure 37 Représentation schématique de l'architecture des classes **Spline**.</figcaption></figure>

La classe **Spline** est composée des champs suivants :

* Un multi-indice **n** correspondant au « nombre » de points de contrôle de la variété spline.

* Un multi-indice **ns** correspondant au « nombre » de sous-variétés de la variété spline.

* Une table de transcodage **alpha** modélisant la fonction de transcodage $\alpha_n$ introduite en section 2.4. A partir d'un multi-indice **i** de type **Dimension::E**, la table **alpha** permet d'accéder à l'indice **i** de la coordonnée *x* du point de contrôle **i** dans le vecteur des coordonnées généralisées **η**.

* Une référence à un objet **FonctionMelange** notée **b**.

La classe **Spline** est dotée des méthodes suivantes :

* Une méthode abstraite **parcours** qui, pour un numéro de sous-variété **j**, retourne un parcoureur des points de contrôle intervenant dans la définition de la sous-variété **j**. C'est à chaque classe dérivée, seule capable d'établir le lien entre points de contrôle et sous-variétés, que revient la responsabilité d'implanter cette méthode. La surcharge de
la méthode abstraite **parcours** assure donc que, dans tous les cas de figure, le parcours
des points de contrôle sera optimal.

□ Une méthode **position** qui, à partir d'un vecteur des coordonnées généralisées $\eta$ renseignant sur la configuration du réseau de points de contrôle, d'un numéro de sous-variété $j$ de type **Dimension::E** et d'un indice paramétrique $\omega$ de type **Dimension::R**, retourne la position dans $\mathbb{R}^3$ de la particule identifiée par le couple $(j, \omega)$ dans la configuration donnée par le vecteur $\eta$. Cette méthode correspond au codage informatique de la fonction de positionnement $\Phi$ introduite en (1.25). La méthode **position** s'appuie naturellement sur la méthode **parcours** pour l'évaluation de la fonction de positionnement.

□ Une méthode **dérivée** qui, à partir d'un vecteur des coordonnées généralisées $\eta$, d'un numéro de sous-variété $j$ de type **Dimension::E**, d'un indice paramétrique $\omega$ de type **Dimension::R** et d'un multi-indice $a$ **Dimension::E** renseignant sur l'ordre de dérivation, retourne la dérivée de la position dans $\mathbb{R}^3$ de la particule identifiée par le couple $(j, \omega)$. Cette méthode permet d'évaluer en l'occurrence des vecteurs tangents aux fibres matérielles de la variété spline, permettant de déduire la normale à un élément surfacique infinitésimal pour l'application d'une force de pression ou pour les besoins du rendu lors de la visualisation d'une variété spline.

Pour résumer, la classe Spline correspond au codage informatique de la structure géométrique d'une variété spline. On peut constater que les points de contrôle n'apparaissent pas dans la liste des champs de la classe. En effet, fidèles à la volonté de proposer un modèle informatique aussi proche que possible du modèle mathématique, je préserve la séparation de la description cinématique du système matériel. Ainsi, à partir de la donnée extérieure d'un jeu de points de contrôle stockés dans un vecteur unidimensionnel $\eta$, la classe Spline est en mesure d'évaluer la position de toute particule ainsi que des directions parallèles aux fibres principales de la variété dans cette configuration proposée.

A ce jour, j'ai implanté essentiellement des splines de nature cubique telles que les splines
cardinales, les splines de Catmull-Rom et les splines de Kochanek. Cette préférence pour les
splines cubiques a été motivée par les choix de Jean-Michel Nourrit lors de la réalisation d'un
modèle géométrique de tricots à base de courbes splines [Nourrit 1999].
# CINEMATIQUE ET DYNAMIQUE DES VARIETES SPLINES

## 1 Cinématique

Afin de simuler le comportement d'objets déformables continus, il me faut dans un premier temps formaliser mathématiquement le mouvement de ces corps modélisés à l'aide de variétés splines. Ce mouvement peut être perçu comme une succession de configurations géométriques instantanées indexées par le temps. Puisque la configuration d'une variété spline est complètement déterminée par les coordonnées de ses points de contrôle, j'introduis des points de contrôle mobiles $\mathbf{p}^i(t)$ pour caractériser le mouvement d'une variété spline. L'équation (1.19) du chapitre 2 devient alors, avec l'introduction du temps :

$$ \mathbf{p}(t, \mathbf{j}, \omega) = \sum_{1 \le i \le n} b_i^j(\omega) \mathbf{p}^i(t) \quad 1 \le \mathbf{j} \le \mathbf{n}, \omega \in [0,1]^d \quad (2.1) $$

Dans le cadre d'une description lagrangienne du mouvement, où chaque particule est identifiée par son numéro de sous-variété $\mathbf{j}$ et son indice paramétrique $\omega$, l'expression $\mathbf{p}(t, \mathbf{j}, \omega)$ de l'équation (2.1) correspond à la définition de la trajectoire particulaire. La cinématique d'une variété spline est synthétisée par l'ensemble de ces trajectoires particulaires. D'apparence indénombrable, cette cinématique peut se réduire à une finitude de « trajectoires généralisées » matérialisées par les trajectoires des points de contrôle. Cet aspect de la modélisation cinématique est développé en détail dans les sections de ce chapitre.

### 1.1 Fonction des configurations

J'introduis naturellement la fonction des configurations $\varepsilon$. Cette trajectoire généralisée, puisqu'elle synthétise l'évolution des coordonnées généralisées, est une courbe (paramétrée par le temps) de l'espace de configurations $\mathcal{U}$. L'équation (2.2) donne l'expression de la fonction des configurations déduite de l'équation (1.23) du chapitre 1.

$$ \varepsilon: \begin{cases} \mathbb{R} \to \mathbb{R}^{3n+1} & \text{où } \varepsilon_0(t)=t \\ t \mapsto \varepsilon(t) & \forall \mathbf{1} \le \mathbf{i} \le \mathbf{n}: \varepsilon_{\alpha_n(i)+(k-1)n}(t) = p_k^i(t) \quad 1 \le k \le 3 \end{cases} \quad (2.2) $$
En utilisant cette définition de la fonction des configurations, l'expression de la trajectoire
particulaire (2.1) d'un point de la variété spline, identifié par son numéro de sous-variété **j** et son
indice paramétrique **ω**, devient :

$$
\mathbf{p}(t, \mathbf{j}, \omega) = \Phi(\varepsilon(t), \mathbf{j}, \omega) \tag{2.3}
$$

La fonction Φ correspond à la fonction de positionnement définie dans le chapitre précédent.
On rappelle tout de même que la fonction de positionnement d'une variété spline est une
application linéaire des coordonnées généralisées dont la représentation matricielle **B**(j,ω) est
de la forme suivante.

$$
\mathbf{B}(\mathbf{j}, \omega) = \begin{pmatrix} 0 & \mathbf{b}^j(\omega) & \cdots & \mathbf{0} \\ 0 & \vdots & \mathbf{b}^j(\omega) & \vdots \\ 0 & \mathbf{0} & \cdots & \mathbf{b}^j(\omega) \end{pmatrix} \tag{2.4}
$$

Le terme **b**<sup>j</sup>(**ω**<sub>1</sub>) est une matrice ligne composée de l'ensemble des valeurs des fonctions de mélange en **ω** :

$$
\mathbf{b}^j(\omega) = (b_{\alpha_n^{-1}(1)}^j(\omega) \cdots b_{\alpha_n^{-1}(n)}^j(\omega)) \quad (2.5)
$$

1.2 Fonction de vitesse

La définition de la trajectoire particulaire (2.1) permet de déduire, par dérivation temporelle, l'expression de la vitesse particulaire d'un point de la variété spline identifié par le couple (j,ω).

$$
\mathbf{v}(t, \mathbf{j}, \omega) = \frac{d\mathbf{p}}{dt}(t, \mathbf{j}, \omega) = \frac{d}{dt} \left( \Phi(\boldsymbol{\varepsilon}(t), \mathbf{j}, \omega) \right) = \frac{\partial \Phi}{\partial \eta}(\boldsymbol{\varepsilon}(t), \mathbf{j}, \omega) \dot{\boldsymbol{\varepsilon}}(t) \quad (2.6)
$$

Comme la fonction de positionnement Φ est une application linéaire des coordonnées généralisées, elle est identique à sa jacobienne. Ainsi, l'expression de la vitesse particulaire devient :

$$
\mathbf{v}(t, \mathbf{j}, \omega) = \mathbf{B}(\mathbf{j}, \omega) \dot{\mathbf{\epsilon}}(t) \quad (2.7)
$$

La **fonction de vitesse** qui, à une configuration géométrique **η**, son taux de variation **ŷ** et une particule de la variété spline identifiée par le couple (**j**, **ω**), fait correspondre la vitesse de cette particule s'exprime alors sous la forme suivante :

$$
\hat{\Phi}: \begin{cases} U \times \hat{U} \times \mathbb{N}^d \times [0.1]^d \to \mathbb{R}^3 \\ (\eta, \hat{\eta}, \mathbf{j}, \omega) \mapsto \hat{\Phi}(\eta, \hat{\eta}, \mathbf{j}, \omega) = \mathbf{B}(\mathbf{j}, \omega) \hat{\eta} \end{cases} \tag{2.8}
$$
Il est important de mentionner que l'évaluation de la fonction de vitesse est indépendante du vecteur des coordonnées généralisées, ne reposant en fait que sur les taux de variation des coordonnées généralisées et de l'identifiant particulière (j, ω).

Evidemment, on peut écrire :

$$v(t, j, \omega) = \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), j, \omega) \quad (2.9)$$

## 2 Dynamique

La dynamique consiste en l'étude des causes de la modification du mouvement qui se manifestent sous la forme d'efforts. Plus encore, la dynamique a pour vocation de mettre en relation le mouvement et ses causes profondes. Cette relation est mathématiquement formalisée par ce que l'on nomme couramment l'équation du mouvement. J'ai choisi d'utiliser les équations de Lagrange pour décrire le mouvement des variétés splines.

Au fil de ce chapitre, je vais commencer par donner une définition de la répartition massique paramétrique d'une variété spline : répartition massique qui, bien que compatible avec la définition géométrique, n'est pas naturelle pour un physicien. Je montrerai alors comment, à partir d'une définition intuitive de la densité de masse par unité de longueur, je peux déduire l'expression de la répartition massique paramétrique d'une variété spline. Par la suite, je donnerai une expression de l'énergie cinétique spécifique aux variétés splines me permettant d'opérer une simplification formelle des équations de Lagrange afin d'aboutir à un système différentiel équivalent aux propriétés structurelles intéressantes en termes d'optimisation. Je terminerai ce chapitre par l'explicitation de quelques efforts externes et internes, en donnant systématiquement leur expression particulière dans le cadre des variétés splines.

### 2.1 Répartition massique

A l'image de la représentation continue des particules au sein de la variété spline, je désire donner une définition d'une répartition massique continue pour la variété. Ainsi, la densité de masse est une fonction du temps et du couple (j, ω) où j est le numéro de sous-variété et ω l'indice paramétrique de la particule.

$$\mu: \begin{cases} \mathbb{N}_{\text{ns}}^d \times [0,1]^d \times \mathbb{R} \rightarrow \mathbb{R} \\ (j, \omega, t) \mapsto \mu(j, \omega, t) \end{cases} \quad (2.10)$$

Ce choix est à distinguer des modélisations dans lesquelles des points caractéristiques de l'objet sont assimilés à des particules massiques. Ces modélisations engendrent une discrétisation de la masse totale, concentrée autour de points caractéristiques.

Mais cette définition de la densité de masse de nature paramétrique n'est pas intuitive. En effet, les physiciens ont pour habitude de considérer des répartitions massiques par unité de longueur, surface ou volume. Mon objectif est donc de déterminer, à partir de la donnée d'une densité de masse par unité métrique ρ(j, ω, t), une densité de masse par unité paramétrique μ(j, ω, t).
Je décide de me limiter d'abord au cas des courbes splines pour mettre en évidence la relation
qui lie la densité par unité de longueur $\rho(j,\omega,t)$ et la densité par unité paramétrique $\mu(j,\omega,t)$.
Pour cela j'introduis la définition de l'abscisse curviligne (2.11) qui permet de mesurer la
longueur d'un segment infinitésimal.

$$
l(j,\omega,t) = \int_{0}^{\omega} \| \hat{\partial}_{u} \mathbf{p}(j,u,t) \| du \quad ; \quad \hat{\partial}_{\omega} l(j,\omega,t) = \| \hat{\partial}_{\omega} \mathbf{p}(j,\omega,t) \| \qquad (2.11)
$$

<figure><img src="image_50.png" /><figcaption>Figure 38 Répartition massique par unité de longueur et répartition massique par unité paramétrique.</figcaption></figure>

En un point ω du segment j de la courbe, pour un élément matériel infinitésimal correspondant
à une variation paramétrique dω , on a la relation suivante entre les densités :

$$
\mu(j, \omega, t) d\omega = \rho(j, \omega, t) \hat{\partial}_{\omega} l(j, \omega, t) d\omega \quad (2.12)
$$

Cette relation est naturellement transposable au cas des surfaces (2.13) et des volumes (2.14).

<figure><img src="image_51.png" /><figcaption>Figure 39 Elément infinitésimal de surface.</figcaption></figure>

$$
\begin{gather*}
\mu(\mathbf{j}, \omega, t) d\omega = \rho(\mathbf{j}, \omega, t) ds(\mathbf{j}, \omega, t) \\
\text{où } ds(\mathbf{j}, \omega, t) = \| \hat{\partial}_{\omega_1} \mathbf{p}(\mathbf{j}, \omega, t) \wedge \hat{\partial}_{\omega_2} \mathbf{p}(\mathbf{j}, \omega, t) \| d\omega \tag{2.13}
\end{gather*}
$$
<figure><img src="image_52.png" /><figcaption>Figure 40 Elément infinitésimal de volume.</figcaption></figure>

$$ 
\begin{gathered}
\mu(\mathbf{j}, \omega, t)d\omega = \rho(\mathbf{j}, \omega, t)dv(\mathbf{j}, \omega, t) \\
\text{où } dv(\mathbf{j}, \omega, t) = (\partial_{\omega_1}\mathbf{p}(\mathbf{j}, \omega, t) \wedge \partial_{\omega_2}\mathbf{p}(\mathbf{j}, \omega, t)) \cdot \partial_{\omega_3}\mathbf{p}(\mathbf{j}, \omega, t)d\omega
\end{gathered}
\tag{2.14}
$$

De cette manière, je suis en mesure de définir une densité massique par unité paramétrique $\mu(\mathbf{j}, \omega, t)$ équivalente à une densité massique par unité métrique $\rho^0(\mathbf{j}, \omega, t)$ valide pour un état de référence caractérisé par un jeu de points de contrôle $\mathbf{p}^0$.

A partir de maintenant, on se place dans le cadre des variétés splines de dimension 3 qui correspond au cadre le plus général. Le résultat pourra ultérieurement être réinterprété en dimensions 1 et 2 de façon aisée.

<figure><img src="image_53.png" /><figcaption>Figure 41 Etat de référence d'une variété spline de dimension 3.</figcaption></figure>
Ainsi, l'état de référence de la variété spline est obtenu à partir d'un réseau régulier de points de contrôle. De forme parallélépipédique, ses arêtes sont alignées aux axes principaux et ont pour longueurs respectives $n_1 a_1, n_2 a_2$ et $n_3 a_3$ (Figure 41). La position d'un point identifié par le couple $(\mathbf{j}, \mathbf{\omega})$ dans cet état de référence est donnée par l'équation ci-dessous :

$$ \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}) = \sum_{1 \le i \le n} b_i^j(\mathbf{\omega}) \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + a_2 \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + a_3 \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix} \quad \mathbf{1} \le \mathbf{j} \le \mathbf{n}\mathbf{s}, \mathbf{\omega} \in [0,1]^d \quad (2.15) $$

Pour cet état de référence, l'équation (2.14) devient :

$$ \mu(\mathbf{j}, \mathbf{\omega}, t)d\mathbf{\omega} = \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}, t)d\mathbf{v}^0(\mathbf{j}, \mathbf{\omega}) \\ \text{où } d\mathbf{v}^0(\mathbf{j}, \mathbf{\omega}) = (\partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}) \wedge \partial_{\omega_2} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega})) \cdot \partial_{\omega_3} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}) d\mathbf{\omega} \quad (2.16) $$

Au vu de la configuration particulière de l'état de référence, donnée par les positions des points de contrôle $\mathbf{p}^0$, on a la relation suivante :

$$ \begin{aligned} d\mathbf{v}^0(\mathbf{j}, \mathbf{\omega}) &= \|\partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega})\| \|\partial_{\omega_2} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega})\| \|\partial_{\omega_3} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega})\| d\mathbf{\omega} \\ &= \partial_{\omega_1} l_1^0(\mathbf{j}, \mathbf{\omega}) \partial_{\omega_2} l_2^0(\mathbf{j}, \mathbf{\omega}) \partial_{\omega_3} l_3^0(\mathbf{j}, \mathbf{\omega}) d\mathbf{\omega} \end{aligned} \quad (2.17) $$

Les trois termes $\| \partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}) \| = (\partial_{\omega_3} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega})^2)^{1/2}$ correspondent aux abscisses curvilignes dans les directions principales pour l'état de référence défini en (2.15). Pour la commodité des calculs et ce, sans perte de généralité, je vais me contenter d'exprimer le terme $\| \partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}) \|$.

$$ \partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}) = \partial_{\omega_1} \sum_{1 \le i \le n} b_i^j(\mathbf{\omega}) \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + a_2 \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + a_3 \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix} \quad (2.18) $$

Je fais l'hypothèse supplémentaire que les fonctions de mélange $b_i^j(\mathbf{\omega})$ sont séparables. Cette hypothèse est peu restrictive dans la mesure où la propriété de séparabilité est largement utilisée pour construire des surfaces ou des volumes splines.

$$ b_i^j(\mathbf{\omega}) = b_{i_1}^{j_1}(\omega_1) b_{i_2}^{j_2}(\omega_2) b_{i_3}^{j_3}(\omega_3) \quad (2.19) $$

Par conséquent, l'équation (2.18) se transforme en :

$$ \partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \mathbf{\omega}) = \sum_{i_1=1}^{n_1} \sum_{i_2=1}^{n_2} \sum_{i_3=1}^{n_3} b_{i_1}^{j_1'}(\omega_1) b_{i_2}^{j_2}(\omega_2) b_{i_3}^{j_3}(\omega_3) \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + a_2 \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + a_3 \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix} \quad (2.20) $$
Après développement du terme de droite, on débouche sur une somme de trois termes :

$$
\begin{align}
\partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \omega) = & \sum_{i_1=1}^{n_1} a_1 i_1 b_{i_1}^{j_1'} (\omega_1) \sum_{i_2=1}^{n_2} b_{i_2}^{j_2} (\omega_2) \sum_{i_3=1}^{n_3} b_{i_3}^{j_3} (\omega_3) \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \nonumber \\
& + \sum_{i_2=1}^{n_2} a_2 i_2 b_{i_2}^{j_2} (\omega_2) \sum_{i_3=1}^{n_3} b_{i_3}^{j_3} (\omega_3) \sum_{i_1=1}^{n_1} b_{i_1}^{j_1'} (\omega_1) \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} \tag{2.21} \\
& + \sum_{i_3=1}^{n_3} a_3 i_3 b_{i_3}^{j_3} (\omega_3) \sum_{i_2=1}^{n_2} b_{i_2}^{j_2} (\omega_2) \sum_{i_1=1}^{n_1} b_{i_1}^{j_1'} (\omega_1) \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \nonumber
\end{align}
$$

On impose une nouvelle condition aux fonctions de mélange. Cette condition qui consiste en une hypothèse de normalité est, tout comme la séparabilité (2.19), largement répandue et ne constitue donc pas une restriction importante :

$$
\forall j, \omega : \sum_{i=1}^{n} b_i^j(\omega) = 1 \Rightarrow \sum_{i=1}^{n} b_i^{j'}(\omega) = 0 \quad (2.22)
$$

L'équation (2.21) se simplifie alors en un seul terme :

$$
\partial_{\omega_1} \mathbf{p}^0(\mathbf{j}, \omega) = \sum_{i_1=1}^{n_1} a_1 i_1 b_{i_1}^{j_1'}(\omega_1) \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \quad (2.23)
$$

L'abscisse curviligne $|| \hat{\partial}_{\omega_1} \mathbf{p}^0(\mathbf{j}, \omega) ||$ s'exprime finalement de la manière suivante :

$$
\| \hat{\partial}_{\omega_1} \mathbf{p}^0(\mathbf{j}, \omega) \| = \sum_{i_1=1}^{n_1} a_1 i_1 b_{i_1}^{j_1'}(\omega_1) \quad (2.24)
$$

Plus généralement, puisque le développement précédent est identique pour les coordonnées $\omega_2$ et $\omega_3$ :

$$
\| \hat{\partial}_{\omega_k} \mathbf{p}^0(\mathbf{j}, \omega) \| = \sum_{i=1}^{n_k} a_k i_k b_i^{j_k'}(\omega_k) \quad (2.25)
$$

Par conséquent l'équation (2.16) devient après calculs :

$$
\mu(\mathbf{j}, \omega, t) d\omega = \rho^0(\mathbf{j}, \omega, t) \sum_{i_1=1}^{n_1} a_1 i_1 b_{i_1}^{j_1'}(\omega_1) \sum_{i_2=1}^{n_2} a_2 i_2 b_{i_2}^{j_2'}(\omega_2) \sum_{i_3=1}^{n_3} a_3 i_3 b_{i_3}^{j_3'}(\omega_3) \quad (2.26)
$$
Si l'on considère une densité massique par unité métrique $\rho^0(j, \omega)$ constante vis à vis du temps, on débouche sur une expression simplifiée de la densité de masse paramétrique, invariante elle aussi vis à vis du temps.

$$ \mu : \left\{ \begin{array}{l} \mathbb{N}_{\text{ns}}^{d} \times [0,1]^{d} \rightarrow \mathbb{R} \\ (j, \omega) \mapsto \mu(j, \omega) \end{array} \right. \qquad (2.27) $$

Cette nouvelle définition paramétrique de la répartition massique suggère que la matière est attachée à sa position paramétrique $\omega$ au sein d'une sous-variété $j$. En d'autres termes, on suppose que la matière ne glisse pas à l'intérieur de la variété spline.

## 2.2 Energie cinétique, équations de Lagrange

Une fois la répartition massique explicitée, je suis en mesure de donner une expression de l'énergie cinétique spécifique aux variétés splines. L'énergie cinétique d'une variété spline est définie de manière classique comme suit :

$$ K(t) = \sum_{1 \le j \le \text{ns}} \frac{1}{2} \int_{[0,1]^d} \mu(j, \omega) v(t, j, \omega)^2 d\omega \qquad (2.28) $$

Je rappelle la définition de la fonction d'énergie cinétique, opérant un découplage entre les caractéristiques structurelles du système et sa configuration.

$$ \kappa : \left\{ \begin{array}{l} \mathbb{U} \times \hat{\mathbb{U}} \rightarrow \mathbb{R} \\ (\eta, \hat{\eta}) \mapsto \kappa(\eta, \hat{\eta}) = \frac{1}{2} \hat{\eta}^T \mathbf{M}(\eta) \hat{\eta} \end{array} \right. \qquad (2.29) $$

La matrice $\mathbf{M}(\eta)$ est la matrice des masses généralisées de la variété spline. L'équation (2.30) fournit l'expression particulière de $\mathbf{M}(\eta)$ pour une variété spline.

$$ \mathbf{M}(\eta) = \sum_{1 \le j \le \text{ns}} \int_{[0,1]^d} \mu(j, \omega) \frac{\partial \Phi}{\partial \eta}(\eta, \omega)^T \frac{\partial \Phi}{\partial \eta}(\eta, \omega) d\omega \qquad (2.30) $$

Comme je l'ai mentionné dans la section consacrée à la cinématique, la fonction de positionnement $\Phi$ d'une variété spline est une application linéaire des coordonnées généralisées. La matrice $\mathbf{B}(j, \omega)$ est donc indépendante des coordonnées généralisées $\eta$, conférant ainsi à la matrice des masses généralisées $\mathbf{M}(\eta)$ la propriété d'invariance. Pour refléter cette propriété, la matrice des masses généralisées sera désormais notée $\mathbf{M}$.
D'autre part, au vu de la structure par blocs de la matrice jacobienne de la fonction de
positionnement Φ, la matrice des masses généralisées hérite de cette structure par blocs
illustrée par l'équation suivante :

$$
\mathbf{M} = \begin{pmatrix} 0 & \cdots & \cdots & 0 \\ \vdots & \mathbf{m} & & \vdots \\ & \mathbf{m} & \ddots & \vdots \\ 0 & & \mathbf{m} & \end{pmatrix} \quad \text{avec} \quad \mathbf{m} \in M_n(\mathbb{R}) \qquad (2.31)
$$

Les trois blocs identiques de la diagonale, désignés par **m**, sont symétriques et s'expriment
ainsi :

$$
\mathbf{m} = \sum_{1 \le j \le n_s} \int_{[0,1]^d} \mu(\mathbf{j}, \omega) \mathbf{b}^i(\omega)^T \mathbf{b}^j(\omega) d\omega \quad (2.32)
$$

Ou, de façon équivalente, si l'on considère les éléments scalaires de la matrice **m** :

$$
m_{ik} = \sum_{1 \le j \le n_s} \int_{[0,1]^d} \mu(\mathbf{j}, \omega) b_{\alpha_n^{-1}(i)}^j(\omega) b_{\alpha_n^{-1}(k)}^j(\omega) d\omega \quad 1 \le i, k \le n \quad (2.33)
$$

La particularité structurelle de la matrice des masses généralisées **M** a d'intéressantes
conséquences sur l'expression du système différentiel équivalent aux équations de Lagrange. Je
rappelle d'abord brièvement l'expression de ce système différentiel. Les équations de Lagrange
ci-dessous,

$$
\frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i} (\varepsilon(t), \dot{\varepsilon}(t)) - \frac{\partial \kappa}{\partial \eta_i} (\varepsilon(t), \dot{\varepsilon}(t)) = Q_i(\varepsilon(t), \dot{\varepsilon}(t)) \quad 1 \le i \le n \quad (2.34)
$$

sont équivalentes au système différentiel suivant :

$$
\mathbf{M}^{**}(\boldsymbol{\varepsilon}(t))\ddot{\boldsymbol{\varepsilon}}(t) = \mathbf{Q}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) - \mathbf{I}\mathbf{n}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) \quad (2.35)
$$

Le vecteur In(ε(t), ˙ε(t)) correspondant au vecteur des termes inertiels. Dans la mesure où la
matrice des masses généralisées M d'une variété spline est invariante, le vecteur des termes
inertiels In(ε(t), ˙ε(t)) est nul. On débouche alors sur une expression simplifiée du système
differentiel conditionnant l'évolution d'une variété spline.

$$
\mathbf{M}^{**}\ddot{\mathbf{e}}(t) = \mathbf{Q}(\mathbf{e}(t), \dot{\mathbf{e}}(t)) \qquad (2.36)
$$

2.3 Efforts

Les efforts interviennent dans les équations de Lagrange par le biais du vecteur des forces
généralisées noté Q(η, ŷ). Je donne dans cette section l'expression de plusieurs natures
d'efforts en exploitant les particularités géométriques des variétés splines.
### 2.3.1 Force ponctuelle générique

Une particularité importante du modèle de variété spline réside dans la faculté d'identifier toute particule du milieu continu que constitue cette variété. Cette propriété permet d'appliquer une force ponctuelle donnée en un point quelconque identifié par son numéro de sous-variété **j** et son indice paramétrique $\omega$. Si $F(\eta, \hat{\eta})$ est un effort ponctuel appliqué en $(j, \omega)$, le vecteur des forces généralisées $Q(\eta, \hat{\eta})$ associé est de la forme :

$$Q_{i+(l-1)n}(\eta, \hat{\eta}) = F(\eta, \hat{\eta})^T \frac{\partial \Phi}{\partial \eta_{i+(l-1)n}}(\eta, j, \omega) \quad \forall 1 \le i \le n, k = 1,2,3 \quad (2.37)$$

De part la définition de la fonction de positionnement d'une variété spline, on est en mesure d'écrire :

$$\frac{\partial \Phi}{\partial \eta_{i+(l-1)n}}(\eta, j, \omega) = \sum_{m=1}^{n} b_{\alpha_n^{-1}(m)}^{j}(\omega) \delta_{mi} \begin{bmatrix} \delta_{1l} \\ \delta_{2l} \\ \delta_{3l} \end{bmatrix} = b_{\alpha_n^{-1}(i)}^{j}(\omega) \begin{bmatrix} \delta_{1l} \\ \delta_{2l} \\ \delta_{3l} \end{bmatrix} \quad (2.38)$$

L'équation (2.37) devient finalement,

$$Q_{i+(l-1)n}(\eta, \hat{\eta}) = b_{\alpha_n^{-1}(i)}^j(\omega) F(\eta, \hat{\eta})^T \begin{bmatrix} \delta_{1l} \\ \delta_{2l} \\ \delta_{3l} \end{bmatrix} = b_{\alpha_n^{-1}(i)}^j(\omega) F_l(\eta, \hat{\eta}) \quad \forall 1 \le i \le n, k = 1,2,3 \quad (2.39)$$

Je mets à profit cette habilité à définir des forces appliquées en des points quelconques de la variété spline pour construire un réseau de ressorts attachés en des points caractéristiques afin de modéliser les efforts internes d'élasticité. Ce choix, difficilement justifiable par son manque d'adéquation avec la réalité du matériau mais aussi par la discrétisation qu'il engendre, consiste en une première approche des phénomènes élastiques. Je propose dans le chapitre 3, la définition d'une énergie continue de déformation afin d'éliminer cette étape de discrétisation indésirable des efforts. Mais ces premiers travaux vers la définition d'un modèle d'efforts continus se limitent pour le moment au cas de courbes splines.

### 2.3.2 Densité d'efforts : amortissement

Afin de réduire les oscillations dues à la réponse des contraintes internes, j'introduis une densité d'efforts d'amortissement définie pour chaque particule de la variété. On ne peut pas à proprement parler de viscosité puisque les efforts considérés sont détachés de toute considération aérodynamique. Il s'agit en fait d'une densité d'efforts proportionnelle à l'opposé du champ des vitesses instantanées. Cette densité d'efforts appliqués en tout point matériel de coordonnées paramétriques $\omega$ au sein de chaque sous-variété **j** a la forme suivante :

$$F(\eta, \hat{\eta}, j, \omega) = -C \hat{\Phi}(\eta, \hat{\eta}, j, \omega) \quad \forall j \in \mathbb{N}_{ns}^d, \forall \omega \in [0,1]^d \quad (2.40)$$

$C$ est le coefficient d'amortissement.
En utilisant l'expression de la fonction de vitesse (2.8), on a aussi :

$$F_l(\eta, \hat{\eta}, j, \omega) = -C \sum_{1 \le k \le n} b_{\alpha_n^{-1}(k)}^j(\omega) \hat{\eta}_{k+(l-1)n}(t) \quad k = 1, 2, 3 \quad (2.41)$$

L'équation (2.38) permet de donner l'expression des coordonnées du vecteur des forces généralisées d'amortissement,

$$ 
\begin{align*}
Q_{i+ln}(\eta, \hat{\eta}) &= \sum_{1 \le j \le ns} \int_{[0,1]^d} F(\eta, \hat{\eta}, j, \omega)^T \frac{\partial \Phi}{\partial \eta_{i+ln}}(\eta, j, \omega) d\omega \\
&= \sum_{1 \le j \le ns} \int_{[0,1]^d} F(\eta, \hat{\eta}, j, \omega)^T \left( b_{\alpha_n^{-1}(i)}^j(\omega) \begin{bmatrix} \delta_{1l} \\ \delta_{2l} \\ \delta_{3l} \end{bmatrix} \right) d\omega \quad \forall 1 \le i \le n, k = 1, 2, 3 \tag{2.42} \\
&= \sum_{1 \le j \le ns} \int_{[0,1]^d} -C \sum_{k=1}^n b_{\alpha_n^{-1}(k)}^j(\omega) \hat{\eta}_{k+(l-1)n}(t) b_{\alpha_n^{-1}(i)}^j(\omega) d\omega \\
&= -C \sum_{k=1}^n \hat{\eta}_{k+(l-1)n} \sum_{1 \le j \le ns} \int_{[0,1]^d} b_{\alpha_n^{-1}(k)}^j(\omega) b_{\alpha_n^{-1}(i)}^j(\omega) d\omega
\end{align*}
$$

Les 3n expressions des coordonnées du vecteur $Q(\eta, \hat{\eta})$ peuvent être synthétisées de la sorte :

$$Q(\eta, \hat{\eta}) = \Theta \hat{\eta} \quad (2.43)$$

$\Theta$ est une matrice carrée de dimension 3n, constituée de trois blocs identiques constants sur sa diagonale.

$$ \mathbf{\Theta} = \begin{pmatrix} \theta \\ \theta \\ \theta \end{pmatrix} \quad \text{et} \quad \theta_{ik} = -C \sum_{1 \le j \le ns} \int_{[0,1]^d} b_{\alpha_n^{-1}(i)}^j(\omega) b_{\alpha_n^{-1}(k)}^j(\omega) d\omega \quad (2.44) $$

### 2.3.3 Energie potentielle : gravitation

Les efforts de gravitation sont déduits de l'énergie potentielle de gravitation définie comme suit :

$$E(\eta) = \sum_{1 \le j \le ns} \int_{[0,1]^d} \mu(j, \omega) \Phi(\eta, j, \omega)^T g d\omega \quad (2.45)$$

le vecteur $g$ est dirigé selon la verticale gravitationnelle et sa norme $g = \|g\|$ est la constante de gravitation locale.

Le vecteur des forces généralisées de gravitation correspond à la dérivée partielle de l'énergie potentielle de gravitation par rapport au vecteur des coordonnées généralisées $\eta$.

$$Q(\eta) = -\frac{\partial E}{\partial \eta^*} = -\sum_{1 \le j \le ns} \int_{[0,1]^d} \mu(j, \omega) \frac{\partial \Phi}{\partial \eta^*}(\eta, j, \omega)^T g d\omega \quad (2.46)$$
On utilise alors l'expression particulière de la jacobienne de la fonction de positionnement (2.4) pour donner une nouvelle formulation du vecteur des forces généralisées de gravitation **Q**.

$$Q = \begin{pmatrix} \Pi \\ \Pi \\ \Pi \end{pmatrix} g \quad \text{avec} \quad \Pi = - \sum_{1 \le j \le ns} \int_{[0,1]^d} \mu(j, \omega) b^i(\omega)^T d\omega \quad (2.47)$$

## 3 Codage informatique

L'introduction de nouveaux éléments mathématiques relatifs à la dynamique des variétés splines a d'importantes conséquences sur l'implantation informatique qui se matérialisent par la construction de nouvelles classes. Ainsi apparaît une classe abstraite **Densite** paramétrée par l'un des types **Dimension**, codant la répartition massique paramétrique définie en (2.10). Fidèle au mécanisme d'héritage utilisé pour les classes précédemment définies, une classe concrète correspondant à une répartition massique particulière est obtenue via l'instantiation pour un type **Dimension** concret d'une classe dérivée de **Densite** (Figure 42).

La classe **Densite** a la charge d'évaluer tous les termes faisant intervenir la répartition massique, à savoir les coefficients de la matrice des masses généralisées **M** (2.31) et ceux de la matrice **Π** des efforts de gravitation (2.47). Au vu de l'expression des coefficients à évaluer, on dote la classe **Densité** de deux méthodes abstraites :

* Une méthode **calculeIntegraleBji** qui, à partir d'une référence à un objet de type **FonctionMelange** et de deux multi-indices **i** et **j** de type **Dimension**::E, calcule le terme intégral (2.48) intervenant dans l'expression des efforts de gravitation.

$$\int_{[0,1]^d} \mu(j, \omega) b_i^j(\omega) d\omega \quad (2.48)$$

* Une méthode **calculeIntegraleBjiBjk** qui, à partir d'une référence à un objet de type **FonctionMelange** et de trois multi-indices **i**, **j** et **k** de type **Dimension**::E, calcule le terme intégral (2.49) intervenant dans l'expression de la matrice des masses généralisées.

$$\int_{[0,1]^d} \mu(j, \omega) b_i^j(\omega) b_k^j(\omega) d\omega \quad (2.49)$$

<figure><img src="image_54.png" /><figcaption>Figure 42 Représentation schématique de l'architecture des classes **Densite**.</figcaption></figure>
La définition de ces deux méthodes dans la classe **Densite** implique l'ajout de méthodes abstraites à la classe **FonctionMelange** de même nature :

□ Une méthode **calculeIntegraleBji** qui, à partir de deux multi-indices i et j de type **Dimension** **::E**, calcule le terme intégral (2.50).

$$ \int_{[0,1]^d} \mathbf{b}_i^j(\omega) d\omega \qquad (2.50) $$

□ Une méthode **calculeIntegraleBjiBjk** qui, à partir de trois multi-indices i, j et k de type **Dimension** **::E**, calcule le terme intégral (2.51).

$$ \int_{[0,1]^d} \mathbf{b}_i^j(\omega) \mathbf{b}_k^j(\omega) d\omega \qquad (2.51) $$

Cette logique permet en l'occurrence de calculer de façon exacte les termes (2.48) et (2.49) pour une variété spline de densité massique constante dont les fonction de mélange sont des polynômes, qui constitue un type de variété spline fréquemment utilisé.

Enfin, la classe **Spline** s'enrichit de nouveaux champs tels que :

□ Une matrice **M** correspondant à la matrice des masses généralisées de la variété spline. Il est important de mentionner que, pour une densité massique invariante et des fonctions de mélange polynomiales, la matrice **M** est constante et peut être évaluée de façon formellement exacte lors de l'étape d'initialisation. De plus, on peut pré-calculer sa décomposition afin d'accélérer la résolution du système linéaire, équivalent aux équations de Lagrange.

□ Une matrice **Theta** correspondant à la matrice d'amortissement définie en (2.44). Sous les mêmes conditions que précédemment, cette matrice est constante et peut donc être pré-calculée de façon formellement exacte durant l'étape d'initialisation.

□ Une matrice **Pi** correspondant à la matrice gravitationnelle définie en (2.47). Si le vecteur **g** est invariant, le vecteur des forces généralisées de gravitation **Q** déduit du produit de **Π** et de **g** peut être pré-calculé lui aussi.

De cette manière, les efforts de gravitation et d'amortissement sont pris en compte de façon efficace au sein de la classe **Spline**. Par contre, Il existe un type d'effort dont je n'ai pas encore parlé : en effet, les efforts issus des contraintes internes au sein d'un corps déformable sont d'une importance primordiale et leur modélisation soulève nombre de problèmes. Voilà pourquoi, j'ai choisi d'y consacrer le chapitre suivant.

# 4 Exemples de variétés splines

## 4.1 Plongeoir

Le système matériel considéré est modélisé à l'aide d'une variété spline cubique d'Hermite de dimension 3 constituée de (9,6,5) points de contrôle et de (6,3,2) sous-variétés splines. L'élasticité interne du plongeoir est garantie par la mise en œuvre d'un réseau de ressorts au nombre de (2,2,2) par sous-variétés, régulièrement répartis.
L'extrémité gauche du plongeoir demeure figée grâce à l'ajout de 8 contraintes de point fixe sur les points de contrôle de la face gauche de la variété spline.

Enfin, le plongeoir est soumis à des efforts de gravité, d'amortissement mais aussi à trois forces ponctuelles appliquées sur le dessus de l'extrémité droite de l'objet durant un court instant, lui conférant un mouvement oscillatoire amorti. Les images de la Figure 43 sont des clichés pris à des instants successifs issus de l'animation répondant à cette définition

<figure><img src="image_55.png" /><figcaption>Figure 43 Images issues d'une animation reproduisant le mouvement oscillatoire amorti d'un plongeoir.</figcaption></figure>

## 4.2 Cube élastique

Le deuxième exemple de système matériel est constitué de deux corps déformables. D'une part, la corde est modélisée à l'aide d'une courbe spline cubique d'Hermite composée de 7 points de contrôle. Deux ressorts par segment spline viennent garantir la cohésion de la corde. D'autre part, le cube est représenté par une variété spline cubique d'Hermite de dimension 3 constituée de (5,5,5) points de contrôle et (2,2,2) sous-variétés splines. De même, le cube est doté de ressorts, au nombre de (2,2) par sous-variété, pour modéliser les efforts internes d'élasticité.
Un contrainte de distance entre une extrémité de la corde et un coin du cube permet de lier ces deux corps déformables au fil de l'animation. Et une contrainte de point fixe imposée à l'autre extrémité en fait un point d'attache.

Le système composé de la corde et du cube est uniquement soumis aux efforts de gravité, d'amortissement et d'élasticité interne (Figure 44).

<figure><img src="image_56.png" /><figcaption>Figure 44 Images issues d'une animation reproduisant le mouvement d'un cube déformable suspendu par un de ses coins au bout d'une corde.</figcaption></figure>

## 4.3 Balançoire

Ce dernier exemple rassemble les trois natures de corps déformables actuellement implantés dans le moteur d'animation dynamique. En effet, la balançoire est composée de :

* 4 cordes modélisées à partir de courbes splines cubiques de Catmull-Rom, composées de 7 points de contrôle et de 4 segments splines. Elles sont dotées de 2 ressorts par segment pour modéliser les efforts internes d'élasticité.
□ Un dossier modélisé par une variété spline cubique d’Hermite de dimension 2, constituée de (5,5) points de contrôle et de (2,2) sous-variétés. Elle est aussi dotée de ressorts, au nombre de (2,2) par sous-variété.

□ Un siège modélisé par une variété spline cubique d’Hermite de dimension 3, constituée de (6,5,4) points de contrôle et de (3,2,1) sous-variétés, dotée de (2,2,2) ressorts par sous-variété.

Des contraintes de distance sont imposées entre les quatre extrémités inférieures des cordes et les quatre coins supérieurs du siège. Des contraintes de même nature sont imposées aux quatre coins du dossier pour qu’il demeure fixé aux cordes. Enfin, quatre contraintes de point fixe sont appliquées aux extrémités supérieures des cordes pour les attacher à la structure qui n’est pas un système matériel mais un objet inerte de l’environnement.

Les éléments de la balançoire sont uniquement soumis aux efforts de gravité, d’amortissement et d’élasticité interne (Figure 45).

<figure><img src="image_57.png" /><figcaption>Figure 45 Images issues d'une animation reproduisant le mouvement d'une balançoire.</figcaption></figure>
# ENERGIE CONTINUE DE DEFORMATION DES VARIETES SPLINES

## 1 Introduction

Jusqu'à présent, pour modéliser l'élasticité interne d'une variété spline, j'utilisais des ressorts attachés en des points caractéristiques constituant un réseau discret au sein de la variété. Bien que les modèles masses-ressorts, toujours populaires, présentent certains avantages [Baraff 1998] [Bourguignon 2000] [Volino 2000], je souhaitais uniformiser le modèle mécanique continu en introduisant une représentation continue des efforts internes d'élasticité. Je me suis alors orienté vers la définition d'une énergie continue de déformation conforme à la théorie de l'élasticité. Dans la mesure où les corps déformables que je cherche à animer peuvent être sujets à de grands déplacements, je me positionne hors du champ d'application des résultats de l'élasticité linéaire, abondamment utilisés en animation dynamique [Gourret 1989] [Collier 1991] [Eischen 1996] [Bro-Nielsen 1996]. En effet, la linéarité des contraintes vis-à-vis des déformations garantit une efficacité calculatoire permettant d'envisager des applications temps-réel [Cotin 1997] mais restreintes à l'animation de corps déformables soumis à de petits déplacements. J'ai alors dirigé mes recherches vers la théorie de l'élasticité non linéaire et plus particulièrement, vers le modèle de St Venant-Kirchhoff. Le contenu de ce chapitre a fait l'objet d'une publication internationale au 12<sup>th</sup> Eurographics Workshop on Computer Animation and Simulation [Nocent 2001c].

Je commence par une synthèse des définitions et résultats majeurs de la théorie de l'élasticité en déplacements finis. Pour de plus amples détails, j' invite le lecteur à compulser les ouvrages de Ciarlet [Ciarlet 1986] et Souchet [Souchet 2001] ainsi que les précieux polycopiés de cours de Garrigues [Garrigues 2000a, 2000b]. J'expose ensuite la manière dont j' applique ces résultats théoriques au cas particulier des variétés splines. Je tiens à mentionner que ce travail s' inscrit dans un projet à long terme visant à tirer parti des nombreux résultats de la théorie l'élasticité pour l' animation dynamique des variétés splines de dimension quelconque. Ainsi, les résultats d' ores et déjà disponibles consistent en une première approche restreinte à la modélisation du comportement purement élastique des courbes splines. Cette restriction provient en fait des considérations pragmatiques dictées par le cadre applicatif qui me concerne, à savoir la
simulation de fils tricotés. Enfin, je présente quelques résultats numériques afin de motiver le choix d'une modélisation continue des contraintes internes plutôt qu'une modélisation discrète à base de ressorts.

## 2 Théorie de l'élasticité

La mécanique des milieux continus permet, entre autres, de décrire l'évolution des déformations d'un corps déformable soumis à des contraintes internes et externes. La relation entre déformations et contraintes internes peut prendre des formes multiples telles que l'élasticité, la viscoélasticité, la plasticité et bien d'autres encore. Comme je l'ai justifié dans l'introduction, je choisis de restreindre, dans cette première approche, mon étude à la modélisation de l'élasticité. La théorie de l'élasticité repose en fait sur les hypothèses suivantes :

* Les déformations subies par le milieu déformable sont supposées réversibles.

* Les contraintes internes générées par les changements de forme du corps sont fonction des déformations.

Pour être en mesure d'utiliser les résultats de la théorie de l'élasticité, il me faut préciser les définitions de notions telles que les **déformations**, les **contraintes internes** et d'enfin expliciter la relation qui les lient : cette relation porte le nom de **loi de comportement**.

### 2.1 Déformations

On distingue classiquement deux approches distinctes en mécanique des milieux continus :

* La description **eulerienne** qui s'appuie sur l'étude du champ de vitesses instantanées exprimant au cours du temps et en tout point **p** de l'espace, la vitesse de la particule occupant la position **p** en cet instant (on ignore d'ailleurs de quelle particule il s'agit)

* La description **lagrangienne** qui s'appuie sur l'étude des trajectoires particulaires en supposant l'existence d'une configuration de référence qui permet une identification de chacune des particules.

Dans le même état d'esprit, une quantité scalaire, vectorielle ou tensorielle est qualifiée d'eulerienne si elle est fonction d'un champ de vitesses instantanées. Et cette même quantité est qualifiée de lagrangienne si elle est fonction de l'identifiant particulaire et du temps.

Conformément aux choix théoriques que j'ai été amené à opérer préalablement (Chapitre 1 : Justification des équations de Lagrange), je privilégie une approche lagrangienne pour la description des déformations.

#### 2.1.1 Déformations en description lagrangienne

Pour une description lagrangienne du mouvement, chaque particule constitutive du système matériel est identifiée par sa position **x**º dans une configuration de référence. Ainsi, la position temporelle **x**(t) de cette particule est fonction de sa position **x**º dans la configuration de référence et du temps t.

$$ \mathbf{x}(t) = \mathbf{p}(\mathbf{x}^0, t) \qquad (2.52) $$
Les coordonnées de la position de référence $\mathbf{x}^0 = (x_1^0, x_2^0, x_3^0) \in \mathbb{R}^3$ et le temps $t \in \mathbb{R}$ sont communément appelés les **variables de Lagrange**, tandis que les coordonnées de la position instantanée $\mathbf{x}(t) = (x_1(t), x_2(t), x_3(t)) \in \mathbb{R}^3$ de la particule sont appelées les **inconnues de Lagrange**.

<figure><img src="image_58.png" /><figcaption>Figure 46 Interprétation géométrique du transformé et du transporté.</figcaption></figure>

L'étude de la déformation d'un milieu continu consiste à comparer les positions particulières dans la configuration instantanée avec les positions des mêmes particules dans la configuration de référence. Cette idée est mathématiquement formalisée par l'équation (2.53) en considérant deux particules identifiées par $\mathbf{x}^0$ et $\mathbf{y}^0$ voisines dans la configuration de référence (Figure 46).

$$ \underbrace{\mathbf{y}(t) - \mathbf{x}(t)}_{\text{transformé de } \mathbf{y}^0 - \mathbf{x}^0} = \mathbf{p}(\mathbf{y}^0, t) - \mathbf{p}(\mathbf{x}^0, t) = \underbrace{\frac{\partial \mathbf{p}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t)(\mathbf{y}^0 - \mathbf{x}^0)}_{\text{transporté de } \mathbf{y}^0 - \mathbf{x}^0} + \|\mathbf{y}^0 - \mathbf{x}^0\| R(\mathbf{y}^0 - \mathbf{x}^0, t) \quad (2.53) $$

le terme $R(\mathbf{y}^0 - \mathbf{x}^0, t)$ correspond au reste du développement de Taylor.

Le **transformé** est le vecteur liant les positions instantanées des particules identifiées par $\mathbf{x}^0$ et $\mathbf{y}^0$. Le **transporté** est le vecteur tangent à l'arc de courbe constitué par les positions instantanées des particules appartenant au segment délimité par $\mathbf{x}^0$ et $\mathbf{y}^0$ dans la configuration de référence. Pour des particules $\mathbf{x}^0$ et $\mathbf{y}^0$ suffisamment proches, le **transformé** et le **transporté** sont confondus, faisant de la jacobienne de la position temporelle l'élément qui relie une direction de la configuration de référence avec la direction « déformée » correspondante de la configuration instantanée.

On définit alors $\mathcal{F}(\mathbf{x}^0, t)$, le **tenseur gradient de la transformation** (i.e. transformation linéaire tangente) de la façon suivante :

$$ \mathcal{F}(\mathbf{x}^0, t) = \frac{\partial \mathbf{p}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) = \operatorname{grad}_{\mathbf{x}^0} \mathbf{p}(\mathbf{x}^0, t) \qquad (2.54) $$

Le tenseur gradient de la transformation $\mathcal{F}(\mathbf{x}^0, t)$ exprime, dans la configuration instantanée, la direction d'une fibre matérielle autour de la particule $\mathbf{x}^0$ à partir de la direction de cette même fibre dans la configuration de référence (Figure 47).
<figure><img src="image_59.png" /><figcaption>Figure 47 Interprétation géométrique du tenseur gradient de la transformation.</figcaption></figure>

La dilatation d'une fibre matérielle ou le glissement de deux directions de fibres matérielles $\boldsymbol{\alpha}(t)$ et $\boldsymbol{\beta}(t)$ est caractérisé par leur produit scalaire de $\boldsymbol{\alpha}(t) \cdot \boldsymbol{\beta}(t)$. Si l'on considère que $\boldsymbol{\alpha}(t) = \mathcal{F}(\mathbf{x}^0, t)\boldsymbol{\alpha}^0$ et $\boldsymbol{\beta}(t) = \mathcal{F}(\mathbf{x}^0, t)\boldsymbol{\beta}^0$, on obtient alors :

$$ 
\begin{align} 
\boldsymbol{\alpha}(t) \cdot \boldsymbol{\beta}(t) &= (\mathcal{F}(\mathbf{x}^0, t)\boldsymbol{\alpha}^0)^T (\mathcal{F}(\mathbf{x}^0, t)\boldsymbol{\beta}^0) \nonumber \\
&= \boldsymbol{\alpha}^{0T} \mathcal{F}(\mathbf{x}^0, t)^T \mathcal{F}(\mathbf{x}^0, t)\boldsymbol{\beta}^0 \tag{2.55} \\
&= \boldsymbol{\alpha}^{0T} C(\mathbf{x}^0, t)\boldsymbol{\beta}^0 \quad \text{avec} \quad C(\mathbf{x}^0, t) = \mathcal{F}(\mathbf{x}^0, t)^T \mathcal{F}(\mathbf{x}^0, t) \nonumber
\end{align} 
$$

le tenseur de Cauchy-Green $C(\mathbf{x}^0, t)$ est le tenseur des dilatations. En effet, il permet d'exprimer la dilatation et le glissement, au voisinage de la particule $\mathbf{x}^0$ dans la configuration instantanée, de deux directions unitaires orthogonales $\boldsymbol{\alpha}^0$ et $\boldsymbol{\beta}^0$ de la configuration de référence.

Dilatation dans la direction unitaire $\boldsymbol{\alpha}^0$ :

$$ \lambda(\mathbf{x}^0, t)(\boldsymbol{\alpha}^0) = \sqrt{\boldsymbol{\alpha}^{0T} C(\mathbf{x}^0, t) \boldsymbol{\alpha}^0} \quad (2.56) $$

Glissement de deux directions unitaires orthogonales $\boldsymbol{\alpha}^0$ et $\boldsymbol{\beta}^0$ :

$$ \theta(\mathbf{x}^0, t)(\boldsymbol{\alpha}^0, \boldsymbol{\beta}^0) = \arccos\left(\frac{\boldsymbol{\alpha}^{0T} C(\mathbf{x}^0, t)\boldsymbol{\beta}^0}{\lambda(\mathbf{x}^0, t)(\boldsymbol{\alpha}^0)\lambda(\mathbf{x}^0, t)(\boldsymbol{\beta}^0)}\right) \quad (2.57) $$

Le tenseur de Cauchy-Green est un **tenseur lagrangien** puisqu'il est fonction de l'identifiant $\mathbf{x}^0$ de la particule dans la configuration de référence et du temps $t$.

Finalement, la déformation autour de la particule $\mathbf{x}^0$ dans la configuration instantanée est déterminée par la différence $\boldsymbol{\alpha}(t) \cdot \boldsymbol{\beta}(t) - \boldsymbol{\alpha}^0 \cdot \boldsymbol{\beta}^0$.
S'appuyant sur l'équation (2.55), on obtient alors :

$$
\begin{align*}
\boldsymbol{\alpha}(t) \cdot \boldsymbol{\beta}(t) - \boldsymbol{\alpha}^0 \cdot \boldsymbol{\beta}^0 &= \boldsymbol{\alpha}^{0T} C(\mathbf{x}^0, t) \boldsymbol{\beta}^0 - \boldsymbol{\alpha}^{0T} \mathcal{I} \boldsymbol{\beta}^0 \\
&= \boldsymbol{\alpha}^{0T} (C(\mathbf{x}^0, t) - \mathcal{I}) \boldsymbol{\beta}^0 \tag{2.58} \\
&= \boldsymbol{\alpha}^{0T} 2\mathcal{E}(\mathbf{x}^0, t) \boldsymbol{\beta}^0 \quad \text{avec} \quad \mathcal{E}(\mathbf{x}^0, t) = \frac{1}{2} (C(\mathbf{x}^0, t) - \mathcal{I})
\end{align*}
$$

$\mathcal{I}$ est le tenseur identité.

le tenseur de **Green-Lagrange** $\mathcal{E}(\mathbf{x}^0, t)$ est le tenseur des déformations. Tout comme le tenseur de Cauchy-Green, le tenseur de Green-Lagrange est un tenseur lagrangien. Ce tenseur est sensible aux variations de forme du système matériel mais il est invariant par transformation rigide (i.e. transformation isométrique). En effet, si l'on considère une fonction de position temporelle isométrique, on a :

$$
\forall \boldsymbol{\alpha}^0, \boldsymbol{\beta}^0 : \boldsymbol{\alpha}(t) \cdot \boldsymbol{\beta}(t) = \boldsymbol{\alpha}^0 \cdot \boldsymbol{\beta}^0 \Rightarrow C(\mathbf{x}^0, t) = \mathcal{I} \Rightarrow \mathcal{E}(\mathbf{x}^0, t) = 0 \quad (2.59)
$$

2.1.2 Elasticité linéaire

Au lieu de caractériser la déformation à partir de la position **p**(*x*⁰, *t*) qu'une particule *x*⁰
occupe à l'instant *t*, on peut choisir de la caractériser par son déplacement **u**(*x*⁰, *t*) depuis sa
position de référence. On définit le **déplacement** d'une particule de la façon suivante :

$$
\mathbf{u}(\mathbf{x}^0, t) = \mathbf{p}(\mathbf{x}^0, t) - \mathbf{x}^0 \quad (2.60)
$$

L'expression du tenseur gradient de transformation $\mathcal{F}(\mathbf{x}^0, t)$ devient alors :

$$
\mathcal{F}(\mathbf{x}^0, t) = \frac{\partial \mathbf{p}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) = \frac{\partial}{\partial \mathbf{x}^0}(\mathbf{x}^0 + \mathbf{u}(\mathbf{x}^0, t)) = \mathcal{I} + \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) \quad (2.61)
$$

En s'appuyant sur l'équation (2.61), on peut donner une nouvelle formulation du tenseur de
Green-Lagrange à partir du déplacement **u**(*x*⁰, *t*).

$$
\begin{align*}
\mathcal{E}(\mathbf{x}^0, t) &= \frac{1}{2} (C(\mathbf{x}^0, t) - \mathcal{I}) \\
&= \frac{1}{2} \left( \left( \mathcal{I} + \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) \right)^T \left( \mathcal{I} + \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) \right) - \mathcal{I} \right) \tag{2.62} \\
&= \frac{1}{2} \left( \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t)^T + \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) + \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t)^T \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) \right)
\end{align*}
$$
Finalement,

$$
\begin{equation} \tag{2.63}
\begin{split}
\mathcal{E}(\mathbf{x}^0, t) &= \epsilon(\mathbf{x}^0, t) + \frac{1}{2} \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t)^T \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) \\
\text{avec } \epsilon(\mathbf{x}^0, t) &= \frac{1}{2} \left( \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t)^T + \frac{\partial \mathbf{u}}{\partial \mathbf{x}^0}(\mathbf{x}^0, t) \right)
\end{split}
\end{equation}
$$

Le tenseur de **Cauchy** $\epsilon(\mathbf{x}^0, t)$ correspond à la partie linéaire du tenseur de Green-Lagrange $\mathcal{E}(\mathbf{x}^0, t)$. Dans le cadre des petits déplacements, le tenseur de Cauchy constitue une approximation tout à fait satisfaisante du tenseur de Green-Lagrange. La nature linéaire du tenseur de Cauchy garantit une plus grande rapidité dans son évaluation, justifiant son succès par rapport au tenseur de Green-Lagrange.

Sous cette hypothèse des petits déplacements où $\mathbf{x}(t) \approx \mathbf{x}^0$, le tenseur $\epsilon(\mathbf{x}^0, t)$ de Cauchy pourtant de nature lagrangienne peut être perçu comme un tenseur eulerien.

Alors que le tenseur de Green-Lagrange est invariant par transformation rigide, le tenseur de Cauchy quant à lui ne l'est pas. Cette propriété restreint donc le champ d'utilisation du tenseur de Cauchy à des systèmes matériels soumis à des petits déplacements autour de la configuration de référence. Des tentatives d'utilisation de ce tenseur pour des systèmes matériels sujets à des déplacements finis se révèlent infructueuses comme l'ont démontré Debunne [Debunne 2000] et Picinbono [Picinbono 2001].

3 Contraintes

Les contraintes correspondent aux efforts appliqués au sein du corps déformable en réponse à sa
déformation. Pour un découpage arbitraire du corps Ω en deux sous-corps Ω₁, Ω₂ séparés par
la frontière ∂Ω, les efforts appliqués par Ω₁ à Ω₂ sur un élément surfacique infinitésimal dS
de la frontière ∂Ω, au voisinage d'un point matériel x(t), peuvent être déduits de la normale n
à cette surface. Ils sont classiquement mesurés par le tenseur des contraintes de Cauchy qui,
pour une position x(t) de la configuration instantanée, permet de déterminer l'effort F(x(t))
appliqué à un petit élément de surface dS autour de x(t) à partir de sa normale n (Figure 48).
Plus synthétiquement, on peut écrire :

$$
F(\mathbf{x}(t)) = \sigma(\mathbf{x}(t)) \mathbf{n} \, dS \quad (2.64)
$$

Le tenseur $\sigma(\mathbf{x}(t))$ des contraintes de Cauchy est de **nature eulerienne** dans la mesure où il dépend d'une position de la configuration instantanée sans considération d'identification particulaire et qu'il s'applique à une normale de la configuration instantanée.
<figure><img src="image_60.png" /><figcaption>Figure 48 Interprétation géométrique du tenseur des contraintes de Cauchy.</figcaption></figure>

Puisque j'ai caractérisé les déformations d'un milieu continu à l'aide de tenseurs lagrangiens, j'ai recours à la transformée de Piola [Ciarlet 1986] [Garrigues 2000b] afin de déduire l'expression d'un tenseur de pseudo-contraintes lagrangien. Ce tenseur $S(\mathbf{x}^0, t)$ déduit de la transformée de Piola du tenseur des contraintes de Cauchy porte le nom de **premier tenseur des contraintes de Piola-Kirchhoff**.

$$S(\mathbf{x}^0, t) = \det(\mathcal{F}(\mathbf{x}^0, t)) \sigma(\mathbf{x}(t)) (\mathcal{F}(\mathbf{x}^0, t))^T^{-1} \quad \text{avec} \quad \mathbf{x}(t) = \mathbf{p}(\mathbf{x}^0, t) \quad (2.65)$$

Afin de satisfaire à la contrainte de symétrie des tenseurs de déformations et de contraintes, on définit le **second tenseur des contraintes de Piola-Kirchhoff** de la manière suivante :

$$\Sigma(\mathbf{x}^0, t) = \mathcal{F}(\mathbf{x}^0, t)^{-1} S(\mathbf{x}^0, t) \quad (2.66)$$

# 4 Loi de comportement

La loi de comportement permet d'établir le lien entre le tenseur des déformations et le tenseur des contraintes internes issues des modifications de forme du milieu continu. Il apparaît clairement que la nature ainsi que les propriétés structurelles du matériau considéré conditionnent fortement l'expression de la loi de comportement. Je vais donc énumérer plusieurs types de matériaux en donnant pour chacun d'eux la forme de la loi de comportement associée. Je m'attarderai plus particulièrement sur les matériaux de St Venant-Kirchhoff dont la loi de comportement, en dépit de sa simplicité, est tout à fait adaptée à mon besoin d'une réponse non linéaire des fils du tricot.

## 4.1 Matériaux élastiques

Un matériau est **élastique** si le tenseur des contraintes $\Sigma(\mathbf{x}^0, t)$ dépend uniquement de la position $\mathbf{x}^0$ de la particule et du tenseur gradient de la transformation $\mathcal{F}(\mathbf{x}^0, t)$. Cela équivaut à dire qu'il existe une application $\hat{\Sigma}: \Omega \times M_3^+(\mathbb{R}) \to M_3(\mathbb{R})$ telle que :

$$\Sigma(\mathbf{x}^0, t) = \hat{\Sigma}(\mathbf{x}^0, \mathcal{F}(\mathbf{x}^0, t)) \quad (2.67)$$

Pour un matériau élastique, le tenseur des contraintes ne dépend donc que des dérivées premières de la fonction de position. Parmi les matériaux élastiques, on distingue une sous-
famille particulière de matériaux qualifiés d'homogènes. Un matériau élastique est **homogène** si le tenseur des contraintes est indépendant du point considéré. De manière équivalente, il existe une application $\hat{\Sigma}: M_3^+(\mathbb{R}) \to M_3(\mathbb{R})$ telle que :

$$ \Sigma(\mathbf{x}^0, t) = \hat{\Sigma}(\mathcal{F}(\mathbf{x}^0, t)) \qquad (2.68) $$

Un principe général en physique affirme que toute quantité observable de nature intrinsèque doit être indépendante de la base orthonormée dans laquelle elle est calculée. Ce principe m' amène à formuler l'axiome suivant.

### Axiome de l'indifférence matérielle [Ciarlet 1986] [Souchet 2001]

Le vecteur des contraintes doit être indépendant de la base orthonormée choisie pour calculer le tenseur des contraintes à l'aide de la réponse du matériau

Une conséquence de cet axiome consiste à dire que l'on peut exprimer le tenseur des contraintes comme une fonction du tenseur de Cauchy-Green $C(\mathbf{x}^0, t)$. En effet, le tenseur $C(\mathbf{x}^0, t)$ correspond au tenseur des dilatations dont les valeurs sont indépendantes de la base orthonormée dans laquelle il est exprimé. On définit alors une application $\tilde{\Sigma}: \Omega \times S_3^>(\mathbb{R}) \to S_3(\mathbb{R})$ vérifiant la relation suivante :

$$ \hat{\Sigma}(\mathbf{x}^0, \mathcal{F}(\mathbf{x}^0, t)) = \tilde{\Sigma}(\mathbf{x}^0, \mathcal{F}(\mathbf{x}^0, t)^T \mathcal{F}(\mathbf{x}^0, t)) = \tilde{\Sigma}(\mathbf{x}^0, C(\mathbf{x}^0, t)) \qquad (2.69) $$

## 4.2 Matériaux élastiques isotropes

Un matériau élastique est **isotrope** si le tenseur des contraintes est invariant par rotation, c'est-à-dire :

$$ \hat{\Sigma}(\mathbf{x}^0, \mathcal{F}(\mathbf{x}^0, t)) = \hat{\Sigma}(\mathbf{x}^0, \mathcal{F}(\mathbf{x}^0, t) R) \quad \forall R \in O_3(\mathbb{R}) \qquad (2.70) $$

Cette condition signifie que la réponse du matériau est identique dans toutes les directions autour de la particule identifiée par $\mathbf{x}^0$. A partir de cette hypothèse, le théorème de représentation de Rivlin-Erickson [Ciarlet 1986] permet de déduire l'expression générale de la loi de comportement associée à un matériau élastique isotrope.

$$ \tilde{\Sigma}(\mathbf{x}^0, C(\mathbf{x}^0, t)) = \beta_1(l_c)\mathcal{I} + \beta_2(l_c)C(\mathbf{x}^0, t) + \beta_3(l_c)C(\mathbf{x}^0, t)^2 \qquad (2.71) \\ \text{où } l_c = (l_1(C(\mathbf{x}^0, t)), l_2(C(\mathbf{x}^0, t)), l_3(C(\mathbf{x}^0, t))) $$

$\beta_1, \beta_2$ et $\beta_3$ sont des fonctions à valeurs réelles, le vecteur $l_c$ est composé des **invariants principaux** du tenseur de dilatation de Cauchy-Green. L'expression de ces invariants principaux est donnée dans l'équation (2.72).

$$ \begin{aligned} l_1(C(\mathbf{x}^0, t)) &= \operatorname{tr} C(\mathbf{x}^0, t) \\ l_2(C(\mathbf{x}^0, t)) &= \frac{1}{2} \left\{ \left( \operatorname{tr} C(\mathbf{x}^0, t) \right)^2 - \operatorname{tr} C(\mathbf{x}^0, t)^2 \right\} \qquad (2.72) \\ l_3(C(\mathbf{x}^0, t)) &= \det C \end{aligned} $$
En exploitant la relation $C(\mathbf{x}^0, t) = \mathcal{I} + 2\mathcal{E}(\mathbf{x}^0, t)$, on donne une expression des invariants principaux de $C(\mathbf{x}^0, t)$ en fonction du tenseur de Green-Lagrange $\mathcal{E}(\mathbf{x}^0, t)$.

$$ 
\begin{align} 
l_1(C(\mathbf{x}^0, t)) &= 3 + 2\operatorname{tr}\mathcal{E}(\mathbf{x}^0, t) \nonumber \\ 
l_2(C(\mathbf{x}^0, t)) &= 3 + 4\operatorname{tr}\mathcal{E}(\mathbf{x}^0, t) + 2\left\{\left(\operatorname{tr}\mathcal{E}(\mathbf{x}^0, t)\right)^2 - \operatorname{tr}\mathcal{E}(\mathbf{x}^0, t)^2\right\} \tag{2.73} \\ 
l_3(C(\mathbf{x}^0, t)) &= \det\left(2\mathcal{E}(\mathbf{x}^0, t) + \mathcal{I}\right) \nonumber 
\end{align} 
$$

Ainsi, la décomposition de la loi de comportement d'un matériau élastique isotrope énoncée en (2.71) devient :

$$ \tilde{\Sigma}(\mathbf{x}^0, C(\mathbf{x}^0, t)) = -\pi(\mathbf{x}^0)\mathcal{I} + \lambda(\mathbf{x}^0)(\operatorname{tr}\mathcal{E}(\mathbf{x}^0, t))\mathcal{I} + 2\mu(\mathbf{x}^0)\mathcal{E}(\mathbf{x}^0, t) + o(\mathcal{E}^2(\mathbf{x}^0, t), \mathbf{x}^0) \quad (2.74) $$

où $\pi$, $\lambda$ et $\mu$ sont des fonctions de la position $\mathbf{x}^0$ dans la configuration de référence, caractéristiques du matériau. Dans le cas particulier d'un **matériau élastique homogène isotrope**, ces fonctions deviennent invariantes. On définit alors une nouvelle loi de comportement $\tilde{\Sigma}: S_3(\mathbb{R}) \to S_3(\mathbb{R})$ fonction du tenseur des déformations $\mathcal{E}(\mathbf{x}^0, t)$ et de deux scalaires $\lambda$ et $\mu$.

$$ \tilde{\Sigma}(\mathcal{E}(\mathbf{x}^0, t)) = \tilde{\Sigma}(\mathcal{I} + 2\mathcal{E}(\mathbf{x}^0, t)) = \lambda \operatorname{tr}(\mathcal{E}(\mathbf{x}^0, t)\mathcal{I}) + 2\mu\mathcal{E}(\mathbf{x}^0, t) + o(\mathcal{E}^2(\mathbf{x}^0, t)) \quad (2.75) $$

Les grandeurs scalaires $\lambda$ et $\mu$ sont les **constantes de Lamé**. Ces deux grandeurs caractérisent les propriétés du matériau et peuvent être déduites de l'expression du **module de Young E** et du **coefficient de Poisson v** par l'équation suivante :

$$ \lambda = \frac{E\nu}{(1 - 2\nu)(1 + \nu)} \qquad \mu = \frac{E}{2(1 + \nu)} \quad (2.76) $$

## 4.3 Matériaux de St Venant-Kirchhoff

Un **matériau de St Venant-Kirchhoff** est un matériau élastique homogène isotrope dont la loi de comportement correspond à la partie linéaire en $\mathcal{E}(\mathbf{x}^0, t)$ de (2.75)

$$ \tilde{\Sigma}(\mathcal{E}(\mathbf{x}^0, t)) = \tilde{\Sigma}(\mathcal{I} + 2\mathcal{E}(\mathbf{x}^0, t)) = \lambda \operatorname{tr}(\mathcal{E}(\mathbf{x}^0, t)\mathcal{I}) + 2\mu\mathcal{E}(\mathbf{x}^0, t) \quad (2.77) $$

A ce titre, les matériaux de St Venant-Kirchhoff sont, dans le cadre de l'élasticité en déplacements finis, les matériaux les plus étudiés.

## 4.4 Matériaux hyper élastiques

Un matériau est **hyper élastique** s'il existe une fonction d'énergie potentielle de déformation dont le tenseur des contraintes dérive. Dans le cadre d'une description lagrangienne des déformations et des contraintes, on considère une énergie potentielle $e(\mathbf{x}^0, t)$ dont le second tenseur de Piola-Kirchoff (2.66) dérive.
Ainsi,

$$ \Sigma(\mathbf{x}^0, t) = \frac{\partial e}{\partial \mathcal{E}}(\mathbf{x}^0, t) \qquad (2.78) $$

$\mathcal{E}$ représente le tenseur des déformations.

On définit alors l'**énergie totale de déformation** en sommant les contributions des énergies de déformations $e(\mathbf{x}^0, t)$ pour toutes les particules $\mathbf{x}^0$ du domaine $\Omega^0$ que constitue le milieu continu dans sa configuration de référence.

$$ e^{\text{totale}}(t) = \int_{\Omega^0} e(\mathbf{x}^0, t) d\mathbf{x}^0 \qquad (2.79) $$

A noter qu'un **matériau de St Venant-Kirchhoff est un matériau hyper élastique**. A partir de l'expression du tenseur des contraintes $\Sigma(\mathbf{x}^0, t)$ donnée en (2.77), on déduit aisément l'expression de l'énergie potentielle de déformation.

$$ e(\mathbf{x}^0, t) = \frac{\lambda}{2} \text{tr}(\mathcal{E}(\mathbf{x}^0, t))^2 + \mu \text{tr}(\mathcal{E}(\mathbf{x}^0, t)^2) \qquad (2.80) $$

L'énergie potentielle de déformation pour un matériau de St Venant-Kirchhoff peut aussi être obtenue à partir du produit contracté du tenseur des déformations $\mathcal{E}(\mathbf{x}^0, t)$ et du tenseur des contraintes $\Sigma(\mathbf{x}^0, t)$.

$$ e(\mathbf{x}^0, t) = \frac{1}{2} \mathcal{E}(\mathbf{x}^0, t) : \Sigma(\mathbf{x}^0, t) \qquad (2.81) $$

La loi de comportement associée à un matériau de St Venant-Kirchhoff, de nature lagrangienne puisqu'elle lie le tenseur de Green-Lagrange et le second tenseur de Piola-Kirchhoff, permet de calculer les contraintes internes imposées à un milieu continu à partir de ses déformations. Elle demeure valable pour des corps matériels sujets à des déplacements finis dans la mesure où le tenseur de Green-Lagrange est invariant par transformation rigide et donc seulement sensible aux changements de forme du corps.

## 4.5 Loi de Hooke

La loi de Hooke est similaire à la loi de comportement d'un matériau de St Venant-Kirchhoff à ceci près que le tenseur de déformations considéré est le tenseur $\epsilon(\mathbf{x}^0, t)$ de Cauchy (2.63). Cette loi de comportement est donc restreinte à l'étude des corps matériels sujets à des petits déplacements. Dans ce cas de figure, où la position instantanée est presque confondue avec la position de référence, le tenseur de Cauchy est considéré comme une quantité eulerienne. Le tenseur de contraintes déduit de cette expression de la déformation n'est autre que le tenseur des contraintes de Cauchy (2.64).
On a alors :

$$
\sigma(\mathbf{x}(t)) = \lambda \operatorname{tr}(\epsilon(\mathbf{x}(t))) \mathcal{I} + 2\mu \epsilon(\mathbf{x}(t)) \quad (2.82)
$$

5 Application aux courbes splines matérielles

Une courbe spline matérielle est un système mécanique curviligne sans épaisseur. Cette
hypothèse mathématique consistant en une idéalisation du phénomène réel conduit à construire
des objets sans aucune existence physique. Afin de pouvoir utiliser les résultats précédents, qui
reposent sur la nature volumique du milieu continu, ils me faut considérer désormais des
courbes splines matérielles dotées d'une épaisseur infinitésimale. Ainsi on considère, autour de
chaque point de la courbe, une section matérielle orthogonale à la courbe dans sa configuration
de référence. On suppose de plus que cette section demeure orthogonale à la courbe en tout
instant *t* et ne subit que des déformations de nature elliptique. Ces hypothèses semblent
raisonnables dans la mesure où l'épaisseur de la courbe, et donc le rayon moyen de la section,
sont supposés négligeables. Autour de chaque point de la courbe de position **p**(*j*,*ω*,*t*), on
construit un référentiel orthonormal (u¹(*j*,*ω*,*t*), u²(*j*,*ω*,*t*), u³(*j*,*ω*,*t*)) de sorte que son
premier axe u¹(*j*,*ω*,*t*) soit tangent à la courbe, c'est-à-dire que les vecteurs u¹(*j*,*ω*,*t*) et
∂_ω **p**(*j*,*ω*,*t*) soient colinéaires (Figure 49).

<figure><img src="image_61.png" /><figcaption>Figure 49 Section d'une courbe spline matérielle avec son référentiel associé</figcaption></figure>

Finalement, la position d'un point materiel au sein de cette courbe spline matérielle avec
épaisseur, identifiée par un numéro de segment *j* et un multi-indice *ω* = (*ω*₁, *ω*₂, *ω*₃) est donnée
par:

$$
\mathbf{x}(t) = \mathbf{p}(j, \omega_1, t) + \omega_2 a_2(j, \omega_1) \mathbf{u}^2(j, \omega_1, t) + \omega_3 a_3(j, \omega_1) \mathbf{u}^3(j, \omega_1, t) \quad (2.83)
$$

Les termes $a_2(j,\omega_1)$ et $a_3(j,\omega_1)$ représentent respectivement les facteurs de dilatation linéaire au sein de la section matérielle dans les directions $\mathbf{u}^2(j,\omega_1,t)$ et $\mathbf{u}^3(j,\omega_1,t)$. Ces facteurs permettent la modélisation des phénomènes de striction dans le plan de la section.
<figure><img src="image_62.png" /><figcaption>Figure 50 Etat de repos d'une courbe spline matérielle.</figcaption></figure>

La position de référence de chaque particule est obtenue à partir de l'état de repos de la courbe spline matérielle. Cet état de repos correspond à une courbe spline rectiligne, parallèle à l'axe Ox dont les segments ont des longueurs au repos fixées (Figure 50).

$$ \mathbf{x}^0 = (\mathbf{p}_1^0(j, \omega_1) \quad \omega_2 \quad \omega_3)^T = \left(\sum_{i=1}^n b_i^j(\omega_1) \mathbf{p}_1^{0i} \quad \omega_2 \quad \omega_3\right)^T \quad (2.84) $$

<figure><img src="image_63.png" /><figcaption>Figure 51 Construction du tenseur gradient de transformation à partir des fonctions paramétriques de position.</figcaption></figure>

Comme dans le modèle de Terzopoulos [Terzopoulos 1987], chaque particule est identifiée par sa coordonnée paramétrique $\omega$ et non directement par sa position $\mathbf{x}^0$ dans la configuration de
référence. Terzopoulos choisit de mesurer la déformation subie par le corps en faisant la différence des tenseurs métriques $G^0(\omega)$ et $G(\omega, t)$, quantités exprimées dans l'espace paramétrique donc sans grandeur physique. Pour construire un tenseur gradient de transformation $\mathcal{F}(\mathbf{x}^0, t)$ dont l'expression soit conforme à celle proposée par la théorie de l'élasticité, il faut être en mesure d'exprimer la dérivée partielle de la position instantanée $\mathbf{x}(t)$ par rapport à la position de référence $\mathbf{x}^0$. Pour cela, j'introduis l'inverse de la dérivée partielle de $\mathbf{x}^0$ par rapport à $\omega$ (Figure 51) dans l'expression de $\mathcal{F}(\mathbf{x}^0, t)$. Cet artifice de calcul permet de retomber dans les rails de la théorie de l'élasticité en définissant le tenseur gradient de transformation comme une fonction de la configuration de référence, permettant ainsi d'utiliser d'une manière classique les caractéristiques du matériau comme paramètres naturelles des contraintes internes.

$$ \mathcal{F}(\mathbf{x}^0, t) = \frac{\partial \mathbf{x}}{\partial \mathbf{x}^0}(t) = \frac{\partial \mathbf{x}}{\partial \omega}(t) \frac{\partial \omega}{\partial \mathbf{x}^0} = \frac{\partial \mathbf{x}}{\partial \omega}(t) \left( \frac{\partial \mathbf{x}^0}{\partial \omega} \right)^{-1} \quad (2.85) $$

D'autre part,

$$ \frac{\partial \mathbf{x}}{\partial \omega}(t) = \left( \frac{\partial l}{\partial \omega_1} \mathbf{u}^1 + \omega_2 \frac{\partial (a_2 \mathbf{u}^2)}{\partial \omega_1} + \omega_3 \frac{\partial (a_3 \mathbf{u}^3)}{\partial \omega_1} \quad a_2 \mathbf{u}^2 \quad a_3 \mathbf{u}^3 \right) \quad \text{où} \quad \partial_{\omega_1} l = \| \partial_{\omega_1} \mathbf{p} \| \quad (2.86) $$

En supposant l'épaisseur infinitésimale, les paramètres $\omega_2$ et $\omega_3$ sont négligeables devant $\partial_{\omega_1} l(j, \omega_1, t)$. A partir de maintenant, on considère que $\omega_2 = \omega_3 = 0$. Par conséquent, les termes définis précédemment ne dépendent désormais plus que d'un seul paramètre scalaire $\omega$ (auparavant noté $\omega_1$). Ainsi, l'équation (2.86) se simplifie de la manière suivante :

$$ \frac{\partial \mathbf{x}}{\partial \omega}(t) = \left( \frac{\partial l}{\partial \omega_1}(j, \omega, t) \mathbf{u}^1(j, \omega, t) \quad a_2(j, \omega, t) \mathbf{u}^2(j, \omega, t) \quad a_3(j, \omega, t) \mathbf{u}^3(j, \omega, t) \right) (2.87) $$

L'équation (2.87) permet d'aboutir à une forme simplifiée de l'expression du tenseur gradient de la transformation $\mathcal{F}(\mathbf{x}^0, t)$, en effet on a :

$$ \begin{gathered} \mathcal{F}(\mathbf{x}^0, t) = \mathcal{U}(j, \omega, t) \mathcal{A}(j, \omega, t) \\ \mathcal{U}(j, \omega, t) = (\mathbf{u}^1 \quad \mathbf{u}^2 \quad \mathbf{u}^3) \text{ et } \mathcal{A}(j, \omega, t) = \begin{pmatrix} \gamma(j, \omega, t) & & \\ & a_2(j, \omega, t) & \\ & & a_3(j, \omega, t) \end{pmatrix} \end{gathered} \quad (2.88) $$

Le scalaire $\gamma(j, \omega, t) = \partial_{\omega_1} l(j, \omega, t) / \partial_{\omega_1} l^0(j, \omega)$ exprime le facteur de dilatation relative dans la direction tangente à la courbe spline matérielle. Naturellement, l'expression du tenseur de Cauchy-Green $C(\mathbf{x}^0, t)$ devient :

$$ C(\mathbf{x}^0, t) = \mathcal{A}(j, \omega, t)^T \mathcal{U}(j, \omega, t)^T \mathcal{U}(j, \omega, t) \mathcal{A}(j, \omega, t) \quad (2.89) $$
Dans la mesure où la base des vecteurs $(\mathbf{u}^1, \mathbf{u}^2, \mathbf{u}^3)$ est orthonormale et le tenseur $\mathcal{A}(j, \omega, t)$ est diagonal, l'équation (2.89) se simplifie en :

$$C(\mathbf{x}^0, t) = \mathcal{A}(j, \omega, t)^T \mathcal{A}(j, \omega, t) = \mathcal{A}(j, \omega, t)^2 \quad (2.90)$$

Finalement, le tenseur de Green-Lagrange $\mathcal{E}(\mathbf{x}^0, t)$ est de la forme suivante.

$$\mathcal{E}(\mathbf{x}^0, t) = \frac{1}{2} \begin{pmatrix} \gamma(j, \omega, t)^2 - 1 & & \\ & a_2(j, \omega, t)^2 - 1 & \\ & & a_3(j, \omega, t)^2 - 1 \end{pmatrix} \quad (2.91)$$

Les deux facteurs de dilatation radiale $a_2(j, \omega, t)$ et $a_3(j, \omega, t)$ sont inconnus. Deux équations supplémentaires sont donc nécessaires afin de pouvoir expliciter complètement le tenseur des déformations $\mathcal{E}(\mathbf{x}^0, t)$ et le tenseur des contraintes $\Sigma(\mathbf{x}^0, t)$ qui en découle.

S'appuyant sur une hypothèse couramment formulée dans la théorie des coques minces, je suppose qu'aucun effort n'est appliqué à la surface de l'objet. La section étant supposée toujours perpendiculaire à la courbe, un vecteur normal $\mathbf{n}$ est forcément de la forme suivante $\mathbf{n} = \alpha_2 \mathbf{u}^2 + \alpha_3 \mathbf{u}^3$. Conformément à l'hypothèse de nullité des efforts radiaux, on a :

$$\sigma(\mathbf{x}(t), t)\mathbf{n} = 0 \quad \forall \alpha_2, \alpha_3 : \mathbf{n} = \alpha_2 \mathbf{u}^2 + \alpha_3 \mathbf{u}^3 \quad (2.92)$$

En raison de l'épaisseur infinitésimale de la courbe spline matérielle, la propriété (2.92) peut être étendue à toute particule intérieure. On donne une expression alternative (2.93) du tenseur des contraintes de Cauchy $\sigma(\mathbf{x}(t))$, dépendant explicitement du second tenseur de Piola-Kirchhoff, afin d'introduire de nouvelles équations faisant intervenir les coefficients de $\Sigma(\mathbf{x}^0, t)$.

$$\sigma(\mathbf{x}(t)) = \det(\mathcal{F}(\mathbf{x}^0, t))^{-1} \mathcal{F}(\mathbf{x}^0, t) \Sigma(\mathbf{x}^0, t) \mathcal{F}(\mathbf{x}^0, t)^T \quad (2.93)$$

L'équation (2.92) devient alors,

$$(\mathcal{F}(\mathbf{x}^0, t) \Sigma(\mathbf{x}^0, t) \mathcal{F}(\mathbf{x}^0, t)^T) \mathbf{n} = 0 \quad \forall \alpha_2, \alpha_3 : \mathbf{n} = \alpha_2 \mathbf{u}^2 + \alpha_3 \mathbf{u}^3 \quad (2.94)$$

En s'appuyant sur la définition du tenseur gradient de transformation d'une variété spline, on développe l'expression de l'équation (2.94) pour aboutir à la formulation suivante :

$$\mathcal{U}(\mathbf{x}^0, t) \mathcal{A}(\mathbf{x}^0, t) \Sigma(\mathbf{x}^0, t) \mathcal{A}(\mathbf{x}^0, t)^T (\mathcal{U}(\mathbf{x}^0, t)^T \mathbf{n}) = 0 \\ \text{avec } \mathcal{U}(\mathbf{x}^0, t)^T \mathbf{n} = \begin{pmatrix} 0 \\ \alpha_2 \\ \alpha_3 \end{pmatrix} \quad \forall \alpha_2, \alpha_3 : \mathbf{n} = \alpha_2 \mathbf{u}^2 + \alpha_3 \mathbf{u}^3 \quad (2.95)$$
Enfin, puisque les tenseurs $\mathcal{A}(\mathbf{x}^0, t)$ et $\Sigma(\mathbf{x}^0, t)$ sont diagonaux, on débouche sur l'expression simplifiée ci-dessous.

$$ \mathcal{A}^2(\mathbf{x}^0, t) \Sigma(\mathbf{x}^0, t) \begin{pmatrix} 0 \\ \alpha_2 \\ \alpha_3 \end{pmatrix} = \begin{pmatrix} 0 \\ \alpha_2 a_2^2 \Sigma_{22} \\ \alpha_3 a_3^2 \Sigma_{33} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \quad \forall \alpha_2, \alpha_3 \quad (2.96) $$

Par conséquent, l'égalité de (2.96) étant vérifiée pour tout $\alpha_2, \alpha_3$, on obtient deux nouvelles équations portant sur les coefficients du second tenseur de Piola-Kirchhoff.

$$ \Sigma(\mathbf{x}^0, t)_{22} = \Sigma(\mathbf{x}^0, t)_{33} = 0 \quad (2.97) $$

Conformément à la loi de comportement des matériaux de St Venant-Kirchhoff (2.77), les deux équations de (2.97) sont équivalentes à :

$$ 2\mu \mathcal{E}(\mathbf{x}^0, t)_{22} = 2\mu \mathcal{E}(\mathbf{x}^0, t)_{33} = -\lambda \operatorname{tr}(\mathcal{E}(\mathbf{x}^0, t)) \quad (2.98) $$

Après calculs, on débouche sur l'égalité suivante :

$$ \mathcal{E}(\mathbf{x}^0, t)_{22} = \mathcal{E}(\mathbf{x}^0, t)_{33} = -\nu \mathcal{E}(\mathbf{x}^0, t)_{11} \quad (2.99) $$

Cette dernière équation permet de donner l'expression définitive du tenseur de Green-Lagrange $\mathcal{E}(\mathbf{x}^0, t)$:

$$ \mathcal{E}(\mathbf{x}^0, t) = \frac{1}{2} (\gamma (j, \omega, t)^2 - 1) \begin{pmatrix} 1 & & \\ & -\nu & \\ & & -\nu \end{pmatrix} \quad (2.100) $$

Le tenseur $\mathcal{E}(\mathbf{x}^0, t)$ présente une caractéristique structurelle très intéressante. En fait, la déformation le long de l'axe $\mathbf{u}^1$ est couplée avec une déformation radiale dans le plan de la section définie par $(\mathbf{u}^2, \mathbf{u}^3)$ pondérée par le coefficient de Poisson $\nu$ qui est, en mécanique, la mesure naturelle de la striction.

Le premier terme du second tenseur des contraintes de Piola-Kirchhoff peut être obtenu grâce à la loi de comportement (2.77).

$$ \begin{align} \Sigma(\mathbf{x}^0, t)_{11} &= \lambda \operatorname{tr}(\mathcal{E}(\mathbf{x}^0, t)) + 2\mu \mathcal{E}(\mathbf{x}^0, t)_{11} = -2\mu \mathcal{E}(\mathbf{x}^0, t)_{22} + 2\mu \mathcal{E}(\mathbf{x}^0, t)_{11} \nonumber \\ &= 2\mu \left[ -(-\nu \mathcal{E}(\mathbf{x}^0, t))_{11} + \mathcal{E}(\mathbf{x}^0, t)_{11} \right] = \frac{2E}{2(1+\nu)} (1+\nu) \mathcal{E}(\mathbf{x}^0, t)_{11} \tag{2.101} \\ &= E \mathcal{E}(\mathbf{x}^0, t)_{11} \nonumber \end{align} $$
En définitive,

$$ \Sigma(\mathbf{x}^0, t) = \frac{E}{2} (\gamma(j, \omega, t)^2 - 1) \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \qquad (2.102) $$

Tout matériau de St Venant-Kirchhoff étant hyper élastique, je donne l'expression de l'énergie de déformation pour un volume unitaire lagrangien (i.e. exprimé dans la configuration de référence).

$$ \begin{aligned} e(\mathbf{x}^0, t) &= \frac{1}{2} \Sigma(\mathbf{x}^0, t) : \mathcal{E}(\mathbf{x}^0, t) \\ &= \frac{E}{2} \mathcal{E}_{11}(\mathbf{x}^0, t)^2 \qquad (2.103) \\ &= \frac{E}{8} (\gamma(j, \omega, t)^2 - 1)^2 \end{aligned} $$

L'énergie totale de déformation est obtenue par sommation.

$$ e^{\text{totale}}(t) = \int e(\mathbf{x}^0, t) d\mathbf{x}^0 = \iiint e(\mathbf{x}^0, t) dx_1^0 dx_2^0 dx_3^0 \qquad (2.104) $$

A l'aide d'un changement de coordonnées approprié relatif à l'état de repos de la courbe spline matérielle,

$$ \begin{gathered} x_1^0 = p_1^0(j, \omega_1), \quad x_2^0 = \omega_2, \quad x_3^0 = \omega_3 \\ dx_1^0 = \partial_{\omega_1} l^0(j, \omega_1) d\omega_1, \quad dx_2^0 = d\omega_2, \quad dx_3^0 = d\omega_3 \end{gathered} \qquad (2.105) $$

On débouche sur l'expression suivante :

$$ \begin{aligned} e^{\text{totale}}(t) &= \sum_{j=1}^{ns} \iiint e(j, \omega_1, t) \partial_{\omega_1} l^0(j, \omega_1) d\omega_1 d\omega_2 d\omega_3 \\ &= \sum_{j=1}^{ns} \iint d\omega_2 d\omega_3 \int e(j, \omega_1, t) \partial_{\omega_1} l^0(j, \omega_1) d\omega_1 \end{aligned} \qquad (2.106) $$

On note s la surface de la section infinitésimale supposée constante

$$ e^{\text{totale}}(t) = \frac{E s}{8} \sum_{j=1}^{ns} \int (\gamma(j, \omega, t)^2 - 1)^2 \partial_{\omega} l^0(j, \omega) d\omega \qquad (2.107) $$
D'une manière différente,

$$
\begin{align}
e^{\text{totale}}(t) &= \frac{E s}{2} \{c_1(t) - c_2(t) + c_3\} \quad \text{avec} \nonumber \\
c_1(t) &= \frac{1}{4} \sum_{j=1}^{ns} \int \frac{\partial_{\omega} l(j, \omega, t)^4}{\partial_{\omega} l^0(j, \omega)^3} d\omega, \quad c_2(t) = \frac{1}{2} \sum_{j=1}^{ns} \int \frac{\partial_{\omega} l(j, \omega, t)^2}{\partial_{\omega} l^0(j, \omega)} d\omega \tag{2.108} \\
\text{et } c_3 &= \frac{1}{4} \sum_{j=1}^{ns} \int \partial_{\omega} l^0(j, \omega) d\omega \nonumber
\end{align}
$$

5.1 Aspects informatiques de l'énergie continue de déformation

Conformément au formalisme lagrangien, les efforts interviennent dans le moteur d'animation
dynamique par le biais de leurs puissances virtuelles dans les mouvements virtuels induits par
chaque coordonnée généralisée ηᵢ. De plus, si les efforts dérivent d'un potentiel, leurs
puissances virtuelles peuvent être directement déduites par dérivation de l'énergie potentielle :

$$
Q_i(t) = - \frac{\partial e^{\text{totale}}}{\partial \eta_i}(t) \qquad (2.109)
$$

Les coordonnées généralisées ηᵢ+(k-1)n correspondent, pour une courbe spline materielle, aux coordonnées pᵢⱼ de ses points de contrôle. En s'appuyant sur la définition de l'énergie totale de déformation donnée en (2.108), on déduit l'expression des coordonnées Qᵢ+(k-1)n(t) du vecteur des forces généralisées correspondant aux contraintes internes.

$$
Q_{i+(k-1)n}(t) = - \frac{\partial e^{\text{totale}}}{\partial p_k^i}(t) = \frac{E s}{2} \sum_{m=1}^{n} p_k^m(t) \left[ \mathcal{B}_{im} - \sum_{p,q=1}^{n} \mathcal{B}_{impq} \mathbf{p}^p(t) \cdot \mathbf{p}^q(t) \right] \quad (2.110)
$$

ou les termes $\mathcal{B}_{impq}$ et $\mathcal{B}_{im}$ correspondent respectivement à

$$
\mathcal{B}_{impq} = \sum_{j=1}^{ns} \int \frac{b_i^j(\omega) b_m^j(\omega) b_p^j(\omega) b_q^j(\omega)}{\partial_{\omega} l^0(j, \omega)^3} d\omega \quad \text{et} \quad \mathcal{B}_{im} = \sum_{j=1}^{ns} \int \frac{b_i^j(\omega) b_m^j(\omega)}{\partial_{\omega} l^0(j, \omega)} d\omega \quad (2.111)
$$

Même si l'on considère des fonctions de mélange polynomiales, les scalaires $\mathcal{B}_{im}$ et $\mathcal{B}_{impq}$ sont le résultat de l'intégration d'une fonction rationnelle qui ne peut être a priori calculée de façon exacte. Mais puisque ces termes intégraux sont constants, donc évalués une seule fois durant l'étape d'initialisation, je suis en mesure d'utiliser une méthode d'intégration numérique d'ordre élevé pour en déduire une valeur approchée de précision souhaitée.

En observant l'équation (2.110), le calcul des puissances virtuelles des efforts internes d'élasticité révèle une complexité d'ordre o(n⁴), n correspondant au nombre de points de contrôle. Par contre, dans un cadre pratique, nombre de modèles splines vérifient la propriété de localité. Cette propriété garantit que la définition d'un segment spline ne dépend que d'un petit nombre noté δ de points de contrôle consécutifs.
On a alors,

$$b_i^j(\omega) = 0 \quad \text{pour } i \notin \mathfrak{T}^j, 1 \le j \le ns, \omega \in [0,1] \qquad (2.112)$$

$\mathfrak{T}^j$ est l'ensemble des indices des $\delta$ points de contrôle intervenant dans la définition du segment spline $j$. Il apparaît alors plus judicieux de permuter les sommations afin d'obtenir une expression alternative des puissances virtuelles des efforts internes d'élasticité.

$$Q_{i+(k-1)n}(t) = \frac{Es}{2} \sum_{j \in \mathfrak{T}^i} \sum_{m \in \mathfrak{T}^j} p_k^m(t) \left[ \mathcal{B}_{im}^j - \sum_{p,q \in \mathfrak{T}^j} \mathcal{B}_{impq}^j \mathbf{p}^p(t) \cdot \mathbf{p}^q(t) \right] \quad (2.113)$$
$$\text{où } \mathfrak{T}^i = \{ j \in [1,ns] : i \in \mathfrak{T}^j \}$$

$$\mathcal{B}_{impq}^j = \int \frac{b_i^{j'}(\omega)b_m^{j'}(\omega)b_p^{j'}(\omega)b_q^{j'}(\omega)}{\partial_{\omega} l^0(j,\omega)^3} d\omega, \mathcal{B}_{im}^j = \int \frac{b_i^{j'}(\omega)b_m^{j'}(\omega)}{\partial_{\omega} l^0(j,\omega)} d\omega \quad (2.114)$$

Le calcul des termes de l'équation (2.113) est d'une complexité d'ordre $o(ns \cdot \delta^4)$ si on le réalise à l'aide d'une boucle sur les $ns$ segments $j$ qui renseigne les $3\delta$ coordonnées $Q_{i+(k-1)n}(t)$ (pour $i \in \mathfrak{T}^j$) en supposant que les $\delta^4$ termes $\mathcal{B}_{impq}^j$ et les $\delta^2$ termes $\mathcal{B}_{im}^j$ (pour $i,m,p,q \in \mathfrak{T}^j$) ont été pré-calculés pour chaque segment durant l'étape d'initialisation. Comme le nombre $\delta$ est une constante du modèle spline choisi, la complexité effective est de l'ordre de $o(ns) = o(n)$.

## 5.2 Loi de ressort équivalente

Afin de comparer l'énergie continue de déformation que je propose avec une succession de ressorts répartis le long de la courbe spline matérielle, je dois être en mesure d'exprimer la loi de ressort dont la réponse est équivalente à celle fournie par l'énergie de déformation pour une elongation donnée. Je rappelle l'expression alternative du tenseur des contraintes de Cauchy $\sigma(\mathbf{x}(t))$ à partir du second tenseur des pseudo contraintes de Piola-Lagrange $\Sigma(\mathbf{x}^0, t)$.

$$\sigma(\mathbf{x}(t)) = \det(\mathcal{F}(\mathbf{x}^0, t))^{-1} \mathcal{F}(\mathbf{x}^0, t) : \Sigma(\mathbf{x}^0, t) : \mathcal{F}(\mathbf{x}^0, t)^T \quad (2.115)$$

Si $\mathbf{n}$ est le vecteur normal à un élément de surface infinitésimal $dS$, la force $\mathbf{F}$ résultante d'une contrainte appliquée à $dS$ est donnée par la relation $\mathbf{F} = dS \sigma \mathbf{n}$. J'applique cette formule dans le cas particulier d'une courbe spline matérielle dont l'aire de la section instantanée est égale à $s' = a_2 a_3 s$ de vecteur normal $\mathbf{n} = \pm \mathbf{u}^1$. Puisque je ne m'intéresse qu'au cas des efforts internes d'élasticité, j'espère trouver une relation entre réponse et vecteur normal de la forme :

$$\lambda \mathbf{u}^1 = \pm s' \sigma \mathbf{u}^1 \qquad (2.116)$$

où $\lambda$ est l'intensité algébrique de la réponse de la contrainte des efforts internes. Des équations (2.88) (2.100) (2.115) et (2.116) on déduit une expression de $\lambda$ qui répond à cette attente, sous la forme d'un polynôme d'ordre 3, fonction du facteur de dilatation relative $\gamma$ :
$$ \lambda = \pm a_2 a_3 s \mathbf{u}^1 \cdot \boldsymbol{\sigma} \mathbf{u}^1 = \pm \frac{E s}{2} \gamma (\gamma^2 - 1) \quad (2.117) $$

Pour illustrer les développements théoriques aboutissant à la formule (2.117), je propose une expérience statique consistant en une mesure de réaction à la traction. Afin de demeurer statique, cette expérience doit garantir que les efforts internes sont homogènes au sein du système matériel. Cette hypothèse est réalisée en choisissant une gravité nulle et un état de départ homothétique à l'état de repos. De cette manière, les points de contrôle de l'état initial sont déduits de l'état repos par la relation suivante :

$$ \mathbf{p}^i(0) = \alpha \mathbf{p}^{0i} \quad (2.118) $$

Cette définition de la position des points de contrôle implique :

$$ \mathbf{p}(j, \omega, 0) = \alpha \mathbf{p}^0(j, \omega) \Rightarrow \gamma(j, \omega, 0) = \alpha \quad 1 \le j \le ns, \omega \in [0,1] \quad (2.119) $$

Afin de m'assurer que la courbe spline matérielle demeure fixe dans son état étiré, j'impose deux contraintes de position à ses extrémités. De cette manière, les efforts visant à garantir la réalisation des contraintes imposées aux faces extrêmes de la courbe doivent exactement s'opposer aux contraintes internes appliquées aux mêmes faces. Ces efforts de maintien sont, dans les équations de Lagrange, représentés par deux multiplicateurs de Lagrange notés $\lambda_0, \lambda_1$. Les valeurs respectives $\lambda_0$ et $\lambda_1$ correspondent exactement aux intensités des efforts de maintien de chaque contrainte (2.120), fournissant un moyen de mesurer les contraintes internes (opposées aux efforts de maintien) pour une elongation donnée.

$$ \lambda_0 = \pm \frac{E s}{2} \gamma (\gamma^2 - 1), \quad \lambda_1 = \mp \frac{E s}{2} \gamma (\gamma^2 - 1) \quad (2.120) $$

Les résultats de cette expérience pour différentes valeurs du facteur d'élongation relative $\gamma - 1$ sont rassemblés et comparés avec la courbe de réponse théorique (2.117) dans le tableau 1.

<figure><img src="image_64.png" /><figcaption>Tableau 1. Résultats de l'expérience de réponse à la traction.</figcaption></figure>
Les précédents développements démontrent clairement que tout segment spline rectiligne doté de cette énergie de déformation se comporte, lors d'une elongation homogène, comme un ressort de loi :

$$f = \frac{1}{2} E s (l/l_0) (l/l_0)^2 - 1) \quad (2.121)$$

Par conséquent, si l'on considère qu'un segment spline demeurent, entre deux particules, rectiligne avec une elongation homogène, les contraintes internes issues de l'énergie de déformation peuvent être approchées en attachant un ressort entre ces deux particules dont la loi correspond à l'équation (2.121). Cette hypothèse me permet d'envisager de comparer le mouvement d'une courbe spline avec énergie continue de déformation avec celui d'une courbe munie de ressorts répartis sur sa longueur

## 5.3 Energie de déformation vs. ressorts

Les tests de comparaison mis en œuvre sont d'une nature relativement simple. En effet, ils consistent en des simulations physiques de courbes splines matérielles dont les seules différences résident dans la nature des efforts internes d'élasticité. Une des cordes utilise une énergie continue de déformation pour modéliser les contraintes internes tandis que les autres sont dotées de ressorts de loi compatible dont le nombre par segments varie selon la courbe spline matérielle.

Le tableau 2 regroupe les résultats de deux expériences illustrées par la Figure 52 et la Figure 53. Le terme « Erreur » correspond à la moyenne par rapport au temps des racines carrés des moyennes des carrés des distances séparant les particules du courbe étudiée et celles de la courbe de référence dotée d'une énergie de déformation. Le « Temps » correspond au temps machine total requis pour calculer une simulation de 5 secondes sur un processeur AMD Athlon 1.2 GHz.

Ces résultats numériques ainsi que les images issues des différentes animations dynamiques révèlent que l'usage de trop peu de ressorts par segment (moins d'une douzaine) entraîne des erreurs inacceptables mais surtout que le modèle s'appuyant sur l'énergie de déformation est d'un coût équivalent à un modèle à base de ressorts dont le nombre garantit seulement une précision moyenne.

<table><caption>Tableau 2 Résultats numériques des deux expériences réalisées avec des cordes.</caption><thead><tr><th rowspan="2">Ressorts par segment</th><th colspan="2">Corde fixée à une extrémité</th><th colspan="2">Corde fixée aux deux extrémités</th></tr><tr><th>Erreur</th><th>Temps</th><th>Erreur</th><th>Temps</th></tr></thead><tbody><tr><td>2</td><td>0.2721 m</td><td>2.120s</td><td>0.0976 m</td><td>4.110s</td></tr><tr><td>5</td><td>0.0459 m</td><td>3.300s</td><td>0.0814 m</td><td>5.310s</td></tr><tr><td>10</td><td>0.0234 m</td><td>5.290s</td><td>0.0825 m</td><td>7.280s</td></tr><tr><td>15</td><td>0.0133 m</td><td>7.260s</td><td>0.0825 m</td><td>9.260s</td></tr><tr><td>20</td><td>0.0107 m</td><td>9.250s</td><td>0.0825 m</td><td>11.240s</td></tr></tbody><tfoot><tr><td>Elasticité continue</td><td>-</td><td>6.230s</td><td>-</td><td>8.250s</td></tr></tfoot></table>
<figure><img src="image_65.png" /><figcaption>Figure 52 Images issues d'une simulation de cordes suspendues par une de leur extrémité. Les numéros renseignent sur le nombre de ressorts pour chacun des 4 segments, le symbole $\infty$ désigne la corde d'élasticité continue.</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>Figure 53 Images issues d'une simulation d'une corde suspendue par ses deux extrémités, la courbe transparente désigne la courbe d'élasticité continue. A gauche, 2 ressorts par segments. A droite, 10 ressorts par segment.</figcaption></figure>
# 6 Conclusion

La mise en œuvre d'un modèle continu d'élasticité interne a surtout été motivé par cette volonté d'aboutir à la réalisation d'un modèle continu unifié pour représenter des corps déformables. Cette nouvelle étape, sans être l'avènement de ce but, constitue tout de même un grand pas vers un modèle purement continu.

D'autre part, le choix de la théorie de l'élasticité en grands déplacements fournit un cadre théorique pour la description de la déformation qui soit :

* Compatible avec les modèles mécaniques utilisés dans l'industrie. Les paramètres de l'énergie de déformation correspondent aux constantes caractéristiques du matériau issues de mesures expérimentales.

* Non restreint aux petits déplacements, garantissant la robustesse du modèle pour représenter des mouvements de grande amplitude.

Je tiens néanmoins à modérer mon propos. En effet, le travail présenté dans ce chapitre, qui s'inscrit dans une démarche à long terme, constitue mes premiers pas dans la représentation des déformations. Bien que limités aux cas des courbes, les résultats obtenus sont satisfaisants à plusieurs titres. Ce modèle se révèle en fait plus proche de la nature du milieu qu'un modèle à base de ressorts sans pour autant engendrer un surcoût de calcul important.

Le développement théorique suffisamment générique peut naturellement s'étendre au cas des surfaces et des volumes splines. Mais l'accroissement de la complexité combinatoire induit par le passage aux dimensions supérieures m'encourage à envisager un modélisation modulaire des lois de comportement me permettant de choisir l'ordre maximum des termes intervenant dans l'expression de la loi. Cette modularité devrait permettre d'établir un compromis entre richesse de représentation et performance.
# REDUCTION DE PARAMETRAGE : APPLICATION AUX TEXTILES TRICOTES

## 1 Introduction

Par expérience, on constate que l'agencement complexe des mailles joue un rôle important dans le comportement dynamique global d'une pièce d'étoffe tricotée. Cette réflexion a amenée l'équipe d'animation dynamique du LERI à développer un modèle de textiles tricotés [Nourrit 1999] à l'échelle mésoscopique (une échelle intermédiaire entre l'échelle microscopique et macroscopique, décrivant le phénomène à l'échelle du fil), contrairement aux modèles dédiés aux textiles tissés dont la structure simple autorise une modélisation à l'échelle macroscopique. Le choix d'un modèle mécanique à l'échelle mésoscopique pour l'étude du mouvement d'un tricot me semble d'une part légitime mais surtout nécessaire. En effet ce modèle, respectueux des propriétés intrinsèques du fil (répartition massique et propriétés élastiques), est par ailleurs seul capable de rendre compte de la complexité structurelle des liages. Ce travail coïncide avec celui entrepris par Meißner et Eberhardt à la même époque [Meißner 1998]. En effet, ils proposent aussi un modèle de tricots à l'échelle mésoscopique en assimilant un fil à une ligne polygonale joignant ce qu'ils appellent les points de déformation qui ne sont rien d'autre que les positions échantillonnées d'application des efforts. Tout comme dans le modèle initial, les efforts d'élasticité interne sont modélisés par des ressorts disposés entre les points de déformation. Ce modèle est utilisé dans un cadre statique pour observer l'état de repos d'une structure tricoté. Cet état de repos est atteint par minimisation d'une énergie modélisant d'une part l'élasticité des fils mais aussi les contraintes visant à préserver la structure tricoté. L'année suivante, Eberhardt et Weber [Eberhardt 1999] envisagent un modèle physique pour l'animation dynamique de tricots en proposant une extension de leur modèle à base de particules, utilisé avec succès pour l'animation de textiles tissés. Ce modèle plus simple permet d'envisager la production d'animations de tricots en des temps raisonnables. Mais reste le problème de l'adéquation du modèle avec la nature fortement non linéaire du comportement des tricots : ils ont alors recours à des courbes de traction-allongement réalisées sur des échantillons de tricots. Cette approche s'éloigne d'une optique de prototypage virtuel dans la mesure où il faut concevoir les échantillons que l'on souhaite simuler. Une autre lacune du modèle proposé réside
dans la réponse au cisaillement et à la flexion. En effet, ne disposant pas de courbes de réponse de cette nature puisque jugées non significatives pour des tricots, ils utilisent des courbes expérimentales issues de tests sur des textiles tissés. Ce choix arbitraire renforce ma conviction d'aborder l'animation dynamique de textiles tricotés par le biais d'un modèle mécanique à l'échelle mésoscopique.

Malheureusement la précision de ce modèle aboutit à des systèmes d'équations de grande taille dont la résolution requiert des temps de calcul prohibitifs, rendant cette approche inadaptée à l'étude du comportement de pièces d'étoffe tricotée de grandes dimensions. Il semble donc primordial d'y apporter des modifications afin de permettre la production d'animations de tricots en des temps de calcul raisonnables. La solution que je propose s'articule autour d'une méthode robuste et générique, devant déboucher à terme sur un modèle à plusieurs niveaux de détail mécanique pour l'animation, permettant d'adapter le nombre de données à gérer lors de la génération du mouvement. Cette idée de niveau de détail mécanique a déjà été abordée dans le cadre des modèles « masses ressorts ». Sa mise en œuvre consiste en un raffinement ou un appauvrissement du réseau des points matériels [Palazzi 1994] [Hutchinson 1996]. Cette méthodologie exploite l'aspect discret des modèles « masses ressorts », elle n'est donc pas envisageable dans mon cadre de travail où j'attache une grande importance au caractère continu de l'approche dynamique. Je présente donc une nouvelle méthode conciliant mon souci de manipuler un modèle continu et la garantie d'une réduction importante de la taille du système à résoudre. Pour cela, je m'appuie sur la méthode de changement de paramétrage qui, à partir d'un nouveau jeu de coordonnées généralisées réduit et d'une fonction de changement de paramétrage, permet de reconstruire le paramétrage initial. Yannick Remion illustre les nombreux intérêts de cette méthode en l'appliquant à la gestion des contraintes par réduction cinématique [Rémion 2000]. Pour ma part, j'exploite cette démarche afin de construire un nouveau paramétrage du tricot dont les configurations accessibles demeurent suffisamment diversifiées tout en garantissant que la taille du nouveau système d'équations soit bien moindre. Je commence par rappeler les étapes importantes de la méthode de changement de paramétrage. Notamment je détaille comment, à partir d'un nouveau jeu de coordonnées généralisées et de la fonction de changement de paramétrage associée, on exprime les termes fondateurs du nouveau système d'équations à partir des termes du système initial. Enfin, j'expose la mise en œuvre de cette méthode dans le cas particulier des tricots en étayant ma démarche de résultats qualitatifs et quantitatifs justifiant l'intérêt de cette approche en termes de performance et de respect du comportement global de l'étoffe. Les résultats présentés dans ce chapitre ont déjà fait l'objet de publications [Nocent 2001a] [Nocent 2001b].

## 2 Réduction de paramétrage

Le méthode de réduction de paramétrage repose sur la construction d'un nouveau jeu de coordonnées généralisées $\tilde{\eta}$ évoluant dans un espace réduit des configurations étendues $\tilde{U}$ de dimension $\tilde{n}+1 \le n+1$. On définit alors une **fonction de réduction de paramétrage** $\gamma: \tilde{U} \to U$ qui, à partir des coordonnées généralisées réduites $\tilde{\eta}$, permet de calculer le vecteur $\eta$ des coordonnées généralisées initiales. On a donc la relation suivante :

$$ \eta = \gamma(\tilde{\eta}) \text{ avec } \tilde{\eta} = (\tilde{\eta}_0 \ \tilde{\eta}_1 \ \dots \ \tilde{\eta}_{\tilde{n}})^T \text{ et } \tilde{\eta}_0 = \eta_0 = t \quad (3.1) $$
Ainsi, l'ensemble des configurations permises par le paramétrage réduit $\tilde{\eta}$ correspond au sous-espace $\gamma(\tilde{U})$ inclus dans l'espace des configurations étendues $\tilde{U}$. Dans la mesure où les paramétrages $\eta$ et $\tilde{\eta}$ décrivent la configuration d'un même système matériel, les termes constitutifs des équations du mouvement associées à chaque paramétrage doivent être en relation. En dérivant l'équation (3.1), on établit une première relation entre les taux de variation $\hat{\eta}$ et $\hat{\tilde{\eta}}$ de chaque paramétrage.

$$ \hat{\eta} = \nabla(\tilde{\eta})\hat{\eta} \quad \text{où} \quad \nabla(\tilde{\eta}) = \frac{\partial\gamma}{\partial\tilde{\eta}}(\tilde{\eta}) \text{ est la matrice jacobienne de } \gamma \quad (3.2) $$

En s'appuyant sur la remarque suivante,

$$ \nabla_{0j}(\tilde{\eta}) = \frac{\partial}{\partial \tilde{\eta}_j} \gamma_0(\tilde{\eta}) = \frac{\partial}{\partial \tilde{\eta}_j} \tilde{\eta}_0 = \delta_{0j} \quad (3.3) $$

la matrice jacobienne $\nabla(\tilde{\eta})$ de $\gamma$ présente des propriétés structurelles intéressantes, synthétisées par l'équation (3.4) :

$$ \nabla(\tilde{\eta}) = \begin{bmatrix} 1 & \mathbf{0} \\ \frac{\partial}{\partial t} \gamma^*(\tilde{\eta}) & \nabla^{**}(\tilde{\eta}) \end{bmatrix} \quad (3.4) $$

Avec l'équation (3.5), je rappelle l'expression du système différentiel équivalent aux équations de Lagrange. Pour le mouvement réel $\tilde{\varepsilon}(t)$ exprimé dans le paramétrage réduit, la matrice $\tilde{M}(\tilde{\varepsilon}(t))$ correspond à la matrice des masses généralisées et le vecteur $\tilde{E}(\tilde{\varepsilon}(t), \dot{\tilde{\varepsilon}}(t))$ représente l'ensemble des efforts appliqués au système comme une différence des forces généralisées $\tilde{Q}(\tilde{\varepsilon}(t), \dot{\tilde{\varepsilon}}(t))$ et des termes inertiels $\tilde{\mathrm{In}}(\tilde{\varepsilon}(t), \dot{\tilde{\varepsilon}}(t))$. Je ne traite pas ici le problème du changement d'expression des contraintes, ce travail ayant déjà été abordé en détail par Yannick Remion pour les besoins de son habilitation à diriger des recherches [Remion 2000].

$$ \tilde{M}(\tilde{\varepsilon}(t)) \ddot{\tilde{\varepsilon}}(t) = \tilde{E}(\tilde{\varepsilon}(t), \dot{\tilde{\varepsilon}}(t)) \quad (3.5) $$

Dans la pratique, la taille du système (3.5) est nettement plus petit que la taille du système différentiel correspondant au paramétrage initial, garantissant ainsi le gain de performance escompté. Afin d'établir les similitudes entre les deux systèmes linéaires, je commence par comparer les expressions de l'énergie cinétique propres aux deux paramétrages.

$$ \begin{align} \kappa(\tilde{\eta}, \hat{\eta}) &= \kappa(\eta = \gamma(\tilde{\eta}), \hat{\eta} = \nabla(\tilde{\eta})\hat{\eta}) \nonumber \\ &= \frac{1}{2} \hat{\eta}^T \mathbf{M}(\eta) \hat{\eta} \tag{3.6} \\ &= \frac{1}{2} \hat{\eta}^T (\nabla(\tilde{\eta})^T \mathbf{M}(\gamma(\tilde{\eta})) \nabla(\tilde{\eta})) \hat{\eta} \nonumber \end{align} $$
Au final, on obtient la relation suivante :

$$ \tilde{M}(\tilde{\eta}) = \nabla(\tilde{\eta})^T M(\gamma(\tilde{\eta})) \nabla(\tilde{\eta}) = \nabla(\tilde{\eta})^T M(\eta) \nabla(\tilde{\eta}) \quad (3.7) $$

Plus précisément,

$$ \tilde{M}^{**}(\tilde{\eta}) = \nabla^{*}(\tilde{\eta})^{T} M(\eta) \nabla^{*}(\tilde{\eta}) = [\mathbf{0} \quad \nabla^{**}(\tilde{\eta})^{T}] M(\eta) \begin{bmatrix} \mathbf{0} \\ \nabla^{**}(\tilde{\eta}) \end{bmatrix} \quad (3.8) $$

soit en définitive,

$$ \tilde{M}^{**}(\tilde{\eta}) = \nabla^{**}(\tilde{\eta})^{T} M^{**}(\eta) \nabla^{**}(\tilde{\eta}) \quad (3.9) $$

Je m'intéresse maintenant au terme des efforts généralisés $\tilde{Q}(\tilde{\eta}, \hat{\tilde{\eta}})$ correspondant au vecteur des puissances virtuelles des efforts donnés dans les $\tilde{n}$ mouvements figés d'entraînement du paramétrage réduit. On cherche en fait à exprimer le vecteur $\tilde{Q}(\tilde{\eta}, \hat{\tilde{\eta}})$ à partir du terme $Q(\eta, \hat{\eta})$, vecteur des puissances virtuelles des efforts donnés dans les $n$ mouvements figés d'entraînement du paramétrage initial.

La définition du terme $Q(\eta, \hat{\eta})$ dans les équations de Lagrange vient de l'expression des puissances virtuelles des efforts donnés appliqués au système dans les différents mouvements virtuels compatibles avec le paramétrage initial $\eta$ et coïncidant en l'instant étudié avec le mouvement « réel » étudié. La puissance $P^*(\eta, \hat{\eta}, \tilde{\eta}^*)$ des efforts donnés appliqués au système matériel dans un mouvement virtuel compatible avec le paramétrage $\eta$ et coïncidant en l'instant étudié avec le mouvement « réel » étudié s'écrit :

$$ P^*(\eta, \hat{\eta}, \tilde{\eta}^*) = \sum_{i=1}^{n} Q_i(\eta, \hat{\eta}) \hat{\eta}_i^* = \hat{\eta}^* \cdot Q(\eta, \hat{\eta}) \quad (3.10) $$

De façon équivalente, la puissance $\tilde{P}^*(\tilde{\eta}, \hat{\tilde{\eta}}, \tilde{\eta}^*)$ de ces efforts dans un mouvement virtuel compatible avec le paramétrage $\tilde{\eta}$ vaut :

$$ \tilde{P}^*(\tilde{\eta}, \hat{\tilde{\eta}}, \tilde{\eta}^*) = \sum_{i=1}^{\tilde{n}} \tilde{Q}_i(\tilde{\eta}, \hat{\tilde{\eta}}) \hat{\eta}_i^* = \hat{\tilde{\eta}}^* \cdot \tilde{Q}(\tilde{\eta}, \hat{\tilde{\eta}}) \quad (3.11) $$

Par construction, un mouvement virtuel compatible avec le paramétrage réduit $\tilde{\eta}$ est forcément compatible avec le paramétrage initial $\eta$. Ainsi, les mouvements virtuels utilisés dans (3.11) sont aussi valides pour (3.10) et l'expression de la puissance virtuelle des efforts donnés dans ces mouvements peut être indifféremment calculée d'après (3.10) ou (3.11). On aboutit donc à l'égalité suivante :

$$ \begin{align} \hat{\eta}^* \cdot \tilde{Q}(\tilde{\eta}, \hat{\tilde{\eta}}) &= (\nabla^{**}(\tilde{\eta}) \hat{\eta}^*) \cdot Q(\gamma(\tilde{\eta}), \nabla(\tilde{\eta}) \hat{\eta}) \nonumber \\ &= (\nabla^{**}(\tilde{\eta}) \hat{\eta}^*)^T Q(\gamma(\tilde{\eta}), \nabla(\tilde{\eta}) \hat{\eta}) \tag{3.12} \\ &= \hat{\eta}^* \cdot (\nabla^{**}(\tilde{\eta})^T Q(\gamma(\tilde{\eta}), \nabla(\tilde{\eta}) \hat{\eta})) \nonumber \end{align} $$
L'égalité (3.12) est valide pour tout mouvement virtuel compatible avec le paramétrage réduit $\tilde{\eta}$. Ainsi, en considérant les mouvements figés d'entraînement, on arrive à la conclusion suivante :

$$ \tilde{Q}(\tilde{\eta}, \hat{\eta}) = \nabla^{**}(\tilde{\eta})^T Q(\eta, \hat{\eta}) \quad (3.13) $$

Il est important de mentionner que les efforts considérés sont identiques dans les deux paramétrages. Mais comme les évolutions cinématiques permises par les deux paramétrages sont différentes, les effets de ces efforts sur le mouvement du système, « mesurés » par leurs puissances virtuelles dans des jeux de mouvements virtuels différents, ne sont pas identiques selon le paramétrage. Par contre si un mouvement réel $\varepsilon(t)$ exprimé dans le paramétrage initial reste toujours compatible avec le paramétrage réduit, son analyse à l'aide de ce paramétrage réduit donne exactement le même mouvement.

Enfin, je termine l'exposé de la méthode de réduction de paramétrage en donnant succinctement l'expression des termes inertiels (3.14) dans le paramétrage réduit $\tilde{\eta}$. Je m'attarde peu sur l'expression de ces termes dans la mesure où, comme je le montrerai dans la section suivante, le paramétrage réduit que j'utilise pour les tricots génère un nouveau système différentiel à matrice des masses généralisées constante, dépourvu par conséquent de termes inertiels. J' invite le lecteur désireux d'obtenir plus de détails concernant les termes inertiels à lire l'annexe II de l'habilitation à diriger des recherches de Yannick Remion.

$$ \begin{gathered} \tilde{\mathrm{In}}(\tilde{\eta}, \hat{\eta}) = \nabla^{**}(\tilde{\eta})^T (\mathrm{In}(\eta, \hat{\eta}) + \mathbf{M}^{**}(\eta)\Psi(\tilde{\eta}, \hat{\eta})) \\ \Psi_k(\tilde{\eta}, \hat{\eta}) = \hat{\eta}^T \mathbf{H}^k(\tilde{\eta}) \hat{\eta} \end{gathered} \quad (3.14) $$

$$ H_{ij}^k(\tilde{\eta}) = \frac{\partial}{\partial \tilde{\eta}_{ij}} \gamma_k(\tilde{\eta}) \quad \text{matrice hessienne} \quad (3.15) $$

## 3 Application à l'animation des textiles tricotés

Conformément au modèle mésoscopique proposé par Jean Michel Nourrit [Nourrit 1999], le paramétrage d'un morceau d'étoffe est donné par les positions des points de contrôle de chaque fil qui le constitue. Mais ce paramétrage, au nombre de degrés de liberté important, est de nature incomplète puisqu'il ne rend pas compte de la nature du liage. La structure des boucles est alors modélisée à l'aide de contraintes de contact entre fils. Ces contraintes, qui implicitement réduisent considérablement le nombre de degrés de liberté de la pièce d'étoffe, engendrent de nouvelles inconnues correspondant aux efforts de liaison garantissant leur réalisation. Ceci a pour conséquence d'accroître la taille déjà très grande du système d'équations à résoudre.

Ainsi, une réduction du paramétrage est tout à fait justifiée puisqu'un grand nombre des inconnues intervenant dans la résolution ne correspondent pas à proprement parler à de nouveaux degrés de liberté mais plutôt à des efforts venant réduire l'ensemble des configurations possibles. Ainsi, au prix d'une hypothèse simplificatrice, j'ai mis en place le
mécanisme de réduction de paramétrage en vue de se rapprocher du nombre véritable de degrés de liberté de l'étoffe tout en s'efforçant de conserver un comportement macroscopique riche. Bien qu'il ne s'agisse pas exactement d'une réduction cinématique (i.e. utilisation des équations de contraintes pour exprimer certains degrés de liberté en fonction d'autres), cette approche a deux avantages majeurs. Elle permet tout d'abord de diminuer le nombre d'inconnues correspondant au système matériel étudié. Mais surtout, elle autorise une prise en compte implicite des nombreuses contraintes de contact entre fils grâce à la réduction de paramétrage choisie. Je développe maintenant en détail ma solution au problème de réduction de paramétrage pour les textiles tricotés.

## 3.1 Mise en œuvre

### 3.1.1 Volume paramétrique englobant

Pour la réduction du paramétrage, j'ai choisi d'utiliser un volume paramétrique englobant l'ensemble des mailles de l'étoffe. Puisque les variations d'épaisseur d'un tricot sont relativement faibles, il apparaît plus pertinent de privilégier, à la place d'un volume spline dans les trois directions, un volume dont la variation paramétrique en épaisseur soit linéaire. Ce choix permet de réduire encore le nombre de degrés de liberté du tricot. En fait, ce volume noté V est obtenu par interpolation linéaire entre deux surfaces paramétriques $S_1$ et $S_2$ « parallèles » à la surface moyenne du liage (Figure 54). Ces deux surfaces $S_1$ et $S_2$ sont des surfaces splines bicubiques séparables définies à partir d'un même nombre $n = (n_1, n_2)$ de points de contrôle notés respectivement $\mathbf{p}^{S_1,i}$ et $\mathbf{p}^{S_2,i}$, réparties en un même nombre de patches $\mathbf{n}s = (ns_1, ns_2)$.

<figure><img src="image_67.png" /><figcaption>Figure 54 Volume paramétrique V composé de (2,2,1) patches, défini à partir de (5,5,2) points de contrôle.</figcaption></figure>

Pour garantir la lisibilité du document, je choisis de reprendre les notations utilisées pour les variétés splines. Les points de contrôle du volume V sont donc indexés à l'aide d'un triplet $i = (i_1, i_2, i_3)$ de la façon suivante :

$$ \forall \mathbf{i} = (i_1, i_2, i_3) \in \mathbb{N}_n^2 \times \{1,2\} : \mathbf{p}^{\mathbf{V},i} = \begin{cases} \mathbf{p}^{S_1, (i_1, i_2)} & \text{si } i_3 = 1 \\ \mathbf{p}^{S_2, (i_1, i_2)} & \text{sinon} \end{cases} \qquad (3.16) $$
De cette manière, le volume englobant $V$ doit être considéré comme une variété spline de dimension 3 composée de $(n_1, n_2, 2)$ (que l'on notera désormais **n**) points de contrôle et de $(ns_1, ns_2, 1)$ (dorénavant noté **ns**) sous-variétés splines. La position temporelle $\mathbf{p}(t, \mathbf{j}, \mathbf{\omega})$ d'un point du volume $V$ de coordonnées paramétriques $\mathbf{\omega} = (\omega_1, \omega_2, \omega_3)$ à l'intérieur de la sous variété d'indice $\mathbf{j} = (j_1, j_2, 1)$ délimitée par les patches de $S_1$ et $S_2$ d'indice $(j_1, j_2)$ est donnée par l'expression suivante :

$$ \mathbf{p}: \left\{ \begin{array}{l} \mathbb{R} \times \mathbb{N}_{\mathrm{ns}}^3 \times [0,1]^3 \to \mathbb{R}^3 \\ (t, \mathbf{j}, \mathbf{\omega}) \mapsto \mathbf{p}(t, \mathbf{j}, \mathbf{\omega}) = \sum_{1 \le i \le n} b_i^j(\mathbf{\omega}) \mathbf{p}^{V,i} \end{array} \right. \qquad (3.17) $$

Où la fonction de mélange $b_i^j$ correspond au poids de la contribution des points de contrôle d'indice $\mathbf{i}$ pour la définition de la sous-variété $\mathbf{j}$, obtenue par le produit de la fonction de mélange séparable $b_{i_1}^{j_1}(\omega_1)b_{i_2}^{j_2}(\omega_2)$ des surfaces $S_1$ et $S_2$. avec une fonction d'interpolation linéaire de la variable $\omega_3$.

$$ b_i^j(\mathbf{\omega}) = b_{i_1}^{j_1}(\omega_1) b_{i_2}^{j_2}(\omega_2) b_{i_3}^{j_3}(\omega_3) \quad \text{avec} \quad b_{i_3}(\omega_3) = \begin{cases} (1-\omega_3) & \text{si } i_3 = 1 \\ \omega_3 & \text{sinon} \end{cases} \qquad (3.18) $$

Les points de contrôle $(\mathbf{p}^{V,i})_{1 \le i \le n}$ du volume englobant $V$ correspondent aux nouveaux degrés du liberté du système matériel. On introduit maintenant la bijection $\alpha^V$, fonction de transcodage permettant d'établir le lien entre l'indexation multi-indicielle des points de contrôle et l'indexation entière du vecteur des coordonnées généralisées.

$$ \alpha^V: \left\{ \begin{array}{l} \mathbb{N}_n^3 \to \{1, \dots, n\} \quad \text{avec} \quad n = \mathrm{Card}(\mathbf{n}) = 2n_1n_2 \\ (i_1, i_2, i_3) \mapsto \alpha^V(i_1, i_2, i_3) = (i_3 - 1)n_1n_2 + \alpha^{S_1}(i_1, i_2) \end{array} \right. \qquad (3.19) $$

D'une manière classique, les coordonnées généralisées $\eta^V$ du volume englobant $V$ sont obtenues grâce à la relation ci-dessous :

$$ \forall \mathbf{1} \le \mathbf{i} \le \mathbf{n}: \eta_{\alpha^V(i)+(k-1)n}^V = p_k^{V,i} \quad 1 \le k \le 3 \qquad (3.20) $$

Par conséquent, le vecteur des coordonnées généralisées $\eta^V$ a la structure suivante :

$$ \eta^V = \left\{ \begin{array}{c} \underbrace{p_1^{S_1,1}, \dots, p_1^{S_1,n}}_{\text{surface } S_1}, \dots, \underbrace{p_1^{S_2,1}, \dots, p_1^{S_2,n}}_{\text{surface } S_2}, \dots, \underbrace{p_2^{S_1,1}, \dots, p_2^{S_1,n}}_{\text{surface } S_1}, \dots, \underbrace{p_2^{S_2,1}, \dots, p_2^{S_2,n}}_{\text{surface } S_2}, \dots, \underbrace{p_3^{S_1,1}, \dots, p_3^{S_1,n}}_{\text{surface } S_1}, \dots, \underbrace{p_3^{S_2,1}, \dots, p_3^{S_2,n}}_{\text{surface } S_2}, \dots, \underbrace{p_3^{S_3,1}, \dots, p_3^{S_3,n}}_{\text{surface } S_3} \\ \multicolumn{1}{c}{\left( \begin{array}{c} \mathbf{1} \\ \mathbf{x} \\ \mathbf{y} \\ \mathbf{z} \end{array} \right)^T} \end{array} \right\} $$

A partir de la définition des coordonnées généralisées $\eta^V$ du volume englobant $V$ (3.20) et de la position temporelle (3.17), je suis en mesure de définir la fonction de positionnement $\Phi$ associée à cette variété spline.
$$
\Phi : \left\{
\begin{array}{l}
\mathbb{U} \times \mathbb{N}_{\mathrm{ns}}^3 \times [0.1]^3 \to \mathbb{R}^3 \\
(\boldsymbol{\eta}, \mathbf{j}, \boldsymbol{\omega}) \mapsto \Phi(\boldsymbol{\eta}, \mathbf{j}, \boldsymbol{\omega})
\end{array}
\right.
\qquad (3.21)
$$
$$
\text{où } \Phi_k(\boldsymbol{\eta}, \mathbf{j}, \boldsymbol{\omega}) = \sum_{1 \le i \le n} b_i^j(\boldsymbol{\omega}) \eta_{\alpha(i)+(k-1)n} = \sum_{1 \le i \le n} b_{\alpha^{-1}(i)}^j(\boldsymbol{\omega}) \eta_{i+(k-1)n} \quad 1 \le k \le 3
$$

3.1.2 Plongement du liage

Avant toute chose, pour alléger les notations déjà chargées, je considère des évolutions ne dépendant pas explicitement du temps. De cette manière, j'omettrai dorénavant de mentionner la coordonnée $\eta_0$ du vecteur des configurations étendues $\eta$.

Je dispose d'ores et déjà d'un jeu réduit de degrés de liberté pour un tricot, correspondant aux
points de contrôle du volume englobant V. Dans un deuxième temps, il me faut définir une
fonction de réduction de paramétrage pour déduire des degrés de liberté réduits les degrés de
liberté initiaux, correspondant aux points de contrôle de chacun des fils. A cette fin, je formule
l'hypothèse simplificatrice suivante : les points de contrôle des fils sont à coordonnées
paramétriques constantes au sein du volume englobant V. Cette hypothèse est à rapprocher de
celle effectuée dans les **Free Form Deformations** [Sederberg 1986] dans la mesure où les fils
sont considérés comme figés dans un volume déformable. Mais, conformément à la démarche
de changement de paramétrage, j'ajoute une dimension dynamique à l'approche de Sederberg
qui constitue l'innovation de cette méthode. En effet, les déformations appliquées au volume V
proviennent de l'application des efforts sur les fils du tricot et non d'une simple transformation
géométrique.

Cette fonction de réduction de paramétrage que je qualifie de **fonction de plongement du liage** permet, à partir d'une configuration $\eta^V$ du volume englobant V, de calculer les positions des points de contrôle de chaque fil constitutif du tricot.

Soit un tricot composé de *r* rangées de mailles et par conséquent de *r* fils, on note $(\mathbf{p}^{f,i})_{1\le i \le nf}$ le vecteur constitué des *nf* points de contrôle du fil $f$ ($1 \le f \le r$).

Les coordonnées généralisées $\eta^f = (\eta_i^f)_{1\le i\le 3nf}$ du fil $f$ s'expriment de la manière suivante :

$$
\eta^f = \left( \underbrace{p_1^{f,1}, \dots, p_1^{f,nf}}_{x}, \underbrace{p_2^{f,1}, \dots, p_2^{f,nf}}_{y}, \dots, \underbrace{p_3^{f,1}, \dots, p_3^{f,nf}}_{z} \right)^T = \begin{pmatrix} (p_1^{f,i})_{1\le i \le nf} \\ (p_2^{f,i})_{1\le i \le nf} \\ (p_3^{f,i})_{1\le i \le nf} \end{pmatrix} \quad (3.22)
$$

En s'appuyant sur ces notations, on donne la définition de la fonction $\gamma^f : \mathbb{R}^{3nV} \rightarrow \mathbb{R}^{3nf}$ qui
permet de déterminer la position des points de contrôle du fil $f$ :

$$
\begin{align*}
\eta_{k+(l-1)n}^f &= (\gamma^f(\boldsymbol{\eta}^V))_{k+(l-1)n} \\
&= \Phi_l(t, \mathbf{j}^{f,k}, \boldsymbol{\omega}^{f,k}) \tag{3.23} \\
&= \sum_{1 \le i \le n} b_{\alpha^{-1}(i)}^{j^{f,k}} (\boldsymbol{\omega}^{f,k}) \eta_{i+(l-1)n}^V
\end{align*}
$$
où $\omega^{f,k} = (\omega_1^{f,k}, \omega_2^{f,k}, \omega_3^{f,k})$ sont les coordonnées paramétriques constantes, au sein de la sous-variété d'indice $\mathbf{j}^{f,k} = (j_1^{f,k}, j_2^{f,k}, 1)$ du volume englobant, du point de contrôle $k$ appartenant au fil $f$.

De par l'équation (3.23), il apparaît que la fonction $\gamma^f : \mathbb{R}^{3nV} \to \mathbb{R}^{3nf}$ est une application linéaire, on note alors $\Gamma^f \in M_{3nF, 3nV}(\mathbb{R})$ sa représentation matricielle. De plus, les coefficients de $\Gamma^f$ associés aux coordonnées $\eta_{i+(l-1)n}^V$ dans l'expression (3.23) sont identiques pour $l=1,2,3$. Cette particularité implique, au vu de la structure des vecteurs $\eta^V$ et $\eta^f$, que la matrice $\Gamma^f$ est constituée de trois blocs identiques :

$$ \eta^f = \Gamma^f \eta^V \quad \text{où} \quad \Gamma^f = \begin{pmatrix} \tilde{\Gamma}^f \\ \tilde{\Gamma}^f \\ \vdots \\ \tilde{\Gamma}^f \end{pmatrix} \quad \text{avec} \quad \tilde{\Gamma}^f \in M(\mathbb{R}^{nf}, \mathbb{R}^{nV}), \quad \tilde{\Gamma}_{pq}^f = b_{\alpha^{-1}(q)}^{i^{f,p}} (\omega^{f,p}), \quad (3.24) $$

Si l'on nomme $\eta$ le vecteur correspondant à la concaténation des coordonnées généralisées de chacun des fils, on obtient l'expression suivante :

$$ \eta = \left( \underbrace{\eta_{1,1}, \dots, \eta_{1,n_1}}_{\text{fil 1}}, \dots, \underbrace{\eta_{1,1}, \dots, \eta_{1,n_1}}_{\text{fil r}}, \dots, \underbrace{\eta_{2,1}, \dots, \eta_{2,n_2}}_{\text{fil 1}}, \dots, \underbrace{\eta_{2,1}, \dots, \eta_{2,n_2}}_{\text{fil r}}, \dots, \underbrace{\eta_{3,1}, \dots, \eta_{3,n_3}}_{\text{fil 1}}, \dots, \underbrace{\eta_{3,1}, \dots, \eta_{3,n_3}}_{\text{fil r}} \right)^T $$

Le vecteur des configurations $\eta$ correspond, en reprenant la terminologie de la section 2, au paramétrage initial du système matériel. Et la fonction de plongement du liage $\Gamma$ définie en (3.25) représente la fonction de réduction de paramétrage.

$$ \eta = \Gamma \eta^V \quad \text{avec} \quad \Gamma = \begin{pmatrix} \tilde{\Gamma}^1 \\ \vdots \\ \tilde{\Gamma}^{n_r} \\ \vdots \\ \tilde{\Gamma}^{n_r} \\ \vdots \\ \tilde{\Gamma}^{n_r} \end{pmatrix} \quad (3.25) $$

## 3.2 Validation qualitative

Il semble clair que cette méthode de réduction de paramétrage n'est viable que si elle préserve suffisamment de degrés de liberté significatifs au système matériel, de sorte que son évolution dans ce nouvel espace cinématique réduit soit assez riche pour permettre des mouvements comparables, à l'échelle macroscopique, à ceux fournis par le paramétrage initial. Ceci est parfaitement envisageable pour le cas des textiles tricotés dans la mesure où l'espace initial des
configurations est largement sous exploité du fait de l'existence de nombreuses contraintes de contact. La réduction de paramétrage s'effectue alors au détriment des petites déformations locales des boucles, possibles dans le paramétrage initial mais considérées négligeables dans le paramétrage réduit. Afin de vérifier cette affirmation, j'entreprend tout d'abord une validation de nature qualitative. J'impose au volume englobant V plusieurs déformations de nature géométrique. Les mailles du tricot sont alors reconstruites au sein du volume déformé par le biais de la fonction de plongement du liage.

Cette méthode permet d'apprécier la pertinence du choix de la fonction de paramétrage (Figure 55, Figure 56, Figure 57 et Figure 58) dans la mesure où les déformations imposées au tricot ne viennent pas détériorer la structure du liage et conduisent à des géométries de boucles convaincantes. Cette première vérification concluante permet d'envisager de poursuivre la démarche en déduisant le système différentiel réduit, dont les inconnues correspondent aux coordonnées généralisées du volume V, à partir des éléments constitutifs du système différentiel correspondant au paramétrage initial.

<figure><img src="image_68.png" /><figcaption>Figure 55 Traction dans le sens des colonnes et des rangées</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>Figure 56 Cisaillement dans le sens des colonnes et des rangées</figcaption></figure>
<figure><img src="image_70.png" /><figcaption>Figure 57 Flexion dans le sens des rangées et des colonnes</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>Figure 58 Torsion dans le sens des rangées et des colonnes</figcaption></figure>

## 3.3 Construction du système différentiel « réduit »

### 3.3.1 Matrice des masses généralisées

Je commence par expliciter la structure de la matrice des masses généralisées pour le paramétrage réduit $\eta^V$. Pour cela, je rappelle l'expression de $\mathbf{M}^f$ ($1 \le f \le r$), la matrice des masses généralisées associée à chaque fil $f$. Les propriétés structurelles de $\mathbf{M}^f$ ont été détaillées dans le chapitre 1 : Géométrie des variétés splines.

$$ \mathbf{M}^f \in M_{3nf}(\mathbb{R}) \quad ; \quad \mathbf{M}^f = \begin{pmatrix} \mathbf{m}^f \\ \mathbf{m}^f \\ \mathbf{m}^f \end{pmatrix} \qquad (3.26) $$

Si on note $\mathbf{M}$ la matrice des masses généralisées associée au paramétrage initial $\eta$, concaténation des coordonnées généralisées des fils, on a :
$$
\mathbf{M} \in M_{3, \sum_{f}^{nf}}(\mathbb{R}) \ ; \ \mathbf{M} = \begin{pmatrix} \mathbf{m} \\ \mathbf{m} \\ \mathbf{m} \end{pmatrix} \quad \text{avec} \quad \mathbf{m} = \begin{pmatrix} \mathbf{m}^1 \\ \vdots \\ \mathbf{m}^r \end{pmatrix} \qquad (3.27)
$$

De par l'expression particulière de la fonction Γ (3.25), la matrice des masses généralisées **M**^V associée au paramétrage réduit **η**^V s'exprime comme suit :

$$
\mathbf{M}^V \in M_{3nV}(\mathbb{R}) \ ; \quad \mathbf{M}^V = \Gamma^T \mathbf{M} \Gamma = \begin{pmatrix} \mathbf{m}^V \\ \mathbf{m}^V \\ \mathbf{m}^V \end{pmatrix} \quad \text{avec} \quad \mathbf{m}^V = \sum_{i=1}^r \tilde{\Gamma}^{iT} \mathbf{m}^i \tilde{\Gamma}^i \quad (3.28)
$$

Comme la fonction de réduction de paramétrage est linéaire, sa jacobienne représentée par la
matrice Γ est constante. Cette invariance de la jacobienne garantit que la matrice **M**^V est elle
aussi constante. Ce résultat permet de calculer par avance la matrice **M**^V ainsi que sa
décomposition lors de l'initialisation afin d'accélérer la résolution du système réduit.

**3.3.2 Puissance virtuelle des efforts**

D'autre part, le vecteur **Q**^V des puissances virtuelles des efforts exprimés dans le paramétrage réduit **η**^V s'obtient de la façon suivante :

$$
\begin{align}
\mathbf{Q}^V &= \Gamma^T \mathbf{Q} \quad \text{avec} \quad \mathbf{Q} = (\mathbf{Q}^{1^T} \cdots \mathbf{Q}^{r^T})^T \tag{3.29} \\
\mathbf{Q}^V &= \sum_{f=1}^r \Gamma^{f^T} \mathbf{Q}^f \nonumber
\end{align}
$$

Le vecteur **Q**^f correspond au vecteur des puissances virtuelles des efforts appliqués au fil f,
exprimées dans le paramétrage initial **η**^f.

**3.4 Exemple numérique**

On considère un volume englobant V constitué de (5,5,2) points de contrôle, englobant un
morceau de tricot de 15 rangées et 15 colonnes de maille jersey. A raison de 8nc+3 points de
contrôle et de 2nc contraintes vectorielles de contact fil à fil par rangée pour nc colonnes de
jersey, on obtient les données suivantes.
<table><thead><tr><th></th><th>Tricot<br>(paramétrage initial)</th><th>Volume englobant<br>(paramétrage réduit)</th></tr></thead><tbody><tr><td>Degrés de liberté par fil</td><td>3×(8×15+3)=369</td><td></td></tr><tr><td>Total des degrés de liberté</td><td>15×369=55535</td><td>3×(5×5×2)=150</td></tr><tr><td>Contraintes scalaires par fil</td><td>3×(2×15)=90</td><td></td></tr><tr><td>Total des contraintes scalaires</td><td>15×90=1350</td><td></td></tr><tr><td>Total des inconnues scalaires</td><td>55535+1350=6886</td><td>150</td></tr></tbody></table>

Lorsque l'on s'intéresse au nombre d'inconnues scalaires intervenant dans chaque paramétrage, on s'aperçoit que la méthode de réduction de paramétrage mise en place induit une diminution du nombre d'inconnues d'un rapport à peu près égal à 45. Cette réduction drastique du nombre d'inconnues à gérer lors de la résolution garantit comme on l'espérait un gain très significatif en terme de temps de calculs. Les 4 images de la Figure 59 illustrent, malgré la réduction importante des degrés de liberté du système, la préservation de la richesse du mouvement.

<figure><img src="image_72.png" /></figure>
<figure><img src="image_73.png" /><figcaption>Figure 59 Images issues d'une animation d'un volume englobant composé de (2,2,1) sous-variétés, englobant un échantillon de 15 colonnes et 15 rangées de maille jersey.</figcaption></figure>

Les deux images de la Figure 60 illustrent une autre particularité de la méthode de réduction de paramétrage. En effet, même si les degrés de liberté du tricot correspondent aux coordonnées des points de contrôle du volume englobant V, les efforts qui lui sont appliqués sont une ré- interprétation dans le paramétrage réduit des efforts appliqués aux fils. Cet aspect permet de rendre compte de déformations non uniformes au sein du volume englobant pourtant régulier.

<figure><img src="image_74.png" /><figcaption>Figure 60 Images issues d'une animation d'un échantillon suspendu, lesté par des masses disposées de façon asymétrique.</figcaption></figure>

# 4 Conclusion

Pour résumer, j'ai présenté dans ce chapitre une application de la méthode de changement de paramétrage visant à mettre en place les bases théoriques d'une nouvelle méthode de niveau de détail mécanique compatible avec mon modèle dynamique continu. Cette méthode ne procède
pas par raffinement ou appauvrissement de la discrétisation du système matériel étudié mais
plutôt par réduction de ses degrés de liberté, conformément au formalisme lagrangien. Cette
méthode permet d'obtenir un système d'équations beaucoup plus petit et par conséquent de
diminuer les temps de calcul pour la production d'animations. Mais le mécanisme de niveau de
détail sur des structures discrètes est en règle générale adaptatif, permettant ainsi de faire varier
dynamiquement la taille du système à résoudre. Pour le moment statique, cette méthode doit
evoluer afin de me permettre d'adapter dynamiquement le nombre de sous-variétés constituant
le volume englobant V, modifiant le nombre de degrés de liberté des boucles selon les
configurations et les efforts appliqués.
<figure><img src="image_75.png" /></figure>
<figure><img src="image_76.png" /></figure>

# CONCLUSION GENERALE
<figure><img src="image_77.png" /></figure>
Ce dernier chapitre est l'occasion de dresser le bilan des travaux réalisés durant ces quelques années de recherche. Comme je l'ai mentionné en introduction, les objectifs de cette thèse étaient doubles de par la diversité des contextes de l'époque.

Tout d'abord, l'importante restructuration des travaux de l'équipe d'animation dynamique du LERI m'encouragea à participer à cet effort commun en entreprenant la rédaction de plusieurs travaux bibliographiques. Ces travaux qui constituent une bonne partie de ce mémoire de thèse sont au nombre de trois.

La première étude bibliographique a pour thème une justification des équations de Lagrange. Ce sujet, bien qu'en apparence très éloigné des considérations informatiques, est néanmoins d'une importance capitale pour tout informaticien désireux de mettre en œuvre une application d'animation s'appuyant sur les lois de la physique. Mais je ne cherche pas à donner une démonstration mathématique rigoureuse de ces équations. Je cherche plutôt à guider le lecteur dans le développement de ce formalisme : je m'emploie à dégager son sens physique en illustrant chaque concept fondateur. Cette étude a pour vocation d'éviter au lecteur de buter sur les nombreuses difficultés ou ambiguïtés des équations de Lagrange et de lui fournir tous les éléments nécessaires à une implantation informatique des équations de Lagrange.

La deuxième étude qui traite des problèmes d'intégration numérique, autre discipline étrangère à l'informatique, est rédigée avec le même souci didactique. La résolution d'équations différentielles de nature quelconque s'appuie sur des méthodes numériques pour déterminer une solution approchée. Mais les méthodes existantes sont nombreuses et ne présentent pas toutes les mêmes avantages ou inconvénients. Une classification se révèle alors précieuse pour aider l'informaticien à choisir la méthode d'intégration la plus appropriée au problème qu'il cherche à résoudre. J'expose alors plusieurs critères de comparaison avant d'introduire les méthodes les plus couramment utilisées. L'intérêt de cette étude réside dans la synthèse de ces méthodes, répertoriées selon leur nature, leurs avantages et inconvénients respectifs.

Enfin la dernière étude de cette première partie consiste en un état de l'art des modèles déformables pour l'animation dynamique. Cette étude, tout à fait classique pour un mémoire de thèse, permet de dresser une cartographie des différentes approches retenues pour l'animation dynamique d'objets non rigides. Cette étude a le mérite de présenter les avantages et les inconvénients et de justifier la nécessité du modèle que je présente dans la deuxième partie de ce mémoire.

La deuxième partie de ce document, initialement motivée par le contexte applicatif de mon travail, développe les étapes de l'élaboration d'un modèle mécanique générique pour représenter un milieu continu déformable. En effet, mon travail de thèse s'inscrit dans la continuité de celui de Jean-Michel Nourrit qui consista à proposer un modèle informatique pour la représentation, l'animation et la visualisation de textiles tricotés. Ce modèle repose principalement sur un modèle de courbes splines qui correspondent aux fils constitutifs du tricot. Ce modèle que nous qualifions de courbe spline matérielle est en fait à l'origine du modèle innovant que je propose. En effet, animé par cette volonté d'unification et de généralisation des travaux de l'équipe, j'entrepris une refonte du modèle de Jean-Michel afin qu'il réponde aux nouvelles exigences du
moteur d'animation dynamique unifié mais aussi qu'il permette de représenter une plus grande
gamme de corps déformables. La première étape de cette unification consiste en la proposition
d'une définition de courbe spline indépendante de la nature des fonctions de mélange.
L'introduction de la notion de segment spline permet de généraliser la propriété de localité
(réalisée par certaines courbes splines) à toutes les courbes splines sans distinction. La deuxième
étape débouche sur une définition de variété spline qui vient unifier les définitions très
similaires de courbes, surfaces et volumes splines. Ce nouveau modèle rassemble une grande
variété de formes d'allures et de dimensions quelconques. Ce modèle de corps déformable se
distingue des autres par sa capacité à préserver la continuité du milieu qu'il représente. D'autre
part, sa simplicité (la position d'un point de la variété spline est obtenue par combinaison
linéaire d'un nombre réduit de points de contrôle pondérée par des fonctions de mélange)
permet de réaliser de nombreuses simplifications des équations du mouvement, débouchant sur
l'expression d'un système différentiel équivalent et formellement exact. La préservation de la
continuité demeure une préoccupation centrale lorsque j'étend l'expression du modèle de
variété spline à la dynamique. Je détaille d'ailleurs la construction d'une répartition massique
continue paramétrique à partir de la définition d'une densité de masse par unité de longueur,
plus intuitive pour l'utilisateur du modèle. La nature même du modèle autorise l'expression de
forces ponctuelles en des points quelconques de la variété du fait de sa nature continue.
Contrairement aux approches couramment répandues, je ne procède pas à une discrétisation du
milieu pour déboucher sur une expression simplifiée des équations du mouvement. La seule
discrétisation persistante que je souhaite éliminer réside dans la modélisation des efforts
internes d'élasticité. Dans un premier temps, j'utilisais des ressorts régulièrement répartis au
sein de la variété pour représenter les contraintes internes. Mais, fidèle à cette volonté d'aboutir
à la définition d'un modèle mécanique continu, je développe mes efforts pour construire une
énergie continue de déformation issue de la théorie de l'élasticité en déplacements finis. Ces
premiers travaux se limitent pour le moment au cas des courbes splines. Mais les résultats
satisfaisants de par leur adéquation à la définition des matériaux et le surcoût négligeable en
temps de calculs m'encouragent à continuer dans cette voie.

Une fois le modèle de variété spline construit, je m'oriente de nouveau vers la simulation de
tricots. J'exploite le modèle générique proposé afin de réaliser certaines optimisations sur le
modèle initial de tricots pour équilibrer le souci de réalisme avec les besoins de performance. Je
mets à profit la méthode de réduction de paramétrage en plongeant les mailles du tricot dans un
volume spline englobant. Le nouveau système matériel que constitue ce volume est doté d'un
nombre réduit de points de contrôle et est dénué de toute contrainte de contact entre fils. Cette
approche permet d'aboutir à un système d'équations de taille très réduite en comparaison avec
celui déduit d'une modélisation à l'échelle du fil.

Nombre des objectifs ont alors été atteints durant cette thèse mais naturellement pas tous. Ainsi,
ce travail mérite d'être prolongé pour atteindre une véritable maturité. Parmi les nombreuses
perspectives que l'on peut envisager, voici celles qui me semblent nécessaires à l'évolution du
modèle de variétés splines :

Le nombre de fonctions de mélange actuellement implantées au sein du moteur d'animation
dynamique est relativement faible et se limite essentiellement aux fonctions cubiques
d'interpolation. L'architecture de classes a été conçue pour que l'ajout de nouvelles fonctions se fasse de façon naturelle. Ne reste plus alors qu'à retrousser ses manches pour entreprendre leur codage et enrichir ainsi la palette de formes envisageables au sein du modèle. Toujours dans le but d'enrichir les formes réalisables, j'envisage d'autres configurations de réseau de contrôle que la configuration parallélépipédique jusqu'à là utilisée. En utilisant, des objets splines circulaires dont les extrémités partagent les mêmes points de contrôle, on peut alors naturellement représenter des objets sphériques, cylindriques ou toriques. Une autre piste consiste à travailler non pas sur le réseau des points de contrôle mais sur le domaine paramétrique qui identifie les particules du corps déformable représenté par une variété spline. Ces approches non encore implantées sont sérieusement envisagées.

Une autre direction de recherche issue du modèle de variété spline est l'enrichissement de la définition de l'énergie continue de déformation dont découlent les efforts internes d'élasticité. En effet, ce travail qui s'inscrit dans une démarche rigoureuse et générale se limite dans ce mémoire au cas des objets linéiques. Ce modèle ayant fait ses preuves, l'étape de validation est satisfaisante et me permet d'envisager son extension au cas des dimensions supérieures. Je suis conscient que cette extension va engendrer un surcoût important du fait de l'explosion combinatoire. Pour cette raison, j'envisage aussi de rendre l'expression de la loi de comportement plus modulaire afin de pouvoir régler l'ordre des termes intervenant dans son expression.

De façon plus pragmatique, les efforts réalisés pour accroître la performance de modèle de textiles tricotés sont satisfaisants mais doivent être approfondis. La méthode de réduction de paramétrage permet de réduire de façon importante le nombre d'inconnues intervenant dans les équations du mouvement. Mais cette réduction est statique et peut se révéler inadapté pour rendre compte d'une déformation importante très localisée d'un échantillon textile. Pour que cette méthode puisse devenir une véritable méthode de niveau de détail mécanique, il nous faut envisager la possibilité d'adapter dynamiquement les degrés de liberté du volume pour pouvoir rendre compte plus finement des déformations locales.
<figure><img src="image_78.png" /></figure>
BIBLIOGRAPHIE

[Arnold 1989] V.I. Arnold, *Mathematical Methods of Classical Mechanics* (Second Edition). Springer Verlag, 1989.

[Bathe 1982] K.L. Bathe, *Finite Element Procedures in Engineering Analysis*. Prentice Hall, Inc., Enlewood Cliffs, New Jersey, 1982.

[Blanc 1994] C. Blanc, *Techniques de Modélisation et de Déformation de Surfaces pour la Synthèse d'Images*. Thèse de doctorat en informatique, Université de Bordeaux I, 1994.

[Baraff 1992] D. Baraff, A. Witkin, *Dynamic simulation of non-penetrating flexible bodies*. SIGGRAPH'92 Conference Proceedings, 26(4) p303-308, 1992.

[Baraff 1996] D. Baraff, *Linear-time dynamics using Lagrange multipliers*. Computer Graphics Proceedings, Annual Conference Series, p137-146, 1996.

[Baraff 1998] D. Baraff, A. Witkin, *Large steps in cloth simulation*. Computer Graphics Proceedings, Annual Conference Series, p43-54, 1998.

[Barnhill 1974] R. Barnhill, R. Riesenfeld, *Computer aided Geometric design*. Academic Press, 1974.

[Barroso Filho 1994] W. Barroso Filho, *La mécanique de Lagrange (Principes et méthodes)*. Editions Karthala, Paris, 1994.

[Barsky 1981] B. Barsky, *The Beta-spline: a local representation based on shape parameters and fundamental geometric measures*. PhD Thesis, University of Utah, 1981.

[Barsky 1983] B. Barsky, J. Beatty, *Local control of bias and tension in Beta-splines*. SIGGRAPH'83 Conference Proceedings, 17(3) p193-218, 1983.

[Baumgarte 1972] J. Baumgarte, *Stabilization of constraints and integrals of motion*. Computer Meth. Appl. Mech. Eng. (1) p1-16, 1972

[Benz 1988] W. Benz, *Application of Smooth Particle Hydrodynamics (SPH) to Astrophysical Problems*. Comp. Phys. Comm., (48) p97-105, 1988.

[Bérest 1997] P. Bérest, *Calcul des variations, application à la mécanique et la physique*. X Ecole Polytechnique. Ellipses, 1997.

[Berger 1987] M. Berger, B. Gostiaux, *Géométrie différentielle: variétés, courbes et surfaces*. Presses Universitaires de France, 1987.

[Bézier 1966] P. Bézier, *Définition numérique de courbes et surfaces*, Automatisme, (11) p 625-632 (partie 1), (12) p 17-21 (partie 2), 1966.
[Bézier 1977] P. Bézier, *Essai de définition numérique des courbes et surfaces expérimentales*. Thèse de doctorat, Université Paris VI, 1977.

[Bourguignon 2000] D. Bourguignon, M.P. Cani, *Controlling Anisotropy in Mass-Spring Systems*. Proceedings of the 11th Eurographics Workshop on Computer Animation and Simulation. Springer-Verlag, 2000.

[B. de Casson 2000] F. Boux de Casson, C. Laugier. *Simulating 2D tearing phenomena for interactive medical surgery simulators*. Computer Animation Proceedings, 2000.

[Breen 1992] D.E. Breen, D.H. House, P.H. Getto, *A physical based particle model of woven cloth*. The Visual Computer, (8) p264-277, 1992.

[Breen 1994] D.E. Breen, D.H. House, M. Wozny, *Predicting the drape of woven cloth using interacting particles*. Computer Graphics Proceedings, Annual Conference Series, p365-372, 1994.

[Bro-Nielsen 1996] M. Bro-Nielsen, S. Cotin, *Real-time volumetric deformable models for surgery simulation using finite elements and condensation*. Eurographics'96 Proceedings, p21-30, 1996.

[Butcher 1964] J.C. Butcher, *On Runge-Kutta processes of high order*. J. Australian Math. Soc., (4) p179-194, 1964.

[Calkin 1996] M.G. Calkin, *Lagrangian & Hamiltonian mechanics*. World Scientific, 1996.

[Cartwright 1992] J.H.E. Cartwright, O. Piro, *The Dynamics of Runge-Kutta Methods*. Int. J. Bifurcation and Chaos (2) p427-449, 1992.

[Catmull 1974] E. Catmull, R. Rom, *A class of local interpolating splines*. Dans [Barnhill 74] p317-326.

[Ciarlet 1986] P.G. Ciarlet, *Elasticité tridimensionnelle*. Masson, Paris, 1986.

[Collier 1991] J. Collier, B. Collier, G. O'Toole, S. Sargand, *Drape prediction by means of finite elements analysis*. Journal of the Textile Institute, 82(1) p96-107, 1991.

[Coquillart 1990] S. Coquillart, *Extended free form deformation: a sculpturing tool for 3D geometric modeling*. SIGGRAPH'90 Conference Proceedings, 24(4) p187-193, 1990.

[Coquillart 1991] S. Coquillart, P. Jancene, *Animated FreeForm Deformation: An Interactive Animation Technique*. SIGGRAPH'91 Conference Proceedings, 25(4) p23-26, 1991.

[Cotin 1997] S. Cotin, *Modèles anatomiques déformables en temps-réel*. Thèse de doctorat en Sciences de l'Ingénieur, Université de Nice-Sophia Antipolis, 1997.
[Cozot 1996] R. Cozot, *Environnement de simulation de systèmes physiques : modèle et langage*. Thèse de doctorat en informatique, Université de Rennes I, 1996.

[Crouzeix 1989] M. Crouzeix, A. Mignot, *Analyse numérique des équations différentielles, 2ᵉ édition*. Masson, Paris, 1989.

[Dacorogna 1992] B. Dacorogna, *Introduction au calcul des variations*. Cahiers mathématiques de l'Ecole Polytechnique Fédérale de Lausanne. Presses polytechniques et universitaires romandes, Lausanne, 1992.

[Debunne 2000] G. Debunne, *Animation multirésolution d'objets déformables en temps-réel, Application à la simulation chirurgicale*. Thèse de doctorat spécialité Imagerie, Vision et Robotique, Institut National Polytechnique de Grenoble, 2000.

[Demailly 1991] J.P. Demailly, *Analyse numérique et équations différentielles*. Editions Presses Universitaires de Grenoble, 1991.

[Desbrun 1997] M. Desbrun, *Modélisation et animation d'objets hautement deformables à l'aide de surfaces implicites*. Thèse de doctorat spécialité Imagerie, Vision et Robotique, Institut National Polytechnique de Grenoble, 1997.

[Dieudonné 1968] J.A. Dieudonné, *Calcul infinitésimal*. Hermann, 1968.

[Eberhardt 1996] B. Eberhardt, A. Weber, W. Strasser, *A Fast, Flexible, Particle-System Model for Cloth Draping*. IEEE Computer Graphics and Applications, 16(5) p52-59, 1996.

[Eberhardt 1999] B. Eberhardt, A. Weber, *A particule system approach to knitted textiles*. Computers & Graphics, (23) p599-606, 1999.

[Eischen 1996] J.W. Eischen, S. Deng, T.G. Clapp, *Finite-Element Modeling and Controll of Flexible Fabric Parts*. IEEE Computer Graphics and Applications, 16(5) p71-80, 1996.

[Enright 1995] W.H. Enright, D.J. Higham, B. Owren, P.W. Sharp, *A survey of the Explicit Runge-Kutta method*. Technical Report, Computer Sciences Department, University of Toronto, 1995

[Farin 1990] G. Farin, *Curves and surfaces for Computer Aided Geometric Design: a practical guide*. Academic Press, second edition, 1990.

[Fehlberg 1974] E. Fehlberg, *Classical seventh, sixth and fifth-order Runge-Kutta-Nyström formulas with step size control for general second order differential equations*. Tech. Report 432, NASA, 1974

[Foley 1989] J. Foley, A. Van Dam, S. Feiner, J. Hughes, *Computer Graphics: principles and practice*, Addison Wesley, 1989.

[Fox 1967] E. A. Fox, *Mechanics*. Rensselaer Polytechnic Institute. A Harper International Edition, jointly published by Harper & Row, New York, Evanston & London and John Weatherhill, Inc., Tokyo, 1967.
[Fox 1987] C. Fox, *An Introduction to the Calculus of Variations*. Dover Publications, New York, 1987.

[Garrigues 2000a] J. Garrigues, *Mécanique des milieux continus*. 2000. http://esm2.imt-mrs.fr/gar/mmc.html

[Garrigues 2000b] J. Garrigues, *Mécanique des milieux continus en déformations finies*. 2000. http://esm2.imt-mrs.fr/gar/gd.html

[Garrigues 2001] J. Garrigues, *La méthode des éléments finis*. 2001. http://esm2.imt-mrs.fr/gar/ef.html

[Germain 1986] P. Germain, *Mécanique X, Tome 1*. Ecole Polytechnique, Ellipses, 1986.

[Giaquinta 1996] M. Giaquinta, S. Hildebrandt, *Calculus of Variations I*. Springer Berlin, 1996.

[Gillard 1997] D. Gillard, *Vers l'animation dynamique d'un corps humain virtuel : utilisation de données biomédicales 3D*. Thèse de doctorat en informatique, Université de Reims Champagne-Ardenne, 1997.

[Gourret 1989] J.P. Gourret, N. Magnenat Thalmann, D. Thalmann, *Simulation of object and human skin deformation in a grasping task*. SIGGRAPH'89 Conference Proceedings, 23(4) p21-30, 1989.

[Hairer 1991] E. Hairer, G. Wanner, *Solving Ordinary Differential Equations II Stiff and Differential Algebraic Problems*. Springer Series in Comput. Mathematics, (14), 1991, Second revised edition 1996.

[Hairer 1993] E. Hairer, *Introduction à l'Analyse Numérique*. 1993. http://www.unige.ch/math/folks/hairer/polycop.html

[Hauth 2001] M. Hauth, O. Etzmuss, *High Performance Solver for the Animation of Deformable Objects using Advanced Numerical Methods*. Computer Graphics forum, 20(3) Conference Issue p320-328, 2001.

[Hughes 1987] T.J. Hughes, *The Finite Element Method : Linear Static and Dynamic Finite Element Analysis*. Prentice-Hall, Inc., Enlewood Cliffs, New Jersey, 1987.

[Hutchinson 1996] D. Hutchinson, M. Preston, T. Hewitt, *Adaptive refinement of mass/spring simulation*. Proceedings of the 7th Eurographics workshop on Computer Animation and Simulation, p31-35, 1996.

[Jackson 1990] K.R. Jackson, S.P. Norsett, *The potential for Parallelism in RK methods*. Technical Report No 239/90, Computer Sciences Department, University of Toronto, 1990.

[Kawabata 1980] S. Kawabata, *The Standardization and Analysis of Hand Evaluation*. Textile Machinery Soc. Of Japan, Osaka, 1980.

[Kochanek 1984] D. Kochanek, R. Bartels, *Interpolating splines with local tension, continuity and bias control*. SIGGRAPH'84 Conference Proceedings, 18(3) p33-41, 1984.
[Kutta 1901] W. Kutta, *Beitrag zur näherungsweisen Integration totaler Differentialgleichungen*. Zeitschr. für Math. u. Phys., (46) p435-453, 1901.

[Lambert 1991] J.D. Lambert, *Numerical Methods for Ordinary Differential Systems: The Initial Value Problem*. John Wiley, 1991.

[Lawson 1966] J.D. Lawson, *An order six Runge-Kutta process with extended region of stability*. SIAM J. Numer. Anal., (3) p593-597, 1966.

[Luciani 1985] A. Luciani, *Un outil informatique de création d'images animées*. Thèse de docteur ingénieur d'électronique, Institut National Polytechnique de Grenoble, 1985.

[Meißner 1998] M. Meißner, B. Eberhardt, *The Art of Knitted Fabrics, Realistic & Physically Based Modelling of Knitted Patterns*. EUROGRAPHICS'98 Proceedings, 17(3) p355-362, 1998.

[Miller 1988] G. Miller, *The motion dynamics of snakes and worms*. SIGGRAPH'88 Conference Proceedings, 22(4) p169-177, 1988.

[Nocent 2001a] O. Nocent, J.M. Nourrit, Y. Remion, *Towards mechanical level of detail for knitwear simulation*. WSCG'2001 Proceedings, p252-p259, 2001.

[Nocent 2001b] O. Nocent, J.M. Nourrit, Y. Remion, *Vers du niveau de détail mécanique pour l'animation dynamique d'objets continus, Application aux textiles tricotés*. Revue internationale de CFAO et d'Informatique Graphique (16) 1 p115-130, 2001.

[Nocent 2001c] O. Nocent, Y. Remion, *Continuous deformation energy for Dynamic Material Splines subject to finite displacements*. Proceedings of the 12th Eurographics workshop on Computer Animation and Simulation, p87-98, 2001.

[Nourrit 1999] J.M. Nourrit, *Modélisation, animation et visualisation de textiles à base de mailles*. Thèse de doctorat en informatique, Université de Reims Champagne-Ardenne, 1999.

[O'Brien 1999] J.F. O'Brien, J.K. Hodgins, *Graphical Modeling and Animation of Brittle Fracture*. Computer Graphics Proceedings, Annual Conference Series, p137-146, 1999.

[Palazzi 1994] L.F. Palazzi, D.R. Forsey, *A Multilevel Approach to Surface Response in Dynamically Deformable Models*. Proceedings of the 5th Eurographics workshop on Computer Animation and Simulation, p21-30, 1994.

[Pentland 1989] A. Pentland, J. Williams, *Good Vibrations: Modal Dynamics for Graphics and Animation*. SIGGRAPH'89 Conference Proceedings, 23(4) p215-222, 1989.

[Picinbono 2001] G. Picinbono, *Modèles géométriques et physiques pour la simulation d'interventions chirurgicales*. Thèse de doctorat en Sciences de l'Ingénieur, Université de Nice-Sophia Antipolis, 2001.
[Piegl 1991] L. Piegl, *On NURBS: a survey*. Computer Graphics and Application, 11(1) p55-71, 1991.

[Platt 1981] S. Platt, N. Badler, *animating facial expressions*. SIGGRAPH'81 Conference Proceedings, 15(3) p245-252, 1981.

[Press 1988] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, *Numerical Recipes in C*. Cambridge University Press, 1988.

[Provot 1995] X. Provot, *Deformation constraints in a mass-spring model to describe rigid cloth behaviour*. Graphics Interface Proceedings, p147-154, 1995.

[Remion 2000] Y. Remion, *Animation dynamique : moteur lagrangien généraliste et applications*. Habilitation à diriger des recherches, Université de Reims Champagne-Ardenne, 2000.

[Remion 2001] Y. Remion, J.M. Nourrit, O. Nocent, *D-dimensional parametric models for dynamic animation of deformable objects*. The Visual Computer Journal, (17) 3 p167-178, 2001.

[Riesenfeld 1973] R. Riesenfeld, *Applications of B-spline Approximation to geometric problems of Computer Aided Design*. PhD Thesis, University of Syracuse, 1973.

[Runge 1895] C. Runge, *Ueber die numerische Auflösung von Differentialgleichungen*. Math. Ann., (46) p167-178, 1895.

[Sederberg 1986] T. Sederberg, S. Parry, *Free-Form Deformation of Solid Geometric Models*, SIGGRAPH'86 Conference Proceedings, (20)4 p151-160, 1986.

[Souchet 2001] R. Souchet, *Leçons sur les grandes déformations*. Cépaduès Editions, Toulouse, 2001.

[Stroud 1974] A. H. Stroud, *Numerical Quadrature and Solution of Ordinary Differential Equations*. Applied Mathematical Sciences 10, Springer-Verlag, 1974.

[Terzopoulos 1987] D. Terzopoulos, J. Platt, A. Barr, K. Fleischer, *Elastically deformable models*. SIGGRAPH'87 Conference Proceeding, 21(4) p205-214, 1987.

[Terzopoulos 1988] D. Terzopoulos, K. Fleischer, *Modeling Inelastic Deformation: Viscoelasticity, Plasticity, Fracture*. SIGGRAPH'88 Conference Proceedings, 22(4) p269-278, 1988.

[Terzopoulos 1989] D. Terzopoulos, J. Platt, K. Fleischer, *Heating and melting deformable models (from goop to glop)*. Graphics Interface Proceedings, p219-226, 1989.

[Van Gelder 1998] A. Van Gelder, *Approximate simulation of elastic membranes by triangulated spring meshes*. Journal of Graphic Tools, 3(2) p21-41, 1998.

[Volino 1995] P. Volino, M. Courchesne, N. Magnenat Thalmann, *Versatile and Efficient Techniques for Simulating Cloth and Other Deformable Objects*.
Computer Graphics Proceedings, Annual Conference Series, p137-144, 1995.

[Volino 2000] P. Volino, N. Magnenat Thalmann, *Implementing Fast Cloth Simulation with Collision Response*. Computer Graphics International Proceedings, 2000.

[Witkin 1990] A. Witkin, W. Welch, *Fast Animation and Control of Nonrigid Structures*. SIGGRAPH'90 Conference Proceedings, 24(4) p243-251, 1990.

[Yngve 2000] G.D. Yngve, J.F. O'Brien, J.K. Hodgins, *Animating Explosions*. Computer Graphics Proceedings, Annual Conference Series, p29-36, 2000.

[Zienkiewicz 1977] O. Zienkiewicz, *The finite element method, 3rd edition*. Mac Graw Hill, London, 1977.

[Zonneveld 1964] J.A. Zonneveld, *Automatic Numerical Integration*. Mathematical Centre Tracts, No 8, Matematisch Centrum, Amsterdam, 1964.
<figure><img src="image_79.png" /></figure>
<figure><img src="image_80.png" /></figure>
# Animation dynamique de corps déformables continus
## Application à la simulation de textiles tricotés

Les progrès réalisés en animation par ordinateur ne cessent de repousser les limites de l'ensemble des formes dont on peut reproduire le mouvement. L'évolution de formes changeantes notamment est d'ores et déjà un centre d'intérêt important au sein de la communauté de l'informatique graphique. L'animation dynamique, qui repose sur les équations de la physique, consiste en une des voies envisageables pour représenter la déformation. L'avantage d'une telle approche réside dans l'adéquation du modèle avec les propriétés intrinsèques de la matière décrites par la physique. D'autre part, l'évolution de ces corps étant conditionnée par les équations du mouvement, l'animation dynamique permet d'automatiser le processus de génération d'animations.

Mais l'animation dynamique requiert la maîtrise de compétences étrangères à l'informatique telles que la mécanique et l'analyse numérique. Afin de tenter de répondre à quelques problèmes cruiaux auxquels l'informaticien doit faire face pour la réalisation d'un logiciel d'animation dynamique, la première partie de ce mémoire rassemble des études bibliographiques traitant d'une justification des équations de Lagrange, d'une classification des méthodes usuelles d'intégration numérique et d'un état de l'art des modèles déformables proposés en animation dynamique.

Nombre de modèles de corps déformables, même s'ils reposent sur des équations régissant le comportement des milieux continus, ont recours à une étape de discrétisation. Cette étape se traduit par une répartition de la masse totale du corps parmi un ensemble discret de points caractéristiques. Dans la deuxième partie de ce mémoire, je développe ma contribution à l'animation dynamique en présentant un modèle innovant de corps déformables dont la répartition massique est continue. Ce modèle, qui s'appuie sur une définition spline, englobe à la fois les courbes, les surfaces et les volumes splines. La définition géométrique de ce modèle couplée avec une description paramétrique continue de la masse permet de réaliser des simplifications formelles des équations du mouvement pour déboucher sur un système différentiel équivalent aux propriétés structurelles intéressantes. Je propose aussi une modélisation continue des efforts élastiques internes à partir d'une énergie de déformation paramétrée par les constantes caractéristiques de la matière. Ce travail s'inscrivant dans une démarche générale se limite pour le moment au cas des corps linéiques. Enfin, j'utilise le modèle des variétés splines pour mettre en œuvre une méthode de réduction de paramétrage afin d'accroître les performances de production d'animations dynamiques de textiles tricotés

Mots clés : Animation dynamique, corps déformables, simulation de textiles tricotés.

# Dynamic animation of continuous deformable objects
## Application to knitted cloth simulation

Recent improvements in computer animation allow simulating the motion of a large variety of shapes. As an example, deformable objects animation is becoming a major research topic for the computer graphics community. Dynamic animation is one of the possible ways to model deformation through time. First, a physics-based modelling approach gives a natural mapping between the parameters of the model and the physical properties of matter. Second, the motion considered as the solution of the underlying mechanical equations can be produced automatically with a given set of initial conditions.

But dynamic animation greatly depends on several disciplines, foreign to the computer sciences area like mechanics, numerical analysis. In order to help every developer who wants to build a dynamic animation engine, the first part of the document gathers three bibliographical studies dealing with a justification of the Lagrange equations, a classification of usual integration methods and a state of the art about physical-based modelling of deformable objects.

Many of the proposed mechanical models for deformable objects, even if they are built upon continuum mechanics, proceed in a discretizing step. This step consists of distributing the total mass upon a discrete set of characteristic points. In the second part of this document, I present an innovative model for deformable objects using a continuous mass repartition. This model, based on a spline definition, encompasses the definitions of spline curves, surfaces and volumes and allows performing formal simplifications on the equations of motion to obtain an equivalent linear differential system exhibiting remarkable structural properties. I also propose a continuous model for inner elastic constraints built upon a deformation energy function. This work, developed in a general way, is currently restricted to spline curves. Finally, I take advantage of the properties of this new mechanical model to propose a parameters reduction method in order to speed up the computation times during the process of knitwear simulation.

Keywords: Dynamic animation, deformable objects, knitted cloth simulation.