Введение

В этом ИДЗ мы проведём полный расчёт балочки с любым количеством промежуточных опор (рис.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. Теперь запускайте на счёт. Результаты можно сравнить с содержимым области вывода (вот такой синий шрифт).

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

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

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

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

Рассмотрим пример. Пусть мы имеем балочку длиной 10м, жёстко защемлённую на левом краю и с шарниром на правом краю. Промежуточные опоры − в точках 2, 5 и 7 м. Она нагружена системой моментов, распределённых и сосредоточенных нагрузок, показанных на рис.9.

Заполним исходные данные. Длину балочки задаём в переменной 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=[].

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

clear all
L=10; % длина, м
LR=[1 2]; % 1-жёсткое защемление, 2-свободное опирание, 3-свободный край
Op=[2 5 7]; % места расположения опор
M=[[2,-10];[6,15]]; % изгибающие моменты, м, кНм
F=[[4,20];[6,-30]]; % сосредоточенные силы, м, кН
q=[[0,5,-20];[8,10,15]]; % распределённая нагрузка, м, кН/м
Gr={'жёсткое защемление','свободное опирание','свободный край'};
fprintf('Длина балочки L=%d м\n',L)
disp('Граничные условия')
fprintf('Cлева: %s\nСправа: %s\n',Gr{LR(1)},Gr{LR(2)})
fprintf('Число промежуточных опор: %d\n',length(Op))
if ~isempty(Op),
  disp('Положение опор, м')
  fprintf('%d   ',Op)
  fprintf('\n')
end
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
Длина балочки L=10 м
Граничные условия
Cлева: жёсткое защемление
Справа: свободное опирание
Число промежуточных опор: 3
Положение опор, м
2   5   7   
Изгибающие моменты
 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

Теперь нарисуем эскиз балочки. Выбираем подходящие коэффициенты масштаба. Рисуем саму балочку, заделку слева, и все приложенные к балочке нагрузки. Для рисования стрелок удобнее всего использовать готовую функцию 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:length(Op), % промежуточные опоры
  plot(Op(k),-2*s002,'ok') % опора
  plot([Op(k) Op(k)-L002 Op(k)+L002 Op(k)],...
    [-2*s002 -6*s002 -6*s002 -2*s002],'k')
  plot([Op(k)-1.5*L002 Op(k)+1.5*L002],...
    [-6*s002 -6*s002],'k')
  text(Op(k)-L002,0.1*s,['\itR_{' num2str(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
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, кН/м')

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

Нахождение начальных параметров и неизвестных реакций опор

Наша балочка является статически неопределимой: неизвестные реакции опор и силовые факторы на концах не могут быть найдены из уравнений статики. Применим для решения задачи метод начальных параметров. Запишем уравнение изогнутой оси балки:

Здесь EJxw0 − прогиб в левом сечении (с точностью до множителя EJx), EJxθ0 − угол поворота левого сечения (также с точностью до множителя EJx), M0 и Q0 − изгибающий момент и перерезывающая сила в левом сечении. Все эти параметры (они называются начальными) − неизвестны. В каждой из сумм суммирование проводится по всем силовым факторам, расположенным слева от текущего сечения. Во второй сумме (сосредоточенных усилий Fk) нужно также учитывать неизвестные реакции опор R1, R2, …. Т.о., в уравнении (5)n+4 неизвестных, где n − количество промежуточных опор. Если мы их найдём, то сможем построить эпюру перемещений по формуле (5) и другие эпюры с помощью производных от (5), которые по (1-3) дают:

Для нахождения этих n+4 неизвестных у нас есть столько же уравнений:

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

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

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

zk=unique([0;L;Op']); % точки переключения
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')
Число интервалов переключения nk=7
Точки переключения, м
0   2   4   5   6   7   8   10   

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

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

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

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

nu=4+length(Op); % число неизвестных и уравнений
fprintf('Общее число неизвестных n=%d\n',nu)
seq1=cell(nu,1); % уравнения
seq2=seq1; % развёрнутые уравнения
seq3=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};
seq3{1}=seq1{1};
seq3{2}=seq1{2};
switch LR(2), % граничное условие справа
  case 1, % жёсткое защемление
    seq1{3}=['EJx*w(' num2str(L) ')'];
    seq1{4}=['EJx*theta(' num2str(L) ')'];
    seq2{3}=['EJx*w0+EJx*theta0*L+'...
      'M0*L^2/2+Q0*L^3/6'];
    seq2{4}='EJx*theta0+M0*L+Q0*L^2/2';
    seq3{3}=['EJx*w0+EJx*theta0*' num2str(L) ...
      '+M0*' num2str(L) '^2/2+Q0*' num2str(L) '^3/6'];
    seq3{4}=['EJx*theta0+M0*' num2str(L) ...
      '+Q0*' num2str(L) '^2/2'];
    A(3,1)=1; 
    A(3,2)=L; 
    A(3,3)=L^2/2; 
    A(3,4)=L^3/6;
    A(4,2)=1; 
    A(4,3)=L; 
    A(4,4)=L^2/2;
    if ~isempty(Op), % есть промежуточные опоры
      for k1=1:length(Op),
        seq2{3}=[seq2{3} '+R' num2str(k1)...
          '*(L-' num2str(Op(k1)) ')^3/6'];
        seq2{4}=[seq2{4} '+R' num2str(k1)...
          '*(L-' num2str(Op(k1)) ')^2/2'];
        seq3{3}=[seq3{3} '+R' num2str(k1) ...
          '*(' num2str(L) '-' num2str(Op(k1)) ')^3/6'];
        seq3{4}=[seq3{4} '+R' num2str(k1) ...
          '*(' num2str(L) '-' num2str(Op(k1)) ')^2/2'];
        A(3,4+k1)=(L-Op(k1))^3/6;
        A(4,4+k1)=(L-Op(k1))^2/2;
      end
    end
    if ~isempty(M),
      for k1=1:size(M,1),
        seq2{3}=[seq2{3} '+M' num2str(k1)...
          '*(L-' num2str(M(k1,1)) ')^2/2'];
        seq2{4}=[seq2{4} '+M' num2str(k1)...
          '*(L-' num2str(M(k1,1)) ')'];
        seq3{3}=[seq3{3} sprintf('%+d',M(k1,2)) ...
          '*(' num2str(L) '-' num2str(M(k1,1)) ')^2/2'];
        seq3{4}=[seq3{4} sprintf('%+d',M(k1,2)) ...
          '*(' num2str(L) '-' num2str(M(k1,1)) ')'];
      end
      b(3)=b(3)+sum(M(:,2).*(L-M(:,1)).^2/2);
      b(4)=b(4)+sum(M(:,2).*(L-M(:,1)));
    end
    if ~isempty(F),
      for k1=1:size(F,1),
        seq2{3}=[seq2{3} '+F' num2str(k1)...
          '*(L-' num2str(F(k1,1)) ')^3/6'];
        seq2{4}=[seq2{4} '+F' num2str(k1)...
          '*(L-' num2str(F(k1,1)) ')^2/2'];
        seq3{3}=[seq3{3} sprintf('%+d',F(k1,2)) ...
          '*(' num2str(L) '-' num2str(F(k1,1)) ')^3/6'];
        seq3{4}=[seq3{4} sprintf('%+d',F(k1,2)) ...
          '*(' num2str(L) '-' num2str(F(k1,1)) ')^2/2'];
      end
      b(3)=b(3)+sum(F(:,2).*(L-F(:,1)).^3/6);
      b(4)=b(4)+sum(F(:,2).*(L-F(:,1)).^2/2);
    end
    if ~isempty(q),
      for k1=1:size(q,1),
        seq2{3}=[seq2{3} '+q' num2str(k1)...
          '*(L-' num2str(q(k1,1)) ')^4/24-q'...
          num2str(k1) '*(L-' num2str(q(k1,2)) ')^4/24'];
        seq2{4}=[seq2{4} '+q' num2str(k1)...
          '*(L-' num2str(q(k1,1)) ')^3/6-q'...
          num2str(k1) '*(L-' num2str(q(k1,2)) ')^3/6'];
        seq3{3}=[seq3{3} sprintf('%+d',q(k1,3)) ...
          '*(' num2str(L) '-' num2str(q(k1,1)) ')^4/24'...
          sprintf('%+d',-q(k1,3)) '*(' ...
          num2str(L) '-' num2str(q(k1,2)) ')^4/24'];
        seq3{4}=[seq3{4}  sprintf('%+d',q(k1,3)) ...
          '*(' num2str(L) '-' num2str(q(k1,1)) ')^3/6'...
          sprintf('%+d',-q(k1,3)) '*(' num2str(L) '-' ...
          num2str(q(k1,2)) ')^3/6'];
      end
      b(3)=b(3)+sum(q(:,3).*(L-q(:,1)).^4/24)-...
        sum(q(:,3).*(L-q(:,2)).^4/24);
      b(4)=b(4)+sum(q(:,3).*(L-q(:,1)).^3/6)-...
        sum(q(:,3).*(L-q(:,2)).^3/6);
    end
  case 2, % шарнир
    seq1{3}=['EJx*w(' num2str(L) ')'];
    seq1{4}=['M(' num2str(L) ')'];
    seq2{3}=['EJx*w0+EJx*theta0*L+'...
      'M0*L^2/2+Q0*L^3/6'];
    seq2{4}='M0+Q0*L';
    seq3{3}=['EJx*w0+EJx*theta0*' num2str(L) ...
      '+M0*' num2str(L) '^2/2+Q0*' num2str(L) '^3/6'];
    seq3{4}=['M0+Q0*' num2str(L)];
    A(3,1)=1; 
    A(3,2)=L; 
    A(3,3)=L^2/2; 
    A(3,4)=L^3/6;
    A(4,3)=1; 
    A(4,4)=L;
    if ~isempty(Op), % есть промежуточные опоры
      for k1=1:length(Op),
        seq2{3}=[seq2{3} '+R' num2str(k1)...
          '*(L-' num2str(Op(k1)) ')^3/6'];
        seq2{4}=[seq2{4} '+R' num2str(k1) ...
          '*(L-' num2str(Op(k1)) ')'];
        seq3{3}=[seq3{3} '+R' num2str(k1) ...
          '*(' num2str(L) '-' num2str(Op(k1)) ')^3/6'];
        seq3{4}=[seq3{4} '+R' num2str(k1) ...
          '*(' num2str(L) '-' num2str(Op(k1)) ')'];
        A(3,4+k1)=(L-Op(k1))^3/6;
        A(4,4+k1)=L-Op(k1);
      end
    end
    if ~isempty(M),
      for k1=1:size(M,1),
        seq2{3}=[seq2{3} '+M' num2str(k1) ...
          '*(L-' num2str(M(k1,1)) ')^2/2'];
        seq2{4}=[seq2{4} '+M' num2str(k1)];
        seq3{3}=[seq3{3} sprintf('%+d',M(k1,2)) ...
          '*(' num2str(L) '-' num2str(M(k1,1)) ')^2/2'];
        seq3{4}=[seq3{4} sprintf('%+d',M(k1,2))];
      end
      b(3)=b(3)+sum(M(:,2).*(L-M(:,1)).^2/2);
      b(4)=b(4)+sum(M(:,2));
    end
    if ~isempty(F),
      for k1=1:size(F,1),
        seq2{3}=[seq2{3} '+F' num2str(k1)...
          '*(L-' num2str(F(k1,1)) ')^3/6'];
        seq2{4}=[seq2{4} '+F' num2str(k1)...
          '*(L-' num2str(F(k1,1)) ')'];
        seq3{3}=[seq3{3} sprintf('%+d',F(k1,2)) ...
          '*(' num2str(L) '-' num2str(F(k1,1)) ')^3/6'];
        seq3{4}=[seq3{4} sprintf('%+d',F(k1,2)) ...
          '*(' num2str(L) '-' num2str(F(k1,1)) ')'];
      end
      b(3)=b(3)+sum(F(:,2).*(L-F(:,1)).^3/6);
      b(4)=b(4)+sum(F(:,2).*(L-F(:,1)));
    end
    if ~isempty(q),
      for k1=1:size(q,1),
        seq2{3}=[seq2{3} '+q' num2str(k1)...
          '*(L-' num2str(q(k1,1)) ')^4/24-q'...
          num2str(k1) '*(L-' num2str(q(k1,2)) ')^4/24'];
        seq2{4}=[seq2{4} '+q' num2str(k1) '*(L-' num2str(q(k1,1)) ...
          ')^2/2-q' num2str(k1) '*(L-' num2str(q(k1,2)) ')^2/2'];
        seq3{3}=[seq3{3} sprintf('%+d',q(k1,3)) '*(' num2str(L) '-' ...
          num2str(q(k1,1)) ')^4/24' sprintf('%+d',-q(k1,3)) '*(' ...
          num2str(L) '-' num2str(q(k1,2)) ')^4/24'];
        seq3{4}=[seq3{4} sprintf('%+d',q(k1,3)) ...
          '*(' num2str(L) '-' num2str(q(k1,1)) ')^2/2'...
          sprintf('%+d',-q(k1,3)) '*(' ...
          num2str(L) '-' num2str(q(k1,2)) ')^2/2'];
      end
      b(3)=b(3)+sum(q(:,3).*(L-q(:,1)).^4/24)-...
        sum(q(:,3).*(L-q(:,2)).^4/24);
      b(4)=b(4)+sum(q(:,3).*(L-q(:,1)).^2/2)-...
        sum(q(:,3).*(L-q(:,2)).^2/2);
    end
  otherwise % свободный край
    seq1{3}=['M(' num2str(L) ')'];
    seq1{4}=['Q(' num2str(L) ')'];
    seq2{3}='M0+Q0*L';
    seq2{4}='Q0';
    seq3{3}=['M0+Q0*' num2str(L)];
    seq3{4}='Q0';
    A(3,3)=1; 
    A(3,4)=L;
    A(4,4)=1;
    if ~isempty(Op), % есть промежуточные опоры
      for k1=1:length(Op),
        seq2{3}=[seq2{3} '+R' num2str(k1)...
          '*(L-' num2str(Op(k1)) ')'];
        seq2{4}=[seq2{4} '+R' num2str(k1)];
        seq3{3}=[seq3{3} '+R' num2str(k1) ...
          '*(' num2str(L) '-' num2str(Op(k1)) ')'];
        seq3{4}=[seq3{4} '+R' num2str(k1)];
        A(3,4+k1)=(L-Op(k1));
        A(4,4+k1)=1;
      end
    end
    if ~isempty(M),
      for k1=1:size(M,1),
        seq2{3}=[seq2{3} '+M' num2str(k1)];
        seq3{3}=[seq3{3} sprintf('%+d',M(k1,2))];
      end
      b(3)=b(3)+sum(M(:,2));
    end
    if ~isempty(F),
      for k1=1:size(F,1),
        seq2{3}=[seq2{3} '+F' num2str(k1)...
          '*(L-' num2str(F(k1,1)) ')'];
        seq2{4}=[seq2{4} '+F' num2str(k1)];
        seq3{3}=[seq3{3} sprintf('%+d',F(k1,2)) ...
          '*(' num2str(L) '-' num2str(F(k1,1)) ')'];
        seq3{4}=[seq3{4} sprintf('%+d',F(k1,2))];
      end
      b(3)=b(3)+sum(F(:,2).*(L-F(:,1)));
      b(4)=b(4)+sum(F(:,2));
    end
    if ~isempty(q),
      for k1=1:size(q,1),
        seq2{3}=[seq2{3} '+q' num2str(k1)...
          '*(L-' num2str(q(k1,1)) ')^2/2-q'...
          num2str(k1) '*(L-' num2str(q(k1,2)) ')^2/2'];
        seq2{4}=[seq2{4} '+q' num2str(k1)...
          '*(L-' num2str(q(k1,1)) ')-q'...
          num2str(k1) '*(L-' num2str(q(k1,2)) ')'];
        seq3{3}=[seq3{3} sprintf('%+d',q(k1,3)) ...
          '*(' num2str(L) '-' num2str(q(k1,1)) ')^2/2'...
          sprintf('%+d',-q(k1,3)) '*(' ...
          num2str(L) '-' num2str(q(k1,2)) ')^2/2'];
        seq3{4}=[seq3{4} sprintf('%+d',q(k1,3)) ...
          '*(' num2str(L) '-' num2str(q(k1,1)) ...
          ')' sprintf('%+d',-q(k1,3)) '*(' ...
          num2str(L) '-' num2str(q(k1,2)) ')'];
      end
      b(3)=b(3)+sum(q(:,3).*(L-q(:,1)).^2/2)-...
        sum(q(:,3).*(L-q(:,2)).^2/2);
      b(4)=b(4)+sum(q(:,3).*(L-q(:,1)))-...
        sum(q(:,3).*(L-q(:,2)));
    end
end
for k=1:length(Op), % есть промежуточные опоры
  seq1{k+4}=['EJx*w(' num2str(Op(k)) ')'];
  seq2{k+4}=['EJx*w0+EJx*theta0*' num2str(Op(k)) ...
    '+M0*' num2str(Op(k)) '^2/2+Q0*' ...
    num2str(Op(k)) '^3/6'];
  seq3{k+4}=seq2{k+4};
  A(k+4,1)=1; 
  A(k+4,2)=Op(k);
  A(k+4,3)=Op(k)^2/2; 
  A(k+4,4)=Op(k)^3/6;
  if k>1, % есть другие опоры левее
    for k1=1:k-1,
      seq2{k+4}=[seq2{k+4} '+R' num2str(k1) ...
        '*(' num2str(Op(k)) '-' num2str(Op(k1)) ')^3/6' ];
      seq3{k+4}=seq2{k+4};
      A(k+4,k1+4)=(Op(k)-Op(k1))^3/6;
    end
  end
  if ~isempty(M)
    nM=find(M(:,1)<Op(k)); % изгибающие моменты слева
    for k1=1:length(nM),
      seq2{k+4}=[seq2{k+4} '+M' num2str(nM(k1)) ...
        '*(' num2str(Op(k)) '-' num2str(M(nM(k1),1)) ')^2/2'];
      seq3{k+4}=[seq3{k+4} sprintf('%+d',M(nM(k1),2)) ...
        '*(' num2str(Op(k)) '-' num2str(M(nM(k1),1)) ')^2/2'];
      b(k+4)=b(k+4)+M(nM(k1),2)*(Op(k)-M(nM(k1),1))^2/2;
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<Op(k)); % сосредоточенные силы слева
    for k1=1:length(nF),
      seq2{k+4}=[seq2{k+4} '+F' num2str(nF(k1)) ...
        '*(' num2str(Op(k)) '-' num2str(F(nF(k1),1)) ')^3/6'];
      seq3{k+4}=[seq3{k+4} sprintf('%+d',F(nF(k1),2)) ...
        '*(' num2str(Op(k)) '-' num2str(F(nF(k1),1)) ')^3/6'];
      b(k+4)=b(k+4)+F(nF(k1),2)*(Op(k)-F(nF(k1),1))^3/6;
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<Op(k)); % начало q слева
    for k1=1:length(nq),
      seq2{k+4}=[seq2{k+4} '+q' num2str(nq(k1)) ...
        '*(' num2str(Op(k)) '-' num2str(q(nq(k1),1)) ')^4/24'];
      seq3{k+4}=[seq3{k+4} sprintf('%+d',q(nq(k1),3)) ...
        '*(' num2str(Op(k)) '-' num2str(q(nq(k1),1)) ')^4/24'];
      b(k+4)=b(k+4)+q(nq(k1),3)*(Op(k)-q(nq(k1),1))^4/24;
    end
    nq=find(q(:,2)<Op(k)); % конец q слева
    for k1=1:length(nq),
      seq2{k+4}=[seq2{k+4} '-q' num2str(nq(k1)) ...
        '*(' num2str(Op(k)) '-' num2str(q(nq(k1),2)) ')^4/24'];
      seq3{k+4}=[seq3{k+4} sprintf('%+d',-q(nq(k1),3)) ...
        '*(' num2str(Op(k)) '-' num2str(q(nq(k1),1)) ')^4/24'];
      b(k+4)=b(k+4)-q(nq(k1),3)*(Op(k)-q(nq(k1),2))^4/24;
    end
  end
end
disp(['Система уравнений для неизвестных '...
  'начальных параметров и реакций опор:'])
for k=1:nu,
  disp([seq1{k} '=0'])
end
disp('Распишем её:')
for k=1:nu,
  disp([seq2{k} '=0'])
end
disp('Подставим значения M, F и q:')
for k=1:nu,
  disp([seq3{k} '=0'])
end
disp('Посчитаем коэффициенты:')
for k=1:nu,
  fprintf(['%d*EJx*w0%+d*EJx*theta0'...
    '%+d*M0%+d*Q0'],A(k,1:4))
  for k1=5:nu,
    fprintf('%+d×R%d',A(k,k1),k1-4)
  end
  fprintf('%+d=0\n',b(k))
end
Общее число неизвестных n=7
Система уравнений для неизвестных начальных параметров и реакций опор:
EJx*w0=0
EJx*theta0=0
EJx*w(10)=0
M(10)=0
EJx*w(2)=0
EJx*w(5)=0
EJx*w(7)=0
Распишем её:
EJx*w0=0
EJx*theta0=0
EJx*w0+EJx*theta0*L+M0*L^2/2+Q0*L^3/6+R1*(L-2)^3/6+R2*(L-5)^3/6+R3*(L-7)^3/6+M1*(L-2)^2/2+M2*(L-6)^2/2+F1*(L-4)^3/6+F2*(L-6)^3/6+q1*(L-0)^4/24-q1*(L-5)^4/24+q2*(L-8)^4/24-q2*(L-10)^4/24=0
M0+Q0*L+R1*(L-2)+R2*(L-5)+R3*(L-7)+M1+M2+F1*(L-4)+F2*(L-6)+q1*(L-0)^2/2-q1*(L-5)^2/2+q2*(L-8)^2/2-q2*(L-10)^2/2=0
EJx*w0+EJx*theta0*2+M0*2^2/2+Q0*2^3/6+q1*(2-0)^4/24=0
EJx*w0+EJx*theta0*5+M0*5^2/2+Q0*5^3/6+R1*(5-2)^3/6+M1*(5-2)^2/2+F1*(5-4)^3/6+q1*(5-0)^4/24=0
EJx*w0+EJx*theta0*7+M0*7^2/2+Q0*7^3/6+R1*(7-2)^3/6+R2*(7-5)^3/6+M1*(7-2)^2/2+M2*(7-6)^2/2+F1*(7-4)^3/6+F2*(7-6)^3/6+q1*(7-0)^4/24-q1*(7-5)^4/24=0
Подставим значения M, F и q:
EJx*w0=0
EJx*theta0=0
EJx*w0+EJx*theta0*10+M0*10^2/2+Q0*10^3/6+R1*(10-2)^3/6+R2*(10-5)^3/6+R3*(10-7)^3/6-10*(10-2)^2/2+15%(10-6)^2/2+20*(10-4)^3/6-30*(10-6)^3/6-20*(10-0)^4/24+20*(10-5)^4/24+15*(10-8)^4/24-15*(10-10)^4/24=0
M0+Q0*10+R1*(10-2)+R2*(10-5)+R3*(10-7)-10+15+20*(10-4)-30*(10-6)-20*(10-0)^2/2+20*(10-5)^2/2+15*(10-8)^2/2-15*(10-10)^2/2=0
EJx*w0+EJx*theta0*2+M0*2^2/2+Q0*2^3/6-20*(2-0)^4/24=0
EJx*w0+EJx*theta0*5+M0*5^2/2+Q0*5^3/6+R1*(5-2)^3/6-10*(5-2)^2/2+20*(5-4)^3/6-20*(5-0)^4/24=0
EJx*w0+EJx*theta0*7+M0*7^2/2+Q0*7^3/6+R1*(7-2)^3/6+R2*(7-5)^3/6-10*(7-2)^2/2+15*(7-6)^2/2+20*(7-4)^3/6-30*(7-6)^3/6-20*(7-0)^4/24+20*(7-0)^4/24=0
Посчитаем коэффициенты:
1*EJx*w0+0*EJx*theta0+0*M0+0*Q0+0*R1+0*R2+0*R3+0=0
0*EJx*w0+1*EJx*theta0+0*M0+0*Q0+0*R1+0*R2+0*R3+0=0
1*EJx*w0+10*EJx*theta0+50*M0+1.666667e+002*Q0+8.533333e+001*R1+2.083333e+001*R2+4.500000e+000*R3-7.602500e+003=0
0*EJx*w0+0*EJx*theta0+1*M0+10*Q0+8*R1+5*R2+3*R3-715=0
1*EJx*w0+2*EJx*theta0+2*M0+1.333333e+000*Q0+0*R1+0*R2+0*R3-1.333333e+001=0
1*EJx*w0+5*EJx*theta0+1.250000e+001*M0+2.083333e+001*Q0+4.500000e+000*R1+0*R2+0*R3-5.625000e+002=0
1*EJx*w0+7*EJx*theta0+2.450000e+001*M0+5.716667e+001*Q0+2.083333e+001*R1+1.333333e+000*R2+0*R3-2020=0

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

disp('Решаем систему уравнений')
np=-A\b;
fprintf(['EJx*w0=%9.5f\nEJx*theta0=%9.5f\n'...
  'M0=%9.5f\nQ0>=%9.5f\n'],np(1:4))
for k=5:nu,
  fprintf('R%d=%9.5f\n',k-4,np(k))
end
Решаем систему уравнений
EJx*w0=  0.00000
EJx*theta0=  0.00000
M0= -8.95349
Q0= 23.43023
R1= 41.39750
R2= 28.11316
R3=  5.96845

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

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

disp('Перемещения на участках:')
z=[];
EJwz=[];
for k=1:nk,
  EJwz1='EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6';
  EJwz2=sprintf(['%d%+d*z%+d*z^2/2'...
    '%+d*z^3/6'],np(1:4));
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  EJwz0=np(1)+np(2)*z0+np(3)*z0.^2/2+np(4)*z0.^3/6; % функции для графика
  if ~isempty(Op),
    nR=find(Op<zk(k+1)); % промежуточные опоры слева
    for k1=1:length(nR),
      EJwz1=[EJwz1 '+R' num2str(nR(k1)) '*(z-' ...
        num2str(Op(nR(k1))) ')^3/6'];
      EJwz2=[EJwz2 sprintf('%+d',np(nR(k1)+4)) '*(z-' ...
        num2str(Op(nR(k1))) ')^3/6'];
      EJwz0=EJwz0+np(nR(k1)+4)*(z0-Op(nR(k1))).^3/6;
    end
  end
  if ~isempty(M)
    nM=find(M(:,1)<zk(k+1)); % изгибающие моменты слева
    for k1=1:length(nM),
      EJwz1=[EJwz1 '+M' num2str(nM(k1))...
        '*(z-' num2str(M(nM(k1),1)) ')^2/2'];
      EJwz2=[EJwz2 sprintf('%+d',M(nM(k1),2))...
        '*(z-' num2str(M(nM(k1),1)) ')^2/2'];
      EJwz0=EJwz0+M(nM(k1),2)*(z0-M(nM(k1),1)).^2/2;
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      EJwz1=[EJwz1 '+F' num2str(nF(k1)) '*(z-' ...
        num2str(F(nF(k1),1)) ')^3/6'];
      EJwz2=[EJwz2 sprintf('%+d',F(nF(k1),2)) '*(z-' ...
        num2str(F(nF(k1),1)) ')^3/6'];
      EJwz0=EJwz0+F(nF(k1),2)*(z0-F(nF(k1),1)).^3/6;
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      EJwz1=[EJwz1 '+q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),1)) ')^4/24'];
      EJwz2=[EJwz2 sprintf('%+d',q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),1)) ')^4/24'];
      EJwz0=EJwz0+q(nq(k1),3)*(z0-q(nq(k1),1)).^4/24;
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      EJwz1=[EJwz1 '-q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),2)) ')^4/24'];
      EJwz2=[EJwz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),2)) ')^4/24'];
      EJwz0=EJwz0-q(nq(k1),3)*(z0-q(nq(k1),2)).^4/24;
    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,'k')
