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 hariso.f で(単純支持された)柱のオイラー座屈を解かせるための c 入力データを作るプログラム。oiraso を実行させると oiraso.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.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 柱のオイラー座屈公式 eul=(mode*pi/ell)**2*eiy write(*,*) ' オイラーの座屈解：',eul c c write(*,*)' 一発目の弧長を入れて下さい' c read(*,*) r write(*,*)' 増分弧長を入れて下さい(目安はオイラー/1000)' write(*,*) ' Euler/1000=',eul/1000.d0 read(*,*) rn r=rn c ### データの出力 ### open(7,file='oiraso.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) 0.d0, 0.d0, 0.d0 20 continue write(7,300) e write(7,300) et write(7,300) sgmy write(7,300) sgmu write(7,300) g write(7,300) b write(7,300) h write(7,*) nhaba write(7,*) ntaka 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=1 write(7,150) nhikak c 載荷条件 c 載荷箇所は右端の1箇所 write(7,150) 1 c 右端の節点 nofe+1 では write(7,150) nofe+1 c Fz に -1 を載荷 write(7,400) 0.d0 write(7,400) 0.d0 write(7,400) -1.d0 write(7,400) 0.d0 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 θxもθyも拘束しなければ、普通は、弱軸回りで座屈すると思う。 c Ix=Iyみたいな断面の場合は、全節点で面外変位を拘束して c (平面問題にして)座屈する方向を決めてやらないと解が出ないかも。 c そり(ねじれ率)を拘束すべきかどうかは、適宜 判断。 c c c 局所系で与えたい境界条件はなし write(7,150) 0 c c 画面出力させたい要素 c 要素の左節点を出力させたいのが0箇所 write(7,150) 0 c 要素の右節点を出力させたいのが1箇所 write(7,150) 1 c 右端の要素 write(7,150) nofe c ファイル出力させたい要素（以下同様） write(7,150) 0 write(7,150) 1 write(7,150) nofe ccccc gnuplotで描かせる軸応力を出力させたい要素 ccccc write(7,150) nofe/2 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 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