Until now, plane skeleton analysis was performed with a program created in Fortran, but in the case of plane skeleton, the total degree of freedom is not so large and it seems that it can be enjoyed in Python, so I decided to work on it.
Just at home, I had a textbook "PC Structural Mechanics-Visualization of Mechanics Behavior-Kenji Miyazawa, Keigaku Shuppan, June 1984", so I rewrote this Basic program to Python. Of course, simultaneous linear equations are solved using numpy.
I think that the calculation part of the program can be written quite compactly, but the input / output part seems to be messy. Is it unavoidable?
The analysis function is only basic. External force is applied to obtain the displacement of each node and the cross-sectional force of the element, and it does not correspond to forced displacement other than zero, temperature change, automatic calculation of inertial force, etc.
First of all, from this basic program, I will try to enhance the functions if I have time in the future. I will also try Python-matplotlib for the section force diagram creation program.
npoin nele nsec npfix nlod
A I E
...(1 to nele)...
node_1 node_2 isec
...(1 to nele)...
x y
...(1 to npoin)...
lp fix_x fix_y fix_r
...(1 to npfix)...
lp fp_x fp_y fp_r
...(1 to nlod)...
| npoin | Number of nodes |
| nele | number of elements |
| nsec | Number of cross-sectional characteristics |
| npfix | Number of constraint nodes |
| nlod | Number of loading nodes |
| A, I, E | Cross-sectional area, moment of inertia of area, elastic modulus |
| node_1, node_2, isec | Node 1, Node 2, Sectional characteristic number |
| x, y | x coordinate, y coordinate |
| lp, fix_x, fix_y, fix_r | Node number, x · y. With or without restraint in the direction of rotation (0: free, 1: complete restraint) |
| lp, fp_x, fp_y, fp_r | Node number, xy-rotational load |
Since the input data is small, the data file is created in the execution script.
a_py_frame.txt
cat << EOT > inp.txt
6 5 2 3 1
1.0 1.0 1.0
1.0 1.0 1.0
1 2 1
3 4 1
5 6 1
2 4 2
4 6 2
0 0
0 3.5
6 0
6 3.5
12 0
12 3.5
1 1 1 1
3 1 1 1
5 1 1 1
2 4 0 0
EOT
python3 py_frame.py inp.txt out.txt
npoin nele nsec npfix nlod
6 5 2 3 1
sec A I E
1 1.0000000e+00 1.0000000e+00 1.0000000e+00
2 1.0000000e+00 1.0000000e+00 1.0000000e+00
node x y fx fy fr kox koy kor
1 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 1 1 1
2 0.0000000e+00 3.5000000e+00 4.0000000e+00 0.0000000e+00 0.0000000e+00 0 0 0
3 6.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 1 1 1
4 6.0000000e+00 3.5000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0 0 0
5 1.2000000e+01 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 1 1 1
6 1.2000000e+01 3.5000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0 0 0
elem i j sec
1 1 2 1
2 3 4 1
3 5 6 1
4 2 4 2
5 4 6 2
node dis-x dis-y dis-r
1 0.0000000e+00 0.0000000e+00 0.0000000e+00
2 1.6079284e+01 2.3039125e+00 -4.5858390e+00
3 0.0000000e+00 0.0000000e+00 0.0000000e+00
4 5.6044784e+00 -1.4855500e+00 -6.2687943e-01
5 0.0000000e+00 0.0000000e+00 0.0000000e+00
6 2.6990174e+00 -8.1836247e-01 -5.5363182e-01
elem N_i S_i M_i N_j S_j M_j
1 -6.5826071e-01 2.2541991e+00 5.2550881e+00 6.5826071e-01 -2.2541991e+00 2.6346087e+00
2 4.2444286e-01 1.2615574e+00 2.3868338e+00 -4.2444286e-01 -1.2615574e+00 2.0286170e+00
3 2.3381785e-01 4.8424351e-01 1.0056067e+00 -2.3381785e-01 -4.8424351e-01 6.8924562e-01
4 1.7458009e+00 -6.5826071e-01 -2.6346087e+00 -1.7458009e+00 6.5826071e-01 -1.3149555e+00
5 4.8424351e-01 -2.3381785e-01 -7.1366148e-01 -4.8424351e-01 2.3381785e-01 -6.8924562e-01
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
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
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
disg = np.linalg.solve(gk, fp)
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()
If the simultaneous equations are large, I think you should use the method at the following site.
http://org-technology.com/posts/solving-linear-equations-QR.html
that's all
Recommended Posts