Exponenta Banner

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

Скрипт Open Live Script

В этом примере показано, как решить уравнение тепла с характеристиками выброса.

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

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

Сгенерируйте и постройте mesh.

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) команда.

Документация Partial Differential Equation Toolbox

Поддержка