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