diff --git a/Part H/Rel_error.m b/Part H/Rel_error.m new file mode 100644 index 0000000..5e17536 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Part H/SE_diff.m b/Part H/SE_diff.m new file mode 100644 index 0000000..1eaddea --- /dev/null +++ b/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; \ No newline at end of file diff --git a/Part H/bisect.m b/Part H/bisect.m new file mode 100644 index 0000000..47ee7ec --- /dev/null +++ b/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{:}); \ No newline at end of file diff --git a/Part H/membrane_solution.m b/Part H/membrane_solution.m new file mode 100644 index 0000000..af4d17d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Part H/part_h.asv b/Part H/part_h.asv new file mode 100644 index 0000000..d344d5f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Part H/part_h.m b/Part H/part_h.m new file mode 100644 index 0000000..30fa740 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/Part H/setDefaults.m b/Part H/setDefaults.m new file mode 100644 index 0000000..51e7ee3 --- /dev/null +++ b/Part H/setDefaults.m @@ -0,0 +1,3 @@ +set (0, 'defaultaxesfontsize', 18) +set (0, 'defaulttextfontsize', 18) +set (0, 'defaultlinelinewidth', 4) \ No newline at end of file diff --git a/Part H/tension_sol.m b/Part H/tension_sol.m new file mode 100644 index 0000000..b0f0132 --- /dev/null +++ b/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);