R Codes For Planar Structures

Explorations of the Finite Element Method (FEM) with the R Programming Language

This page contains codes, written in the R programming language, of the finite element analyses (FEA) of some simple problems in Engineering Mechanics.

Simplicity along with code transparency is the driving force for all of the codes revealed here. The writing of these codes in the R programming language doesn’t in any way suggests that it offers the best application environment for the FEM coding. Exquisitely written texts are available on the subject of FEA, so I have omitted the gory details of the solutions to the problems presented here.

However, most of the problems will be taken from “Logan, D. L. (2007). A First Course in the Finite Element Method, Thomson”.

1. The Linear Spring Element with 1DOF/Node

The three functions below simplify the solution of problems (Examples 1.1 - 1.3) relating to a linear spring element with 1DOF at each node.

Comments:

1. Generates the matrix for a bar element whose axis coincides with the global x-axis. Argument=Deg of freedom and stiffness
2. Generates the expanded matrix, filled with zeros #Argument=Total DOF of the assemblage (TDOF)
3. Generates the expanded matrix for each element, filled with actual values #Argument=TDOF of the assemblage, an unexpanded element matrix and local node numbers of an element.

##Example1.1, (pp 42-44 of Logan , an assemblage with a homogenous boundary condition)

K1=Element_Matrix(2,1000);K2=Element_Matrix(2,2000);K3=Element_Matrix(2,3000);                                                                       #[1]

K1Exp=Expanded_Element_Matrix(4,K1,1,3);K2Exp=Expanded_Element_Matrix(4,K2,3,4);

                                                                 K3Exp=Expanded_Element_Matri2(4,K3,4,2)                         #[2]     

KTotal=K1Exp+K2Exp+K3Exp                                        #[3]

ReducedKtotal=KTotal[3:4,3:4]                                                                                                                        #[4]

ReducedForce=matrix(c(0,5000),ncol=1,byrow=F)                                                                                                        #[5]

 (solve(ReducedKtotal,ReducedForce))                                                                                                                    #[6]

Comments:

1. Form the element matrix
2. Form the expanded matrix for each element
3. Form the global matrix from the assembly of the expanded element matrices
4. Extract the reduced stiffness by specifying the location of the known forces in the form nrow X ncom
5. Specify the columns of the known forces
6. Obtain the unknown displacements

#######################################End of Example 1.1############

##Example1.2, (pp 45-47, an assemblage with a non-homogenous boundary condition)

Objective: To determine the nodal displacements and nodal forces for the assemblage below.

K1=Element_Matrix(2,200);K2=Element_Matrix(2,200);K3=Element_Matrix(2,200); K4=Element_Matrix(2,200);                                                         #[1]

K1Exp=Expanded_Element_Matrix(5,K1,1,2);K2Exp=Expanded_Element_Matrix(5,K2,2,3);K3Exp=Expanded_Element_Matrix(5,K3,3,4); K4Exp=Expanded_Element_Matrix(5,K4,4,5);

                                                                                  #[2]

KTotal=K1Exp+K2Exp+K3Exp+K4Exp                                                    #[3]

ReducedKtotal=KTotal[2:4,2:4])                                                    #[4]

ReducedForce=matrix(c(0,0,4),ncol=1,byrow=F)                                      #[5]

(solve(ReducedKtotal,ReducedForce))                                               #[6]

Comments:

1. Form the element matrix
2. Form the expanded matrix for each element
3. Form the global matrix from the assembly of the expanded element matrices
4. Extract the reduced stiffness by specifying the location of the known forces in the form nrow X ncol
5. Specify the column of the known nodal forces
6. Solve for the unknown nodal displacements

##Example1.3, (Ex 3.6 on pp 133 of Logan, an assemblage with parallel members)

Objective: To determine the nodal displacements and global nodal forces using the direct stiffness method for the assemblage below.

k1=2*(30*10^6)/50;k2=2*(10*10^6)/30;k3=k2 K1=Element_Matrix(2,k1);K2=Element_Matrix(2,k2);K3=Element_Matrix(2,k3); K1Exp=Expanded_Element_Matrix(4,K1,1,2);K2Exp=Expanded_Element_Matrix(4,K2,2,3);K3Exp=Expanded_Element_Matrix(4,K3,2,4); KTotal=K1Exp+K2Exp+K3Exp) ReducedKtotal=KTotal[2,2]) ReducedForce=matrix(c(8000),ncol=1,byrow=F) disp_of_node_2=(solve(ReducedKtotal,ReducedForce));disp_of_node_2 globalForcevector=KTotal%*%matrix(c(0,disp_of_node_2,0,0),nrow=4,byrow=T);

Comments:

The stiffness value (k) for each element is found by using the relation AE/L.

Running the above code, one obtains disp_of_node_2 as 0.003157895 in.

The global nodal forces (in lb) (the last line) are obtained as:

