Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php on line 456
! まずは、beam8p.inpと同じデータを、要素番号と節点番号の ! 並び順を変えて任意の要素数で生成するプログラムを作ってみる。 implicit real*8(a-h,o-z) dimension x(10000),y(10000),z(10000),es66(6,6),ew66(6,6) ! 節点座標：x(節点番号),y(節点番号),z(節点番号) print'("**")' print'("** beam8p.inpと同じデータを要素番号と節点番号の")' print'("** 並び順(と要素数)を変えたデータ")' print'("**")' print'("*NODE, NSET=Nall")' ! 桁幅(x方向)：b(m) !bs=6.d-3 ! 鋼板の桁幅(x方向)：bs(m) bs=9.d-3 ! 鋼板の桁幅(x方向)：bs(m) bw=840.d-3/2 ! 木材の桁幅(x方向)：bw(m) hw=120.d-3 ! 木材の桁高(y方向)：hw(m) h=500.d-3 ! 桁高(y方向)：h(m) hs=h-hw*2 ! 鋼板中央部の高さ(m) ell=3.d0 ! 梁の軸長(z方向)：ell(m) p=0.1d0 !中央集中荷重(MN) ! !bs=2.d-3: bw=8.d-3 hw=4.d-3: h=10.d-3: hs=3.d-3: ell=10.d-3 ! ! x,y,z方向の要素分割数：nx,ny,nz nx=4 ny=6 nz=8 nxs=1 !鋼板の幅(x)方向の要素分割数 nxw=nx-nxs !木材の幅(x)方向の要素分割数 nyw=1 !木材の高さ(y)方向の要素分割数 nys=ny-nyw*2 !鋼板中央部の高さ(y)方向の要素分割数 nsiten=1 !梁両端から支点までの要素数 nmen=(nx+1)*(nyw+1)*2+(nxs+1)*(nys-1) !断面の全節点数 ! ! 1要素のx,y,z方向の長さ：xn,yn,zn xw=bw/nxw !木材要素1要素のx方向の長さ xs=bs/nxs !木材要素1要素のx方向の長さ yw=hw/nyw !木材要素1要素のy方向の長さ ys=hs/nys !鋼板中央部要素1要素のy方向の長さ zn=ell/nz !各要素の1要素のz方向の長さ ! ! print*,'要素数：',(2*nx*nyw+nxs*nys)*nz ! print*,'節点数：',(2*(nx+1)*(nyw+1)+(nxs+1)*(nys-1))*(nz+1) ! ! 鋼板の材料定数 es=206.d9 /1.d6 pois=0.3d0 gs=es/2.d0/(1.d0+pois) call zeromx(6,ew66) es66(1,1)= 1.d0/es ; es66(1,2)=-pois/es ; es66(1,3)=-pois/es es66(2,1)=-pois/es ; es66(2,2)= 1.d0/es ; es66(2,3)=-pois/es es66(3,1)=-pois/es ; es66(3,2)=-pois/es ; es66(3,3)= 1.d0/es es66(4,4)=gs !逆行列を取るのは左上の3*3だけなので es66(5,5)=gs !右下3*3は最初からGを入れておく es66(6,6)=gs call gyaku(es66) ! 木材の材料定数 ezz=6.d9 /1.d6 exx=ezz/25.d0 eyy=ezz/25.d0 poixy=0.016d0 poixz=poixy poiyx=poixy poiyz=poixy poizx=0.4d0 poizy=poizx gyz=ezz/15.d0 gxz=gyz gxy=gyz call zeromx(6,ew66) ew66(1,1)= 1.d0/exx ; ew66(1,2)=-poixy/exx ; ew66(1,3)=-poixz/exx ew66(2,1)=-poiyx/eyy ; ew66(2,2)= 1.d0/eyy ; ew66(2,3)=-poiyz/eyy ew66(3,1)=-poizx/ezz ; ew66(3,2)=-poizy/ezz ; ew66(3,3)= 1.d0/ezz call gyaku(ew66) ew66(4,4)=gyz ew66(5,5)=gxz ew66(6,6)=gxy ! !do i=1,6 ! write(*,*) (ew66(i,j),j=1,6) !end do ! ! ! 節点番号ごとの節点座標 ! 節点番号：nset nset=0 do k=0,nz do j=0,ny if(j<=nyw) then do i=0,nx nset=nset+1 if(i<=nxs) then x(nset)=xs*i y(nset)=yw*j z(nset)=zn*k print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)' & &,nset,',',x(nset),',',y(nset),',',z(nset) else x(nset)=xs*nxs+xw*(i-nxs) y(nset)=yw*j z(nset)=zn*k print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)' & &,nset,',',x(nset),',',y(nset),',',z(nset) end if end do else if(jny-nyw) then print'(i6,a)',nmen*(k-1)+(nx+1)*(nyw+1)+(nxs+1)*(nys-1)+(nx+1)*(j-nyw-nys),',' end if end do end do ! !! ヒンジ支承 !print'("*NSET,NSET=hinzi")' !nowari=nmen*(nsiten+1) !nhazime=nowari-nx !do i=nhazime,nowari ! print'(i6,a)',i,',' !end do !! ローラー支承 !print'("*NSET,NSET=rooraa")' !nowari=nset-nmen*nsiten !nhazime=nowari-nx !do i=nhazime,nowari ! print'(i6,a)',i,',' !end do ! ! 固定端 print'("*NSET,NSET=koteitan")' do i=1,nmen print'(i6,a)',i,',' end do ! ! 上記(NSET=taisyou)の拘束節点の拘束変位 print'("*BOUNDARY")' print'("taisyou,1")' ! print'("hinzi,2,3")' ! print'("rooraa,2")' print'("koteitan,2,3")' ! ! 材料条件 ! print'("*MATERIAL,NAME=kouhan")' print'("*ELASTIC,TYPE=ORTHO")' 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)',& es66(1,1),',',es66(1,2),',',es66(2,2),',',es66(1,3),',',& &es66(2,3),',',es66(3,3),',',es66(6,6),',',es66(5,5),',' print*,es66(4,4),',',293.d0 ! ! print'("*MATERIAL,NAME=mokuzai")' print'("*ELASTIC,TYPE=ORTHO")' 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)',& ew66(1,1),',',ew66(1,2),',',ew66(2,2),',',ew66(1,3),',',& &ew66(2,3),',',ew66(3,3),',',ew66(6,6),',',ew66(5,5),',' print*,ew66(4,4),',',293.d0 ! print'("*SOLID SECTION,ELSET=mokuzai,MATERIAL=mokuzai")' print'("*SOLID SECTION,ELSET=kouhan,MATERIAL=kouhan")' ! ! 載荷条件 print'("*NSET,NSET=saika")' !nhazime=nmen*(nz/2)+1 !nowari=nhazime+nx !do i=nhazime,nowari ! print'(i6,a)',i,',' !end do ! !自由端載荷 print'("*NSET,NSET=saika")' do i=nset-nmen+1,nset print'(i6,a)',i,',' end do ! ! print'("*STEP")' print'("*STATIC,SOLVER=SPOOLES")' print'("*CLOAD")' ! y方向に載荷 ! pn=p/(nx+1) !中央に線載荷 pn=p/nmen !自由端に面載荷 write(*,*)'saika,2,',pn ! print'("*NODE PRINT,NSET=Nall")' print'("U")' print'("*EL PRINT,ELSET=kouhan")' print'("S")' print'("*EL PRINT,ELSET=mokuzai")' print'("S")' print'("*NODE FILE")' print'("U")' print'("*EL FILE,POSITION=AVERAGED AT NODES")' print'("S")' print'("*END STEP")' ! ! ! end ! ! ! ひずみ−応力行列の逆数を取って、応力−ひずみ行列を求める。 ! subroutine gyaku(a) implicit real*8(a-h, o-z) dimension a(6,6), b(3,3) det=a(1,1)*a(2,2)*a(3,3)+a(2,1)*a(3,2)*a(1,3)& & +a(3,1)*a(1,2)*a(2,3)-a(2,1)*a(1,2)*a(3,3)& & -a(3,1)*a(2,2)*a(1,3)-a(1,1)*a(3,2)*a(2,3) b(1,1)=(-1.d0)**(1+1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3)) b(1,2)=(-1.d0)**(1+2)*(a(2,1)*a(3,3)-a(3,1)*a(2,3)) b(1,3)=(-1.d0)**(1+3)*(a(2,1)*a(3,2)-a(3,1)*a(2,2)) b(2,1)=(-1.d0)**(2+1)*(a(1,2)*a(3,3)-a(3,2)*a(1,3)) b(2,2)=(-1.d0)**(2+2)*(a(1,1)*a(3,3)-a(3,1)*a(1,3)) b(2,3)=(-1.d0)**(2+3)*(a(1,1)*a(3,2)-a(3,1)*a(1,2)) b(3,1)=(-1.d0)**(3+1)*(a(1,2)*a(2,3)-a(2,2)*a(1,3)) b(3,2)=(-1.d0)**(3+2)*(a(1,1)*a(2,3)-a(2,1)*a(1,3)) b(3,3)=(-1.d0)**(3+3)*(a(1,1)*a(2,2)-a(2,1)*a(1,2)) do i=1,3 do j=1,3 a(i,j)=b(j,i)/det end do end do return end ! ! subroutine zeromx(n,a) implicit real*8(a-h,o-z) dimension a(n,n) do i=1,n do j=1,n a(i,j)=0.d0 end do end do return end ! !