Week 11 Project 2- MATLAB Module to calculate Forces and Displacements for a 3D Frame Structure(nonlinear analysis)

Objective: MATLAB Module to calculate Forces and Displacements for a 3D Frame Structure(nonlinear analysis)

Solution procedure:

1) Develop a Matlab function to accept Node information: Node numbers, Node coordinates, and information on nodes to perform analysis on a 3D planar frame structure.

2) Based on the Matlab module developed in problem 1 process the node information to calculate the degree of freedom for each node assigned

3) Develop a Matlab module to accept information on elements in the structure and their connectivity info and receive the section and material properties as inputs

4) Based on the Matlab module developed in problem 3 processes the element information and based on their connectivity calculate element degrees of freedom for each element

5) Develop a Matlab module to receive boundary conditions and degrees of freedom specified for the structure

6) Develop a Matlab module to receive the load vector information acting on the structure

7) Based on the stiffness method develop a Matlab module to compute the stiffness matrix for
the structure in global coordinates

8) Using the stiffness matrix developed in problem 7 compute the forces and displacements by implementing the concept of nested functions

1) Develop a Matlab function to accept Node information: Node numbers, Node coordinates, and information on nodes to perform analysis on a 3D planar frame structure

function node_xyz =GetNodeCoordinates(~)
    %Taking input
prompt_1 ='Enter the Node coordinates ';
disp('input in [node_x1 node_y1,node_z1;node_x2 node_y2,node_z2;.......;node_xn,node_yn,node_zn]');
node_xyz =input(prompt_1);
end

2) Based on the Matlab module developed in problem 1 process the node information to calculate the degree of freedom for each node assigned

function [node_coord,dof] = GetNodeInfo(~)
    node_xyz = GetNodeCoordinates();
    %To calculate the degree of freedom for each nodes
    %For 3D frame elements

    %for nodes number
    node_no = node_xyz(:,1);
    %calculate number of nodes

    no_of_nodes = size(node_xyz(:,1),1)
    dof =zeros(no_of_nodes,3);

    for i =1:no_of_nodes
        node_coord(i,:) = [node_xyz(i,2),node_xyz(i,3)]
        dof(i,:) = (3*(node_no(i)-1)^2)+(1:3)
    end
end

3) Develop a Matlab module to accept information on elements in the structure and their connectivity info and receive the section and material properties as inputs

function[node_coord,element_no,element_node,length_element,element_geometry,element_dof]=GetElementInfo
[node_coord,dof] =GetNodeInfo();
prompt ='Input Element Connectivity Data';
disp('Enter vector in this format [Element No,node_i,node_j,b,h]');
elements =input(prompt);
element_no =elements(:,1)
element_node =[elements(:,2),elements(:,3)]
length_element =sqrt((node_coord(element_node(:,2),1)-node_coord(element_node(:,1),1)).^2+(node_coord(element_node(:,2),2)-node_coord(element_node(:,1),2)).^2)
element_geometry = [elements(:,4),elements(:,5)]
element_dof = [dof(element_node(:,1),:),dof(element_node(:,2),:)]
disp(length_element)
end

4) Based on the Matlab module developed in problem 3 processes the element information and based on their connectivity calculate element degrees of freedom for each element

5) Develop a Matlab module to receive boundary conditions and degrees of freedom specified for the structure

function[boundary_conditions] = GetBoundary_Conditions()
    prompt = 'Enter Boundary Conditions for the Structure System';
    disp('Enter in the following format [node_no. dof displacement]');
    boundary_conditions = input(prompt);
end

6) Develop a Matlab module to receive the load vector information acting on the structure

function [node_loads]=GetNodalLoads()
    prompt = 'Enter Nodal Loads';
    disp('Enter in this format[node_no dof of applied load value')
    node_loads = input(prompt);
    
end

7) Based on the stiffness method develop a Matlab module to compute the stiffness matrix for
the structure in global coordinates (Hint: Calculate stiffness matrix for each element, Transform the element stiffness
the matrix in their global coordinates and assemble the element stiffness matrix based on the element degrees of freedom for the structure)

