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