Уравнения теплопередачи для пластины

Пластина имеет плоские размерности один метр на один метр и 1 см толщиной. Поскольку пластина является относительно тонкой по сравнению с плоскими размерностями, температура может быть принята постоянная в направлении толщины; получившаяся проблема 2D.

Конвекция и теплопередача излучения приняты, чтобы произойти между двумя поверхностями пластины и заданной температуры окружающей среды.

Количество тепла, переданное от каждой поверхности пластины на единичную площадь из-за конвекции, задано как

где $$T_a$$ температура окружающей среды, $$T$$ температура в конкретном местоположении X и Y на поверхности пластины, и $$h_c$$ заданный коэффициент конвекции.

Количество тепла, переданное от каждой поверхности пластины на единичную площадь из-за излучения, задано как

где $$ϵ$$ излучаемость поверхности и $$\sigma$$ Stefan-постоянная-Больцмана. Поскольку тепло, переданное из-за излучения, пропорционально четвертой степени температуры поверхности, проблема нелинейна.

УЧП, описывающий температуру в этой тонкой пластине,

где $$\rho$$ существенная плотность, $$C_p$$ удельная теплоемкость, $$t_z$$ толщина пластины, и факторы два составляют теплопередачу от обеих поверхностей пластины.

Удобно переписать это уравнение в форме, ожидаемой Тулбоксом УЧП

Setup задач

Пластина состоит из меди, которая имеет следующие свойства:

k = 400; % thermal conductivity of copper, W/(m-K)
rho = 8960; % density of copper, kg/m^3
specificHeat = 386; % specific heat of copper, J/(kg-K)
thick = .01; % plate thickness in meters
stefanBoltz = 5.670373e-8; % Stefan-Boltzmann constant, W/(m^2-K^4)
hCoeff = 1; % Convection coefficient, W/(m^2-K)
% The ambient temperature is assumed to be 300 degrees-Kelvin.
ta = 300;
emiss = .5; % emissivity of the plate surface

Создайте модель PDE с одной зависимой переменной.

numberOfPDE = 1;
model = createpde(numberOfPDE);

Для квадрата геометрия и mesh легко заданы как показано ниже.

width = 1; 
height = 1;

Задайте квадрат путем предоставления этих 4 x-мест, сопровождаемых 4 y-местами углов.

gdm = [3 4 0 width width 0 0 0 height height]';
g = decsg(gdm, 'S1', ('S1')');

Преобразуйте геометрию DECSG в геометрический объект при выполнении, таким образом, это добавлено к PDEModel

geometryFromEdges(model,g);

Постройте геометрию и отобразите метки ребра для использования в определении граничного условия.

figure; 
pdegplot(model,'EdgeLabels','on'); 
axis([-.1 1.1 -.1 1.1]);
title 'Geometry With Edge Labels Displayed';

Задайте коэффициенты. Выражения для коэффициентов, требуемых Тулбоксом УЧП, могут легко быть идентифицированы путем сравнения уравнения выше скалярным параболическим уравнением в документации Тулбокса УЧП.

c = thick*k;

Из-за граничного условия излучения "a" коэффициент является функцией температуры, u. Это задано как выражение MATLAB, таким образом, это может быть оценено для различных значений u во время анализа.

a = @(~,state) 2*hCoeff + 2*emiss*stefanBoltz*state.u.^3;
f = 2*hCoeff*ta + 2*emiss*stefanBoltz*ta^4;
d = thick*rho*specificHeat;
specifyCoefficients(model,'m',0,'d',0,'c',c,'a',a,'f',f);

Базовый край пластины установлен в 1 000 кельвинов степеней.

Примените граничные условия. Три из ребер пластины изолируются. Поскольку Нейманово граничное условие, равный нуль является значением по умолчанию в формулировке конечного элемента, граничных условиях на этих ребрах, не должно быть установлено явным образом. Условие Дирихле установлено на всех узлах на базовом краю, ребро 1,

applyBoundaryCondition(model,'dirichlet','Edge',1,'u',1000);

Задайте исходное предположение.

setInitialConditions(model,0);

Создайте треугольную mesh на квадрате приблизительно с десятью элементами в каждом направлении.

hmax = .1; % element size
msh = generateMesh(model,'Hmax',hmax);
figure; 
pdeplot(model); 
axis equal
title 'Plate With Triangular Element Mesh'
xlabel 'X-coordinate, meters'
ylabel 'Y-coordinate, meters'

Решение для устойчивого состояния

Поскольку a и f коэффициенты являются функциями температуры (из-за граничных условий излучения), solvepde автоматически выбирает нелинейный решатель, чтобы получить решение.

R = solvepde(model);
u = R.NodalSolution;
figure; 
pdeplot(model,'XYData',u,'Contour','on','ColorMap','jet');
title 'Temperature In The Plate, Steady State Solution'
xlabel 'X-coordinate, meters'
ylabel 'Y-coordinate, meters'
axis equal

p = msh.Nodes;
plotAlongY(p,u,0);
title 'Temperature As a Function of the Y-Coordinate'
xlabel 'Y-coordinate, meters'
ylabel 'Temperature, degrees-Kelvin'

fprintf(['Temperature at the top edge of the plate =' ...
         ' %5.1f degrees-K\n'],u(4));

Temperature at the top edge of the plate = 449.8 degrees-K

Переходное решение

Включайте d коэффициент.

specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f);
endTime = 5000;
tlist = 0:50:endTime;
numNodes = size(p,2);

Установите начальную температуру всех узлов к окружающей среде, 300 K.

u0(1:numNodes) = 300;

Установите начальную температуру на базовом краю E1 к значению постоянного BC, 1000 K.

setInitialConditions(model,1000,'Edge',1);

Установите следующие опции решателя.

model.SolverOptions.RelativeTolerance = 1.0e-3; 
model.SolverOptions.AbsoluteTolerance = 1.0e-4;

Решите задачу при помощи solvepde. Решатель автоматически выбирает параболический решатель, чтобы получить решение.

R = solvepde(model,tlist);
u = R.NodalSolution;
figure; 
plot(tlist,u(3, :)); 
grid on
title ['Temperature Along the Top Edge of ' ...
       'the Plate as a Function of Time']
xlabel 'Time, seconds'
ylabel 'Temperature, degrees-Kelvin'

figure;
pdeplot(model,'XYData',u(:,end),'Contour','on','ColorMap','jet');
title(sprintf(['Temperature In The Plate,' ...
               'Transient Solution( %d seconds)\n'],tlist(1,end)));
xlabel 'X-coordinate, meters'
ylabel 'Y-coordinate, meters'
axis equal;

fprintf(['\nTemperature at the top edge(t = %5.1f secs) = ' ...
         '%5.1f degrees-K\n'],tlist(1,end),u(4,end));

Temperature at the top edge(t = 5000.0 secs) = 441.8 degrees-K

Сводные данные

Графики температуры в пластине от устойчивого состояния и переходного решения в то время окончания очень близки. Таким образом, приблизительно после 5 000 секунд переходное решение достигло значений устойчивого состояния. Температуры из этих двух решений в верхнем краю пластины соглашаются на в одном проценте.