3

I've been writing a finite difference code for the simulation and detection of cracks using laser-induced thermography. The crack is implemented by factors a and b, which are "damping" the heat flow through the air filled crack using a ghost-point-approach. The 2D-Model runs as it is supposed to, stability condition is satisfied, all went fine. It is even well-proved with experimental data. Just copy and paste, and it will work.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       2D-Wärmeleitungsgleichung mit Ghost-Point-Methode und         %%
%%                       Finiter Differenzen                           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Leeren des Workspace und des Editors
clc;
clear all;
format long;
%%
% Abmessungen und Schrittweiten des Bleches im Raum
NX = 121;                            % Schrittzahl in x-Richtung  
NY = 121;                            % Schrittzahl in y-Richtung
XMAX = 30E-3;                        % Abmessung x-Richtung [m]
YMAX = 30E-3;                        % Abmessung y-Richtung [m]
dx = XMAX/(NX-1);                    % Schrittweite in x-Richtung [m]
dy = YMAX/(NY-1);                    % Schrittweite in y-Richtung [m]
x = -XMAX/2:dx:XMAX/2;               % Vektor mit x-Werten
y = -YMAX/2:dy:YMAX/2;               % Vektor mit y-Werten
% Laserparameter
P = 8325;                            % Laserleistung [W]
DIST = 10E-3;                        % Abtaststrecke [m]
SPOTD = 60E-6;                       % Spotdurchmesser [m]
ALPHA = 0.07;                        % Absorptionskoeffizient
% Schrittweiten in der Zeit                          
dt = 0.0002;                         % Zeitschritt [s]
NT = 400;                            % Anzahl der Zeitschritte
% Materialdaten Aluminium
DENS = 2700;                         % Dichte [kg*m^-3]
K_ALU = 180;                         % Wärmeleitfähigkeit Alu [W*(m*K)^-1]
C = 895;                             % spez. Wärmekapazität [J*K^-1 ]
k = K_ALU/(DENS*C);                  % Temperaturleitfähigkeit [m^2*s^-1]
% Materialdaten Luft im Riss
K_AIR = 0.025;                       % Wärmeleitfähigkeit Luft [W*(m*K)^-1]
% Variablen für die Ghost-Point-Methode
delta = 50E-6;                       % Breite Riss [m]
EPS = ((K_ALU)/(K_AIR)-1)*delta;     % Relation K_ALU, K_AIR, delta
a = (3*(EPS)+4*dx)/(EPS+dx);         % Faktor a
b = (dx)/(EPS+dx);                   % Faktor b
% Speicherallokation für die Temperatur-Matrix
T_OLD = zeros(NX,NY);                % Allokation alte Temperaturen
T_NEW = zeros(NX,NY);                % Allokation neue Temperaturen
% Speicherallokation für die Last-Matrix
q = zeros(NX,NY);                    % Allokation der Lasten  
%%
% Anfangsbedingung (Blechtemperatur)
for i=1:NX
    for j=1:NY
            T_OLD(i,j)=30;
    end
end
%%
% Instationärer Wärmestrom (Wärmestromdichte durch Line-Scan)
for i=1:NX
    for j=1:NY
        if ((i>=40) && (i<=80) && (j==61))
            q(i,j)=k*ALPHA*((P)/(DIST*SPOTD))/(K_ALU);  
        else
            q(i,j)=0;
        end
    end
