! 直方体のエレメント、節点をだしてくれるプログラム ! 09/12/10 ! 節点座標:x(節点番号),y(節点番号),z(節点番号) implicit real*8(a-h,o-z) dimension x(10000000),y(10000000),z(1000000) print'("**")' print'("** 直方体に荷重をかけるプログラム")' print'("** プログラム")' print'("**")' print'("*NODE, NSET=Nall")' b=0.006 ! 直方体の横(x方向)(m) tf=0.06 ! 直方体の縦(y方向)(m) ell=1. ! 直方体の奥行き(z方向)(m) bairitu=1.d2 bairitu=1.d0 t=t*bairitu tf=tf*bairitu ell=ell*bairitu ! 要素分割 nz=10 ! 軸(z)方向の分割 nfx=6! 幅(x)方向の分割 nfy=6 ! 厚さ(y)方向の分割 ! 静的載荷問題の場合の荷重 pload pload=1.d0 ! 静的載荷で変位を求める問題か、座屈解析か ! nzakutu=0 !静的縮み nzakutu=1 !座屈 ! そり拘束するかどうか(上下フランジの拘束の有無) ! nsori=0 !そり拘束する nsori=1 !そり拘束しない ! ! ! 3種類の要素の1要素分のx,y,z方向の長さ zn=ell/nz ! z方向の長さ xnf=b/nfx ! x方向の長さ ynf=tf/nfy ! y方向の長さ ! ! print*,'要素数:',nx*ny*nz ! print*,'節点数:',(nx+1)*(ny+1)*(nz+1) ! 節点番号ごとの節点座標 ! 節点番号:nset nset=0 do k=0,nz do j=0,nfy do i=0,nfx nset=nset+1 x(nset)=xnf*i y(nset)=ynf*j z(nset)=zn*k print'(i6,a,1pd11.3,a,1pd11.3,a,1pd11.3)', & nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do ! print*,'nset=',nset ! ! ! ! 要素番号とその要素の8節点の節点番号 print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")' nyou=0 do k=1,nz do j=1,nfy do i=1,nfx nyou=nyou+1 n1=(nfx+1)*(nfy+1)*(k-1)+(nfx+1)*(j-1)+i n2=n1+1 n3=n1+(nfx+1) n4=n3+1 n5=n1+(nfx+1)*(nfy+1) n6=n2+(nfx+1)*(nfy+1) n7=n3+(nfx+1)*(nfy+1) n8=n4+(nfx+1)*(nfy+1) print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)', & nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6 end do end do end do ! ! ! 境界条件 ! 節点ごとに異なる拘束条件を与える場合をふまえて、 ! 拘束条件の名前をつける。 print'("*NSET,NSET=FIX1")' ! その拘束節点の節点番号 print'(i6,a,)',nfx+1,',' ! 上記(NSET=FIX1)の拘束節点の拘束変位 print'("*BOUNDARY")' print'("FIX1,1,3")' ! 2番めの拘束節点 print'("*NSET,NSET=FIX2")' ! その拘束節点の節点番号 ! ! do i=1, (nfx+1)*(nfy+1) ! print'(i6,a,)',i,',' !手前の端面をすべて拘束 ! end do ! do i=2,(nfy+1) print'(i6,a,)',(nfx+1)*i,',' !縁部の1辺を拘束の場合 ! print'(i6,a,)',(nfx+1)*i-nfx/2,',' !中立軸を拘束 end do ! 上記(NSET=FIX2)の拘束節点の拘束変位 print'("*BOUNDARY")' print'("FIX2,1")' print'("FIx2,3")' ! 拘束条件の名前を付ける print'("*NSET,NSET=FIX3")' !奥の端面の右上の点 ! その拘束節点の節点番号 print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1),',' ! 上記(NSET=FIX3)の拘束節点の拘束変位 print'("*BOUNDARY")' print'("FIX3,1,2")' ! 拘束節点の名前を付ける print'("*NSET,NSET=FIX4")' ! 拘束節点の節点番号 ! 奥の端面部分 do i=2,nfy+1 ! 奥の端面の縁の1辺を拘束の場合 print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i,',' ! 奥の端面の中立軸を拘束の場合 !print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i-nfx/2,',' end do ! 上記(NSET=FIX4)の拘束節点の拘束変位 print'("*BOUNDARY")' print'("FIX4,1")' print'("*MATERIAL,NAME=EL")' print'("*ELASTIC")' ! ヤング率とポアソン比 !e=210000. ! ヤング率(MPa) e=203000. e=e*1.d3 poi=0.3 ! ポアソン比 print*, e,',',poi print'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")' ! ! 載荷条件 ! 節点ごとに異なる載荷条件を与える場合をふまえて、 ! 載荷条件の種類に対して適当な名前(LASTとか)を付ける ! ! print'("*NSET,NSET=LAST")' ! 直方体の端面全面に載荷の場合 ! do i=1,(nfx+1)*(nfy+1) ! print'(i6,a,)',(nfx+1)*(nfy+1)*nz+i,',' ! end do ! !直方体端面の1辺(拘束辺)に載荷の場合 ! do i=1,nfy+1 ! print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i,',' ! end do ! !直方体端面の中立軸線に載荷の場合 do i=1,nfy+1 print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i-nfx/2,',' end do ! ! print'("*STEP")' if(nzakutu==0) then print'("*STATIC,SOLVER=SPOOLES")' print'("*CLOAD")' print*,'LAST,3,',(-1.d0/((nfy+1)*(nfx+1))) else if(nzakutu==1) then print'("*BUCKLE")' print*,'5,0.01' print'("*CLOAD")' print*,'LAST,3,',(-1.d0/((nfy+1)*(nfx+1))) yi=b*h**3/12.d0 pi=2.d0*asin(1.d0) euler=e*yi*(pi/2.d0/ell)**2 open(8,file='euler') write(8,*) euler close(8) end if ! ! ! ! print'("*NODE PRINT,NSET=Nall")' print'("U")' print'("*EL PRINT,ELSET=Eall")' print'("S")' print'("*NODE FILE")' print'("U")' print'("*EL FILE,POSITION=AVERAGED AT NODES")' print'("S")' print'("*END STEP")' ! ! ! end