マトリクス構造解析第15回

Processing math: 100%

[Google] [ウィキペディア] [附属図書館] [ シラバス ] (構力I) (構力II) (構造メモ) (構造実験) (フォートラン) (後藤資料)

マトリクス構造解析II 第7回

（マトリクス構造解析第15回）

小さい字は補足説明なので、読み飛ばしてもいいです。

このページのオリジナルの作者は後藤文彦です。

マトリクス構造解析オンライン授業用テキスト
第15回オンライン授業

5要素トラスを骨組として解いてみる

II第5回(13回)で、図のような5要素トラスをトラス要素のプログラム(torasu2sj.f90)で解いたが、今回は、この同じ5要素トラスを、節点を剛結として骨組のプログラム(hone2sj.f90)で解いて比較してみたい。現代の一般的なトラスは、格点部にヒンジはなく、上弦材や下弦材から突き出したガセットプレートと言われるプレート部分に、トラス部材がボルトでしっかりと連結されているが、それでも、変位や部材力はトラスモデルとして計算した結果とよく合うと考えられている。それを確かめてみようという趣旨である。なので、第13回と同様に、ホームセンターに売っている10mm角の 1m程度の角材を想定し、第13回と同じく
$E=10$ GPa $=10000$ N/mm $^{2}$
$A=(10$ mm $)^{2}=100$ mm $^{2}$
$\ell=1$ m $=1000$ mm
$P=500$ N
とする。

まず、以下のプログラムをC:\TDM-GCC-64 の中にダウンロードする。第13回の5要素トラスを、 torasu2sj.f90に設定したときとほぼ同じ設定内容だた、一応、説明する。

hone2sj.f90(Windows日本語版用Shift JIS版)
hone2.f90 (UTF-8版。日本語版以外のWindowsの場合は、こちらの方が文字化けしにくいかもしれない)

要素数、節点数

pi=asin(1.d0)*2.d0
!
nyou=5 !要素数
nset=4 !節点数

修正不要だ。

剛性マトリクス

まず、以下のように書いてあるところを見つける。

! 例題1
! 要素1の剛性マトリクス
el=1000.d0 !1本のL(mm)
ea=10.d3*10.d0**2 !EA(N)
ei=10.d3*10.d0**4/12.d0 !EI(N/mm^2・mm^4)
s(1,2,2)=       ea/el   ; s(1,2,5)=      -ea/el
s(1,5,2)=      -ea/el   ; s(1,5,5)=       ea/el
s(1,1,1)= 12.d0*ei/el**3; s(1,1,3)= -6.d0*ei/el**2
s(1,1,4)=-12.d0*ei/el**3; s(1,1,6)= -6.d0*ei/el**2
s(1,3,1)= -6.d0*ei/el**2; s(1,3,3)=  4.d0*ei/el
s(1,3,4)=  6.d0*ei/el**2; s(1,3,6)=  2.d0*ei/el
s(1,4,1)=-12.d0*ei/el**3; s(1,4,3)=  6.d0*ei/el**2
s(1,4,4)= 12.d0*ei/el**3; s(1,4,6)=  6.d0*ei/el**2
s(1,6,1)= -6.d0*ei/el**2; s(1,6,3)=  2.d0*ei/el
s(1,6,4)=  6.d0*ei/el**2; s(1,6,6)=  4.d0*ei/el

これも修正不要だ。単位は力をN, 長さをmmで与えている。
$\ell=1000$ mmだから el=1000.d0
ヤング率は、10GPa=10kMPa=10kN/mm $^{2}=10\times 10^{3}$ N/mm $^{2}$ だから、10.d3となる。
$A=(10$ mm $)^{2}$ だから
$EA=10\times 10^{3}$ N/mm $^{2}\times(10$ mm $)^{2}$ よって、ea=10.d3*10.d0**2
1辺10mmの正方形の断面2次モーメントは、 $I=(10$ mm $)^{4}/12$ だから、 $EI=10\times 10^{3}$ N/mm $^{2}\times(10$ mm $)^{4}/12$ よって、
ei=10.d3*10.d0**4/12.d0

