どうも食卓塩です。
今回は日刊工業新聞社が2014年に出版した「有限要素法のつくり方！」の批評を行います。

書籍の紹介(Amazonから)
有限要素法(FEM)の基礎を習得した方を対象に10年間以上続いている、NPO法人CAE懇話会で行われている、FEMプログラムを作成する講義および実習「解析塾 FEMプログラミング」を1冊の本で紹介。FEMプログラミングの手順やノウハウを知ることで、自分でFEMプログラムを作ることを学ぶ本。
著者略歴 (「BOOK著者紹介情報」より)
石川/博幸
工学博士。名古屋大学大学院工学研究科修了。S39年生まれ、愛知県出身、メーカーでCAE技術者として勤務

青木/伸輔
技術士(情報工学部門、総合技術監理部門)、計算力学技術者(固体1級)。豊橋技術学大学大学院工学研究科建設工学専攻修士課程修了。S53年生まれ、高知県出身、株式会社中電シーティアイ勤務

日比/学
慶應義塾大学理工学部機械工学科卒業。S38年生まれ、岐阜県出身、ブラザー工業株式会社勤務。NPO法人CAE懇話会の中部CAE懇話会発足より幹事を務める(本データはこの書籍が刊行された当時に掲載されていたものです)

監修者　石川　覚志
工学博士　京都大学大学院工学研究科修了
S36年生まれ、大阪市出身、株式会社IDAJ勤務
非線形有限要素法プログラムのカスタマーサポートを経て、非線形解析の受託解析業務に従事。2001年から、NPO法人CAE懇話会の解析塾において非線形構造解析コースを担当。延べ350名以上の受講生を指導する。2006年より、日本ゴム協会・ゴムの力学研究分科会の書記を担当。

序文：小寺秀俊京大理事・副学長 CAE懇話会副理事長

　はじめに、言い方が雑というか礼儀など放り出した書き方をしてるのは書籍の内容による。敬称は省略させていただく。つける気すら起きないが。
　内容は懇話会とやらのテキストを編集したものとあるが、ターゲットは誰？もう知ってる人向け？それとも初心者？そういう内容の書籍である。
　これが初心者向けなら、まぁ、最初の方は理解できなくもない気がするが、プログラミング編を見るとこの期待はあっさり裏切られ、知ってる人向けであればこいつ等ふざけてるのか、な内容になる。この本はそういう代物である。
　つまりお勧めできない。Amazonの書評にある通りVBAの学習用と割り切った方がいい内容である。
　このようなものを買うのだったら、例えば三好敏郎著「有限要素法入門」を購入することを薦める。

　章立てを以下に示す。

第1章有限要素法概要
第2章シンプルな2次元プログラム
第3章高度な2次元プログラミング
第4章前処理・後処理

　まず第一章から始める。
　そもそもこの部分からして正確さを欠いている。初学者だからと言って雑な説明をすればいいというものではない。この部分を大まかに引用すると次のようになる。

　三次元空間内の微小空間を要素と称する。最終的に三角形になるとすれば…

　はい？何がどうなると「最終的」に二次元空間の話になるんだ。それとも、もっと説明が難しいシェル要素の話でもする気なのか。なら初めから二次元の平面応力状態から始めればいいものを何でこんな説明をしているのか理解できない。一般的な機械構造物の説明をすれば十分だったはずである。
　この部分の説明については、要するに「機械構造物に荷重がかかれば変形をし、応力やひずみが生じます。」これで十分である。この後取り敢えず構成則について書き、エネルギー理論の触りから仮想仕事の式の説明をざっくばらんにすればいいだろう。要するに詳しい式の導出は避けるが、とりあえずこれを使って計算をするんだよ、ということである。
　続いて応力をベクトル量と書くのは正確ではない。吹き出しではなく本文で、ベクトル量を使った方が説明が簡単になるので、ここではベクトル表記を使いますと書いた方がよかっただろう。
　次に構成則について書いてるが、本文では「両者(著者註:応力とひずみのこと)はあるマトリックスDによって関連付けられる」とあるだけである。この後に

\sigma=D\epsilon

の式が出てくる。これだけですか？「このような関係を構成則と言い、金属材料では専らフックの法則が使われている。」位は書けよ。説明が足らないだろうが。
　この後改行無しにひずみエネルギーの説明が始まる。内容が変わるのだから改行を入れておくべきである。
　この章の最大の問題は定式化が怪しいことである。式の説明は正しいのだが、肝心の式がおかしい。

E_i=\frac{1}{2} \int_V \left\{\sigma\right\}^T\left\{\epsilon\right\}dV

　応力とひずみの位置が逆である。内積なので交換則が成り立つから逆にしても間違いではないが、この後の定式化でボロが出る。
　初心者向けならば、ベクトルは基本的に縦に書くことを本文中で書くべきだろう。因みに応力をベクトル表記すると二次元の場合次のようになる。