hold on
plot(Op,zeros(size(Op)),'o')
hold off
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭпюра прогибов')
xlabel('\itz\rm, м')
ylabel('\itEJ_xw\rm(\itz\rm), кНм^3')
Перемещения на участках:
x є [0,2): EJx*w1(z)=EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6+q1*(z-0)^4/24=0+2.007283e-013*z-8.953488e+000*z^2/2+2.343023e+001*z^3/6-20*(z-0)^4/24
x є [2,4): EJx*w2(z)=EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6+R1*(z-2)^3/6+M1*(z-2)^2/2+q1*(z-0)^4/24=0+2.007283e-013*z-8.953488e+000*z^2/2+2.343023e+001*z^3/6+4.139750e+001*(z-2)^3/6-10*(z-2)^2/2-20*(z-0)^4/24
x є [4,5): EJx*w3(z)=EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6+R1*(z-2)^3/6+M1*(z-2)^2/2+F1*(z-4)^3/6+q1*(z-0)^4/24=0+2.007283e-013*z-8.953488e+000*z^2/2+2.343023e+001*z^3/6+4.139750e+001*(z-2)^3/6-10*(z-2)^2/2+20*(z-4)^3/6-20*(z-0)^4/24
x є [5,6): EJx*w4(z)=EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6+R1*(z-2)^3/6+R2*(z-5)^3/6+M1*(z-2)^2/2+F1*(z-4)^3/6+q1*(z-0)^4/24-q1*(z5)^4/24=0+2.007283e-013*z-8.953488e+000*z^2/2+2.343023e+001*z^3/6+4.139750e+001*(z-2)^3/6+2.811316e+001*(z-5)^3/6-10*(z-2)^2/2+20*(z-4)^3/6-20*(z-0)^4/24+20*(z-5)*4/24
x є [6,7): EJx*w5(z)=EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6+R1*(z-2)^3/6+R2*(z-5)^3/6+M1*(z-2)^2/2+M2*(z-6)^2/2+F1*(z-4)^3/6+F2*(z-6)^3/6+q1*(z-0)^4/24-q1*(z-5)^4/24=0+2.007283e-013*z-8.953488e+000*z^2/2+2.343023e+001*z^3/6+4.139750e+001*(z-2)^3/6+2.811316e+001*(z-5)^3/6-10*(z-2)^2/2+15*(z-6)^2/2+20*(z-4)^3/6-30*(z-6)^3/6-20*(z-0)^4/24+20*(z-5)^4/24
x є [7,8): EJx*w6(z)=EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6+R1*(z-2)^3/6+R2*(z-5)^3/6+R3*(z-7)^3/6+M1*(z-2)^2/2+M2*(z-6)^2/2+F1*(z-4)^3/6+F2*(z-6)^3/6+q1*(z-0)^4/24-q1*(z-5)^4/24=0+2.007283e-013*z-8.953488e+000*z^2/2+2.343023e+001*z^3/6+4.139750e+001*(z-2)^3/6+2.811316e+001*(z-5)^3/6+5.968454e+000*(z-7)^3/6-10*(z-2)^2/2+15*(z-6)^2/2+20*(z-4)^3/6-30*(z-6)^3/6-20*(z-0)^4/24+20*(z-5)^4/24
x є [8,10): EJx*w7(z)=EJx*w0+EJx*theta0*z+M0*z^2/2+Q0*z^3/6+R1*(z-2)^3/6+R2*(z-5)^3/6+R3*(z-7)^3/6+M1*(z-2)^2/2+M2*(z-6)^2/2+F1*(z-4)^3/6+F2*(z-6)^3/6+q1*(z-0)^4/24+q2*(z-8)^4/24-q1*(z-5)^4/24=0+2.007283e-013*z-8.953488e+000*z^2/2+2.343023e+001*z^3/6+4.139750e+001*(z-2)^3/6+2.811316e+001*(z-5)^3/6+5.968454e+000*(z-7)^3/6-10*(z-2)^2/2+15*(z-6)^2/2+20*(z-4)^3/6-30*(z-6)^3/6-20*(z-0)^4/24+15*(z-8)^4/24+20*(z-5)^4/24

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

