3D Wärmeleitungsgleichung mit expliziter FDM

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: 07.12.2022, 13:29 Titel: 3D Wärmeleitungsgleichung mit expliziter FDM

Hallo,

ich bin gerade dabei die 3D Wärmeleitungsgleichung für einen Stahlwürfel in Matlab zu programmieren.
Die Randbedigungen sind konstant und sind auf die aussenflächen bezogen.
Zum Zeitpunkt 0s sind dalle Knoten im inneren des Würfels bei 100 K und die Aussenflächen haben eine konstante Temperatur.

Hier ist der verwendete Code:

Code:

clc;
clear all;
clf;
format short g

fontSize = 18;
tv=1;

% Materialeigenschaften
rho=7850; % kg/m3
cp=500; % J/kgK
lambda=50; % W/mK
alpha=lambda/(rho*cp);% m2/s

% Domänengröße
L=0.1;
B=0.1;
H=0.1;

nx=31;
ny=31;
nz=31;
dx=L/(nx-1);
dy=B/(ny-1);
dz=H/(nz-1);

nt=600; % Anzahl der Zeitschritte
t=600; % Gesamte Berechnungszeit
dt=t/nt; % Zeitschrittweite

% CFL Bedingung
CFL=alpha*(dt/(dx^2) + dt/(dy^2) + dt/(dz^2));

%
Tn=zeros(nx,ny,nz);
x=linspace(0,L,nx);
y=linspace(0,H,ny);
z=linspace(0,B,nz);
[X,Y,Z]=meshgrid(x,y,z);

K=zeros(nx,ny,nz);
K([1 end],:,:)=alpha;
K(:,[1 end],:)=alpha;
K(:,:,[1 end])=alpha;

% Initiale Bedingungen
T=zeros(nx,ny,nz);
T(:,:,:)=100;
t=0;

% Randbedigungen
Tu=600; % unten
To=900; % oben
Tl=400; % links
Tr=800; % rechts
Tv=300; % vorn
Th=300; % hinten

% Seitenflächen
Tn(1,1:nx-1,1:nz-1)=Tv; % vorn
Tn(ny,1:nx-1,1:nz-1)=Th; % hinten
Tn(1:ny-1,1:nx-1,1)=Tu; % unten
Tn(1:ny-1,1:nx-1,nz)=To; % oben
Tn(1:ny-1,1,1:nz-1)=Tl; % links
Tn(1:ny-1,nx,1:nz-1)=Tr; % rechts

% Innere Knoten
Tn(2:ny-1,2:nx-1,2:nz-1)=T(2:ny-1,2:nx-1,2:nz-1);
Tc=Tn;

