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

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

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

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

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

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

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

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

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

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

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

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

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

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

$ρ 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}$ толщина пластины, и факторы два составляют теплопередачу от обеих поверхностей пластины.

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

$ρ 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 object. The axes object with title Geometry With Edge Labels Displayed contains 5 objects of type line, text.

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

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'

Figure contains an axes object. The axes object 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 object. The axes object 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 object. The axes object 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 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));

Figure contains an axes object. The axes object with title ['Temperature Along the Top Edge of ' contains an object of type line.

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

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

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

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

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