disp('Углы поворота на участках:')
EJtz=[];
for k=1:nk,
  EJtz1=['EJx*theta0+M0*z+Q0*z^2/2'];
  EJtz2=sprintf('%d%+d*z%+d*z^2/2',np(2:4));
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  EJtz0=np(2)+np(3)*z0+np(4)*z0.^2/2; % функции для графика
  if ~isempty(Op),
    nR=find(Op<zk(k+1)); % промежуточные опоры слева
    for k1=1:length(nR),
      EJtz1=[EJtz1 '+R' num2str(nR(k1)) '*(z-' ...
        num2str(Op(nR(k1))) ')^2/2'];
      EJtz2=[EJtz2 sprintf('%+d',np(nR(k1)+4)) '*(z-' ...
        num2str(Op(nR(k1))) ')^2/2'];
      EJtz0=EJtz0+np(nR(k1)+4)*(z0-Op(nR(k1))).^2/2;
    end
  end
  if ~isempty(M)
    nM=find(M(:,1)<zk(k+1)); % изгибающие моменты слева
    for k1=1:length(nM),
      EJtz1=[EJtz1 '+M' num2str(nM(k1))...
        '*(z-' num2str(M(nM(k1),1)) ')'];
      EJtz2=[EJtz2 sprintf('%+d',M(nM(k1),2))...
        '*(z-' num2str(M(nM(k1),1)) ')'];
      EJtz0=EJtz0+M(nM(k1),2)*(z0-M(nM(k1),1));
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      EJtz1=[EJtz1 '+F' num2str(nF(k1)) '*(z-' ...
        num2str(F(nF(k1),1)) ')^2/2'];
      EJtz2=[EJtz2 sprintf('%+d',F(nF(k1),2)) '*(z-' ...
        num2str(F(nF(k1),1)) ')^2/2'];
      EJtz0=EJtz0+F(nF(k1),2)*(z0-F(nF(k1),1)).^2/2;
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      EJtz1=[EJtz1 '+q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),1)) ')^3/6'];
      EJtz2=[EJtz2 sprintf('%+d',q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),1)) ')^3/6'];
      EJtz0=EJtz0+q(nq(k1),3)*(z0-q(nq(k1),1)).^3/6;
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      EJtz1=[EJtz1 '-q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),2)) ')^3/6'];
      EJtz2=[EJtz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),2)) ')^3/6'];
      EJtz0=EJtz0-q(nq(k1),3)*(z0-q(nq(k1),2)).^3/6;
    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,'k')
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭпюра углов поворота')
xlabel('\itz\rm, м')
ylabel('\itEJ_x\rm\theta(\itz\rm), кНм^2')
Углы поворота на участках:
x є [0,2): EJx*theta1(z)=EJx*theta0+M0*z+Q0*z^2/2+q1*(z-0)^3/6=2.007283e-013-8.953488e+000*z+2.343023e+001*z^2/2-20*(z-0)^3/6
x є [2,4): EJx*theta2(z)=EJx*theta0+M0*z+Q0*z^2/2+R1*(z-2)^2/2+M1*(z-2)+q1*(z-0)^3/6=2.007283e-013-8.953488e+000*z+2.343023e+001*z^2/2+4.139750e+001*(z-2)^2/2-10*(z-2)-20*(z-0)^3/6
x є [4,5): EJx*theta3(z)=EJx*theta0+M0*z+Q0*z^2/2+R1*(z-2)^2/2+M1*(z-2)+q1*(z-0)^3/6+F1*(z-4)^2/2+q1*(z-0)^3/6=2.007283e-013-8.953488e+000*z+2.343023e+001*z^2/2+4.139750e+001*(z-2)^2/2-10*(z-2)+20*(z-4)^2/2-20*(z-0)^3/6
x є [5,6): EJx*theta4(z)=EJx*theta0+M0*z+Q0*z^2/2+R1*(z-2)^2/2+R2*(z-5)^2/2+M1*(z-2)+F1*(z-4)^2/2+q1*(z-0)^3/6-q1*(z-5)^3/6=2.007283e-013-8.953488e+000*z+2.343023e+001*z^2/2+4.139750e+001*(z-2)^2/2+2.811316e+001*(z-5)^2/2-10*(z-2)+20*(z-4)^2/2-20*(z-0)^3/6+20*(z-5)^3/6
x є [6,7): EJx*theta5(z)=EJx*theta0+M0*z+Q0*z^2/2+R1*(z-2)^2/2+R2*(z-5)^2/2+M1*(z-2)+M2*(z-6)+F1*(z-4)^2/2+F2*(z-6)^2/2+q1*(z-0)^3/6-q1*(z-5)^3/6=2.007283e-013-8.953488e+000*z+2.343023e+001*z^2/2+4.139750e+001*(z-2)^2/2+2.811316e+001*(z-5)^2/2-10*(z-2)+15*(z-6)+20*(z-4)^2/2-30*(z-6)^2/2-20*(z-0)^3/6+20*(z-5)^3/6
x є [7,8): EJx*theta6(z)=EJx*theta0+M0*z+Q0*z^2/2+R1*(z-2)^2/2+R2*(z-5)^2/2+R3*(z-7)^2/2+M1*(z-2)+M2*(z-6)+F1*(z-4)^2/2+F2*(z-6)^2/2+q1*(z-0)^3/6-q1*(z-5)^3/6=2.007283e-013-8.953488e+000*z+2.343023e+001*z^2/2+4.139750e+001*(z-2)^2/2+2.811316e+001*(z-5)^2/2+5.968454e+000*(z-7)^2/2-10*(z-2)+15*(z-6)+20*(z-4)^2/2-30*(z-6)^2/2-20*(z-0)^3/6+20*(z-5)^3/6
x є [8,10): EJx*theta7(z)=EJx*theta0+M0*z+Q0*z^2/2+R1*(z-2)^2/2+R2*(z-5)^2/2+R3*(z-7)^2/2+M1*(z-2)+M2*(z-6)+F1*(z-4)^2/2+F2*(z-6)^2/2+q1*(z-0)^3/6+q2*(z-8)^3/6-q1*(z-5)^3/6=2.007283e-013-8.953488e+000*z+2.343023e+001*z^2/2+4.139750e+001*(z-2)^2/2+2.811316e+001*(z-5)^2/2+5.968454e+000*(z-7)^2/2-10*(z-2)+15*(z-6)+20*(z-4)^2/2-30*(z-6)^2/2-20*(z-0)^3/6+15*(z-8)^3/6+20*(z-5)^3/6

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

