3
我一直在写一个有限差分代码，用于使用激光诱导热成像进行模拟和裂纹检测。裂缝由因子a和b实现，这些因子通过使用鬼点方法“阻尼”通过充气裂缝的热流。二维模型按预期运行，稳定条件满足，一切正常。它甚至可以用实验数据证明。只需复制并粘贴即可使用。将我的有限差分模型从2D更改为3D会导致不稳定行为
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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值比预期增加了数量级。我已经尝试了一切。使用一个更简单的负载，评论“裂缝循环”，但没有奏效。
计算的稳定性条件，
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
先谢谢了，我对这个问题非常绝望。
问候亚历
非常感谢，这是一个令人厌恶的简单的错误... – user6539314
另一个尴尬刚刚出现。由于负载是作为瞬态热流施加的，而不是Neumann-B.C，所以我必须处理边界。正如你在脚本中看到的那样，我假设每个时间增量的区域外微分的网格点都是环境温度。实际上这是行不通的，因为k = NZ的点也必须升温，这几乎不会发生。在2D中这样做不是问题，因为在Z方向上有n个渐变。 – user6539314
嘿，学习它几乎是一样的方式;）这是一个有趣的错误:)思考Stack规则，你可能不得不作为一个单独的问题发布这个问题。我不介意，但有些人可能会不同意它在这篇文章中作为答案。也许你可以保留它并在稍后移动它。我会尽快回复！ :) – atru