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

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

オンライン用テキストの準備中:8/5(水)までに用意します。 小さい字は補足説明なので、読み飛ばしてもいいです。

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

目次

トラス例題(単位荷重法)

前回、プログラム(torasu2.f90)の 使い方をやったので、このプログラムを使って、今回は、 図のような、トラスを解いてみたい。 まず、正解を単位荷重法で求めておく。 伸び剛性は$EA$とするが、今回は、想像できるぐらいの大きさ、固さの 材料諸元で解いてみよう。 ホームセンターに売ってるような 10mm角の角材1m(斜材は$\sqrt{2}$m)で、図のようなトラスを作って、 50kgfぐらいの体重で引っ張ったぐらい (図のトラスを90$^{\circ}$回転させて壁に固定して、人がぶら下がったぐらい)の 問題を想定して、
$E=10$GPa$=10000$N/mm$^{2}$
$A=(10$mm$)^{2}=100$mm$^{2}$
$\ell=1$m$=1000$mm
$P=500$N
とする。
$\sum\uparrow=V_{1}+V_{2}=0$
$\sum\rightarrow=H_{1}+P=0\;$ よって、$H_{1}=-P$
$\sum_{1}\circlearrowleft=V_{2}\ell-P\ell=0\;$よって、$V_{2}=P$
よって、$V_{1}=-P$

元の構造の部材力

反力は求まったので、次に各要素の部材力を求める。 まず要素③、⑤、①を切断して、トラスを2つのピースに切り離し、 今回は左側のピースを取り出してつりあいを考える。

$\sum_{1}\circlearrowright=N_{③}\ell=0$よって、$N_{③}=0$
$\sum_{3}\circlearrowleft=N_{①}\ell-P\ell+P\ell=0$よって、$N_{①}=0$
$\sum\uparrow=\frac{N_{⑤}}{\sqrt{2}}-P=0\;$ よって、$N_{⑤}=\sqrt{2}P$

次に、要素④と③を切断して、節点4側のピースを取り出す。 $N_{③}=0$と上で求まっているから、 $N_{④}$しか力はないので、力のつりあいから、$N_{④}=0$

次に、要素①と②を切断して、節点2側のピースを取り出す。 $N_{①}=0$と上で求まっているから、 力のつりあいから、
$\sum\uparrow=P+N_{②}=0\;$よって、 $N_{②}=-P$

$w_{3}$を求める

まず、載荷点の$z$方向変位$w_{3}$を求めよう。 単位荷重法で求めるには、節点変位を求めたい節点の求めたい変位の方向に、 仮想単位荷重を載荷して、各節点力を求める。 元の構造の部材力の$P$に$1$を代入すればいいから、この仮想構造の部材力は、 以下のように求まる。
$\overline{N}_{①}=\overline{N}_{③}=\overline{N}_{④}=0$
$\overline{N}_{②}=-1$
$\overline{N}_{⑤}=\sqrt{2}$
第5回でやったように 公式から、
$$\overline{1}\cdot w_{3} =\int_{全部材}\frac{N\overline{N}}{EA}dz^{\ell}$$
$$=\frac{1}{EA}\left(N_{②}\overline{N}_{②}\ell+N_{⑤}\overline{N}_{⑤}\cdot \sqrt{2}\ell \right)$$ $$=\frac{1}{EA}\left( (-P)(-1)\ell + \sqrt{2}P\cdot\sqrt{2}\cdot\sqrt{2}\ell \right)$$ $$=(1+2\sqrt{2})\frac{P\ell}{EA}\simeq 3.8284271\frac{P\ell}{EA}$$
$\frac{P\ell}{EA}$に、冒頭に示した以下の諸元を代入すると
$E=10$GPa$=10000$N/mm$^{2}$
$A=(10$mm$)^{2}=100$mm$^{2}$
$\ell=1$m$=1000$mm
$P=500$N
$$\frac{P\ell}{EA}= \frac{500\text{N}\times 1000\text{mm}}{10000\text{N}/\text{mm}^{2} \times 100\text{mm}^{2}}=0.5\text{mm}$$
すると
$w_{3}=(1+2\sqrt{2})\times 0.5 \simeq 1.91421356\text{mm}$ つまり、だいたい 2mm ぐらいということだ。

$v_{3}$を求める。

$v_{3}$を求めるには、 節点3の$v_{3}$方向に、単位荷重を与えて各要素の部材力を求める。

この場合、反力は、$V_{2}=1$だけだから、

