La méthode des éléments finis est un moyen classique de résoudre les équations de la mécanique des milieux continus. C'est un outil mathématique permettant de discrétiser un problème variationnel continu 〚24〛. Cette méthode est basée sur une subdivision du domaine d'intérêt en volumes élémentaires. Nous avons choisi d'utiliser des éléments finis de type P₁ pour lesquels l'élément de volume est un tétraèdre (Fig. 1). En chaque point $\mathrm{M}\left(\mathrm{x},\mathrm{y},z\right)$ du tétraèdre $𝖳,$ le vecteur déplacement est interpolé à partir des valeurs des déplacements U_k des sommets P_k. Pour les éléments P₁, les fonctions d'interpolation Λ_k sont linéaires ($\left\{{\Lambda }_k;\mathrm{k}=0,\cdots ,3\right\}$ sont les coordonnées barycentriques de M dans le tétraèdre) (éqs. (6) et (7)) :

En utilisant cette équation pour U, on obtient la formulation de l'énergie de déformation élastique linéaire dans le tétraèdre $𝖳$ (éq. (8)) :

Finalement, on obtient la force ${\mathrm{F}}_p^{𝖳}$ appliquée par le tétraèdre $𝖳$ sur le sommet P_p en dérivant l'énergie élastique par rapport à la position de ce sommet (éq. (9)) :

Nous utilisons cette formulation de l'élasticité linéaire depuis plusieurs années au travers de deux modèles déformables, le modèle pré-calculé 〚25〛 et le modèle masse–tenseur 〚8〛. De plus, nous les avons étendus à l'élasticité linéaire anisotrope 〚9〛, ce qui permet de modéliser des matériaux renforcés par des fibres, qui sont très répandus parmi les tissus biologiques (muscles, tendons...), ou d'autres structures anatomiques, comme les vaisseaux sanguins.

La principale limitation du modèle linéaire est qu'il n'est pas invariant par rotation. Lorsqu'un élément de volume subit une rotation, il l'interprète comme une augmentation de l'énergie élastique, qui se traduit par une distorsion. La Fig. 2 montre la déformation que subit un modèle élastique linéaire lorsqu'il est soumis à une rotation.

Lorsqu'il s'agit d'une rotation globale, on pourrait imaginer résoudre le problème par un changement de repère approprié. Mais, dans la plupart des cas, il n'y a qu'une partie de l'objet simulé qui subit une rotation par rapport au reste du domaine. C'est en particulier le cas dès lors qu'une partie de l'objet est fixée et que le reste est soumis à un ensemble de contraintes entraînant une rotation locale. La Fig. 3 montre un cylindre dont la face inférieure a été fixée. La face supérieure est soumise à une force dirigée vers la droite. Nous avons représenté sur la Fig. 3a les déformations entraînées pour plusieurs valeurs successives de l'intensité de la force. On observe, comme sur la Fig. 2, une dilatation du modèle aux endroits où il subit une rotation. Sur la Fig. 3b, nous avons représenté à l'aide de flèches les trajectoires de certains points sur la face supérieure du cylindre. Nous constatons que ces points suivent, au cours de la déformation, une trajectoire linéaire, qui est responsable de la dilatation observée. Cette dilatation est d'autant plus irréaliste qu'elle n'est pas isotrope, puisqu'elle ne se produit que dans le plan de rotation. On le constate en observant la déformation subie par la face supérieure du cylindre (Fig. 3c).

Le comportement irréaliste du modèle d'élasticité linéaire, lorsqu'il subit des grands déplacements, nous a conduit à considérer un modèle d'élasticité très proche du modèle linéaire, mais qui reste valable pour les grandes déformations, le modèle d'élasticité de St Venant–Kirchhoff.

En conservant les conventions de notations utilisées pour l'élasticité linéaire, nous exprimons l'énergie de déformation d'un tétraèdre $𝖳$ pour le modèle d'élasticité de St Venant–Kirchhoff sous la forme (éq. (13)) :

