% Energy method to approximate buckling load
% for a stiffened thn plate with in-plane compressive Nox along one side
% the Poisson expansion is not restrained by the supports.
% ----------------------------------------------------
%    Author: Baroudi D. 2021
% ---------------------------------------------------
syms x y
syms delta_P   delta_W
   
syms w w0
syms L b  D  EI nu
syms N0x

% ---------------------------------------
% Displacement approximation  (you can use better approximations)
% ----------------------------------
%% w(x, y, w0, L)    =  w0 / ((L^2) * (b^2)) * x * ( L - x) * y * (b - y)  % less good than the trigonometric
 w(x, y, w0, L)    =  w0 *  sin(pi* x /L) * sin(pi* y /b) % best one

d1x_w(x, y, w0, L, b)  = simplify( diff(w, x) )
d2xy_w(x, y, w0, L, b) = simplify( diff(d1x_w, y) )
d2x_w(x, y, w0, L, b)  = simplify( diff(d1x_w, x) )

d1y_w(x, y, w0, L, b)  = simplify( diff(w, y) ) 
d2yx_w(x, y, w0, L, b) = simplify( diff(d1y_w, x) )
d2y_w(x, y, w0, L, b)  = simplify( diff(d1y_w, y) )  
 

% Strain energy (plate alone)
delta_U(L, b, D, w0) = 0.5* D * int( int(d2x_w * d2x_w, x, [0 L]) , y, [0 b]) + ...
                       0.5* D * int( int(d2y_w * d2y_w, x, [0 L]) , y, [0 b]) + ...
                       nu * D * int( int(d2x_w * d2y_w, x, [0 L]) , y, [0 b]) + ...
                  (1 - nu)* D * int( int(d2xy_w * d2xy_w, x, [0 L]) , y, [0 b]) 
              
% Strain energy (stiffner beam alone)
d2x_w_beam = d2x_w(x, b/2, w0, L, b)
%% delta_U_beam(L, b, D, w0) = 0.5* EI * int( int(d2x_w_beam * d2x_w_beam, x, [0 L]) , y, [0 b])
delta_U_beam(L, b, D, w0) = 0.5* EI * int( d2x_w_beam * d2x_w_beam, x, [0 L]  )
          
                
% Work increment of initial stresses (applied N0x along the boundaries x=0
% and x=L
%--------------------------------
delta_W(L, b, w0, N0x) = 0.5 * int( int(N0x * d1x_w * d1x_w, x, [0 L]) , y, [0 b]   )
%% Note that here Nox is negative
% ------------------------------------
% Total increment of potential energy
% -----------------------------------
delta_Pi = delta_U(L, b, D, w0) + delta_W(L, b, w0, N0x) + delta_U_beam(L, b, D, w0)
delta_Pi = simplify( delta_Pi)
% Equations of neutral equilibrium
% ------------------------------------
delta_Pi_w   = simplify ( diff(delta_Pi, w0) )
% ----------------------------------------------
texti = 'Remember Nox is now negative : = - Nox_ref'