Нелинейная теплопередача в тонком диске

Этот пример показывает, как выполнить анализ теплопередачи тонкого диска.

Пластина квадратная, и температура фиксируется вдоль нижнего ребра. Тепло не передается от остальных трёх ребер (т.е. они изолированы). Тепло передается от верхней и нижней граней пластины конвекцией и излучением. Поскольку излучение включено, проблема нелинейна. Одна из целей этого примера состоит в том, чтобы показать, как справиться с нелинейностями в задачах УЧП.

Выполняются как устойчивое состояние, так и переходный анализ. В устойчивом анализе мы заинтересованы в конечной температуре в разных точках диска после того, как он достиг равновесного состояния. В переходном анализе нас интересует температура в диске как функция времени. Один из вопросов, на который можно ответить при этом переходном анализе, заключается в том, сколько времени требуется для достижения пластиной равновесной температуры.

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

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

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

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

$Q_{c} = h_{c} (T - T_{a})$

где $T_{a}$ - температура окружающей среды, $T$ - температура в определенном местоположении x и y на поверхности диска, и $h_{c}$ является заданным коэффициентом конвекции.

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

$Q_{r} = ϵ σ (T^{4} - T_{a}^{4})$

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

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

$ρ C_{p} t_{z} \frac{\partial T}{\partial t} - k t_{z} \nabla^{2} T + 2 Q_{c} + 2 Q_{r} = 0$

где $ρ$ - плотность материала, $C_{p}$ - удельная теплота, $t_{z}$ - толщина пластины, и коэффициенты двойки учитывают теплопередачу от обеих поверхностей пластины.

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

$ρ C_{p} t_{z} \frac{\partial T}{\partial t} - k t_{z} \nabla^{2} T + 2 h_{c} T + 2 ϵ σ T^{4} = 2 h_{c} T_{a} + 2 ϵ σ T_{a}^{4}$

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

Figure contains an axes. The axes with title Geometry With Edge Labels Displayed contains 5 objects of type line, text.

Задайте коэффициенты. Выражения для коэффициентов, требуемых 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 градусов Кельвина.

Примените граничные условия. Три из пластинчатых ребер являются изолированными. Поскольку граничное условие Неймана, равное нулю, является по умолчанию в формулировке конечного элемента, граничные условия на этих ребрах не должны быть заданы явным образом. Условие Дирихле задается на всех узлах на нижнем ребреребре 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'

Figure contains an axes. The axes with title Plate With Triangular Element Mesh contains 2 objects of type line.

Решение в установившемся состоянии

Поскольку коэффициенты 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

Figure contains an axes. The axes with title Temperature In The Plate, Steady State Solution contains 12 objects of type patch, line.

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

Figure contains an axes. The axes with title Temperature As a Function of the Y-Coordinate contains an object of type line.

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 contains an axes. The axes with title Temperature Along the Top Edge of the Plate as a Function of Time contains an object of type line.

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;

Figure contains an axes. The axes with title Temperature In The Plate, Transient Solution( 5000 seconds) contains 12 objects of type patch, line.

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 секунд переходное решение достигло значений устойчивого состояния. Температура от двух решений на верхнем крае пластины соответствует в пределах одного процента.

Документация