Проблема теплопередачи с температурно-зависимыми свойствами

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

Дифференциальное уравнение с частными производными для переходной теплопередачи проводимости:

$ρ C_{p} \frac{\partial T}{\partial t} - \nabla \cdot (k \nabla T) = f$

где $T$ температура, $ρ$ существенная плотность, $C_{p}$ удельная теплоемкость, и $k$ теплопроводность. $f$ тепло, выработанное в теле, которое является нулем в этом примере.

Установившееся решение: постоянная теплопроводность

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

thermalmodelS = createpde('thermal','steadystate');

Создайте 2D геометрию путем рисования одного прямоугольника размер блока и второго прямоугольника размер паза.

r1 = [3 4 -.5 .5 .5 -.5  -.8 -.8 .8 .8];
r2 = [3 4 -.05 .05 .05 -.05  -.4 -.4 .4 .4];
gdm = [r1; r2]';

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

g = decsg(gdm,'R1-R2',['R1'; 'R2']');

Преобразуйте decsg формат в геометрический объект. Включайте геометрию в модель.

geometryFromEdges(thermalmodelS,g);

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

figure
pdegplot(thermalmodelS,'EdgeLabels','on'); 
axis([-.9 .9 -.9 .9]);
title 'Block Geometry With Edge Labels Displayed'

Установите температуру на левом крае до 100 градусов. На правом краю существует предписанный поток тепла из блока. Верхние и нижние ребра и ребра в полости все изолируются, то есть, никакое тепло не передается через эти ребра.

thermalBC(thermalmodelS,'Edge',6,'Temperature',100);
thermalBC(thermalmodelS,'Edge',1,'HeatFlux',-10);

Задайте теплопроводность материала. Во-первых, полагайте, что постоянная теплопроводность, например, равняется тому. Позже, рассмотрите случай, где теплопроводность является функцией температуры.

thermalProperties(thermalmodelS,'ThermalConductivity',1);

Создайте mesh с элементами, не больше, чем 0,2.

generateMesh(thermalmodelS,'Hmax',0.2);
figure 
pdeplot(thermalmodelS); 
axis equal
title 'Block With Finite Element Mesh Displayed'

Вычислите установившееся решение.

R = solve(thermalmodelS);
T = R.Temperature;
figure
pdeplot(thermalmodelS,'XYData',T,'Contour','on','ColorMap','hot'); 
axis equal
title 'Temperature, Steady State Solution'

Переходное решение: постоянная теплопроводность

Создайте переходную тепловую модель и включайте геометрию.

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

r1 = [3 4 -.5 .5 .5 -.5  -.8 -.8 .8 .8];
r2 = [3 4 -.05 .05 .05 -.05  -.4 -.4 .4 .4];
gdm = [r1; r2]';
g = decsg(gdm,'R1-R2',['R1'; 'R2']');
geometryFromEdges(thermalmodelT,g);

Задайте теплопроводность, массовую плотность и удельную теплоемкость материала.

thermalProperties(thermalmodelT,'ThermalConductivity',1,...
                                'MassDensity',1,...
                                'SpecificHeat',1);

Задайте граничные условия. В переходных случаях температура на левом крае является нулем в time=0 и пандусах до 100 градусов более чем.5 секунд. Можно найти функцию помощника transientBCHeatedBlock под matlab/R20XXx/examples/pde/main.

thermalBC(thermalmodelT,'Edge',6,'Temperature',@transientBCHeatedBlock);

На правом краю существует предписанный поток тепла из блока.

thermalBC(thermalmodelT,'Edge',1,'HeatFlux',-10);

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

Создайте mesh с элементами, не больше, чем 0,2.

msh = generateMesh(thermalmodelT,'Hmax',0.2);
figure 
pdeplot(thermalmodelT); 
axis equal
title 'Block With Finite Element Mesh Displayed'

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

tlist = 0:.1:5;
thermalIC(thermalmodelT,0);
R = solve(thermalmodelT,tlist);
T = R.Temperature;

Два графика полезны в понимании результатов этого анализа переходных процессов. Первым является график температуры в итоговое время. Вторым является график температуры в отдельном моменте в блоке, в этом случае около центра правого края, как функция времени. Чтобы идентифицировать узел около центра правого края, удобно задать эту короткую служебную функцию.

getClosestNode = @(p,x,y) min((p(1,:) - x).^2 + (p(2,:) - y).^2);

Вызовите эту функцию, чтобы получить узел около центра правого края.

[~,nid] = getClosestNode( msh.Nodes, .5, 0 );

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

h = figure;
h.Position = [1 1 2 1].*h.Position;
subplot(1,2,1); 
axis equal
pdeplot(thermalmodelT,'XYData',T(:,end),'Contour','on','ColorMap','hot'); 
axis equal
title 'Temperature, Final Time, Transient Solution'
subplot(1,2,2); 
axis equal
plot(tlist, T(nid,:)); 
grid on
title 'Temperature at Right Edge as a Function of Time';
xlabel 'Time, seconds'
ylabel 'Temperature, degrees-Celsius'

Решение для устойчивого состояния: температурно-зависимая теплопроводность

Свойствам материала весьма свойственно быть функциями зависимых переменных. Например, примите, что теплопроводность является простой линейной функцией температуры:

k = @(~,state) 0.3+0.003*state.u;

В этом случае переменная u является температурой. В данном примере примите, что плотность и удельная теплоемкость не являются функциями температуры.

thermalProperties(thermalmodelS,'ThermalConductivity',k);

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

R = solve(thermalmodelS);
T = R.Temperature;
figure
pdeplot(thermalmodelS,'XYData',T,'Contour','on','ColorMap','hot');  
axis equal
title 'Temperature, Steady State Solution'

Переходное решение: температурно-зависимая теплопроводность

Теперь выполните анализ переходных процессов с температурно-зависимой проводимостью.

thermalProperties(thermalmodelT,'ThermalConductivity',k,...
                                'MassDensity',1,...
                                'SpecificHeat',1);

Используйте тот же промежуток tlist = 0:.1:5 что касается линейного случая.

thermalIC(thermalmodelT,0);
R = solve(thermalmodelT,tlist);
T = R.Temperature;

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

h = figure;
h.Position = [1 1 2 1].*h.Position;
subplot(1,2,1); 
axis equal
pdeplot(thermalmodelT,'XYData',T(:,end),'Contour','on','ColorMap','hot'); 
axis equal
title 'Temperature, Final Time, Transient Solution'
subplot(1,2,2); 
axis equal
plot(tlist(1:size(T,2)), T(nid,:)); 
grid on
title 'Temperature at Right Edge as a Function of Time (Nonlinear)';
xlabel 'Time, seconds'
ylabel 'Temperature, degrees-Celsius'

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