将我的有限差分模型从 2D 更改为 3D 会导致不稳定的行为
Changing my finite difference model from 2D to 3D causes unstable behavior
我一直在编写有限差分代码，用于使用激光诱导热成像法模拟和检测裂纹。裂缝由因素 a 和 b 实现，它们是 "damping" 使用鬼点法通过空气填充裂缝的热流。 2D 模型按预期运行，满足稳定性条件，一切顺利。甚至用实验数据也很好地证明了这一点。只需复制粘贴即可。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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
但是从 2D 变为 3D，稳定行为的 dt 值比预期增加了多个数量级。我已经尝试了一切。使用更简单的负载，注释掉 "crack-loop"，但没有任何效果。
计算稳定性条件，
dt <= dx^2/(6*k) = 1.39E-4 instead of 2E-10(!!!)
应该够了。但只需尝试 2E-9，方案就会开始振荡。问题是，我需要裂缝下方的热流。这就是我需要 3D 模型的原因，以防万一。但是这样算下来需要数年的时间才算出10到100毫秒，这正是我需要的范围。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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
提前谢谢你，我对这个问题非常绝望。
问候亚历克斯
您的代码中有一个错误 - 在 3D 版本中，您为 z 维度引入了一个名为 k 的循环变量。该变量会覆盖您之前定义的 k 系数。修复后，它都适用于 3D 中的 dt = 1e-4 s。我只是将用作循环变量的 k 更改为 kj。你可以选择一个更好的名字。实际上，建议为循环变量使用稍长的名称，而不仅仅是 i、j、k... - 例如 2 或 3 个字母而不是一个。
刚刚出现另一个问题。由于负载是作为瞬态热流而不是 Neumann-B.C 施加的，我必须处理边界。正如您在我的脚本中看到的那样，我假设域外微分的网格点是每个时间增量的环境温度。实际上这是行不通的，因为 k=NZ 处的点也必须升温，而这几乎不会发生。在 2D 中再次这样做是没有问题的，因为在 z 方向上没有梯度。那么你对如何修改我的代码有什么建议吗？我考虑用 T_NEW(i,j,k) 代替 T_AMB，这样 T_NEW(i,j,k+1) 等于 T_NEW(i,j, k).这给出了一个合理的情节。但同样，我不知道这在数学上是否正确。下面，有一个关于循环的稍微更正的代码。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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
我一直在编写有限差分代码，用于使用激光诱导热成像法模拟和检测裂纹。裂缝由因素 a 和 b 实现，它们是 "damping" 使用鬼点法通过空气填充裂缝的热流。 2D 模型按预期运行，满足稳定性条件，一切顺利。甚至用实验数据也很好地证明了这一点。只需复制粘贴即可。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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
但是从 2D 变为 3D，稳定行为的 dt 值比预期增加了多个数量级。我已经尝试了一切。使用更简单的负载，注释掉 "crack-loop"，但没有任何效果。
计算稳定性条件，
dt <= dx^2/(6*k) = 1.39E-4 instead of 2E-10(!!!)
应该够了。但只需尝试 2E-9，方案就会开始振荡。问题是，我需要裂缝下方的热流。这就是我需要 3D 模型的原因，以防万一。但是这样算下来需要数年的时间才算出10到100毫秒，这正是我需要的范围。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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
提前谢谢你，我对这个问题非常绝望。
问候亚历克斯
您的代码中有一个错误 - 在 3D 版本中，您为 z 维度引入了一个名为 k 的循环变量。该变量会覆盖您之前定义的 k 系数。修复后，它都适用于 3D 中的 dt = 1e-4 s。我只是将用作循环变量的 k 更改为 kj。你可以选择一个更好的名字。实际上，建议为循环变量使用稍长的名称，而不仅仅是 i、j、k... - 例如 2 或 3 个字母而不是一个。
刚刚出现另一个问题。由于负载是作为瞬态热流而不是 Neumann-B.C 施加的，我必须处理边界。正如您在我的脚本中看到的那样，我假设域外微分的网格点是每个时间增量的环境温度。实际上这是行不通的，因为 k=NZ 处的点也必须升温，而这几乎不会发生。在 2D 中再次这样做是没有问题的，因为在 z 方向上没有梯度。那么你对如何修改我的代码有什么建议吗？我考虑用 T_NEW(i,j,k) 代替 T_AMB，这样 T_NEW(i,j,k+1) 等于 T_NEW(i,j, k).这给出了一个合理的情节。但同样，我不知道这在数学上是否正确。下面，有一个关于循环的稍微更正的代码。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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