トラスの解析には方向余弦という言葉が出てきます。もう一度定義を復習しておきま しょう。 部材の方向を表す単位長さのベクトルを方向余弦と言います。その式は(1a)のように 表すことができます。 部材の長さをℓとして具体的に成分を表記すれば,2次元の場合は(1b)のように,また, 3次元の場合は(1c)のようになります。
2
次に骨組み構造の用語である,「静定構造」,「不静定構造」という言葉について説明し ます。 「静定構造」とは部材が荷重ベクトルをどのように分担しているかベクトルの分解により 一意に決定できる構造です。 例えば上の図を見てください。節点2に赤矢印で示された外力が作用しているものと仮 定します。 節点2を共有する部材は二つしかありません。この外力を各部材の軸方向に分解する ことが一意にできます。 一方,下の図を見てください。節点2を共有する部材が3つあります。 この場合はベクトルの分解は一意に決定することができません。このような構造を不静 定構造といいます。 不静定構造は力を支えるのに余計な部材があると考えられます。 余分に部材があるのは無駄でしょうか。
さて前回の講義では2部材のトラス構造の変位を求める方法を学びました。 今回はさらに部材数が多い一般的なトラス構造の荷重と変位の関係を導きます。 今,図のように4部材からなるトラス構造を考えます。 部材および節点に通しの番号をつけることにします。要素の番号を「要素番号」、節点 の番号と「節点番号」と呼びます。 右上の表は各要素がどの節点から構成されているかを示していますが,この関係を要 素コネクティビティと呼びます。 ここで,各要素の要素剛性行列を図中の下側の表のように表記することにします. 上添え字の〇つきの数字は要素番号,下添字ij は節点i と節点j に関係する項である ことを意味します。
4
系のポテンシャルエネルギーは (3要素分のひずみエネルギー) - (4節点に作用する 外力のなす仕事) です。これを表したのが(2a)式です。 これをひとつの行列にまとめましょう。前回の講義内容に従ってそれぞれの行列を拡 張します。すると(2b)を経て(2c)を得ることができます.[K]は全体剛性行列と呼ばれま す。 この式に対してポテンシャルエネルギー最小の原理を適用します。 ポテンシャルエネルギー最小の原理は「変位はポテンシャルエネルギーが最小となる 値を取る」というものでしたね。 これよりポテンシャルエネルギーを節点変位で偏微分すると0になるという(3)式を得ま す。 (3)は拘束条件が入っていないためこのままでは解は不定となります。拘束条件(4)を (3)に代入することにより初めて解を求めることができます。
平面トラスの解析プログラムの流れを示します. 初めに要素剛性行列を計算しては全体剛性行列にこれを加えるという作業を行います 。続いて荷重ベクトルを作成します。 そして,拘束条件を代入しながら連立方程式を解いて節点の変位を求めます。 最後に求めた節点変位より部材の伸び,ひずみを求めヤング率を乗じることにより応 力を計算します。 全ての部材の応力が材料の降伏応力よりも十分低ければ解析した構造は安全である と判断します。
6
MATLABでは連立方程式はどう解くのかを意識せずに組み込み関数によって簡単に解 を得ることができます. しかしながら連列方程式の解法は数値計算としては最も基本的でありエンジニアの常 識として知っておいた方が良いと思います. ここで連立方程式の数値解法について簡単に説明しておきます。 連立方程式の数値解法は大別して「反復解法」と「直接解法」があります。 反復法は適当な近似解からスタートして反復計算によって解を改善していく方法で,使 用メモリが少ないことが大きな長所です.短所としては条件が悪い場合は解が収束し ないことがあることが挙げられます. 一方,「直接解法」は行同士の足し算,引き算を巧みに行うことによって解を直接求め る方法です. メモリの使用量が大きいのですが解を求めることに失敗することが少ないという長所が あります。 通常は直接解法で解き,非常に係数行列が大きくメモリが足りないときには反復解法 を用いるといのが一般的でしょう. ここでは直接解法の代表的手法としてガウスの消去法を学びます.
ガウスの消去法では「前進消去」過程と「後退代入」過程のふたつの過程からなります . 前進消去過程では,ある行に適当な係数を乗じて他の行に足したり引いたりしても解 は変わらないという性質を利用して,図のように下の行になるにつれ階段状に未知数 の数が減るよう方程式を変形していきます。 後退代入過程では,まず一番下の行より解を求めます.一番下の行は未知数ひとつ の単純な一次方程式になっているので容易に解を求めることができます。 続いて下から2行目(未知数2)に移り,最初に求めた解を代入した後,残った未知数に ついて解を求めます. 以降,1行ずつ上に移動しながら,既知の解を代入し,未知数を求める作業を続けます .
8
例題で実際に前進消去を行ってみましょう. まず,1行目のxを利用して2行目,3行目のxを消去します。このためにまず1行目のxの 係数を1にします。この式を(1)’とします. 続いて2行目よりxを消去するために(1)’式を3倍したものを2行目よりひきます。これを (2)’式とします. さらに3行目からもxを消去するために1)’式を5倍したものを3行目より引きます。.以上 でxの消去は終わりました。 続いて2行目を利用して3行目のyを消去します。すなわち(2)”式のようにyの係数を1と した上でこれを3倍して(3)’行から引きます。 行列は下三角形となり前進消去過程は終了です。
続いて後退代入過程を開始します。 まず一番下の行を用いてzを求めます.zは6となります。 続いて下から2行目に移ります。 下から2行目の式をyについて解くと y=3+zになりま すので,このzに既に解き終わった6の値を代入します。 こうしてy=9を得ます。 最後に一番上の式に移り,xについて解いた式にyとzの値を代入するとx=33を得ま す。
10
それではこのページで与えられる連立方程式に対してガウスの消去法を適用して前進 消去してください。
ガウスの消去法アルゴリズムをここで整理してみましょう. すなわち前進消去では,k=1,2の順に ①対角項の係数を1にするため自分の行をakk で割る. ②自分より下の行についてk列を消去する,という手続きを繰り返します。 次に一番下の行から1つずつ,さかのぼりながら解を求めていきます。. これを適当な言語でプログラミングすれば任意の未知数の連立方程式を解くことがで きます.
15
有限要素法で作成される全体剛性行列には1)対角項からある幅の部分だけ天の川 のように非ゼロ項が入る(バンド特性)、2)対称性(A(I,j)=A(j,i)) がある,という特徴があ ります. 商業用の有限要素法のプログラムでは,この性質を利用して行列のバンド幅内の上半 分(または下半分)のみを記憶することによってメモリを節約しています。 また,バンドの外に対する計算も結果がゼロになるのは明白なので,とばすようにプロ グラムとして計算時間を節約しています。 以上のような特質のために有限要素法は未知数の割には連立方程式を解くのが早く, これがライバルの数値解法に対する強みとなっています.
ここからは計算力学に必要なmatlabの基礎を学びます。初めにインストールした matlabを立ち上げてください。Tool BOXなどオプションを選択した場合は多少プログラ ムの立ち上がりが遅くなります。
MATLABを立ち上げると図のような画面が出てきます。このウィンドウをコマンドウィ ンドウと呼びます。コマンドウィンドウに計算式を入力するとMATLABは実行結果を表 示します。数行程度の簡単な計算ならコマンドウィンドウに計算式を入力して計算すれ ば良いでしょう。 しかし,ある程度行数が多いプログラムになるとエディタを立ち上げてプログラムを入 力して実行する方が便利です。エディター内ではワープロのように自由自在にプログラ ムを編集することができます。 また,作成したプログラムを保存し繰り返し実行することができます。 この講義ではエディタを使ってプログラムを入力し実行することにします。 それではエディタを立ち上げて見ましょう。画面左上にある新規スクリプトというアイコ ンをクリックしてください。新しいウィンドウが開きエディタが表示されます。
エディタが立ち上がったら1行目に1+2と入力してみましょう。続いてエディタの上の 方にある緑色の三角形の「実行」ボタンをクリックしてください。 するとプログラムを保存するダイアログが開かれます。 プログラムを保存する先はデ フォルトはドキュメントの下のMATLABというフォルダです。 別なフォルダでも構いません。ただし,MATLABは日本語を得意としません。日本語 フォルダは避けてください。 ファイル名も任意でかまいませんが英文としましょう。また,ファイルの最後には必 ず,.m(どっとエム)としてください。このドット以下の文字をファイル識別子といい, ファイルの種類を分類するのに重要です。 ファイルを保存するとプログラムが実行されてコマンドウィンドウに結果の ‘3’ が 表示されます。エディタ画面には何も表示されないので注意してください。
19
最初にスカラの四則演算について説明します。四則演算とべき乗はC言語と同じです。 ただし,文末のセミコロンの意味がC言語とは異なります。 MATLABでは文末にセミコロンをつけると計算結果を子マントウィンドウに表示しませ ん。 一方,セミコロンを付けない場合は計算結果がコマンドウィンドウに表示されます。 また,C言語では変数の型宣言が必要ですが,MATLABは変数の型は自動で判定さ れ,意識する必要がありません。
それでは先ほど1行入力したエディタに上書きでかまいませんので,画面のようなプロ グラムを入力して実行ボタンを押してください。 上書きの場合は実行ボタンを押してもファイル保存のダイアログは立ち上がりません。 先ほどのファイル名で自動的に保存されます。 プログラムは4行からなりますが,上の2行が文末にセミコロンがあるのに対して下の2 行はセミコロンがありません。 MATLABはセミコロンのないこの2行の計算結果をコマンドウィンドウに表示します。
21
次にベクトルの定義方法を説明します。横ベクトルは大かっこでくくって定義します。要 素の区切りは空白またはコンマとなります。 また,その要素は小かっこでくくって引用します。たとえばaの一番目の項は小かっこを つかってa(1)のように引用します。 MATLABでは横ベクトルと縦ベクトルを厳密に区別します。縦ベクトルは大かっこの中 にセミコロンをいれて表示します。 また,横ベクトルを転置することによって縦ベクトルを作ることができます。MATLABで はダッシュが転置を表します。
23
MATLABでは等差数列を簡単に定義することができます。 等差数列は初期値,コロン,増分値,コロン,終わりの値,で定義します。 例えば,a イコール 1 コロン 2 コロン 9 と書くと,1から9まで2ずつ増える数列が横ベクトルとして定義されます。 ただし増分が特に1の場合は初期値と終わりの値だけで定ぐすることができます。 b イコール 1コロン 5 と書くと1から5までの整数ベクトルが作られます。 MATLAB特有の関数にsetdiffというものがあります。これは一番目のベクトルと2番目 のベクトルを比べ,2番目にある要素は削除します。トラスの解析では全体の自由度番 号から 拘束されている自由度番号を引いて拘束されていない自由度番号列を作ります。この ときにこのsetdiffを使います。
ベクトルの和や差は我々が数式で扱うのと同じように変数名同士をプラス記号または マイナス記号で結ぶだけで計算することができます。 ただし,ベクトルの乗算,つまり掛け算には注意が必要です。MATLABは縦ベクトルと横 ベクトルを厳密に区別します。 例えば横ベクトル同士の乗算 (a)式はエラーとなります。 一方,(b)式のように横ベクトルとその転置である縦ベクトルの乗算は内積となります。 組み込み関数 dotを用いても同じ結果を得ることができます。 さらに(c)では縦ベクトルと横ベクトルの乗算になります。この場合は演算結果は行列 となります。
25
次に行列について説明します。行列もベクトルと同じように大かっこを使って定義しま す。行の区切りはセミコロンで表されます。 したがって(a)のように記述すると(b)の行列が定義されます。 また,MATLABでは小行列を使って大行列を定義することができます。 例えば(b)で定義した行列Aについて(c)式のように記述すれば,行列Aの横に行列Aの 転置行列が置かれた(d)式のような行列が定義されます。 また,(e)式のように書けば行列Aの下に行列Aが定義された(f)式で表される行列が定 義されます。
MATLABでは組み込み関数を使って有用な行列を定義することができます。 ゼロスはゼロ行列を定義する関数です。ゼロスのnと書くとn行n列の正方ゼロ行列が 定義されます。 また,アイは単位行列,すなわち体格項のみが1の行列を定義します。 なお,正方行列ではなくn行m列の行列のゼロ行列や単位行列を作りたい場合は(c)式 や(d)式のように引数をnとmのふたつにします。
27
次に小行列を大行列に加える方法について説明します。トラスの解析では要素剛性行 列を全体剛性行列に加えます。 全体剛性行列の中の要素に関係する自由度の行と列に要素剛性を加えることが必要 です。 MATLABはこのような作業が得意です。簡単な例で説明しましょう。 まず3行のプログラムを見てください。 1行目で行列bとして3行3列のゼロ行列を定義しています。 次に,2行目で行列aとして2行2列の単位行列を作成しています。 そして3行目では行列bの1,2行と1,2列に行列aを加えるという演算をしています。 この演算結果が(a)から(c )式で表されています。 (c)では行列bの指定された位置に 行列aが加えられていますね。
最後に連立発定式をMATLABで解く方法を説明します。通常は連立方程式を解くプログ ラムはGaussの消去法などの関数を自作する必要があります。 しかしMATLABは’¥’マークを用いて簡単に連立方程式を解くことができます。 たとえば(a)式のような連立方程式を解くときのプログラムをその下3行に示します。 初めに係数行列をAとして定義します。次に右辺ベクトルをbとして定義します。 連立方程式の解は¥マークを使うことにより計算できるのです。
29
それでは本日の課題2,3に取り組んでください。答えはCoursePowerのレポートに記入 してください。