end
%%
% Berechnung der Feldvariablen für jeden Zeitschritt
for it = 0:NT
    clf;                                 % Löscht aktuelle Figure
    T_NEW = T_OLD;                       % setze T_NEW als T_OLD 
    h=surf(x,y,T_OLD','EdgeColor','k');  % Plotting der Feldvariablen
    set(gca,'fontsize',20)
    colormap jet;                        % Farbschema der Farbskala
    colorbar('location','eastoutside'... % Position und Größe Farbschema
             ,'fontsize',20);
    shading interp                       % Interpolation zwichen Schritten
    axis ([-15E-3 15E-3 -15E-3 15E-3])   % Achsenskalierung 

    % Achsbeschriftungen

    title({['LST for crack detection using finite difference 2D Heat-'...
            'Diffusion'];['and ghost point method'] ;['time (\itt) = '...
            ,num2str(it*dt) 's']})

    xlabel('x in [m]','FontSize',20)
    ylabel('y in [m]','FontSize',20)
    zlabel('Temperatur in [°C]')  

    view(2);                             % Darstellung (1D, 2D, oder 3D)
    drawnow;                             % Aktualisiert die Figure
    pause(1E-40)                         % Pause zwischen einzelnen Figures
    refreshdata(h)                       % Aktualisiert die Daten in Figure

    % Explizites Finite-Differenzen-Verfahren (mittels zentralem DQ)

    for i=2:NX-1
    for j=2:NY-1
        if((j==69) && (i>=52) && (i<=68))
            T_OLD(i,j) = T_NEW(i,j)+(k*dt)/(dx^2)*(T_NEW(i+1,j)-...
                         a*T_NEW(i,j)+T_NEW(i-1,j)+b*T_NEW(i,j+1)+...
                         T_NEW(i,j-1))+dt*q(i,j);

        else
            T_OLD(i,j) = T_NEW(i,j)+(k*dt)/(dx^2)*(T_NEW(i+1,j)-...
                         4*T_NEW(i,j)+T_NEW(i-1,j)+T_NEW(i,j+1)+...
                         T_NEW(i,j-1))+dt*q(i,j);

        end
    end
    end
end
%% Programm Ende

But changing from 2D to 3D, the value of dt for stable behavior increases by orders of magnitude more than expected. I've tried everything. Using a simplier load, commenting out the "crack-loop", but nothing worked.

Calculating the stability-condition, 

dt <= dx^2/(6*k) = 1.39E-4   instead of 2E-10(!!!)

should be enough. But just try 2E-9, and the scheme will start oscillating already. The problem is, I need the thermal flow below the crack. This is why I need a 3D model, just in case you're asking. But this way it would take years to calculate just a few 10 to 100 milliseconds, which is the range I need.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       3D-Wärmeleitungsgleichung mit Ghost-Point-Methode und         %%
%%                       Finiter Differenzen                           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Leeren des Workspace und des Editors
clc;
close all;
format long;
%%
% Abmessungen und Schrittweiten des Bleches im Raum
NX = 121;                            % Schrittzahl in x-Richtung  
NY = 121;                            % Schrittzahl in y-Richtung
NZ = 9;                              % Schrittzahl in y-Richtung
XMAX = 30E-3;                        % Abmessung x-Richtung [m]
YMAX = 30E-3;                        % Abmessung y-Richtung [m]
ZMAX = 2E-3;                         % Abmessung z-Richtung [m]
dx = XMAX/(NX-1);                    % Schrittweite in x-Richtung [m]
dy = YMAX/(NY-1);                    % Schrittweite in y-Richtung [m]
dz = ZMAX/(NZ-1);                    % Schrittweite in z-Richtung [m]
x = 0:dx:XMAX;                       % Vektor mit x-Werten
y = 0:dy:YMAX;                       % Vektor mit y-Werten
z = 0:dz:ZMAX;                       % Vektor mit Z-Werten
% Schrittweiten in der Zeit                          
dt = 2E-10;                          % Zeitschritt [s]
NT = 5E11;                           % Anzahl der Zeitschritte
% Laserparameter
P = 8325;                            % Laserleistung [W]
DIST = 10E-3;                        % Abtaststrecke [m]
SPOTD = 60E-6;                       % Spotdurchmesser [m]
% Materialdaten Aluminium
DENS = 2700;                         % Dichte [kg*m^-3]
K_ALU = 180;                         % Wärmeleitfähigkeit Alu [W*(m*K)^-1]
C = 895;                             % spez. Wärmekapazität [J*K^-1 ]
k = K_ALU/(DENS*C);                  % Temperaturleitfähigkeit [m^2*s^-1]
ALPHA = 0.07;                        % Absorptionskoeffizient
% Materialdaten Luft im Riss
K_AIR = 0.025;                       % Wärmeleitfähigkeit Luft [W*(m*K)^-1]
% Variablen für die Ghost-Point-Methode
delta = 50E-6;                       % Breite Riss [m]
EPS = ((K_ALU)/(K_AIR)-1)*delta;     % Relation K_ALU, K_AIR, delta
a = (5*(EPS)+6*dx)/(EPS+dx);         % Faktor a
b = (dx)/(EPS+dx);                   % Faktor b
% Speicherallokation für die Temperatur-Matrix
T_OLD = zeros(NX,NY,NZ);             % Allokation alte Temperaturen
T_NEW = zeros(NX,NY,NZ);             % Allokation neue Temperaturen
T_AMB = 30;                          % Umgebungstemperatur
% Speicherallokation für die Last-Matrix
q = zeros(NX,NY,NZ);                 % Allokation der Lasten  
%%
% Anfangsbedingung (Blechtemperatur)
for i=1:NX
    for j=1:NY
        for k=1:NZ
            T_OLD(i,j,k)=T_AMB;
        end
    end
end
%%
% Instationärer Wärmestrom (Wärmestromdichte durch Line-Scan)
for i=1:NX
    for j=1:NY
        for k=1:NZ
            if ((j>=40) && (j<=80) && (i==60) && (k==9))
                q(i,j,k)=k*ALPHA*((P)/(DIST*SPOTD))/(K_ALU);  
            else
                q(i,j,k)=0;
            end
        end
    end
end
%%
% Berechnung der Feldvariablen für jeden Zeitschritt
for it = 0:NT
    clf;                                 % Löscht aktuelle Figure
    T_NEW = T_OLD;                       % setze T_NEW als T_OLD
    h = slice(x,y,z,T_OLD,...            % Plotting der Feldvariablen
      [],[],[2E-3]); 
    colormap jet;                        % Farbschema der Farbskala
    colorbar('location','eastoutside'... % Position und Größe Farbschema
             ,'fontsize',12);
    shading interp                       % Interpolation zwichen Schritten
    axis ([0 30E-3 0 30E-3 0 2E-3])      % Achsenskalierung
    % alpha(0.5);

    % Achsbeschriftungen
    title({['LST for crack detection using finite difference 3D Heat-'...
            'Diffusion'];['and ghost point method'] ;['time (\itt) = '...
            ,num2str(it*dt) 's']})
    xlabel('x in [m]')
    ylabel('y in [m]')
    zlabel('z in [m]')  

    view(2);                             % Darstellung (1D, 2D, oder 3D)
    drawnow;                             % Aktualisiert die Figure
    pause(1E-40)                         % Pause zwischen einzelnen Figures
    refreshdata(h)                       % Aktualisiert die Daten in Figure

    % Explizites Finite-Differenzen-Verfahren (mittels zentralem DQ)

    for i=2:NX-1
    for j=2:NY-1
    for k=1:NZ
        if((j>=45) && (j<=75) && (i==50) && (k<=9) && (k>=5))
            T_OLD(i,j,k) = T_NEW(i,j,k)+(k*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         a*T_NEW(i,j,k)+b*T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_AMB+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);      
        elseif(k==1)
            T_OLD(i,j,k) = T_NEW(i,j,k)+(k*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_NEW(i,j,k+1)+T_AMB)+...
                         dt*q(i,j,k);
       elseif(k==NZ)
            T_OLD(i,j,k) = T_NEW(i,j,k)+(k*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_AMB+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);
        else     
            T_OLD(i,j,k) = T_NEW(i,j,k)+(k*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_NEW(i,j,k+1)+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);

       end
    end
    end
    end
end
%% Programm Ende

Thank you in advance, im very desperated about this problem.

Greetings Alex

user6539314
  • 33
  • 5

2 Answers2

2

You have a bug in your code - in the 3D version, you introduce a looping variable called k for the z dimension. That variable overwrites your previously defined k coefficient. When fixed, it all works with dt = 1e-4 s in 3D. I just changed the k serving as a looping variable to kj. You can pick a better name. Actually, it is recommended to use a little longer names for looping variables, not just i, j, k... - like 2 or three letters instead one.

atru
  • 4,699
  • 2
  • 18
  • 19
  • Another queation just came up. As the load is applied as a transient heat flow and not as a Neumann-B.C, I have to handle the boundaries. As you can see in my script, I assumed the grid points for differenciation out of the domain to be ambient temperature for every time increment. Actually this doesn't work, since the points at k=NZ have to heat up, too, which nearly doesn't happen. Doing so in 2D again is no Problem, because there's n Gradient in z-direction. – user6539314 Aug 10 '17 at 13:12
  • Heh, learned it pretty much the same way ;) it's a fun bug :) Thinking about Stack rules, you may have to post this question as a separate question. I don't mind, but chances are some people will disagree with it being in this post as an answer. Maybe you can keep it and move it later. I'll try answering soon! :) – atru Aug 10 '17 at 14:49
  • Cant' comment on your answer since my reputation is still low! In the physical situation, what are your actual conditions in z? Ambient temperature or insulated (no-flux) boundary? If it's no-flux (Neumann) then similarly to what you wrote, you should use T_NEW(i,j,k+1) = T_NEW(i,j,k) i.e. the boundary node is equal to the previous one in z. For the bottom boundary you can do the same thing, just different order: T_NEW(i,j,1) = T_NEW(i,j,2). – atru Aug 10 '17 at 16:13
  • Exactly, there's no heat flux on the surface (k=NZ), except the transient heat flux caused by the laser beam conducting into the x,y and z direction. No convection or radiation is considered. What worries me is the difference between the 2D and the 3D heat distribution using same parameters. Maybe it is because the heat capacity of a 3D model is much higher, resulting in a lower heating, because the bulk material has to be heated. This in fact results in raising up the laser power to scale the plots so that they match. – user6539314 Aug 10 '17 at 22:37
  • I didn't check your scheme rigorously, but I've used a similar scheme years ago, looks ok from afar. It seems logical that you would need more heat input from the heat source in a 3D domain as it is larger. Except if your heat source spanned the entire domain height and you would have no flux or periodic boundaries on top and bottom. I see it spans only one node in z and the bottom boundary is ambient. One thing that comes to mind - have you done the grid independence study? Also, no flux boundary is ok if that portion of the system is insulated in real life. – atru Aug 11 '17 at 01:51
  • I just had to look what grid Independence study means, but reading about it, it seems it is just important handling fluids? I thought a numerial problem can't be grid-independent in general, because a coarser mesh affects the accuracy of the solution. So in my mind, every difference in grid spacing results in a slightly different solution. Do you have experiences modelling heat sources? Actually, I have to apply the whole intensity (which refers to an area) of the laser on just one grid line (which is no area, obviously), as the laser diameter is just about 90 micrometer. – user6539314 Aug 16 '17 at 11:39
  • It's good to do it and it's important for any numerical modeling, other name is "grid refinement study", for instance. With refinement your solution will be converging towards the correct solution. You can then choose which grid you want to work with based on your accuracy needs and computational load/efficiency restrictions. In simple problems you can actually reach that. Within reasonable accuracy you can really reach the correct solution. In more complicated problems you not rarely have to leave the whole idea altogether and just work with one grid hoping that it's close enough. – atru Aug 16 '17 at 16:40
  • It's not a must but in a small problem like this it would be beneficial to check it. Just run a couple grids with spacing ideally different by a factor of 2 and plot the solution as a function of spacing. Solution can be mean temperature in the system or difference of interpolated temperature fields (you'll need to interpolate because you'll be comparing matrices with different number of elements). You can also use your finest grid as the "correct" solution and look at the error relative to it. That's a common practice too. – atru Aug 16 '17 at 16:44
