Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Bücher:

Modellbildung und Simulation Mechatronischer Systeme

Fachkräfte:
weitere Angebote

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

2D Wärmeleitungsgleichung in Matlab

 

SteffenB
Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 03.11.22
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 03.11.2022, 11:06     Titel: 2D Wärmeleitungsgleichung in Matlab
  Antworten mit Zitat      
Hallo,

ich versuche die 2D Wärmeleitungsgleichung mit Orts und Zeitabhängigen Randbedingungen mittels implizitem SOR Algorithmus zu lösen.
Die Umsetzung für konstante Werte an den Rändern ist simpel. Bei der Implementierung der Raum und Zeitabhängigkeit habe ich einige Probleme.

Der Code für die konstanten BC's ist folgender:

Code:
%********Transiente Wärmeleitungsgleichung implizite SOR Methode**********
clc
clear all
close all

%Input Parameter
nx=100;            % Anzahl der Knoten in x
ny=100;            % Anzahl der Knoten in y
nt=50;

L=2;               % Länge der Plate
H=0.75;            % Höhe der Plate

x=linspace(0,L,nx);% Knotenvektor in x
y=linspace(0,H,ny);% Knotenvektor in y

dx=x(2)-x(1);      % Abstand zwischen zwei Knoten in x
dy=y(2)-y(1);      % Abstand zwischen zwei Knoten in x

%Fehlerkriterium
Fehler=9e9;        % Fehler
Toleranz=1e-4;     % Toleranz
dt=1e-3;           % Zeitschrittweite
omega=1.1;         % Überrelaxationsparameter

% Randbedigungen an den Kanten des Rechtecks
TL=400;            % Linke Seite
TR=800;            % Rechte Seite
TT=600;            % Oben
TB=900;            % Unten

% Initiale Temperatur
Tini=300;
T=Tini*ones(nx,ny);
T(2:ny-1,1)=TL;
T(2:ny-1,nx)=TR;
T(1,2:nx-1)=TT;
T(ny,2:ny-1)=TB;

% Temperaturen an den Eckpunkten
T(1,1)=(TL+TT)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;

% Aktualisierung der Temperaturwerte
Tini=T;
Talt=T;
Tfinal=T;

% Berechnung der
rho=1;
Lambda=1;
cp=1;
alpha=Lambda/(rho*cp);
k1=alpha*dt/(dx)^2;
k2=alpha*dt/(dy)^2;
m1=1/(1+2*k1+2*k2);
m2=k1*m1;
m3=k2*m1;

ct=1;  
gs=1;
for k=1:nt
    Fehler=9e9;
    time(k)=k*dt;
    while(Fehler>Toleranz)
        for i=2:nx-1
            for j=2:ny-1
                T(i,j)=((Tini(i,j)*m1))+(m2*(T(i-1,j)+Talt(i+1,j)))+(m3*(T(i,j+1)+Talt(i,j-1)));
                Tfinal(i,j)=((1-omega)*Talt(i,j))+(T(i,j)*omega);
            end
        end
        Fehler=max(max(abs(Talt-T)));
        Talt=Tfinal;
        gs=gs+1;
    end
    Tini=Tfinal;
    % Ergebnissplot des Temperaturfeldes
    figure(1)

    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca,'ydir','reverse')

    xlabel('X -Axis')
    ylabel('Y -Axis')
    title(sprintf('Anzahl der Iterationen Implizite SOR Methode= %d',gs));
    pause(0.00000000000000003)
    M(ct)=getframe(gcf);
    ct=ct+1;

end
 


Ich versuche die Temperatur auf den Rändern zum einen als variable Werte in Abhängigkeit von der Position auf der jeweiligen Kante einzubinden (Links, Oben und Rechts). Unten soll ein Raum und Zeitabhängiger Wärmestrom angesetzt werden.

Die räumliche Variabilität wollte ich erstmal in Form von simplen linearen Funktionen einbinden.

Habe mir dazu folgenden Ansatz überlegt:

