Почему наблюдается расхождение результатов с Matlab и Python? ⇐ Python

[Anonymous]

Я хочу решить одномерное уравнение нестационарной теплопередачи.
Определить одномерную геометрию (линию) в направлении Y

Код: Выделить всё

dT/dt = (k/(rho*cp))*d²T/dy²

I.C.: @t=0, температура всех элементов 25°C
B.C.: Геометрия разделена на несколько элементов.

• Нижний элемент подвергается воздействию горячего ролика -

  Код: Выделить всё

  kdT/dy = h(T_roller-T)

• Верхний элемент подвергается воздействию атмосферы –

  Код: Выделить всё

  k*dT/dy = h(T-T_ambient)

Но я получаю разные результаты в файле Matlab и файле Python.
Граничное условие радиации отлично работает для обоих файлов, но ролик в граничном условии контакта дает разные результаты. Результаты различаются на 4-5°C при низкой температуре валика и разница увеличивается примерно до 10-11°C при повышении температуры валика на 100°C.
Код Matlab:

Код: Выделить всё

function pde()
% Constants

Ny = 1000;            % Number of spatial grid points
Nt = 5000;           % Number of time steps
velocity_x = 15/60; % line speed in m/s
% Thermal properties of different layers
layer1_thickness = 0.001;  % Thickness of layer 1 (meters)
layer1_k = 2.15;            % Thermal conductivity of layer 1 (W/m-K)
layer2_thickness = 0.001;  % Thickness of layer 2 (meters)
layer2_k = layer1_k;            % Thermal conductivity of layer 2 (W/m-K)
rho_1 = 1950;           %density of layer 1 (kg/m³)
rho_2 = 1950;           %density of layer 2 (kg/m³)
cp_1 = 2100;            %specific heat capacity of layer 1 (J/kg-K)
cp_2 = 2100;            %specific heat capacity of layer 2 (J/kg-K)
T_roller = 120;         %Temperature of hot roller (°C)
h_roller = 500;         %Convective heat transfer coefficient (W/m2-K)

Ly = layer1_thickness+layer2_thickness;         % Total thickness of the system (meters)
% Ambient temperature and convection properties
T_ambient = 25;    % Ambient temperature (°C)
h_air = 12;        % Convective heat transfer coefficient of air (W/m²-K)

%Discretization in space
ny1 = ceil(layer2_thickness/Ly*Ny);
ny2 = Ny-ny1;

y1 = linspace(0,layer2_thickness,ny1);
y2 = linspace(layer2_thickness,Ly,ny2);
y = [y1,y2(2:end)];

% Initialize temperature matrix
initial_temperature1 = 25;
initial_temperature2 = 25;

radiative_flux = 150e3;
absorption = 50;
performance = 80;
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);

dia_roller = 0.5;                   %diameter of hot roller (m)
contact_angle = 180;                %contact angle/wrap angle for hot_roller in °

y_1 = pi*dia_roller*contact_angle/360;              %in contact with hot roller
y_2 = 1.3;                                          %after hot roller

t1 = y_1/velocity_x;
t2 = y_2/velocity_x;

T = t1+t2;

nt1 = ceil(t1/(t1+t2)*Nt);
nt2 = Nt-nt1;

t11 = linspace(0,t1,nt1);
t12 = linspace(t1,t1+t2,nt2);

%pdepe settings
m = 0;  %for 1-D cartesian coordinates with no symmetry
phase = 1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t11);
u1 = sol(:,:,1);
phase = 2;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t12);
u2 = sol(:,:,1);

plot([t11,t12],[[u1(:,1);u2(:,1)],[u1(:,25);u2(:,25)],[u1(:,50);u2(:,50)]])
grid on
function [c f s] = pdefun(y,t,u,dudy)
if y 

Подробнее здесь: [url]https://stackoverflow.com/questions/78538874/why-is-there-discrepancy-in-the-result-with-matlab-and-python[/url]