implicit real*8(a-h,o-z) dimension x(315500),y(315500),z(315500),e66(6,6),b(3,3)& &,meshmoku(300000,8),meshkou(300000,8) b1=8.4d-1 b2=0.9d-2 h1=1.2d-1 h2=2.6d-1 h=h1*2+h2 ell=3.5d0 p=1.000d-1!荷重 ellhaji=1.500d-1 ellall=ell*2+ellhaji*2 !xyz方向の要素分割数 nx1=8!木材 nx2=7!鋼板 nx=nx1+nx2*2+1 ny1=5!上部分と下部分の高さ方向の分割数 ny2=10!間の鋼板の高さ方向の分割数(3以上) nznaka=120 nzhaji=40 nz=nzhaji*2+nznaka*2 !ヤング率の入力 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 poizy=poizx !せん断弾性係数の入力 gxy=ezz/15.0d0 gxz=gxy gyz=gxy !1つあたりの要素の長さ? xn1=b1/nx1 xn2=b2/nx2 yn1=h1/ny1 yn2=h2/ny2 zn=ell/nznaka znhaji=ellhaji/nzhaji ! 節点番号ごとの節点座標!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) if (k>=0.and.k<=nzhaji) then z(nset)=znhaji*(k) else if (k>nzhaji.and.k<=(nznaka*2+nzhaji)) then z(nset)=zn*(k-nzhaji)+ellhaji else z(nset)=znhaji*(k-nzhaji-nznaka*2)+ellhaji+ell*2.0 end if 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) if (k>=0.and.k<=nzhaji) then z(nset)=znhaji*(k) else if (k>nzhaji.and.k<=(nznaka*2+nzhaji)) then z(nset)=zn*(k-nzhaji)+ellhaji else z(nset)=znhaji*(k-nzhaji-nznaka*2)+ellhaji+ell*2.0 end if 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 if (k>=0.and.k<=nzhaji) then z(nset)=znhaji*(k) else if (k>nzhaji.and.k<=(nznaka*2+nzhaji)) then z(nset)=zn*(k-nzhaji)+ellhaji else z(nset)=znhaji*(k-nzhaji-nznaka*2)+ellhaji+ell*2.0 end if end do end do end do end do !i5~6が整数,aが文字,f9,6が小数.1peとかはEの*乗 !接点番号と座標の出力 print'("*NODE, NSET=Nall")' do i=1,(nz+1)*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1)) print'(i6,a,f9.6,a,f9.6,a,f9.6)',& &i,',',x(i),',',y(i),',',z(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 print'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")' do i=1,2*nx2*(ny1*2+ny2)*nz print'(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 print'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")' do i=1,2*nx1*ny1*nz print'(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 nitimen=((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)) !載荷接点、拘束接点の設定 print'("*NSET,NSET=kousoku")' do i=1,nitimen*(nz+1) if (z(i)>ellhaji-0.000001d0.and.z(i)ellall-ellhaji-0.000001d0.and.z(i)ellall/2-0.000001d0.and.z(i)h-0.0000001d0) then nsaika=i print'(i6,a)',nsaika,',' end if end do ! 境界条件 print'("*BOUNDARY")' print'("kousoku,1,3")' print'("kousoku2,1,2")' ! 木材の設定 print'("*MATERIAL,NAME=ELmoku")' print'("*ELASTIC,TYPE=ORTHOROPIC")' call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,& & gxz,gxy,poiyx,poizx,poizy) print'(f12.4,a,f12.4,a,f12.4,a,& &f12.4,a,f12.4,a,f12.4,a,f12.4,& & a ,f12.4,a)',& &e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& &,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),',' print'(f12.4,a,a)',e66(4,4),',','293' print'("*MATERIAL,NAME=ELkou")' print'("*ELASTIC")' print'("206.d3, 0.3")' print'("*SOLID SECTION,ELSET=Emoku,MATERIAL=ELmoku")' print'("*SOLID SECTION,ELSET=Ekou,MATERIAL=ELkou")' ! 載荷条件 print'("*STEP")' print'("*STATIC,SOLVER=SPOOLES")' print'("*CLOAD")' print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/nx ! datのオプション print'("*NODE PRINT,NSET=saika")' print'("U")' print'("*NODE FILE")' print'("U")' print'("*END STEP")' ! ! ! 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