Code:
for k=1:nt      
for i=1:nx
    for j=1:ny
        if j==1                 % Knoten der Kante unten
            time=k*dt;
            x=i*dx;
            TB=ab*x+bb;
            T(i,j)=TB;
        elseif j==ny         % Knoten der Kante oben
            time=k*dt;
            x=i*dx;
            TT=at*x+bt
            T(i,j)=TT;
        else
        end
    end
end
for j=1:ny
    for i=1:nx
        if i==1                  % Knoten der Kante Links
            time=k*dt;
            y=j*dy;
            TL=aL*y+bL;
            T(i,j)=TL;
        elseif i==nx          % Knoten der Kante Rechts
            time=k*dt;
            y=j*dy;
            TR=aR*y+bR;
            T(i,j)=TR;
        else
        end
    end
end
T(1,1)=(TL+TT)/2;
T(nx,ny)=(TR+TB)/2;   % Update der Temperatur
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;
end


Ich habe es leider noch nicht geschafft meinen Ansatz um ein variablen Wärmestrom unten zu erweitern und in den funktionierenden Code für die konstanten Randbedingungen einzubinden.

Vieleicht hat ja jemand einen Hinweis bzgl. des Problems

VG
Steffen

Schema.png
 Beschreibung:

Download
 Dateiname:  Schema.png
 Dateigröße:  20.24 KB
 Heruntergeladen:  32 mal
Private Nachricht senden Benutzer-Profile anzeigen


Harald
Forum-Meister

Forum-Meister


Beiträge: 24.265
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 03.11.2022, 11:45     Titel:
  Antworten mit Zitat      
Hallo,

reicht es nicht aus, die Variable time in der Berechnung von T in der gewünschten Form zu verwenden?

Generell würde ich zu Vektorisierung raten. Das macht den Code kompakter, lesbarer und beschleunigt in der Regel die Berechnungen. Für die Aktualisierung der Temperaturen in der while-Schleife kannst du filter2 verwenden.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
SteffenB
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 03.11.22
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 09.11.2022, 14:39     Titel:
  Antworten mit Zitat      
Hallo Harald,

vielen Dank für deine Antwort.
Wenn ich die loops für die Randknoten in den Loop implementiere dann rechnet Matlab ohne einen Ende zu finden. Irgendwo habe ich einen Fehler eingebaut.
Die Schleifen laufen nicht wie geplant.

Code:
%********Transiente Wärmeleitungsgleichung implizite SOR Methode**********
clc
clear all
close all

%Input Parameter
nx=100;            % Anzahl der Knoten in x
ny=100;            % Anzahl der Knoten in y
nt=50;

L=2;               % Länge der Plate
H=0.75;            % Höhe der Plate

x=linspace(0,L,nx);% Knotenvektor in x
y=linspace(0,H,ny);% Knotenvektor in y

dx=x(2)-x(1);      % Abstand zwischen zwei Knoten in x
dy=y(2)-y(1);      % Abstand zwischen zwei Knoten in x

%Fehlerkriterium
Fehler=9e9;        % Fehler
Toleranz=1e-4;     % Toleranz
dt=1e-3;           % Zeitschrittweite
omega=1.1;         % Überrelaxationsparameter

% Randbedigungen an den Kanten des Rechtecks
TL=400;            % Linke Seite
TR=800;            % Rechte Seite
TT=600;            % Oben
TB=900;            % Unten

aB=1;
aT=1;
aL=1;
aR=1;
bT=1;
bB=TB;
bT=TT;
bL=TL;
bR=TR;

% Initiale Temperatur
Tini=300;
T=Tini*ones(nx,ny);
T(2:ny-1,1)=TL;
T(2:ny-1,nx)=TR;
T(1,2:nx-1)=TT;
T(ny,2:ny-1)=TB;

% Temperaturen an den Eckpunkten
T(1,1)=(TL+TT)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;

% Aktualisierung der Temperaturwerte
Tini=T;
Talt=T;
Tfinal=T;

