% Coupled electrmechanical dynamics  of MEMS structures
% using linear normal mode summation method and
% parallel plate approximation method
% Written for MEAM 550 in Fall 1998
% -- Suresh G.K.  Ananthasuresh

%
% This script needs the following m-files
% fem2.m  matcut.m veccut.m stfew.m
%

clear all
hold off
clc
clf
axis('normal');
echo off

global GK GM evec L B ic voltage v_time FM Flag freq ...
vdc vac omega gap e_force_elems ncnt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify the following
%
length = 80;  % Length of the beam  in microns
NELEM = 20;   % Number of elements
% Property value specification
H = 10*ones(NELEM,1); % Thicknesses of beam elements in microns
B = 0.5*ones(NELEM,1); % Widths of the beam elements
% NOTE: This is in plane vibration.
%
E = 169.0e3;    % Units: microN/micron^2
rho = 2231e-18; % Units: microKg/micron^3

vdc  = 20.0;
vac = 2.0;
omega = 4e6*2*pi;
gap = 2.0;

% Elements on which electrostatic force is applied.
e_force_elems = [6 7]

% Number of cycles to run the simulation
% This is based on first natural frequency
Ncycles = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NNODE = NELEM  + 1;
nx = [0:length/NELEM:length];
ny = zeros(size(nx));

for ielem=1:NELEM,
    ncnt(ielem,1) = ielem;
    ncnt(ielem,2) = ielem+1;
end


%
% Displacement boundary condition specification
% Fixed-fixed beam
dispID = [1 2 (NNODE-1)*2+1 (NNODE-1)*2+2];
dispVal = [0 0 0 0];

% Lengths of elements
for ie=1:NELEM,
  eye = ncnt(ie,1);
  jay = ncnt(ie,2);
  L(ie) = sqrt ( (nx(jay) - nx(eye))^2 + (ny(jay) - ny(eye))^2 );
end

%
% Call the fem m-file
% K = stiffness matrix
% M = inertia (mass) matrix
%
[K, M] = fem2(B,H,L,E,rho,nx,ny,ncnt,NELEM,NNODE,dispVal,dispID);

% Compute mode shapes and natural frquencies
% evec = eigenvectorr matrix
% eval = diagonal matrix containing the eigenvalues
D = inv(K)*M;
[evec eval] = eig(D);

% Normalize  mode  shapes so that the maximum transverse
% displacement is unity.
sz  =  size(evec,2);
evec(:,1) = evec(:,1) / max(abs(evec(1:2:sz,1)));
evec(:,2) = evec(:,2) / max(abs(evec(1:2:sz,2)));
evec(:,3) = evec(:,3) / max(abs(evec(1:2:sz,3)));
evec(:,4) = evec(:,4) / max(abs(evec(1:2:sz,4)));

% Generalized inertia matrix
GM = evec'*M*evec;

% Generalized stiffness matrix
GK = evec'*K*evec;

% Natural frequencies
freq = sqrt(diag(GK)./diag(GM));

% Compare frequencies with the  formula estimate of frequency

rho_l = rho*B(1)*H(1); % mass per unit length
omega1 = sqrt( (E*B(1)*H(1)^3/12) / (rho_l*length^4) ) * 22.4;
omega2 = sqrt( (E*B(1)*H(1)^3/12) / (rho_l*length^4) ) * 61.7;
omega3 = sqrt( (E*B(1)*H(1)^3/12) / (rho_l*length^4) ) * 121;

echo  on
% Computed & formula estimate
[freq(1) omega1; freq(2) omega2; freq(3) omega3;]
% percentage error in computation
[(omega1-freq(1))/omega1*100; (omega2-freq(2))/omega2*100; ...
 (omega3-freq(3))/omega3*100]
echo off

