c (薄肉とは見なせない)長方形断面のねじり定数、そりねじり定数を c 断面積分で求めるプログラム。 c ねじり応力関数は、 c ティモシェンコ・グーディア(金多 潔 監訳)『弾性論』(コロナ社) c 4版p.325の応力関数を用いた。 c 正規化そり関数は、 c 下関正義・藤沼平一・茅野寛奈:反りを考慮した長方形断面 c アイソパラメトリック要素, ばね論文集, 第39号, 1994, p.33-39 c の p.34 の正規化反り関数を用いた。この正規そり関数の2乗を c 断面積分して そりねじり定数を求めたけど、比較できる c 正解を知らないので、当たってるかどうか分からない。 c とはいえ、長辺が短辺の10倍以上の断面だと薄肉近似(後述) c の解に近づくから、まあ、そうおかしくもないのかも。 c (比較解を知ってる人がいたら教えて) c ねじり定数もそりねじり定数も、分割数を十分に細かくすれば、 c 短辺と長辺を逆に与えても同じ値になる。 c 有効数字3桁以上の精度が欲しいときは、おおよそ c 級数に関しては7項以上、 c 分割数に関しては 100*100 以上が目安でしょうか c c c 03/11/6版 c プログラムを作った人:後藤文彦 c c このプログラムは無保証です。 c このプログラムは改造/再配布して構いません。 c 但し、私的利用に留まらない場合は、もとのプログラムが c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/hari/kukei.f c から取ってきたものと分かるように明示して下さい。 c 改造/再配布後も上記の条件を保持して下さい。 c バグの報告は、 c http://www.str.ce.akita-u.ac.jp/~gotou/jpg/meeru.jpg c にお願いします。 c program kukei implicit real*8 (a-h, o-z) dimension x(5000),y(5000) c ねじりの応力関数、そり関数の級数の項数:kyuu c kyuu=19 5 continue write(*,*) '-------------------------------------------------' write(*,*) 'ねじり応力関数の級数の項数(奇数)を入れて下さい' read(*,*) kyuu write(*,*) '長方形の短辺と長辺を入力して下さい' read(*,*) tanp,tyou a=tanp/2.d0 b=tyou/2.d0 write(*,*) '短辺の分割数と長辺の分割数を偶数で入力して下さい' read(*,*) na,nb na=na/2 nb=nb/2 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) yi=yi+xl**2*f xi=xi+ym**2*f call neziri(kyuu,a,b,xl,ym,pi,p) pj=pj+p*f call sori(kyuu,a,b,xl,ym,pi,w) ww=ww+w**2*f c c (0<x<a,0<y<b),(-a<x<0,0<y<b),(0<x<a,-b<y<0),(-a<x<0,-b<y<0) c を別々に計算する場合(弾性なら符号が変わるだけだからw**2は c どの象限でも同じになる c call sori(kyuu,a,b,xl,ym,pi,w1) c call sori(kyuu,a,b,-xl,ym,pi,w2) c call sori(kyuu,a,b,xl,-ym,pi,w3) c call sori(kyuu,a,b,-xl,-ym,pi,w4) c ww=ww+(w1**2+w2**2+w3**2+w4**2)*f c 100 continue c c 各断面定数は、0<x<a,0<y<b の領域だけで積分したものを c 4倍してるけど、ファイバーモデルなどで各断面分割要素の弾塑性を c 考慮して弾性係数をかけながら積分するときは、ちゃんと c -a<x<a, -b<y<b の領域で断面積分すること。 xi=xi*4.d0 yi=yi*4.d0 pj=64.d0*a**2/pi**3*pj*4.d0 ww=ww*4.d0 c c write(*,500) ' 断面積:',tanp*tyou write(*,500) ' 強軸回りの断面二次モーメント(断面積分):',xi write(*,500) ' 強軸回りの断面二次モーメント( 正 解 ):', & tanp*tyou**3/12.d0 write(*,500) ' 弱軸回りの断面二次モーメント(断面積分):',yi write(*,500) ' 弱軸回りの断面二次モーメント( 正 解 ):', & tyou*tanp**3/12.d0 write(*,500) ' ねじり定数(断面積分): ',pj call timj(a,b,pi,tj) bt33=tyou*tanp**3/3.d0 write(*,500) ' ねじり定数( 正 解 ): ',tj*bt33 call timj2(kyuu,a,b,pi,tj2) write(*,500) ' ねじり定数(積分は正解、級数は有限個):',tj2 write(*,500) ' ねじり定数(薄肉の近似:bt^3/3): ',bt33 write(*,500) ' そりねじり定数(断面積分): ',ww write(*,500) ' そりねじり定数(薄肉近似:(bt)^3/144):', & (tyou*tanp)**3/144.d0 c (bt)^3/144 というのは _|_ 型断面の出っ張り | の厚さを c 0 としたそりねじり定数 goto 5 500 format(a,1pd15.7) end c c ******* ねじり定数 ******** c ねじり応力関数の級数について和を取る c x,y について断面積分はメインでやる subroutine neziri(kyuu,a,b,x,y,pi,p) implicit real*8 (a-h, o-z) p=0.d0 do 10 n=1,kyuu,2 p=p+(-1.d0)**( (n-1)/2 )/n**3 * & ( 1.d0- & cosh(n*pi*y/2.d0/a) / cosh(n*pi*b/2.d0/a) & ) *cos(n*pi*x/2.d0/a) 10 continue return end c c Timoshenko が計算した正解 c 手で積分して、級数も極限を取れる部分は取ってある c 極限を取れない級数は、31項まで subroutine timj(a,b,pi,tj) implicit real*8 (a-h, o-z) tj=0.d0 kyuu=31 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 積分は(Timoshenkoが)手で計算した正解だけど、 c 級数は、指定した項数までしか取らない subroutine timj2(kyuu,a,b,pi,tj2) implicit real*8 (a-h, o-z) tn=0.d0 tt=0.d0 tj2=0.d0 do 10 n=1,kyuu,2 tn=tn+1.d0/n**4 tt=tt+1.d0/n**5*tanh(n*pi*b/2.d0/a) 10 continue tj2=512.d0*a**3*b/pi**4*tn-1024.d0*a**4/pi**5*tt 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