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'

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

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

Найдите решение в 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

Чтобы проигрывать анимацию, используйте movie(Mv,1) команда.

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