$\sum\uparrow=1+\overline{N}_{②}=0\;$ よって、 $\overline{N}_{②}=-1$
その他の要素の部材力はすべて0. 上で解いた 元の構造の要素②の部材力は、${N}_{②}=-P$
公式から、
$$\overline{1}\cdot v_{3} =\int_{全部材}\frac{N\overline{N}}{EA}dz^{\ell}$$
$$=\frac{1}{EA}\left(N_{②}\overline{N}_{②}\ell\right)$$ $$=\frac{1}{EA}\left( (-P)(-1)\cdot \ell \right)$$ $$=\frac{P\ell}{EA}=0.5\text{mm}$$

トラス例題(プログラム)

では、 前回のプログラム(torasu2sj.f90)を 使って、上のトラスを解いてみよう。 torasu2sj.f90を適当な別名(torasu2sj_gotou2.f90とか)で保存し、 書き換えていく。

要素数、節点数

今回は、要素数は5要素、節点数は4節点。

pi=asin(1.d0)*2.d0
nyou=5 !要素数
nset=4 !節点数
剛性マトリクス

今回は、ホームセンターで売ってる10mm角の角材1mぐらいの材料で作った トラスを体重ぐらいで引っ張るということで、
$E=10$GPa$=10000$N/mm$^{2}$
$A=(10$mm$)^{2}=100$mm$^{2}$
$\ell=1$m$=1000$mm
$P=500$N
を与える。 プログラムに与える単位系は、自分で決められるが、 今回は、小さいものを解くので、力はN, 長さはmmの単位系を使う。 そうすると、ヤング率($\frac{力}{{長さ}^{2}}$)は、N/mm$^{2}=$MPaで与えなければならない。 まず、要素①は長さ$\ell=1000$mmだから、

E=10.d3!MPa
ell=1000.d0!mm
a=10.d0**2!mm^2
sk1=E*a/ell
! 要素1の剛性マトリクス
s(1,2,2)= sk1; s(1,2,4)=-sk1
s(1,4,2)=-sk1; s(1,4,4)= sk1

要素②、③、④はすべて長さ$\ell=1000$mmだから、 sk2, sk3, sk4にはsk1を代入すればよい。

! 要素2の剛性マトリクス
sk2=sk1
s(2,2,2)= sk2; s(2,2,4)=-sk2
s(2,4,2)=-sk2; s(2,4,4)= sk2
! 要素3の剛性マトリクス
sk2=sk1
s(3,2,2)= sk3; s(3,2,4)=-sk3
s(3,4,2)=-sk3; s(3,4,4)= sk3
! 要素4の剛性マトリクス
sk4=sk1
s(4,2,2)= sk4; s(4,2,4)=-sk4
s(4,4,2)=-sk4; s(4,4,4)= sk4

要素⑤は長さ$\sqrt{2}\ell=1000\sqrt{2}$mmだから、 sk1$=k_{①}=\frac{EA}{\ell}$とすると、 sk5$=k_{⑤}=\frac{EA}{\sqrt{2}\ell}=k_{①}/\sqrt{2}=$sk1/sqrt(2) とすればよい。

! 要素5の剛性マトリクス
sk5=sk1/sqrt(2.d0)
s(5,2,2)= sk5; s(5,2,4)=-sk5
s(5,4,2)=-sk5; s(5,4,4)= sk5
各要素の節点番号

なるべく一筆書きできるように、 なるべく要素番号と節点番号が左から右に増えるように、 図のように、初期状態の要素を$z$軸上に並べる。 このトラスは、一筆書きで書けるので、1直線上に並ぶが、 複数の直線になっても構わない。

各要素ごとに、初期状態での左節点番号(idari)と右節点番号(migi)を与える。

!各要素の左節点番号、右節点番号
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
各要素の回転角

次に各要素の回転角を、初期状態での左節点を中心に反時計回りの回転角で与える。 $yz$平面の回転角は、$x$軸右ねじ回転なので、反時計回りになる。

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

次は境界条件である。 回転角の下の以下のように書いてあるところを見つける。

!境界条件
xv(3)=0.d0; xw(3)=0.d0
xv(4)=0.d0!; xw(4)=0.d0

今回は、節点1で$v_{1}=w_{1}=0\;$, 節点2で、$v_{2}=0$だから、

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

今回は、節点3の$z$方向に、500Nを与えるので、

!載荷条件
fz(3)=500.d0 !N

とする。ちなみに、$y$方向に与える場合は、fy()を用いる。

実行

プログラムの修正が終わったら、 修正したプログラムは、例えば"torasu2sj_gotou2.f90"みたいに名前を変えて保存する。 プログラムをコンパイルして実行すると、 画像のような画面が表示される。

