# 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. ```matlab 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. ```matlab % Part B Script % From problem statement: % T=0.006 uN/um % P=0.001 MPa [w] = membrane_solution3(0.006,0.001); ``` ![](https://github.uconn.edu/ltd13002/ME3255_Final_Project/blob/master/Part%20B/PartBFigure.png) # 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. ```matlab 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. ```matlab % 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) ``` ![](https://github.uconn.edu/ltd13002/ME3255_Final_Project/blob/master/Part%20D/PartDFigure.png) # 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)'. ```matlab 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. ```matlab % 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 ``` ```matlab 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); ``` ```matlab 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{:}); ``` ```matlab 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. ```matlab % 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') ``` ![](https://github.uconn.edu/ltd13002/ME3255_Final_Project/blob/master/Part%20G/Part%20g.png) ### 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. ```matlab % 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