Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php on line 456
c harisen.f で等曲げを受ける円弧アーチの横ねじれ座屈を解かせるための c 入力データを作るプログラム。enkosen を実行させると enkosen.d という c データファイルができるから、harisen10とかの材料でこれをやると c ポアソン比が1を越えたりするので注意。 c つまりE/G>10とかの時点で等方性材料になり得ない。 poi=0.d0 c せん断補正係数：pak pak=10.d0*(1.d0+poi)/(12.d0+11.d0*poi) c せん断変形の影響を無視した答えが欲しいときは、 c pak=pak*1.d10 みたいに1より極端に大きくする c 軸長：ell ell=10.244d0 ell=10.d0 c 断面の桁幅：b b=0.1d0 b=5.d-2 c 断面の桁高：h h=0.5d0 h=25.d-2 c 断面積：a a=b*h c x軸（強軸）回りの断面二次モーメント：xi xi=b*h**3/12.d0 c Trahair の解と比較するために面内変位を無視するときは c xi を仮想的にyiに対して極端に大きくする c xi=xi*1.d6 c そのようにxiが極端に大きい場合は面内にたわまないから c その場合は曲げ面内(yz面）のせん断変形の影響も微少になるだろう c c y軸（弱軸）回りの断面二次モーメント：yi yi=h*b**3/12.d0 c 長方形断面のねじり定数:cj, そりねじり定数:wi call kukei(b,h,cj,wi) ea=e*a eix=e*xi eiy=e*yi gj=g*cj eiw=e*wi el=ell/float(nofe) c c ウラソフの解析解 call vlasov(mode,ell,eix,eiy,gj,eiw,th,pi,vln,vlp) write(*,*)' 修正Vlasov（負） =',vln print '(a,1pd11.3)',' ML/EIy=',vln*ell/eiy write(*,*)' 修正Vlasov（正） =',vlp print '(a,1pd11.3)',' ML/EIy=',vlp*ell/eiy eix2=eix*1.d6 call vlasov(mode,ell,eix2,eiy,gj,eiw,th,pi,vln,vlp) write(*,*)' Vlasov（負）=',vln print '(a,1pd11.3)',' ML/EIy=',vln*ell/eiy write(*,*)' Vlasov（正）=',vlp print '(a,1pd11.3)',' ML/EIy=',vlp*ell/eiy write(*,*)' 一発目の弧長を入れて下さい' read(*,*) r write(*,*)' 増分弧長を入れて下さい' read(*,*) rn c ### データの出力 ### open(7,file='enkosen.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) e write(7,300) g write(7,300) a write(7,300) xi write(7,300) yi write(7,300) cj write(7,300) wi write(7,300) el write(7,300) pak write(7,300) pak 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 c c 長方形断面のねじり定数とそりねじり定数を求める c 詳細は c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/#kukei c subroutine kukei(tan,tyou,cj,wi) implicit real*8 (a-h, o-z) dimension x(1000),y(1000) c ねじりの応力関数、そり関数の級数の項数：kyuu kyuu=7 a=tan/2.d0 b=tyou/2.d0 c 短辺/2の分割数:na c と長辺/2の分割数:nb（偶数） na=100 nb=100 an=a/float(na) bn=b/float(nb) f=an*bn pi=2.d0*asin(1.d0) do 10 i=1,na x(i)=an*i -an/2.d0 10 continue do 20 j=1,nb y(j)=bn*j -bn/2.d0 20 continue c xi=0.d0 yi=0.d0 pj=0.d0 ww=0.d0 do 100 l=1,na do 100 m=1,nb xl=x(l) ym=y(m) call sori(kyuu,a,b,xl,ym,pi,w) ww=ww+w**2*f 100 continue wi=ww*4.d0 c call timj(kyuu,a,b,pi,tj) bt33=tyou*tan**3/3.d0 cj=tj*bt33 write(*,500) ' ねじり定数（Timoshenkoの長方形断面）： ',cj write(*,500) ' ねじり定数（薄肉の近似：bt^3/3）： ',bt33 write(*,500) ' そりねじり定数（断面積分）： ',wi write(*,500) ' そりねじり定数（薄肉近似：(bt)^3/144）：', & (tyou*tan)**3/144.d0 c (bt)^3/144 というのは _|_ 型断面の出っ張り | の厚さを c 0 としたそりねじり定数 500 format(a,1pd15.7) return end c c ******* ねじり定数 ******** c Timoshenko が計算した正解 c 手で積分して、級数も極限を取れる部分は取ってある c 極限を取れない級数は、19項まで subroutine timj(kyuu,a,b,pi,tj) implicit real*8 (a-h, o-z) tj=0.d0 do 10 n=1,kyuu,2 tj=tj+1.d0/n**5*tanh(n*pi*b/2.d0/a) 10 continue tj=1.d0-192.d0/pi**5*a/b*tj return end c c c ******** そりねじり定数 ********** c 正規化そり関数の級数について和を取る subroutine sori(kyuu,a,b,x,y,pi,w) implicit real*8 (a-h, o-z) w=0.d0 do 10 n=1,kyuu,2 w=w+(-1.d0)**( (n+1)/2 )/ n**3 * & sinh(n*pi*y/2.d0/a) / cosh(n*pi*b/2.d0/a) & *sin(n*pi*x/2.d0/a) 10 continue w=32.d0*a**2/pi**3* w + x*y return end c