\mathbf\sigma=\left\{\matrix{\sigma_x\\\sigma_y\\\sigma_{xy}}\right\}

　横書きする場合は記号の肩にTをつけて次のように書く。

\mathbf\epsilon^T=\left\{\matrix{\epsilon_x&\epsilon_y&\epsilon_{xy}}\right\}^T

話が逸れたが、初学者にはこの程度のことは教えないといけない。中学数学ではベクトルは基本的に横書き、つまり行ベクトルを用いているが、基本的に縦書き、つまり列ベクトルを用いる。
　この後、式の変形を行ってどこかで見た式が登場する。

E_i=\frac{1}{2}\left\{u\right\}^T\left( \int_V [B]^T[D]^T[B]dV\right)\left\{u\right\}

　これを見て違和感を感じない人は、これ以外の書籍を穴が開くほど見るといいだろう。バカな話だが、最初の式で応力ベクトル${\sigma}$に転置を掛けた結果、Dマトリックスにも転置がかかっているのである。細かい事とは言ってはいけない。定式化というは大学者たちの研鑽の結果による、数十年の結果であって気紛れでできたものではない。
　この本では何故か仮想仕事の原理は出てこない。初学者に教えても理解できないと言いたいのかもしれないが、私の意見は異なる。初学者だからこそ、その点はきちんと教えなければならない。勿論定式化に至るまでの導出過程まで教えるのは骨が折れるため、結果の式だけを出し、なぜこれが使われているのかを説明するべきなのである。というかそれが本来の仕事ではないのか。

　仮想仕事の原理の式は次の式で表される。

\int_V \delta\epsilon^T\sigma dV=\int_V\delta u^TfdV \cdots (1)

　ここでfは外荷重一般を指すものとする。この式を大雑把に説明すると、微小な変位がかかった場合、それに見合う微小なひずみが生じるというものである。つまり系は釣り合う。あくまでも大雑把な説明なので、正確なところは専門書を見て欲しい。初めに勉強する人にとって大変だろうが、なぜ、弾性力学の解析に仮想原理の原理を使用するかというと、応力とひずみを関係づける構成則に依存しない、つまり弾性のみならず弾塑性問題の解析やゴムといった超弾性問題にも使えることと、ここから定式化がなされる変位法の実装が容易だからである。

　ここから、次に示す構成則とひずみと変位の関係式を代入すると、よく知られた変位法の式が出来上がる。
とはいっても、解析するべき物体についてはnelem個の要素に分割されたものであるから、その点を考慮しなければならない。なお、「解析するべき物体」と言ってるが、実際は「解析領域」とか「計算領域」ということがある。
　この結果、左辺は要素ごとの式を式の上では合計したものになる。

\epsilon=Bu \\
\sigma=D\epsilon \\
\epsilon^T=(Bu)^T=u^TB^T \\

\sum^{nelem}_{i=1}\int_V \delta u_i^TB^T_iD_iB_iu_idV=\int_V\delta u^TfdV \cdots (2)

式に出てくる$u_i$は最終的に$u$にまとめられるから、次のように変形される。また、変位に関する項は領域積分に関係しないので、これも外に出せる。
　これらを考慮すると次のように変形される。

\delta u^T\left(\sum^{nelem}_{i=1}\int_V B^T_iD_iB_idV\right)u=\delta u^T\int_VfdV \cdots (3)

このうち両辺に出てくる$\delta u$は式から消すことが出来るので、最終的に以下のようになる。

\left(\sum^{nelem}_{i=1}\int_V B_i^TD_iB_idV\right)u=\int_V fdV \cdots (4)

　次にこの式の左辺を次のように置き換える。

Ke_i=\int_V B_i^TD_iB_idV

　そうすると(4)式は次のようになる。

\left(\sum^{nelem}_{i=1}Ke_i\right)u=\int_V fdV \cdots (5)

更に左辺を次のように置き換える。

K=\sum^{nelem}_{i=1}Ke_i

そうすると(5)式はさらに次のようになる。

Ku=\int_V fdV \cdots (6)

ついでに右辺も単純に

F=\int_V fdV

のように置き換えると(6)式は最終的に以下のようになる。

Ku=F \cdots (7)

　要するにフックの法則そのものである。ただ、ばねの場合と異なるのは、Kがマトリックスになり、u・Fがベクトルになっていることである。つまり連立一次方程式が出来上がるわけだが、式の形としては解析領域全体を一つのばねと見立てて変位や荷重の釣り合いを求めるのである。

　だいぶ脱線したが、説明するならば読者や受講者が解ろうが解るまいがきちんと原理となる式や理論を提示し、なぜそれを使うのかを説明した方が余程いい。

　ところで、序文(正確には「はじめに」)は京大の副学長様が書いておられるのですが、本書をちゃんと読んだのでしょうか。自分なら書き直しを命じるところですが。