Début 2017, je postais un article traitant du déplacement de particules sur un réseau hexagonal (ici). Au cœur des effets visuels obtenus sur canvas HTML5 : l’utilisation de forces physiques comme la viscosité et l’attraction d’un ressort.
L’effet d’amorti créé par la combinaison de ces deux forces spécifiques fut à l’origine constaté et mesuré dans les milieux fluides, alors pourquoi ne pas tenter franchement de reproduire les mouvements d’un fluide à partir d’un réseau de particules liées par ces forces ?
Cet article retrace un parcours de recherche d’effets fluides mis en œuvre avec le langage JavaScript et le <canvas> HTML5. En partant d’un système de particules mobiles liées par des forces simples (ressort, viscosité), la démarche s’oriente vers la discrétisation et la résolution des équations de Navier Stokes qui régissent les fluides dans les milieux continus. L’attention est portée sur la mise en place des systèmes de forces plutôt que sur les aspects graphiques.
Motivé comme jamais et convaincu de la pertinence de la démarche, je m’élançais.
1 – Le montage de la grille
Avant toute chose, il fallait dessiner un terrain de jeu pour mes particules. Devant la complexité inhérente à l’utilisation d’un réseau hexagonal tant dans le positionnement que le calcul des forces, j’optais pour une grille carrée.
Entre chaque particule : un ressort. De cette façon, chaque particule qui se situe pas aux frontières de la grille a 4 voisins qui concourent pour lui arracher les membres ou la repousser (pogo, ou, la danse des maillages). Les particules aux limites ont 2 voisins (sud/nord) ou (est/ouest), à l’exception des coins qui n’en ont qu’un.
Le positionnement peut s’écrire ainsi :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
App.setupGrid = function(h) { for (var x = 0.25 * this.width; x <= 0.75 * this.width; x += h) { var col = []; for (var y = 0.25 * this.height; y <= 0.75 * this.height; y += h) { var particle = { x: x, y: y, xLast: x, yLast: y, xSpeed: 0, ySpeed: 0, m: 250 } col.push(particle); } this.grid.push(col); } }; |
Voilà une belle grille remplie de carrés de côté h.
2 – Mise en tension de la grille
Il s’agit maintenant de parcourir cette grille pour appliquer à chacune des particules les forces de ses voisins, sans oublier de garder en mémoire les positions passées des voisins pour les calculs.
Voici en résumé la boucle de mouvement :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 |
App.move = function() { for (var i = 0; i < this.grid.length; i++) { for (var j = 0; j < this.grid[0].length; j++) { var particle = this.grid[i][j]; // Save previous position particle.xLast = particle.x; particle.yLast = particle.y; // Calculate forces from neighbours [...] var k = this.springStiffness, l = this.springLength; if (i > 0) {// Western neighbour particleWest = this.grid[i-1][j]; var toWest = { x: particleWest.xLast - particle.x, y: particleWest.yLast - particle.y }; var distWest = Math.sqrt(Math.pow(toWest.x, 2) + Math.pow(toWest.y, 2)); var FwestAmp = distWest - l;// Force amplitude var FwestExtraAngle = (FwestAmp < 0 ? Math.PI : 0); var FwestAngle = segmentAngleRad(particle.x, particle.y, particleWest.xLast, particleWest.yLast, true) + FwestExtraAngle; var FwestX = FwestAmp * Math.cos(FwestAngle); var FwestY = FwestAmp * Math.sin(FwestAngle); } if (i < this.grid.length - 1) {// Eastern neighbour particleEast = this.grid[i+1][j]; [...] } if (j > 0) { particleNorth = this.grid[i][j-1]; [...] } if (j < this.grid[0].length - 1) { particleSouth = this.grid[i][j+1]; [...] } // Total forces var Fx = k * ( (!!particleWest ? FwestX : 0) + (!!particleEast ? FeastX : 0) + (!!particleNorth ? FnorthX : 0) + (!!particleSouth ? FsouthX : 0) ); var Fy = k * ( (!!particleWest ? FwestY : 0) + (!!particleEast ? FeastY : 0) + (!!particleNorth ? FnorthY : 0) + (!!particleSouth ? FsouthY : 0) ); // Now get acceleration var xAcc = Fx / particle.m; var yAcc = Fy / particle.m; // Simple discrete integration particle.xSpeed += xAcc; particle.ySpeed += yAcc; particle.x += 0.1 * particle.xSpeed; particle.y += 0.1 * particle.ySpeed; } } }; |
On remarque que l’on soustrait ici une longueur l à l’amplitude (ici l’amplitude est homogène à une distance ultérieurement agrémentée de la raideur du ressort, uniforme sur la grille) de la force : c’est la longueur à vide du ressort. Au lancement, la grille devrait donc passer d’un ensemble de carrés de côté h (positionnement initial) à un état d’oscillation jusqu’à ce que ses carrés soient stabilisés et de longueur l.
Mais là, il y eut un imprévu :
SimplGrid 1.6 by Alexandre Andrieux – SimplX
1
Bien sûr ! La grille s’effondre sur elle même et continue dans un état instable dit explosif car les valeurs en jeu sont trop importantes et rien ne les retient, rien ne dissipe l’énergie initiale du système. J’introduisis alors la viscosité, qui vint « ralentir » l’accélération en lui imputant une partie de la vitesse :
|
1 2 |
var xAcc = Fx / particle.m - this.viscosity * particle.xSpeed; var yAcc = Fy / particle.m - this.viscosity * particle.ySpeed; |
En calibrant cette viscosité puis la raideur et la taille des ressorts, ça fonctionnait !
SimplGrid 1.8 by Alexandre Andrieux – SimplX
1
J’y ajoutai finalement une possibilité de cliquer pour attraper le coin de la grille.
Quel émerveillement ! La propagation des forces sur toute la grille. Une grille stable.
(Note : cette technique est utilisée dans les jeux vidéos pour reproduire les mouvements complexes des tissus)
3 – Sous l’eau
Fort de l’obtention d’une grille qui répond calmement aux efforts extérieurs, j’entrepris de préparer un ensemble de contraintes qui rapprocheraient le système d’une vraie piscine 2D.
La surface la plus intéressante ici pour l’étude du fluide est le dessus ; c’est dans un cas réel l’interface avec l’air qui fait office de fracture entre deux milieux et nous permet de visualiser les mouvements surfaciques de l’eau, les vagues. Dans notre cas bidimensionnel, il s’agit donc d’appliquer une condition aux limites de stabilité pour les côtés gauche, droit, et le fond, pour laisser libre de mouvement la ligne supérieure.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
App.setupGrid = function(h) { [...] // Stabilize some particles (fixed contraints) var leftWall = this.grid[0], rightWall = this.grid[this.grid.length - 1], bottomWall = []; for (var c = 0; c < this.grid.length; c++) { // Add last particle on the column (bottom) bottomWall.push(this.grid[c][this.grid[c].length - 1]); } // Entire wall var walls = _.union(leftWall, rightWall, bottomWall); _(walls).each(function(particle) { particle.fixed = true; }); }; |
Les fonctions utilitaires (du namespace « _ ») font partie de la librairie Lodash, très pratique. Ces particules se voient ainsi infliger un booléen qui, dans la boucle de mouvement, sera détecté et les empêchera de subir les agitations de leurs voisins.
J’ajoutai à cela une rotation du socle basée sur la position du curseur, afin d’influer sur les conditions aux limites en direct.
SimplGrid 1.18 by Alexandre Andrieux – SimplX
1
En désactivant la fonction de dessin des noeuds :
SimplGrid 1.19 by Alexandre Andrieux – SimplX
1
Somme toute, ces effets liquides pourraient être utilisés dans une application web.
En réduisant l’espacement des particules et avec masque hexagonal en HTML/CSS, j’obtins ce résultat :
SimplGrid 1.23 by Alexandre Andrieux – SimplX
1
Ici seules 4 lignes de particules sont laissées libres pour des questions de performance. Cela apporte aussi une certaine stabilité, les particules ayant moins de facilité à passer sous la ligne fixe Sud.
Cependant, cela ne les empêche pas de franchir ces lignes. Dans le cas des lignes Est et Ouest, le réalisme visuel en est altéré. Il fallait essayer une autre approche, un angle d’attaque complètement différent. Peut-être que…
4 – La discrétisation des équations de Navier-Stokes
Depuis le temps qu’il me démangeait d’exploiter les équations de Navier-Stokes ! Ce sont elles qui régissent les fluides dans les milieux continus. Pas de temps à perdre, finies les particules. Cette fois-ci il s’agira d’un maillage, c’est-à-dire un ensemble de nœuds immobiles qui servent de support de calcul aux équations.
L’équation de Navier-Stokes sous sa forme tensorielle s’écrit :
Où ρ est la densité, u est le champ de vitesse (flux), ∇ représente l’opérateur nabla (qui servira surtout sous forme de Laplacien), p est la pression, I est la matrice identité, τ est le tenseur des contraintes (d’ordre 2 ici), g représente toutes les forces d’accélération continues (inertie, gravité, électrostatique…) et ⊗ est le produit vectoriel.
Note : en notation scientifique tout ce qui n’est pas en gras est une grandeur scalaire (tenseur d’ordre 0), comme un nombre. Le reste sont des tenseurs d’ordre 1 (vecteurs) au moins, sinon plus (2 pour les matrices).
L’équation sous cette forme englobant un grand nombre de cas, je choisis de limiter mes calculs au cas où le fluide est Newtonien, incompressible, et subit un unique champ de forces externe conservatif : la gravité. Une telle simplification peut me priver des effets dissipatifs de type vortex, de natures complexes. Je garde toutefois la mécanique de diffusion telle qu’on la retrouve dans l’équation de la chaleur, qui représente la dissipation d’énergie par conduction dans un corps homogène. On retrouvera ce phénomène dans le résultat final.
Quelques heures de calculs (souvent faux, héhé) me permirent de transformer ces équations continues en équations discrètes. J’obtins les relations qui lient un nœud du maillage à ses 4 voisins. Au passage je supprimai la gravité qui avait peu d’intérêt en deux dimensions dans un cas simple.
Voici la transcription de ces équations en JavaScript :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
// Get new Vx and Vy based on current pressure field (might be controlled by user click or something) // also based on neighbours's speeds. // Careful, we're using old speeds for North and West as they've already been computed if (!node.isWall && !isNaN(node.Vx)) { var Vx = node.Vx, Vy = node.Vy, VxN = nodeNorth.VxLast || Vx, VyN = nodeNorth.VyLast || Vy, VxW = nodeWest.VxLast || Vx, VyW = nodeWest.VyLast || Vy, VxS = nodeSouth.Vx || Vx, VyS = nodeSouth.Vy || Vy, VxE = nodeEast.Vx || Vx, VyE = nodeEast.Vy || Vy, P = node.P, Psouth = nodeSouth.P || P, Peast = nodeEast.P || P, h = this.h, visc = this.viscosity, rho = this.rho, tau = this.tau; // I have discretized the incompressible Navier-Stokes equation with no gravity or anythin', juste // some gradients, derivatives and laplacians var Ax = VxN + VxS + VxE + VxW - 4 * Vx, Bx = Vx * (VxE - Vx) + Vy * (VxS - Vx), Cx = Peast - P; var deltaVx = tau * ( visc/(h*h) * Ax - 1/h * Bx - 1/(rho*h) * Cx ); var Ay = VyN + VyS + VyE + VyW - 4 * Vy, By = Vx * (VyE - Vy) + Vy * (VyS - Vy), Cy = Psouth - P; var deltaVy = tau * ( visc/(h*h) * Ay - 1/h * By - 1/(rho*h) * Cy ); // Update speed vector! (cautiously) node.Vx += Math.max(-1000, Math.min(1000, deltaVx)); node.Vy += Math.max(-1000, Math.min(1000, deltaVy)); |
S’il faut retenir une chose de ces tribulations, c’est qu’on va pouvoir ici calculer la vitesse à partir de la pression. Ces équations créent un lien direct entre le champ de pression tel qu’appliqué à chaque nœud, et la vitesse du nœud (valeur vectorielle calculée sur le nœud, mais rappelons qu’il ne bouge pas) en fonction de ses voisins. La variable tau ici (qui fait office de pas de temps, non pas le tenseur des contraintes évoqué précédemment) est un facteur fort de régulation du système. C’est cette variable qui va en grande partie dicter la vitesse de propagation des variations de pression.
5 – L’utilisateur met la pression
Mais ces variations de pression, d’où viennent-elles ?
C’est l’utilisateur, par sa toute puissance et son divin clic, qui va localement augmenter la pression et en voir les effets sur tout le maillage ! En somme, toutes les valeurs de la pression p vont être augmentées en fonction de leur distance au curseur lors du clic, selon une loi exponentielle décroissante.
|
1 2 3 4 5 6 7 8 |
// Get pressure from user action // If close enough to target, increase P (thus add some spatial gradient, biatch!) if (App.target.x) { var distToTarget = Math.sqrt(Math.pow(node.x - App.target.x, 2) + Math.pow(node.y - App.target.y, 2)); node.P += 0.1 * this.P0 * Math.exp(-distToTarget / 100); } else { node.P = node.initialP; } |
Après de longs ajustements de la densité, de la viscosité et de l’amplitude de pression appliquée afin d’obtenir, là aussi, un système stable :
Naviering the Stokes 1.4 by Alexandre Andrieux – SimplX
1
Ça marche ! Miracle. J’assistai enfin à la diffusion des variations de pression sur le champ de vitesse.
En cliquant sur le rectangle supérieur, on visualise la modification en intensité du champ de vitesse. Elle est rendue visible par le biais du modèle de couleurs HSL, en particulier sa première dimension (hue) qui crée cet effet arc-en-ciel. Le rectangle inférieur représente la pression appliquée sur le point cliqué et ses alentours.
En cliquant de manière prolongée, on peut toutefois arriver à faire sortir le système de son état stable. Il se met alors à clignoter, signe d’une augmentation erratique de l’intensité des vecteurs vitesse. Épileptiques, s’abstenir !
En modifiant un peu les paramètres de rendu :
Naviering the Stokes 1.6 by Alexandre Andrieux – SimplX
1
6 – Résulats finaux
Si maintenant l’on associe l’intensité du champ de vitesse au rayon des nœuds dessinés et à leur luminosité, on obtient un effet organique de retour au calme, d’absorption de la perturbation que l’on crée :
Naviering the Stokes 1.7 by Alexandre Andrieux – SimplX
1
Ma satisfaction fut complète.
HTML5 Icosahedron Smiles only by Alexandre Andrieux – SimplX
1
C’est ainsi que, sous la bannière SimplX , je pus enfin voir en action les équations simplifiées de Navier-Stokes dans mon navigateur – et à plus grande échelle franchir le cap entre mes travaux sur les particules (le monde discret) et le monde du continu.
Merci qui ? Merci…
Naviering the Stokes – SimplX by Alexandre Andrieux – SimplX
1
Note : Ce dernier pen a eu l’honneur d’être partagé à la communauté par la plateforme Codepen où il a été créé. Tout comme ce pendule double élastique basé sur le système de particules détaillé dans cet article 🙂