disp('Изгибающие моменты на участках:')
Mz=[];
for k=1:nk,
  Mz1='M0+Q0*z';
  Mz2=sprintf('%d%+d*z',np(3:4));
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  Mz0=np(3)+np(4)*z0; % функции для графика
  if ~isempty(Op),
    nR=find(Op<zk(k+1)); % промежуточные опоры слева
    for k1=1:length(nR),
      Mz1=[Mz1 '+R' num2str(nR(k1)) ...
        '*(z-' num2str(Op(nR(k1))) ')'];
      Mz2=[Mz2 sprintf('%+d',np(nR(k1)+4)) '*(z-' ...
        num2str(Op(nR(k1))) ')'];
      Mz0=Mz0+np(nR(k1)+4)*(z0-Op(nR(k1)));
    end
  end
  if ~isempty(M)
    nM=find(M(:,1)<zk(k+1)); % изгибающие моменты слева
    for k1=1:length(nM),
      Mz1=[Mz1 '+M' num2str(nM(k1))];
      Mz2=[Mz2 sprintf('%+d',M(nM(k1),2))];
      Mz0=Mz0+M(nM(k1),2);
    end
  end
  if ~isempty(F)
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      Mz1=[Mz1 '+F' num2str(nF(k1)) '*(z-' ...
        num2str(F(nF(k1),1)) ')'];
      Mz2=[Mz2 sprintf('%+d',F(nF(k1),2)) '*(z-' ...
        num2str(F(nF(k1),1)) ')'];
      Mz0=Mz0+F(nF(k1),2)*(z0-F(nF(k1),1));
    end
  end
  if ~isempty(q)
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      Mz1=[Mz1 '+q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),1)) ')^2/2'];
      Mz2=[Mz2 sprintf('%+d',q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),1)) ')^2/2'];
      Mz0=Mz0+q(nq(k1),3)*(z0-q(nq(k1),1)).^2/2;
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      Mz1=[Mz1 '-q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),2)) ')^2/2'];
      Mz2=[Mz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),2)) ')^2/2'];
      Mz0=Mz0-q(nq(k1),3)*(z0-q(nq(k1),2)).^2/2;
    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)=M0+Q0*z+q1*(z-0)^2/2=-8.953488e+000+2.343023e+001*z-20*(z-0)^2/2
