!ifを使って要素を設定するプログラム
!長方形断面
implicit real*8(a-h,o-z)
dimension x(300000),y(300000),z(300000)&
&,narax(300000),naray(300000),naraz(300000)&
&,mesh(200000,8),e66(6,6)
!材料の寸法
b=1.0d0!幅
h=1.0d1!高さ
ell=3.5d1!軸長
p=(1.000d-1)/2
!ヤング率の入力
ezz=7.694d3!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
!要素分割数
nx=4
ny=8
nz=2
nitimen=(nx+1)*(ny+1)
!1要素あたりの長さ    
xn=b/nx
yn=h/ny
zn=ell/nz
!節点の設定
nset=0
do k=1,nz+1
do j=1,ny+1
do i=1,nx+1
nset=nset+1
x(nset)=xn*(i-1)
y(nset)=yn*(j-1)
z(nset)=zn*(k-1)
narax(nset)=i
naray(nset)=j
naraz(nset)=k
end do
end do
end do
!要素の設定
nyou=0
do i=1,nset
if (narax(i)>0.and.narax(i)<nx+1.and.&
    &naray(i)<ny+1.and.naraz(i)<nz+1) then
nyou=nyou+1
mesh(nyou,1)=i
end if
end do
do j=1,nset
do i=1,nyou
if (naraz(j)==(naraz(mesh(i,1)))&
&.and.narax(j)==narax(mesh(i,1))&
&.and.naray(j)==naray(mesh(i,1))+1) then
mesh(i,2)=j
end if
end do
end do
do j=1,nset
do i=1,nyou
if (naraz(j)==(naraz(mesh(i,1))+1)&
&.and.narax(j)==narax(mesh(i,1))&
&.and.naray(j)==naray(mesh(i,1))+1) then
mesh(i,3)=j
end if
end do
end do
do j=1,nset
do i=1,nyou
if (naraz(j)==(naraz(mesh(i,1))+1)&
&.and.narax(j)==narax(mesh(i,1))&
&.and.naray(j)==naray(mesh(i,1))) then
mesh(i,4)=j
end if
end do
end do
do j=1,nset
do i=1,nyou
if (naraz(j)==(naraz(mesh(i,1)))&
&.and.narax(j)==narax(mesh(i,1))+1&
&.and.naray(j)==naray(mesh(i,1))) then
mesh(i,5)=j
end if
end do
end do
do j=1,nset
do i=1,nyou
if (naraz(j)==(naraz(mesh(i,1)))&
&.and.narax(j)==narax(mesh(i,1))+1&
&.and.naray(j)==naray(mesh(i,1))+1) then
mesh(i,6)=j
end if
end do
end do
do j=1,nset
do i=1,nyou
if (naraz(j)==(naraz(mesh(i,1))+1)&
&.and.narax(j)==narax(mesh(i,1))+1&
&.and.naray(j)==naray(mesh(i,1))+1) then
mesh(i,7)=j
end if
end do
end do
do j=1,nset
do i=1,nyou
if (naraz(j)==(naraz(mesh(i,1))+1)&
&.and.narax(j)==narax(mesh(i,1))+1&
&.and.naray(j)==naray(mesh(i,1))) then
mesh(i,8)=j
end if
end do
end do
print'("*NODE, NSET=Nall")'
do i=1,nset
print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
&i,',',x(i),',',y(i),',',z(i)
end do
print'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")'
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),',',&
&mesh(i,4),',',&
&mesh(i,5),',',&
&mesh(i,6),',',&
&mesh(i,7),',',&
&mesh(i,8)
end do
!載荷接点、拘束接点の設定
print'("*NSET,NSET=kousoku")'
do i=1,nitimen
nkousoku=i
print'(i6,a)',nkousoku
end do
print'("*NSET,NSET=saika")'
do i=1,nitimen
nsaika=nitimen*nz+i
print'(i6,a)',nsaika
end do
!境界条件
print'("*BOUNDARY")'
print'("kousoku,1,3")'
! 木材の設定
print'("*MATERIAL,NAME=ELmoku")'
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'("*MATERIAL,NAME=ELkou")'
!print'("*ELASTIC")'
!print'("206.d3, 0.3")'
print'("*SOLID SECTION,ELSET=Emoku,MATERIAL=ELmoku")'
! 載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/nitimen
! datのオプション
print'("*NODE PRINT,NSET=saika")'
print'("U")'
print'("*NODE FILE")'
print'("U")'
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