Comme on pouvait s'y attendre, cette énergie comporte maintenant, en plus du terme quadratique de l'élasticité linéaire, des termes d'ordre trois et quatre mettant en œuvre respectivement les déplacements de trois et quatre sommets du tétraèdre. Ces termes font apparaître les nouvelles grandeurs $𝖢{}_{𝖳}^{\mathrm{jkl}}$ et $𝖣{}_{𝖳}^{\mathrm{jklm}},$ qui sont respectivement des vecteurs et des scalaires de rigidité. Ils viennent compléter les tenseurs de rigidité $𝖡{}_{𝖳}^{\mathrm{jk}}$ (éq. (8)) et dépendent des fonctions de forme du tétraèdre (éq. (6)) et des coefficients de Lamé du matériau (éq. (14)) :

La force exercée par le tétraèdre sur chacun de ses sommets P_p est alors dérivée de l'énergie élastique $W\left(𝖳\right)$ pour donner une équation polynomiale de degré trois (éq. (15)) :

On remarque que le premier terme de la force élastique, ${\mathrm{F}}_1^p\left(𝖳\right),$ correspond exactement à la force rencontrée dans le cas de l'élasticité linéaire.

L'idée principale du modèle masse–tenseur est de séparer, pour chaque tétraèdre, la force appliquée sur un sommet en deux parties : une force créée par le déplacement du sommet et une somme de forces dues aux déplacements de ses voisins (éq. (16)) :

De cette façon, nous pouvons définir, pour chaque tétraèdre, un ensemble de tenseurs de rigidité locaux pour les sommets ($\left\{𝖡{}_{𝖳}^{\mathrm{pp}};\mathrm{p}=0,\cdots ,3\right\}$) et pour les arêtes ($\left\{𝖡{}_{𝖳}^{\mathrm{pj}};\mathrm{p},\mathrm{j}=0,\cdots ,3;\mathrm{p}\ne \mathrm{j}\right\}$). Après avoir effectué la même opération pour tous les tétraèdres, on peut accumuler, sur les sommets et les arêtes du maillage, les contributions correspondantes pour former les tenseurs de rigidité globaux :

Ces tenseurs de rigidité sont calculés au moment de la construction du modèle et sont stockés sur chaque sommet et arête du maillage.

Le même principe peut être appliqué au terme quadratique ${\mathrm{F}}_2^p\left(𝖳\right)$ et au terme cubique ${\mathrm{F}}_3^p\left(𝖳\right)$ de l'éq. (15). Le premier apporte des vecteurs de rigidité pour les sommets, les arêtes et les triangles, alors que le second définit des scalaires de rigidité pour les sommets, les arêtes, les triangles et les tétraèdres. Le Tableau 1 résume la répartition de toutes ces informations de rigidité sur chaque primitive géométrique du maillage.

Étant donné un maillage tétraédrique d'un solide – dans notre cas une structure anatomique –, nous construisons une structure de données comprenant les notions de sommet, d'arête, de triangle et de tétraèdre, avec toutes les relations de voisinage nécessaires. Pour chaque sommet, nous stockons sa position actuelle P_p, sa position au repos P⁰_p et ses informations de rigidité. Pour chaque arête et chaque triangle nous stockons les informations de rigidité. Enfin, pour chaque tétraèdre, nous stockons les coefficients de Lamé λ et μ, les quatre vecteurs de forme α_k et les informations de rigidité.

Durant la simulation, nous calculons les forces pour chaque sommet, arête, triangle et tétraèdre, puis nous mettons à jour la position des sommets en considérant que l'équation qui gouverne le mouvement de notre modèle élastique est une équation différentielle newtonienne, de la forme (éq. (17)) :

Cette équation est reliée aux équations différentielles rencontrées en mécanique des milieux continus 〚26〛 (éq. (18)) :

