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

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

Пластина является квадратной, и температура фиксируется вдоль базового края. Никакое тепло не передается от других трех ребер (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';

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

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

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