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

以　上