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

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

3D Wärmeleitungsgleichung mit expliziter FDM

 

SteffenB
Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 03.11.22
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 07.12.2022, 13:29     Titel: 3D Wärmeleitungsgleichung mit expliziter FDM
  Antworten mit Zitat      
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


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:  88 mal
start.png
 Beschreibung:
erster Zeitschritt

Download
 Dateiname:  start.png
 Dateigröße:  91.06 KB
 Heruntergeladen:  89 mal
Private Nachricht senden Benutzer-Profile anzeigen


Harald
Forum-Meister

Forum-Meister


Beiträge: 24.432
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 07.12.2022, 13:47     Titel:
  Antworten mit Zitat      
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;

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 ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.432
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 07.12.2022, 13:50     Titel:
  Antworten mit Zitat      
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 ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
SteffenB
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 6
Anmeldedatum: 03.11.22
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 08.12.2022, 12:49     Titel:
  Antworten mit Zitat      
Hallo Harald,

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

VG
Steffen
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 | FAQ | goMatlab RSS Button 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.