マトリクス構造解析オンライン授業用テキスト
第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$
まず、載荷点の$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}$を求めるには、
節点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}$$
各要素を節点の近傍で切断して、部材力を書き込んでみると、
要素にゼロでない部材力が作用しているのは、⑤要素と②要素のみ。
⑤要素は$\sqrt{2}P=500\sqrt{2}$Nの引張、
②要素は$P=500$Nの圧縮。
ということは、これらを$yz$方向の節点力で表してやると、
⑤要素の節点1に、$S_{1}^{⑤}=500$N, $N_{1}^{⑤}=-500$N
⑤要素の節点3に、$S_{3}^{⑤}=-500$N, $N_{3}^{⑤}=500$N
②要素の節点2に、$S_{2}^{②}=-500$N, $N_{2}^{②}=0$
②要素の節点3に、$S_{3}^{②}=500$N, $N_{3}^{②}=0$
では、 前回のプログラム(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 ということだ。
各要素を節点の近傍で切断して、部材力を書き込んでみると、
要素にゼロでない部材力が作用しているのは、⑤要素と②要素のみ。
⑤要素は$\sqrt{2}P=500\sqrt{2}$Nの引張、
②要素は$P=500$Nの圧縮。
ということは、これらを$yz$方向の節点力で表してやると、
②要素の節点2に、$S_{2}^{②}=-500$N, $N_{2}^{②}=0$
②要素の節点3に、$S_{2}^{②}=500$N, $N_{3}^{②}=0$
⑤要素の節点1に、$S_{1}^{⑤}=500$N, $N_{1}^{⑤}=-500$N
⑤要素の節点3に、$S_{3}^{⑤}=-500$N, $N_{3}^{⑤}=500$N
以下の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の人。
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
!
!