x є [2,4): M2(z)=M0+Q0*z+R1*(z-2)+M1+q1*(z-0)^2/2=-8.953488e+000+2.343023e+001*z+4.139750e+001*(z-2)-10-20*(z-0)^2/2
x є [4,5): M2(z)=M0+Q0*z+R1*(z-2)+M1+F1*(z-4)+q1*(z-0)^2/2=-8.953488e+000+2.343023e+001*z+4.139750e+001*(z-2)-10+20*(z-4)-20*(z-0)^2/2
x є [5,6): M4(z)=M0+Q0*z+R1*(z-2)+R2*(z-5)+M1+F1*(z-4)+q1*(z-0)^2/2-q1*(z-5)^2/2=-8.953488e+000+2.343023e+001*z+4.139750e+001*(z-2)+2.811316e+001*(z-5)-10+20*(z-4)-20*(z-0)^2/2+20*(z-5)^2/2
x є [6,7): M5(z)=M0+Q0*z+R1*(z-2)+R2*(z-5)+M1+M2+F1*(z-4)+F2*(z-6)+q1*(z-0)^2/2-q1*(z-5)^2/2=-8.953488e+000+2.343023e+001*z+4.139750e+001*(z-2)+2.811316e+001*(z-5)-10+15+20*(z-4)-30*(z-6)-20*(z-0)^2/2+20*(z-5)^2/2
x є [7,8): M6(z)=M0+Q0*z+R1*(z-2)+R2*(z-5)+R3*(z-7)+M1+M2+F1*(z-4)+F2*(z-6)+q1*(z-0)^2/2-q1*(z-5)^2/2=-8.953488e+000+2.343023e+001*z+4.139750e+001*(z-2)+2.811316e+001*(z-5)+5.968454e+000*(z-7)-10+15+20*(z-4)-30*(z-6)-20*(z-0)^2/2+20*(z-5)^2/2
x є [8,10): M7(z)=M0+Q0*z+R1*(z-2)+R2*(z-5)+R3*(z-7)+M1+M2+F1*(z-4)+F2*(z-6)+q1*(z-0)^2/2+q2*(z-8)^2/2-q1*(z-5)^2/2=-8.953488e+000+2.343023e+001*z+4.139750e+001*(z-2)+2.811316e+001*(z-5)+5.968454e+000*(z-7)-10+15+20*(z-4)-30*(z-6)-20*(z-0)^2/2+15*(z-8)^2/2+20*(z-5)^2/2

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

