Python 軸対称応力解析プログラム
はじめに
軸対称モデルとして扱う構造物として、水平トンネルなど水平方向に回転軸を持つ構造物と、調圧水槽など鉛直方向に回転軸を持つ構造物がある。 これまで、四角形要素の2次元応力解析プログラムを軸対称解析用に書き換えたプログラムを使ってきたが、座標系(座標軸の方向)を意識せず回転軸を鉛直にして計算を行ったところ、荷重に対して変位及び応力が本来あるべき方向と逆方向の解が計算されてきた。 このため、回転軸(z軸)が水平の場合と、回転軸が鉛直の場合で、場合分けして解析するようプログラムの見直しを行った。
はじめから下図Case-2の座標系を反時計回りに90度回転させてz軸を上向きとした座標を与えればいいのだが、実際のメッシュ切では、2次元応力解析用のメッシャーで作成した要素座標を使うのが便利なため、今回の対応策を導入することにしたものである。
プログラムの改良
以下に座標系と要素のイメージ図を示す。
Case-1およびCase-2では、回転軸をz軸、半径方向軸をr軸としている。
Reference
Reference は通常の2次元応力解析用の座標である。x軸を水平右向き、y軸を鉛直上向きに設定している。 この場合、通常、要素番号は、節点番号を反時計回りに指定して要素を定義する。
Case-1
Case-1は、z軸(回転軸)を水平右向き、r軸(半径方向)を鉛直上向きに設定している。 この場合も、Referenceの2次元応力解析と同様に、要素番号は、節点番号を反時計回りに指定して要素を定義する。(そうするようにプログラムが組まれている)
Case-2
Case-2は、z軸(回転軸)を鉛直上向き、r軸(半径方向)を水平右向きに設定している。 この場合、Reference の2次元応力解析と同様に反時計回りの節点番号で要素を定義するすると、Case-1と比べてz軸とr軸が入れ替わっているため、要素の節点番号が逆周りで定義されたことになり、剛性マトリックスや物体力ベクトル、温度荷重ベクトルなど要素座標を用いて計算する行列やベクトルが通常と逆符号で計算されることとなる。
対応策
このため、今回整備した軸対称応力解析プログラムでは、変数 nzdir
を導入し、座標の方向に応じて、剛性マトリックス、物体力ベクトル、温度荷重ベクトルの方向を逆転させることにした。
すなわち、上図Case-1では、nzdir=1
、上図Case-2ではnzdir=-1
とし、メインプログラム上で、求まった剛性マトリックス、物体力ベクトル、温度荷重ベクトルの符号を調整するようにしている。
sm=sm*float(nzdir) # 剛性マトリックス tfe=tfe*float(nzdir) # 温度荷重ベクトル bfe=bfe*float(nzdir) # 物体力(慣性力)ベクトル
入力データ・フォーマット
npoin nele nsec npfix nlod nzdir # Basic values for analysis E po alpha gamma gkz # Material properties ..... (1 to nsec) ..... # node1 node2 node3 node4 isec # Element connectivity, material set number ..... (1 to nele) ..... # (counter-clockwise order of node number) z r deltaT # Coordinate, Temperature change of node ..... (1 to npoin) ..... # node koz kor rdisz rdisr # Restricted node and restrict conditions ..... (1 to npfix) ..... # node fz fr # Loaded node and loading conditions ..... (1 to nlod) ..... # (omit data input if nlod=0)
変数 | 説明 |
---|---|
npoin, nele, nsec | 節点数,要素数,材料特性数 |
npfix, nlod | 拘束節点数,載荷節点数 |
nzdir | z軸方向(水平右向き: 1,鉛直上向き: −1) |
E, po, alpha | 弾性係数,ポアソン比,線膨張係数 |
gamma, gkz | 単位体積重量,z方向加速度(gの比) |
node1, node2, node3, node4, isec | 節点1,節点2,節点3,節点4,材料特性番号 |
z, r, deltaT | 節点z座標,節点r座標,節点温度変化 |
node, koz, kor | 拘束節点番号,z及びr方向拘束の有無(拘束: 1,自由: 0) |
rdisz, rdisr | z及びr方向変位(無拘束でも0を入力) |
node, fz, fr | 載荷重節点番号,z方向荷重,r方向荷重 |
(注)回転軸zの方向が水平(nzdir=1)であっても鉛直(nzdir=-1)であっても、座標入力は (z , r) の順に行う。
プログラム全文
################################ # Axisymmetric Stress Analysis ################################ import numpy as np from scipy.sparse.linalg import spsolve from scipy.sparse import csr_matrix import sys import time def inpdata_ax4(fnameR,nod,nfree): 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 nzdir =int(text[5]) # direction of z-axis (=1 or =-1) # array declanation ae =np.zeros([5,nsec],dtype=np.float64) # Section characteristics node =np.zeros([nod+1,nele],dtype=np.int) # Node-element relationship x =np.zeros([2,npoin],dtype=np.float64) # Coordinates of nodes deltaT=np.zeros(npoin,dtype=np.float64) # Temperature change of node mpfix =np.zeros([nfree,npoin],dtype=np.int) # Ristrict conditions rdis =np.zeros([nfree,npoin],dtype=np.float64) # Ristricted displacement fp =np.zeros(nfree*npoin,dtype=np.float64) # External force vector # section characteristics for i in range(0,nsec): text=f.readline() text=text.strip() text=text.split() ae[0,i]=float(text[0]) #E : Elastic modulus ae[1,i]=float(text[1]) #po : Poisson's ratio ae[2,i]=float(text[2]) #alpha: Thermal expansion coefficient ae[3,i]=float(text[3]) #gamma: Unit weight of material ae[4,i]=float(text[4]) #gkz : Acceleration in z-direction # 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]) #node_3 node[3,i]=int(text[3]) #node_4 node[4,i]=int(text[4]) #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]) # z-coordinate x[1,i] =float(text[1]) # r-coordinate deltaT[i]=float(text[2]) # Temperature change of node # boundary conditions (0:free, 1:restricted) if 0<npfix: 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 z-direction mpfix[1,lp-1]=int(text[2]) #fixed in r-direction rdis[0,lp-1]=float(text[3]) #fixed displacement in z-direction rdis[1,lp-1]=float(text[4]) #fixed displacement in r-direction # load if 0<nlod: for i in range(0,nlod): text=f.readline() text=text.strip() text=text.split() lp=int(text[0]) #loaded node fp[2*lp-2]=float(text[1]) #load in z-direction fp[2*lp-1]=float(text[2]) #load in r-direction f.close() return npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp def prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node): fout=open(fnameW,'w') # print out of input data print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}' .format('npoin','nele','nsec','npfix','nlod','nzdir'),file=fout) print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}' .format(npoin,nele,nsec,npfix,nlod,nzdir),file=fout) print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s}' .format('sec','E','po','alpha','gamma','gkz'),file=fout) for i in range(0,nsec): print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e}' .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i]),file=fout) print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s}' .format('node','z','r','fz','fr','deltaT','koz','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}' .format(lp,x[0,i],x[1,i],fp[2*i],fp[2*i+1],deltaT[i],mpfix[0,i],mpfix[1,i]),file=fout) if 0<npfix: print('{0:>5s} {1:>5s} {2:>5s} {3:>15s} {4:>15s}' .format('node','koz','kor','rdis_z','rdis_r'),file=fout) for i in range(0,npoin): if 0<mpfix[0,i]+mpfix[1,i]: lp=i+1 print('{0:5d} {1:5d} {2:5d} {3:15.7e} {4:15.7e}' .format(lp,mpfix[0,i],mpfix[1,i],rdis[0,i],rdis[1,i]),file=fout) print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}' .format('elem','i','j','k','l','sec'),file=fout) for ne in range(0,nele): print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}' .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne],node[4,ne]),file=fout) fout.close() def prout_ax4(fnameW,npoin,nele,disg,strs): fout=open(fnameW,'a') # displacement print('{0:>5s} {1:>15s} {2:>15s}'.format('node','dis-z','dis-r'),file=fout) for i in range(0,npoin): lp=i+1 print('{0:5d} {1:15.7e} {2:15.7e}'.format(lp,disg[2*lp-2],disg[2*lp-1]),file=fout) # section force print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}' .format('elem','sig_z','sig_r','sig_t','tau_zr','p1','p2','ang'),file=fout) for ne in range(0,nele): sigz =strs[0,ne] sigr =strs[1,ne] sigt =strs[2,ne] tauzr=strs[3,ne] ps1,ps2,ang=pst_ax4(sigz,sigr,tauzr) print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}' .format(ne+1,sigz,sigr,sigt,tauzr,ps1,ps2,ang),file=fout) fout.close() def dmat_ax4(E,po): d=np.zeros([4,4],dtype=np.float64) d[0,0]=1.0-po; d[0,1]=po ; d[0,2]=po ; d[0,3]=0.0 d[1,0]=po ; d[1,1]=1.0-po; d[1,2]=po ; d[1,3]=0.0 d[2,0]=po ; d[2,1]=po ; d[2,2]=1.0-po; d[2,3]=0.0 d[3,0]=0.0 ; d[3,1]=0.0 ; d[3,2]=0.0 ; d[3,3]=0.5*(1.0-2.0*po) d=d*E/(1.0+po)/(1.0-2.0*po) return d def bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4): bm=np.zeros([4,8],dtype=np.float64) #[N] nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 #dN/da,dN/db dn1a=-0.25*(1.0-b) dn2a= 0.25*(1.0-b) dn3a= 0.25*(1.0+b) dn4a=-0.25*(1.0+b) dn1b=-0.25*(1.0-a) dn2b=-0.25*(1.0+a) dn3b= 0.25*(1.0+a) dn4b= 0.25*(1.0-a) #Jacobi matrix and det(J) J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4 J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4 J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4 J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4 detJ=J11*J22-J12*J21 #[B]=[dN/dx][dN/dy] bm[0,0]= J22*dn1a-J12*dn1b bm[0,2]= J22*dn2a-J12*dn2b bm[0,4]= J22*dn3a-J12*dn3b bm[0,6]= J22*dn4a-J12*dn4b bm[1,1]=-J21*dn1a+J11*dn1b bm[1,3]=-J21*dn2a+J11*dn2b bm[1,5]=-J21*dn3a+J11*dn3b bm[1,7]=-J21*dn4a+J11*dn4b bm[2,1]=nm1/rr*detJ bm[2,3]=nm2/rr*detJ bm[2,5]=nm3/rr*detJ bm[2,7]=nm4/rr*detJ bm[3,0]=-J21*dn1a+J11*dn1b bm[3,1]= J22*dn1a-J12*dn1b bm[3,2]=-J21*dn2a+J11*dn2b bm[3,3]= J22*dn2a-J12*dn2b bm[3,4]=-J21*dn3a+J11*dn3b bm[3,5]= J22*dn3a-J12*dn3b bm[3,6]=-J21*dn4a+J11*dn4b bm[3,7]= J22*dn4a-J12*dn4b bm=bm/detJ return bm,detJ def ab_ax4(kk): if kk==0: a=-0.5773502692 b=-0.5773502692 if kk==1: a= 0.5773502692 b=-0.5773502692 if kk==2: a= 0.5773502692 b= 0.5773502692 if kk==3: a=-0.5773502692 b= 0.5773502692 return a,b def sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4): sm=np.zeros([8,8],dtype=np.float64) #Stiffness matrix [sm]=[B]T[D][B]*t*det(J) d=dmat_ax4(E,po) for kk in range(0,4): a,b=ab_ax4(kk) bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4) nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 sm=sm+np.dot(bm.T,np.dot(d,bm))*rr*detJ return sm def calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4): eps0 =np.zeros(4,dtype=np.float64) stress=np.zeros(4,dtype=np.float64) #stress vector {stress}=[D][B]{u} d=dmat_ax4(E,po) for kk in range(0,4): a,b=ab_ax4(kk) bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4) eps=np.dot(bm,wd) # Strain due to temperature change nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4 #Thermal strain eps0[0]=alpha*tem eps0[1]=alpha*tem eps0[2]=alpha*tem eps0[3]=0.0 stress=stress+np.dot(d,(eps-eps0)) stress=0.25*stress return stress def pst_ax4(sigz,sigr,tauzr): ps1=0.5*(sigz+sigr)+np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr) ps2=0.5*(sigz+sigr)-np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr) if sigz==sigr: if tauzr >0.0: ang= 45.0 if tauzr <0.0: ang=135.0 if tauzr==0.0: ang= 0.0 else: ang=0.5*np.arctan(2.0*tauzr/(sigz-sigr)) ang=180.0/np.pi*ang if sigz>sigr and tauzr>=0.0: ang=ang if sigz>sigr and tauzr <0.0: ang=ang+180.0 if sigz<sigr: ang=ang+90.0 return ps1,ps2,ang def tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4): eps0=np.zeros(4,dtype=np.float64) tfe =np.zeros(8,dtype=np.float64) # {tfe=[B]T[D]{eps0} d=dmat_ax4(E,po) for kk in range(0,4): a,b=ab_ax4(kk) bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4) # Strain due to temperature change nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4 #Thermal strain eps0[0]=alpha*tem eps0[1]=alpha*tem eps0[2]=alpha*tem eps0[3]=0.0 rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 tfe=tfe+np.dot(bm.T,np.dot(d,eps0))*rr*detJ return tfe def bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4): ntn=np.zeros([8,8],dtype=np.float64) nm =np.zeros([2,8],dtype=np.float64) bfe=np.zeros(8,dtype=np.float64) for kk in range(0,4): a,b=ab_ax4(kk) nm1=0.25*(1.0-a)*(1.0-b) nm2=0.25*(1.0+a)*(1.0-b) nm3=0.25*(1.0+a)*(1.0+b) nm4=0.25*(1.0-a)*(1.0+b) rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4 dn1a=-0.25*(1.0-b) dn2a= 0.25*(1.0-b) dn3a= 0.25*(1.0+b) dn4a=-0.25*(1.0+b) dn1b=-0.25*(1.0-a) dn2b=-0.25*(1.0+a) dn3b= 0.25*(1.0+a) dn4b= 0.25*(1.0-a) #Jacobi matrix and det(J) J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4 J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4 J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4 J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4 detJ=J11*J22-J12*J21 nm[0,0]=nm1 nm[0,2]=nm2 nm[0,4]=nm3 nm[0,6]=nm4 nm[1,1]=nm[0,0] nm[1,3]=nm[0,2] nm[1,5]=nm[0,4] nm[1,7]=nm[0,6] ntn=ntn+np.dot(nm.T,nm)*gamma*rr*detJ w=np.array([gkz,0,gkz,0,gkz,0,gkz,0],dtype=np.float64) bfe=np.dot(ntn,w) return bfe def main_ax4(): start=time.time() args = sys.argv fnameR=args[1] # input data file fnameW=args[2] # output data file nod=4 # Number of nodes per element nfree=2 # Degree of freedom per node # data input npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_ax4(fnameR,nod,nfree) # print out of input data prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node) # array declanation ir =np.zeros(nod*nfree,dtype=np.int) # Work vector for matrix assembly gk =np.zeros([nfree*npoin,nfree*npoin],dtype=np.float64) # Global stiffness matrix # assembly of stiffness matrix & load vector for ne in range(0,nele): i=node[0,ne]-1 j=node[1,ne]-1 k=node[2,ne]-1 l=node[3,ne]-1 m=node[4,ne]-1 z1=x[0,i]; r1=x[1,i] z2=x[0,j]; r2=x[1,j] z3=x[0,k]; r3=x[1,k] z4=x[0,l]; r4=x[1,l] E =ae[0,m] #E : Elastic modulus po =ae[1,m] #po : Poisson's ratio alpha =ae[2,m] #alpha: Thermal expansion coefficient gamma =ae[3,m] #gamma: Unit weight of material gkz =ae[4,m] #gkz : Acceleration in z-direction dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l] # temperature change sm=sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4) # element stiffness matrix tfe=tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4) # termal load vector bfe=bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4) # body force vector sm=sm*float(nzdir) tfe=tfe*float(nzdir) bfe=bfe*float(nzdir) ir[7]=2*l+1; ir[6]=ir[7]-1 ir[5]=2*k+1; ir[4]=ir[5]-1 ir[3]=2*j+1; ir[2]=ir[3]-1 ir[1]=2*i+1; ir[0]=ir[1]-1 for i in range(0,nod*nfree): it=ir[i] fp[it]=fp[it]+tfe[i]+bfe[i] for j in range(0,nod*nfree): jt=ir[j] gk[it,jt]=gk[it,jt]+sm[i,j] # treatment of boundary conditions for i in range(0,npoin): for j in range(0,nfree): if mpfix[j,i]==1: iz=i*nfree+j fp[iz]=0.0 for i in range(0,npoin): for j in range(0,nfree): if mpfix[j,i]==1: iz=i*nfree+j fp=fp-rdis[j,i]*gk[:,iz] gk[:,iz]=0.0 gk[iz,iz]=1.0 # solution of simultaneous linear equations #disg = np.linalg.solve(gk, fp) gk = csr_matrix(gk) disg = spsolve(gk, fp, use_umfpack=True) # recovery of restricted displacements for i in range(0,npoin): for j in range(0,nfree): if mpfix[j,i]==1: iz=i*nfree+j disg[iz]=rdis[j,i] # calculation of element stress wd =np.zeros(8,dtype=np.float64) strs =np.zeros([4,nele],dtype=np.float64) # Section force vector for ne in range(0,nele): i=node[0,ne]-1 j=node[1,ne]-1 k=node[2,ne]-1 l=node[3,ne]-1 m=node[4,ne]-1 z1=x[0,i]; r1=x[1,i] z2=x[0,j]; r2=x[1,j] z3=x[0,k]; r3=x[1,k] z4=x[0,l]; r4=x[1,l] wd[0]=disg[2*i]; wd[1]=disg[2*i+1] wd[2]=disg[2*j]; wd[3]=disg[2*j+1] wd[4]=disg[2*k]; wd[5]=disg[2*k+1] wd[6]=disg[2*l]; wd[7]=disg[2*l+1] E =ae[0,m] po =ae[1,m] alpha=ae[2,m] dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l] strs[:,ne]=calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4) # print out of result prout_ax4(fnameW,npoin,nele,disg,strs) # information dtime=time.time()-start print('n={0} time={1:.3f}'.format(nfree*npoin,dtime)+' sec') fout=open(fnameW,'a') print('n={0} time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout) fout.close() #============== # Execution #============== if __name__ == '__main__': main_ax4()
以 上