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

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

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

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

где $${}_{\mathrm{Ta}}$$ - температура окружающей среды, $$T$$ - температура в конкретном положении x и y на поверхности пластины, и $${}_{\mathrm{hc}}$$ - заданный коэффициент конвекции.

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

где $$ϵ$$ - эмиссионная способность лица, $$\mathrm{а\; \lambda }$$ - постоянная Штефана - Больцмана. Поскольку теплота, передаваемая из-за излучения, пропорциональна четвертой мощности поверхностной температуры, проблема нелинейна.

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

где $$\mathrm{start}$$- плотность материала, $${}_{\mathrm{Cp}}$$ - удельная теплота, $${\mathrm{tz}}_{}$$- толщина пластины, а коэффициенты двух учитывают теплопередачу от обеих граней пластины.

Это уравнение удобно переписывать в форме, ожидаемой PDE Toolbox

Настройка проблемы

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

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);

Для квадрата геометрия и сетка легко определяются, как показано ниже.

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';

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

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);

Нижний край пластины установлен на 1000 градусов по Кельвину.

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

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

Укажите начальное предположение.

setInitialConditions(model,0);

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

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 К.

u0(1:numNodes) = 300;

Установите начальную температуру на нижней кромке E1 на величину постоянной ВС, 1000 К.

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

Резюме

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