Axially loaded stepped shaft analysis in MATLAB

Question

I have a stepped shaft as per the attached image. Following information available as an input parameters:

1. Young's modulus 123e3N/mm^2.
2. Cross-sectional area 300mm^2 for the length of 400mm
3. Cross-sectional area 400mm^2 for the length of 250mm
4. Axial force of 200kN acts axially on the shaft and the location of load is at 200mm from the one end of the shaft on cross-sectional area of 300mm^2

I need help to make do finite element analysis in MATALB. Please help me in making MATLAB code for this.

Answer 1

%% Clearing workspace

clc
clear
close all

%% Element specifications 

ne = 3; % Number of elements 
nne = 2; % Number of nodes per element 
nn = ne*(nne - 1) + 1; % total number of nodes
ndof = 1; % Number of degress of freedom per node 
sg = nn*ndof; % size of global stiffness matrix 
se = nne*ndof; % size of elemental stiffness matrix 
KG = zeros(sg,sg); % Global stiffness matrix 
Ke = zeros(se,se); % Elemental Stiffness MAtrix 
Fe = zeros(se,1); % Elemental Force Vector 
FG = zeros(sg,1); % Global Force Vector

%% Geometrical parameters

E = 123e3*ones(1,ne); % Young's Modulus in N/mm^2
P = 200e3; % Force in N
F = P;

A = ones(1,ne) ; % Area of cross-section 
A(1)=300; % Area of cross-section of 1st element in mm^2
A(2)=300; % Area of cross-section of 2nd element in mm^2
A(3)=400; % Area of cross-section of 3rd element in mm^2

L = ones(1,ne); % Length of elements in mm
L(1)=200; % Length of 1st element in mm
L(2)=200; % Length of 2nd element in mm
L(3)=250; % Length of 3rd element in mm

%% Assembly of Global Stiffness Matrix

for i = 1:ne
    Ke = (A(i)*E(i)/L(i))*[1 -1;-1 1]; % Element Stiffness Matrix
    for j = 1:se
        for k = 1:se
            KG(i + j - 1, i + k - 1) = KG(i + j - 1, i + k - 1) + Ke(j,k);
        end
    end
end

%% Concentrated Load Vector at end

FG(2,1) = F; % Defining location of concentrated load

%% Application of boundary conditions

KGS = KG;
cdof = [1 4]; % specify fixed degree of freedom number
Lcdof = length(cdof);
for a = 1:Lcdof
    KGS(cdof(a),:) = 0;
    KGS(:,cdof(a)) = 1;
    FG(cdof(a),1) = 0;
end
FGL = length(FG); 
for b = 1:FGL
    if(b > length(FG))
        elseif(FG(b)<0)
            FG(b) = [];
    end
end

%% Solving for displacement

U = linsolve(KGS,FG) 
U1=KGS\FG 

%% Calculation of Reaction Forces 

FR = KG*U1