SE_diff.m

function [pw_se,w] = SE_diff(T,P,n)
% function that calculates the difference between
% strain energy and work done by pressure on the membrane.
% Input:
% T = tension per unit length (uN/um)
% P = pressure (MPa)
% n = # of interior node rows/columns
% Output:
% pw_se = Absolute value of difference between strain energy and work done by pressure
% w = displacement vector for interior nodes

E = 1e6;        % 1 TPa ~= 10^6 MPa
t = 3*10^-4;    % thickness [um]
h = 10/(n+1);   % height [um]
v = 0.31;       % Poisson's Ratio

% Displacement vector w found using Part C
w = membrane_solution(T,P,n);
z = zeros(n + 2);
z(2:end-1,2:end-1) = reshape(w,[n n]);

% Calculate average displacement, wavg, for each element by taking the displacement at each
% corner and then average the found values.
num = n + 1;
wavg = zeros(num);
for i = 1:num
  for j = 1:num
    wavg(i,j) = mean([z(i,j),z(i+1,j),z(i,j+1),z(i+1,j+1)]);
  end
end

% final work done by pressure
pw = sum(sum(wavg.*h^2.*P))

% to find= change in displacement, find the change in displacement on
% the x-axis, dwdx, and the change in displacement on the y-axis, dwdy, and
% average the found values.
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

% Using dwdx and dwdy, calculate the strain energy, se.
se = (E*t*h^2)/(2*(1-v^2)) * sum(sum((1/4).*dwdx.^4+(1/4).*dwdy.^4+(1/4).*(dwdx.*dwdy).^2));

% Final value of difference between strain energy and work done by pressure
pw_se = pw - se;