% Berechnung der
rho=1;
Lambda=1;
cp=1;
alpha=Lambda/(rho*cp);
k1=alpha*dt/(dx)^2;
k2=alpha*dt/(dy)^2;
m1=1/(1+2*k1+2*k2);
m2=k1*m1;
m3=k2*m1;

ct=1;  
gs=1;
for k=1:nt
    Fehler=9e9;
    time(k)=k*dt;
    while(Fehler>Toleranz)
        for i=1:nx
            for j=1:ny
                if j==1
                    time=k*dt;
                    x=i*dx;
                    TB=aB*x+bB;
                    T(i,j)=TB;
                elseif j==ny
                    time=k*dt;
                    x=i*dx;
                    TT=aT*x+bT;
                    T(i,j)=TT;
                else
                end
            end
        end
        for j=1:ny
            for i=1:nx
                if i==1
                    time=k*dt;
                    y=j*dy;
                    TL=aL*y+bL;
                    T(i,j)=TL;
                elseif i==nx
                    time=k*dt;
                    y=j*dy;
                    TR=aR*y+bR;
                    T(i,j)=TR;
                else
                end
            end
        end
        T(1,1)=(TL+TT)/2;
        T(nx,ny)=(TR+TB)/2;
        T(1,ny)=(TT+TR)/2;
        T(nx,1)=(TL+TB)/2;
        for i=2:nx-1
            for j=2:ny-1
                T(i,j)=((Tini(i,j)*m1))+(m2*(T(i-1,j)+Talt(i+1,j)))+(m3*(T(i,j+1)+Talt(i,j-1)));
                Tfinal(i,j)=((1-omega)*Talt(i,j))+(T(i,j)*omega);
            end
        end
        Fehler=max(max(abs(Talt-T)));
        Talt=Tfinal;
        gs=gs+1;
    end
    Tini=Tfinal;
    % Ergebnissplot des Temperaturfeldes
    figure(1)

    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca,'ydir','reverse')

    xlabel('X -Axis')
    ylabel('Y -Axis')
    title(sprintf('Anzahl der Iterationen Implizite SOR Methode= %d',gs));
    pause(0.00000000000000003)
    M(ct)=getframe(gcf);
    ct=ct+1;

end
time'


Viele Grüße
Steffen
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.265
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 09.11.2022, 20:48     Titel:
  Antworten mit Zitat      
Hallo,

die Initialisierung würde ich nicht in die while-Schleife packen.
Bei so länglichem Code kann ich dir nur empfehlen, das zu machen, was ich auch machen würde: debuggen und schauen, ob das Problem liegt. Und natürlich wie gesagt: vektorisieren. Dann ist das ganze lesbarer.

Das einzige, was mir auffällt: sollte
T(i,j)=((Tini(i,j)*m1)) ... nicht
T(i,j)=((Talt(i,j)*m1)) .. sein?

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
SteffenB
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 03.11.22
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 15.11.2022, 14:18     Titel:
  Antworten mit Zitat      
Hallo Harald,

ich habe den Part mit der Erstellung der Temperaturwerte an den Rändern separat getestet.
Er enthält Randbedingungen die einen linear ansteigenden Verlauf aufweisen.
Der Code ist folgender:

Code:
% Randbedigungen an den Kanten des Rechtecks
nx=6;              % Anzahl der Knoten in x
ny=6;              % Anzahl der Knoten in y
nt=2;
dt=1e-3;           % Zeitschrittweite

L=0.01;            % Länge der Plate
H=0.01;            % Höhe der Plate

x=linspace(0,L,nx);% Knotenvektor in x
y=linspace(0,H,ny);% Knotenvektor in y

dx=x(2)-x(1);      % Abstand zwischen zwei Knoten in x
dy=y(2)-y(1);      % Abstand zwischen zwei Knoten in x

% Initiale Temperatur
Tini=00;
T=Tini*ones(nx,ny);
T0=T;

TL=75;             % Linke Seite
TR=115;            % Rechte Seite
TT=100;            % Oben
TB=100;            % Unten

% Initialisierung an den Rändern
for j=1:ny
    if j==1
        for i=2:nx-1
            T(i,1)=TL;
        end
    elseif j==ny
        for i=2:nx-1
            T(i,ny)=TR;
        end
    else
    end
