Введение

В этом индивидуальном домашнем задании мы будем рассчитывать стержень сложного поперечного сечения под действием сжимающей силы, приложенной не в центре тяжести сечения (рис. 1). При таком приложении силы на стержень, кроме сжимающей силы, действуют и изгибающие моменты. Поэтому, возможно, в вашем вузе это ИДЗ называется "совместное действие изгиба и растяжения-сжатия" или как-нибудь похоже.

В этом ИДЗ обычно задаётся какой-либо сложный профиль, точка приложения силы P и допускаемые напряжения на растяжение [σ+] и сжатие [σ]. Требуется найти положение нейтральной оси, разделяющей области растянутых и сжатых волокон сечения, допускаемую силу [P] и ядро сечения. Иногда вместо допускаемых напряжений [σ+] и [σ] задаётся величина силы P, а требуется определить максимальные растягивающие и сжимающие напряжения в сечении и точки, в которых они достигаются. Мы будем решать обе эти задачи.

Данное пособие позволит вам упростить выполнение этого ИДЗ. Как и любой помощник, оно не избавляет вас от необходимости думать. Используя это пособие, вы получите техническую помощь, избавитесь от досадных ошибок вычислений, но понимать существо проблемы вы всё равно должны. Но не пугайтесь: если вы смогли найти эту страницу в Internet, то разобраться в выполнении этого задания сможете наверняка.

