!オンサイト木橋の等曲げのモデルを作るプログラム
implicit real*8(a-h,o-z)
dimension x(315500),y(315500),z(315500),e66(6,6),b(3,3)
print'("*NODE, NSET=Nall")'
!材料の寸法の設定
b1=1.2d-1!角材の幅
b2=9.0d-3!鋼板の幅
h1=1.2d-1!角材の高さ
h2=5.0d-1!鋼板の高さ
h3=h2-h1*2!間の鋼板の高さ
ell1=1.80d0/2.0d0!支点から荷重までの距離
ell2=4.2d0/2.0d0!荷重間距離の半分
ellhaji=1.5d-1!出ている部分の距離
! x,y,z方向の要素分割数
nx1=2!木材の分割数
nx2=3!鋼板の分割数
ny1=5!上と下部分の高さ方向の分割数
ny2=3!間の鋼板の高さ方向の分割数
nz1=20!支点から荷重までのz方向の距離の分割数
nz2=30!荷重から中央までの分割数
nzhaji=2!出ている部分の分割数
!材料の数
nk=7!角材の数
nmk=2!鋼板の数
!ヤング率の入力
ezz=6.921d3!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
!
! 1要素のx,y,z方向の長さ：xn,yn,zn
xn1=b1/nx1!木材の1要素あたりの幅方向長さ
xn2=b2/nx2!鋼板の1要素あたりの幅方向長さ
yn1=h1/ny1!上と下部分の1要素あたりの高さ方向の長さ 
yn2=h3/ny2!間の鋼板の1要素あたりの高さ方向の長さ
zn1=ell1/nz1!支点と荷重の間の1要素あたりの長さ
zn2=ell2/nz2!荷重と荷重の間の1要素あたりの長さ
znhaji=ellhaji/nzhaji!出ている部分の1要素あたりの長さ
! 節点番号ごとの節点座標
! 下と上部分の節点　鋼板
do m=1,2
do k=1,ny1+1
do l=1,nmk
do j=1,nx2+1
do i=1,nz1*2+nz2*2+nzhaji*2+1
nset=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)+&
(l-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2)&
&+(k-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)+&
&(m-1)*((nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*(ny1+1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1))
x(nset)=xn2*(j-1)+(l-1)*(b2+b1*nk)
y(nset)=yn1*(k-1)+(m-1)*(h1+h3)
if (i>=0.and.i<=nzhaji) then
z(nset)=znhaji*(i-1)
else if (i>nzhaji.and.i<=(nz1+1+nzhaji)) then
z(nset)=zn1*(i-1-nzhaji)+ellhaji
else if (i>(nzhaji+nz1+1).and.i<nz1+nz2*2+nzhaji+1) then
z(nset)=zn2*(i-1-nzhaji-nz1)+ellhaji+ell1
else if (i>=(nzhaji+nz1+nz2*2+1).and.i<=(nz1*2+nz2*2+nzhaji+1)) then
z(nset)=zn1*(i-1-nzhaji-nz1-nz2*2)+ellhaji+ell1+ell2*2
else
z(nset)=znhaji*(i-nz1*2-nz2*2-nzhaji-1)+ellhaji+ell1*2+ell2*2
end if
print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
end do
!下と上部分の節点 木材
do m=1,2
do k=1,ny1+1
do j=1,nx1*nk-1
do i=1,nz1*2+nz2*2+nzhaji*2+1
nset=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)&
&+(k-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)+&
&(m-1)*((nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*(ny1+1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1))+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)
x(nset)=xn1*j+b2
y(nset)=yn1*(k-1)+(m-1)*(h1+h3)
if (i>=0.and.i<=nzhaji) then
z(nset)=znhaji*(i-1)
else if (i>nzhaji.and.i<=(nz1+1+nzhaji)) then
z(nset)=zn1*(i-1-nzhaji)+ellhaji
else if (i>(nzhaji+nz1+1).and.i<nz1+nz2*2+nzhaji+1) then
z(nset)=zn2*(i-1-nzhaji-nz1)+ellhaji+ell1
else if (i>=(nzhaji+nz1+nz2*2+1).and.i<=(nz1*2+nz2*2+nzhaji+1)) then
z(nset)=zn1*(i-1-nzhaji-nz1-nz2*2)+ellhaji+ell1+ell2*2
else
z(nset)=znhaji*(i-nz1*2-nz2*2-nzhaji-1)+ellhaji+ell1*2+ell2*2
end if
print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!  真ん中部分の鋼材の節点        
nset=(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*(ny1+1)
do k=1,ny2-1
do l=1,nmk
do j=1,nx2+1
do i=1,nz1*2+nz2*2+nzhaji*2+1
nset=nset+1
x(nset)=xn2*(j-1)+(l-1)*(b2+b1*nk)
y(nset)=yn2*k+h1
if (i>=0.and.i<=nzhaji) then
z(nset)=znhaji*(i-1)
else if (i>nzhaji.and.i<=(nz1+1+nzhaji)) then
z(nset)=zn1*(i-1-nzhaji)+ellhaji
else if (i>(nzhaji+nz1+1).and.i<nz1+nz2*2+nzhaji+1) then
z(nset)=zn2*(i-1-nzhaji-nz1)+ellhaji+ell1
else if (i>=(nzhaji+nz1+nz2*2+1).and.i<=(nz1*2+nz2*2+nzhaji+1)) then
z(nset)=zn1*(i-1-nzhaji-nz1-nz2*2)+ellhaji+ell1+ell2*2
else
z(nset)=znhaji*(i-nz1*2-nz2*2-nzhaji-1)+ellhaji+ell1*2+ell2*2
end if
print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!
! 要素番号とその要素の8節点の節点番号
print'("*ELEMENT, TYPE=C3D8, ELSET=Eall1")'
!下と上部分木材の要素
nyou=0
do m=1,2
do k=1,ny1
do j=1,nx1*nk
do i=1,nz1*2+nz2*2+nzhaji*2
nyou=nyou+1
n1=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)+(k-1)*&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)&
&+(m-1)*((nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1))&
&+(nz1*2+nz2*2+nzhaji*2+1)*nx2
n2=n1+1
n3=n1+(nz1*2+nz2*2+nzhaji*2+1)
n4=n3+1
n5=n1+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
n6=n2+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
n7=n3+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
n8=n4+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4
end do
end do
end do
end do
print'("*ELEMENT, TYPE=C3D8, ELSET=Eall2")'
!鋼板の下と上部分
do m=1,2
do k=1,ny1
do l=1,nmk
do j=1,nx2
do i=1,nz1*2+nz2*2+nzhaji*2
nyou=nyou+1
n1=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)&
&+(l-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx2+nx1*nk)+&
&(k-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)&
&+(m-1)*((nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)&
&*(ny1+1)+(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*&
&nmk*(ny2-1))
n2=n1+1
n3=n1+(nz1*2+nz2*2+nzhaji*2+1)
n4=n3+1
n5=n1+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
n6=n2+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
n7=n3+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
n8=n4+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4
end do
end do
end do
end do
end do
!鋼板の真ん中部分最初(例外)
do l=1,nmk
do j=1,nx2
do i=1,nz1*2+nz2*2+nzhaji*2
nyou=nyou+1
n1=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*&
&ny1+(l-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx2+nx1*nk)
n2=n1+1
n3=n1+(nz1*2+nz2*2+nzhaji*2+1)
n4=n3+1
n5=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+&
&(l-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)
n6=n5+1
n7=n5+(nz1*2+nz2*2+nzhaji*2+1)
n8=n7+1
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4
end do
end do
end do
!鋼板の真ん中部分
do k=1,ny2-2
do l=1,nmk
do j=1,nx2
do i=1,nz1*2+nz2*2+nzhaji*2
nyou=nyou+1
n1=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)&
&+(l-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)+&
&(k-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)
n2=n1+1
n3=n1+(nz1*2+nz2*2+nzhaji*2+1)
n4=n3+1
n5=n1+(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk
n6=n2+(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk
n7=n3+(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk
n8=n4+(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4
end do
end do
end do
end do
!
!鋼板真ん中部分最後(例外)
do l=1,nmk
do j=1,nx2
do i=1,nz1*2+nz2*2+nzhaji*2
nyou=nyou+1
n1=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)+&
&(l-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)&
&+(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-2)
n2=n1+1
n3=n1+(nz1*2+nz2*2+nzhaji*2+1)
n4=n3+1
n5=i+(j-1)*(nz1*2+nz2*2+nzhaji*2+1)+&
&(l-1)*(nz1*2+nz2*2+nzhaji*2+1)*(nx2+nx1*nk)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2*nmk+nx1*nk+1)*(ny1+1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1)
n6=n5+1
n7=n5+(nz1*2+nz2*2+nzhaji*2+1)
n8=n7+1
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&nyou,',',n1,',',n5,',',n6,',',n2,',',n3,',',n7,',',n8,',',n4
end do
end do
end do
! 拘束節点の節点番号
print'("*NSET,NSET=kousoku1")'
nx3=(nx2*nmk+nx1*nk+1)
do i=1,nx3
nkousoku1=1+nzhaji+(nz1*2+nz2*2+nzhaji*2+1)*(i-1)
print'(i6,a)',nkousoku1,','
end do
print'("*NSET,NSET=kousoku2")'
do i=1,nx3
nkousoku2=1+nzhaji+nz1+nz2*2+(nz1*2+nz2*2+nzhaji*2+1)*(i-1)
print'(i6,a)',nkousoku2,','
end do
! 載荷点の節点番号
print'("*NSET,NSET=saika")'
do m=1,2
do i=1,nx3
nsaika=(nz1*2+nz2*2+nzhaji*2+1)*(i-1)+&
&(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*(ny1+1)&
&+(nz1*2+nz2*2+nzhaji*2+1)*(nx2+1)*nmk*(ny2-1)&
&+(nz1*2+nz2*2+nzhaji*2+1)*(nx1*nk+nx2*nmk+1)*ny1&
&+nzhaji+nz1+1+(m-1)*nz2*2
print'(i6,a)',nsaika,','
end do
end do
!載荷点下の真ん中の節点
print'("*NSET,NSET=saikasita")'
do i=1,ny1+1
nsaikasita=(nz1*2+nz2*2+nzhaji*2+1)*(nx2+nx1*nk/2+1)&
&+(nz1+nzhaji+nz2+1)+&
&(i-1)*(nz1*2+nzhaji*2+nz2*2+1)*(nx2*nmk+nx1*nk+1)
print'(i6,a)',nsaikasita,','
end do
! 載荷断面の節点番号
print'("*NSET,NSET=mannnaka")'
nx4=(nx1*nk+nx2*nmk+1)*(ny1+1)*2+(nx2+1)*nmk*(ny2-1)
do i=1,nx4
ntyuuou=(i-1)*(nz1*2+nz2*2+nzhaji*2+1)&
&+nzhaji+nz1+nz2+1
print'(i6,a)',ntyuuou,','
end do
! 境界条件
print'("*BOUNDARY")'
print'("kousoku1,1,3")'
print'("kousoku2,1,2")'
! 材料の設定　ヤング率
! 鋼板の設定
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'("206.d3, 0.3")'
print'("*SOLID SECTION,ELSET=Eall2,MATERIAL=EL1")'
! 木材の設定
print'("*MATERIAL,NAME=EL2")'
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,ELSET=Eall1,MATERIAL=EL2")'
! 載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',1.000d-1/nx3
! datのオプション
print'("*NODE PRINT,NSET=saikasita")'
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
