Python Code

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()