Python Code
Version 3.11.1, Raw code pasted as text
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import math
#
#x1 = np.arange(9.0).reshape((3, 3))
#x2 = np.arange(9.0).reshape((3, 3))
#print(x1)
#print('\n')
#print(x2)
#print('\n')
#print(x1 * x2)
#Simulation Parameters
members = 3
nodes = 4
dof = 3
conmats = np.array([[4, 5, 6, 7, 8, 9],[1, 2, 3, 10, 11, 12], [7, 8, 9, 10, 11, 12]])
#Properties
E = 30 * 10**6 #psi
I = 650 #in^4
A = 10 #in^2
L = 10 #in
E1 = E
E2 = E
E3 = 2*E
A1 = A
A2 = 0.5*A
A3 = A
I1 = I
I2 = 0.5*I
I3 = I
L1 = L
L2 = 2*L
L3 = 1.25*L
Ek = np.array([E1,E2,E3])
Ik = np.array([I1,I2,I3])
Ak = np.array([A1,A2,A3])
Lk = np.array([L1,L2,L3])
thetak = np.array([0.5*np.pi, 2.0/3.0*np.pi,0])
lamdak = np.cos(thetak)
muk = np.sin(thetak)
rk = np.array([0.0,0.0,0.0])
i=0
while( i < thetak.size):
rk[i] = Ak[i] / Ik[i] * Lk[i]**2
i+=1
i=0
##print("Rk")
##print(rk)
##print("Ik")
##print(Ik)
##print("Ak")
##print(Ak)
##print("Lk")
##print(Lk)
print("Lamdak")
print(lamdak)
print("Muk")
print(muk)
#Stiffness Matrix for Elastic Frame Elements
k = 0
sizeG = int(dof*nodes)
sizeK = dof * 2 + 1
Kg = np.zeros((sizeG, sizeG), dtype = float)
Kk = np.empty((sizeK, sizeK), dtype = float)
Kmem = np.array([[Kk],[Kk],[Kk]])
#populate global stiffness matrix
while (k < members):
multK = (Ek[k] * Ik[k] / (Lk[k]**3))
Kk = [[0, conmats[k][0], conmats[k][1], conmats[k][2], conmats[k][3], conmats[k][4], conmats[k][5] ],
#Row 1
[conmats[k][0], (rk[k] * (lamdak[k]**2) + 12 * muk[k]**2) * multK, ((rk[k] - 12) * lamdak[k] * muk[k]) * multK,
(0-6 * Lk[k] * muk[k]) * multK, (0-rk[k]*lamdak[k]**2 - 12 * muk[k]**2) * multK, (0 - (rk[k] - 12) *lamdak[k]*muk[k]) * multK, (0-6*Lk[k]*muk[k]) * multK],
#Row 2
[ conmats[k][1], ((rk[k]-12)*lamdak[k]*muk[k]) * multK, (12 * lamdak[k]**2 + rk[k] * muk[k]**2)*multK,
(6 * Lk[k] * lamdak[k]) * multK, (0 - (rk[k] - 12) * lamdak[k] * muk[k]) * multK, (0 - (12 * lamdak[k]**2 + rk[k] * muk[k]**2))*multK, (6 * Lk[k] * lamdak[k]) * multK ],
#Row 3
[ conmats[k][2], (0-6 * Lk[k] * muk[k]) *multK, (6 * Lk[k] * lamdak[k]) * multK, (4 * Lk[k]**2) * multK, (6 * Lk[k]*muk[k]) * multK,
(0-6 * Lk[k] * lamdak[k]) * multK, (2 * Lk[k]**2) * multK],
#Row 4
[ conmats[k][3], (0 - rk[k] * lamdak[k]**2 - 12 * muk[k]**2) * multK, (0 - (rk[k] - 12) * lamdak[k] * muk[k]) * multK, (6 * Lk[k] * muk[k])*multK,
(rk[k] * lamdak[k]**2 + 12 * muk[k]**2) *multK, ((rk[k] - 12) *lamdak[k] *muk[k]) * multK, (6 * Lk[k] * muk[k])*multK],
#Row 5
[conmats[k][4], (0- (rk[k] - 12) * lamdak[k]* muk[k]) * multK, (0 - 12 * lamdak[k]**2 - rk[k] * muk[k]**2) *multK, (0 - 6 *Lk[k] *lamdak[k]) *multK,
((rk[k] - 12) * lamdak[k] * muk[k]) * multK, (12 * lamdak[k]**2 + rk[k] * muk[k]**2) *multK, (0 - 6 * Lk[k] * lamdak[k]) * multK],
#Row 6
[conmats[k][5], (0 - 6 * Lk[k] * muk[k]) * multK, (6 * Lk[k] * lamdak[k]) * multK, 2 * Lk[k]**2 * multK, 6* Lk[k]* muk[k] * multK,
(0 - 6 * Lk[k] * lamdak[k]) * multK, (4 * Lk[k]**2) * multK]
]
## ii = 0
## while(ii<7):
## print(Kk[ii])
## ii+=1
## ii = 0
Kmem[k] = Kk
jj = 1
while (jj < sizeK):
row = int(Kk[jj][0])
ll = 1
while(ll < sizeK):
column = int(Kk[0][ll])
Kg[row-1][column-1] += Kk[jj][ll]
ll += 1
ll = 1
jj+=1
jj = 1
k +=1
ii = 0
print("--------------GLOBAL-----------------")
while(ii<12):
print(Kg[ii])
ii+=1
ii = 0
#Assemble Load Vectors
#Reference Frames and in-frame vector loads
frame1 = np.array([[lamdak[0], muk[0],0,0,0,0],
[0-muk[0], lamdak[0],0,0,0,0],
[0,0,1,0,0,0],
[0,0,0,lamdak[0],muk[0],0],
[0,0,0,0-muk[0],lamdak[0],0],
[0,0,0,0,0,1]])
q1 = np.array([[0],[0],[0],[0],[0],[0]])
frame2 = np.array([[lamdak[1], muk[1],0,0,0,0],
[0-muk[1], lamdak[1],0,0,0,0],
[0,0,1,0,0,0],
[0,0,0,lamdak[1],muk[1],0],
[0,0,0,0-muk[1],lamdak[1],0],
[0,0,0,0,0,1]])
q2 = np.array([[0],[0],[0],[0],[0],[0]])
frame3 = np.array([[lamdak[2], muk[2],0,0,0,0],
[0-muk[2], lamdak[2],0,0,0,0],
[0,0,1,0,0,0],
[0,0,0,lamdak[2],muk[2],0],
[0,0,0,0-muk[2],lamdak[2],0],
[0,0,0,0,0,1]])
q3 = np.array([[3000],[-3125],[-6510.42],[0],[-3125],[6510.42]])
Q1 = np.zeros((6,2), float)
Q2 = np.zeros((6,2), float)
Q3 = np.zeros((6,2), float)
ii=0
while(ii<2*dof):
Q1[ii][0]= conmats[0][ii]
Q1[ii][1]= np.matmul(frame1, q1)[ii]
Q2[ii][0]= conmats[1][ii]
Q2[ii][1]= np.matmul(frame2, q2)[ii]
Q3[ii][0]= conmats[2][ii]
Q3[ii][1]= np.matmul(frame3, q3)[ii]
ii+=1
ii=0
#Assemble Global Force Vector
Qg = np.zeros((dof*nodes,1), float)
ii = 0
while(ii < 2*dof):
Qg[int(Q1[ii][0]) - 1] += Q1[ii][1]
Qg[int(Q2[ii][0]) - 1] += Q2[ii][1]
Qg[int(Q3[ii][0]) - 1] += Q3[ii][1]
ii+=1
ii=0
print(Q1)
print(Q2)
print(Q3)
print(Qg)
#segment global stiffness matrix and force vectors into restrained and unrestrained segments
Kff = np.zeros((7,7), float)
Kss = np.zeros((5,5), float)
Ksf = np.zeros((5,7), float)
Kfs = np.zeros((7,5), float) # not needed
for x in range(0,5):
for y in range(0,5):
Kss[x][y] = Kg[x][y]
for x in range(0,7):
for y in range(0,7):
Kff[x][y] = Kg[x+5][y+5]
for x in range(0,5):
for y in range(0,7):
Ksf[x][y] = Kg[x][y+5]
Kfs = np.transpose(Ksf)
Qf = np.zeros((7,1),float)
x=0
for x in range(0,7):
Qf[x] = Qg[x+5]
print("Qf")
print(Qf)
##ii=0
##print("Kss")
##while(ii<5):
## print(Kss[ii])
## ii+=1
##ii = 0
##print("Kff")
##while(ii<7):
## print(Kff[ii])
## ii+=1
##ii = 0
##print("Ksf")
##
##while(ii<5):
## print(Ksf[ii])
## ii+=1
##ii = 0
##print("Kfs")
##while(ii<7):
## print(Kfs[ii])
## ii+=1
##ii = 0
#Compute Unknown Forces and Displacements
#displacements
qf = np.matmul(np.linalg.inv(Kff),Qf)
print("qf")
print(qf)
Qs = np.matmul(Ksf, qf)
print("Qs")
print(Qs)
#while (i < members):
#
#
# i += 1
qe1 = np.zeros((6,1), float)
qe1[0][0]= 0
qe1[1][0]= 0
qe1[2][0]= qf[0][0]
qe1[3][0]= qf[1][0]
qe1[4][0]= qf[2][0]
qe1[5][0]= qf[3][0]
qe2 = np.zeros((6,1), float)
qe2[0][0]=0
qe2[1][0]=0
qe2[2][0]=0
qe2[3][0]= qf[4][0]
qe2[4][0]= qf[5][0]
qe2[5][0]= qf[6][0]
qe3 = np.zeros((6,1), float)
qe3[0][0]= qf[1][0]
qe3[1][0]= qf[2][0]
qe3[2][0]= qf[3][0]
qe3[3][0]= qf[4][0]
qe3[4][0]= qf[5][0]
qe3[5][0]= qf[6][0]
#Calculate in-element forces
k=0
lmultK = float(Ek[k] * Ik[k] / Lk[k]**3)
inLoadCoordTrans = np.array([
[rk[k]*lamdak[k], rk[k]*muk[k], 0, 0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0],
[0-12*muk[k], 12*lamdak[k], 6*Lk[k], 12*muk[k], 0-12*lamdak[k], 6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 4*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 2*Lk[k]**2],
[0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0, rk[k]*lamdak[k], rk[k]*muk[k], 0],
[12*muk[k], 0-12*lamdak[k], 0-6*Lk[k], 0-12*muk[k], 12*lamdak[k], 0-6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 2*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 4*Lk[k]**2]
])
inLoad1 = np.zeros((6,1), float)
k=0
lmultK = float(Ek[k]*Ik[k]/Lk[k]**3)
inLoadCoordTrans = np.array([
[rk[k]*lamdak[k], rk[k]*muk[k], 0, 0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0],
[0-12*muk[k], 12*lamdak[k], 6*Lk[k], 12*muk[k], 0-12*lamdak[k], 6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 4*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 2*Lk[k]**2],
[0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0, rk[k]*lamdak[k], rk[k]*muk[k], 0],
[12*muk[k], 0-12*lamdak[k], 0-6*Lk[k], 0-12*muk[k], 12*lamdak[k], 0-6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 2*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 4*Lk[k]**2]
])
inLoad1 = lmultK * np.matmul(inLoadCoordTrans, qe1)
inLoad2 = np.zeros((6,1), float)
k=1
lmultK = float(Ek[k]*Ik[k]/Lk[k]**3)
inLoadCoordTrans = np.array([
[rk[k]*lamdak[k], rk[k]*muk[k], 0, 0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0],
[0-12*muk[k], 12*lamdak[k], 6*Lk[k], 12*muk[k], 0-12*lamdak[k], 6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 4*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 2*Lk[k]**2],
[0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0, rk[k]*lamdak[k], rk[k]*muk[k], 0],
[12*muk[k], 0-12*lamdak[k], 0-6*Lk[k], 0-12*muk[k], 12*lamdak[k], 0-6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 2*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 4*Lk[k]**2]
])
inLoad2 = lmultK * np.matmul(inLoadCoordTrans, qe2)
inLoad3 = np.zeros((6,1), float)
k=2
lmultK = float(Ek[k]*Ik[k]/Lk[k]**3)
inLoadCoordTrans = np.array([
[rk[k]*lamdak[k], rk[k]*muk[k], 0, 0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0],
[0-12*muk[k], 12*lamdak[k], 6*Lk[k], 12*muk[k], 0-12*lamdak[k], 6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 4*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 2*Lk[k]**2],
[0-rk[k]*lamdak[k], 0-rk[k]*muk[k], 0, rk[k]*lamdak[k], rk[k]*muk[k], 0],
[12*muk[k], 0-12*lamdak[k], 0-6*Lk[k], 0-12*muk[k], 12*lamdak[k], 0-6*Lk[k]],
[0-6*Lk[k]*muk[k], 6*Lk[k]*lamdak[k], 2*Lk[k]**2, 6*Lk[k]*muk[k], 0-6*Lk[k]*lamdak[k], 4*Lk[k]**2]
])
inLoad3 = lmultK * np.matmul(inLoadCoordTrans, qe3)
print("Loads Member 1:")
print(inLoad1)
print("Loads Member 2:")
print(inLoad2)
print("Loads Member 3:")
print(inLoad3)
x = np.linspace(0, Lk[k])
y = np.linspace(-math.sqrt(Ak[k]), math.sqrt(Ak[k]))
K = np.empty((dof*2, dof*2), dtype = float)
m = 0
n = 0
##k=2
##multK = (Lk[k]/L)
##x[0]=0
##y[0]=0
##strainmat3 = np.array([0-lamdak[k]/(10*multK)-y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2))*muk[k],
## y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2)*lamdak[k])-muk[k]/(10*multK),
## y[n]*(6*x[m]/(10*multK**3)-4/(10*multK**2))*10*multK,
## lamdak[k]/(10*multK)-y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*muk[k],
## y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*lamdak[k]+muk[k]/(10*multK),
## y[n]*(6*x[m]/(10*multK**3)-(3/(10*multK**2)))*10*multK])
##strain1 = np.matmul(strainmat3, qe3)
##print('*****')
##print(x)
##print(y)
##print(strain1)
m=0
n=0
strainlog1 = np.zeros((50,50), float)
while(m < 50 ):
while(n < 50):
k=0
multK = (Lk[k]/L)
strainmat1 = np.array([0-lamdak[k]/(10*multK)-y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2))*muk[k],
y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2)*lamdak[k])-muk[k]/(10*multK),
y[n]*(6*x[m]/(10*multK**3)-4/(10*multK**2))*10*multK,
lamdak[k]/(10*multK)-y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*muk[k],
y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*lamdak[k]+muk[k]/(10*multK),
y[n]*(6*x[m]/(10*multK**3)-(3/(10*multK**2)))*10*multK])
strain1 = np.matmul(strainmat1, qe1)
strainlog1[m][n] = float(strain1)
n+=1
n=0
m+=1
print(strainlog1)
m=0
n=0
strainlog2 = np.zeros((50,50), float)
x = np.linspace(0, Lk[k])
y = np.linspace(-math.sqrt(Ak[k]), math.sqrt(Ak[k]))
while(m < 50 ):
while(n < 50):
k=1
multK = (Lk[k]/L)
strainmat2 = np.array([0-lamdak[k]/(10*multK)-y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2))*muk[k],
y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2)*lamdak[k])-muk[k]/(10*multK),
y[n]*(6*x[m]/(10*multK**3)-4/(10*multK**2))*10*multK,
lamdak[k]/(10*multK)-y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*muk[k],
y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*lamdak[k]+muk[k]/(10*multK),
y[n]*(6*x[m]/(10*multK**3)-(3/(10*multK**2)))*10*multK])
strain2 = np.matmul(strainmat2, qe2)
strainlog2[m][n] = float(strain2)
n+=1
n=0
m+=1
print(strainlog2)
m=0
n=0
strainlog3 = np.zeros((50,50), float)
x = np.linspace(0, Lk[k])
y = np.linspace(-math.sqrt(Ak[k]), math.sqrt(Ak[k]))
while(m < 50 ):
while(n < 50):
k=2
multK = (Lk[k]/L)
strainmat3 = np.array([0-lamdak[k]/(10*multK)-y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2))*muk[k],
y[n]*(12*x[m]/(10*multK**3)-6/(10*multK**2)*lamdak[k])-muk[k]/(10*multK),
y[n]*(6*x[m]/(10*multK**3)-4/(10*multK**2))*10*multK,
lamdak[k]/(10*multK)-y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*muk[k],
y[n]*(6/(10*multK**2)-12*x[m]/(10*multK**3))*lamdak[k]+muk[k]/(10*multK),
y[n]*(6*x[m]/(10*multK**3)-(3/(10*multK**2)))*10*multK])
Bsd3 = np.array([[-1/Lk[k], 0, 0, 1/Lk[k], 0, 0],
[0, -y[n] * (6/Lk[k]**2 -12*x[m]/Lk[k]**3), -y[n] * (4/Lk[k] - 6 * x[m] / Lk[k]**2),
0, -y[n] * ( 12 * x[m] / Lk[k] **3 - 6 / Lk[k] **2), -y[n] * (2/Lk[k] - 6 * x[m] / Lk[k]**2)]])
#strain3 = np.matmul(Bsd3, inLoad3)
#print(strain3)
strain3 = np.matmul(strainmat3, qe3)
strainlog3[m][n] = float(strain3)
n+=1
n=0
m+=1
print(strainlog3)
k=0
stresslog1 = Ek[k]*strainlog1
k=1
stresslog2 = Ek[k]*strainlog2
k=2
stresslog3 = Ek[k]*strainlog3
print("elem1")
print(stresslog1.max())
print(stresslog1.min())
print("elem2")
print(stresslog2.max())
print(stresslog2.min())
print("elem3")
print(stresslog3.max())
print(stresslog3.min())
figure1 = plt.matshow(np.transpose(stresslog1))
plt.title('max: ' + str(int(stresslog1.max()))+"psi" + '---min:' + str(int(stresslog1.min()))+ "psi")
figure2 = plt.matshow(np.transpose(stresslog2))
plt.title('max: ' + str(int(stresslog2.max()))+"psi" + '---min:' + str(int(stresslog2.min()))+ "psi")
figure3 = plt.matshow(np.transpose(stresslog3))
plt.title('max: ' + str(int(stresslog3.max()))+"psi" + '---min:' + str(int(stresslog3.min()))+ "psi")
plt.style.use('_mpl-gallery-nogrid')
# make data
X, Y = np.meshgrid(np.linspace(0, 10), np.linspace(-math.sqrt(Ak[k]), math.sqrt(Ak[k])))
Z = 3.773*X*Y-277.96
levels = np.linspace(Z.min(), Z.max(), 7)
# plot
fig, ax = plt.subplots()
ax.contourf(X, Y, Z, levels=levels)
plt.show()