Изменение моей конечно-разностной модели с 2D на 3D вызывает нестабильное поведение

Я писал конечно-разностный код для моделирования и обнаружения трещин с помощью лазерной термографии. Трещина реализуется факторами a и b, которые «демпфируют» тепловой поток через заполненную воздухом трещину с использованием метода призрачной точки. 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 для стабильного поведения увеличивается на порядки больше, чем ожидалось. Я пробовал все. Используя более простую загрузку, закомментировав "трещину-петлю", но ничего не получилось.

Вычисление условия устойчивости,

dt <= dx^2/(6*k) = 1.39E-4   instead of 2E-10(!!!)

должно быть достаточно. Но попробуйте 2Е-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 09.08.2017 источник

Ответы (2)

arrow_upward
2
arrow_downward

У вас есть ошибка в вашем коде — в 3D-версии вы вводите циклическую переменную с именем k для измерения z. Эта переменная перезаписывает ваш ранее определенный коэффициент k. При фиксе все работает с dt=1e-4 с в 3D. Я просто изменил k, служащий переменной цикла, на kj. Вы можете выбрать лучшее имя. На самом деле рекомендуется использовать более длинные имена для циклических переменных, а не только i, j, k... - например, две или три буквы вместо одной.

atru 10.08.2017

comment

Хех, выучил это почти так же;) это забавная ошибка :) Думая о правилах стека, вам, возможно, придется опубликовать этот вопрос как отдельный вопрос. Я не возражаю, но есть вероятность, что некоторые люди не согласятся с тем, что это ответ в этом посте. Может быть, вы можете сохранить его и переместить его позже. Постараюсь ответить в ближайшее время! :) - atru; 10.08.2017

comment

Не могу прокомментировать ваш ответ, так как моя репутация все еще низкая! В физической ситуации, каковы ваши фактические условия в z? Температура окружающей среды или изолированная (без потока) граница? Если это no-flux (Neumann), то аналогично тому, что вы написали, вы должны использовать T_NEW(i,j,k+1) = T_NEW(i,j,k), т.е. граничный узел равен предыдущему по z. Для нижней границы можно сделать то же самое, только в другом порядке: T_NEW(i,j,1) = T_NEW(i,j,2). - atru; 10.08.2017

comment

Именно, на поверхности нет никакого теплового потока (k=NZ), за исключением нестационарного теплового потока, вызванного лазерным лучом, проводящим в направлениях x, y и z. Конвекция и излучение не учитываются. Что меня беспокоит, так это разница между 2D- и 3D-распределением тепла с использованием одних и тех же параметров. Возможно, это связано с тем, что теплоемкость 3D-модели намного выше, что приводит к меньшему нагреву, потому что сыпучий материал должен быть нагрет. Фактически это приводит к увеличению мощности лазера для масштабирования графиков, чтобы они совпадали. - user6539314; 11.08.2017

comment

Я не проверял вашу схему строго, но я использовал подобную схему несколько лет назад, издалека выглядит нормально. Кажется логичным, что вам потребуется больше тепла от источника тепла в трехмерном домене, поскольку он больше. За исключением случаев, когда ваш источник тепла охватывал всю высоту домена, и у вас не было бы потока или периодических границ сверху и снизу. Я вижу, что он охватывает только один узел в z, а нижняя граница является внешней. Одна вещь, которая приходит на ум — вы проводили исследование независимости сети? Кроме того, отсутствие границы потока недопустимо, если в реальной жизни эта часть системы изолирована. - atru; 11.08.2017

comment

Мне просто нужно было посмотреть, что означает сетка Independence Study, но, читая об этом, кажется, это просто важно для работы с жидкостями? Я думал, что числовая задача вообще не может быть независимой от сетки, потому что более грубая сетка влияет на точность решения. Так что, на мой взгляд, каждое различие в шаге сетки приводит к немного другому решению. Есть ли у вас опыт моделирования источников тепла? На самом деле, я должен применить всю интенсивность (которая относится к площади) лазера только на одной линии сетки (которая, очевидно, не является площадью), поскольку диаметр лазера составляет всего около 90 микрометров. - user6539314; 16.08.2017

comment

Это хорошо, и это важно для любого численного моделирования, другое название которого, например, исследование с измельчением сетки. С уточнением ваше решение будет приближаться к правильному решению. Затем вы можете выбрать, с какой сеткой вы хотите работать, исходя из ваших требований к точности и ограничений вычислительной нагрузки/эффективности. В простых задачах вы действительно можете достичь этого. С разумной точностью вы действительно можете найти правильное решение. В более сложных задачах нередко приходится вообще отказываться от всей идеи и работать только с одной сеткой, надеясь, что она достаточно близка. - atru; 16.08.2017

comment

Это не обязательно, но в такой небольшой проблеме было бы полезно проверить это. Просто запустите пару сеток с расстоянием, идеально отличающимся в 2 раза, и постройте решение как функцию расстояния. Решением может быть средняя температура в системе или разность интерполированных температурных полей (вам потребуется интерполяция, потому что вы будете сравнивать матрицы с разным количеством элементов). Вы также можете использовать свою лучшую сетку в качестве правильного решения и посмотреть на ошибку относительно нее. Это тоже обычная практика. - atru; 16.08.2017

arrow_upward
0
arrow_downward

Только что возник еще один вопрос. Поскольку нагрузка применяется как переходный тепловой поток, а не как Нейман-БК, мне приходится обрабатывать границы. Как вы можете видеть в моем сценарии, я предположил, что точки сетки для дифференциации вне домена являются температурой окружающей среды для каждого приращения времени. На самом деле это не работает, так как точки при k=NZ тоже должны нагреваться, чего почти не происходит. Сделать это в 2D снова не проблема, потому что в z-направлении нет градиента. Итак, у вас есть совет, как изменить мой код? Я подумал о замене T_AMB на T_NEW(i,j,k), чтобы 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

user6539314 10.08.2017

Изменение моей конечно-разностной модели с 2D на 3D вызывает нестабильное поведение

Ответы (2)

Похожие вопросы