end
for i=1:nx
    if i==1
        for j=2:ny-1
            T(1,j)=TT;
        end
    elseif i==nx
        for j=2:ny-1
            T(nx,j)=TB;
        end
    else
    end
end
T1=T;
% Temperaturen an den Eckpunkten
T(1,1)=(TL+TT)/2;   % Links oben
T(nx,ny)=(TR+TB)/2; % Rechts unten
T(1,ny)=(TT+TR)/2;  % Rechts oben
T(nx,1)=(TL+TB)/2;  % Links unten

aR=1e04;
bR=TR;
aL=2e4;
bL=TL;

aT=1e04;
bT=TT;
aB=1e04;
bB=TB;

% RB in Form von räumlichen Temperaturfunktionen an den Rändern
% Links und Rechts
for j=1:ny
    if j==1
        for i=2:nx-1
            yy=i*dy;
            TL=aL*yy+bL;
            T(i,1)=TL;
        end
    elseif j==ny
        for i=2:nx-1
            yy=i*dy;
            TR=aR*yy+bR;
            T(i,ny)=TR;
        end
    else
    end
end
T2=T;

% Oben und Unten
for i=1:nx
    if i==1
        for j=2:ny-1
            xx=j*dx;
            TT=aT*xx+bT;
            T(1,j)=TT;
        end
    elseif i==ny
        for j=2:ny-1
            xx=j*dx;
            TB=aB*xx+bB;
            T(nx,j)=TB;
        end
    else
    end
end
% Temperaturen an den Eckpunkten
T(1,1)=(TL+TT)/2;   % Links oben
T(nx,ny)=(TR+TB)/2; % Rechts unten
T(1,ny)=(TT+TR)/2;  % Rechts oben
T(nx,1)=(TL+TB)/2;  % Links unten

T3=T;
figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')

T0
T1
T2
T3


[img][/img]