%Global Stiffness Matrix
for k=1:no_elements
    ele_nodes=element_node(k,:);
    node_xy=node_coord(ele_nodes,:);
    A = element_geometry(k,1)*element_geometry(k,2);
    I=element_geometry(k,1)*element_geometry(k,2).^3/12;
    length=length_element(k);
    k_element=Frame_Stiffness_Matrix1(node_xy,A,I,length);
    eldofs=element_dof(k,:);
    global_stiffness_matrix(eldofs,eldofs)=global_stiffness_matrix(eldofs,eldofs)+k_element;

8) Using the stiffness matrix developed in problem 7 compute the forces and displacements by implementing the concept of nested functions

displacement_matrix(free_dof)=global_stiffness_matrix(free_dof,free_dof)(force_matrix(free_dof)-global_stiffness_matrix(free_dof,dof_specified)*displacement_matrix(dof_specified));
force_matrix(dof_specified)=global_stiffness_matrix(dof_specified,:)*displacement_matrix;
disp('Displacement Vector');
disp(displacement_matrix);
disp('Force Vector');
disp(force_matrix);

end

function [global_stiffness_matrix]=ComputeGlobalStiffnessMatrix()
%Get Element Info
[node_coord,element_no,element_node,length_element,element_geometry,element_dof]=GetElementInfo();
%Get Boundary Conditions Info
boundary_conditions=GetBoundary_Conditions();
%Get Nodal loads Info
node_loads = GetNodalLoads();


%Calculating Number Of Elements And Nodes
no_elements = size(element_no,1);
no_nodes = size(node_coord,1);

%Intializing Matrix with Zeros
global_stiffness_matrix= zeros(3*no_nodes,3*no_nodes);
displacement_matrix= zeros(3*no_nodes,1);
force_matrix = zeros(3*no_nodes,1);

%Boundary Conditions
all_dof= 1:3*no_nodes;
dof_specified =[]
for i=1:no_nodes
    n_dof=3*(boundary_conditions(i,1)-1)+boundary_conditions(i,2);
    dof_specified=[dof_specified n_dof];
    displacement_matrix(n_dof)=boundary_conditions(i,3);
end
free_dof=all_dof;
free_dof(dof_specified)=[]

%Compute Nodal Loads
for j =1:size(node_loads,1)
    force_matrix(3*(node_loads(j,1)-1)+node_loads(j,2))=node_loads(j,3);
end

%Global Stiffness Matrix
for k=1:no_elements
    ele_nodes=element_node(k,:);
    node_xyz=node_coord(ele_nodes,:);
    A = element_geometry(k,1)*element_geometry(k,2);
    I=element_geometry(k,1)*element_geometry(k,2).^3/12;
    length=length_element(k);
    k_element=Frame_Stiffness_Matrix1(node_xy,A,I,length);
    eldofs=element_dof(k,:);
    global_stiffness_matrix(eldofs,eldofs)=global_stiffness_matrix(eldofs,eldofs)+k_element;
    
end
displacement_matrix(free_dof)=global_stiffness_matrix(free_dof,free_dof)(force_matrix(free_dof)-global_stiffness_matrix(free_dof,dof_specified)*displacement_matrix(dof_specified));
force_matrix(dof_specified)=global_stiffness_matrix(dof_specified,:)*displacement_matrix;
disp('Displacement Vector');
disp(displacement_matrix);
disp('Force Vector');
disp(force_matrix);

end

Explanation of the above code:- 

function [global_stiffness_matrix]=ComputeGlobalStiffnessMatrix()
%Get Element Info
[node_coord,element_no,element_node,length_element,element_geometry,element_dof]=GetElementInfo(); - taking data from the functions
%Get Boundary Conditions Info
boundary_conditions=GetBoundary_Conditions();
%Get Nodal loads Info
node_loads = GetNodalLoads();

%Calculating Number Of Elements And Nodes
no_elements = size(element_no,1);
no_nodes = size(node_coord,1);

%Intializing Matrix with Zeros
global_stiffness_matrix= zeros(3*no_nodes,3*no_nodes);
displacement_matrix= zeros(3*no_nodes,1);
force_matrix = zeros(3*no_nodes,1);