Мы будем выполнять это ИДЗ с помощью системы инженерных и научных расчётов MATLAB. Запустите свой MATLAB. Из любой области ввода (она отмечена вот таким зелёным шрифтом скопируйте содержимое и вставьте его в командное окно или редактор-отладчик MATLAB. Теперь запускайте на счёт. Результаты можно сравнить с содержимым области вывода (вот такой синий шрифт).

Ввод исходных данных

Исходными данными для выполнения этого ИДЗ является сечение стержня, точка приложения силы, величина силы и допускаемые напряжения на растяжение и сжатие. Для задания сечения мы будем использовать специальный инструментарий − Partial Differential Equation Toolbox (PDE). Вообще-то он предназначен для других целей (решение дифференциальных уравнений в частных производных методом конечных элементов). Но некоторые его функции носят универсальный характер и могут быть использованы для решения нашей задачи − описания сечения сложного профиля.

Область решения в PDE строится из примитивов (простейших областей). Примитивами являются круг, многоугольник, прямоугольник и эллипс. Из них можно составлять сложные области, используя логические операции + (или), * (и) и − (и не). Операции + и * имеют одинаковый приоритет, а операция − более низкий. Таким образом, чтобы задать геометрию области, нужно задать:

  1. простейшие области (примитивы),
  2. формулу для построения области из примитивов с помощью логических операций +, *, −.

Примитивы задаются в виде матрицы. Число столбцов этой матрицы равно числу примитивов, каждый столбец её содержит данные по одному примитиву. Будем обозначать эту матрицу идентификатором gd.

Рассмотрим, как задаются различные примитивы. Чтобы задать круг, нужно в соответствующий столбец матрицы gd занести числа в таком порядке: первое число − это 1 (признак круга); 2-е и 3-е числа − это x и y координаты центра; 4-е число − радиус круга. Т.о., положение круга на плоскости будет полностью определено. Аналогично задаются и другие примитивы.

Чтобы задать многоугольник, нужно в соответствующий столбец занести числа в таком порядке: первое число − это 2 (признак многоугольника); 2-е число − это n − число вершин; следующие n чисел содержат x координаты вершин; следующие n чисел − это y координаты вершин.

Прямоугольник задаётся так же, как и многоугольник, но первое число − это 3 (признак прямоугольника).

Для задания эллипса нужно в соответствующий столбец занести числа в таком порядке: первое число − 4 (признак эллипса); 2-е и 3-е числа − это x и y координаты центра; 4-е и 5-е числа − это длины полуосей; 6-е число − угол поворота эллипса (в радианах, положительный поворот − против часовой стрелки).

Если столбцы полученной матрицы gd имеют различную длину, нужно дополнить более короткие столбцы нулями.

Чтобы сформировать из этих примитивов область, нужно задать в программе формулу, с помощью которой такая область будет формироваться из столбцов матрицы gd. В этой формуле столбцы матрицы gd будут обозначены какими-либо буквами, поэтому нужно также задать соответствие между буквами в формуле и столбцами матрицы. Такое соответствие задаётся строкой, каждый символ которой соответствует одному столбцу. Можно, например, задать строку вида 'abcdefgh' (или более длинную, если это необходимо). Тогда в формуле для построения области буквой a будет обозначен 1-й столбец (примитив), буквой b2-й, и т.д.

Формирование области выполняется с помощью команды decsg. В результате работы этой процедуры возвращаются 2 массива, которые мы обозначим: dl и bt. Массив dl − это разложенная матрица геометрии (the decomposed geometry matrix). Она имеет не более 12 строк и ne столбцов, где ne − число элементарных участков границы. Элементарный участок границы − это часть какого-либо примитива с монотонным изменением x и y. Так, например, круглая или эллиптическая области имеют 4 элементарных участка границы, разделяемые самой правой, левой, верхней и нижней точками. Каждый столбец массива dl соответствует одному участку границы и содержит такие числа. 1-е число − это тип линии, которая образует данную границу: 1 − окружность, 2 − многоугольник, 3 − прямоугольник, 4 − эллипс. 2-е и 3-е числа − это начальная и конечная x координаты участка границы. 4-е и 5-е числа − это y координаты начальной и конечной точки участка границы. 6-е и 7-е числа − это номера примитивов, которые находятся слева и справа от данного участка границы. Эти номера берутся в соответствии с нумерацией примитивов в массиве gd, число 0 соответствует внешней границе. Остальные числа зависят от того, какая линия образует границу (т.е. от 1-го числа). Для окружности 8-е и 9-е числа − это координаты центра окружности, 10-е число − радиус окружности. Для прямолинейной границы ничего не задаётся (граница уже полностью определена). Для эллиптической границы 8-е и 9-е числа содержат координаты центра эллипса, 10-е и 11-е − величины полуосей, и 12-е − угол поворота. В массиве bt возвращается булевская таблица соответствия примитивов и построенной области (более подробно о bt можно узнать по команде help decsg).

Внутренние границы удаляются командой csgdel. Эта команда возвращает новые разложенную матрицу геометрии и таблицу соответствия, с удалёнными внутренними границами. С помощью команды pdegplot можно нарисовать полученную область.

Рассмотрим пример. Пусть область представляет из себя квадрат со стороной 40 см, скруглённый по верхнему краю дугой окружности радиуса 40 см, с добавкой полукруга слева и с вырезанной частью эллиптического очертания справа. Центр полукруга находится посередине левой стороны, а центр эллипса − на расстоянии 10 см наружу от середины правой стороны. Эллипс с полуосями 20 см и 10 см повёрнут на 30° против часовой стрелки. Начало координат выберем в центре полукруга (рис. 2).

Наша область образована следующими примитивами:

  1. квадрат со стороной 0.4;
  2. круг радиусом 0.4 с центром в точке (0.2; −0.2); он отрезает от квадрата верхнюю часть;
  3. круг радиусом 0.1 с центром в начале координат; он добавляется к области;
  4. эллипс с центром в точке (0.5, 0), с полуосями 0.2, 0.1, повёрнутый на угол 30°; он вырезается из области.

Зададим эти примитивы в виде матрицы gd. В 1-м столбце зададим данные для прямоугольника, во 2-м и 3-м − для кругов, и в 4-м − для эллипса. Уравняем длины всех столбцов.

Затем сформируем из этих примитивов область. Обозначим столбцы матрицы буквами a, b, c, d. В программе это − идентификатор ns. Тогда область может быть получена с помощью формулы Ω=a*b+cd. Действительно, из квадрата a кругом b вырезается область, затем к ней добавляется круг c, и из полученной области вырезается эллипс d. Формула для построения области − идентификатор fo.

Ещё нам нужно задать точку приложения силы M (в см), величину силы P в кН и допускаемые напряжения на растяжение и сжатие в Мпа (идентификаторы sdp и sdm). Если величина силы P у вас не задаётся, задайте любое ненулевое значение (можно оставить то, что есть). Сжимающая сила − отрицательная.

Введём эти данные и напечатаем их. Нарисуем на ней красным крестиком точку приложения силы. Выровняем масштабы по осям координат. Надпишем заголовок и метки осей.

clear all
gd=[2;4;0;40;40;0;-20;-20;20;20]; % многоугольник
gd=[gd,[1;20;-20;40;0;0;0;0;0;0]]; % верхняя окружность
gd=[gd,[1;0;0;10;0;0;0;0;0;0]]; % левая окружность
gd=[gd,[4;50;0;20;10;pi/6;0;0;0;0]]; % эллипс
ns='abcd'; % строка соответствия столбцов
fo='(a*b+c)-d'; % формула операций
M=[10;-10]; % точка приложения силы, см
P=-100; % величина силы, кН
sdp=30; % допускаемое напряжение на растяжение, МПа
sdm=300; % допускаемое напряжение на сжатие, МПа
fprintf('Сила приложена в точке M(%d,%d) (см)\n',M);
fprintf('Величина силы P=%d (кН)\n',P);
disp('Допускаемые напряжения:');
fprintf('На растяжение: [Sigma+]=%d (МПа)\n',sdp);
fprintf('На сжатие: [Sigma-]=%d (МПа)\n',sdm);
[dl,bt]=decsg(gd,fo,ns); % формируем область
[dl1,bt1]=csgdel(dl,bt); % удалили внутренние границы
ne=size(dl1,2); % число участков границы
figure % фигура
pdegplot(dl1) % нарисовали область
hold on
plot(M(1),M(2),'rx')
hold off
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
da=daspect; % текущие масштабы осей
da(1:2)=min(da(1:2)); % одинаковые масштабы
daspect(da); % установили одинаковые масштабы
title('\bfЭскиз сечения')
xlabel('\itx') % ось OX
ylabel('\ity') % ось OY
Сила приложена в точке M(10,-10) (см)
Величина силы P=-100 (кН)
Допускаемые напряжения:
На растяжение: [Sigma+]=30 (МПа)
На сжатие: [Sigma-]=300 (МПа)

Проверьте, правильно ли ввелись ваши исходные данные. Если да, то можете переходить к следующему пункту.

Характеристики сечения

Чтобы выяснить, насколько наше сжатие является внецентренным, нужно найти центр тяжести сечения. Далее нам нужны будут и другие характеристики сечения: площадь, моменты инерции, главные центральные оси инерции. Мы уже делали ранее это ИДЗ, но там мы брали характеристики отдельных профилей из таблиц. Здесь же у нас − сложное сечение, и таблиц данных для него нет.

Поступим следующим образом: разобьём сложное сечение на множество простых, для которых мы можем вычислить центры тяжести, площади и моменты инерции. В PDE Toolbox есть команды, позволяющие разбить сложную область на треугольники. Такая разбивка осуществляется командой initmesh, которая возвращает 3 выходных параметра: p, e и t. Смысл их следующий:

Полученную сетку можно один или несколько раз измельчить при помощи команды refinemesh. Эту команду можно использовать, например, в цикле для достижения нужной точности вычислений.

Изобразить сетку разбиения можно командой pdemesh.

Сформируем треугольную сетку и измельчим её. Напечатаем количество узлов, элементов, граничных линий. Нарисуем полученную сетку. Выравняем масштабы по осям, надпишем заголовок, подпишем координатные оси.

[p,e,t]=initmesh(dl1); % формируем FEM-сетку
[p,e,t]=refinemesh(dl1,p,e,t); % измельчаем
np=size(p,2); % число узлов
nel=size(t,2); % число элементов
nb=size(e,2); % число граничных линий
fprintf('Число узлов np=%d\n',np)
fprintf('Число треугольных элементов nel=%d\n',nel)
fprintf('Число граничных линий nb=%d\n',nb)
figure % фигура
pdemesh(p,e,t) % рисуем сетку
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
da=daspect; % текущие масштабы осей
da(1:2)=min(da(1:2)); % одинаковые масштабы
daspect(da); % установили одинаковые масштабы
title('\bfТреугольная сетка')
xlabel('\itx') % ось OX
ylabel('\ity') % ось OY
Число узлов np=514
Число треугольных элементов nel=940
Число граничных линий nb=86

Переходим к определению характеристик отдельных треугольников и всего сечения в целом. Из курса аналитической геометрии вы знаете, что центр тяжести однородного треугольника с вершинами в точках M1(x1,y1), M2(x2,y2) и M3(x3,y3) находится в точке пересечения его медиан:

а площадь может быть вычислена по формуле

причём, если вершины пронумерованы против часовой стрелки, то модуль в формуле (2) не нужен. Нетрудно также вывести формулы для вычисления моментов инерции однородного плоского треугольника относительно его центра тяжести (1). Эти формулы имеют вид:

где F − площадь (2). Для вывода формул (3-5) можно использовать, например, Symbolic Toolbox, вычислив соответствующие двойные интегралы.

Вычислим величины (1-5) для всех треугольников. Далее находим общую площадь:

где Fk − площадь k-го треугольника, а ne − общее их число. Вычисляем координаты центра тяжести фигуры:

где (xC)k, (yC)k − координаты центра тяжести k-го треугольника.

Теперь вычисляем моменты инерции всего сечения относительно общего центра тяжести по формулам:

Печатаем полученные величины.

x1=p(1,t(1,:))';
x2=p(1,t(2,:))';
x3=p(1,t(3,:))';
y1=p(2,t(1,:))';
y2=p(2,t(2,:))';
y3=p(2,t(3,:))';
xc=mean([x1 x2 x3],2);
yc=mean([y1 y2 y3],2);
F=0.5*(x2.*y3+x3.*y1+x1.*y2-x2.*y1-x1.*y3-x3.*y2);
Jx=F.*((y1-y2).^2+(y2-y3).^2+(y3-y1).^2)/36;
Jy=F.*((x1-x2).^2+(x2-x3).^2+(x3-x1).^2)/36;
Jxy=F.*((x1-x2).*(y1-y2)+(x2-x3).*(y2-y3)+...
  (x3-x1).*(y3-y1))/36;
F0=sum(F);
Xc0=sum(xc.*F)/F0;
Yc0=sum(yc.*F)/F0;
Jx0=sum(Jx+F.*(yc-Yc0).^2);
Jy0=sum(Jy+F.*(xc-Xc0).^2);
Jxy0=sum(Jxy+F.*(xc-Xc0).*(yc-Yc0));
disp('Характеристики сечения:')
disp('Центр тяжести:')
fprintf('Xc=%14.7f (см)\nYc=%14.7f (см)\n',Xc0,Yc0)
disp('Моменты инерции относительно центра тяжести:')
fprintf('Jx=%18.7f (см^4)\nJy=%18.7f (см^4)\nJxy=%18.7f (см^4)\n',...
  Jx0,Jy0,Jxy0)
Характеристики сечения:
Центр тяжести:
Xc=    16.5140457 (см)
Yc=    -0.4631797 (см)
Моменты инерции относительно центра тяжести:
Jx=    188268.3522014 (см^4)
Jy=    240485.8997257 (см^4)
Jxy=      6464.7342617 (см^4)

Определяем положение главных центральных осей инерции. Угол поворота находится из соотношения

Главные центральные моменты инерции находятся так:

При этом будет Juv = 0. Новые координаты (u,v) связаны со старыми (x,y) обычными формулами поворота и параллельного переноса:

В частности, по формулам (16) находим новые координаты точки M приложения силы P. Печатаем полученные результаты. Рисуем эскиз сечения и положение на нём главных центральных осей инерции.

if abs(Jxy0)<1e-8,
  alpha=0;
elseif abs(Jx0-Jy0)<1e-8,
  alpha=pi/2;
else
  tg2a=2*Jxy0/(Jy0-Jx0);
  alpha=atan(tg2a)/2;
end
cosa=cos(alpha);
sina=sin(alpha);
disp('Главные центральные оси инерции:')
fprintf('alpha=%5.2f (рад) = %5.2f (град)\n',...
  alpha,alpha*180/pi)
Ju=Jx0*cosa^2+Jy0*sina^2-Jxy0*2*sina*cosa;
Jv=Jy0*cosa^2+Jx0*sina^2+Jxy0*2*sina*cosa;
Juv=Jxy0*(cosa^2-sina^2)-(Jy0-Jx0)*sina*cosa;
disp('Главные центральные моменты инерции:')
fprintf('Ju=%18.7f (см^4)\nJv=%18.7f (см^4)\nJuv=%18.7f (см^4)\n',...
  Ju,Jv,Juv)
Mp=[[cosa sina];[-sina cosa]]*(M-[Xc0;Yc0]);
disp('Координаты точки приложения силы в новой системе:')
fprintf('M(%14.7f, %14.7f) (см)\n',Mp)
um1=100*cos(alpha);
um2=100*sin(alpha);
vm1=100*cos(alpha+pi/2);
vm2=100*sin(alpha+pi/2);
figure % фигура
pdegplot(dl1) % нарисовали область
hold on
plot(M(1),M(2),'rx')
plot([um1 -um1]+Xc0,[um2 -um2]+Yc0,'-m')
plot([vm1 -vm1]+Xc0,[vm2 -vm2]+Yc0,'-m')
hold off
xlim([min(p(1,:)),max(p(1,:))]);
ylim([min(p(2,:)),max(p(2,:))]);
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
da=daspect; % текущие масштабы осей
da(1:2)=min(da(1:2)); % одинаковые масштабы
daspect(da); % установили одинаковые масштабы
title('\bfГлавные центральные оси инерции')
xlabel('\itx') % ось OX
ylabel('\ity') % ось OY
Главные центральные оси инерции:
alpha= 0.12 (рад) =  6.95 (град)
Главные центральные моменты инерции:
Ju=    187479.8982916 (см^4)
Jv=    241274.3536355 (см^4)
Juv=         0.0000000 (см^4)
Координаты точки приложения силы в новой системе:
M(    -7.6207093,     -8.6780484) (см)

Нейтральная ось, опасные точки и допускаемая сила

В системе главных центральных осей инерции (uOv) напряжения в каждой точке сечения складываются из напряжений от силы P и вызываемых ею изгибающих моментов:

где uP, vP − координаты точки приложения силы, iu, iv − радиусы инерции относительно соответствующих осей. Нейтральной осью сечения называется линия, на которой напряжения равны нулю. Из (17) видно, что это − прямая. Её уравнение в системе (uOv) имеет вид:

Найдём коэффициенты этого уравнения и запишем уравнение нейтральной оси. Для рисования переведём уравнение в старую систему (xOy) по формулам (15).

Теперь находим максимальные и минимальные напряжения в сечении. У нас есть координаты всех узлов в старой системе (xOy). Находим по (16) новые координаты, затем по (17) вычисляем напряжения в каждой точке сечения, и выбираем из них максимальное и минимальное. Печатаем их и точки, где они достигаются. Рисуем эскиз нашего сечения, на котором показываем положение нейтральной оси и опасных точек. Нейтральную ось рисуем красной линией, точку с максимальным напряжением − красным кружком, а с минимальным − красным квадратом.

iu2=Ju/F0;
iv2=Jv/F0;
a=Mp(1)/iv2;
b=Mp(2)/iu2;
disp('Квадраты радиусов инерции:')
fprintf('iu^2=%18.7f (см^2)\niv^2=%18.7f (см^2)\n',iu2,iv2)
disp('Уравнение нейтральной оси:')
disp('(up/iv^2)*u+(vp/iu^2)*v=-1')
disp('Вычисляем коэффициенты:')
fprintf('%18.7f*u%+18.7f*v=-1\n',a,b)
if abs(b)<eps,
  vn=[-100 100];
  un=(-1-b*vn)/a;
else
  un=[-100 100];
  vn=(-1-a*un)/b;
end
xn=un+cosa-vn*sina+Xc0;
yn=un*sina+vn*cosa+Yc0;
u=(p(1,:)'-Xc0)*cosa+(p(2,:)'-Yc0)*sina;
v=-(p(1,:)'-Xc0)*sina+(p(2,:)'-Yc0)*cosa;
Sigma=10*P/F0*(1+a*u+b*v);
[smax,ismax]=max(Sigma);
[smin,ismin]=min(Sigma);
fprintf(['Максимальное напряжение SigmaMax=%12.7f (МПа)\nдостигается в точке '...
  'Mmax(%12.7f,%12.7f) (см)\n'],smax,p(:,ismax))
fprintf(['Минимальное напряжение SigmaMin=%12.7f (МПа)\nдостигается в точке '...
  'Mmin(%12.7f,%12.7f) (см)\n'],smin,p(:,ismin))
figure % фигура
pdegplot(dl1) % нарисовали область
hold on
plot(M(1),M(2),'rx')
plot(xn,yn,'-r')
plot(p(1,ismax),p(2,ismax),'ro')
plot(p(1,ismin),p(2,ismin),'rs')
hold off
xlim([min(p(1,:)),max(p(1,:))]);
ylim([min(p(2,:)),max(p(2,:))]);
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
da=daspect; % текущие масштабы осей
da(1:2)=min(da(1:2)); % одинаковые масштабы
daspect(da); % установили одинаковые масштабы
title('\bfНейтральная ось и опасные точки')
xlabel('\itx') % ось OX
ylabel('\ity') % ось OY
Квадраты радиусов инерции:
iu^2=       118.4589820 (см^2)
iv^2=       152.4489535 (см^2)
Уравнение нейтральной оси:
(up/iv^2)*u+(vp/iu^2)*v=-1
Вычисляем коэффициенты:
        -0.0499886*u        -0.0732578*v=-1
Максимальное напряжение SigmaMax=   0.7270292 (МПа)
достигается в точке Mmax(  38.4699445,  15.4804333) (см)
Минимальное напряжение SigmaMin=  -2.0294431 (МПа)
достигается в точке Mmin(   0.0000000, -20.0000000) (см)

Исходя из допускаемых напряжений, определяем, что более опасно: растяжение или сжатие? Находим допускаемую силу и проверяем максимальное и минимальное напряжения.

[op,iop]=max([smax/sdp, -smin/sdm]);
if iop==1,
  disp('Более опасно растяжение')
else
  disp('Более опасно сжатие')
end
Pdop=P/op;
fprintf(['Максимально допускаемая сжимающая сила '...
  '[P]=%18.7f (МПа)\n'],Pdop)
disp('При этом максимальное и минимальное напряжения:')
fprintf('SigmaMax=%18.7f (МПа)\nSigmaMin=%18.7f (МПа)\n',...
  [smax smin]/P*Pdop)
Более опасно растяжение
Максимально допускаемая сжимающая сила P=     -4126.3818464 (МПа)
При этом максимальное и минимальное напряжения:
SigmaMax=        30.0000000 (МПа)
SigmaMin=       -83.7425705 (МПа)

Ядро сечения

Ядром называется такая область внутри сечения, что, если сжимающая сила будет находиться внутри этого ядра, то напряжения во всём сечении будут одного знака (для сжимающей силы − отрицательными). Чтобы найти ядро, нужно нейтральную ось (18) проводить по касательной к сечению (точнее, к его выпуклой оболочке). При каждом положении такой касательной мы можем найти коэффициенты уравнения (18), а по ним − соответствующую точку приложения силы, которая будет находиться на линии, ограничивающей ядро.

У нас есть координаты всех узлов (массив p) и номера граничных точек (массив e). Построим по ним массив, в котором будут находиться координаты всех точек границы (последовательно одна за другой − массивы xbor и ybor), а затем найдём выпуклую оболочку полученного многоугольника (команда convhull). Координаты выпуклого многоугольника, ограничивающего наше сечение − это xk и yk. Переведём их в новую систему по формулам (16): получим uk и vk. Для каждого отрезка прямой (две последовательные точки) находим параметры прямой (18), а затем − соответствующую точку ядра. Переводим координаты точек ядра в старую систему по формулам (15). Рисуем ядро сечения.

mye=e(1:2,:); % номера узлов граничных точек
nbor=mye(:,1); % начинаем формировать последовательность
mye=mye(:,2:end); % исключили 1-ю линию
while ~isempty(mye), % пока ещё есть линии
  ii=find(nbor(end)==mye(1,:)); % номер следующей линии
  if ~isempty(ii),
    nbor=[nbor;mye(2,ii)]; % присоединили к списку
    mye=mye(:,setdiff([1:size(mye,2)],ii)); % исключили
  else
    ii=find(nbor(end)==mye(2,:)); % номер следующей линии
    if ~isempty(ii),
      nbor=[nbor;mye(1,ii)]; % присоединили к списку
      mye=mye(:,setdiff([1:size(mye,2)],ii)); % исключили
    else % есть несколько контуров
      nbor=[nbor;mye(:,1)]; % присоединили к списку 1-й
      mye=mye(:,2:end); % исключили
    end
  end
end
xbor=p(1,nbor);
ybor=p(2,nbor);
K=convhull(xbor,ybor); % выпуклая оболочка
xk=xbor(K);
yk=ybor(K);
uk=(xk-Xc0)*cosa+(yk-Yc0)*sina;
vk=-(xk-Xc0)*sina+(yk-Yc0)*cosa;
u1=uk(1:end-1);
u2=uk(2:end);
v1=vk(1:end-1);
v2=vk(2:end);
uker=iv2*(v2-v1)./(u1.*v2-u2.*v1);
vker=iu2*(u2-u1)./(u2.*v1-u1.*v2);
uker=[uker,uker(1)];
vker=[vker,vker(1)];
xker=uker*cosa-vker*sina+Xc0;
yker=uker*sina+vker*cosa+Yc0;
figure % фигура
pdegplot(dl1) % нарисовали область
hold on
plot(xker,yker,'r')
hold off
xlim([min(p(1,:)),max(p(1,:))]);
ylim([min(p(2,:)),max(p(2,:))]);
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
da=daspect; % текущие масштабы осей
da(1:2)=min(da(1:2)); % одинаковые масштабы
daspect(da); % установили одинаковые масштабы
title('\bfЯдро сечения')
xlabel('\itx') % ось OX
ylabel('\ity') % ось OY

Что делать дальше

Теперь вы можете скопировать содержимое этой страницы (через ClipBoard), перебросить его, например, в Office-Word, выбросить всё лишнее и распечатать. Осталось повторить теорию, и можно идти сдавать задание преподавателю.

Литература

  1. Справочник по сопротивлению материалов / Писаренко Г.С., Яковлев А.П., Матвеев В.В.: Отв. Ред. Писаренко Г.С. - 2-е изд., перераб. И доп. - Киев: Наукова думка, 1988. - 7736 с. - ISSN 5-12-000299-4.