3

我一直在编写有限差分代码，用于使用激光诱导热成像来模拟和检测裂缝。裂纹由因素 a 和 b 实现，它们使用鬼点方法“阻尼”通过空气填充裂纹的热流。二维模型按预期运行，满足稳定条件，一切正常。甚至用实验数据很好地证明了这一点。只需复制和粘贴，它就会工作。

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

提前谢谢你，我对这个问题非常绝望。

问候亚历克斯

4

2 回答 2

2

您的代码中有一个错误 - 在 3D 版本中，您为 z 维度引入了一个名为 k 的循环变量。该变量会覆盖您之前定义的 k 系数。修复后，它在 3D 中都适用于 dt = 1e-4 s。我只是将作为循环变量的 k 更改为 kj。你可以选择一个更好的名字。实际上，建议为循环变量使用稍长的名称，而不仅仅是 i、j、k……——比如 2 或 3 个字母而不是 1 个。

于 2017-08-10T03:56:52.610 回答
0

另一个问题刚刚出现。由于负载是作为瞬态热流而不是作为 Neumann-BC 施加的，因此我必须处理边界。正如您在我的脚本中看到的那样，我假设域外的微分网格点是每次增量时的环境温度。实际上这是行不通的，因为 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
于 2017-08-10T13:16:25.387 回答