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