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

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

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

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

$$Q_c=h_c(T-T_a)$$,

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

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

$$Q_r=ϵ\sigma (T^4-T_a^4)$$,

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

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

где $$\rho$$ - плотность материала пластины, $$C_p$$ его удельная жара, $$t_z$$ его толщина диска, $$k$$ является его теплопроводностью, и факторы двойки учитывают теплопередачу от обеих его граней.

Определение параметров 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

Извлечение коэффициентов УЧП

Определите УЧП в символьной форме с температурой диска как зависимой переменной $$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, как требуется Partial Differential Equation Toolbox. Для этого сначала извлеките коэффициенты символического УЧП как структуру символьных выражений с помощью 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 тип данных, который делает их допустимыми входами для Partial Differential Equation Toolbox.

coeffs = pdeCoefficientsToDouble(symCoeffs)

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

Задайте модель PDE, геометрию и коэффициенты

Теперь, используя Partial Differential Equation Toolbox, решите задачу УЧП, используя анализ конечных элементов на основе этих коэффициентов.

Во-первых, создайте модель PDE с одной зависимой переменной.

numberOfPDE = 1;
model = createpde(numberOfPDE);

Задайте геометрию для модели PDE - в этом случае размерность квадрата.

width = 1; 
height = 1;

Задайте матрицу описания геометрии. Создайте квадратную геометрию с помощью decsg (Partial Differential Equation Toolbox) функция. Для прямоугольной геометрии первая строка содержит 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');

Затем создайте треугольный mesh в модели 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 К. Задайте это с помощью условия Дирихле на всех узлах на нижнем ребре (ребро E1).

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

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

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

Определите временной интервал, чтобы найти переходное решение УЧП задачи.

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 секунд.

Отобразите температурный профиль диска через 10 000 секунд.

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;

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

u(3,end)

ans = 449.8031