Jusqu'à présent, l'analyse du squelette d'avion était faite avec un programme créé en Fortran, mais dans le cas du squelette d'avion, le degré de liberté total n'est pas si grand et il semble que cela puisse être apprécié avec Python, j'ai donc décidé de travailler dessus.
J'avais chez moi un manuel "Dynamique structurelle de l'ordinateur personnel - Visualisation du comportement dynamique - Kenji Miyazawa, Keigaku Shuppan, juin 1984", alors j'ai réécrit ce programme de base en Python. Bien sûr, les équations linéaires simultanées sont résolues en utilisant numpy.
Je pense que la partie calcul du programme peut être écrite de manière assez compacte, mais la partie entrée / sortie semble être compliquée. Est-ce inévitable?
La fonction d'analyse n'est que basique. Il applique une force externe pour obtenir le déplacement de chaque nœud et la force de section transversale de l'élément, et ne supporte pas de déplacement forcé autre que zéro, changement de température, calcul automatique de la force d'inertie, etc.
Tout d'abord, à partir de ce programme de base, je vais essayer d'améliorer les fonctions si j'ai le temps dans le futur. Aussi, j'essaierai Python-matplotlib pour le programme de création de diagramme de force de section transversale.
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 | Nombre de nœuds |
| nele | nombre d'éléments |
| nsec | Nombre de caractéristiques transversales |
| npfix | Nombre de nœuds de contrainte |
| nlod | Nombre de nœuds de chargement |
| A, I, E | Section transversale, moment secondaire transversal, coefficient élastique |
| node_1, node_2, isec | Node 1, Node 2, Numéro de caractéristique de section |
| x, y | coordonnée x, coordonnée y |
| lp, fix_x, fix_y, fix_r | Numéro de nœud, x · y. Présence ou absence de contrainte dans le sens de rotation (0: libre, 1: contrainte complète) |
| lp, fp_x, fp_y, fp_r | Numéro de nœud, xy / charge rotationnelle |
Comme les données d'entrée sont petites, le fichier de données est créé dans le script d'exécution.
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()
Si les équations simultanées sont grandes, je pense que vous devriez utiliser la méthode sur le site suivant.
http://org-technology.com/posts/solving-linear-equations-QR.html
c'est tout
Recommended Posts