! 要素2~4の剛性マトリクス
do k=2,4
do i=1,6
do j=1,6
   s(k,i,j)=s(1,i,j)
end do
end do
end do

要素②から要素④までは、要素①と同じでいいから、これも修正不要。

! 要素5の剛性マトリクス
el=1000.d0*sqrt(2.d0) !1本のL(mm)
s(5,2,2)=       ea/el   ; s(5,2,5)=      -ea/el
s(5,5,2)=      -ea/el   ; s(5,5,5)=       ea/el
s(5,1,1)= 12.d0*ei/el**3; s(5,1,3)= -6.d0*ei/el**2
s(5,1,4)=-12.d0*ei/el**3; s(5,1,6)= -6.d0*ei/el**2
s(5,3,1)= -6.d0*ei/el**2; s(5,3,3)=  4.d0*ei/el
s(5,3,4)=  6.d0*ei/el**2; s(5,3,6)=  2.d0*ei/el
s(5,4,1)=-12.d0*ei/el**3; s(5,4,3)=  6.d0*ei/el**2
s(5,4,4)= 12.d0*ei/el**3; s(5,4,6)=  6.d0*ei/el**2
s(5,6,1)= -6.d0*ei/el**2; s(5,6,3)=  2.d0*ei/el
s(5,6,4)=  6.d0*ei/el**2; s(5,6,6)=  4.d0*ei/el

要素⑤は、 $\ell=1000\sqrt{2}$ mmとすればいいから、これも修正不要。

各要素の節点番号

第13回のときと同様に設定。

idari(1)=1; migi(1)=2
idari(2)=2; migi(2)=3
idari(3)=3; migi(3)=4
idari(4)=4; migi(4)=1  !ここだけ修正
idari(5)=1; migi(5)=3

各要素の回転角

第13回のときと同様に設定。

!各要素の回転角（上の左右節点番号がz軸に横たわる状態からの）
th(1)=0.d0
th(2)=pi/2.d0
th(3)=pi
th(4)=pi+pi/2.d0
th(5)=pi/4.d0

境界条件

第13回のときと同じ境界条件を与える。節点1と節点2は、回転は自由だから、xt()は拘束しない。

!境界条件
xv(1)=0.d0
xw(1)=0.d0
xv(2)=0.d0

載荷条件

第13回のときと同じ載荷条件を与える。

!載荷条件
!cx(2)=1.d0 ! 例題1
fz(3)=500!N

実行

まずは、 $\ell=1000$ mmとして梁要素のプログラム(hone2sj.f90)を実行してみる。第13回のときのトラス要素(torasu2sj.f90)のプログラムの結果も上に並べてみるが、変位については、手計算とほぼ一致するトラス要素が $v_{3}=0.49999..., w_{3}=1.91421...$ に対して、梁要素が $v_{3}=0.49988..., w_{3}=1.91360...$ と有効桁3桁ぐらいは合っている。節点力も、トラス要素は理論値通りの鉛直・水平に500Nずつとなるところが、梁要素だと、499.8Nぐらいだ。まあ、トラスモデルとほぼ変わらない。節点は剛結なので、モーメントも発生しているが、最も大きい値でも 48Nmm とかだ。 50mmの片持ち梁に1N( $\simeq$ 100gf)を載せたぐらいのモーメントだから、 500Nの荷重から考えて、大したことはない。これなら、トラスモデルでも十分に実用的予測に使えるのではないか。

