07_final_project

Part A

Problem Statement

Create a function called membrane_solution3 that uses the central finite difference approximation to solve for the displacements of a membrane with 3 x 3 interior nodes. The function should take the pressure (P) and tension (T) as inputs and output a column vector, w.

Approach

Using the central difference approximation, the matrix K is created using the coefficients from the diagonal and off-diagonals and the matrix y is a column vector of ones multiplied by h^2*(P/T). The vector w is found by using the backslash key, which is a built-in Matlab function for matrix division.

function [w] = membrane_solution3(T,P) % Set up initial matrix
    % This function outputs a central finite difference approximation
    % The inputes include the following:
    % T = given tension (microNewton/micrometer)
    % P = given pressure (MPa)
    % The output is a vector w

    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);

    % Solve for unknown matrix, w
    y = -(10/4)^2*(P/T)*ones(9,1); %output vector
    w = k\y; %solves for w in microm

    % The following will create a mesh grid for a given pressure and
    % tension
    [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)
end

Part B

Problem Statement

Using the function from part A as well as a given tension of 0.006 uN/um and a given pressure of 0.001 MPa, solve for w and plot the result in 3 dimensions

Approach

The function created in part A is called using the values from the problem statement. Surf(x,y,z) plots the surface of the membrane in 3 dimensions.

[w] = membrane_solution3(0.006,0.001);

Part C

Problem Statement

Create a function called membrane_solution that uses the central finite difference approximation to solve for the displacements of a membrane with n x n interior nodes. The function should take the pressure (P), tension (T), and number of interior nodes (n) as inputs and output the column vector, w.

Approach

Similar to part A, the central difference approximation is used to create the K matrix. However, this time the dimensions of the matrix are dependent on the number of interior nodes, which is an input to the function. The w matrix is once again solved using the backslash key.

function [w] = membrane_solution(T,P,n)% Set up initial matrix
    % This functions outputs a central finite differnce approximation
    % for n by n interior nodes
    % The inputs are the following:
    % T = given tension (microNewton/micrometer)
    % P = given pressure (MPa)
    % n = number of interior nodes
    % The output is a vector w that represents the membrane solution for n
    % by n 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);

    % Solve for unknown matrix, w
    y = -(10/(n+1))^2*(P/T)*ones(n^2,1); %output vector
    w = k\y; %solves for w in micron

    % The following creates a mesh grid for n by n interior nodes
    [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)
end

Part D

Problem Statement

Solve for the vector w using the function from part C with a P value of 0.001 MPa, T value of 0.006 un/um, and 10 interior nodes. The 3 dimensional coordinates will then be plotted.

Approach

The function written in part C is called using the input values seem in the problem statement. The 3D surface plot of the membrane is created.

[w] = membrane_solution(0.006,0.001,10)

Part E

Problem Statement

Create a function that calculates the difference between the strain energy and the work done by pressure for a membrane with n x n elements. The function should output a single value for this difference (pw_se) as well as the w vector from part C.

Approach

Using the output vector from part C (w), a matrix of the interior and exterior nodes is created. From here, dx and dy matrices are made, which are found by taking the difference between the columns and rows of the initial matrix. From here, the dw_dx and dw_dy matrices are created and plugged into the equation given in the project description. From here, the difference between the right and left sides of the equation could be found.

function [pw_se,w] = SE_diff(T,P,n)
    % This function calculates the strain energy and work done by pressure
    % for n by n elements
    % The w vector from the function "membrane_solution" is called and used
    % to solve for the average w across the each element
    % Call function from part c to obtain vector w
    w = membrane_solution(T,P,n);
    % Create All Nodes Matrix (Interior & Exterior)
    m = zeros(n+2);
    c = vec2mat(w,n); % Vector w rearranged to n x n matrix.
    m(2:n+1,2:n+1) = c; % insert matrix c into middle of m matrix
    % m = symmetric matrix

    h = 10/(n+1); % h = delta_x = delta_y (micron)

    % Compute w bar
    w_bar = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
            w_bar(i,j) = (m(i,j)+m(i,j+1)+m(i+1,j)+m(i+1,j+1))/4;
        end
    end

    % Left side of the equation
    left_side = zeros(n+1,n+1);
    for  i = 1:(n+1)^2
        left_side(i) = P * h^2*w_bar(i);
    end
    L = sum(sum(left_side));

    % Compute dw_dx
    dx = zeros(n+2,n+1);
    for i = 2:n+2
        dx(:,i-1)= (m(:,i)-m(:,i-1))/(h);
    end

    dw_dx = zeros(n+1,n+1);
    for i = 1:n+1 % rows
        for j = 1:n+1 % columns
        dw_dx(i,j) = (dx(i,j)+dx(i+1,j))/2;
        end
    end

    % Compute dw_dy
    dy = zeros(n+1,n+2);
    for i = 2:n+2
        dy(i-1,:)= (m(i-1,:)-m(i,:))/(h);
    end

    dw_dy = zeros(n+1,n+1);
    for i = 1:n+1 % rows
        for j = 1:n+1 % columns
        dw_dy(i,j) = (dy(i,j)+dy(i,j+1))/2;
        end
    end

    right_side = zeros(n+1,n+1);

    for i = 1:(n+1)^2
        right_side(i) = (0.25*(dw_dx(i))^4+0.25*(dw_dy(i))^4+0.5*(dw_dx(i)*dw_dy(i))^2);
    end
    % Compute the strain energy
    R = sum(sum(right_side));
    R = (1000*10^3*0.3*10^-3*h^2)/(2*(1-0.31^2))*R;

    % Compute the difference between the strain energy and the work done by
    % the pressure
    pw_se = R - L;
