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

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

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

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

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

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

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

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

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

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

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

Задайте параметры УЧП

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

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]

Затем замените символьными переменными, которые представляют параметры УЧП их числовыми значениями.

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

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

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

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

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

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

На основе графика переходное решение начинает достигать своего значения устойчивого состояния после 6 000 секунд. Температура равновесия верхнего края приближается к 450 K после 6 000 секунд.

Покажите температурный профиль пластины после 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