! 両端支持に任意の面積で荷重をかけ解析するプログラム     

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='kai2.inp')

! 桁幅(x方向)：b(m)
! 桁高(y方向)：h(m)
! 梁の軸長(z方向)：ell(m)
      psaika=1.d3*1.d-6 !(MN)
      b=1.d-1
      h=1.d-1
      ell=1.0d0
      yban=1.d-3 !支承版の厚さ(m)
!
! x,y,z方向の要素分割数：nx,ny,nz
      nx=4
      ny=8
      nz=200
       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)*(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=ISO")')
 write(7,'("206.d3, 0.3, 293d0")')
! 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
!
!
