% This program solves the beam differential equation:
%
%    EIg'''' - Tg'' = -eV^2w/2g^2 (1 + 0.65g/w)
% 
% for a fixed/fixed beam over a ground plane using the finite difference
% method. Note that the fringing term has been included in the above equation.
% There are two Modes to this program.  The Normal Mode is the familiar
% single voltage/single deflection ouput mode.  The new Pull-In Mode
% will automatically iteratively calculate the exact Pull-In Voltage for
% the input'ed case to within 2 decimal points precison without the user
% having to manually iterate in the Normal Mode as was the case in the past.
% Peter Osterberg, MIT
% 4/1/94

% This is now updated for obtaining more precise deflection profiles
% of a fixed-fixed beam under electrostatic load, especially near
% pullin.  The code has been tightened to eliminate redundancies
% and allows for more user input into convergence criteria.
% Raj Gupta, MIT
% 1st October 1994

% A stress stiffening option is now added to this model.  A uniform
% tension is added to any pre-existing tension due to residual stress
% by calculating the strain incurred by deflecting the beam.  I've
% also changed the Newton's solver so that it's correction is not directly
% added to the previous solution, but instead some fraction of it is.
% In particular, one-half of this correction value is added to avoid the
% oscillatory behavior found in some of the more stiffer nonlinear
% cases.  This hack gives the same voltage-deflection profile after
% convergence, except at a cost of time.
%
% It is interesting to note that with the stress stiffening option
% on, the scripts obtain two stable equilibrium states for voltages
% less than pullin.  I am not sure if this is an artifact of the
% numerical solver, or a reality.  Also interesting is that if we
% start off with the larger of the two coverged deflected states, it 
% is easier to approach solutions near the pullin voltage, making it 
% easier to eventually find the exact pullin voltage.
% Raj Gupta, MIT
% 1st February 1995

% Pull-in calculations by Mode 1 show that pull-in voltages do 
% go to zero at almost exactly the theoretical buckling limit 
% for ideal fixed-fixed beams.  However, the Hilton Head analytic 
% model is very poor in this regime close to buckling, and
% predicts a much higher pull-in voltages than obtained directly
% by the numerical simulation of Mode 1.
% Raj Gupta, MIT
% 25th February 1995

% Because a moment induced on a stepped support causes rotation
% of the support, I have calculated the effects due to stress-gradients,
% net residual stress, stress-stiffening, and electrostatic loading.  
% Note, this moment does nothing, if the support is ideally fixed.
% Raj Gupta, MIT
% 4th October 1995 & 6th January 1996

% I have realized that the stress stiffening obtion does not work
% using Newton's method for calculating the differential equation
% at voltages near pull-in.  This problem can be solved using
% other finite-difference based techniques, but will not be 
% implemented in this script.
% Raj Gupta, MIT
% June 1996

% Minor changes made to improve convergence.
% Raj Gupta, MIT
% June 1996


case1=input('New Case ? (1=yes or 0=no): ');

if case1==1

  
  ask=input('Enter discretization number (50 to 100 range is suggested): ');
  if (~isempty(ask)) d = ask; end;
  N=d+4;
  ask=1e-6*input('Enter length in um: ');
  if (~isempty(ask)) L = ask; end;
  ask=1e-6*input('Enter width in um: ');
  if (~isempty(ask)) w = ask; end;
  ask=1e-6*input('Enter thickness in um: ');
  if (~isempty(ask)) t = ask; end;
  ask=1e-6*input('Enter gap in um: ');
  if (~isempty(ask)) g = ask; end;
  ask=1e9*input('Enter elastic modulus in GPa: ');
  if (~isempty(ask)) E = ask; end;
  ask=input('Enter poisson ratio: ');
  if (~isempty(ask)) nu = ask; end;
  ask=1e6*input('Enter biaxial residual stress in MPa: ');
  if (~isempty(ask)) s = ask; end;
  ask=1e-6*input('Enter step height in um: ');
  if (~isempty(ask)) step = ask; end;
  if step ~= 0
    ask = 1e12*input('Enter transverse stress gradient (MPa/um): ');
    if (~isempty(ask)) dsdy = ask; end;
  else
    dsdy = 0;
  end

  half = d + 2;
  r = g*ones(N,1);
  e = 8.8541878e-12;
  rho = 2331;
  h = L/(2*(d-1));
  % Two widths allow for trapezoidal geometry, 
  % w1 = top width, w = bottom width
  w1 = w;
  I = t^3/36 * (w^2+4*w*w1+w1^2)/(w+w1);
  T = s*(1-nu)*w*t;
  a = T/(E*I);
  Bckl_stress_in_MPa = -8*pi^2*E*I/(L^2*(1-nu)*t*(w+w1))/1e6
  % This could use modification with the new 'I'
  k1=32*E*t^3/L^4;
  % This also could use modification
  k2=8*s*(1-nu)*t/L^2;
  k=k1+k2;
  V_pull_in_estimate = (8*k*g^3/(27*e))^0.5
  g1 = 2.79;
  g2 = 0.97;
  g3 = 0.42;
  u=((18*s*(1-nu)*(w+w1)^2/(E*t^2*(w^2+w1^2+4*w*w1)))^0.5)*(L/2);
  V_pull_in_HH=real(((g1*g^3*s*(1-nu)*t)/(L^2*e*(1+g3*g/w)* ...
	(1+(2*(1-cosh(g2*u)))/(g2*u*sinh(g2*u)))))^0.5)
  alpha = 0.8;
  SLB = pi^2/3 + s*(1-nu)*L^2/(E*t^2);
  Cn = (g/t)^2;
  V_pull_in_EM=3.474541868/L^2*(E*t^3*g^3/e*(1+3*s*(1-nu)*L^2/(pi^2*E*t^2)))^(1/2)
%  V_pull_in_HH_ss=V_pull_in_HH*(1 + 0.2664*Cn/SLB - 0.02867*Cn/SLB^2 ...
%	- 0.004352*Cn^2/SLB - 0.00008376*(Cn/SLB)^2)
  V_pull_in_HH_ss=V_pull_in_HH*(1 + 0.206524*Cn/SLB - 0.008042074*Cn/SLB^2 ...
	- 0.002744049*Cn^2/SLB + 0.000097048*(Cn/SLB)^2)


  % Note the following estimates of resonant frequency in Hz are based on stress = 0

  ResFreq_1st_E = (sqrt((E*w*t^3/12)/(rho*w*t*L^4))*22.4)/(2*pi)
  ResFreq_2nd_E = ResFreq_1st_E * 61.7/22.4;
  ResFreq_3rd_E = ResFreq_1st_E * 121/22.4;

  % This estimate of the fundamental frequency in Hz includes effects of stress.
  % (from Tilmans' Ph.D. thesis, page 165) 

  ResFreq_1st_E_s = (1.0279*sqrt(E/(rho*L^2)*((t/L)^2+0.295*s*(1-nu)/E)))

  Mo = dsdy*I + T*step;
  MoEI = Mo/(E*I);

  y=g*ones(N,1);
  A=zeros(N,N);
  C=zeros(N,1);
  y2p=zeros(N,1);
  y4p=zeros(N,1);

  % Boundary conditions:
  %	y(3) = 0	=>	no deflection at support
  %	y'(3) = ?	=>	determined later in script as slope at support due to 
  %					  moments this is normally zero, if support is truly "fixed".
  %	y'(N-2) = 0	=>	zero slope at beam center, by symmetry
  %	y'''(N-2) = 0	=>	no shear force at beam center - this is the derivative of the
  %						  moment, which is proportional to the radius of curvature.
  
  if (s <= Bckl_stress_in_MPa*1e6)
    for i=1:N
      y(i,1)=(1 + 1e-2*sign(Mo+1e-100)*(s/Bckl_stress_in_MPa*1e-6)*(1-cos(2*pi*(i-3)*h/L)))*g;
    end
  end
  
  A(1,3) = 1;
  A(2,2) = -1;
  A(2,4) = 1;
  A(N-1,N-3) = -1;
  A(N-1,N-1) = 1;
  A(N,N-4) = -1;
  A(N,N-3) = 2;
  A(N,N-1) = -2;
  A(N,N) = 1;
  for i = 3:N-2
    A(i,i-2) =  1;
    A(i,i-1) = -4-a*h^2;
    A(i,i+1) = -4-a*h^2;
    A(i,i+2) =  1;
  end

