ltd13002 / ME3255_Final_Project

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

ltd13002 / ME3255_Final_Project

README.md

ME3255_Final_Project

Part A

Problem

Approach

Part B

Problem

Part C

Problem

Approach

Part D

Problem

Part E

Problem

Approach

Part F

Problem

Approach

Part G

Problem

Approach

Part H

Problem

Approach

About

Releases

Contributors 3

Languages

ltd13002 / ME3255_Final_Project

Clone with HTTPS

Launching GitHub Desktop

Launching GitHub Desktop

Launching Xcode

Launching Visual Studio

Latest commit

Git stats

Files

README.md

ME3255_Final_Project

Part A

Problem

Approach

Part B

Problem

Part C

Problem

Approach

Part D

Problem

Part E

Problem

Approach

Part F

Problem

Approach

Part G

Problem

Approach

Part H

Problem

Approach

About

Resources

Releases

Contributors 3

Languages