Part H

Part H/Rel_error.m

@@ -0,0 +1,12 @@
+function re = Rel_error (T)
+% Rel_error: calculates relative error of a vector
+% re = Rel_error (T)
+% input:
+% T = vector of numbers
+% output:
+% re = relative error of vector
+
+re = zeros(1,length(T)-1);
+for i = 2:length(T)
+    re(i-1)= abs(T(i)-T(i-1))/T(i-1);
+end

Part H/SE_diff.m

@@ -0,0 +1,38 @@
+function [pw_se,w]=SE_diff(T,P,n)
+% SE_diff: calculates difference between strain energy and work done by pressure in
+% membrane
+% [pw_se,w]=SE_diff(T,P,n)
+% input:
+% T = Tension (microNewton/micrometer)
+% P = Pressure (MPa)
+% n = number of rows and columns of interior nodes
+% output:
+% pw_se = difference between strain energy and work done by pressure in
+% membrane
+% w = vector of displacement of interior nodes
+
+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 H/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 H/membrane_solution.m

@@ -0,0 +1,30 @@
+function [w] = membrane_solution(T,P,n)
+% membrane_solution: dispalacement of node for membrane with nxn interior nodes
+% [w] = membrane_solution(T,P,n)
+% input:
+% T = Tension (microNewton/micrometer)
+% P = Pressure (MPa)
+% n = number of rows and columns of interior nodes
+% output:
+% w = vector of displacement 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 H/part_h.asv

@@ -0,0 +1,9 @@
+P=0.001;
+n = [20:5:40];
+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

Part H/part_h.m

@@ -0,0 +1,13 @@
+P=0.001;
+n = 5:5:45;
+a = zeros(1,length(n));
+for i = 1:length(n)
+    T = tension_sol(P,n(i));
+    w = membrane_solution(T,P,n(i));
+    wmax = max(w);
+    x = wmax';
+    y = P';
+    Z = x.^3;
+    a(i) = Z\y;
+end
+Rel_error(a)

Part H/setDefaults.m

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

Part H/tension_sol.m

@@ -0,0 +1,13 @@
+function [T,ea] = tension_sol(P,n)
+% tension_sol: outputs tension of a membrane given the pressure and number
+% of nodes
+% [T,ea] = tension_sol(P,n)
+% input:
+% P = Pressure (MPa)
+% n = number of rows and columns of interior nodes
+% output:
+% T = Tension (microNewton/micrometer)
+% ea = approximate relative error (%)
+
+y =@(T) SE_diff(T,P,n);
+[T,fx,ea,iter]=bisect(y,.01,1);