Fluide stabile

Stabil-fluidelor este o metodă de simulare fluid ( NMF ) , descris de Jos Stam în 1999 ca parte a conferinței SIGGRAPH99.

Folosește o grilă euleriană  : un spațiu compus exclusiv din fluid în care mișcarea fluidului este descrisă de un câmp de vectori de viteză (numit și viteză ).

Formularea matematică

Ecuațiile Navier-Stokes

Modelul Stable-Fluid utilizează ecuațiile Navier-Stokes în cazul fluidelor incompresibile.

• Ecuația de incompresibilitate ∇→⋅v→=0{\ displaystyle {\ overrightarrow {\ nabla}} \ cdot {\ vec {v}} = 0}
• Ecuația echilibrului impulsului ∂v→∂t=-(v→⋅∇→)v→-1ρ∇→p+ν∇2v→+f→{\ displaystyle {\ frac {\ partial {\ vec {v}}} {\ partial t}} = - \ left ({\ vec {v}} \ cdot {\ overrightarrow {\ nabla}} \ right) {\ vec {v}} - {\ frac {1} {\ rho}} {\ overrightarrow {\ nabla}} p + \ nu \ nabla ^ {2} {\ vec {v}} + {\ vec {f}} }

Aceste ecuații pot fi interpretate după cum urmează:

• Termenul reprezintă vectorul vitezei într-un punct dat al fluidului.v→{\ displaystyle {\ vec {vb}}}
• Termenii accelerațiilor corespund: unde al doilea membru se numește advecție .
  ∂v→∂t+(v→⋅∇→)v→{\ displaystyle {\ frac {\ partial {\ vec {v}}} {\ partial t}} + \ left ({\ vec {v}} \ cdot {\ overrightarrow {\ nabla}} \ right) {\ vec {vb}}}
  
• Forțele de vâscozitate corespund termenului:
  ν∇2v→{\ displaystyle \ nu \ nabla ^ {2} {\ vec {vb}}}
• Forțele de compresie corespund termenului:
  1ρ∇→p{\ displaystyle {\ frac {1} {\ rho}} {\ overrightarrow {\ nabla}} p}
• La forțele externe corespund termenului .f→{\ displaystyle {\ vec {f}}}
• Incompresibilitatea are ca rezultat relația:
  ∇→⋅tu=0{\ displaystyle {\ overrightarrow {\ nabla}} \ cdot u = 0}

Teorema descompunerii Helmholtz-Hodge

De acolo, încercăm să simplificăm ecuațiile Navier și Stokes utilizând teorema descompunerii Helmholtz-Hodge . În acest fel includem termenul de proiecție . P{\ displaystyle \ mathbb {P}}

P∂v→∂t=P(-(v→⋅∇→)v→-1ρ∇→p+ν∇2v→+f→){\ displaystyle \ mathbb {P} {\ frac {\ partial {\ vec {v}}} {\ partial t}} = \ mathbb {P} \ left (- \ left ({\ vec {v}} \ cdot {\ overrightarrow {\ nabla}} \ right) {\ vec {v}} - {\ frac {1} {\ rho}} {\ overrightarrow {\ nabla}} p + \ nu \ nabla ^ {2} {\ vec {v}} + {\ vec {f}} \ right)}

În plus, această teoremă ne permite să eliminăm mai multe elemente din ecuație, cum ar fi ecuația de incompresibilitate și termenul referitor la presiune:

P(∇→v)=0{\ displaystyle \ mathbb {P} ({\ overrightarrow {\ nabla}} v) = 0} P(∇→p)=0{\ displaystyle \ mathbb {P} ({\ overrightarrow {\ nabla}} p) = 0} P(v→)=v→{\ displaystyle \ mathbb {P} ({\ vec {v}}) = {\ vec {v}}}

Astfel ajungem la următoarea ecuație:

∂v→∂t=P(-(v→⋅∇→)v→+ν∇2v→+f→){\ displaystyle {\ frac {\ partial {\ vec {v}}} {\ partial t}} = \ mathbb {P} \ left (- \ left ({\ vec {v}} \ cdot {\ overrightarrow {\ nabla}} \ right) {\ vec {v}} + \ nu \ nabla ^ {2} {\ vec {v}} + {\ vec {f}} \ right)}

Modelul digital

