2D Wärmeleitungsgleichung in Matlab

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

SteffenB

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 03.11.22

	Wohnort: ---

	Version: ---

Verfasst am: 03.11.2022, 11:06 Titel: 2D Wärmeleitungsgleichung in Matlab

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

Funktion ohne Link?

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

Funktion ohne Link?

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:

149 mal

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 03.11.2022, 11:45 Titel:

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 ;)

SteffenB

Themenstarter

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 03.11.22

	Wohnort: ---

	Version: ---

Verfasst am: 09.11.2022, 14:39 Titel:

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'

Funktion ohne Link?

Viele Grüße
Steffen

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 09.11.2022, 20:48 Titel:

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 ;)

SteffenB

Themenstarter

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 03.11.22

	Wohnort: ---

	Version: ---

Verfasst am: 15.11.2022, 14:18 Titel:

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

Funktion ohne Link?

[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

Funktion ohne Link?

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:

144 mal

Randbedingungen als f(x bzw. y).png

Beschreibung:

Download

Dateiname:

Randbedingungen als f(x bzw. y).png

Dateigröße:

31.02 KB

Heruntergeladen:

161 mal

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 16.11.2022, 20:40 Titel:

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 ;)


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

RSS

Hosted by:

Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.