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