GLobal stiffness matrix in Matlab for my structure

axes
clear
clc
close all

format short g

ndofn = 3; % number of degrees of freedom per node

% each number corresponds to a force for the corresponding DOF
Force = [0 0 0  0 0 0  0 0 0  0 -1000 0  0 -1000 0  0 -1000 0  0-1000 0  0 0 0  0 0 0  0 0 0]';

% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2  pi/2  0  0  0  0  0  0    -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2  0     0  0  0  0  0 -pi/2 -pi/2];

% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0  -240 0   0   0   0   0  -240 0]; %in

% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1  B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9;
     2 3 4 5 6 8 9 10];

%for straight
syms L x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
       0 1 0;
       0 240 1];
HxB_T = HxB'; % transpose

% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
    0 1/(G*Ay) 0;
    0 0 1/(E*Iz)];

fux_straight = HxB_T*fuq_straight*HxB;

fBA_straight = int(fux_straight, x, 0, L);

% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in

% Stiffness matrix of cantilever straight element
kB_straight = eval(inv(fBA_straight));

%for columns
syms L1 x A1 E G Ay1 Iz1 real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
        0 1 0;
        0 240 1];
HxB_Tc = HxBc'; % transpose

% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A1*E) 0 0;
             0 1/(G*Ay1) 0;
             0 0 1/(E*Iz1)];

fux_column = HxB_Tc*fuq_column*HxBc;

fBA_column = int(fux_column, x, 0, L1);

% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in

% Stiffness matrix of cantilever straight element
kB_column = eval(inv(fBA_column));

%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);

% rotation operator
lambda_QB_T = [c s 0;
              -s c 0;
              0 0 1];
lambda_QB = lambda_QB_T';

% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB =   [1 0 0;
         0 1 0;
      -240 240 1];
HQB_T = HQB';

% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
        0 1/(G*Ay) 0;
        0 0 1/(E*Iz)];

TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member

theta = pi/2;
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);

% Stiffness matrix of cantilever curved element
kB_curved = eval(inv(fBA_curved));

%global matrice
nnodes = length(t)+1;

% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn,nnodes*ndofn);
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn,2*ndofn);

%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);

for ii=1:length(t)
    HAB = [1 0 0;
           0 1 0;
        -dY(ii) dX(ii) 1];    % translation operator for straight members
    
    cA = cos(alphaA(ii));
    sA = sin(alphaA(ii));
    cB = cos(alphaB(ii));
    sB = sin(alphaB(ii));
    
    lA = [cA -sA 0;
          sA  cA 0;
          0   0  1];
    
    lB = [cB -sB 0;
          sB  cB 0;
          0   0  1];
    
    angle = pi/2;
    
    TAB = [  cos(angle)    sin(angle)    0;    % rotation/translation operator for curved members
        -sin(angle)    cos(angle)    0     % local ref systems in A and B are off by angle
        -r*(1-cos(angle)) r*sin(angle)  1];
    
    disp(['element= ',num2str(ii)])
    disp(['   '])
    
    if ii==1 || ii==6   % for straight components
        
        KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_straight*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_straight*lB';      % KAB
        KG(4:6,4:6) = lB*kB_straight*lB';      % KBB
        
    elseif ii==2 || ii==5     % for all columns
        
        KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA';      % KAA
        KG(4:6,1:3) = -lB*kB_column*HAB'*lA';      % KBA
        KG(1:3,4:6) = -lA*HAB*kB_column*lB';      % KAB
        KG(4:6,4:6) = lB*kB_column*lB';      % KBB
        
    else        % for curved members
        
        KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA';    % KAA
        KG(4:6,1:3) = -lB*kB_curved*TAB'*lA';    % KBA
        KG(1:3,4:6) = -lA*TAB*kB_curved*lB';    % KAB
        KG(4:6,4:6) = lB*kB_curved*lB';    % KBB
        
    end
    
    node1 = t(1,ii);
    node2 = t(2,ii);
    
    AA = node1*3-2 : node1*3;  % AA contains the dofs for the node A of the member
    BB = node2*3-2 : node2*3;  % BB contains the dofs for the node B of the member
    
    % Adding the element stiffness matrix coefficients in the global
    % stiffness matrix of the entire structure
    
    KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
    KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
    KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
    KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
    
end

% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];

KK = KKG([1,3:6,10:12,16:19,21],[1,3:6,10:12,16:19,21]);

displ = KK\Force;

%Displacements of unconstrained nodes
ux = displ([4,7,10,13,16,19,22,25]);
uy = displ([11,14,17,20]);
phi = displ([6,9,12,15,18,21,24]);

% Table of nodal displacements for the structure
disp(['     Node       '   '  ux [in]       '   'uy [in]         '   'phi [rad]'])
disp([[2;4;6;7] , ux ,  [0;uy;0] ,    phi])

% Displacement components for all nodes
displALL_x = [ux(1),ux(2),0,ux(3),0,ux(4),ux(5)];
displALL_y = [0,uy(1),0,uy(2),0,uy(3),0];

n_coord = [000 120 120 195 270 270 390;...
    035 075 000 085 000 075 035].*12;

xcoord = n_coord(1,:);
ycoord = n_coord(2,:);

scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on

I am new to matlab. I have to analyse the structure in the image below by solving global stiffness matrix.I have code this.But since i am new i dont know where i am making mistakes and I am not able to run the results.

Answer 1

[Gautam]

Hello Vruddhi,

I understand that you are facing errors while executing the MATLAB code. There are a few modifications you can make to the code to mitigate the errors

• Line 114 initializes the matrix ‘KG’ as

KG = zeros(2*ndofn,2*ndofn);

This creates a matrix of type ‘double’. Line 148 attempts to replace a portion of the array KG with a symbolic matrix

KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA'; 

Since KG is of type ‘double’, MATLAB will throw an error saying that it is unable to convert expression containing symbolic variables into double array. To resolve the error, you can consider converting the matrix KG into a symbolic array for direct replacement

KG = sym(zeros(2*ndofn,2*ndofn)); 

• Similar is the case with lines 111 and 178. Converting the matrix of type ‘double’ to type ‘sym’ would resolve the error

KKG = sym(zeros(nnodes*ndofn,nnodes*ndofn)); 

• Lines 179, 180 and 181 use the arrays AA and BB for indexing the matrix KKG which is of size 27x27. However, the array BB contains the values ‘28’, ‘29’ and ‘30’. This would give the ‘Index out of bounds’ error. To resolve this, make sure that the matrix KKG is declared with correct dimensions and the indexing is correct
• Similarly, in line 185, the matrix “Force” matrix has size 29x1 but the code tries to access the 30th element.

Force([1:3,28:29],:) = []; 

This will again give ‘index out of bound’ error. Make sure the matrix is of appropriate size and that you’re accessing the correct elements.

Following these steps, the code would give results without any errors.

Thank you

Gautam Murthy