În ceea ce privește rezoluția digitală, diferitele componente (difuzie, advecție, proiecție, forțe externe) nu sunt calculate și apoi adunate împreună, spre deosebire de forma sa matematică. Mai degrabă, considerăm că fiecare componentă, dintr-un câmp de vectori de intrare, efectuează modificări pe câmp și dă rezultatul componentei următoare.

Definim operatori pentru advecție ( ), difuzie ( ), aplicarea forțelor externe ( ) și proiecție ( ) și operator echivalent cu formularea matematică pentru o etapă de timp. LA{\ displaystyle \ mathbb {A}}D{\ displaystyle \ mathbb {D}}F{\ displaystyle \ mathbb {F}}P{\ displaystyle \ mathbb {P}}S{\ displaystyle \ mathbb {S}}

S(tu)=P∘F∘D∘LA(tu){\ displaystyle \ mathbb {S} (u) = \ mathbb {P} \ circ \ mathbb {F} \ circ \ mathbb {D} \ circ \ mathbb {A} (u)}

Astfel, la fiecare pas de timp, câmpul vectorial u trece succesiv prin operatorul de advecție apoi pe cele de difuzie, forțe externe și în cele din urmă prin operatorul de proiecție.

Implementarea metodei

Luați în considerare un caz 2D (mică diferență pentru 3D) și luați în considerare o rețea regulată de celule N × N, în care fiecare celulă este definită de un vector de viteză (în centrul celulei) și o presiune (pentru partea de proiecție).

Datele vor fi reprezentate în tabele în care i și j corespund poziției celulei curente din grilă. veu,j{\ displaystyle v ^ {i, j}}

Difuzie

Pasul de difuzie corespunde rezoluției vâscozității. Viteza fluidului va influența celulele vecine. Pentru aceasta, va fi utilizat un rezolvator liniar.

Dacă corespunde vitezei uneia dintre celule și , , și la celulele vecine, atunci Solver trebuie să rezolve sistemul în cazul în care pentru fiecare celulă: veu,j{\ displaystyle v ^ {i, j}}veu+1,j{\ displaystyle v ^ {i + 1, j}}veu-1,j{\ displaystyle v ^ {i-1, j}}veu,j+1{\ displaystyle v ^ {i, j + 1}}veu,j-1{\ displaystyle v ^ {i, j-1}}

(1-4LA)⋅veu,j→=LA⋅veu+1,j→+LA⋅veu-1,j→+LA⋅veu,j+1→+LA⋅veu,j-1→{\ displaystyle (1-4 \ mathrm {A}) \ cdot {\ vec {v ^ {i, j}}} = \ mathrm {A} \ cdot {\ vec {v ^ {i + 1, j}} } + \ mathrm {A} \ cdot {\ vec {v ^ {i-1, j}}} + \ mathrm {A} \ cdot {\ vec {v ^ {i, j + 1}}} + \ mathrm {A} \ cdot {\ vec {v ^ {i, j-1}}}}

sau

• A = d t ⋅ν⋅N × N
• d t  : pasul de timp.
• ν  : coeficient de difuzie / vâscozitate.
• N × N  : numărul de celule ⇒ corespunde 1 / Suprafață.

Pentru rezoluția sistemului, poate fi utilizat solverul Gauss-Seidel . În acest exemplu, lucrăm într-un tabel în care sunt folosiți doar vecinii direcți, celelalte celule nu apar și au un coeficient de 0.

Pseudo-cod, algoritm de difuzie int N = 100 // largeur et hauteur du tableau int kmax = 100 // nombre d'itérations à effectuer (solveur) float diff = 0.01 // coefficient de diffusion float u[N,N] // composante horizontale des vecteurs vélocité float v[N,N] // composante verticale des vecteurs vélocité float u0[N,N], v0[N,N] // tableaux temporaires relatifs à u et v // solveur linéaire Gauss-Seidel void linear_solver (int b, float * x, float * x0, float a, float c, int kmax ) { int i, j, k; for ( k=0 ; k<kmax ; k++ ) { FOR_EACH_CELL x[i,j] = (x0[i,j] + a*(x[i-1,j]+x[i+1,j]+x[i,j-1]+x[i,j+1]))/c; END_FOR set_bnd (b, x ); // conditions limites } } // opérateur de diffusion void diffuse_step() { SWAP ( u0, u ); SWAP ( v0, v ); float a=_fm->dt()*diff*N*N; lin_solve (1, u, u0, a, 1+4*a,kmax ); lin_solve (2, v, v0, a, 1+4*a,kmax ); }  

