% Energy method to approximate buckling load
% in lateral torsional buckling
% ----------------------------------------------------
%  Canteliver beam - lateral buckling under tip load P
%    at the center of gravity G 
%    Author: Baroudi D. 2021
% ---------------------------------------------------
syms EI  L P  v v0 x
syms delta_P delta_F
syms P1 P2 
 
syms EI_y GI_t   
syms theta0  w0   theta w

% ---------------------------------------
% Displacement approximation  (you can use better approximations)
% ----------------------------------
%% w(x, w0, L)    =  w0 *  x.^2 / L^2   % less good approx.
w(x, w0, L)    =  w0 *  ( 1 - cos(pi * x/ (2*L) ) )  % better one
d1_w(x, w0, L) = simplify( diff(w, x) )
d2_w(x, w0, L) = simplify( diff(d1_w, x) )
 
% twist angle approximation (you can use better aprox.)
% -----------------------------
phi(x, L)   =  theta0 * x /L 
d1_phi(x, theta0, L) = simplify( diff(phi, x) )
d2_phi(x, theta0, L) = simplify( diff(d1_phi, x) );
%-------------------------------------

% Strain energy
delta_U_B(L, EI_y, w0) = 0.5* EI_y * int(d2_w * d2_w, [0 L])
delta_U_T(L, GI_t, theta0) = 0.5 * GI_t * int(d1_phi * d1_phi, [0 L])
% --- work of initial stresses (bending M0z)
% ------------------------------------------

% Work increment of initial stresses
%--------------------------------
% Initial bending moment: M0z
% ----------------------
 M0z(x, P, L) = -P * (L -x)
d_1_Mz0_Phi(x, P, L, theta0) = simplify( diff( M0z(x, P, L)*phi(x, L) , x) )
delta_W(P, L, w0, theta0) =  int( d_1_Mz0_Phi(x, P, L, theta0) * d1_w(x, w0, L), [0 L])

% ------------------------------------
% total increment of potential energy
% -----------------------------------

delta_Pi = delta_U_B(L, EI_y, w0) + delta_U_T(L, GI_t, theta0) ...
            + delta_W(P, L, w0, theta0) 

% Equations of neutral equilibrium
% ------------------------------------
delta_Pi_w   = simplify ( diff(delta_Pi, w0) )
delta_Pi_phi = simplify ( diff(delta_Pi, theta0) )
% ----------------------------------------------