Das funktioniert soweit ganz gut (Siehe Bild "Randbedigungen als f(x bzw. y)".

Wenn ich jetzt die Loops einbinde dann scheint etwas noch nicht zu funktionieren.
Ich habe die Methode abgeändert und das implizite Crank-Nicholson Verfahren verwendet.

Code:
%********Transiente Wärmeleitungsgleichung implizite C-N Methode**********
clc
clear all
close all
solver=1;

%Input Parameter
nx=100;            % Anzahl der Knoten in x
ny=100;            % Anzahl der Knoten in y
t=1;

L=1;               % Länge der Plate
H=1;               % Höhe der Plate

x=linspace(0,L,nx);% Knotenvektor in x
y=linspace(0,H,ny);% Knotenvektor in y

dx=L/(nx-1);      % Abstand zwischen zwei Knoten in x
dy=L/(ny-1);      % Abstand zwischen zwei Knoten in x
dt=0.001;         % Zeitschrittweite
nt=t/dt;

%Fehlerkriterium
Fehler=9e9;        % Fehler
Toleranz=1e-4;     % Toleranz

% Randbedigungen an den Kanten des Rechtecks
TL=400;            % Linke Seite
TR=800;            % Rechte Seite
TT=600;            % Oben
TB=900;            % Unten

% Randbedigungen vereinfacht, nur Raumabhängige für den Anfabg
aR=0;
bR=TR;

aL=0;
bL=TL;

aT=0;
bT=TT;

aB=0;
bB=TB;

% Initiale Temperatur
Tini=300;
T=Tini*ones(nx,ny);
T(2:ny-1,1)=TL;
T(2:ny-1,nx)=TR;
T(1,2:nx-1)=TT;
T(ny,2:ny-1)=TB;
Tsor=zeros(nx,ny);

% Temperaturen an den Eckpunkten
T(1,1)=(TL+TT)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;

% Aktualisierung der Temperaturwerte
Talt=T;
Tprev = Talt;

% Berechnung der
rho=1;
Lambda=1;
cp=1;
alpha=Lambda/(rho*cp);

k1=alpha*dt/((dx)^2);
k2=alpha*dt/((dy)^2);
m1=(1-k1-k2)/(1+k1+k2);
m2=k2/(2*(1+k1+k2));
m3=k1/(2*(1+k1+k2));

ct=1;  
cn=1;
time=zeros(nt);

for k=1:nt
    Fehler=9e9;
    time(k)=k*dt;
    while(Fehler>Toleranz)
        % Temperaturen an den Randpunkten
        % RB in Form von räumlichen Temperaturfunktionen an den Rändern
        % Links und Rechts
        for j=1:ny
            if j==1
                for i=2:nx-1
                    yy=i*dy;
                    TL=aL*yy+bL;
                    T(i,1)=TL;
                end
            elseif j==ny
                for i=2:nx-1
                    yy=i*dy;
                    TR=aR*yy+bR;
                    T(i,ny)=TR;
                end
            else
            end
        end
        % Oben und Unten
        for i=1:nx
            if i==1
                for j=2:ny-1
                    xx=j*dx;
                    TT=aT*xx+bT;
                    T(1,j)=TT;
                end
            elseif i==ny
                for j=2:ny-1
                    xx=j*dx;
                    TB=aB*xx+bB;
                    T(nx,j)=TB;
                end
            else
            end
        end
        % Temperaturen an den Eckpunkten
        T(1,1)=(TL+TT)/2;   % Links oben
        T(nx,ny)=(TR+TB)/2; % Rechts unten
        T(1,ny)=(TT+TR)/2;  % Rechts oben
        T(nx,1)=(TL+TB)/2;  % Links unten
        % Aktualisierung der Temperaturwerte
        Talt=T;
        Tprev = Talt;
        % Temperaturen in den inneren Knoten
        for i=2:nx-1
            for j=2:ny-1
                V=(Tprev(i+1,j)+Tprev(i-1,j)+Tprev(i,j+1)+Tprev(i,j-1));
                H=(Talt(i+1,j)+Talt(i-1,j)+Talt(i,j+1)+Talt(i,j-1));
                T(i,j)=(Tprev(i,j)*m1)+(m3*H)+(m2*V);
            end
        end
        Fehler=max(max(abs(Talt-T)));
        Talt=T;
        cn=cn+1;
    end
    Tprev=T;
    Fehler=9e9;
    Rechenzeit=toc;
    % Ergebnissplot des Temperaturfeldes
    figure(1)

    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca,'ydir','reverse')

    xlabel('X -Axis')
    ylabel('Y -Axis')
    title(sprintf('Anzahl der Iterationen Implizite CN Methode= %d',cn));
    pause(0.00000000000000003)
    M(ct)=getframe(gcf);
    ct=ct+1;
end


Irgendwas scheint bei der Übergabe nicht richtig zu funktionieren bzw. ich glaube das die Schleifen nicht so durchlaufen werden wie gewünscht.

VG
Steffen

Resultat CN.png
 Beschreibung:

Download
 Dateiname:  Resultat CN.png
 Dateigröße:  19.49 KB
 Heruntergeladen:  20 mal
Randbedingungen als f(x bzw. y).png
 Beschreibung:

Download
 Dateiname:  Randbedingungen als f(x bzw. y).png
 Dateigröße:  31.02 KB
 Heruntergeladen:  29 mal
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.265
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 16.11.2022, 20:40     Titel:
  Antworten mit Zitat      
Hallo,

ich kann mich bei meinen Empfehlungen nur wiederholen:
1. Debugggen um den Fehler zu finden
2. Vektorisieren. Code wird lesbarer, übersichtlicher, und vielleicht erledigt sich damit auch der Fehler.
3. Überlegen, ob die Randbedingungen wirklich in der while-Schleife gesetzt werden sollen.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
.





 Impressum  | Nutzungsbedingungen  | Datenschutz  | Werbung/Mediadaten | Studentenversion | FAQ | goMatlab RSS Button RSS


Copyright © 2007 - 2023 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks
Partner: LabVIEWforum.de

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.