Hi David,
There are a few issues in the code you have provided:
- The plotting command plot(u(n,1)) is incorrect because it will just plot a single value rather than a series of values. To plot the axial deformation at the free-end as a function of the number of elements, you need to collect the last displacement value (u(end)) for each iteration over a range of element counts.
- The code will need a for loop to iterate through different numbers of elements and store each corresponding displacement in an array. The plot needs to be generated outside of the loop.
- The code needs to store the displacement of the free-end for each number of elements in a separate array to be used for plotting.
Here is the modified MATLAB code:
clc;
clear;
%============================= Inputs ===================================
elements = 10; % Number of elements desired, nodes at left and right ends of element
Lo = 32; % Length of Tapered Bar
Ao = 12; % Initial Cross-sectional Area at X=0
E = 11*(10^3); % Modulus of Elasticity
F_applied = -2.5; % Applied force at the free-end
% Arrays to store results for plotting
element_counts = 1:elements;
free_end_displacements = zeros(1, elements);
for e = 1:elements
X = Lo/e; % Length of individual element
n = e + 1;
F = zeros(n,1);
F(n) = F_applied; % Force Vector Matrix
%============================= Global Stiffness =========================
ks = [1,-1; -1,1]; % Local Stiffness matrix
K = zeros(n, n);
A = zeros(1, e); % Initialize cross-sectional area array
for i = 1:e
A(i) = (Ao*(1-0.5*(i-1)*X/Lo) + Ao*(1-0.5*i*X/Lo)) / 2; % Generating avg. cross-section area of each element
end
k = A*E/X; % Generating local stiffness matrices k
for i = 1:e
K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + k(i)*ks; % Generating Global Stiffness Matrix K
end
fixed_nodes = 1; % Fixed Node at left end of bar
free_nodes = setxor(1:n, fixed_nodes);
Kred = K(free_nodes, free_nodes); % Reduced Global Stiffness Matrix after partitioning
Fred = F(free_nodes, 1); % Reduced Force Matrix after partitioning
%============================= Solve Displacements =====================
ured = Kred\Fred; % Reduced vector matrix after partitioning
u = zeros(n, 1); % displacement vector matrix
u(free_nodes) = ured; % displacement of free-end is last value in matrix
% Store the displacement of the free-end for the current number of elements
free_end_displacements(e) = u(end);
end
% Now plot the axial deformation at the free-end vs. the number of elements
plot(element_counts, free_end_displacements);
grid on;
title('Axial Deformation vs. Number of Elements');
xlabel('Number of Elements');
ylabel('Axial Deformation (inches)');
Hope this helps!