では、10mm角の部材を、断面はそのままでもっと短くして $\ell=200$ mm にして解いてみる。上がトラス要素、下が梁要素である。トラス要素の方は、恐らく手計算にほぼ一致する答えが出ているのだと思うが、実際の剛結トラスに近いであろう梁要素の結果は、やや違いが出てきた。 $v_{3}$ はそれほどでもないが、 $w_{3}$ は、トラス要素が0.3828mmに対して、梁要素が1.461mmだから、 4倍ぐらいになっている。実物のトラスだったら剛結した方が剛性が大きくなると思うのに、梁モデルでは剛結した方が変位が大きく出るのは不思議な感じだ。トラス要素モデルというのは曲げ剛性 $EI$ が入ってないから、（現実にはあり得ないけど）部材は全く曲がらずに（言わば、曲げ剛性が無限大で）、伸び剛性 $EA$ による部材の伸び縮みだけで変位が生じている。つまりトラスモデルは、部材が全く曲がらないように理想化されたモデルだということに注意する必要がある。それに対して、節点が剛結された梁要素の方は曲げ剛性 $EI$ による曲げ変形も入っている。節点力も、理想的トラスモデルでは、500Nとなるところが、490とか482とかだから、まあ、数パーセントぐらいは違ってくる。節点のモーメントは、大きいところで943Nmmと約1000Nmmぐらいだが、 1000mmの片持ち梁に1N( $\simeq$ 100gf)を載せたぐらいのモーメントだから 500Nの荷重から考えると、まだ大したことはないかもしれない。このように、節点が剛結のトラスがトラスモデルとして扱えるかどうかは、部材の細長さに関係する。後期に鋼構造設計学の授業で「細長比」というパラメーターについて習うと思うが、鋼構造では部材の細長さを表す細長比はとても重要である。

先頭、目次、

プログラム課題7：

第13回でやった以下の5要素トラス（自分の学籍番号のもの）を、節点を剛結とみなして梁要素のプログラム(hone2sj.f90)で解き、節点変位や節点力がトラス要素(torasu2sj.f90)で解いたものとどれくらい変わるか考察せよ。トラスの諸元は第13回と同じとするが、長さ $\ell$ については、短くして(torasu2sj.f90とhone2sj.f90それぞれで) 解いたものと比較し、部材の細長さと、節点が剛結のトラスをトラスモデルとして扱うことの妥当性について考察せよ。
長さ $\ell$ の水平、鉛直の4本の部材が正方形を作り、長さ $\sqrt{2}\ell$ の斜材が1本入っている。
$E=10$ GPa $=10000$ N/mm $^{2}$
$A=(10$ mm $)^{2}=100$ mm $^{2}$
$\ell=200\sim 1000$ mm
$P=500$ N
とする。
WebClassの「プログラム課題7」から、 hone2sj.f90を修正したプログラム（ellは変えて計算するので）の代表的な 1つ（例えば hone2sj_ell500.f90 とか）と、その出力ファイル（例えば kekka_hone500.out とか）をそれぞれアップロードせよ。また、最後のコメント欄に、節点が剛結のトラスをトラスモデルとして扱う妥当性についての考察を書くこと。うまく実行できなかった場合は、どういう症状で実行できないのかの説明を書き込むこと。

学籍番号の末尾が1か7の人。

学籍番号の末尾が2か8の人。

学籍番号の末尾が3か9の人。

学籍番号の末尾が4か6の人。

学籍番号の末尾が5か0の人。

hone2.f90プログラムソース

