Введение

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

Граничные условия на каждом краю могут быть:

  1. жёсткая заделка;
  2. шарнир;
  3. свободный край.

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

Выберем систему координат так, как показано на рис.2.

Начало координат O поместим на левом краю (в жёсткой заделке), ось Oz направим вдоль оси балки, а оси Ox и Oy − вдоль главных центральных осей инерции. Все силовые факторы считаем действующими в плоскости yOz, как показано на рис.2.

Будем использовать правило знаков плюс-плюс-плюс-плюс:

В соответствии с [1] выберем положительное направление прогиба w(z) вверх, в сторону положительного направления оси Oy (рис.3).

Тогда положительные значения углов поворота θ(z) будут соответствовать возрастанию прогиба w(z), а отрицательные − убыванию (рис.4).

Изгибающий момент − это вторая производная от прогиба (с точностью до множителя) и первая производная от угла поворота θ(z) (опять-таки с точностью до множителя); поэтому положительное значение момента M(z) соответствует увеличению угла поворота θ(z), т.е. изгибу балочки выпуклостью вниз, а отрицательный M(z) − изгибу выпуклостью вверх (рис.5).

При построении эпюр мы будем разрезать балку в данном сечении z, отбрасывать левую часть и заменять её эквивалентной системой сил и моментов. Положительное значение M(z) (выпуклостью вниз) при этом даст момент, направленный по часовой стрелке (рис.6).

Поэтому в исходных данных сосредоточенные моменты будем задавать положительными, если они направлены по часовой стрелке.

Теперь рассмотрим правило знаков для перерезывающих сил. В соответствии с (3) положительной будем считать такую силу Q(z), которая соответствует возрастанию изгибающего момента M(z) при увеличении z. Наглядно представить себе увеличение вогнутости трудно, поэтому применим другое правило для определения знака Q(z). Заменим отрезанную левую часть такой силой, которая соответствует увеличению M(z) (рис.7). Т.к. момент равен произведению силы на плечо, то положительное значение сосредоточенной силы соответствует направлению её вверх. Такая сила стремится повернуть элемент балки по часовой стрелке.

И, наконец, выведем правило знаков для распределённой нагрузки q(z). Положительная q(z) соответствует возрастанию перерезывающей силы Q(z). На рис.8 показано положительное направление q(z): вверх. Именно такое направление q(z) соответствует возрастанию Q(z).

Итак, подытожим всё вышесказанное. При задании исходных данных будем считать:

При построении эпюр будем руководствоваться формулами (1-4). Считаем:

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

Краткие теоретические сведения

В соответствии с [1] выведем дифференциальное уравнение изогнутой оси балки. Это уравнение имеет вид: вторая производная от изгибающего момента (т.е. EJxwIV(z)) равна сумме всех распределённых нагрузок в данном сечении на единицу длины. Такими нагрузками будут внешние силовые факторы и реакция упругого основания. Коэффициент жёсткости упругого основания c − это сила, с которой действует упругое основание на единицу площади нижней поверхности балки при единичном прогибе основания. Поэтому размерность c будет Н/м3. Если обозначить ширину нижней поверхности балки через b, то

− это сила, с которой действует упругое основание на единицу длины балки при единичном прогибе. Размерность α будет Н/м2. Сила реакции упругого основания на единицу длины тогда будет равна αw, и направлена в сторону, противоположную прогибу. Т.о., дифференциальное уравнение изогнутой оси балки на упругом основании имеет вид:

где q(z) − обобщённая внешняя нагрузка, учитывающая все силы и моменты. Дифференциальное уравнение (6) дополняется 4 граничными условиями: по 2 на левом и правом концах, в зависимости от условий закрепления.

Решение уравнения (6) удобно записывать через Крылова:

Эти функции обладают свойством:

Если ввести в рассмотрение приведённую длину:

то общее решение дифференциального уравнения (6) записывается в виде:

Здесь EJxw0 − прогиб в левом сечении (с точностью до множителя EJx), EJxθ0 − угол поворота левого сечения (также с точностью до множителя EJx), M0 и Q0 − изгибающий момент и перерезывающая сила в левом сечении. Эти 4 параметра (они называются начальными) находятся из граничных условий. В каждой из сумм суммирование проводится по всем силовым факторам, расположенным слева от текущего сечения.

Последовательное дифференцирование выражения (10) даёт угол поворота:

изгибающий момент:

и перерезывающую силу:

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

В данном методическом пособии можно использовать такие нагрузки:

Если в вашем вузе преподаватели задают студентам другие виды нагружения (распределённые моменты, линейную нагрузку и т.п.) − напишите мне, и мы вместе доработаем это пособие.

Исходными данными для выполнения этого ИДЗ являются длина балочки L, граничные условия, данные по сечению, физические характеристики материала балки и упругого основания, и нагрузка на неё: значения M, F, q и точки (интервалы) их приложения.

Рассмотрим пример. Пусть мы имеем двутавровую балку сортамента 30а длиной 10 м, жёстко защемлённую на левом краю и с шарниром на правом краю. Она нагружена системой моментов, распределённых и сосредоточенных нагрузок, показанных на рис.9. Жёсткость упругого основания 106 Н/м3, модуль упругости 2×105 Мпа.

Заполним исходные данные. Длину балочки задаём в переменной L (в метрах).

Граничные условия на левом и правом краях зададим в массиве из двух чисел LR. Первое число − это граничное условие на левом краю, а второе − на правом. Каждое из этих чисел должно быть равно или 1, или 2, или 3 в зависимости от того, какое граничное условие мы имеем на соответствующем краю:

  1. жёсткая заделка;
  2. шарнир;
  3. свободный край.

Места расположения промежуточных опор задаём в массиве Op (в метрах). Если промежуточных опор нет, задаём пустой массив: Op=[].

Сосредоточенные моменты задаём в двумерном массиве M размером 2× количество моментов. Каждая строка этого массива соответствует одному моменту: 1-е число − координата приложения в метрах, 2-е − величина момента в кНм. Если моментов нет, задаём пустой массив: M=[]. Не забывайте про знаки: положительные моменты направлены по часовой стрелке.

Точно так же задаём сосредоточенные силы: в двумерном массиве F размером 2× количество сосредоточенных сил. 1-е число каждой строки − координата приложения силы в метрах, 2-е − величина нагрузки в кН. Если сосредоточенных сил нет, задаём пустой массив: F=[]. Положительные силы направлены вверх.

Для распределённой нагрузки нужно задать и начало, и конец её приложения. Поэтому распределённые нагрузки задаём в двумерном массиве q размером 3× количество распределённых нагрузок. 1-е число каждой строки − это координата начала приложения нагрузки в метрах, 2-е − координата конца приложения нагрузки в метрах; и 3-е − величина распределённой нагрузки в кН/м (положительная − вверх). Если распределённых нагрузок нет, задаём пустой массив: q=[].

Далее задаём модуль упругости E в Мпа и жёсткость упругого основания c в Н/м3. Если в вашем вузе преподаватели задают сразу коэффициент α размерности Н/м2, задайте его дальше в программе сами (сейчас он вычисляется). Нам осталось задать профиль сечения. Базы данных различных профилей перепишите отсюда (zip-архив, 6kb). В этом архиве также находится функция для рисования сложного профиля, которая нам дальше пригодится. Базы данных в этом архиве − обычные текстовые файлы, процедура рисования − m-файл, т.е. также обычный текстовый файл. Распакуйте переписанные файлы в какой-либо каталог, доступный системе MATLAB. Номер профиля двутавра задаём в переменной Num. Если сортамент включает букву "а", её нужно задать, как одну десятую.

Введём эти данные и напечатаем их.

