Help generating a plot that shows the axial deformation at the free-end as a function of the number of elements, using Finite Element Analysis

Hi, I need help generating a plot in MATLAB that shows the axial deformation at the free-end of a tapered bar as a function of the number of elements, with the matlab code automating this and not hardwiring in any values. Below is the example bar and my code for solving the displacement at each node. I am a complete novice when it comes to Finite Element Analysis and MATLAB, but I know the code works and is correct for the determination of displacements at each node with any number of elements as it converges to the exact solution of -0.000840178 inches at the free end (last node) as the number of elements increase. I've tried for hours, but haven't got anything to work. I'm assuming it must need some "if" and/or "for" statement to iterate each free-end displacement for any number of elements from 1:e, but I don't know how to generate this and plot each free-end displacement with its relative number of elements used. Thank you for your help.

clc
 
%============================= Inputs ===================================
e = 10; %Number of elements desired, nodes at left and right ends of element
Lo=32; %Length of Tapered Bar
X=Lo/e; %Length of individual element
Ao=12; %Initial Cross-sectional Area at X=0


E = 11*(10^3); %Modulus of Elasticity


n = e + 1; %Number of Nodes
F = zeros(n,1); 
F(n) = -2.5; %Force Vector Matrix

%============================= Global Stiffness =========================
ks = [1,-1; -1,1]; %Local Stiffness matrix
K = zeros(n, n);
 
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


for i=1:e
    k(i) = A(i)*E/X; %Generating local stiffness matirces k
end


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 matix after partitioning
u=zeros(n,1); %displacment vector matirx
u(free_nodes)=ured %displacement of free-end is last value in matrix


plot(u(n,1)) %Incorrect plot where y-axis is free-end deformation and x-axis is number of elements

    title('Axial Deformation vs. Number of Elements')
    xlabel('Number of Elements')
    ylabel('Axial Deformation')

Answer 1

[Jaynik]

Hi David,

There are a few issues in the code you have provided:

1. 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.
2. 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.
3. 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!