D'après la théorie des éléments finis, la matrice de masse M et la matrice d'amortissement C sont creuses et dépendent des propriétés de chaque tétraèdre. Dans notre cas, nous considérons que M et C sont des matrices diagonales, c'est-à-dire que la masse et les effets d'amortissement sont concentrés sur les sommets du maillage. Cette simplification, appelée mass-lumping, découple les mouvements de chacun des nœuds et permet donc d'écrire l'éq. (18) comme un ensemble d'équations différentielles indépendantes du type (17), une pour chaque sommet.

De plus, nous avons choisi un schéma d'intégration explicite, dans lequel les forces élastiques sont estimées à l'instant t pour calculer la nouvelle position du sommet à l'instant t+1 (éq. (19)) :

Cette intégration explicite permet une remise à jour du modèle avec une fréquence compatible avec les besoins de la simulation. Une approche fondée sur un schéma semi-implicite permettrait un comportement plus rigide, mais demanderait la résolution d'un système linéaire à chaque itération, ce qui est trop coûteux en temps de calcul, compte tenu de la taille des maillages.

Un des gestes de base de la simulation de chirurgie consiste à découper des tissus mous. Avec notre modèle déformable, ce geste peut être effectué de manière très efficace. Nous simulons l'action d'un bistouri électrique sur un tissu mou en supprimant successivement tous les tétraèdres qui sont en contact avec l'extrémité de l'instrument. Lorsque l'on supprime un tétraèdre, environ 280 soustractions de réels sont effectuées pour retrancher la contribution du tétraèdre aux informations de rigidité des sommets, arêtes et triangles voisins. En mettant à jour localement les informations de rigidité, le tissu garde exactement les mêmes propriétés que si nous avions retiré le tétraèdre correspondant sur le maillage au repos. Grâce à la continuité volumique de la modélisation par les éléments finis, la déformation du tissu reste réaliste durant la découpe.

Les tissus vivants, qui sont composés essentiellement d'eau, peuvent être considérés comme quasiment incompressibles. Cette propriété est particulièrement difficile à modéliser et conduit la plupart du temps à des problèmes d'instabilité. C'est le cas avec le modèle de St Venant–Kirchhoff, car le matériau simulé atteint un comportement incompressible lorsque la constante de Lamé λ tend vers l'infini. Mais prendre une valeur trop grande pour λ nous forcerait à diminuer le pas de temps et donc à augmenter le temps de convergence. Une autre raison qui nous a poussé à ajouter une contrainte d'incompressibilité externe au modèle élastique est liée au modèle lui-même : le principal avantage du modèle de St Venant-Kirchhoff est d'utiliser le tenseur de déformation E, qui est invariant par rotation. Mais il est aussi invariant par symétrie, ce qui peut conduire au retournement de certains tétraèdres lorsqu'ils sont soumis à de fortes contraintes. Ce retournement sera rendu impossible par une contrainte géométrique de conservation du volume.

Nous avons donc choisi de pénaliser la variation de volume de chaque tétraèdre en appliquant, sur chacun de ses sommets, une force dirigée selon la normale à la face opposée. La norme de cette force étant proportionnelle au carré de la variation relative de volume et le sens étant donné par le signe de cette variation (Fig. 4) :

Ces forces agissent comme une augmentation de la pression à l'intérieur du tétraèdre. Cette méthode est relativement proche des multiplicateurs de Lagrange, qui sont souvent utilisés pour résoudre les problèmes de minimisation d'énergie sous contrainte.

Afin de mettre en évidence l'apport du modèle masse–tenseur «grands déplacements», nous avons renouvelé l'expérience proposée au § 2 sur le cylindre (Fig. 3). On peut ainsi apprécier l'avantage du modèle non linéaire, qui se déforme de manière beaucoup plus réaliste. En effet, on constate sur la Fig. 5a que les points peuvent maintenant suivre des trajectoires non rectilignes, permettant au cylindre de se courber naturellement. Si on compare ce résultat avec le modèle linéaire (Fig. 5b et c), on se rend compte de l'erreur commise par ce dernier dès que l'on sort de l'hypothèse des petits déplacements. Le modèle non linéaire permet au cylindre de se plier sans que sa longueur augmente ni que sa section se déforme. On peut aussi apprécier l'avantage qu'il présente en termes de conservation du volume.