disp('Перерезывающие силы на участках:')
Qz=[];
for k=1:nk,
  Qz1='Q0';
  Qz2=sprintf('%d',np(4));
  z0=linspace(zk(k),zk(k+1)); % аргументы для графика
  Qz0=np(4)*ones(size(z0)); % функции для графика
  if ~isempty(Op),
    nR=find(Op<zk(k+1)); % промежуточные опоры слева
    for k1=1:length(nR),
      Qz1=[Qz1 '+R' num2str(nR(k1))];
      Qz2=[Qz2 sprintf('%+d',np(nR(k1)+4))];
      Qz0=Qz0+np(nR(k1)+4);
    end
  end
  if ~isempty(F),
    nF=find(F(:,1)<zk(k+1)); % сосредоточенные силы слева
    for k1=1:length(nF),
      Qz1=[Qz1 '+F' num2str(nF(k1))];
      Qz2=[Qz2 sprintf('%+d',F(nF(k1),2))];
      Qz0=Qz0+F(nF(k1),2);
    end
  end
  if ~isempty(q),
    nq=find(q(:,1)<zk(k+1)); % начало q слева
    for k1=1:length(nq),
      Qz1=[Qz1 '+q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),1)) ')'];
      Qz2=[Qz2 sprintf('%+d',q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),1)) ')'];
      Qz0=Qz0+q(nq(k1),3)*(z0-q(nq(k1),1));
    end
    nq=find(q(:,2)<zk(k+1)); % конец q слева
    for k1=1:length(nq),
      Qz1=[Qz1 '-q' num2str(nq(k1)) ...
        '*(z-' num2str(q(nq(k1),2)) ')'];
      Qz2=[Qz2 sprintf('%+d',-q(nq(k1),3)) ...
        '*(z-' num2str(q(nq(k1),2)) ')'];
      Qz0=Qz0-q(nq(k1),3)*(z0-q(nq(k1),2));
    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)=Q0+q1*(z-0)=2.343023e+001-20*(z-0)
