In diesem Programm wird die Randbedingung eingeführt, ohne die Größe der Steifigkeitsmatrix zu ändern (siehe die folgende Formel). Da hier nur Null für die bekannte Verschiebung verarbeitet werden kann, beispielsweise wenn die Verschiebung $ \ delta_i = 0 $ ist, setzen Sie $ k_ {ii} = 1 $ und setzen Sie die anderen $ i $ Zeilen und $ i $ Spalten. Das Element ist Null. Wenn hier $ f_i $ nicht Null ist, muss $ f_i = 0 $ gesetzt werden. Dieser Prozess wurde nicht in das vorherige Programm aufgenommen, daher wurde er diesmal korrigiert.
Ursprüngliche Steifheitsgleichung
Starrheitsgleichung nach Einführung der Randbedingung ($ \ delta_i = 0, \ quad f_i = 0
\vdots & \vdots & \vdots & \vdots & \vdots \
0 & \ldots & 1 & \ldots & 0 \
\vdots & \vdots & \vdots & \vdots & \vdots \
k_{n1} & \ldots & 0 & \ldots & k_{nn}
\end{bmatrix}
\begin{Bmatrix}
\delta_1 \
\vdots \
\delta_i \
\vdots \
\delta_n
\end{Bmatrix}
=\begin{Bmatrix}
f_1 \
\vdots \
f_i \
\vdots \
f_n
\end{Bmatrix}
\end{equation}
$$
py_frame.py
import numpy as np
import sys
def STIFF(i,node,x,ae):
ek=np.zeros([6,6],dtype=np.float64) # local stiffness matrix
tt=np.zeros([6,6],dtype=np.float64) # transformation matrix
mpf=node[0,i]
mps=node[1,i]
x1=x[0,mpf-1]
y1=x[1,mpf-1]
x2=x[0,mps-1]
y2=x[1,mps-1]
xx=x2-x1
yy=y2-y1
el=np.sqrt(xx**2+yy**2)
mb=node[2,i]
aa=ae[0,mb-1]
am=ae[1,mb-1]
ee=ae[2,mb-1]
EA=ee*aa
EI=ee*am
ek[0,0]= EA/el; ek[0,1]= 0.0; ek[0,2]= 0.0; ek[0,3]=-EA/el; ek[0,4]= 0.0; ek[0,5]= 0.0
ek[1,0]= 0.0; ek[1,1]= 12*EI/el**3; ek[1,2]= 6*EI/el**2; ek[1,3]= 0.0; ek[1,4]=-12*EI/el**3; ek[1,5]= 6*EI/el**2
ek[2,0]= 0.0; ek[2,1]= 6*EI/el**2; ek[2,2]= 4*EI/el ; ek[2,3]= 0.0; ek[2,4]= -6*EI/el**2; ek[2,5]= 2*EI/el
ek[3,0]=-EA/el; ek[3,1]= 0.0; ek[3,2]= 0.0; ek[3,3]= EA/el; ek[3,4]= 0.0; ek[3,5]= 0.0
ek[4,0]= 0.0; ek[4,1]=-12*EI/el**3; ek[4,2]=-6*EI/el**2; ek[4,3]= 0.0; ek[4,4]= 12*EI/el**3; ek[4,5]=-6*EI/el**2
ek[5,0]= 0.0; ek[5,1]= 6*EI/el**2; ek[5,2]= 2*EI/el ; ek[5,3]= 0.0; ek[5,4]= -6*EI/el**2; ek[5,5]= 4*EI/el
s=yy/el
c=xx/el
tt[0,0]= c; tt[0,1]= s; tt[0,2]=0.0; tt[0,3]=0.0; tt[0,4]=0.0; tt[0,5]=0.0
tt[1,0]= -s; tt[1,1]= c; tt[1,2]=0.0; tt[1,3]=0.0; tt[1,4]=0.0; tt[1,5]=0.0
tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]=1.0; tt[2,3]=0.0; tt[2,4]=0.0; tt[2,5]=0.0
tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]=0.0; tt[3,3]= c; tt[3,4]= s; tt[3,5]=0.0
tt[4,0]=0.0; tt[4,1]=0.0; tt[4,2]=0.0; tt[4,3]= -s; tt[4,4]= c; tt[4,5]=0.0
tt[5,0]=0.0; tt[5,1]=0.0; tt[5,2]=0.0; tt[5,3]=0.0; tt[5,4]=0.0; tt[5,5]=1.0
return ek,tt,mpf,mps
# Main routine
args = sys.argv
fnameR=args[1] # input data file
fnameW=args[2] # output data file
f=open(fnameR,'r')
text=f.readline()
text=text.strip()
text=text.split()
npoin=int(text[0]) # Number of nodes
nele =int(text[1]) # Number of elements
nsec =int(text[2]) # Number of sections
npfix=int(text[3]) # Number of restricted nodes
nlod =int(text[4]) # Number of loaded nodes
n=3*npoin
x =np.zeros([2,npoin],dtype=np.float64) # Coordinates of nodes
ae =np.zeros([3,nsec],dtype=np.float64) # Section characteristics
node =np.zeros([3,nele],dtype=np.int) # Node-element relationship
fp =np.zeros(3*npoin,dtype=np.float64) # External force vector
fw =np.zeros(3*npoin,dtype=np.float64) # Work vector for External force vector
mpfix=np.zeros([3,npoin],dtype=np.int) # Ristrict conditions
ir =np.zeros(6,dtype=np.int) # Work vector for matrix assembly
gk =np.zeros([n,n],dtype=np.float64) # Global stiffness matrix
fsec =np.zeros([6,nele],dtype=np.float64) # Section force vector
work =np.zeros(6,dtype=np.float64) # work vector for section force calculation
# section characteristics
for i in range(0,nsec):
text=f.readline()
text=text.strip()
text=text.split()
ae[0,i]=float(text[0]) #section area
ae[1,i]=float(text[1]) #moment of inertia
ae[2,i]=float(text[2]) #elastic modulus
# element-node
for i in range(0,nele):
text=f.readline()
text=text.strip()
text=text.split()
node[0,i]=int(text[0]) #node_1
node[1,i]=int(text[1]) #node_2
node[2,i]=int(text[2]) #section characteristic number
# node coordinates
for i in range(0,npoin):
text=f.readline()
text=text.strip()
text=text.split()
x[0,i]=float(text[0]) # x-coordinate
x[1,i]=float(text[1]) # y-coordinate
# boundary conditions (0:free, 1:restricted)
for i in range(0,npfix):
text=f.readline()
text=text.strip()
text=text.split()
lp=int(text[0]) #fixed node
mpfix[0,lp-1]=int(text[1]) #fixed in x-direction
mpfix[1,lp-1]=int(text[2]) #fixed in y-direction
mpfix[2,lp-1]=int(text[3]) #fixed in rotation
# load
for i in range(0,nlod):
text=f.readline()
text=text.strip()
text=text.split()
lp=int(text[0]) #loaded node
fp[3*lp-3]=float(text[1]) #load in x-direction
fp[3*lp-2]=float(text[2]) #load in y-direction
fp[3*lp-1]=float(text[3]) #load in rotation
f.close()
for ne in range(0,nele):
ek,tt,mpf,mps=STIFF(ne,node,x,ae)
ck=np.dot(np.dot(tt.T,ek),tt)
ir[5]=3*mps-1
ir[4]=ir[5]-1
ir[3]=ir[4]-1
ir[2]=3*mpf-1
ir[1]=ir[2]-1
ir[0]=ir[1]-1
for i in range(0,6):
for j in range(0,6):
it=ir[i]
jt=ir[j]
gk[it,jt]=gk[it,jt]+ck[i,j]
# boundary condition
fw=fp
for i in range(0,npoin):
for j in range(0,3):
if mpfix[j,i]==1:
iz=i*3+j
for k in range(0,n):
gk[iz,k]=0.0
gk[k,iz]=0.0
gk[iz,iz]=1.0
fw[iz]=0.0
disg = np.linalg.solve(gk, fw)
for ne in range(0,nele):
ek,tt,mpf,mps=STIFF(ne,node,x,ae)
work[0]=disg[3*mpf-3]; work[1]=disg[3*mpf-2]; work[2]=disg[3*mpf-1]
work[3]=disg[3*mps-3]; work[4]=disg[3*mps-2]; work[5]=disg[3*mps-1]
fsec[:,ne]=np.dot(ek,np.dot(tt,work))
fout=open(fnameW,'w')
# print out of input data
print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('sec','A','I','E'),file=fout)
for i in range(0,nsec):
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(i+1,ae[0,i],ae[1,i],ae[2,i]),file=fout)
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s} {8:>5s}'
.format('node','x','y','fx','fy','fr','kox','koy','kor'),file=fout)
for i in range(0,npoin):
lp=i+1
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d} {8:5d}'
.format(lp,x[0,i],x[1,i],fp[3*lp-3],fp[3*lp-2],fp[3*lp-1],mpfix[0,i],mpfix[1,i],mpfix[2,i]),file=fout)
print('{0:>5s} {1:>5s} {2:>5s} {3:>5s}'.format('elem','i','j','sec'),file=fout)
for ne in range(0,nele):
print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
# displacement
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('node','dis-x','dis-y','dis-r'),file=fout)
for i in range(0,npoin):
lp=i+1
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(lp,disg[3*lp-3],disg[3*lp-2],disg[3*lp-1]),file=fout)
# section force
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
.format('elem','N_i','S_i','M_i','N_j','S_j','M_j'),file=fout)
for ne in range(0,nele):
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
.format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
fout.close()
das ist alles
Recommended Posts