.

Part G/SE_diff.m

@@ -0,0 +1,26 @@
+function [pw_se,w]=SE_diff(T,P,n)
+E = 1000000; %TPa Units may need to be changed
+v = .31; %Poissons ratio
+t = .0003; %um
+h = 10/(n+1); %nm
+w = membrane_solution(T,P,n);
+z = zeros(n+2);
+z(2:end-1,2:end-1) = reshape(w,[n n]);
+num = n + 1;
+wbar = zeros(num);
+for i = 1:num
+    for j = 1:num
+        wbar(i,j) = mean([z(i,j),z(i+1,j),z(i,j+1),z(i+1,j+1)]);
+    end
+end
+pw = sum(sum(wbar.*h^2.*P));
+dwdx = zeros(num);
+dwdy = zeros(num);
+for i = 1:num
+    for j = 1:num
+        dwdx(i,j) = mean([z(i+1,j)-z(i,j),z(i+1,j+1)-z(i,j+1)])./h;
+        dwdy(i,j) = mean([z(i,j+1)-z(i,j),z(i+1,j+1)-z(i+1,j)])./h;
+    end
+end
+se = E*t*h^2/(2*(1-v^2))*sum(sum(0.25.*dwdx.^4+.25.*dwdy.^4+0.5.*(dwdx.*dwdy).^2));
+pw_se = pw-se;

Part G/bisect.m

@@ -0,0 +1,37 @@
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+  xrold = xr;
+  xr = (xl + xu)/2;
+  iter = iter + 1;
+  if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+  test = func(xl,varargin{:})*func(xr,varargin{:});
+  if test < 0
+    xu = xr;
+  elseif test > 0
+    xl = xr;
+  else
+    ea = 0;
+  end
+  if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});

Part G/membrane_solution.m

@@ -0,0 +1,25 @@
+function [w] = membrane_solution(T,P,n)
+    % T = Tension (microNewton/micrometer)
+    % P = Pressure (MPa)
+    % n = # of interior nodes
+
+    od = ones(n^2-1,1);
+    od(n:n:end) = 0;
+    k = -4*diag(ones(n^2,1))+diag(ones((n^2)-n,1),n)+diag(ones((n^2)-n,1),-n)+diag(od,1)+diag(od,-1);
+
+    y = -(10/(n+1))^2*(P/T)*ones(n^2,1);
+    w = k\y;
+
+    % Solves for displacement (micrometers)
+    % Output w is a vector
+    % Solution represents a 2D data set w(x,y)
+
+    [x,y] = meshgrid(0:10/(n+1):10,0:10/(n+1):10);
+    z = zeros(size(x));
+    z(2:end-1,2:end-1) = reshape(w,[n n]);
+    surf(x,y,z)
+    title('Membrane Displacement')
+    zlabel('Displacement (micrometer)')
+
+    % Membrane displacement is shown on chart
+end

Part G/part_g.asv

@@ -1,5 +1,20 @@
 P = linspace(.001,.01,10);
 n = 20;
+T = zeros(1,length(P));
+wmax = zeros(1,length(P));
 for i = 1:length(P)
     T(i) = tension_sol(P(i),n);
-
+    w = membrane_solution(T(i),P(i),n);
+    wmax(i) = max(w);
+end
+clf
+setDefaults
+plot(wmax,P,'o')
+%P_out = interp1(wmax,P,linspace(.1,.3),'pchip');
+hold on 
+%plot(linspace(.1,.3),P_out)
+x = wmax;
+y = P;
+Z = x.^3;
+a = Z/
+hold off

Part G/part_g.m

@@ -7,4 +7,11 @@ for i = 1:length(P)
     w = membrane_solution(T(i),P(i),n);
     wmax(i) = max(w);
 end
-plot(wmax,P,'o')
+clf
+setDefaults
+x = wmax;
+y = P';
+Z = [x.^3]';
+a = Z\y;
+x_fcn = linspace(min(x),max(x));
+plot(x,y,'o',x_fcn,a.*x_fcn)

Part G/setDefaults.m

@@ -0,0 +1,3 @@
+set (0, 'defaultaxesfontsize', 18)
+set (0, 'defaulttextfontsize', 18)
+set (0, 'defaultlinelinewidth', 4)

Part G/tension_sol.m

@@ -0,0 +1,3 @@
+function [T,ea] = tension_sol(P,n)
+y =@(T) SE_diff(T,P,n);
+[T,fx,ea,iter]=bisect(y,.01,1);