0

Another queation just came up. As the load is applied as a transient heat flow and not as a Neumann-B.C, I have to handle the boundaries. As you can see in my script, I assumed the grid points for differenciation out of the domain to be ambient temperature for every time increment. Actually this doesn't work, since the points at k=NZ have to heat up, too, which nearly doesn't happen. Doing so in 2D again is no problem, because there's no gradient in z-direction. So do you have some advice how to modify my code? I thought about substituting T_AMB by T_NEW(i,j,k), so that T_NEW(i,j,k+1) equals T_NEW(i,j,k). which gives a reasonable plot. But again, I don't know if this is mathematically correct. Below, theres a slightly corrected code regarding to the loops.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       3D-Wärmeleitungsgleichung mit Ghost-Point-Methode und         %%
%%                       Finiter Differenzen                           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Leeren des Workspace und des Editors
clc;
close all;
format long;
%%
% Abmessungen und Schrittweiten des Bleches im Raum
NX = 121;                            % Schrittzahl in x-Richtung  
NY = 121;                            % Schrittzahl in y-Richtung
NZ = 33;                             % Schrittzahl in y-Richtung
XMAX = 30E-3;                        % Abmessung x-Richtung [m]
YMAX = 30E-3;                        % Abmessung y-Richtung [m]
ZMAX = 8E-3;                         % Abmessung z-Richtung [m]
dx = XMAX/(NX-1);                    % Schrittweite in x-Richtung [m]
dy = YMAX/(NY-1);                    % Schrittweite in y-Richtung [m]
dz = ZMAX/(NZ-1);                    % Schrittweite in z-Richtung [m]
x = 0:dx:XMAX;                       % Vektor mit x-Werten
y = 0:dy:YMAX;                       % Vektor mit y-Werten
z = 0:dz:ZMAX;                       % Vektor mit Z-Werten
% Schrittweiten in der Zeit                          
dt = 1E-4;                           % Zeitschritt [s]
NT = 5E11;                           % Anzahl der Zeitschritte
% Laserparameter
P = 160000;                           % Laserleistung [W]
DIST = 10E-3;                        % Abtaststrecke [m]
SPOTD = 60E-6;                       % Spotdurchmesser [m]
% Materialdaten Aluminium
DENS = 2700;                         % Dichte [kg*m^-3]
K_ALU = 180;                         % Wärmeleitfähigkeit Alu [W*(m*K)^-1]
C = 895;                             % spez. Wärmekapazität [J*K^-1 ]
kappa = K_ALU/(DENS*C);              % Temperaturleitfähigkeit [m^2*s^-1]
ALPHA = 0.07;                        % Absorptionskoeffizient
% Materialdaten Luft im Riss
K_AIR = 0.025;                       % Wärmeleitfähigkeit Luft [W*(m*K)^-1]
% Variablen für die Ghost-Point-Methode
delta = 10E-6;                       % Breite Riss [m]
EPS = ((K_ALU)/(K_AIR)-1)*delta;     % Relation K_ALU, K_AIR, delta
a = (5*(EPS)+6*dx)/(EPS+dx);         % Faktor a
b = (dx)/(EPS+dx);                   % Faktor b
% Speicherallokation für die Temperatur-Matrix
T_OLD = zeros(NX,NY,NZ);             % Allokation alte Temperaturen
T_NEW = zeros(NX,NY,NZ);             % Allokation neue Temperaturen
T_AMB = 30;                          % Umgebungstemperatur
% Speicherallokation für die Last-Matrix
q = zeros(NX,NY,NZ);                 % Allokation der Lasten  
%%
% Anfangsbedingung (Blechtemperatur)
for i=1:NX
    for j=1:NY
        for k=1:NZ
            T_OLD(i,j,k)= T_AMB;
        end
    end
