Probleme bei der Implementierung der Lattice-Boltzmann-Methode

Ich möchte eine grundlegende 2D-Implementierung der Lattice Boltzmann-Methode (in Javascript) erstellen, um den Gas-/Flüssigkeitsfluss zu simulieren, aber ich bin auf einige Probleme gestoßen und kann die Ursache dafür anscheinend nicht finden.

Das Problem ist, dass beim Ausführen des Kollisionsschritts (ich habe Streaming noch nicht vollständig implementiert) die Simulation explodiert und ins Unendliche abläuft und ich meinen Code sehr gründlich auf Fehler überprüft habe.

Ich beginne mit der Initialisierung meiner Variablen:

function Grid(width, height){
        this.width = width;
        this.height = height;
        this.f = [];
        this.e = [
            new vec2(0,0),
            new vec2(0,1),
            new vec2(1,0),
            new vec2(0,-1),
            new vec2(-1,0),
            new vec2(1,1),
            new vec2(1,-1),
            new vec2(-1,-1),
            new vec2(-1,1)
        ];
        this.w  = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36];
        for(var x=0; x<width; x++){
            this.f[x] = [];
            for(var y=0; y<height; y++){
                this.f[x][y] = [];
                for(var i=0; i<9; i++){
                    this.f[x][y][i] = 20+Math.random();
                }
            }
        }
    }

Hier ist this.f die Verteilungsfunktion mit 9 Richtungen: 0 für keine Bewegung, 1-4 für gerade Bewegung und 5-9 für die schrägen Geschwindigkeiten. Die Verteilungsfunktion hat Werte im Bereich von 20-21. This.e enthält alle Richtungen als Vektoren (ich habe meine eigene einfache vec2-Klasse erstellt) und this.w enthält die Gewichte für jede Richtung. Die erste Frage, die ich stellen möchte, lautet: Habe ich die Distribution mit den richtigen Werten initialisiert? Ich konnte nicht viele Informationen über die Initialisierungswerte finden.

Zweitens könnten meine Fehler im Kollisionsschritt liegen:

    collide: function(){
        for(var x=0; x<width; x++){
                for(var y=0; y<height; y++){
                    var rho = 0;
                    for(var i=0; i<9; i++){
                        rho += this.f[x][y][i];
                    }

                    var u = new vec2(0,0);
                    for(var i=0; i<9; i++){
                        u  = u.add( this.e[i].scale( this.f[x][y][i] ) );
                    }
                    u = u.scale( 1/rho );
                    for(var i=0; i<9; i++){
                        var eu = u.dot(this.e[i]);
                        var uu = u.dot(u);
                        var fiEq = this.w[i]*rho*(1+3*eu+4.5*eu*eu-1.5*uu);
                        this.f[x][y][i] = this.f[x][y][i] - .005*(this.f[x][y][i]-fiEq);
                    }
                }
            }
        }

Zuerst berechne ich die makroskopische Dichte rho, dann berechne ich mit Vektoraddition die makroskopische Geschwindigkeit u. Ich verwende diese in der folgenden Formel, um die Gleichgewichtsverteilung zu berechnen: und dann entspanne ich die Verteilung in Richtung Gleichgewicht. Als zweite Frage möchte ich fragen, ob ich den Kollisionsschritt richtig implementiert habe.

Ein Live-Beispiel ist hier zu sehen

Wenn mehr Informationen benötigt werden, um zu helfen, lassen Sie es mich wissen und wenn Sie mir helfen könnten, würde ich es wirklich schätzen!