Finite Differenzen Wärmegleichung



  • Hat jemand einen Link zu einer guten Implementirung der Wärmeleitungsgleichung mittels finiter Differenzen in C++ parat?

    Ich dachte iegentlich, das wäre ein Standardproblem mit zig Varianten im Netz, finde aber wirklich nichts gutes.



  • Wenn möglich, wäre 2D oder 3D optimal.

    Vielen Dank fürs Lesen



  • Das ist immernoch sehr allgemein. Was soll denn der Code können? Welches Koordinatensystem (kartesisch, zylinder...)? Ist dein Rechengitter strukturiert? Welches Dateiformat hat es? Welche Zeitdiskretisierungen brauchst du? Implizit/Explizit? Welche Diskretisierung für den Diffusionsterm? Gibts nen Advektionsterm? Wovon ist der Quellterm abhängig? Welche Randbedingungen brauchst du?



  • Hey!

    Das ganze sol nicht praxisrelevant eingesetzt werden, sondern nur mir als Grundlage für einige Spielereien und Erweiterungen dienen.

    Ein normales rechteckiges Gitter oder Einheitsquadrat als Gebiet in kartesischen Koordianten sowie ein 5-Punktstern Laplace reichen mir vollkommen aus. Advektionsterm muss nicht sein, Randbedingung kann beliebig sein.

    Mir reicht also das einfachste vom einfachen.



  • Wenn man einige Vereinfachungen annimmt, ist es schnell gemacht. Annahmen:

    • Der Diffusionskoeffizient ist konstant
    • Das Gitter ist äquidistant und kartesisch
    • Kein Quellterm

    Den Code habe ich gerade aus dem Kopf geschrieben. Hoffe, es haben sich keine Fehler eingeschlichen.

    #include <vector>
    #include <string>
    #include <iostream>
    using namespace std;
    
    typedef vector<vector<double> > grid; // ein 2D Array (mehr oder weniger)
    
    // Funktion, um Temperaturen in Datei zu schreiben
    void write_grid(const grid& g, const std::string& fname)
    {
        ofstream write(fname.c_str());
        for (int i = 0; i != g.size(); ++i)
        {
            for (int j = 0; j != g[i].size(); ++j)
                write << g[i][j] << ' ';
            write << '\n';
        }
    }
    
    int main()
    {
        const int X = 20; // Dimensionen des Gitters
        const int Y = 20;
    
        const double dx = 0.02;  // Abstand der Gitterpunkte zueinander in m. Das Rechengebiet ist damit 4x4 cm gross. dx = dy weil aequidistant
        const double dt = 0.1;   // Zeitschritt in s
        const double a = 115e-6; // Waermeleitzahl m2/s
    
        const double CFL = a*dt / (dx*dx); // Berechne CFL Kriterium
        cout << "CFL Zahl = " << CFL;
        if (CFL > 0.5) // Fehler. Verfahren nicht stabil
            return -1;
    
        // Rechengitter erstellen
        grid T(X, vector<double>(Y, 300.0)); // ein 20x20 Rechengitter, das an allen Gitterpunkten eine Temperatur von 300 K
    
        T[10][10] = 10000.0; // in der Mitte ist es 10000 K heiss (Anfangsbedingung)
    
        write_grid(T, "start.txt");
    
        //Randbedingung: Raender links und rechts haben konstante Temperatur von 400 K. Da die Werte konstant bleiben, kann man es ausserhalb der Zeitschleife machen
        for (int x = 0; x != X; ++x)
        {
            T[x][0] = 400.0;
            T[x][Y - 1] = 400.;
        }
    
        grid T_neu = T; // T_neu ist die Temperatur im nachsten Zeitschritt
    
        for (double t = 0; t <= 30; t += dt) // simuliere 30 s lang
        {
            // laufe ueber alle Punkte des Rechengebietes ausser dem Rand
            for (int y = 1; y != Y - 1; ++y)
                for (int x = 1; x != X - 1; ++x)
                    // der 5-Punkte Stern. Euler vorwarts (explizit) in der Zeit, zentraler Differenzenquotient  (2. Ordnung) fuer Diffusionsterm
                    T_neu[x][y] = dt*a*( (T[x-1][y]-2.0*T[x][y]+T[x+1][y])/(dx*dx) + (T[x][y-1]-2.0*T[x][y]+T[x][y+1])/(dx*dx) ) + T[x][y];
    
            // Randbedingung oben und unten: adiabate Wand (Nullgradient in Wandnormalenrichtung erzwingen)
            for (int y = 1; y != Y - 1; ++y)
            {
                T_neu[0][y] = T_neu[1][y];
                T_neu[X - 1][y] = T_neu[X - 2][y];
            }
    
            T_neu.swap(T); // fuer den neachsten Zeitschritt wird T_neu zu T
        }
    
        write_grid(T, "end.txt");
    }
    

    Meinst du sowas?



  • Edit: Es ist die DFL Zahl, nicht CFL...


Log in to reply