ME3255_Final_Project
Part A
function [w] = membrane_solution3(T,P)
% T = Tension (microNewton/micrometer)
% P = Pressure (MPa)
od = ones(8,1);
od(3:3:end) = 0;
k = -4*diag(ones(9,1))+diag(ones(9-3,1),3)+diag(ones(9-3,1),-3)+diag(od,1)+diag(od,-1);
y = -(10/4)^2*(P/T)*ones(9,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/4:10,0:10/4:10);
z = zeros(size(x));
z(2:end-1,2:end-1) = reshape(w,[3 3]);
surf(x,y,z)
title('Membrane Displacement')
zlabel('Displacement (micrometer)')
% Membrane displacement is shown on chart
end
Part B
% Part B Script
[w] = membrane_solution3(0.006,0.001);
Part C
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 D
% Part D Script
[w] = membrane_solution(0.006,0.001,10)
Part E
function [pw_se,w]=SE_diff(T,P,n)
E = 1; %TPa Units may need to be changed
v = .31; %Poissons ratio
t = .3; %nm
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)]);
dwdy(i,j) = mean([z(i,j+1)-z(i,j),z(i+1,j+1)-z(i+1,j)]);
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 F
n=[3,20:5:40];
P=0.001; %MPa
T = zeros(1,length(n));
ea = zeros(1,length(n));
for i = 1:length(n)
[T(i), ea(i)] = tension_sol(P,n(i));
end
function [T,ea] = tension_sol(P,n)
y =@(T) SE_diff(T,P,n);
[T,fx,ea,iter]=bisect(y,.01,1);
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{:});
function re = Rel_error (T)
re = zeros(1,length(T)-1);
for i = 2:length(T)
re(i-1)= abs(T(i)-T(i-1))/T(i-1);
end
number of nodes | Tension (uN/um) | rel. error |
---|---|---|
3 | 0.0489 | n/a |
20 | 0.0599 | 22.6% |
25 | 0.0601 | 0.27% |
30 | 0.0602 | 0.15% |
35 | 0.0602 | 0.09% |
40 | 0.0603 | 0.06% |
Part G
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
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.^3)
title('Pressure vs Maximum Deflection')
xlabel('Maximum Deflection (um)')
ylabel('Pressure (MPa)')
print('Part g','-dpng')
![part g](./Part G/Part g.png)