!片持ち梁のねじり(CCX)の非線形座屈解析
! 円柱(B32)、弾性体、等方性材料

implicit real*8(a-h,o-z)
dimension x(1000000),y(1000000),z(1000000),e66(6,6),b(3,3)&
&,mesh(30000000,8),u(1000000),v(1000000),w(1000000)
!プログラムはメートル、MPa
zakutu=0.8634147d-3!N・m:起こしたいモードの座屈荷重
nstep=1000
to=zakutu/0.8    !MN:ちょっと大きめに荷重をかける
husei=1000       !線形座屈の節点変位をこれで割る
!寸法
d=13.d-3!直径(m)
!円形なら上のd、矩形なら下のb1,h
b1=1.d-3!幅(m)
h=3.d-3!高さ(m)
ell=3.0d0!長さ(m)
p=0.0d-6!荷重 先端の引張力(MN)
!z方向の分割数
nz=1000!軸方向
!1つあたりの要素の長さ
pi=asin(1.d0)*2.d0
zn=ell/nz!軸方向
!材料定数
young=206.d3!ヤング率(MN)
poi=0.3!ポアソン比
!#datファイルのデータの読み飛ばし
do i=1,2056 !初めの読み飛ばす行数(モードによって変わる)
read(*,'()')
enddo
!datファイルの節点変位の読み込み
!do i=1,(nr*nc+1)*(nz+1) !節点数
do i=1,nz*2+1 !節点数
read(*,*) a,u(i),v(i),w(i)
enddo
! 節点番号ごとの節点座標!i~mで始まるのがdoで使える 
nset=0
do k=0,nz*2
nset=nset+1
x(nset)=u(nset)/husei
y(nset)=v(nset)/husei
z(nset)=k*zn/2.d0+w(nset)/husei
end do
nsettenzenbu=nset!全節点数
!接点番号と座標の出力
print'("*NODE, NSET=Nall")'
do i=1,nsettenzenbu
print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
    &i,',',x(i),',',y(i),',',z(i)
end do
!要素番号とその要素の3接点の節点番号
nyou=0
do k=0,nz-1
nyou=nyou+1
mesh(nyou,1)=k*2+1
mesh(nyou,2)=mesh(nyou,1)+1
mesh(nyou,3)=mesh(nyou,2)+1
end do
!要素の出力
print'("*ELEMENT, TYPE=B32, ELSET=E1")'
do i=1,nyou
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
    &i,',',&
    &mesh(i,1),',',&
    &mesh(i,2),',',&
    &mesh(i,3)
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!固定端の節点群を書きだす
print'("*NSET,NSET=kotei")'
print'("1")'
!載荷点の節点群を書きだす
print'("*NSET,NSET=saika")'
print'(i6)',nsettenzenbu
!梁中央の節点群を書きだす
print'("*NSET,NSET=tyuuou")'
print'(i6)',(nsettenzenbu+1)/2
!境界条件
print'("*BOUNDARY")'
print'("kotei,1,6")'
print'(A,a,i1,a,i1)','saika',',',1,',',2
print'(A,a,i1,a,i1)','saika',',',4,',',5
print'(A,a,i1)','tyuuou',',',1!梁中央
print'(A,a,i1)','tyuuou',',',4!梁中央
! 材料の設定
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'(f20.9,a,f10.6)',young,',',poi
!print'("*BEAM SECTION,ELSET=E1,MATERIAL=EL1,SECTION=RECT")'!矩形断面
!print'(f10.6,a,f10.6)',b1,',',h!矩形断面
print'("*BEAM SECTION,ELSET=E1,MATERIAL=EL1,SECTION=CIRC")'!円形断面
print'(f10.6,a,f10.6)',d,',',d!円形断面
print'(f10.6,a,f10.6,a,f10.6)',0.d0,',',1.d0,',',0.d0
! 出力ステップ数
print'("*TIME POINTS,NAME=T1")'
do i=0,nstep
ts=real(i)/nstep
print'(f12.4,a)',ts,','
enddo
! 載荷条件
print'("*STEP,NLGEOM,INC=10000")'
print'("*STATIC")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',6,',',to
print'(A,a,i6,a,1pe13.5)','saika',',',3,',',p
! datのオプション
print'("*NODE PRINT,NSET=Nall,TIMEPOINTS=T1")'
print'("U")'
!print'("*EL PRINT,ELSET=E1,TIMEPOINTS=T1")'
!print'("S")'
print'("*NODE FILE,TIMEPOINTS=T1")'
print'("U")'
print'("*EL FILE,TIMEPOINTS=T1")'
print'("S")'
print'("*END STEP")'
!
end

