Pythonで平面骨組解析プログラムを作る
これまで平面骨組解析はFortranで作成したプログラムで行っていましたが，平面骨組の場合，総自由度もそれほど大きくなく，Pythonでも楽しめそうなので，取り組んでみることにしました．
ちょうど家に教科書「パソコン構造力学ー力学挙動の可視化ー宮澤健二著 啓学出版 1984年6月」があったので，このBasicプログラムをPythonに書き換えました．もちろん連立一次方程式は，numpyを用いて解いています．
プログラムの計算部分は結構コンパクトに書けていると思いますが，入出力部分がゴチャゴチャした感じです．しょうがないのですかね．
解析機能は，基本的なもののみです．外力を与え各節点変位および要素の断面力を求めるものであり，ゼロ以外の強制変位，温度変化，慣性力の自動計算などには対応していません．
まず，この基本プログラムから，今後時間があれば機能強化を図っていこうと思います．
また断面力図作成プログラムもPython-matplotlibでやってみようと思います．
プログラムと入出力事例
入力データフォーマット
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 | 節点数 |
| nele | 要素数 |
| nsec | 断面特性数 |
| npfix | 拘束節点数 |
| nlod | 載荷節点数 |
| A, I, E | 断面積，断面2次モーメント，弾性係数 |
| node_1, node_2, isec | 節点１，節点２，断面特性番号 |
| x, y | ｘ座標，ｙ座標 |
| lp, fix_x, fix_y, fix_r | 節点番号，ｘ・ｙ．回転方向の拘束の有無 （0: 自由，1：完全拘束） |
| lp, fp_x, fp_y, fp_r | 節点番号，ｘ・ｙ・回転方向荷重 |
実行用スクリプト
入力データは少ないので，実行用スクリプトの中でデータファイルを作成しています．
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()
参考サイト
連立方程式が大規模になる場合は，以下のサイトの方法を用いると良いと思います．
http://org-technology.com/posts/solving-linear-equations-QR.html
以 上