Animacija fluida u okruženju s preprekama

Programsko rješenje se sastoji od dvije temeljne klase (uz par pomoćnih). Klasa Solver je što se simulacije fluida tiče temelj programa. Sve metode koje se bave simulacijom ponašanja fluida, računanjem tlakova i brzina, iterativnim postupkom, provjeravanjem konvergencije i slično smještene su isključivo na ovom mjestu. Klasa Visualizer služi samo za vizualizaciju rezultata dobivenih korištenjem podataka i metoda u ovoj klasi. Ovdje je pojašnjenje metoda klase Solver.

Iteriranje glavne petlje algoritma

Jedina metoda koju je potrebno pozvati nakon što su zadane vrijednosti svih potrebnih parametara je advanceTimeStep(). Ova metoda vraća logičku vrijednost true u slučaju da je korak simulacije završio konvergencijom ili false ako do konvergencije nije došlo. Sama metoda samo poziva druge metode u slijedu koji u potpunosti prati algoritam SIMPLE. Izvorni kod je dan u nastavku.

      public bool advanceTimeStep ()
      {
            switchVelocities ();
           
            for (int iteracija = 0; iteracija < maxIteracija; iteracija++)
            {
                  calculateNewVelocities ();
                  calculatePCorrection ();
                  applyPCorrection ();
                  if (maxPCorrection < convergenceCondition)
                  {
                        applyVelocityCorrections ();
                        return true;
                  }
            }
      return false;
      }

U kodu su vidljivi koraci SIMPLE algoritma. Unutar for petlje se prvo računaju novi iznosi brzina (drugi korak algoritma), zatim se iteracijskim postupkom koji je u potpunosti unutar metode calculatePCorrection() računaju korekcije tlaka koje se primjenjuju na stare vrijednosti tlaka da bi se dobile nove. Nakon toga se provjerava konvergencija rezultata. Ako je konvergencija postignuta, petlja se prekida i u pozivajući program se vraća vrijednost koja označava uspješno završen korak. Ako do konvergencije nije došlo, cijeli se postupak nastavlja. Petlja će se izvršiti najviše maxIteracija puta. Ovaj je broj postavljen dovoljno velik da se može smatrati da ako rezultat nije dostignut u tolikom broju pokušaja da neće ni biti dostignut (rješenje je divergiralo), ali ipak dovoljno nizak da program ne upadne u dugotrajnu besmislenu blokadu. Ako nakon maxIteracija konvergencija nije dostignuta, petlja se prekida i u pozivajući program se vraća odgovarajuća vrijednost.

Proračun vrijednosti brzina za novi vremenski korak

Računanje novih brzina, kako je vidljivo u metodi advanceTimeStep(), vrši metoda calculateNewVelocities(). Kako ova metoda ima vrlo jednostavan kod koji samo za svaki element mreže poziva metodu calculateOneElementNewVelocity(), bitno je obratiti pažnju na ovu drugu metodu. Njen pseudokod se može prikazati otprilike kako slijedi:

     calculateOneElementNewVelocity (i, j)
          postavi vrijednosti karakterističnih susjednih veličina
          za svaki susjedni element koji je prepreka
               prilagodi karakteristične susjedne veličine
          izračunaj nove brzine prema diskretnom izrazu momentne jednadžbe
     kraj metode

Pri tome su karakteristične susjedne veličine sve one vrijednosti tlaka i brzina koje su potrebne za rješavanje momentne jednadžbe. Unutar ove metode te veličine imaju imena dana prema stranama svijeta kako je prikazano na slici. Ovakav način zapisa je razumljiviji od korištenja indeksa polja, pogotovo uzme li se u obzir kompleksnost momentne jednadžbe u diskretnom obliku.

Proračun korekcija tlaka

Sljedeća metoda koja se poziva u advanceTimeStep() je  metoda calculatePCorrection(). Njen pseudokod se može prikazati kako slijedi:

    bool calculatePCorrection()
        za svaki element mreže
            postavi vrijednost korekcije tlaka na nula
            postavi vrijednost stare korekcije tlaka na nula
        ponavljaj do konvergencije ili najviše do maxIteracija puta
            postavi najveću korekciju tlaka na nula
            postavi najveću razliku korekcija tlaka na nula
            za svaki element mreže
                izračunaj korekciju tlaka
                ako je razlika nove i stare korekcije tekućeg elementa veća
                od dosad najveće razlike korekcije tlaka
                    postavi najveću razliku korekcije tlaka na razliku nove
                    i stare korekcije tekućeg elementa
                ako je korekcija tlaka tekućeg elementa veća od dosad
                najveće korekcije tlaka
                    postavi najveću korekciju tlaka na korekciju tlaka
                    tekućeg elementa
            ako je najveća razlika korekcije tlaka manja od uvjeta
            konvergencije
                korekcije tlaka su dostigle konačnu vrijednost; izađi iz
                metode
            inače
                za svaki element mreže
                    postavi staru korekciju na vrijednost nove korekcije
    kraj metode

Ova je metoda posebno bitna zbog iterativnog postupka koji dosad nije objašnjen. Radi se o jednostavnom postupku. Nove vrijednosti korekcija tlaka se izračunavaju sve dok najveća razlika korekcije između dvije uzastopne iteracije ne padne na prihvatljivu razinu. Unutar cijelog postupka, koji ne ovisi o ovome koraku, vrlo je bitno da se pamti najveća korekcija tlaka na mreži u apsolutnom iznosu. Ta je veličina, naime, mjerilo konvergencije algoritma SIMPLE.

Ova metoda za svaki element mreže poziva zasebnu metodu koja zapravo rješava Poissonovu jednadžbu. Njen je naziv calculateOneElementNewPressure(), a pseudokod se može prikazati na sljedeći način:

    calculateOneElementNewPressure(i, j)
        postavi vrijednosti karakterističnih susjednih veličina
        za svaki susjedni element koji je prepreka
            prilagodi karakteristične susjedne veličine
        izračunaj nove korekcije prema diskretnom izrazu Poissonove jednadžbe
    kraj metode

Karakteristične susjedne veličine se označavaju jednako kao kod rješavanja momentne jednadžbe, s tim da ovdje postoje dvije dodatne varijable za tlak: westPressure i southPressure.

Ostale metode

U klasi Solver postoje još neke dodatne pomoćne metode. Tako metoda switchVelocities() prije ulaska u glavnu petlju algoritma SIMPLE zapisuje tekuće vrijednosti brzina u polje starih brzina. Da bi se postupak ubrzao, umjesto prepisivanja svih vrijednosti jednostavno se zamijene samo imena polja, tj. zamijene se reference na ove objekte. Metoda applyPCorrection() vrši jednostavno zbrajanje tekućih vrijednosti tlaka i njihovih korekcija kako je predviđeno u četvrtom koraku algoritma SIMPLE. Slično, ali s brzinama, radi i applyVelocityCorrections(). Također, postoje metode za postavljanje različitih veličina na njihove početne vrijednosti.