単位荷重法では、$v_{3}=0.5$mm, $w_{3}=(1+2\sqrt{2})\times 0.5 \simeq 1.91421356$mmと求まっていたが、 節点番号3のところを見ると、ほぼ合っている。

 節点番号:           3 v=  0.49999999999999989      w=   1.9142135623730945     

ちなみに、節点番号4のところで、

 節点番号:           4 v=   3.5162371365825888E-016 w=   1.9142135623730947     

$v_{4}$に値が入っているが、"E-016"の部分は、$\times 10^{-16}$という 意味なので、$v_{4}$はほぼ 0 ということだ。

先頭目次

プログラム課題5:

以下の5要素トラス(自分の学籍番号のもの)の 載荷点の節点変位($v, w$)を単位荷重法で求めよ。 長さ$\ell$の水平、鉛直の4本の部材が正方形を作り、 長さ$\sqrt{2}\ell$の斜材が1本入っている。
$E=10$GPa$=10000$N/mm$^{2}$
$A=(10$mm$)^{2}=100$mm$^{2}$
$\ell=1$m$=1000$mm
$P=500$N
とする。 torasu2sj.f90で解いて、 手計算の答えとどれくらい合うか考察せよ。 WebClassの「プログラム課題5」から、 単位荷重法の手計算を撮影した画像と プログラム本体(例えば torasu2sj_gotou2.f90 とか)と、 出力ファイル(例えば kekka2.out とか)を それぞれアップロードせよ。 どうしてもエラーが取れないなどの場合は、 kekka2.outの代わりに、エラー画面等 を写真で撮影したものをアップロードする。 また、最後のコメント欄に、手計算とプログラムの計算とを比較した考察や、 うまく実行できなかった場合は、 どういう症状で実行できないのかの説明を書き込むこと。 要素番号は、初期状態で左から①②③④⑤の順番に並べれば、 一筆書きできるようにしているつもりだが、 一筆書きできないように並べたとしても、適切に節点番号を指定すれば、 解けるはずである。

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

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

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

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

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

torasu2.f90プログラムソース


implicit real*8(a-h, o-z)
dimension f(18),& !節点力ベクトル
& d(18),& !節点変位ベクトル
& x(18),& !境界条件(拘束節点に0, その他に1が入る)
& ss(18,18),& !全体剛性行列
& xv(9),xw(9),& !節点ごとの変位境界条件
& fy(9),fz(9),& !節点ごとの荷重条件
& idari(9),migi(9),& !各要素の左節点番号、右節点番号
& s(9,4,4),& !要素剛性行列
& th(9),& !各要素の回転角
& t(4,4),& !座標変換行列(1要素ごとに計算)
& s44(4,4),& !要素剛性行列(1要素ごとに移し替える入れ物)
& ts(4,4),& !TK(1要素ごとに計算)
& tst(4,4) !TKT(1要素ごとに計算)
dimension sn(4) !節点力の出力のために追加

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

! 初期化
do n=1,nyou
do i=1,4
do j=1,4
s(n,i,j)=0.d0
end do
end do
end do

do i=1,18
d(i)=0.d0 !; f(i)=0.d0; x(i)=0.d0
do j=1,18
ss(i,j)=0.d0
end do
end do

do i=1,9
xv(i)=1.d0; xw(i)=1.d0 !境界条件は1で初期化
fy(i)=0.d0; fz(i)=0.d0
end do

E=10.d3!MPa
ell=1000.d0!mm
a=10.d0**2!mm^2
sk1=E*a/ell
! 要素1の剛性マトリクス
s(1,2,2)= sk1; s(1,2,4)=-sk1
s(1,4,2)=-sk1; s(1,4,4)= sk1
! 要素2の剛性マトリクス
sk2=sk1
s(2,2,2)= sk2; s(2,2,4)=-sk2
s(2,4,2)=-sk2; s(2,4,4)= sk2
sk3=sk1
! 要素3の剛性マトリクス
s(3,2,2)= sk3; s(3,2,4)=-sk3
s(3,4,2)=-sk3; s(3,4,4)= sk3
! 要素4の剛性マトリクス
sk4=sk1
s(4,2,2)= sk4; s(4,2,4)=-sk4
s(4,4,2)=-sk4; s(4,4,4)= sk4
! 要素5の剛性マトリクス
sk5=sk1/sqrt(2.d0)
s(5,2,2)= sk5; s(5,2,4)=-sk5
s(5,4,2)=-sk5; s(5,4,4)= sk5