end
%%
% Instationärer Wärmestrom (Wärmestromdichte durch Line-Scan)
for i=1:NX
    for j=1:NY
        for k=1:NZ
            if ((j>=40) && (j<=80) && (i==60) && (k==33))
                q(i,j,k)=kappa*ALPHA*((P)/(DIST*SPOTD))/(K_ALU);  
            else
                q(i,j,k)=0;
            end
        end
    end
end
%%
% Berechnung der Feldvariablen für jeden Zeitschritt
for it = 0:NT
    clf;                                 % Löscht aktuelle Figure
    T_NEW = T_OLD;                       % setze T_NEW als T_OLD
    h = slice(x,y,z,T_OLD,...            % Plotting der Feldvariablen
      [],[],[8E-3]); 
    colormap jet;                        % Farbschema der Farbskala
    colorbar('location','eastoutside'... % Position und Größe Farbschema
             ,'fontsize',12);
    shading interp                       % Interpolation zwichen Schritten
    axis ([0 30E-3 0 30E-3 0 8E-3])      % Achsenskalierung
    %alpha(0.5);

    % Achsbeschriftungen
    title({['LST for crack detection using finite difference 3D Heat-'...
            'Diffusion'];['and ghost point method'] ;['time (\itt) = '...
            ,num2str(it*dt) 's']})
    xlabel('x in [m]')
    ylabel('y in [m]')
    zlabel('z in [m]')  

    view(2);                             % Darstellung (1D, 2D, oder 3D)
    drawnow;                             % Aktualisiert die Figure
    pause(1E-40)                         % Pause zwischen einzelnen Figures
    refreshdata(h)                       % Aktualisiert die Daten in Figure

    % Explizites Finite-Differenzen-Verfahren (mittels zentralem DQ)

    for i=2:NX-1
    for j=2:NY-1
    for k=1:NZ
        if((j>=52) && (j<=68) && (i==65) && (k==NZ))
            T_OLD(i,j,k) = T_NEW(i,j,k)+(kappa*dt)/(dx^2)*(b*T_NEW(i+1,j,k)-...
                         a*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_AMB+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);   
        elseif((j>=52) && (j<=68) && (i==65) && (k<NZ) && (k>=15))   
            T_OLD(i,j,k) = T_NEW(i,j,k)+(kappa*dt)/(dx^2)*(b*T_NEW(i+1,j,k)-...
                         a*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_NEW(i,j,k+1)+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);  
        elseif(k==1)
            T_OLD(i,j,k) = T_NEW(i,j,k)+(kappa*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_NEW(i,j,k+1)+T_AMB)+...
                         dt*q(i,j,k);
        elseif((k==NZ) && (j<52))
            T_OLD(i,j,k) = T_NEW(i,j,k)+(kappa*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_AMB+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);
        elseif((k==NZ) && (j>68))
            T_OLD(i,j,k) = T_NEW(i,j,k)+(kappa*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_AMB+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);                    
        elseif((k==NZ) && (j>=52) && (j<=68) && (i~=65))
            T_OLD(i,j,k) = T_NEW(i,j,k)+(kappa*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_AMB+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);                            
        else     
            T_OLD(i,j,k) = T_NEW(i,j,k)+(kappa*dt)/(dx^2)*(T_NEW(i+1,j,k)-...
                         6*T_NEW(i,j,k)+T_NEW(i-1,j,k)+T_NEW(i,j+1,k)+...
                         T_NEW(i,j-1,k)+T_NEW(i,j,k+1)+T_NEW(i,j,k-1))+...
                         dt*q(i,j,k);

       end
    end
    end
    end
end
%% Programm Ende
user6539314
  • 33
  • 5