end

• At 10 nodes, pw_se = 53.4902 pJ

Part F

Problem Statement

Using a root finding method, calculate the tension in the membrane based on a given pressure of 0.001 MPa and testing every 5 nodes from 20 to 40. Display the results in a table and show that the error in the tension is decreasing as the number of nodes is increased.

Approach

Using the bisect method, the tension in the membrane based on the conditions given in the problem statement is found. The relative error is found by calculating the difference between the error from the initial number of nodes to all of the subsequent node values.

% This script calls for the result of the SE_diff function and uses the
% bisect root finding method to solve for the tension given a pressure (P)
% and desired range of nodes (n)
[pw_se,w] = SE_diff(T,P,n);
[root,fx,ea,iter] = bisect_final_project(@(T) SE_diff(T,0.001,40),.001,1,.1)ff(T,0.001,25),.001,1,.1)

• The same code was run five times, changing the value of n from 20 to 40 for intervals of 5.

number of nodes	Tension (uN/um)	rel. error
3	0.0489	n/a
20	0.0599	18.4%
25	0.0601	0.30%
30	0.0602	0.17%
35	0.0603	0.16%
40	0.0603	0.00%

Part G

Problem Statement

For a given pressure range of 0.001 to 0.01 MPa with 10 steps, determine the membrane tension using a root finding method at each pressure. Plot the results of pressure vs. maximum deflection and use a cubic best-fit line to determine the coefficient A, which is the best-fit constant for the line.

Approach

Using the given pressure range, the tension is found at each different pressure value using the bisect method. From here, the deflection values at the calculated tensions and given pressures are found by calling the membrane_solution function from part C. The best fit line and coefficient are found using the built-in Matlab functions polyfit and polyval. These return a 4 x 1 matrix of coefficients, but since there is only one x term in the function (x^3), the first value in the matrix is determined to be the coefficient A.

% This creates a plot for the max_w vs a given pressure
clear
P = linspace(0.001,0.01,10); % Assign the range of pressures used to find T
% In order to find the tension the bisect method is used
for i = 1:length(P)
    func = @(T) SE_diff(T,P(i),10);
    [root(i),fx,ea,iter] = bisect_final_project(func,0.001,1,.1);
end
% Each value of w is calculated using each root and pressure value
for i = 1:length(root)
    w = membrane_solution(root(i),P(i),10);
    w1(:,i) = w;
end
% in order to get the w_max we take the maximum w value from each column of the w vector
w_max = max(w1);
coefficients = polyfit(w_max,P,3);
Y = polyval(coefficients,w_max);
% plot the w_max vs. the Pressure and include a cubic best fit curve.
plot(w_max,P,w_max,Y,'or')
xlabel('Max Deflection (micron)')
ylabel('Pressure (MPa)')
title('Pressure vs. Maximum Deflection')

Output - Based on 10 interior nodes:

• A = 0.4443

Part H - Extra Credit

Problem Statement

Show that the constant A is converging as the number of interior nodes is increased and display it in a table similar to the one seen in part F.

Approach

Running the code in Part G while increasing the number of nodes by 5 each time will yield the A constant. The relative error is found by calculating the difference between the error from the initial number of nodes to all of the subsequent node values.

number of nodes	Value of A	rel. error
3	0.3510	N/A
20	0.5622	37.57%
25	0.6422	12.46%
30	0.4932	30.21%
35	0.5366	8.09%
40	0.5339	0.51%