c hari.f で等曲げを受ける円弧アーチの横ねじれ座屈を解かせるための
c 入力データを作るプログラム。enko を実行させると enko.d という
c データファイルができるから、hari<enko.d と実行させればよい。
c
program enko
implicit real*8 (a-h,o-z)
dimension alp(256),gmm(256),phi(256)
write(*,*) '要素数を入れて下さい'
read(*,*) nofe
write(*,*)'座屈モードを入れて下さい'
read(*,*) mode
write(*,*)'外乱を入れる場合は、その倍率を入れて下さい'
write(*,*)'(0 を入れた場合は座屈点を捜します)'
read(*,*) gmode
write(*,*)'モード計算のための出力が要りますか 要1／不要0 '
read(*,*) modeo
write(*,*)'正の曲げなら+1、負の曲げなら-1を入力'
read(*,*) porn
write(*,*)'中心角を「度」で入れて下さい'
read(*,*) th
pi=2.d0*asin(1.d0)
th=th*pi/180.d0
alp(1)=th/2.d0
alp(2)=alp(1)-3.d0*th/2.d0/float(nofe)
do 5 i=3,nofe-1
alp(i)=alp(i-1)-th/float(nofe)
5 continue
alp(nofe)=-th/2.d0
c 断面定数
c x軸回りの曲げ試験で測定されるヤング率：ex
ex=200.d6
c y軸回りの曲げ試験で測定されるヤング率：ey
ey=ex
c ねじり試験で測定されるせん断弾性係数：g
g=77.2d6
c 軸長：ell
ell=10.244
c 断面積：a
a=9.288d-2
c x軸（強軸）回りの断面二次モーメント：xi
xi=1.1363d-4
c 座屈前の面内変位の影響を無視したVlasovの解と比較するときは
c xi=xi*1.d6 などとして仮想的に面内剛性を極端に大きくする
c
c y軸（弱軸）回りの断面二次モーメント：yi
yi=3.871d-5
c ねじり定数：ej
ej=5.89d-7
c そりねじり定数：wi
wi=5.5869d-7
c
eix=ex*xi
eiy=ey*yi
ea=a*sqrt(ex*ey)
eiw=wi*sqrt(ex*ey)
gj=g*ej
el=ell/float(nofe)
call vlasov(mode,ell,eix,eiy,gj,eiw,th,pi,vln,vlp)
write(*,*)' 修正Vlasov（負） =',vln
write(*,*)' 修正Vlasov（正） =',vlp
eix2=eix*1.d6
call vlasov(mode,ell,eix2,eiy,gj,eiw,th,pi,vln,vlp)
write(*,*)' Vlasov（負）=',vln
write(*,*)' Vlasov（正）=',vlp
write(*,*)' 一発目の弧長を入れて下さい'
read(*,*) r
write(*,*)' 増分弧長を入れて下さい'
read(*,*) rn
c ### データの出力 ###
open(7,file='enko.d')
write(7,150) nofe
write(7,150) mode
write(7,300) gmode
write(7,150) modeo
do 20 i=1,nofe
write(7,200) alp(i), 0.d0, 0.d0
20 continue
write(7,300) ea
write(7,300) eix
write(7,300) eiy
write(7,300) gj
write(7,300) eiw
write(7,300) el
c 比較対照解の種類：nhikak
c nhikak=1:圧縮単純梁のオイラー座屈の解
c nhikak=2：等曲げ円弧アーチの横ねじれ座屈のウラソフの解
c nhikak=3: 片持ち梁の横ねじれ座屈のトラヘアの解（面内変位無視）
nhikak=2
write(7,150) nhikak
c 載荷条件
c 載荷箇所は両端の二箇所
write(7,150) 2
c 左端の節点1では
write(7,150) 1
c Mxに -1*porn を載荷
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) -1.d0*porn
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) 0.d0
c 右端の節点 nofe+1 では
write(7,150) nofe+1
c Mx に +1*porn を載荷
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) +1.d0*porn
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) 0.d0
c 境界条件
c 両端の二箇所
write(7,150) 2
c 左端では
write(7,150) 1
c u,v,w,θzを拘束
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) 1.d0
write(7,400) 1.d0
write(7,400) 0.d0
write(7,400) 1.d0
c 右端では
write(7,150) nofe+1
c u,v,θzを拘束
write(7,400) 0.d0
write(7,400) 0.d0
write(7,400) 1.d0
write(7,400) 1.d0
write(7,400) 1.d0
write(7,400) 0.d0
write(7,400) 1.d0
c ただ、これだと面外回転の回転軸が鉛直軸回りになって
c しまって、Vlasovとかの座屈モードと一致しない。で、
c 局所系で与えたい境界条件
c 両端の二箇所
write(7,150) 2
c 左端の要素で
write(7,150) 1
c 左側の回転変位を
write(7,500) 0,1,0,0
c 右端の要素で
write(7,150) nofe
c 右側の回転変位を
write(7,500) 0,0,0,1
c 局所系で与える（アーチの半径軸回りに回転するようにする）
c
c 画面出力させたい要素
c 要素の左節点を出力させたいのが一箇所
write(7,150) 1
c 真ん中（より一つ右の要素）
write(7,150) int(float(nofe)/2.d0)+1
c 要素の右節点を出力させたいのはなし
write(7,150) 0
c ファイル出力させたい要素（以下同様）
write(7,150) 1
write(7,150) int(float(nofe)/2.d0)+1
write(7,150) 0
c 弧長
write(7,100) r
write(7,100) rn
close(7)
150 format(i3)
100 format(f15.5)
200 format(1p3d17.9)
300 format(1pd13.5)
400 format(1pd9.1)
500 format(4i2)
stop
end
c
c
c 修正Vlasov の解析解を求めるサブルーチン
subroutine vlasov(n,blen,eix,eiy,gj,eiw,th,pi,vln,vlp)
implicit real*8(a-h,o-z)
if(th.lt.1.d-6) th=1.d-6
r=blen/th
f=(n*pi/blen)**2*eiw+gj
a=eiy/eix-1.d0+(1.d0/eix-eiy/eix/eix)*f
b=(eiy+(1.d0-2.d0*eiy/eix)*f)/r
c=((n*pi/blen)**2-1.d0/r/r)*eiy*f
d=b**2-4.d0*a*c
vln=(-b+sqrt(d))/2.d0/a
vlp=(-b-sqrt(d))/2.d0/a
return
end