!円筒の圧縮(CCX)

implicit none
!#--入力
real :: r=32.85d-3 !  m :半径
real :: h=120.d-3  !  m :高さ
real :: t=1.d-4    !  m :厚さ
real :: p=1.d0     ! MN :荷重
integer :: nc=64   !    :円周方向分割数
!ncは4の倍数(固定端を滑れるようにしたから)
integer :: nh=100  !    :高さ方向分割数
!------
real :: x(1000000),y(1000000),z(1000000)
integer :: mesh(1000000,8)
real :: pi,th,dh
integer :: i,j,k,ndan,nset,nsetzen,nyou,nyouzen

pi=2.d0*asin(1.d0) ! rad :円周率
th=2.d0*pi/nc      ! rad :一辺の中心角
dh=h/nh            !   m :一要素の高さ

!---節点座標の計算
nset=0
do j=0,nh*2
   if (mod(j,2)==1) then
      do i=0,nc-1
         nset=nset+1
         x(nset)=r*cos(th*i)
         y(nset)=r*sin(th*i)
         z(nset)=dh/2*j
      enddo
   else
      do i=0,nc*2-1
         nset=nset+1
         x(nset)=r*cos(th/2*i)
         y(nset)=r*sin(th/2*i)
         z(nset)=dh/2*j
      enddo
   endif
enddo
nsetzen=nset
!---節点の出力
print'("*NODE, NSET=Nall")'
do i=1,nset
   print'(i6,a,f17.9,a,f17.9,a,f17.9)',&
        &i,',',x(i),',',y(i),',',z(i)
enddo
!---要素を構成する節点
ndan=nc*3
nyou=0
do j=0,nh-1
   do i=0,nc-2
      nyou=nyou+1
      mesh(nyou,1)=1+i*2+ndan*j
      mesh(nyou,2)=mesh(nyou,1)+2
      mesh(nyou,3)=mesh(nyou,2)+ndan
      mesh(nyou,4)=mesh(nyou,3)-2
      mesh(nyou,5)=mesh(nyou,1)+1
      mesh(nyou,6)=i+nc*2+ndan*j+2
      mesh(nyou,7)=mesh(nyou,3)-1
      mesh(nyou,8)=mesh(nyou,6)-1
   enddo
   i=nc-1
   nyou=nyou+1
   mesh(nyou,1)=1+i*2+ndan*j
   mesh(nyou,2)=mesh(nyou,1)-i*2
   mesh(nyou,3)=mesh(nyou,2)+ndan
   mesh(nyou,4)=mesh(nyou,3)+i*2
   mesh(nyou,5)=mesh(nyou,1)+1
   mesh(nyou,6)=1+nc*2+ndan*j
   mesh(nyou,7)=mesh(nyou,3)+nc*2-1
   mesh(nyou,8)=mesh(nyou,6)+i
enddo
nyouzen=nyou
!---要素の出力
print'("*ELEMENT, TYPE=S8, ELSET=Eall")'
do i=1,nyouzen
   print'(i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8)',&
       &i,',',&
       &mesh(i,1),',',&
       &mesh(i,2),',',&
       &mesh(i,3),',',&
       &mesh(i,4),',',&
       &mesh(i,5),',',&
       &mesh(i,6),',',&
       &mesh(i,7),',',&
       &mesh(i,8)
end do
!---固定端の節点群を出力
print'("*NSET,NSET=koteiz")'
do i=1,nc*2
   print'(i8,a)',i
enddo
print'("*NSET,NSET=koteiy")'
print'(i8,a,i8)',1,',',nc+1
print'("*NSET,NSET=koteix")'
print'(i8,a,i8)',nc/2+1,',',nc*3/2+1
!---載荷点の節点群を出力
print'("*NSET,NSET=saika")'
do i=nsetzen-nc*2+1,nsetzen
   print'(i7,a)',i
enddo
!---境界条件
print'("*BOUNDARY")'
print'("koteiz,1,3")'!下面を固定する場合
!print'("koteiz,3")'!下面を滑らせる場合
!print'("*BOUNDARY")'
!print'("koteiy,2")'
!print'("*BOUNDARY")'
!print'("koteix,1")'
print'("*BOUNDARY")'!上面を固定する場合
print'("saika,1,2")'
!#---材料定数と厚さ
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'("69d3, 0.3d0")'!アルミ缶69GPa
print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL1")'
print'(1pe13.5)',t
!---載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*BUCKLE")'
print'("100,0.01")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',3,',',-p/(nc*2)
!---出力の設定
print'("*NODE PRINT,NSET=Nall")'
print'("U")'
!print'("*EL PRINT,ELSET=Eall")'
!print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE")'
print'("S")'
print'("*END STEP")'
end