x є [2,4): Q2(z)=Q0+R1+q1*(z-0)=2.343023e+001+4.139750e+001-20*(z-0)
x є [4,5): Q3(z)=Q0+R1+F1+q1*(z-0)=2.343023e+001+4.139750e+001+20-20*(z-0)
x є [5,6): Q4(z)=Q0+R1+R2+F1+q1*(z-0)-q1*(z-5)=2.343023e+001+4.139750e+001+2.811316e+001+20-20*(z-0)+20*(z-5)
x є [6,7): Q5(z)=Q0+R1+R2+F1+F2+q1*(z-0)-q1*(z-5)=2.343023e+001+4.139750e+001+2.811316e+001+20-30-20*(z-0)+20*(z-5)
x є [7,8): Q6(z)=Q0+R1+R2+R3+F1+F2+q1*(z-0)-q1*(z-5)=2.343023e+001+4.139750e+001+2.811316e+001+5.968454e+000+20-30-20*(z-0)+20*(z-5)
x є [8,10): Q7(z)=Q0+R1+R2+R3+F1+F2+q1*(z-0)+q2*(z-8)-q1*(z-5)=2.343023e+001+4.139750e+001+2.811316e+001+5.968454e+000+20-30-20*(z-0)+15*(z-8)+20*(z-5)

