%% Расчет балки на упругом основании методом начальных параметров
% Расчетная схема: балка лежит на упругом основании свободно, под действием
% сосредоточенной силы Р, приложенной на расстоянии а от левого конца,
% учитывается нагрузка от собственного веса.
clear all;
%% Параметры балки
E = 2e11;                   %Па, модуль упругости
b = 0.135;                  %м, ширина полки, двутавр №30
I = 7080e-8;                %м4, момент инерции отн. оси х
EI = E*I;                   %изгибная жесткость
W = 472e-6;                 %м3, момент сопротивления отн. оси х
F = 46.5e-4;                %м2, площадь сечения
m = 36.5;                   %кг, вес 1 м балки (погонная масса)
L = 6;                      %м, длина балки
%% Параметры упругого основания
c = 4000000;                %H/м3, коэффициент постели пластичной глины
k_osn = c * b;              %H/м2, коэффициент упругости основания

%% Параметры нагрузки
q = 0;%-m*9.81;     %Н/м, погонный вес балки
P = -10e3;       %H, сосредоточенная сила, приложена на расстоянии а от левого конца 
a = L/3;

%% Волновое число и длина полуволны краевого эффекта
k = (k_osn/(4*EI))^0.25;    %волновое число
lambda = pi/k;              %длина полуволны краевого эффекта

%% Определение вектора начальных параметров
% Для свободно лежащей балки с незагруженными концами 
% момент и поперечная сила на концах равны нулю, 
% прогибы и углы поворота на концах не известны.    
global EI;
global k;           %глобальные переменные

ksi_a = k * a;      %значение безразмерной координаты ksi в точке z=a
ksi_L = k * (L-a);  %значение безразмерной координаты ksi в точке z=L

Matr = V(ksi_L) * V(ksi_a); % *вывод формул см. занятие*
P_vect = [0; 0; 0; P];
Vect = V(ksi_L) * ( q*Vq(ksi_a) + P_vect ) + q*Vq(ksi_L);

%вектор неизвестных начальных параметров x=[v_0; fi_0; v_L; fi_L]
A = [Matr(1,1)  Matr(1,2)   -1  0;
     Matr(2,1)  Matr(2,2)    0 -1;
     Matr(3,1)  Matr(3,2)    0  0;
     Matr(4,1)  Matr(4,2)    0  0];
 b = -Vect;
 
 x = A^-1*b    %x=[v_0; fi_0; v_L; fi_L]

%% Определение внутренних усилий и перемещений 
% Начальные параметры для первого участка (от левого конца балки до силы Р)
% 0<= ksi<= ksi_a
V_0 = [x(1); x(2); 0; 0]; %[прогиб, угол, момент, сила] при z=0

%% Первый участок
% 
N1 = 20;                     % число точек  
z1 = linspace(0, a, N1);     % массив координат 

for i=1:N1
    
    % координата
    z(i) = z1(i);
    ksi = k*z1(i);
    Vtmp = V(ksi);      %значение матрицы функций Крылова  
    Vqtmp = Vq(ksi);    %значение вектора частного решения 
    
    %прогиб
    v(i) = Vtmp(1,:)* V_0 + q*Vqtmp(1);  
    
    %угол поворота сечения
    fi(i)= Vtmp(2,:)* V_0 + q*Vqtmp(2);  
    
    %изгибающий момент
    M(i) = Vtmp(3,:)* V_0 + q*Vqtmp(3); 
    
    %поперечная сила
    Q(i) = Vtmp(4,:)* V_0 + q*Vqtmp(4);  
    
end;

%% Второй участок 
% Начальные параметры для второго участка (от силы Р до правого конца балки):
% 0 <= ksi<= ksi_L - ksi_a
V_0 = [v(N1); fi(N1); M(N1); (Q(N1) + P)]; %[прогиб, угол, момент, сила] при z=a

N2 = 20;                     % число точек  
z2 = linspace(a, L, N2);     % массив координат 

for i=1:N2
    
    ii=i + N1;
    
    %глобальная координата
    z(ii) = z2(i);
    ksi = k*( z2(i) - a );  %в начале участка ksi должно быть равно 0
    
    Vtmp = V(ksi);      %текущее значение матрицы функций Крылова  
    Vqtmp = Vq(ksi);    %текущее значение вектора частного решения  
    
    %прогиб
    v(ii) = Vtmp(1,:)* V_0 + q*Vqtmp(1);  
    
    %угол поворота сечения
    fi(ii)= Vtmp(2,:)* V_0 + q*Vqtmp(2);  
    
    %изгибающий момент
    M(ii) = Vtmp(3,:)* V_0 + q*Vqtmp(3); 
    
    %поперечная сила
    Q(ii) = Vtmp(4,:)* V_0 + q*Vqtmp(4);  
    
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

%% Расчет на прочность
% *Выполнить самостоятельно*
maxM = max(M)
max_sigma = maxM/W



