%lab5.m - solves the 4th order beam bending equation by the finite
%difference method for steady-state and beam not touching the substrate
%- fixed-fixed beam boundary conditions
% For lab 5 in EE/ME 550

clear all
%beam geometry
Lb = 100; %um - beam length
tb = 2; %um - beam thickness
wb = 10; %um - beam width
gap = 1; %um - gap under beam

%beam properties
sigma0 = 20; %MPa - residual stress in the film
V = 100; %V

%material properties
E = 165e3; %MPa - Young's modulus
nu = 0.2; %Poisson's ratio (unitless)

%mesh properties
N = 1000; %number of elements in the beam - must be at least four
h = Lb/(2*N); %step size in x
tol = 1e-6; %ratio of change in deflection to deflection magnitude

%a few derived constants
Epl = 1/(1/E - nu^2/E*(((wb/Lb))/(0.18 + (wb/Lb)))^(1.77*(Lb/tb)^-0.061)); %MPa - plate modulus
I = wb*tb^3/12; %um^4 - moment of inertia
Area = wb*tb; %um^2 - cross-sectional area
rig = Epl*I; %uN*um^2 - beam flexural rigidity
P = sigma0*(1-nu)*Area; %uN - force at the supports due to residual stress


%analyze the beam 
%mesh the beam and initialize some variables
wu = zeros(N+1,1);
A = sparse(N,N);
B = sparse(N,N);
b = zeros(N,1);
x = [0:h:Lb/2]';
fe = zeros(N,1);

%calculate some constants which will be used in the equations
c1 = h^4/rig; %um^2/uN
g1 = -gap; %um
c2 = P*h^2/rig; %unitless

%set up the matrices
A(1,1:3) = [7 -4 1];
B(1,1:2) = [-2 1];
A(2,1:4) = [-4 6 -4 1];
B(2,1:3) = [1 -2 1];
A(N-1,N-3:N) = [1 -4 7 -4];
B(N-1,N-2:N) = [1 -2 1];
A(N,N-2:N) = [2 -8 6];
B(N,N-1:N) = [2 -2];
for counter = 3:N-2
    A(counter,counter-2:counter+2) = [1 -4 6 -4 1];
    B(counter,counter-1:counter+1) = [1 -2 1];
end
B = c2*B;

%Now, we will iterate on the solution until the beam stops deflecting
finished = 0;
wlast = 2*wu+4;

while finished==0
    fe = electf(wu(2:N+1),gap,wb,V);
    bss = b + c1*fe;
    wu(2:N+1) = (A - B)\bss;
    if (max(abs(wu(2:N) - wlast(2:N))./abs(wu(2:N)))<tol)|min(wu+gap)<0
        finished = 1;
    end
    wlast = wu;
end
%Check whether pull-in occurred, and display the maximum deflection if it
%didn't.
if min(wu+gap)<0
    display('Pull-In')
else
    display('Maximum Deflection is')
    display(num2str(max(abs(wu))))
    
    %Plot the beam's deflected shape.
    figure(1)
    clf
    plot(x,wu)
    xlabel('Beam x-coordinate (\mum)')
    ylabel('Deflection (\mum)')
    title('Deflection of Half of the Beam')
end