Этот пример показывает, как решить уравнение тепла с исходным термином.
Основное уравнение тепла с членом источника модуля является
Это уравнение решается в квадратной области с прерывистым начальным условием и нулевой температурой на контурах.
Создайте переходную тепловую модель.
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) команда.