implicit real*8(a-h, o-z)
dimension f(27),& !節点力ベクトル
& d(27),& !節点変位ベクトル
& x(27),& !境界条件（拘束節点に0, その他に1が入る)
& ss(27,27),& !全体剛性行列
& xv(9),xw(9),xt(9),& !節点ごとの変位境界条件:v,w,θ
& fy(9),fz(9),cx(9),& !節点ごとの荷重条件:S,N,M
& idari(9),migi(9),& !各要素の左節点番号、右節点番号
& s(9,6,6),& !要素剛性行列
& th(9),& !各要素の回転角
& t(6,6),& !座標変換行列（1要素ごとに計算）
& s66(6,6),& !要素剛性行列（1要素ごとに移し替える入れ物）
& ts(6,6),& !TK（1要素ごとに計算）
& tst(6,6) !TKT（1要素ごとに計算）
dimension snm(6) !節点力の出力のために追加
!
pi=asin(1.d0)*2.d0
!
nyou=5 !要素数
nset=4 !節点数
!
! 初期化
do n=1,nyou
do i=1,6
do j=1,6
s(n,i,j)=0.d0
end do
end do
end do
!
do i=1,27
d(i)=0.d0
do j=1,27
ss(i,j)=0.d0
end do
end do
!
do i=1,9
xv(i)=1.d0; xw(i)=1.d0; xt(i)=1.d0 !境界条件は1で初期化
fy(i)=0.d0; fz(i)=0.d0; cx(i)=0.d0
end do
!
! 例題1
! 要素1の剛性マトリクス
el=1000.d0 !1本のL(mm)
ea=10.d3*10.d0**2 !EA(N)
ei=10.d3*10.d0**4/12.d0 !EI(N/mm^2・mm^4)
s(1,2,2)=       ea/el   ; s(1,2,5)=      -ea/el
s(1,5,2)=      -ea/el   ; s(1,5,5)=       ea/el
s(1,1,1)= 12.d0*ei/el**3; s(1,1,3)= -6.d0*ei/el**2
s(1,1,4)=-12.d0*ei/el**3; s(1,1,6)= -6.d0*ei/el**2
s(1,3,1)= -6.d0*ei/el**2; s(1,3,3)=  4.d0*ei/el
s(1,3,4)=  6.d0*ei/el**2; s(1,3,6)=  2.d0*ei/el
s(1,4,1)=-12.d0*ei/el**3; s(1,4,3)=  6.d0*ei/el**2
s(1,4,4)= 12.d0*ei/el**3; s(1,4,6)=  6.d0*ei/el**2
s(1,6,1)= -6.d0*ei/el**2; s(1,6,3)=  2.d0*ei/el
s(1,6,4)=  6.d0*ei/el**2; s(1,6,6)=  4.d0*ei/el
! 要素2~4の剛性マトリクス
do k=2,4
do i=1,6
do j=1,6
   s(k,i,j)=s(1,i,j)