ct=1;
cn=1;
% iteration der Inneren Temperaturwerte des Würfels
for n=1:nt
t=n*dt;
for i=2:nx-1
for j=2:ny-1
for k=2:nz-1
Tn(i,j,k)=Tc(i,j,k)+dt*(K(i,j,k)/rho/cp)*((Tc(i+1,j,k)-2*Tc(i,j,k)+Tc(i-1,j,k))/dx/dx+(Tc(i,j+1,k)-2*Tc(i,j,k) +Tc(i,j-1,k))/dy/dy+(Tc(i,j,k+1)-2*Tc(i,j,k)+Tc(i,j,k-1))/dz/dz);
end
end
end
Tc=Tn;
Tbar=mean(Tn(:));
cn=cn+1;
% Visualisierung
if (mod(t,dt)==0)
clf;
%
f = figure(1);
f.Position(3:4) = [1080 840];
%XY
subplot(2,2,1)
slice(X,Y,Z,Tn,[],[],[0,H],'cubic');
set(gca,'FontSize',fontSize)
colormap(jet)
xlabel('X -Axis','FontSize', fontSize)
ylabel('Y -Axis','FontSize', fontSize)
zlabel('Z -Axis','FontSize', fontSize)
shading(gca,'interp')
axis([0 L 0 B 0 H]);
view(30,30);
% XZ
subplot(2,2,2)
slice(X,Y,Z,Tn,[],[0,B],[],'cubic');
set(gca,'FontSize',fontSize)
colormap(jet)
%set ( gca, 'ydir', 'reverse' )
xlabel('X -Axis','FontSize', fontSize)
ylabel('Y -Axis','FontSize', fontSize)
zlabel('Z -Axis','FontSize', fontSize)
shading(gca,'interp')
axis([0 L 0 B 0 H]);
view(30,30);
% YZ
subplot(2,2,3)
slice(X,Y,Z,Tn,[0,L],[],[],'cubic');
set(gca,'FontSize',fontSize)
colormap(jet)
%set ( gca, 'ydir', 'reverse' )
xlabel('X -Axis','FontSize', fontSize)
ylabel('Y -Axis','FontSize', fontSize)
zlabel('Z -Axis','FontSize', fontSize)
shading(gca,'interp')
axis([0 L 0 B 0 H]);
view(30,30);
% 3D
subplot(2,2,4)
slice(X,Y,Z,Tn,[0,L/2,L],[B/2,B],[0],'linear');
set(gca,'FontSize',fontSize)
set(gca,'YDir','normal')
colormap(jet)
set ( gca, 'xdir', 'reverse' )
set ( gca, 'ydir', 'reverse' )
xlabel('X -Axis','FontSize', fontSize)
ylabel('Y -Axis','FontSize', fontSize)
zlabel('Z -Axis','FontSize', fontSize)
axis([0 L 0 B 0 H])
view(-150,30);
hold on
shading(gca,'interp')
p = get(subplot(2,2,4),'Position');
cb=colorbar('Position', [0.93 0.25 0.025 0.6]);
set(cb,'FontSize',fontSize);
title(sprintf('xy t = %.2f',t),'FontSize', fontSize);
sgtitle(sprintf('Durchschnittstemperatur = %.1f Zeit = %.2f minutes',Tbar,t))
caxis([0, 900]);
pause(0.1);
end
end

Funktion ohne Link?

Leider ändert sich die Temperaturverteilung im Würfel überhaupt nicht.
Ich habe irgendwo einen Fehler gemacht, kann diesen aber nicht finden.

VG
Steffen

end.png

Beschreibung:

letzter Zeitschritt

Download

Dateiname:

end.png

Dateigröße:

93.09 KB

Heruntergeladen:

140 mal

start.png

Beschreibung:

erster Zeitschritt

Download

Dateiname:

start.png

Dateigröße:

91.06 KB

Heruntergeladen:

145 mal

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 07.12.2022, 13:47 Titel:

Hallo,

wo zumindest ein Fehler ist, kann ich dir sagen. Wenn man die Zeile in der innersten for-Schleife so umschreibt:

Code:

dT = dt*(K(i,j,k)/rho/cp)*((Tc(i+1,j,k)-2*Tc(i,j,k)+Tc(i-1,j,k))/dx/dx+(Tc(i,j+1,k)-2*Tc(i,j,k) +Tc(i,j-1,k))/dy/dy+(Tc(i,j,k+1)-2*Tc(i,j,k)+Tc(i,j,k-1))/dz/dz);
if dT ~= 0
dT
end
Tn(i,j,k)=Tc(i,j,k)+dT;

Funktion ohne Link?

dann erfolgt keine Ausgabe, d.h. dT ist immer 0. Warum das so ist, müsstest du nun herausfinden.

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

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 07.12.2022, 13:50 Titel:

Oh... nein, musst du nicht.

Die Werte von K, die du in der besagten Zeile verwendest, sind die im Inneren des Würfels. Die Nicht-Nullwerte von K sind aber am Rand des Würfels.

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: 08.12.2022, 12:49 Titel:

Hallo Harald,

danke für den Hinweis.
Jetzt läuft alles wie geplant.

VG
Steffen


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.