!オンサイト木橋の3点曲げのモデルを作るプログラム implicit real*8(a-h,o-z) dimension x(315500),y(315500),z(315500),e66(6,6),b(3,3) print'("*NODE, NSET=Nall")' !材料の寸法の設定 b1=1.2d-1!角材の幅 b2=9.0d-3!鋼板の幅 h1=1.2d-1!角材の高さ h2=5.0d-1!鋼板の高さ h3=h2-h1*2!間の鋼板の高さ ell=6.0d0/2.0!スパンの半分の長さ ellhaji=1.500d-1!出ている部分の長さ p=1.000d-1!荷重 ! x,y,z方向の要素分割数 nx1=2!木材の分割数(偶数) nx2=3!鋼板の分割数 ny1=5!上と下部分の高さ方向の分割数 ny2=10!間の鋼板の高さ方向の分割数 nz=50!軸方向の半分の分割数 nzhaji=2!出ている部分の分割数 !材料の数 nk=7!角材の数 nmk=2!鋼板の数 !ヤング率の入力 ezz=6.921d3!z軸方向のヤング率 exx=ezz/25.0d0!x軸方向のヤング率 eyy=ezz/25.0d0!y軸方向のヤング率 !ポアソン比の入力 poixy=0.016d0 poixz=poixy poiyx=poixy poiyz=poixy poizx=0.4d0 poizy=poizx !せん断弾性係数の入力 gxy=ezz/15.0d0 gxz=gxy gyz=gxy !mk=2!鋼板の数 ! 1要素のx,y,z方向の長さ:xn,yn,zn xn1=b1/nx1!木材の1要素あたりの幅方向長さ xn2=b2/nx2!鋼板の1要素あたりの幅方向長さ yn1=h1/ny1!上と下部分の1要素あたりの高さ方向の長さ yn2=h3/ny2!間の鋼板の1要素あたりの高さ方向の長さ zn=ell/nz!支点間の1要素あたりの長さ znhaji=ellhaji/nzhaji!出ている部分の1要素あたりの長さ ! 節点番号ごとの節点座標 ! 下と上部分の節点 鋼板 do m=1,2 do k=1,ny1+1 do l=1,nmk do j=1,nx2+1 do i=1,nz*2+nzhaji*2+1 nset=i+(j-1)*(nz*2+nzhaji*2+1)+(l-1)*(nz*2+nzhaji*2+1)*(nx1*nk+nx2)& &+(k-1)*(nz*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)+& &(m-1)*((nz*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*(ny1+1)+& &(nz*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1)) x(nset)=xn2*(j-1)+(l-1)*(b2+b1*nk) y(nset)=yn1*(k-1)+(m-1)*(h1+h3) if (i>=0.and.i<=nzhaji) then z(nset)=znhaji*(i-1) else if (i>nzhaji.and.i<=(nz*2+1+nzhaji)) then z(nset)=zn*(i-1-nzhaji)+ellhaji else z(nset)=znhaji*(i-1-nzhaji-nz*2)+ellhaji+ell*2.0 end if print'(i6,a,f9.6,a,f9.6,a,f9.6)',& &nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do end do end do !下と上部分の節点 木材 do m=1,2 do k=1,ny1+1 do j=1,nx1*nk-1 do i=1,nz*2+nzhaji*2+1 nset=i+(j-1)*(nz*2+nzhaji*2+1)& &+(k-1)*(nz*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)+& &(m-1)*((nz*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*(ny1+1)+& &(nz*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1))+& &(nz*2+nzhaji*2+1)*(nx2+1) x(nset)=xn1*j+b2 y(nset)=yn1*(k-1)+(m-1)*(h1+h3) if (i>=0.and.i<=nzhaji) then z(nset)=znhaji*(i-1) else if (i>nzhaji.and.i<=(nz*2+1+nzhaji)) then z(nset)=zn*(i-1-nzhaji)+ellhaji else z(nset)=znhaji*(i-1-nzhaji-nz*2)+ellhaji+ell*2.0 end if print'(i6,a,f9.6,a,f9.6,a,f9.6)',& &nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do end do ! 真ん中部分の鋼材の節点 nset=(nz*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*(ny1+1) do k=1,ny2-1 do l=1,nmk do j=1,nx2+1 do i=1,nz*2+nzhaji*2+1 nset=nset+1 x(nset)=xn2*(j-1)+(l-1)*(b2+b1*nk) y(nset)=yn2*k+h1 if (i>=0.and.i<=nzhaji) then z(nset)=znhaji*(i-1) else if (i>nzhaji.and.i<=(nz*2+1+nzhaji)) then z(nset)=zn*(i-1-nzhaji)+ellhaji else z(nset)=znhaji*(i-1-nzhaji-nz*2)+ellhaji+ell*2.0 end if print'(i6,a,f9.6,a,f9.6,a,f9.6)',& &nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do end do ! ! 要素番号とその要素の8節点の節点番号 print'("*ELEMENT, TYPE=C3D8, ELSET=Eall1")' !下と上部分木材の要素 nyou=0 do m=1,2 do k=1,ny1 do j=1,nx1*nk do i=1,nz*2+nzhaji*2 nyou=nyou+1 n1=i+(j-1)*(nz*2+nzhaji*2+1)+(k-1)*(nz*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)& &+(m-1)*((nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+& &(nz*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1))+(nz*2+nzhaji*2+1)*nx2 n2=n1+1 n3=n1+(nz*2+nzhaji*2+1) n4=n3+1 n5=n1+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) n6=n2+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) n7=n3+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) n8=n4+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4 end do end do end do end do print'("*ELEMENT, TYPE=C3D8, ELSET=Eall2")' !鋼板の下と上部分 do m=1,2 do k=1,ny1 do l=1,nmk do j=1,nx2 do i=1,nz*2+nzhaji*2 nyou=nyou+1 n1=i+(j-1)*(nz*2+nzhaji*2+1)& &+(l-1)*(nz*2+nzhaji*2+1)*(nx2+nx1*nk)+(k-1)*(nz*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)& &+(m-1)*((nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+(nz*2+nzhaji*2+1)*(nx2+1)*& &nmk*(ny2-1)) n2=n1+1 n3=n1+(nz*2+nzhaji*2+1) n4=n3+1 n5=n1+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) n6=n2+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) n7=n3+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) n8=n4+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1) print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4 end do end do end do end do end do !鋼板の真ん中部分最初(例外) do l=1,nmk do j=1,nx2 do i=1,nz*2+nzhaji*2 nyou=nyou+1 n1=i+(j-1)*(nz*2+nzhaji*2+1)+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*& &ny1+(l-1)*(nz*2+nzhaji*2+1)*(nx2+nx1*nk) n2=n1+1 n3=n1+(nz*2+nzhaji*2+1) n4=n3+1 n5=i+(j-1)*(nz*2+nzhaji*2+1)+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+& &(l-1)*(nz*2+nzhaji*2+1)*(nx2+1) n6=n5+1 n7=n5+(nz*2+nzhaji*2+1) n8=n7+1 print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4 end do end do end do !鋼板の真ん中部分 do k=1,ny2-2 do l=1,nmk do j=1,nx2 do i=1,nz*2+nzhaji*2 nyou=nyou+1 n1=i+(j-1)*(nz*2+nzhaji*2+1)+(l-1)*(nz*2+nzhaji*2+1)*(nx2+1)+& &(k-1)*(nz*2+nzhaji*2+1)*(nx2+1)*nmk+& &(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1) n2=n1+1 n3=n1+(nz*2+nzhaji*2+1) n4=n3+1 n5=n1+(nz*2+nzhaji*2+1)*(nx2+1)*nmk n6=n2+(nz*2+nzhaji*2+1)*(nx2+1)*nmk n7=n3+(nz*2+nzhaji*2+1)*(nx2+1)*nmk n8=n4+(nz*2+nzhaji*2+1)*(nx2+1)*nmk print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4 end do end do end do end do ! !鋼板真ん中部分最後(例外) do l=1,nmk do j=1,nx2 do i=1,nz*2+nzhaji*2 nyou=nyou+1 n1=i+(j-1)*(nz*2+nzhaji*2+1)+(l-1)*(nz*2+nzhaji*2+1)*(nx2+1)& &+(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+& &(nz*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-2) n2=n1+1 n3=n1+(nz*2+nzhaji*2+1) n4=n3+1 n5=i+(j-1)*(nz*2+nzhaji*2+1)+(l-1)*(nz*2+nzhaji*2+1)*(nx2+nx1*nk)+& &(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+& &(nz*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1) n6=n5+1 n7=n5+(nz*2+nzhaji*2+1) n8=n7+1 print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4 end do end do end do ! 拘束節点の節点番号 print'("*NSET,NSET=kousoku1")' nx3=(nx2*nmk+nx1*nk+1) do i=1,nx3 nkousoku1=1+nzhaji+(nz*2+nzhaji*2+1)*(i-1) print'(i6,a)',nkousoku1,',' end do print'("*NSET,NSET=kousoku2")' do i=1,nx3 nkousoku2=1+nzhaji+nz*2+(nz*2+nzhaji*2+1)*(i-1) print'(i6,a)',nkousoku2,',' end do ! 載荷点の節点番号 print'("*NSET,NSET=saika")' do i=1,nx3 nsaika=(nz*2+nzhaji*2+1)*(i-1)+& &(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+& &(nz*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1)+& &(nz*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*ny1& &+nzhaji+nz+1 print'(i6,a)',nsaika,',' end do ! 載荷点の下の節点番号 print'("*NSET,NSET=saikasita")' do i=1,ny1+1 nsaikasita=(nz*2+nzhaji*2+1)*(nx1*nk/2+nx2)+& &(nz+nzhaji+1)+(nx2*nmk+nx1*nk+1)*& &(nz*2+nzhaji*2+1)*(i-1) print'(i6,a)',nsaikasita,',' end do ! 境界条件 print'("*BOUNDARY")' print'("kousoku1,1,3")' print'("kousoku2,1,2")' ! 材料の設定 ヤング率 ! 鋼板の設定 print'("*MATERIAL,NAME=EL1")' print'("*ELASTIC")' print'("206.d3, 0.3")' print'("*SOLID SECTION,ELSET=Eall2,MATERIAL=EL1")' ! 木材の設定 print'("*MATERIAL,NAME=EL2")' 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'("*SOLID SECTION,ELSET=Eall1,MATERIAL=EL2")' ! 載荷条件 print'("*STEP")' print'("*STATIC,SOLVER=SPOOLES")' print'("*CLOAD")' print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/nx3 ! 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