end do
end do
end do
! 要素5の剛性マトリクス
el=1000.d0*sqrt(2.d0) !1本のL(mm)
s(5,2,2)=       ea/el   ; s(5,2,5)=      -ea/el
s(5,5,2)=      -ea/el   ; s(5,5,5)=       ea/el
s(5,1,1)= 12.d0*ei/el**3; s(5,1,3)= -6.d0*ei/el**2
s(5,1,4)=-12.d0*ei/el**3; s(5,1,6)= -6.d0*ei/el**2
s(5,3,1)= -6.d0*ei/el**2; s(5,3,3)=  4.d0*ei/el
s(5,3,4)=  6.d0*ei/el**2; s(5,3,6)=  2.d0*ei/el
s(5,4,1)=-12.d0*ei/el**3; s(5,4,3)=  6.d0*ei/el**2
s(5,4,4)= 12.d0*ei/el**3; s(5,4,6)=  6.d0*ei/el**2
s(5,6,1)= -6.d0*ei/el**2; s(5,6,3)=  2.d0*ei/el
s(5,6,4)=  6.d0*ei/el**2; s(5,6,6)=  4.d0*ei/el
!
!do i=1,6
!print'(9f10.2)',(s(1,i,j),j=1,9)
!end do
!print*
!!
!
!
!各要素の左節点番号、右節点番号
idari(1)=1; migi(1)=2
idari(2)=2; migi(2)=3
idari(3)=3; migi(3)=4
idari(4)=1; migi(4)=4
idari(5)=1; migi(5)=3
!
!各要素の回転角（上の左右節点番号がz軸に横たわる状態からの）
th(1)=0.d0
th(2)=pi+pi/2.d0
th(3)=pi
th(4)=pi+pi/2.d0
th(5)=pi+pi*3.d0/4.d0
!
!
!境界条件
xv(4)=0.d0!; xw(4)=0.d0; xt(4)=0.d0
xv(3)=0.d0; xw(3)=0.d0!; xt(3)=0.d0
!
!
do i=1,9
x(3*i-2)=xv(i)
x(3*i-1)=xw(i)
x(3*i)=xt(i)
end do
!
!載荷条件
!cx(2)=1.d0 ! 例題1
fz(2)=200!N
!
do i=1,9
f(3*i-2)=fy(i)
f(3*i-1)=fz(i)
f(3*i)=cx(i)
end do
!
!
!要素ごとのTKTの計算
do n=1,nyou
!
! 座標変換マトリクス
!
call zahyou(t,th(n))
!
do i=1,6
do j=1,6
s66(i,j)=s(n,i,j)
end do
end do
!
call mxtmx(t,s66,ts)
!
call mxmx(ts,t,tst)
!
!do i=1,6
!print'(9f10.2)',(tst(i,j),j=1,9)
!end do
!print*
!!!
!!重ねあわせ
i1=3*idari(n)-2
i2=3*idari(n)-1
i3=3*idari(n)
m1=3*migi(n)-2
m2=3*migi(n)-1
m3=3*migi(n)
!
ss(i1,i1)=ss(i1,i1)+tst(1,1)
ss(i1,i2)=ss(i1,i2)+tst(1,2)
ss(i1,i3)=ss(i1,i3)+tst(1,3)
ss(i1,m1)=ss(i1,m1)+tst(1,4)
ss(i1,m2)=ss(i1,m2)+tst(1,5)
ss(i1,m3)=ss(i1,m3)+tst(1,6)
!
ss(i2,i1)=ss(i2,i1)+tst(2,1)
ss(i2,i2)=ss(i2,i2)+tst(2,2)
ss(i2,i3)=ss(i2,i3)+tst(2,3)
ss(i2,m1)=ss(i2,m1)+tst(2,4)
ss(i2,m2)=ss(i2,m2)+tst(2,5)
ss(i2,m3)=ss(i2,m3)+tst(2,6)
!
ss(i3,i1)=ss(i3,i1)+tst(3,1)
ss(i3,i2)=ss(i3,i2)+tst(3,2)
ss(i3,i3)=ss(i3,i3)+tst(3,3)
ss(i3,m1)=ss(i3,m1)+tst(3,4)
ss(i3,m2)=ss(i3,m2)+tst(3,5)
ss(i3,m3)=ss(i3,m3)+tst(3,6)
!         !
!
ss(m1,i1)=ss(m1,i1)+tst(4,1)
ss(m1,i2)=ss(m1,i2)+tst(4,2)
ss(m1,i3)=ss(m1,i3)+tst(4,3)
ss(m1,m1)=ss(m1,m1)+tst(4,4)
ss(m1,m2)=ss(m1,m2)+tst(4,5)
ss(m1,m3)=ss(m1,m3)+tst(4,6)
!
ss(m2,i1)=ss(m2,i1)+tst(5,1)
ss(m2,i2)=ss(m2,i2)+tst(5,2)
ss(m2,i3)=ss(m2,i3)+tst(5,3)
ss(m2,m1)=ss(m2,m1)+tst(5,4)
ss(m2,m2)=ss(m2,m2)+tst(5,5)
ss(m2,m3)=ss(m2,m3)+tst(5,6)
!
ss(m3,i1)=ss(m3,i1)+tst(6,1)
ss(m3,i2)=ss(m3,i2)+tst(6,2)
ss(m3,i3)=ss(m3,i3)+tst(6,3)
ss(m3,m1)=ss(m3,m1)+tst(6,4)
ss(m3,m2)=ss(m3,m2)+tst(6,5)
ss(m3,m3)=ss(m3,m3)+tst(6,6)
!
!
end do !n要素について
!!
!!
!do i=1,8
!print'(8f10.3)', (ss(i,j),j=1,8)
!end do
!print*
!!
!　境界条件を入れる
do i=1,nset*3
do j=1,nset*3
  ss(i,j)=x(i)*ss(i,j)
  ss(j,i)=x(i)*ss(j,i)
end do
end do
do i=1,nset*3
if(x(i)<1.d-3) then
    ss(i,i)=1.d0