[1,] -3789.474 [2,] 8000.000 [3,] -2105.263 [4,] -2105.263

2. The Linear Spring Element with 2DOF/Node

The functions below simplify the solution of problems (Examples 1.4 - 1.6) relating to a linear bar element with 2DOF at each node.

Element_Matrix_Inc=function(stiff,ele_DOF,theta){m=cos(theta*pi/180);n=sin(theta*pi/180);stiff*matrix(c(m^2,m*n,-m^2,-m*n,m*n,n^2,-m*n,-n^2,-m^2,-m*n,m^2,m*n,-m*n,-n^2,m*n,n^2),nrow=ele_DOF,byrow=T)};

Expanded_Element_Matrix_Inc=function(TDOF,eMatrix,i,j){r1=(i-1)+i;r2=(i-1)+(i+1);r3=(j-2)+(j+1);r4=(j-2)+(j+2) bigMatrix=matrix(vector(l=TDOF*TDOF),nrow=TDOF,byrow=T);bigMatrix[c(r1,r2,r3,r4),c(r1,r2,r3,r4)]=eMatrix; return (bigMatrix)}

Axial_Stress_Inc=function(YoungMod,elem_len,theta,disp){m=cos(theta*pi/180);n=sin(theta*pi/180);disp_vec=matrix(disp,nrow=length(disp),byrow=T);(YoungMod/elem_len)*(matrix(c(-m,-n,m,n),nrow=1,byrow=T)%*%disp_vec)}

Comments:

[1] The function named Element_Matrix_Inc generates the local stiffness matrix. It requires theta, which of course is the angle of the local coordinate of the element. This matrix is coded based on Eq. 3.4.23 on page 80 of Logan (2007).

[2] The function named Expanded_Element_Matrix_Inc yields the expanded matrix for each element.

[3] The function named Axial_Stress_Inc calculates the axial stress for a given Young's modulus, an element's length, the angle of a local coordinate and the nodal displacements.

##Example1.4, (page 84 of Logan, a plane truss with three members)

Objective: To determine the nodal displacements and the stresses in each member.

Note: Element1=1-2; Element 2=1-3; Element 3= 1-4.

k1=30*(10^6)*2/(120);k2=30*(10^6)*2/(120*sqrt(2));k3=30*(10^6)*2/(120);

Element_Matrix_Inc1=round(Element_Matrix_Inc(k1,4,90),0);

Element_Matrix_Inc2=round(Element_Matrix_Inc(k2,4,45),0);

Element_Matrix_Inc3=round(Element_Matrix_Inc(k3,4,0),0);

K1Exp=Expanded_Element_Matrix_Inc(8,Element_Matrix_Inc1,1,2);

K2Exp=Expanded_Element_Matrix_Inc(8,Element_Matrix_Inc2,1,3);

K3Exp=Expanded_Element_Matrix_Inc(8,Element_Matrix_Inc3,1,4);

KTotal=K1Exp+K2Exp+K3Exp);ReducedKtotal=KTotal[c(1,2),c(1,2)]);

ReducedForce=matrix(c(0,-10000),ncol=1,byrow=F);

disp_of_nodes_1x1y=(solve(ReducedKtotal,ReducedForce));disp_of_nodes_1x1y

globalForcevector=KTotal%*%matrix(c(disp_of_nodes_1x1y[1],disp_of_nodes_1x1y[2],0,0,0,0,0,0),nrow=8,byrow=T));

#######Nodal displacement vectors for each element

Nodal_disp_elem1=c(disp_of_nodes_1x1y[1],disp_of_nodes_1x1y[2],0,0);Nodal_disp_elem2=c(disp_of_nodes_1x1y[1],disp_of_nodes_1x1y[2],0,0);

Nodal_disp_elem3=c(disp_of_nodes_1x1y[1],disp_of_nodes_1x1y[2],0,0)

####Stress Calculation

AxialStress1=Axial_Stress_Inc(30*10^6,120,90,Nodal_disp_elem1);

AxialStress2=Axial_Stress_Inc(30*10^6,120*sqrt(2),45,Nodal_disp_elem2);

AxialStress3=Axial_Stress_Inc(30*10^6,120,0,Nodal_disp_elem3);

stresses=cbind(c(AxialStress1,AxialStress2,AxialStress3));stresses

Comments:

[1] Using the above code, we obtain the displacements (with the unit in) of node 1, given by the variable

disp_of_nodes_1x1y (x- and y-directions) respectively as:

[1,] 0.00414214 [2,] -0.01585786

[2] The stresses (with the unit psi ) in elements 1, 2 and 3 are also yielded, by the variable named stresses in the above code, respectively, as:

[3] The axial force (in kN) developed in each element due to the 20 kN load at node 2 are obtained as:

[1,] -15972.710 [2,] 12036.301 [3,] 9612.617