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

clear all;
%Построение упругой линии балки при прямом изгибе методом начальных параметров.
%Шарнирно-опертая балка под действием центральной сосредоточенной силы.
 
%Исходные данные:
l = 2.4;                % м, длина первого и второго участка балки
a = 0.8;                % м, длина третьего участка балки
I = 17.5*10^-5;        % м^4, момент инерции сечения
E = 2*10^11;            % Па, модуль Юнга для стали
EI = E*I;               % жесткость при изгибе
q = 2300;               % H/м, максимальная допускаемая нагрузка
 
%Начальные параметры для первого участка:
v_0 = 0;                % м, прогиб               
fi_0 = 0;               % рад,  угол поворота 
M_0 = q*(l^2+a^2);      % Нм, момент
Q_0 = (q*(l^2+a^2))/l;  % Н, сила
q_0 = 2300;             % Н/м, интенсивность распределенной нагрузки
 
%Первый участок
N1 = 20;                    % число точек  
z1 = linspace(0, l, N1);     % координата 
 
for i=1:N1
    
    % координата
    z(i) = z1(i);
    
    %поперечная сила
    Q(i) = Q_0 - q_0*z1(i);  
    
    %изгибающий момент
    M(i) = M_0 - Q_0*z1(i) + q*(z1(i))^2/2; 
    
    %угол поворота сечения
    fi(i)= fi_0 + M_0*z1(i)/EI + Q_0*(z1(i))^2/(2*EI) + q_0*(z1(i))^3/(6*EI);  
    
    %прогиб
    v(i) = v_0 + fi_0*z1(i) + M_0*z1(i)*z1(i)/(2*EI) + Q_0*(z1(i))^3/(6*EI) + q_0*(z1(i))^4/(24*EI);  
    
end;
 
%Начальные параметры для второго участка:
v_0 = v(N1);            % м, прогиб               
fi_0 = fi(N1);          % рад,  угол поворота 
M_0 = M(N1);            % Нм, момент
Q_0 = q*a^2/l;        % Н, сила
q_0 = 2300;                % Н/м, интенсивность распределенной нагрузки
 
%Второй участок
N2 = 20;                     % число точек  
z2 = linspace(l, 2*l, N2);   % координата 
 
for i=1:N2
    
    ii=i + N1;
    
    %глобальная координата
    z(ii) = z2(i);
    
    %поперечная сила
    Q(ii) = Q_0 + q_0*(z2(i));  
    
    %изгибающий момент
    M(ii) = q*l^2/2 -q*l^2*(z2(i))/l  + q_0*(z2(i))^2/2; 
    
    %угол поворота сечения
    fi(ii)= fi_0 + M_0*z2(i)/EI + Q_0*(z(i))^2/(2*EI) + q_0*(z2(i))^3/(6*EI);  
    
    %прогиб
    v(ii) = v_0 + fi_0*z2(i) + M_0*z2(i)*z2(i)/(2*EI) + Q_0*(z2(i))^3/(6*EI) + q_0*(z2(i))^4/(24*EI);  
    
end;
 
%Начальные параметры для третьего участка:
v_0 = v(N2);            % м, прогиб               
fi_0 = fi(N2);          % рад,  угол поворота 
M_0 = q*a^2;            % Нм, момент
Q_0 = (q*(l^2+a^2))/l ;        % Н, сила
q_0 = 0;                % Н/м, интенсивность распределенной нагрузки
 
%Третий участок
N3 = 20;                     % число точек  
z3 = linspace(2*l, 2*l+a, N3);   % координата 
 
 for i=1 : N3
    iii = i + N2;
    %глобальная координата
    z(iii) = z3(i);
    
    %поперечная сила
    Q(iii) = q_0*(z3(i))^2/l;  
    
    %изгибающий момент
    M(iii) = M_0; 
    
    %угол поворота сечения
    fi(iii)= fi_0 + M_0*z3(i)/EI + Q_0*(z3(i))^2/(2*EI) + q_0*(z3(i))^3/(6*EI);  
    
    %прогиб
    v(iii) = v_0 + fi_0*z3(i) + M_0*z3(i)*z(i)/(2*EI) + Q_0*(z3(i))^3/(6*EI) + q_0*(z3(i))^4/(24*EI);  
end;
 
subplot(4,1,1), area(z, Q, 'FaceColor', [0.75 1 0.75]), xlabel('z, m'); ylabel('Q, H'), title([' Поперечная сила Q ']), grid on,
subplot(4,1,2), area(z, M, 'FaceColor', [0.75 1 0.75]); xlabel('z, m'); ylabel('M, Hm'), title([' Изгибающий момент М ']),grid on,
subplot(4,1,3), area(z, fi, 'FaceColor', [0.75 1 0.75]); xlabel('z, m'); ylabel('\phi, rad'), title([' Угол поворота \phi ']),grid on,
subplot(4,1,4), area(z, v, 'FaceColor', [0.75 1 0.75]); xlabel('z, m'); ylabel('v, m'), title([' Прогиб v ']),grid on