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

Создайте статическую тепловую модель.

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

Создайте 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);

Задайте граничные условия. В переходных случаях температура на левом краю равна нулю в момент времени = 0 и наклоняется на 100 степени в течение 5 секунд. Вы можете найти функцию helper 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'