TRUSS ANALYSIS v4.1 was first created in 14 June 2001 using QBasic and updated
in 09 December 2023 using Phyton by Ernel Solis Filipinas
(ernel_filipinas@yahoo.com).
https://colab.research.google.com/drive/1xNp5Hd96YDruSFZOiwcobsjLhk-DI_Aa?usp=sharing
(For educational
reference only)
-----------------------------------------------------copy
below-----------------------------------------------------
import numpy as np
## INPUT DATA
#From Cañete (2021) Example 3-3 of pagse 98-104.
# Define joint coordinate [x,y] in mm
jdat1 = np.array([
[0,0],
[1000,0],
[500,1000]
])
# Define node/joint type [] where 1-Hinge, 2-Roller x, 3-Roller y, 4-Free Joint
jdat3 = np.array([[1],[2],[4]])
# Define external horizontal load in kN (+ve going right, -ve going left)
jdat4 = np.array([[0],[0],[10]])
# Define external vertical load in kN (+ve going up, -ve going down)
jdat5 = np.array([[0],[0],[0]])
# Define element/member incidences (joint end1, joint end2)
mdat1 = np.array([[1,2],[1,3],[2,3]])
# Define Area in mm2
mdat3 = np.array([[2500],[2500],[2500]])
# Define Modulus of Elasticity (GPa)
mdat4 = np.array([[200],[200],[200]])
#----------------------------------------------------------------------------------
#----------------------------------------------------------------------------------
## CALCULATION FUNCTION
# Step 0---------------------------------------------------------------------------
def title_disclaimer():
print("2D TRUSS ANALYSIS v4.1 first created in 14 June 2001 using QBasic")
print("Updated in 09 December 2023 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,jdat3,jdat4,jdat5,mdat1, mdat3, mdat4):
print("Step 1: Checking number of joint and number of member-----------------------\n")
noj = jdat1.shape[0]
nojdat3 = jdat3.shape[0]
nojdat4 = jdat4.shape[0]
nojdat5 = jdat5.shape[0]
print ("Number of joints =" ,noj)
if noj == nojdat3:
if noj == nojdat4:
if noj == nojdat5:
print("The number joints and joint datum is OK")
else:
print("Number of joints not equal to number of vertical loads")
else:
print("Number of joints not equal to number of horizontal loads ")
else:
print("Number of joints not equal to number of joint type")
nom = mdat1.shape[0]
nonmdat3 = mdat3.shape[0]
nonmdat4 = mdat4.shape[0]
print ("Number of members " ,nom)
if nom == nonmdat3:
if nom == nonmdat4:
print("The number members and member datum is OK\n")
else:
print("Number of members not equal to number of member areas\n")
else:
print("Number of members not equal to number of member modulus\n")
#Step 2---------------------------------------------------------------------------
# Define functions for NDF, TDOF
def DOF(jdat1,jdat3):
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 jdat3[i] == 1:
NDF[i,0] = 0
NDF[i,1] = 0
elif jdat3[i] == 2:
NDF[i,0] = TDOF + 1
NDF[i,1] = 0
TDOF+=1
elif jdat3[i] == 3:
NDF[i,0] = 0
NDF[i,1] = TDOF + 1
TDOF+=1
elif jdat3[i] == 4:
NDF[i,0] = TDOF + 1
NDF[i,1] = TDOF + 2
TDOF+=2
else:
print("Please check the joint data")
print("NDF Array x = ", NDF[i,0],"NDF Array y = ", NDF[i,1])
print ("TDOF = ",TDOF,"\n")
return NDF, TDOF
#Step 3---------------------------------------------------------------------------
# Define functions for Qsys
def Q(jdat1,jdat4,jdat5, NDF,TDOF):
print("Step 3: Compute for Qsys-----------------------\n")
noj = jdat1.shape[0]
Qsys = np.zeros(shape=(TDOF),dtype = float)
for i in range (0,noj):
for j in range (0,2):
if NDF[i,j] > 0:
if j==0:
Qsys[NDF[i,0]-1] = jdat4[i].item ()
print("Qsys",NDF[i,0],"=", Qsys[NDF[i,0]-1],"kN")
if j==1:
Qsys[NDF[i,1]-1] = jdat5[i].item()
print("Qsys",NDF[i,1],"=", Qsys[NDF[i,1]-1],"kN")
print("")
return Qsys
#Step 4---------------------------------------------------------------------------
# Define functions for Mcode
def MCODE(mdat1):
print("Step 4: Compute for Mcode-----------------------\n")
nom = mdat1.shape[0]
mcode = np.zeros(shape=(nom,4),dtype = int)
MDIFF = 0
for i in range(0,nom):
print("Member:",i+1)
print("Member Incidence",mdat1[i,0],mdat1[i,1])
mcode[i,0] = NDF[mdat1[i,0]-1,0]
mcode[i,1] = NDF[mdat1[i,0]-1,1]
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 5---------------------------------------------------------------------------
# Define functions for Ksys
def create_stiffness_matrix(jdat1, mdat1, mdat3, mdat4,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,4))
# 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,"mm")
# 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)
#AEL
AEL = mdat4[i] * mdat3[i] / L
mprop[i,3] = AEL.item()
print("AE/L:",AEL,"kN mm^2")
# Local stiffness matrix (4x4)
K = AEL * np.array([[ l*l , l*m ,-l*l ,-l*m],
[ l*m , m*m ,-l*m ,-m*m],
[-l*l ,-l*m , l*l , l*m],
[-l*m ,-m*m , l*m , m*m]])
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]
#Allocating space for the solution vector
Dsys = np.zeros_like(y, dtype=np.double);
#Here we perform the back-substitution.
#Initializing with the last row.
Dsys[-1] = y[-1] / U[-1, -1]
#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(n-2, -1, -1):
Dsys[i] = (y[i] - np.dot(U[i,i:], Dsys[i:])) / U[i,i]
print("Dsys:",'\n',Dsys,'\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]
f = np.zeros((nom))
R = np.zeros((noj,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("D",j,":",D[j],"mm")
print("AEL =",mprop[i,3])
f[i] = mprop[i,3]*((D[2]-D[0])*mprop[i,1]+(D[3]-D[1])*mprop[i,2])
print("member forces =",f[i],"kN")
Fx1 = -f[i]*mprop[i,1]
Fy1 = -f[i]*mprop[i,2]
Fx2 = f[i]*mprop[i,1]
Fy2 = f[i]*mprop[i,2]
print("Fx1=",Fx1,"kN")
print("Fy1=",Fy1,"kN")
print("Fy2=",Fx2,"kN")
print("Fy2=",Fy2,"kN\n")
R[mdat1[i,0]-1,0] = R[mdat1[i,0]-1,0] + Fx1
R[mdat1[i,0]-1,1] = R[mdat1[i,0]-1,1] + Fy1
R[mdat1[i,1]-1,0] = R[mdat1[i,1]-1,0] + Fx2
R[mdat1[i,1]-1,1] = R[mdat1[i,1]-1,1] + Fy2
for m in range(0,noj):
print ("Joint: ", m+1)
for n in range(0,2):
if NDF[m,n]!=0:
R[m,n] = 0
print("Reaction x",m+1,":",R[m,0],"kN")
print("Reaction y",m+1,":",R[m,1],"kN")
print ("\n")
return f,R
## RUN---------------------------------------------------------------------------
# Call for the title and disclaimer
title_disclaimer()
# Step 1 - Compute for correctness of in terms of numbers
check(jdat1,jdat3,jdat4,jdat5,mdat1, mdat3, mdat4)
# Step 2 - Compute for NDF and TDOF
NDF,TDOF = DOF(jdat1,jdat3)
# Step 3 - Compute for Qsys
Qsys = Q(jdat1,jdat4,jdat5, NDF,TDOF)
# Step 4 - Compute for Mcode
mcode, HBW = MCODE(mdat1)
# Step 5 - Compute for stiffness matrix
Ksys , mprop = create_stiffness_matrix(jdat1, mdat1, mdat3, mdat4,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
f,R = Forces(jdat1,mdat1,mcode,mprop,NDF)
No comments:
Post a Comment