Пример работы в Matlab толстые оболочки

Решаем на заказ! Пример работы в Matlab толстые оболочки

Заказать работу
Узнать стоимость

ПРИМЕРЫ РЕШЕНИЙ / ПЛОХИЕ советы

Пример нашей работы в Матлаб. Можем решить подобную для Вас. Пишите в ВК или Телеграм.
Это часть диплома. Сначала идет код, дальше 
скрины -- результаты  работы программы.
Код полностью рабочий. Если правильно скопируете -- будет работать у Вас.

 

% Стационарная задача теплопроводности в цилиндрических

% координатах в осесимметричной постановке для составной

% цилиндрическо-сферической толстой оболочки

 

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 рисунка: область решения, сетка МКЭ, распределение температур в виде сверху и в аксонометрии. Если вы сами запустите этот скрипт на счёт, то последнюю фигуру (аксонометрию) можно будет покрутить мышкой и посмотреть, как меняется температура на стыке цилиндрической и сферической частей.