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...