LoginSignup
2
4

More than 5 years have passed since last update.

@damyarou

Pythonで平面骨組解析 (2) 修正プログラム

Posted at

概要

このプログラムでは，剛性マトリックスの大きさを変えずに境界条件を導入している（下式参照）
ここで，既知変位はゼロのみ処理可能としているため，例えば，変位 $\delta_i=0$ であれば，$k_{ii}=1$ とし，それ以外の $i$ 行および $i$ 列の要素をゼロとしている．
ここで，もし $f_i$ がゼロでない場合は，$f_i=0$ とする必要がある．
前回のプログラムでは，この処理が織り込まれていなかったので，今回修正した．

元の剛性方程式
$$
\begin{equation}
\begin{bmatrix}
k_{11} & \ldots & k_{1i} & \ldots & k_{1n} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
k_{1i} & \ldots & k_{ii} & \ldots & k_{1n} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
k_{n1} & \ldots & k_{ni} & \ldots & k_{nn}
\end{bmatrix}
\begin{Bmatrix}
\delta_1 \\
\vdots \\
\delta_i \\
\vdots \\
\delta_n
\end{Bmatrix}
=\begin{Bmatrix}
f_1 \\
\vdots \\
f_i \\
\vdots \\
f_n
\end{Bmatrix}
\end{equation}
$$

境界条件導入後の剛性方程式 ($\delta_i=0, \quad f_i=0$)
$$
\begin{equation}
\begin{bmatrix}
k_{11} & \ldots & 0 & \ldots & k_{1n} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
0 & \ldots & 1 & \ldots & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
k_{n1} & \ldots & 0 & \ldots & k_{nn}
\end{bmatrix}
\begin{Bmatrix}
\delta_1 \\
\vdots \\
\delta_i \\
\vdots \\
\delta_n
\end{Bmatrix}
=\begin{Bmatrix}
f_1 \\
\vdots \\
f_i \\
\vdots \\
f_n
\end{Bmatrix}
\end{equation}
$$

修正プログラム

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
fw   =np.zeros(3*npoin,dtype=np.float64)   # Work vector for 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
fw=fp
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
            fw[iz]=0.0

disg = np.linalg.solve(gk, fw)

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()

以　上

2
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
4