implicit real*8(a-h,o-z) dimension x(300000),y(300000),z(300000),e66(6,6),b(3,3)& &,mesh(3000000,8) character filename*128,filename2*48 !プログラムはメートル、MPa、 ell=2.0d0 h2=-2.0d-3!高さ(m) do i3=100,490,10 write(filename, '("hako", i4.4, ".inp")') i3! open (7,file=filename, status='replace') !寸法 b1=1.0d-3!幅(m) b2=8.0d-3!幅(m) h1=1.0d-3!高さ(m) h2=h2+1.0d-2!長さ p=1.0d-3 !木材のとき !ヤング率 !ezz=7.694d3 !exx=ezz/25.0d0!x軸方向のヤング率 !eyy=ezz/25.0d0!y軸方向のヤング率 !ezz2=206.0d3!z軸方向のヤング率 !ポアソン比の入力 !poixy=0.016d0 !poixz=poixy !poiyx=poixy !poiyz=poixy !poizx=0.4d0 !せん断弾性係数の入力 !gxy=ezz/15.0d0 !gxz=gxy !gyz=gxy !xyz方向の要素分割数 nx1=2 nx2=6 nx3=2 ny1=2 ny2=6 ny3=2 nz=1000 !1つあたりの要素の長さ? xn1=b1/nx1 xn2=b2/nx2 xn3=b1/nx3 yn1=h1/ny1 yn2=h2/ny2 yn3=h1/ny3 zn=ell/nz ! 節点番号ごとの節点座標!i~mで始まるのがdoで使える !左下 nitimen=((nx1+nx2+nx3+1)*(ny1+ny3+2)+(nx1+nx3+2)*(ny2-1)) nitidan1=(nx1+nx2+nx3+1) do k=0,nz!z do j=0,ny1! do i=0,nx1!x方向分割数を担当 nset=1+i+j*nitidan1+k*nitimen x(nset)=i*xn1 y(nset)=j*yn1 z(nset)=k*zn end do end do end do !真ん中下 do k=0,nz!z do j=0,ny1! do i=0,nx2-2!x方向分割数を担当 nset=nx1+2+i+j*nitidan1+k*nitimen x(nset)=b1+(1+i)*xn2 y(nset)=j*yn1 z(nset)=k*zn end do end do end do !右下 do k=0,nz!z do j=0,ny1! do i=0,nx3!x方向分割数を担当 nset=nx1+nx2+1+i+j*nitidan1+k*nitimen x(nset)=b1+b2+i*xn3 y(nset)=j*yn1 z(nset)=k*zn end do end do end do !左 nitidan2=nx1+nx3+2 do k=0,nz!z do j=0,ny2-2! do i=0,nx1!x方向分割数を担当 nset=((nx1+nx2+nx3+1)*(ny1+1))+1+i+j*nitidan2+k*nitimen x(nset)=i*xn1 y(nset)=h1+((1+j)*yn2) z(nset)=k*zn end do end do end do !右 do k=0,nz!z do j=0,ny2-2! do i=0,nx3!x方向分割数を担当 nset=((nx1+nx2+nx3+1)*(ny1+1))+nx1+2+i+j*nitidan2+k*nitimen x(nset)=b1+b2+i*xn3 y(nset)=h1+((1+j)*yn2) z(nset)=k*zn end do end do end do !左上 do k=0,nz!z do j=0,ny3! do i=0,nx1!x方向分割数を担当 nset=nitidan1*(ny1+1)+nitidan2*(ny2-1)+1+i+j*nitidan1+k*nitimen x(nset)=i*xn1 y(nset)=h1+h2+j*yn3 z(nset)=k*zn end do end do end do !真ん中上 do k=0,nz!z do j=0,ny3! do i=0,nx2-2!x方向分割数を担当 nset=nitidan1*(ny1+1)+nitidan2*(ny2-1)+nx1+2+i+j*nitidan1+k*nitimen x(nset)=b1+(1+i)*xn2 y(nset)=h1+h2+j*yn3 z(nset)=k*zn end do end do end do !右上 do k=0,nz!z do j=0,ny3! do i=0,nx3!x方向分割数を担当 nset=nitidan1*(ny1+1)+nitidan2*(ny2-1)+nx1+nx2+1+i+j*nitidan1+k*nitimen x(nset)=b1+b2+i*xn3 y(nset)=h1+h2+j*yn3 z(nset)=k*zn end do end do end do !i5~6が整数,aが文字,f9,6が小数.1peとかはEの*乗 !接点番号と座標の出力 nsettenzenbu=nitimen*(nz+1) write(7,'("*NODE, NSET=Nall")') do i=1,nsettenzenbu write(7,'(i6,a,f9.6,a,f9.6,a,f9.6)'),& &i,',',x(i),',',y(i),',',z(i) nsettensuu=i end do !要素番号とその要素の8接点の節点番号 !左下 nyou=0 do k=0,nz-1 do j=0,ny1 do i=1,nx1 nyou=nyou+1 mesh(nyou,1)=i+j*nitidan1+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan1 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan1 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !真ん中下 do k=0,nz-1 do j=0,ny1-1 do i=1,nx2 nyou=nyou+1 mesh(nyou,1)=nx1+i+j*nitidan1+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan1 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan1 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !右下 do k=0,nz-1 do j=0,ny1-1 do i=1,nx3 nyou=nyou+1 mesh(nyou,1)=nx1+nx2+i+j*nitidan1+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan1 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan1 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !左 do k=0,nz-1 do j=0,ny2-2 do i=1,nx1 nyou=nyou+1 mesh(nyou,1)=nitidan1*(ny1+1)+i+j*nitidan2+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan2 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan2 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !右 do k=0,nz-1 do j=0,ny2-2 do i=1,nx3 nyou=nyou+1 mesh(nyou,1)=nitidan1*ny1+nx1+nx2+i+j*nitidan2+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan2 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan2 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !左上 do k=0,nz-1 do j=0,ny3-1 do i=1,nx1 nyou=nyou+1 mesh(nyou,1)=nitidan1*(ny1+1)+nitidan2*(ny2-1)+i+j*nitidan1+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan1 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan1 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !真ん中上 do k=0,nz-1 do j=0,ny3-1 do i=1,nx2 nyou=nyou+1 mesh(nyou,1)=nitidan1*(ny1+1)+nitidan2*(ny2-1)+nx1+i+j*nitidan1+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan1 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan1 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !右上 do k=0,nz-1 do j=0,ny3 do i=1,nx3 nyou=nyou+1 mesh(nyou,1)=nitidan1*(ny1+1)+nitidan2*(ny2-2)+nx1+1+i+j*nitidan1+k*nitimen mesh(nyou,2)=mesh(nyou,1)+nitidan1 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou,4)=mesh(nyou,1)+nitimen mesh(nyou,5)=mesh(nyou,1)+1 mesh(nyou,6)=mesh(nyou,5)+nitidan1 mesh(nyou,7)=mesh(nyou,6)+nitimen mesh(nyou,8)=mesh(nyou,5)+nitimen end do end do end do !要素の出力 write(7,'("*ELEMENT, TYPE=C3D8, ELSET=E1")') do i=1,nyou write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)'),& &i,',',& &mesh(i,1),',',& &mesh(i,2),',',& &mesh(i,3),',',& &mesh(i,4),',',& &mesh(i,5),',',& &mesh(i,6),',',& &mesh(i,7),',',& &mesh(i,8) end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !固定端の節点群を書きだす write(7,'("*NSET,NSET=kotei")') do i=1,nitimen nkotei=i write(7,'(i6,a)'),nkotei end do !載荷点の節点群を書きだす write(7,'("*NSET,NSET=saika")') do i=1,nitimen nsaika=nitimen*nz+i write(7,'(i6,a)'),nsaika end do !境界条件 write(7,'("*BOUNDARY")') write(7,'("kotei,1,3")') ! 木材の設定 write(7,'("*MATERIAL,NAME=EL1")') write(7,'("*ELASTIC")') write(7,'("206.d3, 0.3")') write(7,'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1")') ! 載荷条件 write(7,'("*STEP")') write(7,'("*STATIC,SOLVER=SPOOLES")') !print'("*BUCKLE")' !print'("5,0.01")' write(7,'("*CLOAD")') write(7,'(A,a,i6,a,1pe13.5)'),'saika',',',2,',',p/nitimen ! datのオプション write(7,'("*NODE PRINT,NSET=saika")') write(7,'("U")') write(7,'("*EL PRINT,ELSET=Eall")') write(7,'("S")') write(7,'("*NODE FILE")') write(7,'("U")') write(7,'("*EL FILE,POSITION=AVERAGED AT NODES")') write(7,'("S")') write(7,'("*END STEP")') ! end do !jobの設定 do i=100,490,10 write(filename2, '("ccx_2.5 hako", i4.4,">hako", i4.4,".kekka" )') i,i open (8,file='job') write(8,*) filename2 end do close(8) end !######## 直交異方性の定数を出すプログラム ######## subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,& &gxz,gxy,poiyx,poizx,poizy) implicit real*8(a-h,o-z) dimension e66(6,6),b(3,3) !応力-ひずみ行列:e66 do i=1,6 do j=1,6 e66(i,j)=0.d0 end do end do ! e66(1,1)=1.d0/exx e66(1,2)=-poixy/exx e66(1,3)=-poixz/exx e66(2,1)=-poiyx/eyy e66(2,2)=1.d0/eyy e66(2,3)=-poiyz/eyy e66(3,1)=-poizx/ezz e66(3,2)=-poizy/ezz e66(3,3)=1.d0/ezz ! !****** 逆行列を取って、応力-ひずみ行列に変換する****** det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)& &+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)& &-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3) b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3)) b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3)) b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2)) b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3)) b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3)) b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2)) b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3)) b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3)) b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2)) do i=1,3 do j=1,3 e66(i,j)=b(j,i)/det end do end do ! !****** せん断補正係数を入力 ****** !write(7,*,'xy面のせん断弾性係数を入れて下さい' !read*, gxy e66(4,4)=gyz !write(7,*,'xz面のせん断弾性係数を入れて下さい' !read*, gxz e66(5,5)=gxz !write(7,*,'yz面のせん断弾性係数を入れて下さい' !read*, gyz e66(6,6)=gxy ! ! return end