% plot mode shapes
plot([0 nx(2:NNODE-1) length], [0 evec(1:2:sz,1)' 0],'-r', ...
     [0 nx(2:NNODE-1) length], [0 evec(1:2:sz,2)' 0],'--g', ...
     [0 nx(2:NNODE-1) length], [0 evec(1:2:sz,3)' 0],':k',  ...
     [0 nx(2:NNODE-1) length], [0 evec(1:2:sz,4)' 0],'-.m');
  hold on
plot([0 nx(2:NNODE-1) length], [0; evec(1:2:sz,1); 0],'or', ...
     [0 nx(2:NNODE-1) length], [0; evec(1:2:sz,2); 0],'og', ...
     [0 nx(2:NNODE-1) length], [0; evec(1:2:sz,3); 0],'ok',  ...
     [0 nx(2:NNODE-1) length], [0; evec(1:2:sz,4); 0],'om');
legend('Mode1', 'Mode2', 'Mode3', 'Mode4'); grid
xlabel('Length along the  beam,  x');
ylabel('Transverse displacement of the modes');

echo on
% Frequencies in Hertz
[freq(1) freq(2) freq(3) freq(4)]/2/pi

% Press any key  to continue
pause
echo  off

fun = 'stfew';
T0 = 0.0;
Tfinal = Ncycles * 2*pi/freq(1);


% Take only first few modes for simulation
FM = 4; %  Number  of modes considered

Y0 = zeros(1,2*FM);

%
% Time interval for obtaining results from  Runge-Kutta integration
% periodically = deltaT
% Nsteps = number of intervals for Runge-Kutta integration
%
Nsteps = 100;
deltaT = Tfinal/Nsteps;
Tbegin = 0;

ic = 0; % Index used to  keep track of voltage as a function  of time
irk = 0;
t_start_cpu = cputime;

for t = deltaT:deltaT:Tfinal,
    irk = irk + 1;
    Tend = t; 
    [Ts, Ys] = ode23(fun,[Tbegin Tend], Y0);
    last = size(Ts,1);
    for i = 1:2*FM,
        temp = Ys(last,i);
        Y(irk,i) = temp;
    end
    T(irk) = t;
    Tbegin = Tend;
    Y0 = Ys(last,:);
    t
end

% End CPU clock time
t_run = cputime - t_start_cpu;


% Plot deformed beam at different time instants
% Take FM modes
idvec = 1:2:2*FM;
x_soln = Y(:,idvec);
x_soln = evec(:,1:FM)*x_soln';

clf

pmax = size(T,2);
for ip=1:1:pmax,
   plot([0 nx(2:NELEM) length], [0; x_soln(1:2:size(K,2),ip); 0],'-y');
   hold on
end

echo on
% Presss any key to continue
pause

% Plot first four modal coordinate variation with time
clf
subplot(2,2,1);
plot(T, Y(:,1)); ylabel('Mode1'); xlabel('Time (sec)');
subplot(2,2,2);
plot(T, Y(:,3)); ylabel('Mode2'); xlabel('Time (sec)');
subplot(2,2,3);
plot(T, Y(:,5)); ylabel('Mode3'); xlabel('Time (sec)');
subplot(2,2,4);
plot(T, Y(:,7)); ylabel('Mode4'); xlabel('Time (sec)'); 

% Press any key to continue
pause


% Plot disp variation with time of the mid node 
clf
mid = NELEM - 1;
plot(T, x_soln(mid,:),'-y'); hold on; plot(T, x_soln(mid,:),'.c')
grid;
xlabel('Time (sec)');
ylabel('Disp_max');

% Press any key to continue
pause

% Mid node Velocity

x_dot = evec(:,1:FM)*Y(:,2:2:2*FM)';
clf
plot(T, x_dot(mid,:),'-y'); hold on; plot(T, x_dot(mid,:),'.c')
grid;
xlabel('Time (sec)');
ylabel('Vel_max');

% Press any key to continue
pause

%  Plot the  voltage applied as a  function of time
clf
plot(v_time, voltage);
xlabel('Time (sec)');
ylabel('Voltage (V)');