La Fig. 6 présente la comparaison du modèle masse–tenseur non linéaire avec son homologue linéaire. Les modèles déformables ont été fixés au centre de leur face antérieure : cette zone correspond au hile hépatique, c'est-à-dire à l'endroit où les différents arbres vasculaires (artériel, veineux et biliaire) pénètrent dans le parenchyme hépatique. Nous utilisons un coefficient μ de l'ordre de ${10}^4-{10}^5\mathrm{kg}{\mathrm{cm}}^{-2},$ alors que les valeurs de λ sont issues du compromis entre un comportement le plus incompressible possible (λ très grand) et un temps de calcul correct. En général, cela nous mène à des valeurs de λ comprises entre $2\mu$ et $10\mu .$ Nous avons alors appliqué des forces sur le lobe droit du foie. Si on utilise le modèle linéaire, la partie droite du foie subit une importante (et irréaliste) augmentation de volume, alors qu'avec le modèle non linéaire, le lobe droit est capable de subir une rotation partielle, ce qui donne une impression de déformation beaucoup plus proche de la réalité.

Lorsque l'on ajoute la contrainte d'incompressibilité décrite au paragraphe précédent, on obtient des variations de volume encore plus faibles (voir Tableau 2). L'ajout de cette contrainte a aussi pour effet de stabiliser le comportement dans les zones de plus forte déformation.

Enfin, nous présentons un exemple de simulation d'un geste classique de chirurgie laparoscopique sur le foie. Le chirurgien opère à l'aide de deux instruments. Celui de droite est une pince avec laquelle il va saisir une zone du lobe droit du foie afin de le mettre sous tension, pendant qu'il va découper le parenchyme hépatique avec l'instrument de gauche. Au fur et à mesure de la découpe, le chirurgien écarte la partie du foie qu'il veut supprimer afin de faciliter la suite de la découpe. Cette partie peut alors subir un important mouvement de rotation qui sera fidèlement simulé par notre modèle non linéaire (Fig. 7).

Évidemment, le temps de calcul des forces pour ce modèle non linéaire est nettement plus long que dans le cas linéaire. D'abord, parce qu'il nous faut parcourir, en plus des sommets et des arêtes, tous les triangles et tous les tétraèdres du maillage. Ensuite, parce que les expressions des forces sont nettement plus compliquées. Après avoir factorisé au maximum les calculs, nous avons obtenu une fréquence de simulation cinq fois inférieure à celle du modèle linéaire, soit environ 9 Hz sur un pentium PIII 600 MHz pour un modèle de foie comportant 6 342 tétraèdres (ce qui correspond à 1 394 sommets), alors qu'avec le modèle linéaire nous atteignons la fréquence de 45 Hz. Ce temps de calcul permet néanmoins d'obtenir une fréquence de 30 Hz avec un maillage comportant environ 2 000 tétraèdres, ce qui est suffisant pour atteindre le temps réel visuel avec des objets relativement complexes et pour fournir un retour d'effort réaliste en utilisant la méthode d'extrapolation des forces décrite dans 〚9〛. Une architecture multiprocesseurs, disponible aujourd'hui à un prix raisonnable, permettrait déjà de traiter des maillages plus complexes ou d'atteindre une fréquence de simulation supérieure. En effet, en raison de l'utilisation d'un schéma d'intégration explicite, le calcul des forces internes et l'intégration du mouvement se font sommet par sommet, ce qui permet de paralléliser ces calcul en décomposant le maillage en sous-domaines.