Bisher wurde die Flugzeugskelettanalyse mit einem in Fortran erstellten Programm durchgeführt, aber im Fall des Flugzeugskeletts ist der Gesamtfreiheitsgrad nicht so groß und es scheint, dass er mit Python genossen werden kann. Deshalb habe ich beschlossen, daran zu arbeiten.
Ich hatte zu Hause ein Lehrbuch "Strukturdynamik von Personalcomputern - Visualisierung des dynamischen Verhaltens - Kenji Miyazawa, Keigaku Shuppan, Juni 1984", also schrieb ich dieses Basisprogramm in Python um. Natürlich werden simultane lineare Gleichungen mit numpy gelöst.
Ich denke, dass der Berechnungsteil des Programms ziemlich kompakt geschrieben werden kann, aber der Eingabe- / Ausgabeteil scheint chaotisch zu sein. Ist es unvermeidlich?
Die Analysefunktion ist nur grundlegend. Eine externe Kraft wird angewendet, um die Verschiebung jedes Knotens und die Querschnittskraft des Elements zu erhalten, und sie unterstützt keine erzwungene Verschiebung außer Null, Temperaturänderung, automatische Berechnung der Trägheitskraft usw.
Zunächst werde ich von diesem Basisprogramm aus versuchen, die Funktionen zu verbessern, wenn ich in Zukunft Zeit habe. Außerdem werde ich Python-matplotlib für das Programm zur Erstellung von Querschnittskraftdiagrammen ausprobieren.
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 | Anzahl der Knoten |
| nele | Anzahl der Elemente |
| ns | Anzahl der Querschnittseigenschaften |
| npfix | Anzahl der Einschränkungsknoten |
| nlod | Anzahl der Ladeknoten |
| A, I, E | Querschnittsfläche, Querschnittssekundärmoment, Elastizitätskoeffizient |
| Knoten_1, Knoten_2, isec | Knoten 1, Knoten 2, Schnittmerkmalnummer |
| x, y | x-Koordinate, y-Koordinate |
| lp, fix_x, fix_y, fix_r | Knotennummer, x · y. Mit oder ohne Zurückhaltung in Drehrichtung (0: frei, 1: vollständige Rückhaltung) |
| lp, fp_x, fp_y, fp_r | Knotennummer, xy / Rotationslast |
Da die Eingabedaten klein sind, wird die Datendatei im Ausführungsskript erstellt.
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()
Wenn die simultanen Gleichungen groß sind, sollten Sie die Methode an der folgenden Stelle anwenden.
http://org-technology.com/posts/solving-linear-equations-QR.html
das ist alles
Recommended Posts