%Boundary Conditions
all_dof= 1:3*no_nodes;
dof_specified =[]
for i=1:no_nodes
n_dof=3*(boundary_conditions(i,1)-1)+boundary_conditions(i,2);
dof_specified=[dof_specified n_dof];
displacement_matrix(n_dof)=boundary_conditions(i,3);
end
free_dof=all_dof;
free_dof(dof_specified)=[]

%Compute Nodal Loads
for j =1:size(node_loads,1)
force_matrix(3*(node_loads(j,1)-1)+node_loads(j,2))=node_loads(j,3);
end

%Global Stiffness Matrix
for k=1:no_elements - for loop to iteration
ele_nodes=element_node(k,:); - assigning the elements nodes
node_xy=node_coord(ele_nodes,:); - assigning the coordinates for node x and node y
A = element_geometry(k,1)*element_geometry(k,2);
I=element_geometry(k,1)*element_geometry(k,2).^3/12;
length=length_element(k);
k_element=Frame_Stiffness_Matrix1(node_xy,A,I,length);
eldofs=element_dof(k,:);
global_stiffness_matrix(eldofs,eldofs)=global_stiffness_matrix(eldofs,eldofs)+k_element;

end
displacement_matrix(free_dof)=global_stiffness_matrix(free_dof,free_dof)(force_matrix(free_dof)-global_stiffness_matrix(free_dof,dof_specified)*displacement_matrix(dof_specified));
force_matrix(dof_specified)=global_stiffness_matrix(dof_specified,:)*displacement_matrix;
disp('Displacement Vector');
disp(displacement_matrix);- display the dispacement amtrix
disp('Force Vector');- display the force vector name
disp(force_matrix);- display the force matrix

end

9) Solve the following problem using the Matlab module developed in 1-8 to compute forces and displacements. Verify the Matlab answers with the answers of the previous challenge

function [global_stiffness_matrix]=ComputeGlobalStiffnessMatrix()
%Get Element Info
[node_coord,element_no,element_node,length_element,element_geometry,element_dof]=GetElementInfo();
%Get Boundary Conditions Info
boundary_conditions=GetBoundary_Conditions();
%Get Nodal loads Info
node_loads = GetNodalLoads();


%Calculating Number Of Elements And Nodes
no_elements = size(element_no,1);
no_nodes = size(node_coord,1);

%Intializing Matrix with Zeros
global_stiffness_matrix= zeros(3*no_nodes,3*no_nodes);
displacement_matrix= zeros(3*no_nodes,1);
force_matrix = zeros(3*no_nodes,1);

%Boundary Conditions
all_dof= 1:3*no_nodes;
dof_specified =[]
for i=1:no_nodes
    n_dof=3*(boundary_conditions(i,1)-1)+boundary_conditions(i,2);
    dof_specified=[dof_specified n_dof];
    displacement_matrix(n_dof)=boundary_conditions(i,3);
end
free_dof=all_dof;
free_dof(dof_specified)=[]

%Compute Nodal Loads
for j =1:size(node_loads,1)
    force_matrix(3*(node_loads(j,1)-1)+node_loads(j,2))=node_loads(j,3);
end

%Global Stiffness Matrix
for k=1:no_elements
    ele_nodes=element_node(k,:);
    node_xyz=node_coord(ele_nodes,:);
    A = element_geometry(k,1)*element_geometry(k,2);
    I=element_geometry(k,1)*element_geometry(k,2).^3/12;
    length=length_element(k);
    k_element=Frame_Stiffness_Matrix1(node_xy,A,I,length);
    eldofs=element_dof(k,:);
    global_stiffness_matrix(eldofs,eldofs)=global_stiffness_matrix(eldofs,eldofs)+k_element;
    
end
displacement_matrix(free_dof)=global_stiffness_matrix(free_dof,free_dof)(force_matrix(free_dof)-global_stiffness_matrix(free_dof,dof_specified)*displacement_matrix(dof_specified));
force_matrix(dof_specified)=global_stiffness_matrix(dof_specified,:)*displacement_matrix;
disp('Displacement Vector');
disp(displacement_matrix);
disp('Force Vector');
disp(force_matrix);

end

The output of the above code:-