end

% OLD VALUE OF y(i)'s WILL BE USED UNLESS THEY ARE GARBAGE
% THIS IS ESPECIALLY IMPORTANT IF WE'RE IN MODE = 1 AND
% CLOSE TO THE PULLIN VOLTAGE.

if (y(half)<0 | y(half)>g)
  y = g*ones(N,1);
end

mode=input('Normal Mode or Pull-In Mode ? (1=Normal or 0=Pull-In): ');

if mode==1

% BEGIN MODE 1 = NORMAL

% NOTE: Stress-stiffening option will not work for V close to VPI.

ask=input('Include effects of stress stiffening ? (1=yes or 0=no): ');
if (~isempty(ask)) ss = ask; end;

V=input('Enter voltage: ');

% NEWTON ITERATIONS START HERE FOR MODE 1

b = -e*V^2*w/(2*E*I);
c = 0.65*b/w;
phi = 0;
r(half) = 1;
count = 0;
while abs(r(half)/y(half)) > 1e-8
  if (ss == 1)
    newL = 0;
    for i=3:N-3
	newL = newL + (1 + (((y(i+1)-y(i-1)))/(2*h))^2)^0.5 * h;
    end
    if (step ~= 0)
      strain = (2*newL-L+2*phi*step)/L;
    else
      strain = (2*newL-L)/L;
    end
    a = (T + (1-2*nu)*(strain*E)*w*t)/(E*I);
    Mo = (dsdy + a*step*E)*I;
    MoEI = Mo/(E*I);
    for i=3:N-2
      A(i,i-2) =  1;
      A(i,i-1) = -4-a*h^2;
      A(i,i+1) = -4-a*h^2;
      A(i,i+2) =  1;
    end
  end
  for i=3:N-2
    y2p(i) = (y(i-1)-2*y(i)+y(i+1))/h^2;
    y4p(i) = (y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
    A(i,i)   =  6+2*a*h^2+2*h^4/y(i)*(y4p(i)-a*y2p(i)-c/(2*y(i)));
    C(i) = h^4*(b/y(i)^2-y4p(i)+a*y2p(i)+c/y(i));
  end
  if (step ~= 0)
    phi = (-MoEI + (y(4)-2*y(3)+y(2))/(h^2)) * step;
    y1p3 = (y(4)-y(2))/(2*h);
    C(2) = 2*h*(phi-y1p3);
  end
  r=A\C;
  if count > -1
    y=y+(r/20);
  else
    y=y+r;
  end
  y(half);
  count = count+1;
end

%END MODE 1 = NORMAL

else

phi =0;
Z = input('Enter newton iteration number (odd number > 5 suggested): ');
%Z=5;
Q = input('Set the limits of binary search for pullin (1=yes & 0=no): ');
if Q == 1
  VMIN = input('Enter lower limit for pullin: ');
  VMAX = input('Enter upper limit for pullin: ');
else
  if (s < 0)
    VMIN = 0;
  else
    VMIN = V_pull_in_estimate;
  end
  VMAX = V_pull_in_EM;
end

V = VMIN;

%BEGIN BINARY SEARCH FOR PULLIN VOLTAGE

while abs(VMAX-VMIN) > 1e-8

% y(i) is initialized to be g for all i''s upon inception of
% this loop, if this is a 'NEW CASE'; otherwise the previously
% iterated value of y is brought in.

y_old = y;
b = -e*V^2*w/(2*E*I);
c = 0.65*b/w;

%NEWTON ITERATIONS FOR PULLIN

r(half) = 0.5;
err = 1;
count = 1;
lim = 2*floor((log(VMAX/(VMAX-VMIN)))) + Z;
while ((abs(r(half)) < err & count <= lim) | count < Z)
  if (count == 2)
      err = abs(r(half));
  end
  for i=3:N-2
	y2p(i)=(y(i-1)-2*y(i)+y(i+1))/h^2;
	y4p(i)=(y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
	A(i,i-2)=1;
	A(i,i-1)=-4-a*h^2;
	A(i,i)=6+2*a*h^2+2*h^4/y(i)*(y4p(i)-a*y2p(i)-c/(2*y(i)));
	A(i,i+1)=-4-a*h^2;
	A(i,i+2)=1;
	C(i)=h^4*(b/y(i)^2-y4p(i)+a*y2p(i)+c/y(i));
  end
  if (step ~= 0)
    phi = (-MoEI + (y(4)-2*y(3)+y(2))/(h^2)) * step;
    y1p3 = (y(4)-y(2))/(2*h);
    C(2) = 2*h*(phi-y1p3);
  end
  r = A\C;
  if count > 10
    y = y+(r/5);
  else
    y = y+r;
  end
  count = count + 1;
end

%END NEWTON ITERATIONS FOR PULLIN

V_max = VMAX
V_min = VMIN
iterations = count-1
V_pull_in = V
Center_gap = y(half)*1e6

if ((abs(r(half)) < 1e-4*g) & (y(half) < 2*g) & (y(half) > 0.575*g))
  VMIN = V;
  V = (VMAX+4*VMIN)/5;
else
  VMAX = V;
  V = (VMAX+4*VMIN)/5;
  y = y_old;
end

end

%END BINARY SEARCH

end

plot(y*1e6);
ymax = 1.05*max(y(half),y(3));
ymin = 0.95*min(y(half),y(3));
axis([3 half ymin*1e6 ymax*1e6]);
xlabel('Discretization in x');
ylabel('Gap (um)');