!木材の片持ち梁を先端に集中荷重をかけて解くプログラム
implicit real*8(a-h,o-z)
dimension x(1000000),y(1000000),z(1000000),e66(6,6),b(3,3)
! 桁幅(x方向)：b(m)
! 桁高(y方向)：h(m)
! 梁の軸長(z方向)：ell(m)
b1=5.0d-3
h=1.0d-3
ell=100.0d-3
ellhaji=100.0d-4
! x,y,z方向の要素分割数：nx,ny,nz
nxhan=3
nx=nxhan*2
nyhan=3
ny=nyhan*2
nz=5
nzhaji=2
nza=nz+nzhaji
p=1.0d-6
! 1要素のx,y,z方向の長さ：xn,yn,zn
xn=b1/nx
yn=h/ny
zn=ell/nz
znhaji=ellhaji/nzhaji!出ている部分の1要素あたりの長さ
!ヤング率の入力
ezz=5.690d3!z軸方向のヤング率
exx=ezz/25.0d0!x軸方向のヤング率
eyy=ezz/25.0d0!y軸方向のヤング率
!ポアソン比の入力
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
poizy=poizx
!せん断弾性係数の入力
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!
!print*,'要素数：',nx*ny*nz
!print*,'節点数：',(nx+1)*(ny+1)*(nz+1)
!節点番号ごとの節点座標
!節点番号：nset
print'("*NODE, NSET=Nall")'
nset=0
do k=1,nza+1
do j=0,ny
do i=0,nx
nset=nset+1
x(nset)=xn*i
y(nset)=yn*j
if (k>=0.and.k<=nz) then
z(nset)=zn*(k-1)
else
z(nset)=znhaji*(k-1-nz)+ell
end if
print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
&         nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
! 要素番号とその要素の8節点の節点番号
print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")'
nyou=0
do k=1,nza
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)
print'(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
!
! 拘束節点の節点番号3 z変位が0の節点
print'("*NSET,NSET=kotei3")'
nx1ny1=(nx+1)*(ny+1)
do i=1,nx1ny1
print'(i6,a)',i,','
end do
! 拘束節点の節点番号1 x変位が0の節点
print'("*NSET,NSET=kotei1")'
do i=1,ny+1
print'(i6,a)',1+nx/2+(i-1)*(nx+1),','
end do
! 拘束節点の節点番号2 y変位が0の節点
print'("*NSET,NSET=kotei2")'
do i=1,nx+1
print'(i6,a)',(nx+1)*ny/2+i,','
end do
! 節点ごとに異なる載荷条件を与える場合をふまえて、
print'("*NSET,NSET=saika")'
!
do k=1,nza
do i=1,nx+1
print'(i6,a)',i+(nx+1)*ny/2+k*(nx+1)*(ny+1),','
!print'(i6,a)',i+(nx+1)*ny/2+(k-1)*(nx+1)*(ny+1),','!拘束点に載荷しない版
end do
end do
!
! 境界条件
print'("*BOUNDARY")'
print'("kotei3,3")'
print'("kotei2,2")'
print'("kotei1,2")'
! 節点ごとに異なる拘束条件を与える場合は、
! 材料の設定　ヤング率
print'("*MATERIAL,NAME=EL")'
print'("*ELASTIC,TYPE=ORTHOROPIC")'
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy)
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)',&
&e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& 
&,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
print'(f12.4,a,a)',e66(4,4),',','293'
print'("*SOLID SECTION,MATERIAL=EL,ELSET=Eall")'
!
! 載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
! 自由端(つまり上で定義したLAST)の全節点(nx1ny1個)の
! y(つまり2)方向に大きさ1の荷重を分配して載荷したいなら、
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',(-1.0d0)*p/(nx+1)*(nza+1)
!print'(A,a,i6,a,1pe13.5)','saika',',',2,',',(-1.0d0)*p/(nx+1)*nza
! 節点ごとに異なる載荷条件を与える場合は、
print'("*NODE PRINT,NSET=saika")'
print'("U")'
!print'("*EL PRINT,ELSET=Eall")'
!print'("S")'
print'("*NODE FILE")'
print'("U")'
!print'("*EL FILE,POSITION=AVERAGED AT NODES")'
!print'("S")'
print'("*END STEP")'
end
!######## 直交異方性の定数を出すプログラム ########
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3)
!応力－ひずみ行列：e66
do i=1,6
do j=1,6
e66(i,j)=0.d0
end do
end do
!  
e66(1,1)=1.d0/exx
e66(1,2)=-poixy/exx
e66(1,3)=-poixz/exx
e66(2,1)=-poiyx/eyy
e66(2,2)=1.d0/eyy
e66(2,3)=-poiyz/eyy
e66(3,1)=-poizx/ezz
e66(3,2)=-poizy/ezz
e66(3,3)=1.d0/ezz
!
!****** 逆行列を取って、応力－ひずみ行列に変換する******
det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)&
&+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)&
&-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
do i=1,3
do j=1,3
e66(i,j)=b(j,i)/det
end do
end do
!
!****** せん断補正係数を入力 ******
!print*,'xy面のせん断弾性係数を入れて下さい'
!read*, gxy
e66(4,4)=gyz
!print*,'xz面のせん断弾性係数を入れて下さい'
!read*, gxz
e66(5,5)=gxz
!print*,'yz面のせん断弾性係数を入れて下さい'
!read*, gyz
e66(6,6)=gxy
!
!
return
end
