exponenta event banner

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

В этом примере показано, как решить тепловое уравнение с помощью члена источника.

Основное уравнение тепла с единичным термином источника

∂u∂t-Δu=1

Это уравнение решается в квадратной области с прерывистым начальным условием и нулевыми температурами на границах.

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

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

Создание квадратной геометрии с центром в x = 0 и y = 0 со сторонами длины 2. Включите окружность радиуса 0,4, концентричную квадрату.

R1 = [3;4;-1;1;1;-1;-1;-1;1;1];
C1 = [1;0;0;0.4];
C1 = [C1;zeros(length(R1) - length(C1),1)];
gd = [R1,C1];
sf = 'R1+C1';
ns = char('R1','C1')';
g = decsg(gd,sf,ns);

Добавьте геометрию к модели.

geometryFromEdges(thermalmodel,g);

Укажите тепловые свойства материала.

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

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

internalHeatSource(thermalmodel,1);

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

figure
pdegplot(thermalmodel,'EdgeLabels','on','FaceLabels','on')
axis([-1.1 1.1 -1.1 1.1]);
axis equal
title 'Geometry With Edge and Subdomain Labels'

Figure contains an axes. The axes with title Geometry With Edge and Subdomain Labels contains 11 objects of type line, text.

Установите нулевые температуры на всех четырех внешних краях квадрата.

thermalBC(thermalmodel,'Edge',1:4,'Temperature',0);

Прерывистое начальное значение равно 1 внутри круга и нулю снаружи. Укажите нулевую начальную температуру везде.

thermalIC(thermalmodel,0);

Задайте ненулевую начальную температуру внутри окружности (Грань 2).

thermalIC(thermalmodel,1,'Face',2);

Создание и печать сетки.

msh = generateMesh(thermalmodel);
figure;
pdemesh(thermalmodel); 
axis equal

Figure contains an axes. The axes contains 2 objects of type line.

Найдите решение в 20 точках времени от 0 до 0,1.

nframes = 20;
tlist = linspace(0,0.1,nframes);

thermalmodel.SolverOptions.ReportStatistics ='on';
result = solve(thermalmodel,tlist);
99 successful steps
0 failed attempts
200 function evaluations
1 partial derivatives
20 LU decompositions
199 solutions of linear systems
T = result.Temperature;

Постройте график решения.

figure
Tmax = max(max(T));
Tmin = min(min(T));
for j = 1:nframes
    pdeplot(thermalmodel,'XYData',T(:,j),'ZData',T(:,j));
    caxis([Tmin Tmax]);
    axis([-1 1 -1 1 0 1]);
    Mv(j) = getframe;
end

Figure contains an axes. The axes contains an object of type patch.

Для воспроизведения анимации используйте movie(Mv,1) команда.