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

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

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

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

$${}_{\mathrm{Qc}}{=}_{}\mathrm{hc}{(}_{}$$T-Ta),

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

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

$${}_{}\mathrm{Qr=ϵ\sigma }{(}^{\mathrm{}}{}_{}^{\mathrm{}}Т4-Та4$$),

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

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

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

Определение параметров PDE

Настройте проблему PDE путем определения значений параметров. Пластина состоит из меди, которая обладает следующими свойствами.

kThermal = 400;             % thermal conductivity of copper, W/(m-K)
rhoCopper = 8960;           % density of copper, kg/m^3
specificHeat = 386;         % specific heat of copper, J/(kg-K)
thick = 0.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)
tAmbient = 300;             % the ambient temperature
emiss = 0.5;                % emissivity of the plate surface

Извлечение коэффициентов PDE

Определите PDE в символьной форме с температурой пластины в качестве зависимой переменной $$T(t,x,$$y).

syms T(t,x,y)
syms eps sig tz hc Ta rho Cp k
Qc = hc*(T - Ta);
Qr = eps*sig*(T^4 - Ta^4);
pdeeq = (rho*Cp*tz*diff(T,t) - k*tz*laplacian(T,[x,y]) + 2*Qc + 2*Qr)

pdeeq(t, x, y) = 
$$2\hspace{0.17em}\mathrm{eps}\hspace{0.17em}\mathrm{sig}\hspace{0.17em}\left({T\left(t,x,y\right)}^4-{\mathrm{Ta}}^4\right)-k\hspace{0.17em}\mathrm{tz}\hspace{0.17em}\left(\frac{{\partial }^2}{\partial x^2}\mathrm{ }T\left(t,x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }T\left(t,x,y\right)\right)-2\hspace{0.17em}\mathrm{hc}\hspace{0.17em}\left(\mathrm{Ta}-T\left(t,x,y\right)\right)+\mathrm{Cp}\hspace{0.17em}\rho \hspace{0.17em}\mathrm{tz}\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }T\left(t,x,y\right)$$

Теперь создайте коэффициенты для использования в качестве входных данных в модели PDE в соответствии с требованиями панели инструментов дифференциальных уравнений в частных производных. Для этого сначала извлеките коэффициенты символьного PDE как структуру символьных выражений с помощью pdeCoefficients функция.

symCoeffs = pdeCoefficients(pdeeq,T,'Symbolic',true)

symCoeffs = struct with fields:
    m: [1x1 sym]
    a: [1x1 sym]
    c: [2x2 sym]
    f: [1x1 sym]
    d: [1x1 sym]

Затем замените символьные переменные, представляющие параметры PDE, их числовыми значениями.

symVars = [eps sig tz hc Ta rho Cp k];
symVals = [emiss stefanBoltz thick hCoeff tAmbient rhoCopper specificHeat kThermal];
symCoeffs = subs(symCoeffs,symVars,symVals);

Наконец, так как поля symCoeffs являются символическими объектами, используйте pdeCoefficientsToDouble для преобразования коэффициентов в double тип данных, который делает их действительными входами для панели инструментов дифференциального уравнения в частных производных.

coeffs = pdeCoefficientsToDouble(symCoeffs)

coeffs = struct with fields:
    a: @makeCoefficient/coefficientFunction
    c: [4x1 double]
    m: 0
    d: 3.4586e+04
    f: 1.0593e+03

Задание модели PDE, геометрии и коэффициентов

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

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

numberOfPDE = 1;
model = createpde(numberOfPDE);

Укажите геометрию модели PDE - в данном случае размер квадрата.

width = 1; 
height = 1;

Определите матрицу описания геометрии. Создайте квадратную геометрию с помощью decsg(Панель инструментов дифференциальных уравнений в частных производных). Для прямоугольной геометрии первая строка содержит 3, и вторая строка содержит 4. Следующие четыре строки содержат $$$$координаты X начальных точек ребер, а четыре строки после этого содержат $$$$координаты Y начальных точек ребер.

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

Преобразуйте геометрию DECSG в объект геометрии и включите его в модель PDE.

geometryFromEdges(model,g);

Постройте график геометрии и отобразите метки кромок.

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

Затем создайте треугольную сетку в модели PDE с размером сетки приблизительно 0,1 в каждом направлении.

hmax = 0.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');

Укажите коэффициенты в модели PDE.

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

Поиск временного решения

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

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

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

setInitialConditions(model,0);
setInitialConditions(model,1000,'edge',1);

Определите временную область для поиска временного решения проблемы PDE.

endTime = 10000;
tlist = 0:50:endTime;

Задайте допуск опций решателя.

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 (s)'
ylabel 'Temperature (K)'

Основываясь на графике, переходное решение начинает достигать своего устойчивого состояния через 6000 секунд. Равновесная температура верхнего края достигает 450 К через 6000 секунд.

Показать температурный профиль пластины через 10000 секунд.

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;

Показать температуру на верхнем краю в 10000 секунд.

u(3,end)

ans = 449.8031