!応急橋の片持ち梁を解くモデルを作るプログラム
implicit real*8(a-h,o-z)
dimension x(500000),y(500000),z(500000),e66(6,6),b(3,3)
character filename*128,filename2*48
h2=2.00d-1!h2を変化させる時
do i2=250,900,50!h2を変化させる時
!do i3=2,402,25!nzを増やす時用
write(filename, '("oukyukata", i3.3, ".inp")') i2!h2を変化させる時
!write(filename, '("oukyuz", i3.3, ".inp")') i3!nzを増やす時用
open  (7,file=filename, status='replace')
! 部材の寸法と数
b1=1.20d-1!木材一つの幅
b2=2.40d-2!合板一つの幅（厚み）
h1=1.20d-1!木材一つの高さ
h2=h2+5.00d-2!合板一つの高さ
h3=h2-h1*2!真ん中の合板部分の高さ
ell=5.000d0/2.0!スパンの半分
nk=6!角材の数
ng=5!合板の数
!材料定数
!ヤング率の入力
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
!
! x,y,z方向の要素分割数：nx,ny,nz
nx1=2!角材のx方向の分割数
nx2=2!合板のx方向の分割数
ny1=5!角材のx方向の分割数(合板の上と下)
ny2=10!3以上にすること合板のx方向の分割数(真ん中部分)
nz=400!z方向の分割数
!nz=i3!nzを変化させる時用
! 1要素のx,y,z方向の長さ：xn,yn,zn
xn1=b1/nx1!角材の1要素のx方向の長さ
xn2=b2/nx2!合板の1要素のx方向の長さ
yn1=h1/ny1!木材の1要素のy方向の長さ
yn2=h3/ny2!合板の真ん中ぶぶんの1要素のy方向の長さ
zn=ell/nz!全体の1要素のz方向の長さ
! 節点番号ごとの節点座標
write(7,'("*NODE, NSET=Nall")')
!下部分と上部分の角材の節点の設定
do m=1,2
do k=1,ny1+1
do l=1,nk
do j=1,nx1+1
do i=1,nz+1
nset=i+(j-1)*(nz+1)+(k-1)*(nz+1)&
&*(nk*nx1+ng*nx2+1)&
&+(l-1)*(nx1+nx2)*(nz+1)&
&+(m-1)*((nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)+&
&(nz+1)*(nx2+1)*ng*(ny2-1))
x(nset)=xn1*(j-1)+(b1+b2)*(l-1)
y(nset)=yn1*(k-1)+(m-1)*(h1+h3)
z(nset)=zn*(i-1)
write(7,'(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 l=1,ng
do j=1,nx2-1
do i=1,nz+1
nset=i+(nx1+1)*(nz+1)+(j-1)*(nz+1)+&
&(k-1)*(nz+1)*(nk*nx1+ng*nx2+1)&
&+(l-1)*(nx1+nx2)*(nz+1)&
&+(m-1)*((nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)+&
&(nz+1)*(nx2+1)*ng*(ny2-1))
x(nset)=xn2*j+b1+(b1+b2)*(l-1)
y(nset)=yn1*(k-1)+(m-1)*(h1+h3)
z(nset)=zn*(i-1)
write(7,'(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 k=1,ny2-1
do l=1,ng
do j=1,nx2+1
do i=1,nz+1
nset=i+(nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)&
&+(j-1)*(nz+1)+(k-1)*(nz+1)*(nx2+1)*ng&
&+(l-1)*(nz+1)*(nx2+1)
x(nset)=xn2*(j-1)+b1+(b1+b2)*(l-1)
y(nset)=yn2*k+h1
z(nset)=zn*(i-1)
write(7,'(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
!要素番号とその要素の節点番号の設定
write(7,'("*ELEMENT, TYPE=C3D8, ELSET=Eall")')
!下と上の要素
nyou=0
do m=1,2
do k=1,ny1
do j=1,nx1*nk+nx2*ng
do i=1,nz
nyou=nyou+1
n1=i+(j-1)*(nz+1)+(k-1)*(nz+1)*(nx1*nk+nx2*ng+1)&
&+(m-1)*((nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)+(nz+1)*(nx2+1)*ng*(ny2-1))
n2=n1+1
n3=n1+(nz+1)
n4=n3+1
n5=n1+(nz+1)*(nx1*nk+nx2*ng+1)
n6=n2+(nz+1)*(nx1*nk+nx2*ng+1)
n7=n3+(nz+1)*(nx1*nk+nx2*ng+1)
n8=n4+(nz+1)*(nx1*nk+nx2*ng+1)
write(7,'(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,ng
do j=1,nx2
do i=1,nz
nyou=nyou+1
n1=i+(nz+1)*(nx1*nk+nx2*ng+1)*ny1+(nz+1)*nx1+&
&(j-1)*(nz+1)+(l-1)*(nx1+nx2)*(nz+1)
n2=n1+1
n3=n1+nz+1
n4=n3+1
n5=i+(nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)&
&+(j-1)*(nz+1)+(l-1)*(nx2+1)*(nz+1)
n6=n5+1
n7=n5+nz+1
n8=n7+1
write(7,'(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,ng
do j=1,nx2
do i=1,nz
nyou=nyou+1
n1=i+(nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)&
&+(j-1)*(nz+1)+(l-1)*(nx2+1)*(nz+1)&
&+(k-1)*(nz+1)*(nx2+1)*ng
n2=n1+1
n3=n1+nz+1
n4=n3+1
n5=n1+(nz+1)*(nx2+1)*ng
n6=n2+(nz+1)*(nx2+1)*ng
n7=n3+(nz+1)*(nx2+1)*ng
n8=n4+(nz+1)*(nx2+1)*ng
write(7,'(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,ng
do j=1,nx2
do i=1,nz
nyou=nyou+1
n1=i+(nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)&
&+(nz+1)*(nx2+1)*ng*(ny2-2)&
&+(j-1)*(nz+1)+(l-1)*(nx2+1)*(nz+1)
n2=n1+1
n3=n1+nz+1
n4=n3+1
n5=i+(nz+1)*(nx1*nk+nx2*ng+1)*(ny1+1)&
+(nz+1)*(nx2+1)*ng*(ny2-1)+(nz+1)*nx1&
&+(j-1)*(nz+1)+(l-1)*(nx2+nx1)*(nz+1)
n6=n5+1
n7=n5+nz+1
n8=n7+1
write(7,'(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
! 境界条件
! 節点ごとに異なる拘束条件を与える場合をふまえて、
! 拘束条件の種類に対して適当な名前(FIXとか)を付ける
write(7,'("*NSET,NSET=kousoku")')
! 拘束節点の節点番号
do i=1,(nx1*nk+nx2*ng+1)*(ny1+1)*2+(nx2+1)*ng*(ny2-1)
nkousokusettenn=(nz+1)*(i-1)+1
write(7,'(i6,a)')nkousokusettenn,','
end do
! 上記(NSET=FIX)の拘束節点の拘束変位
write(7,'("*BOUNDARY")')
write(7,'("kousoku,1,3")')
! 材料の設定　ヤング率
write(7,'("*MATERIAL,NAME=EL")')
write(7,'("*ELASTIC,TYPE=ORTHOROPIC")')
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy)
write(7,'(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),','
write(7,'(f12.4,a,a)')e66(4,4),',','293'
write(7,'("*SOLID SECTION,MATERIAL=EL,ELSET=Eall")')
!
! 載荷条件
!
! 載荷節点の節点番号
write(7,'("*NSET,NSET=saika")')
do i=1,(nx1*nk+nx2*ng+1)*(ny1+1)*2+(nx2+1)*ng*(ny2-1)
nsaikasettenn=(nz+1)*(i-1)+nz+1
write(7,'(i6,a)')nsaikasettenn,','
end do
!
nx3=(nx1*nk+nx2*ng+1)*(ny1+1)*2+(nx2+1)*ng*(ny2-1)
!
write(7,'("*STEP")')
write(7,'("*STATIC,SOLVER=SPOOLES")')
write(7,'("*CLOAD")')
!
write(7,'(A,a,i6,a,1pe13.5)')'saika',',',2,',',2.008d-2/nx3
write(7,'("*NODE PRINT,NSET=saika")')
!write(7,'("*NODE PRINT,NSET=saika1")')
write(7,'("U")')!たわみを出力
write(7,'("*EL PRINT,ELSET=Eall")')
!write(7,'("E")')
!write(7,'("S")')
write(7,'("*NODE FILE")')
write(7,'("U")')
write(7,'("*EL FILE,POSITION=AVERAGED AT NODES")')
write(7,'("S")')
write(7,'("*END STEP")')
close(7)
end do
do i=250,900,50
write(filename2, '("ccx_2.3 oukyukata", i3.3,">oukyukata", i3.3,".kekka" )') i,i!h2を変化させる時
open  (8,file='jobkata')
write(8,*) filename2
end do
close(8)
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
