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