第2章の残りの部分の批評に入る。
2.7節
ここの順番は反力→ひずみ→応力の順となっているが、ひずみ→応力→反力の順とすべきである。実際の構造解析では応力やひずみを求めることが多いからである。
2.7.1の反力の計算だが、全体剛性マトリックス$\mathbf K$と、これも全体変位ベクトル$\mathbf U$から、$\mathbf K\mathbf U=\mathbf F$を用いて全体荷重ベクトル$\mathbf F$を求めているが、間違いではないにしても、一般的な方法ではない。コラムにある方法をそっくりそのまま説明したほうがよかった。
そのコラムにある方法については後述する。
2.7.2はひずみの計算の説明だが、一行目にあるように端的にいって$\mathbf \epsilon=\mathbf B\mathbf u$から計算しますの一言で終わる。式だけの行がないから、よく見ないと見落とすことになる。次に節点変位ベクトルが必要だが、コネクティビティを用いて、全体変位ベクトル$\mathbf U$から引用すればいいと書けば終わる。というか、基本的にひずみの計算はこれだけであるが、この時点でコネクティビティについて理解していない前提で書くというのが抑々おかしい。
式を書く際、文中に埋める方法と独立した行にする方法があるが、この本ではそういう方法をとったり取らなかったりするので読み辛いのである。グダグダと文章だけで書けばいいというものではない。
プログラム例もいいものとは言えない。全くの初心者向けか「さいごのゆうげんようそほうぷろぐらむ」を作成する連中向けにはいいかもしれないが。
原著のこの部分は次のようになっている。初心者向けを謳っているからだろうが、ひずみと応力の計算ルーチンを別にしているが、2次元要素なら一体とすることも可能である。この辺りは要素の定式化にかかわるので、応力とひずみの計算部分を分離するのが悪いとは言えない。
' ----------------------
' 要素ひずみプロシージャ
' ----------------------
Sub make_strain_element()
Dim Ue(DOF_TRIA3) As Double ' 要素内節点変位ベクトル
Dim e As Integer ' 要素番号カウンター
Dim n As Integer ' 節点番号カウンター
Dim r As Integer, c As Integer ' マトリックスの行と列のカウンター
For e = 1 To ELEMENTS ' 要素ごとに要素ひずみを計算
For n = 1 To NODES_TRIA3 ' 要素内節点変位を計算
Ue(n * 2 - 1) = U(connectivity(e, n) * 2 - 1) ' x成分
Ue(n * 2) = U(connectivity(e, n) * 2) ' y成分
Next n
For r = 1 To COMPONENTS ' 要素のひずみを計算 (B×Ue)
strain_element(e, r) = 0#
For c = 1 To DOF_TRIA3
strain_element(e, r) = strain_element(e, r) + B(e, r, c) * Ue(c)
Next c
Next r
Next e
End Sub
初心者向け説明用ならこれでいいのだろうが、そこで止まっているから、これを読んだ人はそのあと苦労することになるのが目に見えている。この本ではループカウンタに$r,c$を使っている。変数$e$も同じで、意味を持たせるなら$ielem$くらいにはして欲しかった。変に手を抜いてどうする。こんなことするをくらいなら、ループカウンタは伝統的な$i,j,k$で十分である。
コメントも不要である。冒頭のコメントはルーチンの説明になっていない。名前を見ればわかるので、コメント自体不要である。精々「$B\times Ue$から要素のひずみを求める」で十分だろう。
何故なら元の式を熟知していればこの部分は簡単に読み下せるからである。特にループカウンターや$e$のループ初めのコメントは不要である。こんなコメントを読むくらいなら文学作品でも読んだ方がいいだろう。教育目的というならコメントを読ませるのではなく、プログラムを読ませる方に持って行った方がいいだろう。
サブルーチン全体のコメントはそれが何をするかを端的に書くべきである。
この部分を直すと直接$\mathbf \epsilon=\mathbf B\mathbf u$を用いて計算するのはいい方法とは言えない。何度でも繰り返すが、計算機に適した方法に直せばプログラムは容易になる。この部分に$\mathbf A$マトリックスを用いると、ひずみの式は次のように変換される。
初めに$\mathbf A$マトリックスを次のようにおく。
\mathbf A=\left[ \matrix{
\frac {\partial N_1} {\partial r_1} & \frac {\partial N_2} {\partial r_1} & \frac {\partial N_3} {\partial r_1} \\
\frac {\partial N_1} {\partial r_2} & \frac {\partial N_2} {\partial r_2} & \frac {\partial N_3} {\partial r_2}
} \right]=\left[ \matrix{
\mathbf A_1 \\ \mathbf A_2
} \right]
これを求めておけば実質的にひずみ変位マトリックスの計算が終わるが、これに符号行列$\mathbf \Lambda$
\mathbf \Lambda=\left[ \matrix{
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 1 & 1 & 0
} \right]
を用いて次のような式を立てる。
\mathbf B=\mathbf \Lambda\left[\matrix{
\left[\matrix{
\mathbf A_1 \\ \mathbf A_2
} \right] & 0 \\
0 & \left[\matrix{
\mathbf A_1 \\ \mathbf A_2
} \right]
} \right]
これを展開すると次のようになる。
\mathbf B=\left[\matrix{
\mathbf A_1 & 0 \\
0 & \mathbf A_2 \\
\mathbf A_2 & \mathbf A_1
} \right]
どこかで見たような式になるが、これが実はあちこちの書籍で見かける$\mathbf B$マトリックスの具体形に他ならない。
これに要素の変位ベクトル$\mathbf U_e$をかけた式がひずみになる。式中の$u,v$は要素変位の$x,y$方向の変位ベクトルになる。
\mathbf \epsilon = \mathbf B\mathbf u=
\left[\matrix{
\mathbf A_1 & 0 \\
0 & \mathbf A_2 \\
\mathbf A_2 & \mathbf A_1
} \right]
\left[\matrix{
\mathbf u \\ \mathbf v
} \right]=
\left[\matrix{
\mathbf A_1\mathbf u \\
\mathbf A_2\mathbf v \\
\mathbf A_1\mathbf v+\mathbf A_2\mathbf u
} \right]
最後の項がベクトルの内積になっていることに注意。
この部分を見ると初心者向けでないように思える方がいるかもしれないが、重要なことは計算しやすい形に式を変形することだけである。そうすればメモリも速度も節約できる。
これをプログラムにするとこのようになる。変位$dx,dy$を求め、それをAマトリックスの各要素にかけているだけである。上の式そのままに計算していることがわかるだろう。
For j = 1 To 3
dx = disp((conn(i, j) - 1) * 2 + 1)
dy = disp((conn(i, j) - 1) * 2 + 2)
strain(i, 1) = strain(i, 1) + A_mat(i, 1, j) * dx
strain(i, 2) = strain(i, 2) + A_mat(i, 2, j) * dy
strain(i, 3) = strain(i, 3) + A_mat(i, 2, j) * dx + A_mat(i, 1, j) * dy
Next j
2.7.3で応力の計算が出てくる。説明は本文の通りだが、
\mathbf \sigma = \mathbf D \mathbf \epsilon
の式位は載せて欲しかった。
この部分のプログラムは次のようになっているが、コメントとや変数については先にも指摘した通りで、それ以外は誰が書いてもこのような形に落ち着く。
' --------------------
' 要素応力プロシージャ
' --------------------
Sub make_stress_element()
Dim e As Integer ' 要素番号カウンター
Dim r As Integer, c As Integer ' マトリックスの行と列のカウンター
For e = 1 To ELEMENTS ' 要素ごとに要素応力を計算 (D×strain_element)
For r = 1 To COMPONENTS
stress_element(e, r) = 0#
For c = 1 To COMPONENTS
stress_element(e, r) = stress_element(e, r) _
+ D(r, c) * strain_element(e, c)
Next c
Next r
Next e
End Sub
なお、自作プログラムのこの部分は応力とひずみを一気に求めている。最後の方で、Z方向の応力とひずみを求めている。
Sub calc_result()
Dim i As Integer, j As Integer, k As Integer
Dim dx As Double, dy As Double
ReDim strain(ntelem, 4)
ReDim stress(ntelem, 4)
For i = 1 To ntelem
For j = 1 To 3
dx = disp((conn(i, j) - 1) * 2 + 1)
dy = disp((conn(i, j) - 1) * 2 + 2)
strain(i, 1) = strain(i, 1) + A_mat(i, 1, j) * dx
strain(i, 2) = strain(i, 2) + A_mat(i, 2, j) * dy
strain(i, 3) = strain(i, 3) + A_mat(i, 2, j) * dx + A_mat(i, 1, j) * dy
Next j
For j = 1 To 3
For k = 1 To 3
stress(i, k) = stress(i, k) + D_mat(k, j) * strain(i, j)
Next k
Next j
If prob = ProbPSN Then stress(i, 4) = poisson * (stress(i, 1) + stress(i, 2))
If prob = ProbPSS Then strain(i, 4) = -poisson / young * (stress(i, 1) + stress(i, 2))
Next i
End Sub
応力とひずみの成分が4個になっているが、平面問題では応力とひずみどちらかが4個になるからである。なお、軸対称問題では両方とも4個必要である。コラムでこの部分の説明をしているのだが、だったら本文にプログラム例諸共書けばいいと思うのは私だけじゃないだろう。
なお、VBAの組み込み手続き、Redimを使うと必要な大きさの配列を割り当てるが、同時に配列の内容を0にするので初期化の必要はなくなる。この書籍では配列は全部あらかじめ決められた数を指定しているが、実際のプログラムでは入力を読むまでデーターの個数は不定なので、プログラムの実行中で必要な容量を割り当てる必要が出てくる。このような技法を動的割り当てと言っているが、そんなに難しいものではないし、使えば必要なメモリだけ用意できるので無駄がなくなる。
このようなことを説明しなくてどうする。
最後に等価節点力を求める。本書では単に反力と言っているが、実際には要素ごとの節点力を等価節点力といい、反力を求める際にはこれらを自由度ごとに足し合わせて求めることになる。この場合、外力も求まることになる。
これもAマトリックスを用いれば容易に求めることができる。
三角形一次要素の要素剛性マトリックスは次の式で計算できた。
\mathbf Ke=\Delta t\mathbf B^T\mathbf D\mathbf B
等価節点力$f_e$は次のようになる。
\mathbf f_e=\mathbf K_e\mathbf u_e
$\mathbf \epsilon=\mathbf B\mathbf u_e,\mathbf \sigma=\mathbf D\mathbf \epsilon$から$\mathbf f_e$は次のようになる。
\mathbf f_e=\Delta t\mathbf B^T\mathbf \sigma
ここで、力の方向成分をベクトル化した$\mathbf f_x,\mathbf f_y$用いると、上の式は次のようになる。
\mathbf f_e=\left[\matrix{
\mathbf f_x \\ \mathbf f_y
}\right]
\mathbf \epsilon =
\mathbf B^T
\left[\matrix{
\epsilon_x \\ \epsilon_y \\ \epsilon_{xy}
}\right]
=\left[\matrix{
\mathbf A_1^T & 0 & \mathbf A_2^T \\
0 & \mathbf A_2^T & \mathbf A_1^T \\
}\right]
\left[\matrix{
\epsilon_x \\ \epsilon_y \\ \epsilon_{xy}
}\right]
となるから、等価節点力は次のようになる。$\mathbf A_1^T,\mathbf A_2^T$は縦ベクトルになることに注意。また、$\mathbf f_x,\mathbf f_y$も要素の荷重成分ごとのベクトルになる。
\mathbf f_e=\Delta \left[\matrix{
\mathbf f_x \\
\mathbf f_y
}\right] =
\left[\matrix{
\mathbf A_1^T \epsilon_x + \mathbf A_2^T \epsilon_{xy} \\
\mathbf A_2^T \epsilon_y + \mathbf A_1^T \epsilon_{xy}
}\right]
最後に自作プログラムでの例を示す。等価節点力は直接全体荷重ベクトルに足しこまれる。
Sub calc_eqiv_force()
Dim i As Integer, j As Integer, k As Integer, pos As Integer
ReDim eqiv_force(ntnode * 2)
For i = 1 To ntelem
For j = 1 To 3
pos = (conn(i, j) - 1) * 2
eqiv_force(pos + 1) = eqiv_force(pos + 1) + (A_mat(i, 1, j) * stress(i, 1) + A_mat(i, 2, j) * stress(i, 3)) * area(i) * thick
eqiv_force(pos + 2) = eqiv_force(pos + 2) + (A_mat(i, 2, j) * stress(i, 2) + A_mat(i, 1, j) * stress(i, 3)) * area(i) * thick
Next j
Next i
End Sub
この後、2章のプログラム全体が掲載しているが、プログラム自身はCAE懇話会のサイトからも引っ張ってくるのが可能なので、興味のある方は見るといいかも知れない。
この書籍がどこを向いているのか分からないと書いたが、教えるべき内容がプログラミング例を通した有限要素法入門であることを考えると、内容が不足していることは以前にも書いた。まず、基本的な原理について書くことや、式に沿った説明を行うことである。初心者相手に変分原理の説明までする必要はない(死屍累々になるのが目に見えている)が、「何故それを使うのか」の説明はあってしかるべきである。
次にこの分野はコンピューターサイエンスに片足程度は突っ込んでいるはずである。しかし、この章を通してみると冗長な部分が目立つ。初心者向けの説明だからという理由が聞こえてきそうだが、この書籍で行うことはここで指摘したように「以下に計算機向けに式を変形するか」という話である。応力を求める部分は単に式をそのまま書けばいいが、ひずみ変位マトリックスやひずみを求める部分は「こうすればもっとプログラムが簡単になりますよ」という例である。そうすればミスだって減るだろうし、トータルの開発効率も上がる筈である。そういうことが全く反映されていないのはなぜなのかと言いたくなる。
参考文献
林晴比古「BASICによるプログラミングスタイルブック」ソフトバンククリエイティブ (1985/12) BASICによるプログラム書法の本。
SUN FORTRANリファレンスマニュアル