! 両端支持に任意の面積で荷重をかけ解析するプログラム implicit real*8(a-h,o-z) dimension x(10000),y(10000),z(10000),es66(6,6),ew66(6,6) ! コンパイルして実行するとharioas.inpができる。 open(7,file='ccxbunpu.inp') ! 桁幅(x方向):b(m) ! 桁高(y方向):h(m) ! 梁の軸長(z方向):ell(m) psaika=1.d3*1.d-6 !(MN) b=1.d-2 h=1.d-2 ell=8.d-2 yban=1.d-2 !支承版の厚さ(m) ! ! x,y,z方向の要素分割数:nx,ny,nz nx=4 ny=2 nz=20 nsiten=2 !端部から支点の始まりまでの軸方向要素数 nshaba=4 !支点台座の幅の軸方向要素数 nphaba=4 !等分布荷重の載荷幅の軸方向要素数 ! ! ! 鋼板の材料定数 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 ! ! ! ! write(7,'("*NODE, NSET=Nall")') ! 1要素のx,y,z方向の長さ:xn,yn,zn xn=b/nx yn=h/ny zn=ell/nz ! ! 節点番号ごとの節点座標 ! 節点番号:nset nset=0 do k=0,nz do j=0,ny do i=0,nx nset=nset+1 x(nset)=xn*i y(nset)=yn*j z(nset)=zn*k write(7,'(i6,a,f9.3,a,f9.3,a,f9.3)')& & nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do ! nx1ny1=(nx+1)*(ny+1) !断面の節点数 ! ! ヒンジ支承板の始点: nhinzi nhinzi=nx1ny1*(nsiten+1)-nx ! ローラー支承板の始点: nrooraa nrooraa=nset-nx1ny1*(nshaba+nsiten)-nx ! 載荷部(軸方向中央)線の始点:nsaika nsaika=nx1ny1*(nz/2-nphaba/2)+1 ! 支承板の板面の節点数:nbanmen nbanmen=(nx+1)*(2*nshaba+1) ! 支承版をのぞく全節点数:nzen nzen=(nx+1)*(ny+1)*(nz+1) ! ! ヒンジ支承板下部の節点番号ごとの節点座標 ! 節点番号:nset do k=0,nshaba do i=0,nx nset=nset+1 x(nset)=xn*i y(nset)=h+yban z(nset)=zn*(nsiten+k) write(7,'(i6,a,f9.3,a,f9.3,a,f9.3)')& & nset,',',x(nset),',',y(nset),',',z(nset) end do end do ! ローラ支承板下部の節点番号ごとの節点座標 ! 節点番号:nset do k=0,nshaba do i=0,nx nset=nset+1 x(nset)=xn*i y(nset)=h+yban z(nset)=zn*(nz-nsiten-nshaba+k) write(7,'(i6,a,f9.3,a,f9.3,a,f9.3)')& & nset,',',x(nset),',',y(nset),',',z(nset) end do end do ! 要素番号とその要素の8節点の節点番号 nyou=0 do k=1,nz do j=1,ny do i=1,nx nyou=nyou+1 n1=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i n2=n1+1 n3=n1+(nx+1) n4=n3+1 n5=n1+(nx+1)*(ny+1) n6=n2+(nx+1)*(ny+1) n7=n3+(nx+1)*(ny+1) n8=n4+(nx+1)*(ny+1) write(7,'("*ELEMENT, TYPE=C3D8, ELSET=mokuzai")') write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)')& & nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6 end do end do end do ! ! ! ヒンジ支承板部要素の ! 要素番号とその要素の8節点の節点番号 do k=1,nshaba do i=1,nx nyou=nyou+1 n1=nhinzi+(k-1)*nx1ny1+(i-1) n2=n1+1 n3=nzen+i+(k-1)*(nx+1) n4=n3+1 n5=n1+(nx+1)*(ny+1) n6=n2+(nx+1)*(ny+1) n7=n3+(nx+1) n8=n4+(nx+1) write(7,'("*ELEMENT, TYPE=C3D8, ELSET=kouhan")') write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)')& & nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6 end do end do ! ! ローラー支承板部要素の ! 要素番号とその要素の8節点の節点番号 do k=1,nshaba do i=1,nx nyou=nyou+1 n1=nrooraa+(k-1)*nx1ny1+(i-1) n2=n1+1 n3=nzen+i+(nshaba+k)*(nx+1) n4=n3+1 n5=n1+(nx+1)*(ny+1) n6=n2+(nx+1)*(ny+1) n7=n3+(nx+1) n8=n4+(nx+1) write(7,'("*ELEMENT, TYPE=C3D8, ELSET=kouhan")') write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)')& & nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6 end do end do ! 境界条件 nhinha=nzen+(nx+1)*nshaba/2+1 ! ヒンジ支承 write(7,'("*NSET,NSET=hinzi")') do i=nhinha,nhinha+nx write(7,'(i6,a)'),i,',' end do nrooha=nzen+(nx+1)*(nshaba+1)+(nx+1)*nshaba/2+1 ! ローラー支承 write(7,'("*NSET,NSET=rooraa")') do i=nrooha,nrooha+nx write(7,'(i6,a)'),i,',' end do ! 載荷点 write(7,'("*NSET,NSET=saika")') if(nphaba==0) then do i=nx1ny1*nz/2+1,nx1ny1*nz/2+1+nx write(7,'(i6,a)'),i,',' nsaikaten=nx+1 end do else do k=0,nphaba do i=0,nx is=nx1ny1*k+nsaika+i write(7,'(i6,a)'),is,',' nsaikaten=(nx+1)*(2*nphaba+1) end do end do end if ! ! 拘束について ! 1でx軸を、2でy軸を、3でz軸を固定する。 ! 4~6でそれぞれの回転も固定できる。 write(7,'("*BOUNDARY")') write(7,'("hinzi,1,3")') write(7,'("rooraa,1,2")') ! 材料条件 ! write(7,'("*MATERIAL,NAME=kouhan")') write(7,'("*ELASTIC,TYPE=ORTHO")') write(7,'(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),',' write(7,*) es66(4,4),',',293.d0 ! ! write(7,'("*MATERIAL,NAME=mokuzai")') write(7,'("*ELASTIC,TYPE=ORTHO")') write(7,'(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),',' write(7,*) ew66(4,4),',',293.d0 ! write(7,'("*SOLID SECTION,ELSET=mokuzai,MATERIAL=mokuzai")') write(7,'("*SOLID SECTION,ELSET=kouhan,MATERIAL=kouhan")') ! 載荷条件 ! 載荷する荷重の節点 write(7,'("*STEP")') write(7,'("*STATIC,SOLVER=SPOOLES")') ! CLOADで集中荷重を設定 write(7,'("*CLOAD")') !載荷する面積 menseki=(nx+1)*(2*oo+1) ! 荷重の大きさ ! 載荷する荷重の方向 y軸なので2 ! 載荷する荷重の大きさを面積で割ることで、全体に等分布に載荷される。 write(7,*) 'saika,2,', psaika/nsaikaten write(7,'("*NODE PRINT,NSET=Nall")') write(7,'("U")') write(7,'("*EL PRINT,ELSET=mokuzai")') write(7,'("S")') write(7,'("*EL PRINT,ELSET=kouhan")') 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 ! ! ひずみ−応力行列の逆数を取って、応力−ひずみ行列を求める。 ! 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 ! !