Это часть диплома. Сначала идет код, дальше скрины -- результаты работы программы.
Код полностью рабочий. Если правильно скопируете -- будет работать у Вас.
% Стационарная задача теплопроводности в цилиндрических
% координатах в осесимметричной постановке для составной
% цилиндрическо-сферической толстой оболочки
clear all % очистили память
close all % закрыли все окна
clc % очистили командное окно
% Задаём исходные данные
% данные по области решения
Rext = 1; % наружный радиус, м
Rint = 0.85; % внутренний радиус, м
Hcyl = 3; % длина цилиндрической части, м
% физические константы материала
Kti = 22.3; % внутренняя теплопроводность титана, Вт/(м*К)
% Взято из таблицы
% http://thermalinfo.ru/svojstva-materialov/metally-i-splavy/teploprovodnost-i-teplofizicheskie-svojstva-titana
% В данной задаче эту константу вообще можно не задавать,
% т.к. она сокращается и в диф. уравнении,
% и в граничных условиях.
% данные для граничных условий
Ucyl = 200; % температура на внешней поверхности
% цилиндрической оболочки, К
kUsph = 300; % коэффициент изменения температуры
% в сферической оболочке, К/м
% другие исходные данные
dx = 0.1; % расстояние от области до границы рисунка, м
MyEps = 1e-6; % малая для сравнения действительных чисел
% конец задания исходных данных
% Используются следующие имена файлов:
% ShellArea.png - область решения
% ShellGrid.png - сетка МКЭ
% ShellSol2D.png - решение, вид сверху
% ShellSol3D.png - решение, вид в аксонометрии
HeatModel = createpde('thermal','steadystate'); % создали
% пустую модель для стационарной задачи теплопрповодности
% строим область решения
C1 = [1;0;0;Rext;0;0;0;0;0;0;0;0]; % внешний круг с центром
% в точке O(0,0) и радиусом Rext
C2 = [1;0;0;Rint;0;0;0;0;0;0;0;0]; % внутренний круг с центром
% в точке O(0,0) и радиусом Rint;
% (C1-C2) - это кольцо
R1 = [3;4;0;Rext+dx;Rext+dx;0;0;0;Rext+dx;Rext+dx;0;0];
% квадрат x=[0;Rext+dx]; y=[0;Rext+dx] для вырезания части кольца;
% (C1-C2)*R1 - это правая верхняя четверть кольца
P1 = [2;5;-Hcyl;0;(Rext-Rint)/2;0;-Hcyl;...
Rint;Rint;(Rext+Rint)/2;Rext;Rext]; % цилиндрическая часть
% оболочки с небольшим заходом в 1/4 кольца
gd = [C1 C2 R1 P1]; % матрица из примитивов
ns = char('C1','C2','R1','P1')'; % идентификаторы столбцов
fo = '(C1-C2)*R1+P1'; % формула построения области
[dl1,bt1] = decsg(gd,fo,ns); % формируем область
[dl,bt] = csgdel(dl1,bt1); % удалили внутренние границы
geometryFromEdges(HeatModel,dl); % добавили геометрию в модель
% рисуем область решения и номера граничных линий
figure('Position',[40 40 1200 500]) % новое окно фигуры
pdegplot(HeatModel,'EdgeLabels','on') % рисуем
axis equal % одинаковый масштаб осей
box on % ограничивающий прямоугольник
xlim([-Hcyl-dx Rext+dx]) % границы по оси Ox
ylim([-dx Rext+dx]) % границы по оси Oy
xlabel('itzrm, м') % метка оси Ox
ylabel('itrrm, м') % метка оси Oy
title('Область решения') % заголовок
set(get(gcf,'CurrentAxes'),'FontSize',14,...
'FontName','Times New Roman') % поменяли шрифт
print('ShellArea','-dpng') % сохранили рисунок
% узнаём и печатаем номера граничных линий
mmgg = HeatModel.Geometry.geom; % матрица граничный линий
% для правильного задания граничных условий нужны номера линий:
% 1) на наружной поверхности цилиндрической части;
% 2) на наружной поверхности сферической части;
% 3) все остальные
NumAll = [1:size(mmgg,2)]; % номера всех границ
NumUpCyl = find(all(abs(mmgg(4:5,:)-Rext)<MyEps));
% на наружной поверхности цилиндрической части номера тех
% столбцов, для которых обе y-е координаты (т.е. 4-й и 5-й
% элементы столбца) равны Rext
NumUpSph = find((abs(mmgg(1,:)-1)<MyEps)&...
(abs(mmgg(10,:)-Rext)<MyEps)); % на наружной поверхности
% сферической части номера тех столбцов, для которых 1-я
% координата равна 1 (окружность) и 10-я координата =Rext
NumOther = setdiff(setdiff(NumAll,NumUpCyl),NumUpSph);
% номера остальных линий
disp('Номера граничных линий:')
fprintf('на верхней стороне цилиндрической части:')
fprintf(' %d',NumUpCyl)
fprintf('nна верхней стороне сферической части:')
fprintf(' %d',NumUpSph)
fprintf('nна остальных сторонах:')
fprintf(' %d',NumOther)
fprintf('n')
generateMesh(HeatModel,'Hmax',(Rext-Rint)/5,...
'GeometricOrder','linear');
% создаём сетку МКЭ из 3-узловых элементов с линейным изменением
% температуры и максимальным размером элемента (Rext-Rint)/5,
% чтобы по толщине было 5 конечных элементов
disp('Создана сетка МКЭ:')
fprintf('Количество узлов Nodes = %dn',...
size(HeatModel.Mesh.Nodes,2))
fprintf('Количество элементов Elements = %dn',...
size(HeatModel.Mesh.Elements,2))
% рисуем сетку МКЭ
figure('Position',[40 40 1200 500]) % новое окно фигуры
pdemesh(HeatModel) % рисуем сетку МКЭ
axis equal % одинаковый масштаб осей
box on % ограничивающий прямоугольник
xlim([-Hcyl-dx Rext+dx]) % границы по оси Ox
ylim([-dx Rext+dx]) % границы по оси Oy
xlabel('itzrm, м') % метка оси Ox
ylabel('itrrm, м') % метка оси Oy
title('Сетка МКЭ') % заголовок
set(get(gcf,'CurrentAxes'),'FontSize',14,...
'FontName','Times New Roman') % поменяли шрифт
print('ShellGrid','-dpng') % сохранили рисунок
% для стационарной осесимметричной задачи
% в цилиндрических координатах коэффициент теплопроводности
% должен быть умножен на координату y = r:
% https://www.mathworks.com/help/pde/examples/heat-distribution-in-a-circular-cylindrical-rod.html
kFunc = @(region,state) Kti*region.y;
% задаём функцию граничного условия на внешней сфере
bSph = @(region,state) Ucyl+kUsph*region.x;
thermalProperties(HeatModel,'ThermalConductivity',kFunc);
% добавили в модель коэффициент теплопроводности
% задаём граничные условия
thermalBC(HeatModel,'Edge',NumUpCyl,...
'Temperature',Ucyl); % на верхней поверхности цилиндра
thermalBC(HeatModel,'Edge',NumUpSph,...
'Temperature',bSph); % на верхней поверхности сферы
thermalBC(HeatModel,'Edge',NumOther,...
'HeatFlux',0); % на остальных границах, в т.ч. и на оси
% симметрии - теплоизоляция
HeatRes = solve(HeatModel); % решаем задачу
U = HeatRes.Temperature; % стационарное распределение температур
Umax = max(U); % max температура
Umin = min(U); % min температура
disp('Результат решения:')
fprintf('Максимальная температура Umax = %8.4f Кn',Umax)
fprintf('Минимальная температура Umin = %8.4f Кn',Umin)
% рисуем распределение температур в виде сверху
figure('Position',[40 40 1300 500]) % новое окно фигуры
pdeplot(HeatModel,'xydata',U,...
'mesh','off','colorbar','off') % рисуем температуры
view(2) % вид сверху
hold on % добавляем на рисунок другой рисунок
pdegplot(HeatModel) % добавили саму модель
caxis([Umin Umax]) % границы вдоль Oz
hold off % больше добавлять другие рисунки не будем
axis equal % одинаковый масштаб осей Ox, Oy
box on % ограничивающий прямоугольник
colormap('jet') % яркая палитра цветов
h = colorbar; % дескриптор шкалы
h.Limits = [Umin Umax]; % границы шкалы
colorbar % нарисовали шкалу
xlim([-Hcyl-dx Rext+dx]) % границы по оси Ox
ylim([-dx Rext+dx]) % границы по оси Oy
xlabel('itzrm, м') % метка оси Ox
ylabel('itrrm, м') % метка оси Oy
title('Решение стационарной задачи теплопроводности') % заголовок
set(get(gcf,'CurrentAxes'),'FontSize',14,...
'FontName','Times New Roman') % поменяли шрифт
print('ShellSol2D','-dpng') % сохранили рисунок
% рисуем распределение температур в аксонометрии
figure('Position',[40 40 1200 800]) % новое окно фигуры
p = HeatModel.Mesh.Nodes; % узлы
t = [HeatModel.Mesh.Elements;...
ones(1,size(HeatModel.Mesh.Elements,2))];
% элементы с добавкой строки единиц
pdesurf(p,t,U) % рисуем температуры
view(3) % вид в аксонометрии
ax = gca; % дескриптор текущего рисунка
ax.BoxStyle = 'full'; % стиль окаймляющего параллелепипеда
box on % нарисовали окаймляющий параллелепипед
colormap('jet') % яркая палитра цветов
xlim([-Hcyl-dx Rext+dx]) % границы по оси Ox
ylim([-dx Rext+dx]) % границы по оси Oy
zlim([Umin Umax]) % границы по оси Oz
xlabel('itzrm, м') % метка оси Ox
ylabel('itrrm, м') % метка оси Oy
zlabel('itUrm, К') % метка оси Oz
pbaspect([Rext+Hcyl+2*dx Rext+2*dx 1]) % масштабы осей
title('Решение стационарной задачи теплопроводности') % заголовок
set(get(gcf,'CurrentAxes'),'FontSize',14,...
'FontName','Times New Roman') % поменяли шрифт
print('ShellSol3D','-dpng') % сохранили рисунок
=======================================================================
При запуске этого скрипта на счёт в командном окне MATLAB печатаются результаты:
=======================================================================
Номера граничных линий:
на верхней стороне цилиндрической части: 1
на верхней стороне сферической части: 6
на остальных сторонах: 2 3 4 5 7
Создана сетка МКЭ:
Количество узлов Nodes = 1039
Количество элементов Elements = 1770
Результат решения:
Максимальная температура Umax = 500.0000 К
Минимальная температура Umin = 200.0000 К
=======================================================================
И рисуются 4 рисунка: область решения, сетка МКЭ, распределение температур в виде сверху и в аксонометрии. Если вы сами запустите этот скрипт на счёт, то последнюю фигуру (аксонометрию) можно будет покрутить мышкой и посмотреть, как меняется температура на стыке цилиндрической и сферической частей.
