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

Решение для устойчивого состояния

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

Создайте тепловую модель для установившегося анализа.

thermalModelS = createpde('thermal');

2D модель является прямоугольной полосой, y-размерность которой расширяет от оси симметрии до наружной поверхности, и x-размерность расширяет по фактической длине стержня (от-1.5 м до 1,5 м). Геометрия и mesh для этого прямоугольного раздела легко заданы путем определения местоположений X и Y этих четырех углов как показано ниже.

g = decsg([3 4 -1.5 1.5 1.5 -1.5 0 0 .2 .2]');

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

geometryFromEdges(thermalModelS,g);

Стержень состоит из материала со следующими тепловыми свойствами.

k = 40; % thermal conductivity, W/(m-degree C)
rho = 7800; % density, kg/m^3
cp = 500; % specific heat, W-s/(kg-degree C)
q = 20000; % heat source, W/m^3

Тулбокс УЧП позволяет определение непостоянных коэффициентов как функция пространственных координат и решения.

kFunc = @(region,state) k*region.y;
cFunc = @(region,state) cp*region.y;
qFunc = @(region,state) q*region.y;

Для установившегося анализа задайте теплопроводность материала.

thermalProperties(thermalModelS,'ThermalConductivity',kFunc);

Задайте внутренний источник тепла.

internalHeatSource(thermalModelS,qFunc);

При определении граничных условий ниже, необходимо знать числа ребра для граничных ребер геометрии. Удобный способ получить эти числа ребра состоит в том, чтобы построить геометрию с помощью pdegplot с опцией edgeLabels установите на 'on'.

figure
pdegplot(thermalModelS,'EdgeLabels','on');
axis equal
xlim([-2 2]);
title 'Rod Section Geometry With Edge Labels Displayed';

Задайте граничные условия. Ребро 1, который является ребром в $$y$$ равняйтесь нулю, приезжает ось симметрии, таким образом, нет никакого тепла, переданного в направлении, нормальном к этому ребру. Этот контур моделируется значением по умолчанию как изолированный контур. Ребро 2 сохранено при постоянной температуре T = 100 C. Граничные условия для ребер 3 и 4 являются функциями y.

thermalBC(thermalModelS,'Edge',2,'Temperature',100);

outerCC = @(region,~) 50*region.y;
thermalBC(thermalModelS,'Edge',3,...
                       'ConvectionCoefficient',outerCC,...
                       'AmbientTemperature',100);
                   
leftHF = @(region,~) 5000*region.y;
thermalBC(thermalModelS,'Edge',4,'HeatFlux',leftHF);

Сгенерируйте mesh.

generateMesh(thermalModelS,'Hmax',0.1);
figure; 
pdeplot(thermalModelS); 
axis equal
title 'Rod Section With Triangular Element Mesh'

Решите модель и постройте результат.

result = solve(thermalModelS);
T = result.Temperature;
figure; 
pdeplot(thermalModelS,'XYData',T,'Contour','on'); 
axis equal
title 'Steady State Temperature';

Переходное решение

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

thermalModelT = createpde('thermal','transient');

g = decsg([3 4 -1.5 1.5 1.5 -1.5 0 0 .2 .2]');
geometryFromEdges(thermalModelT,g);

generateMesh(thermalModelT,'Hmax',0.1);

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

thermalProperties(thermalModelT,'ThermalConductivity',kFunc,...
                                'MassDensity',rho,...
                                'SpecificHeat',cFunc);

Задайте внутренний источник тепла и граничные условия.

internalHeatSource(thermalModelT,qFunc);

thermalBC(thermalModelT,'Edge',2,'Temperature',100);
thermalBC(thermalModelT,'Edge',3,...
                        'ConvectionCoefficient',outerCC,...
                        'AmbientTemperature',100);
thermalBC(thermalModelT,'Edge',4,'HeatFlux',leftHF);

Вычислите переходное решение в течение времен решения от $$t=0$$ к $$t=\mathrm{}\mathrm{}\mathrm{}\mathrm{}20000$$ секунды. Начальная температура в стержне является нулем.

tfinal = 20000;
tlist = 0:100:tfinal;
thermalIC(thermalModelT,0);
thermalModelT.SolverOptions.ReportStatistics = 'on';

result = solve(thermalModelT,tlist);

114 successful steps
1 failed attempts
232 function evaluations
1 partial derivatives
22 LU decompositions
231 solutions of linear systems

T = result.Temperature;

Постройте решение в t = 20000.

figure; 
pdeplot(thermalModelT,'XYData',T(:,end),'Contour','on'); 
axis equal
title(sprintf('Transient Temperature at Final Time (%g seconds)',tfinal));

Решение для устойчивого состояния и переходное решение в 20 000 секунд находятся в близком соглашении. Это видно путем сравнения двух фигур.

Третья фигура ниже показов температура в левом конце стержня как функция времени. Наружная поверхность стержня отсоединена среде с постоянной температурой 100 градусов по Цельсию. Следовательно, когда температура поверхности стержня меньше 100, стержень нагревается средой и, когда больше, чем 100, охлажденный. Когда температура стержня меньше 100 градусов, наружная поверхность немного теплее, чем внутренняя ось, но когда температура является приблизительно 100 градусами, наружная поверхность становится более холодной, чем внутренняя часть стержня.

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

p = thermalModelT.Mesh.Nodes;
nodesLeftEnd  = find(p(1,:) < -1.5+eps);
nodeCenter = nodesLeftEnd(p(2,nodesLeftEnd) < eps);
nodeOuter = nodesLeftEnd(p(2,nodesLeftEnd) > 0.2-eps);

figure; 
plot(tlist,T(nodeCenter,:)); 
hold all 
plot(tlist,T(nodeOuter,:),'--');
title 'Temperature at Left End as a Function of Time'
xlabel 'Time, seconds'
ylabel 'Temperature, degrees-C'
grid on;
legend('Center Axis','Outer Surface','Location','SouthEast');