Advecție

În timpul etapei de advecție , deci considerăm că există o particulă la vectorul vitezei, efectuează un backtracking ( backtracking ) pentru a determina poziția particulei în no-time-ul anterior, în cele din urmă efectuăm o interpolare pentru a cunoaște viteza particule. Acesta din urmă va înlocui valoarea vectorului viteză. Prin urmare, încercăm să păstrăm energia conținută în particulă.

Pseudo-cod, algoritm de advecție advection (int b, float * d, float * d0, float * u, float * v) { int i, j, i0, j0, i1, j1; float x, y, s0, t0, s1, t1, dt0; dt0 = _fm->dt()*_n; //coefficient FOR_EACH_CELL // position de la particule au pas de temps précédent x = i-dt0*u[IX(i,j)]; y = j-dt0*v[IX(i,j)]; // condition limites, les bords de la grille if (x<0.5f) x=0.5f; if (x>_n+0.5f) x=_n+0.5f; i0=(int)x; i1=i0+1; if (y<0.5f) y=0.5f; if (y>_n+0.5f) y=_n+0.5f; j0=(int)y; j1=j0+1; // interpolation s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)]+t1*d0[IX(i0,j1)])+ s1*(t0*d0[IX(i1,j0)]+t1*d0[IX(i1,j1)]); END_FOR // conditions limites, bord de la grille set_bnd (b, d ); } advection_step() { SWAP ( u0, u ); SWAP ( v0, v ); diffuse_vel_x (1, u, u0, _visc,this->_algoVel_diffuse_k); diffuse_vel_y (2, v, v0, _visc,this->_algoVel_diffuse_k); }  

Proiecție

Proiecția face posibilă restabilirea conservării masei și, prin urmare, a incompresibilității fluidului. Pentru aceasta, trebuie să determinăm cantitățile de material de intrare și de ieșire la nivelul fiecărei celule, apoi să corectăm viteza astfel încât să avem cât mai mult material de intrare ca de ieșire.

Divergența este diferența dintre fluxurile de intrare și de ieșire. Pentru aceasta luăm componenta orizontală a vectorilor vecini situate orizontal și componenta verticală a vectorilor vecini situată vertical.

deuveuj=-0,5(vXeu+1,j-vXeu-1,j)⋅dy-0,5(vyeu,j+1-vyeu,j-1)⋅dX{\ displaystyle div ^ {ij} = - 0,5 (vx ^ {i + 1, j} -vx ^ {i-1, j}) \ cdot \ mathrm {d} y-0,5 (vy ^ {i, j + 1} -vy ^ {i, j-1}) \ cdot \ mathrm {d} x} vXeu,j{\ displaystyle vx ^ {i, j}} : componentă orizontală a vectorului viteză veu,j→{\ displaystyle {\ vec {v ^ {i, j}}}} vyeu,j{\ displaystyle vy ^ {i, j}} : componentă verticală a vectorului viteză veu,j→{\ displaystyle {\ vec {v ^ {i, j}}}} dX{\ displaystyle \ mathrm {d} x} : distanța dintre vectori șiveu+1,j{\ displaystyle v ^ {i + 1, j}}veu-1,j{\ displaystyle v ^ {i-1, j}} dy{\ displaystyle \ mathrm {d} y} : distanța dintre vectori șiveu,j+1{\ displaystyle v ^ {i, j + 1}}veu,j-1{\ displaystyle v ^ {i, j-1}}

În cazul unei rețele de celule N × N avem dX=dy=1NU{\ displaystyle \ mathrm {d} x = \ mathrm {d} y = {\ frac {1} {\ mathrm {N}}}}

Pseudo-cod, algoritm de proiecție project (float * u, float * v, float * p, float * div, int kmax ) { int i, j; //Calcul de la divergence FOR_EACH_CELL div[i,j] = -0.5f*(u[i+1,j]-u[i-1,j]+v[i,j+1]-v[i,j-1])/_n; p[i,j] = 0; END_FOR set_bnd (0, div ); set_bnd (0, p ); //Résolution du système lin_solve (0, p, div, 1, 4, kmax); //Mise à jour des vélocités FOR_EACH_CELL u[i,j] -= 0.5f*_n*(p[i+1,j]-p[i-1,j]); v[i,j] -= 0.5f*_n*(p[i,j+)]-p[i,j-1]); END_FOR set_bnd (1, u ); set_bnd (2, v ); }  

Condiții la limite

Extensii ale metodei

Grila-eșalonată

Într-o grilă eșalonată, variabilele scalare (cum ar fi presiunea) sunt stocate la centrele celulelor, în timp ce viteza este situată la marginile (sau fețele) celulei.

Folosirea unei rețele eșalonate evită decuplarea parului impar între presiune și viteză. Putem cunoaște exact volumul de fluid care intră și iese din celulă, ceea ce nu este cazul dacă vitezele sunt în centru. Deci, Grila Staggered permite calcule mai precise și mai puțină disipare digitală .

Dezavantajul este că valorile sunt stocate în locuri diferite, ceea ce face mai dificilă controlul volumului.

Înapoi pe drumul cel bun

Backtracking-ul ( backtracking-ul ) este metoda de găsire a poziției unei particule într-un eulerian fluid în noile timpuri anterioare. În advecția de bază, este pur și simplu o chestiune de a lua viteza, înmulțind cu pasul de timp și cu dimensiunile unei celule pentru a găsi poziția anterioară a particulei.

Acestea fiind spuse, fluidul urmează o cale curbiliniară și pur și simplu utilizarea vitezei nu este suficientă pentru a găsi cu precizie vechea poziție a particulei. Corecțiile pot fi făcute pentru a îmbunătăți precizia. Există Runge-Kutta , MacCormack și BFECC ( Compensarea și corectarea erorilor înapoi și înainte ) .

Prima figură a diagramei opuse arată cazul revenirii pe o pistă simplă: putem observa că revenirea pe cale recuperează o viteză (în gri) care nu aparține liniei caracteristice bune. MacCormack vă permite să vă apropiați de linia caracteristică (ultima figură din diagramă, vectorul în gri).

Flux pe o suprafață arbitrară și Distorsiune

În cazul în care se alege să lucreze pe un flux de fluid care traversează o suprafață, se regăsește în cazul în care celulele care formează sistemul nu sunt pătrate. Aceste deformări trebuie luate în considerare pentru a evita fenomenele nerealiste.

Suprafata libera

Noțiunea de suprafață liberă se referă la modelarea suprafeței apei. Există diferite metode, cum ar fi setul de nivel sau urmărirea rapidă și robustă a suprafeței fluidului.

Celule triunghiulare și tetraedre

În cazurile anterioare, celulele pătrate (sau cel puțin patrulaterele) au fost utilizate pentru carcasele 2D și cuburi pentru carcasa 3D . Cu toate acestea, este posibil să se utilizeze celule triunghiulare sau tetraedre.

Vorticitate

Este vorba de înlocuirea etapei de proiecție cu o etapă în care se urmărește păstrarea vorticității . Vorticitatea corespunde fenomenelor vortexurilor. Deci, în loc să verificăm dacă o celulă conține cantitatea potrivită de fluid, verificăm incompresibilitatea prin fluxurile din jurul punctelor.

linkuri externe

• Simularea fluxului unui fluid - Humbert Florent - 2008
• Simulare rapidă a dinamicii fluidelor pe GPU - GPU Gems - NVIDIA
• Simularea fluidelor pe suprafețe - Ensiwiki
• Simulare de fluide pentru efecte vizuale - Magnus Wrenninge - 2003

Note

1. „Fluid stabil” - Jos Stam - 1999
2. Dinamica fluidelor în timp real pentru jocuri
3. "Calculul numeric al fluxului vâscoasă, incompresibil, de fluid, cu suprafață liberă" - Harlow, FH și Welch, JE - 1965,
4. „Transfer numeric de căldură” - SV Patankar
5. „O metodă MacCormack stabilă necondiționat”
6. „Utilizarea BFECC pentru simularea fluidelor”
7. „Fluxuri pe suprafețe de topologie arbitrară”
8. „metodă de setare a nivelului pentru interfața fluidă”
9. „Urmărirea rapidă și robustă a suprafețelor fluide”
10. „Simulare de lichid invizibil și incompresibil pe ochiuri triunghiulare”
11. "Stabile, circulație-conservare, fluide simple"