end if
end do
!!
!!
!do i=1,9
!print'(9f10.2)',(ss(i,j),j=1,9)
!end do
!
!print*,'f=',(f(j),j=1,18)
!
!
call gausu(nset*3,ss,d,f)
!
!print*
!do i=1,8
!print'(8f10.3)', (ss(i,j),j=1,8)
!end do
!
do n=1,nset
print*,'節点番号：',n,'v=',d(n*3-2),'w=',d(n*3-1), 'th=',d(n*3)
end do
!追加
    print'(A)','                                                    曲げモーメントの単位に注意!'
   !要素ごとのTKTの計算
   do n=1,nyou
   !
   ! 座標変換マトリクス
   !
   call zahyou(t,th(n))
   !
   do i=1,6
   do j=1,6
   s66(i,j)=s(n,i,j)
   end do
   end do
   !
   call mxtmx(t,s66,ts)
   !
   call mxmx(ts,t,tst)
   !
    snm=0.0
    do i=1,6
      do j=1,3
         snm(i)=snm(i)+tst(i,j)*d((idari(n)-1)*3+j)
      enddo
      do j=4,6
         snm(i)=snm(i)+tst(i,j)*d((migi(n)-1)*3+j-3)
      enddo
    enddo
    print'(A,I2,A,I2,A,f15.10,A,I2,A,f15.10,A,I2,A,f15.10)','部材番号： ',n,&
        &'   S(',idari(n),')=',snm(1),'  ,N(',idari(n),')=',snm(2),'  ,M(',idari(n),')=',snm(3)
    print'(A,A,I2,A,f15.10,A,I2,A,f15.10,A,I2,A,f15.10)','             ',&
        &'   S(',migi(n),')=',snm(4),'  ,N(',migi(n),')=',snm(5),'  ,M(',migi(n),')=',snm(6)
   !
   end do !n要素について
!ここまで追加
!print*,'d=',(d(j),j=1,18)
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
! 以上がメインプログラム
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
!
subroutine mxmx(a,b,c)
implicit real*8(a-h, o-z)
dimension a(6,6), b(6,6), c(6,6)
!
!
do i=1,6
do j=1,6
c(i,j)=0.
 do k=1,6
  c(i,j)=c(i,j)+a(i,k)*b(k,j)
 end do
end do
end do
!
return
end
!
subroutine zahyou(t,theta)
implicit real*8(a-h, o-z)
dimension t(6,6)
do i=1,3
do j=1,3
 t(i+3,j)=0.d0
 t(i,j+3)=0.d0
end do
end do
t(1,1)= cos(theta)
t(1,2)= sin(theta)
t(2,1)=-sin(theta)
t(2,2)= cos(theta)
t(3,3)=1.d0
t(4,4)=t(1,1)
t(4,5)=t(1,2)
t(5,4)=t(2,1)
t(5,5)=t(2,2)
t(6,6)=1.d0
return
end
!
subroutine mxtmx(a,b,c)
implicit real*8(a-h, o-z)
dimension a(6,6), b(6,6), c(6,6)
!
do i=1,6
do j=1,6
c(i,j)=0.
 do k=1,6
  c(i,j)=c(i,j)+a(k,i)*b(k,j)
 end do
end do
end do
!
!
return
end
!
!
subroutine gausu(n,a,x,b)
implicit real*8(a-h, o-z)
!dimension a(n,n),x(n),b(n)
dimension a(27,27),x(27),b(27)
!
!ガウスの消去法の参考としたのは、
!名取亮「すうがくぶっくす12 線形計算」(朝倉書店)p.10-15
!!
!print*
!do i=1,8
!print'(8f10.3)', (a(i,j),j=1,8)
!end do
!
!
do k=1,n-1  !a(k,k)を消去
do i=k+1,n  !k+1行からn行まで
do j=k+1,n    !k+1列からn列まで
   a(i,j)=a(i,j)-a(k,j)*a(i,k)/a(k,k)
end do
   b(i)=b(i)-b(k)*a(i,k)/a(k,k)
end do
end do
!
! 後退代入
x(n)=b(n)/a(n,n)
do k=n,1,-1
   akjxj=0.d0
do j=k+1,n
   akjxj=akjxj+a(k,j)*x(j)
end do
   x(k)=(b(k)-akjxj)/a(k,k)
end do
!
return
end
!
!