Подбор сечения по условиям прочности

Зададим ещё некоторые исходные данные: модуль упругости E и допускаемое напряжение [σ]. Они задаются соответственно в переменных E и Sd (в МПа). Найдём максимальный (по модулю) изгибающий момент Mmax и сечение, в котором он достигается (опасное сечение). Из соотношения

находим минимально допустимый момент сопротивления сечения:

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

Находим нужное двутавровое сечение. Печатаем его характеристики и пересчитанное по формуле (9) нормальное напряжение в опасном сечении. Если вам задали подобрать не двутавр, а, например, швеллер, посмотрите в ИДЗ Геометрические характеристики плоских сечений, как достать из соответствующей базы данных нужные характеристики.

E=2e5; % модуль упругости, МПа
Sd=160; % допускаемое напряжение, МПа
fprintf('Модуль упругости E=%d МПа\n',E);
fprintf('Допускаемое напряжение [sigma]=%d МПа\n',Sd);
[Mmax,iMmax]=max(abs(Mz)); % максимальный изгибающий момент
fprintf(['Максимальный по модулю изгибающий момент' ...
  '|Mmax|=%7.2f кНм\n'],Mmax)
fprintf(['достигается при z=%6.2f м\n'],z(iMmax))
Wmin=(Mmax)/Sd*1e3; % минимально допустимый Wx, см^3
fprintf(['Минимально допустимый момент сопротивления '...
  'Wmin=|Mmax|/[sigma]*10^5=%7.2f см^3\n'],Wmin);
load dwudata.txt % загружаем базу данных двутавров
dwu=find(dwudata(:,10)>Wmin); % двутавры, которые подходят
dwutavr=dwudata(dwu(1),:); % выбираем минимально возможный
Numdwu=dwutavr(1); % номер двутавра
if (Numdwu==round(Numdwu)),
  disp(['Выбираем друтавр No ' num2str(Numdwu)]);
else
  disp(['Выбираем друтавр No ' num2str(round(Numdwu)) 'a']);
end
disp('Его данные:')
Jx=dwutavr(9);
Wx=dwutavr(10);
Sx=dwutavr(11);
fprintf('Момент инерции Jx=%d см^4\n',Jx)
fprintf('Момент cопротивления Wx=%d см^3\n',Wx)
fprintf('Статический момент Sx=%d см^3\n',Sx)
Smax=Mmax/Wx*1000;
fprintf(['Максимальное нормальное напряжение в опасном сечении ' ...
  'SigmaMax=|Mmax|/Wx*10^3=%7.2f МПа\n'],Smax)
Модуль упругости E=200000 МПа
Допускаемое напряжение [sigma]=160 МПа
Максимальный по модулю изгибающий момент|Mmax|=  20.33 кНм
достигается при z=  6.00 м
Минимально допустимый момент сопротивления Wmin=|Mmax|/[sigma]*10^5= 127.07 см^3
Выбираем друтавр No 18
Его данные:
Момент инерции Jx=1290 см^4
Момент cопротивления Wx=143 см^3
Статический момент Sx=8.140000e+001 см^3
Максимальное нормальное напряжение в опасном сечении SigmaMax=|Mmax|/Wx*10^3= 142.18 МПа

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

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

[Qmax,iQmax]=max(abs(Qz)); % максимальная перерезывающая сила
fprintf(['Максимальная по модулю перерезывающая сила ' ...
  '|Qmax|=%7.2f кН\n'],Qmax)
fprintf(['достигается при z=%6.2f м\n'],z(iQmax))
fprintf('при z=%6.2f м: |Qmax|=%7.2f кН\n',...
  z(iMmax),abs(Qz(iMmax)))
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)
Максимальная по модулю перерезывающая сила |Qmax|=  24.83 кН
достигается при z=  2.00 м
при z=  6.00 м: |Qmax|=  17.06 кН
Касательные напряжения:
при z=  2.00: TauMax=  30.72 МПа
при z=  6.00: TauMax=  21.11 МПа

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

Нарисуем распределение нормальных и касательных напряжений по сечению. Нормальные напряжения распределены линейно, а касательные − по параболе. Мы строим эпюру распределения касательных напряжений приближённо: заменяем двутавр набором прямоугольников. Вычисляем по формуле (11) напряжения в крайних волокнах тонкой вертикальной стойки (переменная 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 Numdwu 0 0 -GeoDat(3)/2 -GeoDat(1)/2 0],GeoDat);
hold on
plot([0,-sign(Mz(iMmax))*Smax/10,sign(Mz(iMmax))*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, МПа')

Построение эпюр действительных прогибов и углов поворота

После того, как мы подобрали сечение, построим эпюры действительных перемещений и углов поворота (модуль упругости E и момент инерции сечения Jx у нас есть). Аналитические выражения мы уже не записываем, т.к. они отличаются от полученных ранее только множителем EJx. Находим максимальный прогиб и сечение, в котором он достигается.

figure
plot(z,EJwz/(E*Jx)*1e8,'k')
hold on
plot(Op,zeros(size(Op)),'o')
hold off
set(get(gcf,'CurrentAxes'),...
  'FontName','Times New Roman Cyr','FontSize',12)
title('\bfЭпюра прогибов')
xlabel('\itz\rm, м')
ylabel('\itw\rm(\itz\rm), мм')
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\rm(\itz\rm), радиан')
[wmax,iwmax]=max(abs(EJwz/(E*Jx)*1e8));
fprintf(['Максимальный по модулю прогиб'...
  ' |Wmax|=%8.4f мм\n'],wmax)
fprintf('достигается при z=%6.2f м\n',z(iwmax))
Максимальный по модулю прогиб |Wmax|=  3.9260 мм
достигается при z=  8.61 м


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

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

Литература

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