!各要素の左節点番号、右節点番号
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(3)=0.d0; xw(3)=0.d0
xv(4)=0.d0!; xw(4)=0.d0

do i=1,9
x(2*i-1)=xv(i)
x(2*i)=xw(i)
end do
!
!載荷条件
fz(2)=200.d0 !N 
!
do i=1,9
f(2*i-1)=fy(i)
f(2*i)=fz(i)
end do
!
!
!要素ごとのTKTの計算
do n=1,nyou
!
! 座標変換マトリクス
!
call zahyou(t,th(n))
!
do i=1,4
do j=1,4
s44(i,j)=s(n,i,j)
end do
end do
!
call mxtmx(t,s44,ts)
!
call mxmx(ts,t,tst)
!!
!!重ねあわせ
i1=2*idari(n)-1
i2=2*idari(n)
m1=2*migi(n)-1
m2=2*migi(n)
!
ss(i1,i1)=ss(i1,i1)+tst(1,1)
ss(i1,i2)=ss(i1,i2)+tst(1,2)
ss(i1,m1)=ss(i1,m1)+tst(1,3)
ss(i1,m2)=ss(i1,m2)+tst(1,4)
!         !
ss(i2,i1)=ss(i2,i1)+tst(2,1)
ss(i2,i2)=ss(i2,i2)+tst(2,2)
ss(i2,m1)=ss(i2,m1)+tst(2,3)
ss(i2,m2)=ss(i2,m2)+tst(2,4)
!         !
ss(m1,i1)=ss(m1,i1)+tst(3,1)
ss(m1,i2)=ss(m1,i2)+tst(3,2)
ss(m1,m1)=ss(m1,m1)+tst(3,3)
ss(m1,m2)=ss(m1,m2)+tst(3,4)
!         !
ss(m2,i1)=ss(m2,i1)+tst(4,1)
ss(m2,i2)=ss(m2,i2)+tst(4,2)
ss(m2,m1)=ss(m2,m1)+tst(4,3)
ss(m2,m2)=ss(m2,m2)+tst(4,4)
!
!
end do !n要素について
!!
!!
!do i=1,8
!print'(8f10.3)', (ss(i,j),j=1,8)
!end do
!print*
!!
! 境界条件を入れる
do i=1,nset*2
do j=1,nset*2
  ss(i,j)=x(i)*ss(i,j)
  ss(j,i)=x(i)*ss(j,i)
end do
end do
do i=1,nset*2
if(x(i)<1.d-3) then
    ss(i,i)=1.d0
end if
end do
!!
!!
!
!print*,'f=',(f(j),j=1,18)
!
!
call gausu(nset*2,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*2-1),'w=',d(n*2)
end do
!追加
   !要素ごとのTKTの計算
   do n=1,nyou
   !
   ! 座標変換マトリクス
   !
   call zahyou(t,th(n))
   !
   do i=1,4
   do j=1,4
   s44(i,j)=s(n,i,j)
   end do
   end do
   !
   call mxtmx(t,s44,ts)
   !
   call mxmx(ts,t,tst)
   !!
    sn=0.0
    do i=1,4
      do j=1,2
         sn(i)=sn(i)+tst(i,j)*d((idari(n)-1)*2+j)
      enddo
      do j=3,4
         sn(i)=sn(i)+tst(i,j)*d((migi(n)-1)*2+j-2)
      enddo
    enddo
    print'(A,I2,A,I2,A,f15.10,A,I2,A,f15.10)','部材番号: ',n,&
        &'  S(',idari(n),')=',sn(1),'  ,N(',idari(n),')=',sn(2)
    print'(A,A,I2,A,f15.10,A,I2,A,f15.10)','             ',&
        &'  S(',migi(n),')=',sn(3),'  ,N(',migi(n),')=',sn(4)
   !
   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(4,4), b(4,4), c(4,4)
!
!
do i=1,4
do j=1,4
c(i,j)=0.
 do k=1,4
  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(4,4)
do i=1,2
do j=1,2
 t(i+2,j)=0.d0
 t(i,j+2)=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)=t(1,1)
t(3,4)=t(1,2)
t(4,3)=t(2,1)
t(4,4)=t(2,2)
return
end
!
subroutine mxtmx(a,b,c)
implicit real*8(a-h, o-z)
dimension a(4,4), b(4,4), c(4,4)
!
do i=1,4
do j=1,4
c(i,j)=0.
 do k=1,4
  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(18,18),x(18),b(18)
!
!ガウスの消去法の参考としたのは、
!名取亮「すうがくぶっくす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
!
!