%% Расчет балки на упругом основании методом начальных параметров
% Расчетная схема: полубесконечная балка лежит на упругом основании свободно, под действием
% сосредоточенной силы Р, приложенной на расстоянии а от левого конца,
% учитывается распределенная нагрузка от собственного веса.
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 = -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

%вектор неизвестных параметров x=[v_0; fi_0; С1; С2]
Va = V(ksi_a);  %матрица функций Крылова в точке ksi_a

A = [Va(1,1)  Va(1,2)   -1          0;
     Va(2,1)  Va(2,2)    k         -k;
     Va(3,1)  Va(3,2)    0          2*k*k*EI;
     Va(4,1)  Va(4,2)   -2*k^3*EI  -2*k^3*EI];
 
 b = -q*Vq(ksi_a) + [q/(4*k^4*EI); 0; 0; P];
 
 x = A^-1*b;    %x=[v_0; fi_0; С1; С2]

%% Определение внутренних усилий и перемещений 
% Начальные параметры для первого участка (от левого конца балки до силы Р)
% 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;

%% Второй участок 
% Постоянные интегрирования для второго участка (от силы Р до L = 2*lambda вправо):

C1 = x(3);
C2 = x(4);

L = 2*lambda;   % параметр длины балки для построения графиков

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

for i=1:N2
    
    ii=i + N1;
    
    %глобальная координата
    z(ii) = z2(i) + a;
    
    %обозначим 
    ks = k*z2(i);
    
    %прогиб
    v(ii) = C1*exp(-ks)*cos(ks) + C2*exp(-ks)*sin(ks) + q/(4*k^4*EI);  
    
    %угол поворота сечения
    fi(ii)= -k*exp(-ks)*( (C1-C2)*cos(ks) + (C1+C2)*sin(ks) ); 
    
    %изгибающий момент
    M(ii) = EI*2*k^2*exp(-ks)*( C1*sin(ks) - C2*cos(ks) );
    
    %поперечная сила
    Q(ii) = EI*2*k^3*exp(-ks)*( (C1+C2)*cos(ks) - (C1-C2)*sin(ks) );
    
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;



