Skip to content
72c658286e
Go to file
Code

Latest commit

 

Git stats

Files

Permalink
Failed to load latest commit information.
Type
Name
Latest commit message
Commit time
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

README.md

ME3255_Final_Project

Part A

Problem

Create a central finite difference approximation function to solve for the membrane displacement gradient with 3-by-3 interior nodes '[w]=membrane_solution3(T,P)'. The inputs are tension (T) and pressure (P), and the output is column vector w.

function [w] = membrane_solution3(T,P)
% membrane_solution3: dispalacement of node for membrane with 3x3 interior
% nodes
% [w] = membrane_solution3(T,P)
% input:
% T = Tension (microNewton/micrometer)
% P = Pressure (MPa)
% output:
% w = vector of displacement of interior nodes

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)
% 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 (3 x 3)')
    zlabel('Z Position (micron)')
    ylabel('Y Position (micron)')
    xlabel('X Position (micron)')

% Membrane displacement is shown on chart

end
    
end 

Approach

  • A matrix is set up using finite element analysis of the interior nodes of a 5x5 matrix
  • A vector for the y direction is then set up for the membrane
  • Using linear algebra, a vector for the displacement of the nodes is created
  • The vector can then be transformed to represent the actual surface as a matrix

Part B

Problem

Use membranesolution3 function to solve for w and plot the result with surf(X,Y,W). Pressure P=0.001 MPa and tension T=0.006 uN/um.

% Part B Script
% From problem statement:
% T=0.006 uN/um
% P=0.001 MPa
[w] = membrane_solution3(0.006,0.001);

Part C

Problem

Create a central finite difference approximation function to solve for the membrane displacement gradient with n-by-n interior nodes '[w]=membrane_solution(T,P,n)'. The inputs are tension (T), pressure (P), and number of nodes (n), and the output is column vector w.

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

Approach

  • This problem uses the same steps as with the membrane_solution3 function, except it allows for the user to change the number of interior nodes
  • This required alterations of the all the indexing and sizing in the function along with a few calculations

Part D

Problem

Use membranesolution function to solve for w and plot the result with surf(X,Y,W). Pressure P=0.001 MPa, tension T=0.006 uN/um, and number of nodes n=10.

% Part D Script
% From problem statement:
% T=0.006 uN/um
% P=0.001 MPa
% n = 10 nodes
[w] = membrane_solution(0.006,0.001,10)

Part E

Problem

Create a function that calculates the difference in strain energy and work done by pressure for n-by-n elements '[pw_se,w]=SE_diff(T,P,n)'.

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 = 1e6; %MPa Units may need to be changed
v = .31; %Poissons ratio
t = .0003; %um
h = 10/(n+1); %um
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;

Approach

  • Using the membrane_solution function, a vector of the displacements is formed
  • Next, the average displacement for each element is calculated. For each elements, the dispalcement at all four courners is taken and then averaged
  • Using these values, the work done by pressure can be calculated
  • For the change in dispalcement over the x and y coordinate system, the values of the change on the left and right (y-axis) or top and bottom (x-axis) are taken and averaged
  • Using these values, the strain energy can be calculated

Part F

Problem

Use a root-finding method to calculate the tension in the membrane given a pressure, P=0.001 MPa, and n=[20:5:40] interior nodes. Show the decrease in tension error in a table.

% Part F script
% From Problem Statement:
 n=[3,20:5:40];
 P=0.001; %MPa
 
 % Sets vector length so it doesn't change every iteration in for loop
 T = zeros(1,length(n));
 ea = zeros(1,length(n));
 
 % Uses tension_sol function to solve for the tension for each input
 for i = 1:length(n)
[T(i), ea(i)] = tension_sol(P,n(i));
 end
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);
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)
 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
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%

Approach

  • This problem uses the bisect method for locating roots and the SE_diff function for calculating the difference in work and strain energy of the membrane
  • The script runs through all iterations of different amounts of nodes, zeroing the SE_diff function output and saving the values for tension in the T variable as a vector

Part G

Problem

Plot the pressure vs maximum deflection for pressure range 0.001-0.01 MPa with 10 steps using a root-finding method to determine tension, T, at each pressure.

% Part G Script
% From Problem Statement
P = linspace(.001,.01,10);
n = 20;

% Sets vector length so it doesn't change every iteration in for loop
T = zeros(1,length(P));
wmax = zeros(1,length(P));

% Uses tension_sol and membrane_solution functions to solve for the maximum displacement for all iterations
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

% Sets up figure for plot
clf
setDefaults

% Standard Linear Regression with equation P(x) = A*w^2
x = wmax';
y = P';
Z=x.^3;
a=Z\y;
x_fcn=linspace(min(x),max(x));

% Plots regression line with actual points
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')

Approach

  • This script uses the tension_sol function as described in the part above, running though all iterations of pressure
  • Additionally, the max deflection for each pressure and tension is calculated at each iteration
  • From there, a general linear regression is calculated with the formula P(x) = A*w^3
  • The results are plotted

Part H

Problem

Show that the constant A is converging as the number of nodes is increased.

% Part G Script
% From Problem Statement
P=0.001;
n = 5:5:50;

% Sets vector length so it doesn't change every iteration in for loop
a = zeros(1,length(n));

% Uses tension_sol and membrane_solution functions along with a regression line to solve for the constant A in the regression line
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

% Calculates the relative error between the constant A values
Rel_error(a)
number of nodes vaalue of const. A rel. error
5 0.4423 n/a
10 0.5389 21.84%
15 0.5349 0.73%
20 0.5483 2.5%
25 0.5454 0.52%
30 0.5503 0.89%
35 0.5485 0.31%
40 0.5510 0.45%
45 0.5499 0.2%

Approach

  • Part H uses the same basis as Part G, but with varrying n (number of nodes) valeus instead
  • Additionally, the outputs and a few calculations were deleted to speed up the process since they were not needed here
  • In the chart, while there are some fluctuations, there is a clear decrease in the relative error of A, leading to the assumption that it is converging as the number of nodes increases

About

ME 3255 Final Project

Resources

Releases

No releases published

Languages

You can’t perform that action at this time.