2D BEAM ANALYSIS v4 first created as FRAME ANALYSIS on 24 July 2009 using Excel VBA and updated on 10 July 2024 using Phyton by Ernel Solis Filipinas (ernel_filipinas@yahoo.com).
import math
import numpy as np
import matplotlib.pyplot as plt
import sys
## INPUT DATA
# Define joint coordinate [x,y] in m
jdat1 = np.array([
[0,0],
[1,0],
[2,0],
[3,0],
[4,0],
[5,0],
[6,0],
[7,0],
[8,0],
[9,0],
[10,0]
])
# Define node/joint type [] where 1-Hinge, 2-Roller x, 3-Roller y, 4-Free, 5-Fixed Joints
jdat2 = np.array([[5],[4],[4],[4],[4],[4],[4],[4],[4],[4],[5]])
# Define Moment Joint Load in kNm
# External Force (+ve counter counter clockwise, -ve clockwise)
# Fixed End force (inverse of External Force)
jdat5 = np.array([[0],[0],[0],[0],[0],[0],[0],[0],[0],[0],[0]])
# Define Vertical Joint Load in kN
# External Force (+ve going up, -ve going down)
# Fixed End force (inverse of External Force)
jdat6 = np.array([[0],[0],[0],[0],[0],[-4],[0],[0],[0],[0],[0]])
# Define element/member incidences (joint end1, joint end2)
mdat1 = np.array([[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11]])
# Define Modulus of Elasticity (kPa)
mdat2 = np.array([[200e6],[200e6],[200e6],[200e6],[200e6],[200e6],[200e6],[200e6],[200e6],[200e6]])
# Define Moment of Inertia (m4)
mdat3 = np.array([[5.208e-7],[5.208e-7],[5.208e-7],[5.208e-7],[5.208e-7],[5.208e-7],[5.208e-7],[5.208e-7],[5.208e-7],[5.208e-7]])
# Define Fixed End Moment for beam kNm # Make sure it is as per member incidence
#(+ve counter counter clockwise, -ve clockwise)
mdat4 = np.array([[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0]])
# Define Fixed End Shear for beam kNm
#(+ve going up, -ve going down)
mdat5 = np.array([[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0],[0,0]])
#----------------------------------------------------------------------------------
#----------------------------------------------------------------------------------
#FEM Link
#https://blog.payrollschedule.net/#google_vignette
#https://engineering.purdue.edu/~ce474/Docs/Fixed%20End%20Moments.pdf
#https://mathalino.com/reviewer/strength-materials/fixed-end-moments-fully-restrained-beam#google_vignette
# Scale
dscale = 5000 #Deformation Scale
ascale = 200 #Arrow Scale
#----------------------------------------------------------------------------------
#----------------------------------------------------------------------------------
## CALCULATION FUNCTION
# Step 0---------------------------------------------------------------------------
def title_disclaimer():
print("2D BEAM ANALYSIS v4 first created as FRAME ANALYSIS in 24 January 2009 using Excel VBA")
print("10 July 2024 using Phyton")
print("by Ernel Solis Filipinas. \ne: ernel_filipinas@yahoo.com\n\n")
print("Disclaimer:\n")
print("No Liability is accepted by its software authors for any direct, ")
print("indirect, consequential or incidental loss or damage arising out of ")
print("the software use or any mistakes and negligence in developing this software.")
print("The organisation or person using the software bears all risks and ")
print("responsibility for the quality and performance of the software. \n")
print("Reference:\n")
print("Espejo, Isagani S. (2000). The Matrix Method - Computer-Aided Structural Analysis")
print("Cañete, Alberto (2021). Principles of Matrix Structural Analysis\n")
#Step 1---------------------------------------------------------------------------
# Define functions checking number of joint and number of member
def check(jdat1,jdat2,jdat5,jdat6,mdat1, mdat2, mdat3, mdat4, mdat5):
print("Step 1: Checking number of joint and number of member-----------------------\n")
noj = jdat1.shape[0]
nojdat2 = jdat2.shape[0]
nojdat5 = jdat5.shape[0]
nojdat6 = jdat6.shape[0]
print ("Number of joints =" ,noj)
if noj == nojdat2 and noj == nojdat5 and noj == nojdat6:
print("The number joints and joint datum is OK")
else:
sys.exit("Number of joints not equal for jdat data")
nom = mdat1.shape[0]
nonmdat2 = mdat2.shape[0]
nonmdat3 = mdat3.shape[0]
nonmdat4 = mdat4.shape[0]
nonmdat5 = mdat5.shape[0]
print ("Number of members " ,nom)
if nom == nonmdat2 and nom==nonmdat3 and nom==nonmdat4 and nom==nonmdat5:
print("The number members and member datum is OK\n")
else:
sys.exit("Number of members not equal to nomdat data")
#Step 2---------------------------------------------------------------------------
# Define functions for NDF, TDOF
def DOF(jdat1,jdat2):
print("Step 2: Compute for NDF & TDOF-----------------------\n")
TDOF = 0
noj = jdat1.shape[0]
NDF = np.zeros(shape=(noj, 2),dtype=int)
for i in range (0,noj):
print("Joint: ", i+1)
if jdat2[i] == 1: #Hinge
NDF[i,0] = 0
NDF[i,1] = TDOF + 1
TDOF+=1
elif jdat2[i] == 2: #Roller x
NDF[i,0] = 0
NDF[i,1] = TDOF + 1
TDOF+=1
elif jdat2[i] == 3: #Roller y
NDF[i,0] = TDOF + 1
NDF[i,1] = TDOF + 2
TDOF+=2
elif jdat2[i] == 4: #Free
NDF[i,0] = TDOF + 1
NDF[i,1] = TDOF + 2
TDOF+=2
elif jdat2[i] == 5: #Fixed
NDF[i,0] = 0
NDF[i,1] = 0
else:
print("Please check the joint data")
print("NDF Array Y = ", NDF[i,0],"NDF Array M = ", NDF[i,1])
print ("TDOF = ",TDOF,"\n")
return NDF, TDOF
#Step 3---------------------------------------------------------------------------
# Define functions for Mcode
def MCODE(mdat1,NDF, TDOF):
print("Step 4: Compute for Mcode-----------------------\n")
nom = mdat1.shape[0]
mcode = np.zeros(shape=(nom,4),dtype = int) # for beam only
MDIFF = 0
for i in range(0,nom):
print("Member:",i+1)
print("Member Incidence",mdat1[i,0],mdat1[i,1])
#mcode[i,] = 0 # For beam only
mcode[i,0] = NDF[mdat1[i,0]-1,0]
mcode[i,1] = NDF[mdat1[i,0]-1,1]
#mcode[i,] = 0 # For beam only
mcode[i,2] = NDF[mdat1[i,1]-1,0]
mcode[i,3] = NDF[mdat1[i,1]-1,1]
print("mcode",mcode[i,0],mcode[i,1],mcode[i,2],mcode[i,3])
MAX = 0
MIN = TDOF
for j in range (0,4):
if mcode[i,j] == 0:
continue
elif mcode[i,j] > MAX:
MAX = mcode[i,j]
if mcode[i,j] < MIN:
MIN = mcode[i,j]
print("MAX",MAX)
print("MIN",MIN)
if (MAX-MIN) >= MDIFF:
MDIFF = MAX - MIN
print("MDIFF",MDIFF,"\n")
HBW = MDIFF+1
print("HBW",HBW,"\n")
return mcode, HBW
#Step 4---------------------------------------------------------------------------
# Define functions for Qsys
def Q(jdat1,jdat5,jdat6, NDF,TDOF,MCODE):
print("Step 3: Compute for Qsys-----------------------\n")
noj = jdat1.shape[0]
Qsys = np.zeros(shape=(TDOF),dtype = float)
nom = mdat1.shape[0]
F = np.zeros(shape=(noj),dtype = float)
for i in range (0,noj):
for j in range (0,2):
if NDF[i,j] > 0:
if j==0:
p = NDF[i,0]-1
r = jdat6[i]
Qsys[p] = r.item()
print("Qsys",NDF[i,0],"=", Qsys[p],"kN")
if j==1:
q = NDF[i,1]-1
s = jdat5[i]
Qsys[q] = s.item()
print("Qsys",NDF[i,1],"=", Qsys[q],"kNm")
return Qsys
#Step 5---------------------------------------------------------------------------
# Define functions for Ksys
def create_stiffness_matrix(jdat1, mdat1, mdat2, mdat3,TDOF,mcode):
print("Step 5: Compute for Ksys-----------------------\n")
# Initialize matrix
nom = mdat1.shape[0]
K = np.zeros((4,4))
Ksys = np.zeros((TDOF,TDOF))
mprop = np.zeros((nom,6))
# Extract node coordinates=
#x=jdat1[:,0] # sliced array with only x node coordinates
#y=jdat1[:,1] # sliced array with only y node coordinates
for i in range(0,nom):
print ("Member: ", i+1)
x1 = jdat1[mdat1[i,0]-1,0]
y1 = jdat1[mdat1[i,0]-1,1]
x2 = jdat1[mdat1[i,1]-1,0]
y2 = jdat1[mdat1[i,1]-1,1]
print ("joint",mdat1[i,0],":","x1:",x1,"y1:",y1)
print ("joint",mdat1[i,1],":","x2:",x2,"y2:",y2)
# Element length
L = np.sqrt((x2- x1)**2 + (y2-y1)**2)
mprop[i,0] = L
print("L:",L,"m")
# Cosine and sine of element angle
l = (x2 - x1) / L
m = (y2 - y1) / L
mprop[i,1] = l
mprop[i,2] = m
print("l:",l)
print("m:",m)
#EI
EIL = mdat2[i] * mdat3[i]/L
mprop[i,3] = EIL.item()
mprop[i,4] = (6*EIL.item())/(L**1)
mprop[i,5] = (12*EIL.item())/(L**2)
print("EI/L:",EIL,"kNm/rad")
k11 = 12/(L**2)
k12 = 6/(L**1)
# Local stiffness matrix (4x4)
K = EIL * np.block([[ k11 , k12 ,-k11 , k12],
[ k12 , 4 ,-k12 , 2],
[-k11 ,-k12 , k11 ,-k12],
[ k12 , 2 ,-k12 , 4]
])
print(K)
print("mcode",mcode[i,0],mcode[i,1],mcode[i,2],mcode[i,3])
for j in range(0,4):
for k in range(0,4):
if mcode[i,j] != 0 and mcode[i,k] != 0:
x = mcode[i,j]
y = mcode[i,k]
print("x,y:",x,y)
print("K[",j+1,k+1,"] =",K[j,k])
Ksys[x-1,y-1] = Ksys[x-1,y-1] + K[j,k]
print("Ksys : \n",Ksys,"\n")
# Return global stiffness matrix
return Ksys , mprop
#Step 6---------------------------------------------------------------------------
# Define functions for the Gobal Displacement
def cholesky(Ksys):
print("Step 6: Compute for Gloabal Displacement-----------------------\n")
print("Ksys : \n",Ksys)
n = len(Ksys)
L = np.zeros((n, n))
for i in range(n):
for j in range(i+1):
if i == j:
L[i, j] = np.sqrt(Ksys[i, i] - np.sum(L[i, :]**2))
else:
L[i, j] = (Ksys[i, j] - np.sum(L[i, :] * L[j, :])) / L[j, j]
print("L:",'\n',L)
return L
def forward_substitution(L, Qsys):
print("Qsys:",'\n',Qsys)
#Get number of rows
n = L.shape[0]
#Allocating space for the solution vector
y = np.zeros_like(Qsys, dtype=np.double);
#Here we perform the forward-substitution.
#Initializing with the first row.
y[0] = Qsys[0] / L[0, 0]
#Looping over rows in reverse (from the bottom up),
#starting with the second to last row, because the
#last row solve was completed in the last step.
for i in range(1, n):
y[i] = (Qsys[i] - np.dot(L[i,:i], y[:i])) / L[i,i]
print("y:",'\n',y)
return y
def back_substitution(U, y):
print("U:",'\n',U)
# Number of rows
n = U.shape[0]
# Allocate space for the solution vector
Dsys = np.zeros_like(y, dtype=np.double)
# Initialize with the last row
Dsys[-1] = y[-1] / U[-1, -1]
# Loop over rows in reverse order
for i in range(n - 2, -1, -1):
Dsys[i] = (y[i] - np.dot(U[i, i + 1:], Dsys[i + 1:])) / U[i, i]
print("Dsys:",'\n',Dsys,'\n')
noj = jdat1.shape[0]
for m in range(0,noj):
print ("Joint: ", m+1)
for n in range(0,2):
if NDF[m,n]!=0:
if n==0:
print("Displacement z",m+1,":",Dsys[NDF[m,n]-1],"m")
if n==1:
print("Rotation y",m+1,":",Dsys[NDF[m,n]-1],"rad\n")
if NDF[m,n]==0:
if n==0:
print("Displacement z",m+1,":",0,"m")
if n==1:
print("roation y",m+1,":",0,"rad\n")
return Dsys
# Step 7---------------------------------------------------------------------------
# Define functions for Member Forces & Reactions
def Forces(jdat1,mdat1,mcode,mprop,NDF):
print("Step 7: Compute for Member Forces & Reactions-----------------------\n")
noj = jdat1.shape[0]
nom = mdat1.shape[0]
MDOF = np.zeros((nom,2))
FEM = np.zeros((nom,2))
M = np.zeros((nom,2))
RDOF = np.zeros((nom,2))
FER = np.zeros((nom,2))
R = np.zeros((nom,2))
D = np.zeros((4))
for i in range (0,nom):
print ("Member: ", i+1)
for j in range (0,4):
if mcode[i,j] == 0:
D[j]=0
else:
D[j] = Dsys[mcode[i,j]-1]
print("D1",D[0],"m")
print("D2",D[1],"rad")
print("D3",D[2],"m")
print("D4",D[3],"rad")
print("EI/L =",mprop[i,3])
MDOF[i,0] = mprop[i,3]*(4*D[1]+2*D[3])+mprop[i,4]*(D[0]-D[2])
MDOF[i,1] = mprop[i,3]*(2*D[1]+4*D[3])+mprop[i,4]*(D[0]-D[2])
print("member Moment force end 1 due to DOF = ", MDOF[i,0],"kNm")
print("member Moment force end 2 due to DOF = ", MDOF[i,1],"kNm")
FEM[i,0] = mdat4[i,0]
FEM[i,1] = mdat4[i,1]
print("member Moment force end 1 due to FEM = ", FEM[i,0],"kNm")
print("member Moment force end 2 due to FEM = ", FEM[i,1],"kNm")
M[i,0] = MDOF[i,0]+FEM[i,0]
M[i,1] = MDOF[i,1]+FEM[i,1]
print("member Moment force end 1 Final = ", M[i,0],"kNm")
print("member Moment force end 2 Final = ", M[i,1],"kNm\n")
RDOF[i,0] = mprop[i,4]*(D[1]+D[3])+mprop[i,5]*(D[0]-D[2])
RDOF[i,1] = -mprop[i,4]*(D[1]+D[3])+mprop[i,5]*(D[2]-D[0])
print("member Vertical force end 1 due to DOF = ", RDOF[i,0],"kN")
print("member Vertical force end 2 due to DOF = ", RDOF[i,1],"kN")
FER[i,0] = mdat5[i,0]
FER[i,1] = mdat5[i,1]
print("member Vertica force end 1 due to FEM = ", FER[i,0],"km")
print("member Vertica force end 2 due to FEM = ", FER[i,1],"km")
R[i,0] = RDOF[i,0]+FER[i,0]
R[i,1] = RDOF[i,1]+FER[i,1]
print("member Vertica force end 1 Final = ", R[i,0],"kN")
print("member Vertica force end 2 Final = ", R[i,1],"kN\n")
return M, R
# Step 8---------------------------------------------------------------------------
# Define functions for Graph
def TGraph(jdat1,jdat2,mdat1,mcode):
noj = jdat1.shape[0]
nom = mdat1.shape[0]
dist = np.zeros((4))
Qsys_x = np.zeros(shape=(TDOF),dtype = float)
Qsys_y = np.zeros(shape=(TDOF),dtype = float)
for i in range(0,nom):
x1 = jdat1[mdat1[i,0]-1,0]
y1 = jdat1[mdat1[i,0]-1,1]
x2 = jdat1[mdat1[i,1]-1,0]
y2 = jdat1[mdat1[i,1]-1,1]
# x axis values
xp =np.array([x1,x2])
# corresponding y axis values
yp =np.array([y1,y2])
# plotting the points
plt.plot(xp,yp, color='green', linestyle='solid', linewidth = 2, marker = 'none', ms = 4)
for j in range (0,4):
if mcode[i,j] == 0:
dist[j]=0
else:
dist[j] = Dsys[mcode[i,j]-1]
dx1 = x1 + dist[0]*dscale
dy1 = y1 + dist[1]*dscale
dx2 = x2 + dist[2]*dscale
dy2 = y2 + dist[3]*dscale
# x axis values for displacement
dxp =np.array([dx1,dx2])
# corresponding y axis values jfor displacement
dyp =np.array([dy1,dy2])
# plotting the points
plt.plot(xp,yp, color='green', linestyle='solid', linewidth = 2, marker = ',', ms = 4)
#Adding support
for m in range (0,noj):
Sx1 = jdat1[m,0]
Sx2 = jdat1[m,1]
if jdat2[m] == 1:
plt.plot(Sx1,Sx2, color='Blue', marker = '^', ms = 10)
elif jdat2[m] == 2:
plt.plot(Sx1,Sx2, color='Blue', marker = 'o', ms = 10)
elif jdat2[m] == 3:
plt.plot(Sx1,Sx2, color='Blue', marker = 'o', ms = 10)
elif jdat2[m] == 4:
plt.plot(Sx1,Sx2, color='Blue', marker = ',', ms = 10)
elif jdat2[m] == 5:
plt.plot(Sx1,Sx2, color='Blue', marker = 's', ms = 10)
plt.annotate(m+1, (jdat1[m,0]+ascale/10, jdat1[m,1]+ascale/10))
# naming the x axis
plt.xlabel('x - axis')
# naming the y axis
plt.ylabel('y - axis')
# giving a title to my graph
plt.title('2D Beam Analysis with Graph \n by Ernel S. Filipinas in 2024')
# function to show the plot
plt.show()
## RUN---------------------------------------------------------------------------
# Call for the title and disclaimer
title_disclaimer()
# Step 1 - Compute for correctness of in terms of numbers
check(jdat1,jdat2,jdat5,jdat6,mdat1, mdat2, mdat3, mdat4, mdat5)
# Step 2 - Compute for NDF and TDOF
NDF,TDOF = DOF(jdat1,jdat2)
# Step 3 - Compute for Mcode
mcode, HBW = MCODE(mdat1,NDF, TDOF)
# Step 4 - Compute for Qsys
Qsys = Q(jdat1,jdat5,jdat6, NDF,TDOF,MCODE)
# Step 5 - Compute for stiffness matrix
Ksys , mprop = create_stiffness_matrix(jdat1, mdat1, mdat2, mdat3,TDOF,mcode)
# Step 6 - Compute for Gloabal Displacement
L = cholesky(Ksys)
y = forward_substitution(L, Qsys)
U = np.transpose(L)
Dsys = back_substitution(U, y)
# Step 7: Compute for Member Forces and Reaction
M,R = Forces(jdat1,mdat1,mcode,mprop,NDF)
