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

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

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

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

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

Создайте сетку с элементами размером не более 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);

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

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

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

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

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

Создайте сетку с элементами размером не более 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'