clear all
L=10; % длина, м
LR=[1 2]; % 1-жёсткое защемление, 2-свободное опирание, 3-свободный край
M=[[2,-10];[6,15]]; % изгибающие моменты, м, кНм
F=[[4,20];[6,-30]]; % сосредоточенные силы, м, кН
q=[[0,5,-20];[8,10,15]]; % распределённая нагрузка, м, кН/м
c=1e6; % жёсткость упругого основания, Н/м^3
Num=30.1; % номер двутавра
E=2e5; % модуль упругости, МПа
Gr={'жёсткое защемление','свободное опирание','свободный край'};
fprintf('Длина балочки L=%d м\n',L)
disp('Граничные условия')
fprintf('Cлева: %s\nСправа: %s\n',Gr{LR(1)},Gr{LR(2)})
fprintf('Жёсткость упругого основания c=%d Н/м^3\n',c)
fprintf('Модуль упругости E=%d МПа\n',E)
if ~isempty(M),
  disp('Изгибающие моменты')
  disp(' k      zk         Mk')
  fprintf('%2.0f    %5.2f     %6.2f\n',[(1:size(M,1));M'])
end
if ~isempty(F),
  disp('Сосредоточенные силы')
  disp(' k      zk         Fk')
  fprintf('%2.0f    %5.2f     %6.2f\n',[(1:size(F,1));F'])
end
if ~isempty(q),
  disp('Распределённая нагрузка')
  disp([' k      zk1      ' ...
    'zk2        qk'])
  fprintf('%2.0f    %5.2f    %5.2f     %6.2f\n',[(1:size(q,1));q'])
end
load dwudata.txt % загружаем базу данных двутавров
dwutavr=dwudata(find(dwudata(:,1)==Num),:); % выбираем нужный
Jx=dwutavr(9);
Wx=dwutavr(10);
Sx=dwutavr(11);
b=dwutavr(3);
if (Num==round(Num)),
  disp(['Выбираем друтавр номер ' num2str(Num)]);
else
  disp(['Выбираем друтавр номер ' num2str(round(Num)) 'a']);
end
disp('Его данные:')
fprintf('Момент инерции Jx=%d см^4\n',Jx)
fprintf('Момент cопротивления Wx=%d см^3\n',Wx)
fprintf('Статический момент Sx=%d см^3\n',Sx)
fprintf('Ширина профиля b=%d мм\n',b)
Длина балочки L=10 м
Граничные условия
Cлева: жёсткое защемление
Справа: свободное опирание
Жёсткость упругого основания c=1000000 Н/м^3
Модуль упругости E=200000 МПа
Изгибающие моменты
 k      zk         Mk
 1     2.00     -10.00
 2     6.00      15.00
Сосредоточенные силы
 k      zk         Fk
 1     4.00      20.00
 2     6.00     -30.00
Распределённая нагрузка
 k      zk1      zk2        qk
 1     0.00     5.00     -20.00
 2     8.00    10.00      15.00
Выбираем друтавр номер 30a
Его данные:
Момент инерции Jx=7780 см^4
Момент cопротивления Wx=518 см^3
Статический момент Sx=292 см^3
Ширина профиля b=145 мм

Теперь нарисуем эскиз балочки. Выбираем подходящие коэффициенты масштаба. Рисуем саму балочку, граничные условия, и все приложенные к балочке нагрузки. Для рисования стрелок удобнее всего использовать готовую функцию arrow, которую можно свободно переписать на сайте компании Mathworks или здесь (zip-архив, 17kb). Перепишите её и распакуйте в какой-либо каталог, доступный вашей системе MATLAB. Теперь запускайте на счёт эту область ввода.

s=1; % max размер по вертикали
if ~isempty(M),
  s=max([s;abs(M(:,2))]);
end
if ~isempty(F),
  s=max([s;abs(F(:,2))]);
end
if ~isempty(q),
  s=max([s;abs(q(:,3))]);
end
s002=0.02*s; % масштаб по вертикали
L002=0.02*L; % масштаб по горизонтали
figure
plot([0,L,L,0],[-s002,-s002,s002,s002],'k') % балочка
xlim([-L002,L+L002]) % пределы по осям
ylim([-1.2*s,1.2*s])
hold on
switch LR(1), % рисуем граничное условие слева
  case 1, % жёсткое защемление
    plot([0,0],[-0.1*s,0.1*s],'k') % заделка слева
    for k=0:9,
      plot([0,-L002],[(5-k)*s002,(3-k)*s002],'k');
    end
  case 2, % шарнир
    plot(0,-2*s002,'ok') % левая опора
    plot([0 -L002 +L002 0],...
      [-2*s002 -6*s002 -6*s002 -2*s002],'k')
    plot([-1.5*L002 1.5*L002],...
      [-6*s002 -6*s002],'k')
end
switch LR(2), % рисуем граничное условие справа
  case 1, % жёсткое защемление
    plot([L,L],[-0.1*s,0.1*s],'k') % заделка справа
    for k=0:9,
      plot([L,L+L002],[(5-k)*s002,(3-k)*s002],'k');
    end
  case 2, % шарнир
    plot(L,-2*s002,'ok') % правая опора
    plot([L L-L002 L+L002 L],...
      [-2*s002 -6*s002 -6*s002 -2*s002],'k')
    plot([L-1.5*L002 L+1.5*L002],...
      [-6*s002 -6*s002],'k')
end
for k=1:size(M,1), % изгибающие моменты
  x=M(k,1);
  v=M(k,2);
  plot([x-0.1*L,x,x,x+0.1*L],[-v,-v,v,v],'k');
  arrow([x-0.05*L,-v],[x-0.1*L,-v]);
  arrow([x+0.05*L,v],[x+0.1*L,v]);
  text(x-L002,abs(v)+0.1*s,['\itM_{' num2str(k) '}'])
end
for k=1:size(F,1), % поперечные силы
  x=F(k,1);
  v=F(k,2);
  plot([x,x],[-v,0],'k');
  arrow([x,-0.05*s*sign(v)],[x,0]);
  text(x-L002,-(abs(v)+0.1*s)*sign(v),['\itF_{' num2str(k) '}'])
end
for k=1:size(q,1), % распределённые нагрузки
  x1=q(k,1);
  x2=q(k,2);
  v=q(k,3);
  plot([x1,x1,x2,x2],[0,-v,-v,0],'k');
  arrow([x1,-0.05*s*sign(v)],[x1,0]);
  arrow([x2,-0.05*s*sign(v)],[x2,0]);
  text((x1+x2)/2,-(abs(v)+0.1*s)*sign(v),['\itq_{' num2str(k) '}'])
end
text(0.38*L,-s,['\itc\rm=' num2str(c) ' \itH\rm/\itм^3'])
hold off
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭскиз балочки')
xlabel('\itz\rm, м')
ylabel('\itM\rm, кНм, \itF\rm, кН, \itq\rm, кН/м')

Проверьте, правильно ли вы задали исходные данные. Если да, то идём дальше.

Нахождение начальных параметров

Наша балочка является статически неопределимой: неизвестные начальные параметры не могут быть найдены из уравнений статики. Они находятся из граничных условий. Для нахождения этих 4 неизвестных начальных параметров у нас есть столько же уравнений: по 2 каких-либо параметра на каждом краю балочки равны нулю − всего 4 уравнения.

В зависимости от вида граничных условий будут равны нулю:

  1. в жёсткой заделке − перемещение и угол поворота;
  2. при шарнирном закреплении − перемещение и изгибающий момент;
  3. на свободном краю − изгибающий момент и перерезывающая сила.

Наша задача состоит в том, чтобы сформировать эту систему уравнений и решить её. Выполним эту работу с помощью MATLAB. Вначале проделаем небольшую подготовительную работу. Найдём коэффициент α, приведённую длину балочки Leq по формулам (5) и (9), а затем − точки переключения аналитических выражений в формулах (10-13). Это будут: начало и конец балочки и точки приложения всех силовых факторов (для распределённой нагрузки q берём и начало, и конец приложения).

alpha=c*b*1e-3; % коэффициент жёсткости, Н/м^2
Leq=(4*E*Jx*1e-2/alpha)^0.25; % приведённая длина, м
fprintf('Коэффициент жёсткости alpha=c*b=%d Н/м^2\n',alpha)
fprintf(['Приведённая длина Leq>='...
  '(4*EJx/alpha)^(1/4)=%12.5f м\n'],Leq)
fprintf('Отношение L/Leq=%12.8f\n',L/Leq)
zk=unique([0;L]); % точки переключения
if ~isempty(M)
  zk=unique([zk;M(:,1)]);
end
if ~isempty(F)
  zk=unique([zk;F(:,1)]);
end
if ~isempty(q)
  zk=unique([zk;q(:,1);q(:,2)]);
end
nk=length(zk)-1; % число интервалов
fprintf('Число интервалов переключения nk=%d\n',nk)
disp('Точки переключения, м')
fprintf('%d   ',zk)
fprintf('\n')
Коэффициент жёсткости alpha=c*b=145000 Н/м^2
Приведённая длина Leq=(4*EJx/alpha)^(1/4)=     4.55172 м
Отношение L/Leq=  2.19697242
Число интервалов переключения nk=6
Точки переключения, м
0   2   4   5   6   8   10   

Составим систему уравнений. Нам интересно не только найти готовые коэффициенты, но и записать аналитические выражения для самих уравнений. Мы запишем:

  1. в переменной seq1 − исходное выражение: что мы будем приравнивать нулю;
  2. в переменной seq2 − расписанные по формулам (10-13) выражения для соответствующих величин.

Коэффициенты системы и правые части записываем в матрице A и векторе b.

Этот фрагмент программы − очень длинный, и, возможно, он будет считаться долго (около минуты).

nu=4;
fprintf('Общее число неизвестных n=%d\n',nu)
seq1=cell(nu,1); % уравнения
seq2=seq1; % развёрнутые уравнения
A=zeros(nu); % матрица системы уравнений
b=zeros(nu,1); % вектор правых частей
switch LR(1), % граничное условие слева
  case 1, % жёсткое защемление
    seq1{1}='EJx*w0';
    seq1{2}='EJx*theta0';
    A(1,1)=1;
    A(2,2)=1;
  case 2, % шарнир
    seq1{1}='EJx*w0';
    seq1{2}='M0';
    A(1,1)=1;
    A(2,3)=1;
  otherwise % свободный край
    seq1{1}='M0';
    seq1{2}='Q0';
    A(1,3)=1;
    A(2,4)=1;
end
seq2{1}=seq1{1};
seq2{2}=seq1{2};
switch LR(2), % граничное условие справа
  case 1, % жёсткое защемление
    seq1{3}=['EJx*w(' num2str(L) ')'];
    seq1{4}=['EJx*theta(' num2str(L) ')'];
    seq2{3}=['EJx*w0*Y1(L/Leq)+EJx*theta0*Leq*Y2(L/Leq)+'...
      'M0*Leq^2*Y3(L/Leq)+Q0*Leq^3*Y4(L/Leq)'];
    seq2{4}=['-4*EJx*w0/Leq*Y4(L/Leq)+EJx*theta0*Y1(L/Leq)+'...
      'M0*Leq*Y2(L/Leq)+Q0*Leq^2*Y3(L/Leq)'];
    Y1=cosh(L/Leq)*cos(L/Leq); % Y1(L/Leq)
    Y2=0.5*(cosh(L/Leq)*sin(L/Leq)+...
      sinh(L/Leq)*cos(L/Leq)); % Y2(L/Leq)
    Y3=0.5*sinh(L/Leq)*sin(L/Leq); % Y3end
    Y4=0.25*(cosh(L/Leq)*sin(L/Leq)-...
      sinh(L/Leq)*cos(L/Leq)); % Y4(L/Leq)
    A(3,1)=Y1;
    A(3,2)=Leq*Y2;
    A(3,3)=Leq^2*Y3;
    A(3,4)=Leq^3*Y4;
    A(4,1)=-4/Leq*Y4;
    A(4,2)=Y1;
    A(4,3)=Leq*Y2;
    A(4,4)=Leq^2*Y3;
    if ~isempty(M),
      for k1=1:size(M,1),
        seq2{3}=[seq2{3} '+Leq^2*M'...
          num2str(k1) '*Y3((L-'...
          num2str(M(k1,1)) ')/Leq)'];
        seq2{4}=[seq2{4} '+Leq*M' num2str(k1)...
          '*Y2((L-' num2str(M(k1,1)) ')/Leq)'];
      end
      Y2=0.5*(cosh((L-M(:,1))/Leq).*sin((L-M(:,1))/Leq)+...
        sinh((L-M(:,1))/Leq).*cos((L-M(:,1))/Leq));
      Y3=0.5*sinh((L-M(:,1))/Leq).*sin((L-M(:,1))/Leq);
      b(3)=b(3)+Leq^2*sum(M(:,2).*Y3);
      b(4)=b(4)+Leq*sum(M(:,2).*Y2);
    end
    if ~isempty(F),
      for k1=1:size(F,1),
        seq2{3}=[seq2{3} '+Leq^3*F' num2str(k1)...
          '*Y4((L-' num2str(F(k1,1)) ')/Leq)'];
        seq2{4}=[seq2{4} '+Leq^2*F' num2str(k1)...
          '*Y3((L-' num2str(F(k1,1)) ')/Leq)'];
      end
      Y3=0.5*sinh((L-F(:,1))/Leq).*sin((L-F(:,1))/Leq);
      Y4=0.25*(cosh((L-F(:,1))/Leq).*sin((L-F(:,1))/Leq)-...
        sinh((L-F(:,1))/Leq).*cos((L-F(:,1))/Leq));
      b(3)=b(3)+Leq^3*sum(F(:,2).*Y4);
      b(4)=b(4)+Leq^2*sum(F(:,2).*Y3);
    end
    if ~isempty(q),
      for k1=1:size(q,1),
        seq2{3}=[seq2{3} '-Leq^4/4*q' num2str(k1)...
          '*(Y1((L-' num2str(q(k1,1)) ')/Leq)-1)+'...
          'Leq^4/4*q' num2str(k1) '*(Y1((L-'...
          num2str(q(k1,2)) ')/Leq)-1)'];
        seq2{4}=[seq2{4} '+Leq^3*q' num2str(k1)...
          '*Y4((L-' num2str(q(k1,1)) ')/Leq)-Leq^3*q'...
          num2str(k1) '*Y4((L-' num2str(q(k1,2)) ')/Leq)'];
      end
      Y1=cosh((L-q(:,1))/Leq).*cos((L-q(:,1))/Leq);
      Y4=0.25*(cosh((L-q(:,1))/Leq).*sin((L-q(:,1))/Leq)-...
        sinh((L-q(:,1))/Leq).*cos((L-q(:,1))/Leq));
      b(3)=b(3)-Leq^4/4*sum(q(:,3).*(Y1-1));
      b(4)=b(4)+Leq^3*sum(q(:,3).*Y4);
      Y1=cosh((L-q(:,2))/Leq).*cos((L-q(:,2))/Leq);
      Y4=0.25*(cosh((L-q(:,2))/Leq).*sin((L-q(:,2))/Leq)-...
        sinh((L-q(:,2))/Leq).*cos((L-q(:,2))/Leq));
      b(3)=b(3)+Leq^4/4*sum(q(:,3).*(Y1-1));
      b(4)=b(4)-Leq^3*sum(q(:,3).*Y4);
    end
  case 2, % шарнир
    seq1{3}=['EJx*w(' num2str(L) ')'];
    seq1{4}=['M(' num2str(L) ')'];
    seq2{3}=['EJx*w0*Y1(L/Leq)+EJx*theta0*Leq*Y2(L/Leq)+'...
      'M0*Leq^2*Y3(L/Leq)+Q0*Leq^3*Y4(L/Leq)'];
    seq2{4}=['-4*EJx*w0/Leq^2*Y3(L/Leq)'...
      '-4*EJx*theta0/Leq*Y4(L/Leq)+'...
      'M0*Y1(L/Leq)+Q0*Leq*Y2(L/Leq)'];
    Y1=cosh(L/Leq)*cos(L/Leq); % Y1(L/Leq)
    Y2=0.5*(cosh(L/Leq)*sin(L/Leq)+...
      sinh(L/Leq)*cos(L/Leq)); % Y2(L/Leq)
    Y3=0.5*sinh(L/Leq)*sin(L/Leq); % Y3end
    Y4=0.25*(cosh(L/Leq)*sin(L/Leq)-...
      sinh(L/Leq)*cos(L/Leq)); % Y4(L/Leq)
    A(3,1)=Y1;
    A(3,2)=Leq*Y2;
    A(3,3)=Leq^2*Y3;
    A(3,4)=Leq^3*Y4;
    A(4,1)=-4/Leq^2*Y3;
    A(4,2)=-4/Leq*Y4;
    A(4,3)=Y1;
    A(4,4)=Leq*Y2;
    if ~isempty(M),
      for k1=1:size(M,1),
        seq2{3}=[seq2{3} '+Leq^2*M' num2str(k1) '*Y3((L-'...
          num2str(M(k1,1)) ')/Leq)'];
        seq2{4}=[seq2{4} '+M' num2str(k1) '*Y3((L-' ...
          num2str(M(k1,1)) ')/Leq)'];
      end
      Y1=cosh((L-M(:,1))/Leq).*cos((L-M(:,1))/Leq);
      Y3=0.5*sinh((L-M(:,1))/Leq).*sin((L-M(:,1))/Leq);
      b(3)=b(3)+Leq^2*sum(M(:,2).*Y3);
      b(4)=b(4)+sum(M(:,2).*Y1);
    end
    if ~isempty(F),
      for k1=1:size(F,1),
        seq2{3}=[seq2{3} '+Leq^3*F' num2str(k1) '*Y4((L-'...
          num2str(F(k1,1)) ')/Leq)'];
        seq2{4}=[seq2{4} '+Leq*F' num2str(k1) '*Y2((L-'...
          num2str(F(k1,1)) ')/Leq)'];
      end
      Y2=0.5*(cosh((L-F(:,1))/Leq).*sin((L-F(:,1))/Leq)+...
        sinh((L-F(:,1))/Leq).*cos((L-F(:,1))/Leq));
      Y4=0.25*(cosh((L-F(:,1))/Leq).*sin((L-F(:,1))/Leq)-...
        sinh((L-F(:,1))/Leq).*cos((L-F(:,1))/Leq));
      b(3)=b(3)+Leq^3*sum(F(:,2).*Y4);
      b(4)=b(4)+Leq*sum(F(:,2).*Y2);
    end
    if ~isempty(q),
      for k1=1:size(q,1),
        seq2{3}=[seq2{3} '-Leq^4/4*q' num2str(k1)...
          '*(Y1((L-' num2str(q(k1,1)) ')/Leq)-1)+Leq^4/4'...
          '*q' num2str(k1) '*(Y1((L-'...
          num2str(q(k1,2)) ')/Leq)-1)'];
        seq2{4}=[seq2{4} '+Leq^4*q'...
          num2str(k1) '*Y1((L-'...
          num2str(q(k1,1)) ')/Leq)-Leq^4*q'...
          num2str(k1) '*(Y3((L-'...
          num2str(q(k1,2)) ')/Leq)'];
      end
      Y1=cosh((L-q(:,1))/Leq).*cos((L-q(:,1))/Leq);
      Y3=0.5*sinh((L-q(:,1))/Leq).*sin((L-q(:,1))/Leq);
      b(3)=b(3)-Leq^4/4*sum(q(:,3).*(Y1-1));
      b(4)=b(4)+Leq^2*sum(q(:,3).*Y3);
      Y1=cosh((L-q(:,2))/Leq).*cos((L-q(:,2))/Leq);
      Y3=0.5*sinh((L-q(:,2))/Leq).*sin((L-q(:,2))/Leq);
      b(3)=b(3)+Leq^4/4*sum(q(:,3).*(Y1-1));
      b(4)=b(4)-Leq^2*sum(q(:,3).*Y3);
    end
  otherwise % свободный край
    seq1{3}=['M(' num2str(L) ')'];
    seq1{4}=['Q(' num2str(L) ')'];
    seq2{3}=['-4*EJx*w0/Leq^2*Y3(L/Leq)-'...
      '4*EJx*theta0/Leq*Y4(L/Leq)+'...
      'M0*Y1(L/Leq)+Q0*Leq*Y2(L/Leq)'];
    seq2{4}=['-4*EJx*w0/Leq^3*Y2(L/Leq)-'...
      '4*EJx*theta0/Leq^2*Y3(L/Leq)-'...
      '4*M0/Leq*Y4(L/Leq)+Q0*Y1(L/Leq)'];
    Y1=cosh(L/Leq)*cos(L/Leq); % Y1(L/Leq)
    Y2=0.5*(cosh(L/Leq)*sin(L/Leq)+...
      sinh(L/Leq)*cos(L/Leq)); % Y2(L/Leq)
    Y3=0.5*sinh(L/Leq)*sin(L/Leq); % Y3end
    Y4=0.25*(cosh(L/Leq)*sin(L/Leq)-...
      sinh(L/Leq)*cos(L/Leq)); % Y4(L/Leq)
    A(3,1)=-4/Leq^2*Y3;
    A(3,2)=-4/Leq*Y4;
    A(3,3)=Y1;
    A(3,4)=Leq*Y2;
    A(4,1)=-4/Leq^3*Y2;
    A(4,2)=-4/Leq^2*Y3;
    A(4,3)=-4/Leq*Y4;
    A(4,4)=Y1;
    if ~isempty(M),
      for k1=1:size(M,1),
        seq2{3}=[seq2{3} '+M' num2str(k1)...
          '*Y1((L-' num2str(M(k1,1)) ')/Leq)'];
        seq2{4}=[seq2{4} '-4/Leq*M' num2str(k1)...
          '*Y4((L-' num2str(M(k1,1)) ')/Leq)'];
      end
      Y1=cosh((L-M(:,1))/Leq).*cos((L-M(:,1))/Leq);
      Y4=0.25*(cosh((L-M(:,1))/Leq).*sin((L-M(:,1))/Leq)-...
        sinh((L-M(:,1))/Leq).*cos((L-M(:,1))/Leq));
      b(3)=b(3)+sum(M(:,2).*Y1);
      b(4)=b(4)-4/Leq*sum(M(:,2).*Y4);
    end
    if ~isempty(F),
      for k1=1:size(F,1),
        seq2{3}=[seq2{3} '+Leq*F'...
          num2str(k1) '*Y2((L-' num2str(F(k1,1)) ')/Leq)'];
        seq2{4}=[seq2{4} '+F' num2str(k1)...
          '*Y1((L-' num2str(F(k1,1)) ')/Leq)'];
      end
      Y1=cosh((L-F(:,1))/Leq).*cos((L-F(:,1))/Leq);
      Y2=0.5*(cosh((L-F(:,1))/Leq).*sin((L-F(:,1))/Leq)+...
        sinh((L-F(:,1))/Leq).*cos((L-F(:,1))/Leq));
      b(3)=b(3)+Leq*sum(F(:,2).*Y2);
      b(4)=b(4)+sum(F(:,2).*Y1);
    end
    if ~isempty(q),
      for k1=1:size(q,1),
        seq2{3}=[seq2{3} '+Leq^2*q'...
          num2str(k1) '*Y3((L-' num2str(q(k1,1)) ')/Leq)-'...
          'Leq^2*q' num2str(k1) '*Y3((L-' num2str(q(k1,2)) ')/Leq)'];
        seq2{4}=[seq2{4} '+Leq*q' num2str(k1)...
          '*Y2((L-' num2str(q(k1,1)) ...
          ')/Leq)-Leq*q' num2str(k1) '*Y2((L-'...
          num2str(q(k1,2)) ')/Leq)'];
      end
      Y2=0.5*(cosh((L-q(:,1))/Leq).*sin((L-q(:,1))/Leq)+...
        sinh((L-q(:,1))/Leq).*cos((L-q(:,1))/Leq));
      Y3=0.5*sinh((L-q(:,1))/Leq).*sin((L-q(:,1))/Leq);
      b(3)=b(3)+Leq^2*sum(q(:,3).*Y3);
      b(4)=b(4)+Leq*sum(q(:,3).*Y2);
      Y2=0.5*(cosh((L-q(:,2))/Leq).*sin((L-q(:,2))/Leq)+...
        sinh((L-q(:,2))/Leq).*cos((L-q(:,2))/Leq));
      Y3=0.5*sinh((L-q(:,2))/Leq).*sin((L-q(:,2))/Leq);
      b(3)=b(3)-Leq^2*sum(q(:,3).*Y3);
      b(4)=b(4)-Leq*sum(q(:,3).*Y2);
    end
end
disp(['Система уравнений для неизвестных начальных параметров:'])
for k=1:nu,
  disp([seq1{k} '=0'])
end
disp('Распишем её:')
for k=1:nu,
  disp([seq2{k} '=0'])
end
disp('Посчитаем коэффициенты:')
for k=1:nu,
  fprintf('%d*EJx*w0%+d*EJx*theta0%+d*M0%+d*Q0',A(k,1:4))
  fprintf('%+d=0\n',b(k))
end
Общее число неизвестных n=4
Система уравнений для неизвестных начальных параметров:
EJx*w0=0
EJx*theta0=0
EJx*w(10)=0
M(10)=0
Распишем её:
EJx*w0=0
EJx*theta0=0
EJx*w0*Y1(L/Leq)+EJx*theta0*Leq*Y2(L/Leq)+M0*Leq^2*Y3(L/Leq)+Q0*Leq^3*Y4(L/Leq)+Leq^2*M1*Y3((L-2)/Leq)+Leq^2*M2*Y3((L-6)/Leq)+Leq^3*F1*Y4((L-4)/Leq)+Leq^3*F2*Y4((L-6)/Leq)-Leq^4/4*q1*(Y1((L-0)/Leq)-1)+Leq^4/4*q1*(Y1((L-5)/Leq)-1)-Leq^4/4*q2(Y1((L-8)/Leq)-1)+Leq^4/4*q2(Y1((L-10)/Leq)-1)=0
-4*EJx*w0/Leq^2*Y3(L/Leq)-4*EJx*theta0/Leq*Y4(L/Leq)+M0*Y1(L/Leq)+Q0*Leq*Y2(L/Leq)+M1*Y3((L-2)/Leq)+M2*Y3((L-6)/Leq)+Leq*F1*Y2((L-4)/Leq)+Leq*F2*Y2((L-6)/Leq)+Leq^4*q1*Y1((L-0)/Leq)-Leq^4*q1(Y3((L-5)/Leq)+Leq^4*q2*Y1((L-8)/Leq)-Leq^4*q2*(Y3((L-10)/Leq)=0
Посчитаем коэффициенты:
1*EJx*w0+0*EJx*theta0+0*M0+0*Q0+0=0
0*EJx*w0+1*EJx*theta0+0*M0+0*Q0+0=0
-2.669130e+000*EJx*w0+2.472370e+000*EJx*theta0+3.729565e+001*M0+1.483940e+002*Q0-7.122234e+003=0
-3.475494e-001*EJx*w0-1.382849e+000*EJx*theta0-2.669130e+000*M0+2.472370e+000*Q0-4.605556e+002=0

Решаем полученную систему уравнений − находим начальные параметры EJxw0, EJxθ0, M0 и Q0.

disp('Решаем систему уравнений')
np=-A\b;
fprintf(['EJx*w0=%9.5f кНм^3\nEJx*theta0=%9.5f кНм^2\n'...
  'M0=%9.5f кНм\nQ0=%9.5f кН\n'],np(1:4))
Решаем систему уравнений
EJx*w0=  0.00000 кНм^3
EJx*theta0=  0.00000 кНм^2
M0=-103.90286 кНм
Q0= 74.10919 кН

Построение эпюр

Теперь, когда мы вычислили все необходимые данные для формул (10-13), строим эпюры перемещений, углов поворота, изгибающих моментов и перерезывающих сил. Начинаем с перемещений. Вначале записываем аналитические выражения на каждом участке, а затем строим графики перемещений. На множитель EJx мы делим.

disp('Перемещения на участках:')
z=[];
EJwz=[];
for k=1:nk,
  EJwz1=['EJx*w0*Y1(z/Leq)+EJx*theta0*Leq*Y2(z/Leq)+'...
    'M0*Leq^2*Y3(z/Leq)+Q0*Leq^3*Y4(z/Leq)'];
  EJwz2=sprintf(['%d*Y1(z/Leq)%+d*Leq*Y2(z/Leq)%+d*'...
    'Leq^2*Y3(z/Leq)%+d*Leq^3*Y4(z/Leq)'],np);
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  Y1=cosh(z0/Leq).*cos(z0/Leq); % Y1(z/Leq)
  Y2=0.5*(cosh(z0/Leq).*sin(z0/Leq)+...
    sinh(z0/Leq).*cos(z0/Leq)); % Y2(z/Leq)
  Y3=0.5*sinh(z0/Leq).*sin(z0/Leq); % Y3(z/Leq)
  Y4=0.25*(cosh(z0/Leq).*sin(z0/Leq)-...
    sinh(z0/Leq).*cos(z0/Leq)); % Y4(z/Leq)
  EJwz0=np(1)*Y1+np(2)*Leq*Y2+np(3)*Leq^2*Y3+np(4)*Leq^3*Y4;
  if ~isempty(M)
    nM=find(M(:,1)<zk(k+1)); % изгибающие моменты слева
    for k1=1:length(nM),
      Y3=0.5*sinh((z0-M(nM(k1),1))/Leq).*...
        sin((z0-M(nM(k1),1))/Leq);
      EJwz1=[EJwz1 '+Leq^2*M' num2str(nM(k1)) '*Y3((z-'...
        num2str(M(nM(k1),1)) ')/Leq)'];
      EJwz2=[EJwz2 sprintf('%+d',M(nM(k1),2))...
        '*Leq*2*Y3((z-'... num2str(M(nM(k1),1)) ')/Leq)'];
      EJwz0=EJwz0+M(nM(k1),2)*Leq^2*Y3;
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      Y4=0.25*(cosh((z0-F(nF(k1),1))/Leq).*sin((z0-F(nF(k1),1))/Leq)-...
        sinh((z0-F(nF(k1),1))/Leq).*cos((z0-F(nF(k1),1))/Leq)); % Y4(z/Leq)
      EJwz1=[EJwz1 '+Leq^3*F' num2str(nF(k1)) '*Y4((z-' ...
        num2str(F(nF(k1),1)) ')/Leq)'];
      EJwz2=[EJwz2 sprintf('%+d',F(nF(k1),2))...
        '*Leq^3*Y4((z-' num2str(F(nF(k1),1)) ')/Leq)'];
      EJwz0=EJwz0+F(nF(k1),2)*Leq^3*Y4;
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      Y1=cosh((z0-q(nq(k1),1))/Leq).*cos((z0-q(nq(k1),1))/Leq);
      EJwz1=[EJwz1 '-Leq^4/4*q' num2str(nq(k1)) ...
        '*(Y1((z-' num2str(q(nq(k1),1)) ')/Leq)-1)'];
      EJwz2=[EJwz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*Leq^4/4*(Y1((z-' num2str(q(nq(k1),1)) ')/Leq)-1)'];
      EJwz0=EJwz0-Leq^4/4*q(nq(k1),3)*(Y1-1);
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      Y1=cosh((z0-q(nq(k1),2))/Leq).*cos((z0-q(nq(k1),2))/Leq);
      EJwz1=[EJwz1 '+Leq^4/4*q' num2str(nq(k1)) ...
        '*(Y1((z-' num2str(q(nq(k1),2)) ')/Leq)-1)'];
      EJwz2=[EJwz2 sprintf('%+d',q(nq(k1),3)) ...
        '*Leq^4/4*(Y1((z-' num2str(q(nq(k1),2)) ')/Leq)-1)'];
      EJwz0=EJwz0+Leq^4/4*q(nq(k1),3)*(Y1-1);
    end
  end
  fprintf(['x є [%d,%d): EJx*w%d(z)=%s=%s\n'],...
    zk(k),zk(k+1),k,EJwz1,EJwz2);
  z=[z z0]; % накапливаем аргументы
  EJwz=[EJwz EJwz0]; % накапливаем функцию
end
figure
plot(z,EJwz/(E*Jx)*1e8,'k')
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭпюра прогибов')
xlabel('\itz\rm, м')
ylabel('\itw\rm(\itz\rm), мм')
Перемещения на участках:
x є [0,2): EJx*w1(z)=EJx*w0*Y1(z/Leq)+EJx*theta0*Leq*Y2(z/Leq)+M0*Leq^2*Y3(z/Leq)+Q0*Leq^3*Y4(z/Leq)-Leq^4/4*q1*(Y1((z-0)/Leq)-1)=0*Y1(z/Leq)+1.333744e-013*Leq*Y2(z/Leq)-1.039029e+002*Leq^2*Y3(z/Leq)+7.410919e+001*Leq^3*Y4(z/Leq)+20*Leq^4/4*(Y1((z-0)/Leq)-1)
x є [2,4): EJx*w2(z)=EJx*w0*Y1(z/Leq)+EJx*theta0*Leq*Y2(z/Leq)+M0*Leq^2*Y3(z/Leq)+Q0*Leq^3*Y4(z/Leq)+Leq^2*M1*Y3((z-2)/Leq)-Leq^4/4*q1*(Y1((z-0)/Leq)-1)=0*Y1(z/Leq)+1.333744e-013*Leq*Y2(z/Leq)-1.039029e+002*Leq^2*Y3(z/Leq)+7.410919e+001*Leq^3*Y4(z/Leq)-10*Leq^2*Y3((z-2)/Leq)+20*Leq^4/4*(Y1((z-0)/Leq)-1)
x є [4,5): EJx*w3(z)=EJx*w0*Y1(z/Leq)+EJx*theta0*Leq*Y2(z/Leq)+M0*Leq^2*Y3(z/Leq)+Q0*Leq^3*Y4(z/Leq)+Leq^2*M1*Y3((z-2)/Leq)+Leq^3*F1*Y4((z-4)/Leq)-Leq^4/4*q1*(Y1((z-0)/Leq)-1)=0*Y1(z/Leq)+1.333744e-013*Leq*Y2(z/Leq)-1.039029e+002*Leq^2*Y3(z/Leq)+7.410919e+001*Leq^3*Y4(z/Leq)-10*Leq^2*Y3((z-2)/Leq)+20*Leq^3*Y4((z-4)/Leq)+20*Leq^4/4*(Y1((z-0)/Leq)-1)
x є [5,6): EJx*w4(z)=EJx*w0*Y1(z/Leq)+EJx*theta0*Leq*Y2(z/Leq)+M0*Leq^2*Y3(z/Leq)+Q0*Leq^3*Y4(z/Leq)+Leq^2*M1*Y3((z-2)/Leq)+Leq^3*F1*Y4((z-4)/Leq)-Leq^4/4*q1(Y1((z-0)/Leq)-1)+Leq^4/4*q1*(Y1((z-5)/Leq)-1)=0*Y1(z/Leq)+1.333744e-013*Leq*Y2(z/Leq)-1.039029e+002*Leq^2*Y3(z/Leq)+7.410919e+001*Leq^3*Y4(z/Leq)-10*Leq^2*Y3((z-2)/Leq)+20*Leq^3*Y4((z-4)/Leq)+20*Leq^4/4*(Y1((z-0)/Leq)-1)-20*Leq^4/4*(Y1((z-5)/Leq)-1)
x є [6,8): EJx*w5(z)=EJx*w0*Y1(z/Leq)+EJx*theta0*Leq*Y2(z/Leq)+M0*Leq^2*Y3(z/Leq)+Q0*Leq^3*Y4(z/Leq)+Leq^2*M1*Y3((z-2)/Leq)+Leq^2*M2*Y3((z-6)/Leq)+Leq^3*F1*Y4((z-4)/Leq)+Leq^3*F2*Y4((z-6)/Leq)-Leq^4/4*q1*(Y1((z-0)/Leq)-1)+Leq^4/4*q1*(Y1((z-5)/Leq)-1)=0*Y1(z/Leq)+1.333744e-013*Leq*Y2(z/Leq)-1.039029e+002*Leq^2*Y3(z/Leq)+7.410919e+001*Leq^3*Y4(z/Leq)-10*Leq^2*Y3((z-2)/Leq)+15*Leq^2*Y3((z-6)/Leq)+20*Leq^3*Y4((z-4)/Leq)-30*Leq^3*Y4((z-6)/Leq)+20*Leq^4/4*(Y1((z-0)/Leq)-1)-20*Leq^4/4*(Y1((z-5)/Leq)-1)
x є [8,10): EJx*w6(z)=EJx*w0*Y1(z/Leq)+EJx*theta0*Leq*Y2(z/Leq)+M0*Leq^2*Y3(z/Leq)+Q0*Leq^3*Y4(z/Leq)+Leq^2*M1*Y3((z-2)/Leq)+Leq^2*M2*Y3((z-6)/Leq)+Leq^3*F1*Y4((z-4)/Leq)+Leq^3*F2*Y4((z-6)/Leq)-Leq^4/4*q1*(Y1((z-0)/Leq)-1)-Leq^4/4*q2*(Y1((z-8)/Leq)-1)+Leq^4/4*q1*(Y1((z-5)/Leq)-1)=0*Y1(z/Leq)+1.333744e-013*Leq*Y2(z/Leq)-1.039029e+002*Leq^2*Y3(z/Leq)+7.410919e+001*Leq^3*Y4(z/Leq)-10*Leq^2*Y3((z-2)/Leq)+15*Leq^2*Y3((z-6)/Leq)+20*Leq^3*Y4((z-4)/Leq)-30*Leq^3*Y4((z-6)/Leq)+20*Leq^4/4*(Y1((z-0)/Leq)-1)-15*Leq^4/4*(Y1((z-8)/Leq)-1)-20*Leq^4/4*(Y1((z-5)/Leq)-1)

Строим эпюры углов поворота (опять разделили на множитель EJx).

disp('Углы поворота на участках:')
EJtz=[];
for k=1:nk,
  EJtz1=['-4/Leq*EJx*w0*Y4(z/Leq)+EJx*theta0*Y1(z/Leq)+'...
    'M0*Leq*Y2(z/Leq)+Q0*Leq^2*Y3(z/Leq)'];
  EJtz2=sprintf(['%d*(-4/Leq)*Y4(z/Leq)%+d*Y1(z/Leq)'...
    '%+d*Leq*Y2(z/Leq)%+d*Leq^2Y3(z/Leq)'],np);
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  Y1=cosh(z0/Leq).*cos(z0/Leq); % Y1(z/Leq)
  Y2=0.5*(cosh(z0/Leq).*sin(z0/Leq)+...
    sinh(z0/Leq).*cos(z0/Leq)); % Y2(z/Leq)
  Y3=0.5*sinh(z0/Leq).*sin(z0/Leq); % Y3(z/Leq)
  Y4=0.25*(cosh(z0/Leq).*sin(z0/Leq)-...
    sinh(z0/Leq).*cos(z0/Leq)); % Y4(z/Leq)
  EJtz0=-4/Leq*np(1)*Y4+np(2)*Y1+np(3)*Leq*Y2+np(4)*Leq^2*Y3;
  if ~isempty(M)
    nM=find(M(:,1)<zk(k+1)); % изгибающие моменты слева
    for k1=1:length(nM),
      Y2=0.5*(cosh((z0-M(nM(k1),1))/Leq).*sin((z0-M(nM(k1),1))/Leq)+...
        sinh((z0-M(nM(k1),1))/Leq).*cos((z0-M(nM(k1),1))/Leq));
      EJtz1=[EJtz1 '+Leq*M' num2str(nM(k1)) '*Y2((z-' num2str(M(nM(k1),1))...
        ')/Leq)'];
      EJtz2=[EJtz2 sprintf('%+d',M(nM(k1),2))...
        '*Leq*Y2((z-' num2str(M(nM(k1),1)) ')/Leq)'];
      EJtz0=EJtz0+M(nM(k1),2)*Leq*Y2;
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      Y3=0.5*sinh((z0-F(nF(k1),1))/Leq).*...
        sin((z0-F(nF(k1),1))/Leq);
      EJtz1=[EJtz1 '+Leq^2*F' num2str(nF(k1)) '*Y3((z-' ...
        num2str(F(nF(k1),1)) ')/Leq)'];
      EJtz2=[EJtz2 sprintf('%+d',F(nF(k1),2)) ...
        '*Leq^2*Y3((z-' num2str(F(nF(k1),1)) ')/Leq)'];
      EJtz0=EJtz0+F(nF(k1),2)*Leq^2*Y3;
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      Y4=0.25*(cosh((z0-q(nq(k1),1))/Leq).*sin((z0-q(nq(k1),1))/Leq)-...
        sinh((z0-q(nq(k1),1))/Leq).*cos((z0-q(nq(k1),1))/Leq));
      EJtz1=[EJtz1 '+Leq^3*q' num2str(nq(k1)) '*Y4((z-' ...
        num2str(q(nq(k1),1)) ')/Leq)'];
      EJtz2=[EJtz2 sprintf('%+d',q(nq(k1),3)) ...
        '*Leq^3*Y4((z-' num2str(q(nq(k1),1)) ')/Leq)'];
      EJtz0=EJtz0+Leq^3*q(nq(k1),3)*Y4;
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      Y4=0.25*(cosh((z0-q(nq(k1),2))/Leq).*sin((z0-q(nq(k1),2))/Leq)-...
        sinh((z0-q(nq(k1),2))/Leq).*cos((z0-q(nq(k1),2))/Leq));
      EJtz1=[EJtz1 '-Leq^3*q' num2str(nq(k1)) '*Y4((z-' ...
        num2str(q(nq(k1),2)) ')/Leq)'];
      EJtz2=[EJtz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*Leq^3*Y4((z-' num2str(q(nq(k1),2)) ')/Leq)'];
      EJtz0=EJtz0-Leq^3*q(nq(k1),3)*Y4;
    end
  end
  fprintf(['x є [%d,%d): EJx*theta%d(z)=%s=%s\n'],...
    zk(k),zk(k+1),k,EJtz1,EJtz2);
  EJtz=[EJtz EJtz0]; % накапливаем функцию
end
figure
plot(z,EJtz/(E*Jx)*1e5,'k')
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭпюра углов поворота')
xlabel('\itz\rm, м')
ylabel(' \theta(\itz\rm)')
Углы поворота на участках:
x є [0,2): EJx*theta1(z)=-4/Leq*EJx*w0*Y4(z/Leq)+EJx*theta0*Y1(z/Leq)+M0*Leq*Y2(z/Leq)+Q0*Leq^2*Y3(z/Leq)+Leq^3*q1*Y4((z-0)/Leq)=0*(-4/Leq)*Y4(z/Leq)+1.333744e-013*Y1(z/Leq)-1.039029e+002*Leq*Y2(z/Leq)+7.410919e+001*Leq^2*Y3(z/Leq)-20*Leq^3*Y4((z-0)/Leq)
x є [2,4): EJx*theta2(z)=-4/Leq*EJx*w0*Y4(z/Leq)+EJx*theta0*Y1(z/Leq)+M0*Leq*Y2(z/Leq)+Q0*Leq^2*Y3(z/Leq)+Leq*M1*Y2((z-2)/Leq)+Leq^3*q1*Y4((z-0)/Leq)=0*(-4/Leq)*Y4(z/Leq)+1.333744e-013*Y1(z/Leq)-1.039029e+002*Leq*Y2(z/Leq)+7.410919e+001*Leq^2*Y3(z/Leq)-10*Leq*Y2((z-2)/Leq)-20*Leq^3*Y4((z-0)/Leq)
x є [4,5): EJx*theta3(z)=-4/Leq*EJx*w0*Y4(z/Leq)+EJx*theta0*Y1(z/Leq)+M0*Leq*Y2(z/Leq)+Q0*Leq^2*Y3(z/Leq)+Leq*M1*Y2((z-2)/Leq)+Leq^2*F1*Y3((z-4)/Leq)+Leq^3*q1*Y4((z-0)/Leq)=0*(-4/Leq)*Y4(z/Leq)+1.333744e-013*Y1(z/Leq)-1.039029e+002*Leq*Y2(z/Leq)+7.410919e+001*Leq^2*Y3(z/Leq)-10*Leq*Y2((z-2)/Leq)+20*Leq^2*Y3((z-4)/Leq-20*Leq^3*Y4((z-0)/Leq)
x є [5,6): EJx*theta4(z)=-4/Leq*EJx*w0*Y4(z/Leq)+EJx*theta0*Y1(z/Leq)+M0*Leq*Y2(z/Leq)+Q0*Leq^2*Y3(z/Leq)+Leq*M1*Y2((z-2)/Leq)+Leq^2*F1*Y3((z-4)/Leq)+Leq^3*q1*Y4((z-0)/Leq)-Leq^3*q1*Y4((z-5)/Leq)=0*(-4/Leq)*Y4(z/Leq)+1.333744e-013*Y1(z/Leq)-1.039029e+002*Leq*Y2(z/Leq)+7.410919e+001*Leq^2*Y3(z/Leq)-10*Leq*Y2((z-2)/Leq)+20*Leq^2*Y3((z-4)/Leq)-20*Leq^3*Y4((z-0)/Leq)+20*Leq^3*Y4((z-5)/Leq)
x є [6,8): EJx*theta5(z)=-4/Leq*EJx*w0*Y4(z/Leq)+EJx*theta0*Y1(z/Leq)+M0*Leq*Y2(z/Leq)+Q0*Leq^2*Y3(z/Leq)+Leq*M1*Y2((z-2)/Leq)+Leq*M2*Y2((z-6)/Leq)+Leq^2*F1*Y3((z-4)/Leq)+Leq^2*F2*Y3((z-6)/Leq)+Leq^3*q1*Y4((z-0)/Leq)-Leq^3*q1*Y4((z-5)/Leq)=0*(-4/Leq)*Y4(z/Leq)+1.333744e-013*Y1(z/Leq)-1.039029e+002*Leq*Y2(z/Leq)+7.410919e+001*Leq^2*Y3(z/Leq)-10*Leq*Y2((z-2)/Leq)+15*Leq*Y2((z-6)/Leq)+20*Leq^2*Y3((z-4)/Leq)-30*Leq^2*Y3((z-6)/Leq)-20*Leq^3*Y4((z-0)/Leq)+20*Leq^3*Y4((z-5)/Leq)
x є [8,10): EJx*theta6(z)=-4/Leq*EJx*w0*Y4(z/Leq)+EJx*theta0*Y1(z/Leq)+M0*Leq*Y2(z/Leq)+Q0*Leq^2*Y3(z/Leq)+Leq*M1*Y2((z-2)/Leq)+Leq*M2*Y2((z-6)/Leq)+Leq^2*F1*Y3((z-4)/Leq)+Leq^2*F2*Y3((z-6)/Leq)+Leq^3*q1*Y4((z-0)/Leq)+Leq^3*q2*Y4((z-8)/Leq)-Leq^3*q1*Y4((z-5)/Leq)=0*(-4/Leq)*Y4(z/Leq)+1.333744e-013*Y1(z/Leq)-1.039029e+002*Leq*Y2(z/Leq)+7.410919e+001*Leq^2*Y3(z/Leq)-10*Leq*Y2((z-2)/Leq)+15*Leq*Y2((z-6)/Leq)+20*Leq^2*Y3((z-4)/Leq)-30*Leq^2*Y3((z-6)/Leq)-20*Leq^3*Y4((z-0)/Leq)+15*Leq^3*Y4((z-8)/Leq)+20*Leq^3*Y4((z-5)/Leq)

Эпюра изгибающих моментов:

disp('Изгибающие моменты на участках:')
Mz=[];
for k=1:nk,
  Mz1=['-4/Leq^2*EJx*w0*Y3(z/Leq)-4/Leq*EJx*theta0*Y4(z/Leq)+'...
    'M0*Y1(z/Leq)+Q0*Leq*Y2(z/Leq)'];
  Mz2=sprintf(['%d*(-4/Leq^2)*Y3(z/Leq)%+d*(-4/Leq)*Y4(z/Leq)'...
    '%+d*Y1(z/Leq)%+d*Leq*Y2(z/Leq)'],np);
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  Y1=cosh(z0/Leq).*cos(z0/Leq); % Y1(z/Leq)
  Y2=0.5*(cosh(z0/Leq).*sin(z0/Leq)+...
    sinh(z0/Leq).*cos(z0/Leq)); % Y2(z/Leq)
  Y3=0.5*sinh(z0/Leq).*sin(z0/Leq); % Y3(z/Leq)
  Y4=0.25*(cosh(z0/Leq).*sin(z0/Leq)-...
    sinh(z0/Leq).*cos(z0/Leq)); % Y4(z/Leq)
  Mz0=-4/Leq^2*np(1)*Y3-4/Leq*np(2)*Y4+np(3)*Y1+np(4)*Leq*Y2;
  if ~isempty(M)
    nM=find(M(:,1)<zk(k+1)); % изгибающие моменты слева
    for k1=1:length(nM),
      Y1=cosh((z0-M(nM(k1),1))/Leq).*cos((z0-M(nM(k1),1))/Leq);
      Mz1=[Mz1 '+M' num2str(nM(k1))...
        '*Y1((z-' num2str(M(nM(k1),1)) ')/Leq)'];
      Mz2=[Mz2 sprintf('%+d',M(nM(k1),2))...
        '*Y1((z-' ... num2str(M(nM(k1),1)) ')/Leq)'];
      Mz0=Mz0+M(nM(k1),2)*Y1;
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      Y2=0.5*(cosh((z0-F(nF(k1),1))/Leq).*sin((z0-F(nF(k1),1))/Leq)+...
        sinh((z0-F(nF(k1),1))/Leq).*cos((z0-F(nF(k1),1))/Leq));
      Mz1=[Mz1 '+Leq*F' num2str(nF(k1)) ...
        '*Y2((z-' ... num2str(F(nF(k1),1)) ')/Leq)'];
      Mz2=[Mz2 sprintf('%+d',F(nF(k1),2)) '*Leq*Y2((z-' ...
        num2str(F(nF(k1),1)) ')/Leq)'];
      Mz0=Mz0+F(nF(k1),2)*Leq*Y2;
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      Y3=0.5*sinh((z0-q(nq(k1),1))/Leq).*sin((z0-q(nq(k1),1))/Leq);
      Mz1=[Mz1 '+Leq^2*q' ... num2str(nq(k1)) '*Y3((z-' ...
        num2str(q(nq(k1),1)) ')/Leq)'];
      Mz2=[Mz2 sprintf('%+d',q(nq(k1),3)) ...
        '*Leq^2*Y3((z-' ... num2str(q(nq(k1),1)) ')/Leq)'];
      Mz0=Mz0+Leq^2*q(nq(k1),3)*Y3;
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      Y3=0.5*sinh((z0-q(nq(k1),2))/Leq).*sin((z0-q(nq(k1),2))/Leq);
      Mz1=[Mz1 '-Leq^2*q' ... num2str(nq(k1)) '*Y3((z-' ...
        num2str(q(nq(k1),2)) ')/Leq)'];
      Mz2=[Mz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*Leq^2*Y3((z-' ... num2str(q(nq(k1),2)) ')/Leq)'];
      Mz0=Mz0-Leq^2*q(nq(k1),3)*Y3;
    end
  end
  fprintf('x є [%d,%d): M%d(z)=%s=%s\n',zk(k),zk(k+1),k,Mz1,Mz2);
  Mz=[Mz Mz0]; % накапливаем функцию
end
figure
plot(z,Mz,'k')
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭпюра изгибающих моментов')
xlabel('\itz\rm, м')
ylabel('\itM\rm(\itz\rm), кНм')
Изгибающие моменты на участках:
x є [0,2): M1(z)=-4/Leq^2*EJx*w0*Y3(z/Leq)-4/Leq*EJx*theta0*Y4(z/Leq)+M0*Y1/Leq)+Q0*Leq*Y2(z/Leq)+Leq^2*q1*Y3((z-0)/Leq)=0*(-4/Leq^2)*Y3(z/Leq)+1.333744e-013*(-4/Leq)*Y4(z/Leq)-1.039029e+002*Y1(z/Leq)+7.410919e+001*Leq*Y2(z/Leq)-20*Leq^2*Y3((z-0)/Leq)
x є [2,4): M2(z)=-4/Leq^2*EJx*w0*Y3(z/Leq)-4/Leq*EJx*theta0*Y4(z/Leq)+M0*Y1/Leq)+Q0*Leq*Y2(z/Leq)+M1*Y1((z-2)/Leq)+Leq^2*q1*Y3((z-0)/Leq)=0*(-4/Leq^2)*Y3(z/Leq)+1.333744e-013*(-4/Leq)*Y4(z/Leq)-1.039029e+002*Y1(z/Leq)+7.410919e+001*Leq*Y2(z/Leq)-10*Y1((z-2)/Leq)-20*Leq^2*Y3((z-0)/Leq)
x є [4,5): M3(z)=-4/Leq^2*EJx*w0*Y3(z/Leq)-4/Leq*EJx*theta0*Y4(z/Leq)+M0*Y1/Leq)+Q0*Leq*Y2(z/Leq)+M1*Y1((z-2)/Leq)+Leq*F1*Y2((z-4)/Leq)+Leq^2*q1*Y3((z-0)/Leq)=0*(-4/Leq^2)*Y3(z/Leq)+1.333744e-013*(-4/Leq)*Y4(z/Leq)-1.039029e+002*Y1(z/Leq)+7.410919e+001*Leq*Y2(z/Leq)-10*Y1((z-2)/Leq)+20*Leq*Y2((z-4)/Leq)-20*Leq^2*Y3((z-0)/Leq)
x є [5,6): M4(z)=-4/Leq^2*EJx*w0*Y3(z/Leq)-4/Leq*EJx*theta0*Y4(z/Leq)+M0*Y1/Leq)+Q0*Leq*Y2(z/Leq)+M1*Y1((z-2)/Leq)+Leq*F1*Y2((z-4)/Leq)+Leq^2*q1*Y3((z-0)/Leq)-Leq^2*q1*Y3((z-5)/Leq)=0*(-4/Leq^2)*Y3(z/Leq)+1.333744e-013*(-4/Leq)*Y4(z/Leq)-1.039029e+002*Y1(z/Leq)+7.410919e+001*Leq*Y2(z/Leq)-10*Y1((z-2)/Leq)+20*Leq*Y2((z-4)/Leq)-20*Leq^2*Y3((z-0)/Leq)+20*Leq^2*Y3((z-5)/Leq)
x є [6,8): M5(z)=-4/Leq^2*EJx*w0*Y3(z/Leq)-4/Leq*EJx*theta0*Y4(z/Leq)+M0*Y1/Leq)+Q0*Leq*Y2(z/Leq)+M1*Y1((z-2)/Leq)+M2*Y1((z-6)/Leq)+Leq*F1*Y2((z-4)/Leq)+Leq*F2*Y2((z-6)/Leq)+Leq^2*q1*Y3((z-0)/Leq)-Leq^2*q1*Y3((z-5)/Leq)=0*(-4/Leq^2)*Y3(z/Leq)+1.333744e-013*(-4/Leq)*Y4(z/Leq)-1.039029e+002*Y1(z/Leq)+7.410919e+001*Leq*Y2(z/Leq)-10*Y1((z-2)/Leq)+15*Y1((z-6)/Leq)+20*Leq*Y2((z-4)/Leq)-30*Leq*Y2((z-6)/Leq)-20*Leq^2*Y3((z-0)/Leq)+20*Leq^2*Y3((z-5)/Leq)
x є [8,10): M6(z)=-4/Leq^2*EJx*w0*Y3(z/Leq)-4/Leq*EJx*theta0*Y4(z/Leq)+M0*Y1/Leq)+Q0*Leq*Y2(z/Leq)+M1*Y1((z-2)/Leq)+M2*Y1((z-6)/Leq)+Leq*F1*Y2((z-4)/Leq)+Leq*F2*Y2((z-6)/Leq)+Leq^2*q1*Y3((z-0)/Leq)+Leq^2*q2*Y3((z-8)/Leq)-Leq^2*q1*Y3((z-5)/Leq)=0*(-4/Leq^2)*Y3(z/Leq)+1.333744e-013*(-4/Leq)*Y4(z/Leq)-1.039029e+002*Y1(z/Leq)+7.410919e+001*Leq*Y2(z/Leq)-10*Y1((z-2)/Leq)+15*Y1((z-6)/Leq)+20*Leq*Y2((z-4)/Leq)-30*Leq*Y2((z-6)/Leq)-20*Leq^2*Y3((z-0)/Leq)+15*Leq^2*Y3((z-8)/Leq)+20*Leq^2*Y3((z-5)/Leq)

И, наконец, эпюра перерезывающих сил:

disp('Перерезывающие силы на участках:')
Qz=[];
for k=1:nk,
  Qz1=['-4/Leq^3*EJx*w0*Y2(z/Leq)'...
    '-4/Leq^2*EJx*theta0*Y3(z/Leq)-4/Leq'...
    '*M0*Y4(z/Leq)+Q0*Y1(z/Leq)'];
  Qz2=sprintf(['%d*(-4/Leq^3)*Y2(z/Leq)%+d'...
    '*(-4/Leq^3)*Y3(z/Leq)%+d*(-4/Leq)*Y4(z/Leq)'...
    '%+d*Y1(z/Leq)'],np);
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  Y1=cosh(z0/Leq).*cos(z0/Leq); % Y1(z/Leq)
  Y2=0.5*(cosh(z0/Leq).*sin(z0/Leq)+...
    sinh(z0/Leq).*cos(z0/Leq)); % Y2(z/Leq)
  Y3=0.5*sinh(z0/Leq).*sin(z0/Leq); % Y3(z/Leq)
  Y4=0.25*(cosh(z0/Leq).*sin(z0/Leq)-...
    sinh(z0/Leq).*cos(z0/Leq)); % Y4(z/Leq)
  Qz0=-4/Leq^3*np(1)*Y3-4/Leq^2*np(2)*Y4-4/Leq*np(3)*Y1+np(4)*Y2;
  if ~isempty(M)
    nM=find(M(:,1)<zk(k+1)); % изгибающие моменты слева
    for k1=1:length(nM),
      Y4=0.25*(cosh((z0-M(nM(k1),1))/Leq).*sin((z0-M(nM(k1),1))/Leq)-...
        sinh((z0-M(nM(k1),1))/Leq).*cos((z0-M(nM(k1),1))/Leq));
      Qz1=[Qz1 '-4/Leq*M' num2str(nM(k1))...
        '*Y4((z-' num2str(M(nM(k1),1)) ')/Leq)'];
      Qz2=[Qz2 sprintf('%+d',M(nM(k1),2)) '*(-4/Leq)'...
        '*Y4((z-' num2str(M(nM(k1),1)) ')/Leq)'];
      Qz0=Qz0-4/Leq*M(nM(k1),2)*Y4;
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      Y1=cosh((z0-F(nF(k1),1))/Leq).*cos((z0-F(nF(k1),1))/Leq);
      Qz1=[Qz1 '+F' num2str(nF(k1)) ...
        '*Y1((z-' num2str(F(nF(k1),1)) ')/Leq)'];
      Qz2=[Qz2 sprintf('%+d',F(nF(k1),2)) ...
        '*Y1((z-' num2str(F(nF(k1),1)) ')/Leq)'];
      Qz0=Qz0+F(nF(k1),2)*Y1;
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      Y2=0.5*(cosh((z0-q(nq(k1),1))/Leq).*sin((z0-q(nq(k1),1))/Leq)+...
        sinh((z0-q(nq(k1),1))/Leq).*cos((z0-q(nq(k1),1))/Leq));
      Qz1=[Qz1 '+Leq*q' num2str(nq(k1)) ...
        '*Y2((z-' num2str(q(nq(k1),1)) ')/Leq)'];
      Qz2=[Qz2 sprintf('%+d',q(nq(k1),3)) ...
        '*Leq*Y1((z-' num2str(q(nq(k1),1)) ')/Leq)'];
      Qz0=Qz0+Leq*q(nq(k1),3)*Y2;
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      Y2=0.5*(cosh((z0-q(nq(k1),2))/Leq).*sin((z0-q(nq(k1),2))/Leq)+...
        sinh((z0-q(nq(k1),2))/Leq).*cos((z0-q(nq(k1),2))/Leq));
      Qz1=[Qz1 '-Leq*q' num2str(nq(k1)) ...
        '*Y2((z-' num2str(q(nq(k1),2)) ')/Leq)'];
      Qz2=[Qz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*Y2((z-' num2str(q(nq(k1),2)) ')/Leq)'];
      Qz0=Qz0-Leq*q(nq(k1),3)*Y2;
    end
  end
  fprintf('x є [%d,%d): Q%d(z)=%s=%s\n',zk(k),zk(k+1),k,Qz1,Qz2);
  Qz=[Qz Qz0]; % накапливаем функцию
end
figure
plot(z,Qz,'k')
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭпюра перерезывающих сил')
xlabel('\itz\rm, м')
ylabel('\itQ\rm(\itz\rm), кН')
Перерезывающие силы на участках:
x є [0,2): Q1(z)=-4/Leq^3*EJx*w0*Y2(z/Leq)-4/Leq^2*EJx*theta0*Y3(z/Leq)-4/Leq*M0*Y4(z/Leq)+Q0*Y1(z/Leq)+Leq*q1*Y2((z-0)/Leq)=0*(-4/Leq^3)*Y2(z/Leq)+1.333744e-013*(-4/Leq^2)Y3(z/Leq)-1.039029e+002*(-4/Leq)*Y4(z/Leq)+7.410919e+001*Y1(z/Leq)-20*Leq*Y1((z-0)/Leq)
x є [2,4): Q2(z)=-4/Leq^3*EJx*w0*Y2(z/Leq)-4/Leq^2*EJx*theta0*Y3(z/Leq)-4/Leq*M0*Y4(z/Leq)+Q0*Y1(z/Leq)-4/Leq*M1*Y4((z-2)/Leq)+Leq*q1*Y2((z-0)/Leq)=0*(-4/Leq^3)*Y2(z/Leq)+1.333744e-013*(-4/Leq^2)Y3(z/Leq)-1.039029e+002*(-4/Leq)*Y4(z/Leq)+7.410919e+001*Y1(z/Leq)-10*(-4/Leq)*Y4((z-2)/Leq)-20*Leq*Y1((z-0)/Leq)
x є [4,5): Q3(z)=-4/Leq^3*EJx*w0*Y2(z/Leq)-4/Leq^2*EJx*theta0*Y3(z/Leq)-4/Leq*M0*Y4(z/Leq)+Q0*Y1(z/Leq)-4/Leq*M1*Y4((z-2)/Leq)+F1*Y1((z-4)/Leq)+Leq*q1*Y2((z-0)/Leq)=0*(-4/Leq^3)*Y2(z/Leq)+1.333744e-013*(-4/Leq^2)Y3(z/Leq)-1.039029e+002*(-4/Leq)*Y4(z/Leq)+7.410919e+001*Y1(z/Leq)-10*(-4/Leq)*Y4((z-2)/Leq)+20*Y1((z-4)/Leq)-20*Leq*Y1((z-0)/Leq)
x є [5,6): Q4(z)=-4/Leq^3*EJx*w0*Y2(z/Leq)-4/Leq^2*EJx*theta0*Y3(z/Leq)-4/Leq*M0*Y4(z/Leq)+Q0*Y1(z/Leq)-4/Leq*M1*Y4((z-2)/Leq)+F1*Y1((z-4)/Leq)+Leq*q1*Y2((z-0)/Leq)-Leq*q1*Y2((z-5)/Leq)=0*(-4/Leq^3)*Y2(z/Leq)+1.333744e-013*(-4/Leq^2)Y3(z/Leq)-1.039029e+002*(-4/Leq)*Y4(z/Leq)+7.410919e+001*Y1(z/Leq)-10*(-4/Leq)*Y4((z-2)/Leq)+20*Y1((z-4)/Leq)-20*Leq*Y1((z-0)/Leq)+20*Y2((z-5)/Leq)
x є [6,8): Q5(z)=-4/Leq^3*EJx*w0*Y2(z/Leq)-4/Leq^2*EJx*theta0*Y3(z/Leq)-4/Leq*M0*Y4(z/Leq)+Q0*Y1(z/Leq)-4/Leq*M1*Y4((z-2)/Leq)-4/Leq*M2*Y4((z-6)/Leq)+F1*Y1((z-4)/Leq)+F2*Y1((z-6)/Leq)+Leq*q1*Y2((z-0)/Leq)-Leq*q1*Y2((z-5)/Leq)=0*(-4/Leq^3)*Y2(z/Leq)+1.333744e-013*(-4/Leq^2)Y3(z/Leq)-1.039029e+002*(-4/Leq)*Y4(z/Leq)+7.410919e+001*Y1(z/Leq)-10*(-4/Leq)*Y4((z-2)/Leq)+15*(-4/Leq)*Y4((z-6)/Leq)+20*Y1((z-4)/Leq)-30*Y1((z-6)/Leq)-20*Leq*Y1((z-0)/Leq)+20*Y2((z-5)/Leq)
x є [8,10): Q6(z)=-4/Leq^3*EJx*w0*Y2(z/Leq)-4/Leq^2*EJx*theta0*Y3(z/Leq)-4/Leq*M0*Y4(z/Leq)+Q0*Y1(z/Leq)-4/Leq*M1*Y4((z-2)/Leq)-4/Leq*M2*Y4((z-6)/Leq)+F1*Y1((z-4)/Leq)+F2*Y1((z-6)/Leq)+Leq*q1*Y2((z-0)/Leq)+Leq*q2*Y2((z-8)/Leq)-Leq*q1*Y2((z-5)/Leq)=0*(-4/Leq^3)*Y2(z/Leq)+1.333744e-013*(-4/Leq^2)Y3(z/Leq)-1.039029e+002*(-4/Leq)*Y4(z/Leq)+7.410919e+001*Y1(z/Leq)-10*(-4/Leq)*Y4((z-2)/Leq)+15*(-4/Leq)*Y4((z-6)/Leq)+20*Y1((z-4)/Leq)-30*Y1((z-6)/Leq)-20*Leq*Y1((z-0)/Leq)+15*Leq*Y1((z-8)/Leq)+20*Y2((z-5)/Leq)

Проверка по условиям прочности

Найдём максимальный (по модулю) изгибающий момент Mmax и сечение, в котором он достигается (опасное сечение), а также максимальную по модулю перерезывающую силу Qmax и сечение, в котором она достигается.

[Mmax,iMmax]=max(abs(Mz)); % максимальный изгибающий момент
fprintf(['Максимальный по модулю изгибающий момент' ...
  '|Mmax|=%7.2f кНм\n'],abs(Mmax))
fprintf('достигается при z=%6.2f м\n',z(iMmax))
[Qmax,iQmax]=max(abs(Qz)); % максимальная перерезывающая сила
fprintf(['Максимальная по модулю перерезывающая сила ' ...
  '|Qmax|=%7.2f кН\n'],abs(Qmax))
fprintf('достигается при z=%6.2f м\n',z(iQmax))
fprintf('при z=%6.2f м: |Q|=%7.2f кН\n',...
  z(iMmax),abs(Qz(iMmax)))
Максимальный по модулю изгибающий момент|Mmax|= 103.90 кНм
достигается при z=  0.00 м
Максимальная по модулю перерезывающая сила |Qmax|= 138.70 кН
достигается при z= 10.00 м
при z=  0.00 м: |Q|=  91.31 кН

Из соотношения

находим максимальное нормальное напряжение в опасном сечении. Проверим теперь касательные напряжения. В каждом сечении они подсчитываются по формуле Журавского:

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

Smax=abs(Mmax)/Wx*1000;
fprintf(['Максимальные нормальные напряжения в опасном сечении ' ...
  'SigmaMax=|Mmax|/Wx*10^3=%7.2f МПа\n'],Smax)
Tmax=abs(Qmax)*Sx/(dwutavr(4)*Jx)*100; % Tmax, МПа
TMmax=abs(Qz(iMmax))*Sx/(dwutavr(4)*Jx)*100;
disp('Касатаельные напряжения:')
fprintf('при z=%6.2f: TauMax=%7.2f МПа\n',z(iQmax),Tmax)
fprintf('при z=%6.2f: TauMax=%7.2f МПа\n',z(iMmax),TMmax)
Максимальные нормальные напряжения в опасном сечении SigmaMax=|Mmax|/Wx*10^3= 200.58 МПа
Касатаельные напряжения:
при z= 10.00: TauMax=  80.09 МПа
при z=  0.00: TauMax=  52.72 МПа

Как правило, касательные напряжения значительно меньше нормальных в одном и том же сечении, к тому же они достигаются на разных волокнах: нормальные − на крайних, а касательные − в середине сечения. Поэтому опасными является обычно нормальные напряжения.

Нарисуем распределение нормальных и касательных напряжений по сечению. Нормальные напряжения распределены линейно, а касательные − по параболе. Мы строим эпюру распределения касательных напряжений приближённо: заменяем двутавр набором прямоугольников. Вычисляем по формуле (15) напряжения в крайних волокнах тонкой вертикальной стойки (переменная t2) и во внутренних волокнах широкой горизонтальной полки (переменная t1). На стойке строим параболу, а на короткой полке ограничимся прямолинейным отрезком. Рисуем на одном рисунке сечение (процедуру для его рисования мы переписали!), распределение нормальных и касательных напряжений в опасном сечении (там, где достигается Mmax), и распределение касательных напряжений в том сечении, где достигается Qmax.

t1=Qmax*(dwutavr(3)*dwutavr(5)*(dwutavr(2)-dwutavr(5))/2)/...
  (Jx*dwutavr(3))/10; % tau на полке, МПа
t2=Qmax*(dwutavr(3)*dwutavr(5)*(dwutavr(2)-dwutavr(5))/2)/...
  (Jx*dwutavr(4))/10; % tau под полкой, МПа
a=(Tmax-t2)/(dwutavr(2)/2-dwutavr(5))^2*100; % строим параболу
y=linspace(-(dwutavr(2)/2-dwutavr(5)),dwutavr(2)/2-dwutavr(5))/10;
tau=Tmax-a*y.^2;
y=[-dwutavr(2)/20, y(1), y, y(end), dwutavr(2)/20];
tau=[0 t1 tau t1 0];
GeoDat=zeros(1,7); % массив для геометрических данных
GeoDat(1)=dwutavr(2)/10; % h
GeoDat(3:7)=dwutavr(3:7)/10; % b,s,t,R,r
PlotFig([2 Num 0 0 -GeoDat(3)/2 -GeoDat(1)/2 0],GeoDat);
hold on
plot([0,-sign(Mmax)*Smax/10,sign(Mmax)*Smax/10,0],...
  [-GeoDat(1)/2,-GeoDat(1)/2,GeoDat(1)/2,GeoDat(1)/2],'r')
plot(tau/10,y,'g')
plot(tau/10/Tmax*TMmax,y,'r')
hold off
title('\bfНапряжения\rm \sigma_{\itx}\rm, \tau_{\itx}')
xlabel(' \sigma_{\itx}\rm/10, \tau_{\itx}\rm/10, МПа')

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

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

Литература

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