!onsi21220.f90改良版 !箱断面の板厚を自動で厚くして計算(断面は変わるが、断面の中立軸は一定) implicit real*8(a-h,o-z) dimension x(1000000),y(1315500),z(1315500),e66(6,6),b(3,3)& &,meshmoku(300000,8),meshkou(300000,8) character filename*128,filename2*48 ell=10.0d0 a=50.0d-2 b2=0.0 h1=0.0 p=1.000d-3!荷重 !ここを変えたら下の!jobの設定 のi3も同じく変えて do i3=100,180,10!1から100まで1刻み write(filename, '("hako15-ell", i3.3, ".inp")') i3! open (7,file=filename, status='replace') b2=b2+1.0d-2 h1=(2.0d0*a*b2)/(a-3.0d0*b2)!板厚の関係式 b1=a-b2 h2=a-h1 !xyz方向の要素分割数 nx1=10!木材 nx2=8!鋼板 ny1=10!上部分と下部分の高さ方向の分割数 ny2=8!間の鋼板の高さ方向の分割数(3以上) nz=500 !ヤング率の入力 !木材 ezz=206.0d3 exx=ezz!x軸方向のヤング率 eyy=ezz!y軸方向のヤング率 ezz2=206.0d3!z軸方向のヤング率 !ポアソン比の入力 poixy=0.3d0 poixz=poixy poiyx=poixy poiyz=poixy poizx=0.3d0 poizy=poizx !せん断弾性係数の入力 gxy=ezz/(4.0d0*(1.0d0+poixy)) gxz=gxy gyz=gxy !1つあたりの要素の長さ? xn1=b1/nx1 xn2=b2/nx2 yn1=h1/ny1 yn2=h2/ny2 zn=ell/nz ! 節点番号ごとの節点座標!i~mで始まるのがdoで使える !下と上部分の節点 鋼板! do m=0,1 do k=0,nz do j=0,ny1 do l=0,1 do i=0,nx2 nset=1+i+l*(nx2+nx1)+j*(nx2*2+nx1+1)& &+k*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1))& &+m*((nx2*2+nx1+1)*(ny1+1)+(nx2+1)*2*(ny2-1)) !上と下あるのでmを動かして表す !kは面の一面すべての節点を掛けてz方向に動かす x(nset)=i*xn2+l*(b1+b2) !iが0~xn2までいってi=0,l=1になる y(nset)=j*yn1+m*(h1+h2) z(nset)=k*zn end do end do end do end do end do !下と上部分の節点 木材 do m=0,1 do k=0,nz do j=0,ny1 do i=0,nx1-2 nset=nx2+1+1+i+j*(nx2*2+nx1+1)& &+k*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1))& &+m*((nx2*2+nx1+1)*(ny1+1)+(nx2+1)*2*(ny2-1)) !上と下2つあるのでmを動かして表す x(nset)=i*xn1+b2+xn1 !iが0~xn2までいってi=0,l=1になる y(nset)=j*yn1+m*(h1+h2) z(nset)=k*zn end do end do end do end do !真ん中部分の鋼材の節点 do k=0,nz do j=0,ny2-2 do l=0,1 do i=0,nx2 nset=(1+nx1+2*nx2)*(ny1+1)+1+i+l*(1+nx2)+j*2*(nx2+1)& &+k*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1)) x(nset)=i*xn2+l*(b1+b2) y(nset)=j*yn2+h1+yn2 z(nset)=k*zn end do end do end do end do !i5~6が整数,aが文字,f9,6が小数.1peとかはEの*乗 !接点番号と座標の出力 write(7,'("*NODE, NSET=Nall")') do i=1,(nz+1)*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1)) 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 m=0,1 do k=0,nz-1 do j=0,ny1-1 do i=1,nx1 nyou=nyou+1 meshmoku(nyou,1)=nx2+i+(2*nx2+nx1+1)*j& &+k*((2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))& &+m*((2*nx2+nx1+1)*(ny1+1)+2*(nx2+1)*(ny2-1)) meshmoku(nyou,2)=meshmoku(nyou,1)+(2*nx2+nx1+1) meshmoku(nyou,3)=meshmoku(nyou,2)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshmoku(nyou,4)=meshmoku(nyou,1)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshmoku(nyou,5)=meshmoku(nyou,1)+1 meshmoku(nyou,6)=meshmoku(nyou,5)+(2*nx2+nx1+1) meshmoku(nyou,7)=meshmoku(nyou,6)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshmoku(nyou,8)=meshmoku(nyou,5)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) end do end do end do end do !上部分と下部分の鋼板 nyou=0 do l=0,1 do m=0,1 do k=0,nz-1 do j=0,ny1-1 do i=1,nx2 nyou=nyou+1 meshkou(nyou,1)=i+(2*nx2+nx1+1)*j& &+k*((2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))& &+m*((2*nx2+nx1+1)*(ny1+1)+2*(nx2+1)*(ny2-1))& &+l*(nx2+nx1) meshkou(nyou,2)=meshkou(nyou,1)+(2*nx2+nx1+1) meshkou(nyou,3)=meshkou(nyou,2)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshkou(nyou,4)=meshkou(nyou,1)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshkou(nyou,5)=meshkou(nyou,1)+1 meshkou(nyou,6)=meshkou(nyou,5)+(2*nx2+nx1+1) meshkou(nyou,7)=meshkou(nyou,6)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshkou(nyou,8)=meshkou(nyou,5)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) end do end do end do end do end do !真ん中の鋼板(最初の一段目) do m=0,1 do k=0,nz-1 do i=1,nx2 nyou=nyou+1 meshkou(nyou,1)=i+(nx2*2+nx1+1)*ny1& &+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))& &+m*(nx2+nx1) meshkou(nyou,2)=(2*nx2+nx1+1)*(ny1+1)+i& &+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))& &+m*(nx2+1) meshkou(nyou,3)=meshkou(nyou,2)+(nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshkou(nyou,4)=meshkou(nyou,1)+(nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1) meshkou(nyou,5)=meshkou(nyou,1)+1 meshkou(nyou,6)=meshkou(nyou,2)+1 meshkou(nyou,7)=meshkou(nyou,3)+1 meshkou(nyou,8)=meshkou(nyou,4)+1 end do end do end do !真ん中の鋼板(2段目から最後の前まで) do m=0,1 do k=0,nz-1 do j=0,ny2-3!ny2>3 do i=1,nx2 nyou=nyou+1 meshkou(nyou,1)=(2*nx2+nx1+1)*(ny1+1)+i& &+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))& &+m*(nx2+1)+j*(nx2+1+nx2+1) meshkou(nyou,2)=meshkou(nyou,1)+nx2+1+nx2+1 meshkou(nyou,3)=meshkou(nyou,2)+(nx2+1)*(ny2-1)*2+(nx2*2+1+nx1)*(ny1+1)*2 meshkou(nyou,4)=meshkou(nyou,1)+(nx2+1)*(ny2-1)*2+(nx2*2+1+nx1)*(ny1+1)*2 meshkou(nyou,5)=meshkou(nyou,1)+1 meshkou(nyou,6)=meshkou(nyou,2)+1 meshkou(nyou,7)=meshkou(nyou,3)+1 meshkou(nyou,8)=meshkou(nyou,4)+1 end do end do end do end do ! !真ん中の鋼板(最後の段) do m=0,1 do k=0,nz-1 do i=1,nx2 nyou=nyou+1 meshkou(nyou,1)=(nx2*2+nx1+1)*(ny1+1)+(nx2+1)*(ny2-2)*2+i& &+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))& &+m*(nx2+1) meshkou(nyou,2)=(nx2*2+nx1+1)*(ny1+1)+(nx2+1)*(ny2-1)*2+i& &+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))& &+m*(nx2+nx1) meshkou(nyou,3)=meshkou(nyou,2)+(nx2*2+nx1+1)*2*(ny1+1)+(nx2+1)*(ny2-1)*2 meshkou(nyou,4)=meshkou(nyou,1)+(nx2*2+nx1+1)*2*(ny1+1)+(nx2+1)*(ny2-1)*2 meshkou(nyou,5)=meshkou(nyou,1)+1 meshkou(nyou,6)=meshkou(nyou,2)+1 meshkou(nyou,7)=meshkou(nyou,3)+1 meshkou(nyou,8)=meshkou(nyou,4)+1 end do end do end do write(7,'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")') do i=1,2*nx2*(ny1*2+ny2)*nz write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)'),& &i,',',& &meshkou(i,1),',',& &meshkou(i,2),',',& &meshkou(i,3),',',& &meshkou(i,4),',',& &meshkou(i,5),',',& &meshkou(i,6),',',& &meshkou(i,7),',',& &meshkou(i,8) end do nyousosuukou=2*nx2*(ny1*2+ny2)*nz write(7,'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")') do i=1,2*nx1*ny1*nz write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)'),& &i+nyousosuukou,',',& &meshmoku(i,1),',',& &meshmoku(i,2),',',& &meshmoku(i,3),',',& &meshmoku(i,4),',',& &meshmoku(i,5),',',& &meshmoku(i,6),',',& &meshmoku(i,7),',',& &meshmoku(i,8) end do !載荷接点、拘束接点の設定 write(7,'("*NSET,NSET=kousoku")') do i=1,(nx2*2+nx1+1)*2*(ny1+1)+(nx2+1)*(ny2-1)*2 nkousoku=i write(7,'(i6,a)'),nkousoku,',' end do write(7,'("*NSET,NSET=saika")') nsaikasettensuu=0 do i=1,nsettensuu !節点数 if (ell-0.00002-0.00001d0.and.y(i)<0.00001d0& &.and.b2+0.0001d0hako15-ell", i3.3,".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 ! !****** せん断補正係数を入力 ****** !print*,'xy面のせん断弾性係数を入れて下さい' !read*, gxy e66(4,4)=gyz !print*,'xz面のせん断弾性係数を入れて下さい' !read*, gxz e66(5,5)=gxz !print*,'yz面のせん断弾性係数を入れて下さい' !read*, gyz e66(6,6)=gxy ! ! return end