有限要素法とPythonプログラム（三角形一次要素を用いた２次元応力解析プログラム）

（情報）初期のプログラムから以下の点を変更しています．

• 疎行列処理導入による高速化
• 関数名を小文字にした
• データ入力部を関数にした
• 計算用関数（データ入力・プリントアウト部を除く）の引数を全てスカラーとした
• 無意味な配列初期化の部分をカットした
• 配列初期化時の形状をタプルで指定するようにした（何故かこれまでリストで渡していた）
• メインルーチンも関数化した．これは下の記事を見てハット思い実行しました．これによりグローバル変数となっているものを関数の引数で呼び込んだりというおかしなことがなくなったと思います．（pythonで小さなツールを作る時のtips）

まえがき

はじめに

現在，有限要素法は，微分方程式の数値解法の一つとして，構造解析分野はもとより，浸透流（地下水流動）解析，熱伝導解析，流体解析など様々な分野で用いられています．
この記事では，特に筆者が土木分野での設計用ツールとして慣れ親しんできた，静力学分野での有限要素法の理論とPythonによる有限要素法プログラムについて，筆者なりの解説を試みたいと思います．

なお筆者は，土木構造物の設計技術者であり，教育機関の人間でもなければプログラミングの専門家でもありません．この記事に記載されている事項はあくまでも筆者の経験に基づく個人的な見解です．

ちなみに，２次元応力解析については，四角形（アイソパラメトリック）要素のほうが．三角形一次要素より性能がいいという記載をどこかで見たことがある気もします．また，筆者も四角形要素も使うには使うのですが，実は四角形要素用の手持ちのメッシャーの調子が今一つなことと（気分的にはこちらのほうが大きい），三角形要素のほうが基本的かな，と思うところがあり，ここでは三角形一次要素のプログラムを紹介することにしました．

Pythonでのプログラミングについて

Pythonを使うメリット

まずは，Pythonでプログラムを作って運用することのメリットについて，筆者なりの意見を述べてみたいと思います．

• プログラミング言語として：いずれかのプログラミング言語の経験があれば，特に行列・ベクトル演算が理解できていれば，比較的スムースに理解できる，使いこなせる言語だと，経験的に感じています．
• 数値解析ツールとして：ライブラリが充実している．連立一次方程式の解法，非線形方程式の解法，高速フーリエ変換，統計処理関数，特殊関数など，通常使いそうなものはほぼ揃っており，これらを活用することにより，自前のルーチンを書く必要がなくなり，コードをすっきりさせることができるとともに，本来の問題解決のプログラミングに集中できます．
• 総合的な技術ツールとして：作図ツールとしてmatplotlib，分析支援としてpandas，分析・解析作業の組立・実行のプラットフォームとしてJupyter notebookなど，関連ツールも充実している．またPythonの中から，ImageMagickなど他のアプリケーションを呼び出して実行できるし，Jupyter notebook上でも Shell script を実行できるなど，Pythonとその関連ツールを使うことにより，１つのプラットフォーム上で複数のアプリケーションを整然と利用できます．

以上，Pythonおよびその関連ツールの利用により，多くのことを自由にストレスなく行えることから，Pythonを使うメリットはあると思います．また有限要素法プログラムを考えた場合，基本的に行列とベクトル演算のプログラムとなるので，Pythonによる記述には向いていると考えています．実用上も，Pythonはインタープリタであることから，他のコンパイラ系言語に比べ速度的に不利な面もありますが，ライブラリの実行速度は十分に早く，小規模な設計計算であれば，データの準備から後処理まで考えれば総合的には所要時間的にもそれほど不利ではないように感じています．

自分のプログラムを持つメリット

筆者はこれまで現場の設計技術者という立場で設計を行ってきました．
基本的に大規模構造物の数値解析は，解析専門業者に外部委託するか，社内の専門家にお願いするのですが，現場で設計をやっていると，「あそこの設計を確認したいのだけど，手計算ではラフすぎるし，もう少し精度を上げたい」というようなことがよくあります．そんなとき，平面骨組や２次元応力解析のプログラムを持っていると，すぐに自前で計算できるし，思う存分モデルや条件を変えることができます．また走行式クレーンなどの移動荷重に対する設計など，やることは単純なのだけど多くのケースをこなさなければならない場合など，shell script でケースを回し，出力をPythonで重ね合わせるなどの加工も自由自在です．一言で言うと，「複雑なことはできないが，単純作業の重ね合わせなども含め，即応性・多様性・出力加工の自由度がある」というのが自前プログラムを持つメリットでしょうか．

なお，ついて回る問題としては，「その結果正しいの？」ということです．そこは色々なケースで使い込んで必要あれば直していくという対応しかないと思っています．まあ，そもそも自分が与えた荷重に対し，断面力図や主応力コンターが（定性的に，数値のオーダー的に）正しいかどうかなどは見ればわかるというのが技術者でしょう．ちなみに筆者が使っている平面骨組のプログラムは，言語はBasic, C, Fortran, Pythonと変わってきていますが，30年来使ってきているものです．平面骨組については，何度か若い人たちが使っている市販プログラムと結果を比べたことがありますが，メッシュの切り方や荷重の与え方のテクニックの影響もありますが，筆者の自作版ほうが高精度でした．

ということで，筆者の立場としては，「自前のプログラムを持って，疑問に感じたことはどんどん確認しよう」ということを勧めます．

参考文献

プログラミングの参考にした図書を以下に示します．30年程度以上前の文献ですが，これらをベースに作ったプログラムを書き換えながら使ってきました．プログラムの全体構成は，FORTRAN プログラムが掲載されている文献 (2) をベースにしています．また Python プログラム作成にあたっては，変数名，全体剛性方程式の組立方法など，N88-Basic プログラムが掲載されている文献 (1) を参考に設定しています．

1. 宮沢健二：パソコン構造力学--力学挙動の可視化--，啓学出版株式会社，1984年6月
2. 三本木茂夫・吉村信敏：コンピュータによる構造工学講座I-1-B 日本鋼構造協会編 有限要素法による構造解析プログラム 考え方と解説，培風館，昭和54年5月
3. H.C. Martin / G.F. Carey共著，鷲津久一郎/山本善之共訳：有限要素法の基礎と応用，培風館，昭和56年10月
4. 山田嘉昭：有限要素法の基礎と応用シリーズ２ マトリックス法材料力学，培風館，昭和55年5月
5. S.P. Timoshenko, J.N. Goodier : Theory of Elasticity Third edition, McGRAW-HILL BOOK CAMPANY, 21st printing 1985

関連リンク

• Python：FEMでの疎行列計算利用による高速化（３次元骨組構造解析での事例）
• Pythonによる２次元応力解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）
• Pythonによる平面骨組構造解析プログラムの理論を少し詳しく解説する（改訂版）

静力学分野における有限要素法の一般的な理論

前段が長くなりましたが，そろそろ有限要素法に関する記述に移っていこうと思います．

理論展開およびプログラムには，以下の条件が入っています．

• 等方性線形弾性材料の微小変形を仮定している
• 物体力として慣性力を考慮する
• 初期ひずみとして温度変化によるひずみを考慮する

仮想仕事の式による要素剛性方程式の標準的な形式の導出

静力学分野における有限要素解析は，節点荷重と節点変位の関係を示す連立一次方程式を作成し，これを解く作業であるということができる．
以下に，仮想仕事の原理から，標準的な要素剛性方程式，すなわち節点荷重と節点変位の関係を示す連立一次方程式を導く過程を説明する．

一般に，外力を受ける物体が静的釣り合い状態にある場合，物体内部の応力による仮想仕事と表面力や物体力などの外力による仮想仕事は等しい関係にある．この関係を式で表すと以下の通りとなる．

\begin{equation}
\int_V \{\boldsymbol{\delta \epsilon}\}\{\boldsymbol{\sigma}\}dV=\int_A \{\boldsymbol{\delta u}\}\{\boldsymbol{S}\}dA+\int_V \{\boldsymbol{\delta u}\}\{\boldsymbol{F}\}dV
\end{equation}

\begin{align}
&\{\boldsymbol{\delta \epsilon}\} & &\text{：応力方向の仮想ひずみ}\\
&\{\boldsymbol{\sigma}\}          & &\text{：物体内部応力}\\
&\{\boldsymbol{\delta u}\}        & &\text{：表面力方向の仮想変位}\\
&\{\boldsymbol{S}\}               & &\text{：表面力}\\
&\{\boldsymbol{F}\}               & &\text{：物体力（慣性力など）}\\
&V                                & &\text{：物体の体積}\\
&A                                & &\text{：表面力が作用する領域面積}\\
\end{align}

以降，一般的な静力学の設計に適用することを考慮し，物体力として慣性力（重力含む），初期ひずみとして温度ひずみを考えることとする．

ここで，物体を１つの有限要素とし，要素内任意点のひずみ${\boldsymbol{\epsilon}}$および要素内任意点の変位${\boldsymbol{u}}$を，以下のように節点変位${\boldsymbol{u_{nd}}}$で表せたとする．

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} \\
&\{\boldsymbol{u}\}=[\boldsymbol{N}]\{\boldsymbol{u_{nd}}\}
\end{align}

\begin{align}
&\{\boldsymbol{\epsilon}\} & &\text{：要素内任意点のひずみ}\\
&\{\boldsymbol{u}\}        & &\text{：要素内任意点の変位}\\
&[\boldsymbol{B}]          & &\text{：ひずみ-変位関係行列}\\
&[\boldsymbol{N}]          & &\text{：内挿関数行列}\\
&\{\boldsymbol{u_{nd}}\}   & &\text{：節点変位}
\end{align}

また以下に示す線形弾性体の応力-ひずみ関係を導入する．

\begin{equation}
\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
\end{equation}

\begin{align}
&\{\boldsymbol{\sigma}\}     & &\text{：要素応力}\\
&[\boldsymbol{D_e}]          & &\text{：応力-ひずみ関係行列}\\
&\{\boldsymbol{\epsilon}\}   & &\text{：外力による変位から求まるひずみ}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{：初期ひずみ（温度ひずみ）}
\end{align}

次に要素内任意点の物体力（慣性力）は，要素の単位体積重量，節点加速度を用いて以下のように表せるものとする．

\begin{equation}
\{\boldsymbol{F}\}=\gamma\cdot[\boldsymbol{N}]\{\boldsymbol{w_{nd}}\}
\end{equation}

\begin{align}
&\{\boldsymbol{F}\}      & &\text{：慣性力}\\
&[\boldsymbol{N}]        & &\text{：内挿関数}\\
&\{\boldsymbol{w_{nd}}\} & &\text{：節点加速度（ただし重力加速度との比）}\\
&\gamma                  & &\text{：要素の単位体積重量}
\end{align}

(note) ここで，上記のように定義することにより，単位体積重量 $\gamma$ は，$\gamma=m\cdot g$ （$m$：要素の質量，$g$：重力加速度）と表され，また節点加速度 ${\boldsymbol{w_{nd}}}$ を重力加速度との比 $w_{nd}=a :/: g$ （$a$：要素に働く加速度）とすることにより，次元として，$F=m\cdot g \times a :/: g=m\cdot a$ というNewtonの第二法則が成り立つことになる．

上記関係を用いることにより，仮想仕事式の左辺（内力による仮想仕事）は以下のようになる．

\begin{align}
\int_V \{\boldsymbol{\delta \epsilon}\}\{\boldsymbol{\sigma}\}dV
=&\{\boldsymbol{\delta u_{nd}}\}^T \int_v[\boldsymbol{B}]^T\{\boldsymbol{\sigma}\}dV \\
=&\{\boldsymbol{\delta u_{nd}}\}^T \left(\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV\right) \{\boldsymbol{u_{nd}}\}
-\{\boldsymbol{\delta u_{nd}}\}^T \left(\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV\right) \\
=&\{\boldsymbol{\delta u_{nd}}\}^T[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}-\{\boldsymbol{\delta u_{nd}}\}^T\{\boldsymbol{f_t}\}
\end{align}

また，仮想仕事式の右辺（外力による仮想仕事）は以下のようになる．

\begin{align}
&\int_A \{\boldsymbol{\delta u}\}\{\boldsymbol{S}\}dA=\{\boldsymbol{\delta u_{nd}}\}^T \left(\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA\right)=\{\boldsymbol{\delta u_{nd}}\}^T\{\boldsymbol{f}\} \\
&\int_V \{\boldsymbol{\delta u}\}\{\boldsymbol{F}\}dV=\{\boldsymbol{\delta u_{nd}}\}^T \left(\gamma \int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\}=\{\boldsymbol{\delta u_{nd}}\}^T\{\boldsymbol{f_b}\}
\end{align}

以上により，次に示すような，一般的な剛性方程式，すなわち節点荷重と節点変位の関係を表す連立一次方程式が得られる．

\begin{equation}
\{\boldsymbol{f}\}+\{\boldsymbol{f_t}\}+\{\boldsymbol{f_b}\}=[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}
\end{equation}

\begin{align}
&[\boldsymbol{k}]=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV  & & \text{（要素剛性行列）}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{（節点外力ベクトル）}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{（節点初期ひずみベクトル）}\\
&\{\boldsymbol{f_b}\}=\left(\gamma \int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{（節点物体力ベクトル）}
\end{align}

上式は，変位を求めるための式であり，初期ひずみは荷重項として扱っている．要素内部応力は，前出式

\begin{equation}
\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
\end{equation}

を用いて計算する必要があることに注意する．ここに，${\boldsymbol{\epsilon}}$ は，剛性方程式を解いて求めた変位から求めた要素のひずみベクトルである．

有限要素法においては，解きたい問題に合わせ，トラス要素，梁要素，三角形要素，四角形要素などを選定し，その要素の特性を折り込みながら，上述した剛性方程式の中身を具体化していくことになる．

内挿関数の求め方（概念的な説明）

次に，色々な要素で共通の考え方となる，内挿関数の求め方の概念的な説明を行っておく．

以下に示すように，要素内任意点の変位 ${\boldsymbol{u}}$ を節点変位 ${\boldsymbol{u_{nd}}}$ で表すための内挿関数 $[\boldsymbol{N}]$ を求めることを考える．

\begin{equation}
\{\boldsymbol{u}\}=[\boldsymbol{N}]\{\boldsymbol{u_{nd}}\}
\end{equation}

要素内の任意点における変位 ${\boldsymbol{u}}$ が，要素節点と同じ自由度を有する列ベクトル ${\boldsymbol{\alpha}}$ を用いて以下のように表せるものとする．ここに $[\boldsymbol{X}]$ は，要素内座標を含む行列である．

\begin{equation}
\{\boldsymbol{u}\}=[\boldsymbol{X}]\{\boldsymbol{\alpha}\}
\end{equation}

上記関係を節点での値 ${\boldsymbol{u_{nd}}}$ に適用した場合，以下の関係が成り立つ．

\begin{equation}
\{\boldsymbol{u_{nd}}\}=[\boldsymbol{X_{nd}}]\{\boldsymbol{\alpha}\}
\end{equation}

上式を ${\boldsymbol{\alpha}}$ について解けば，

\begin{equation}
\{\boldsymbol{\alpha}\}=[\boldsymbol{C}]\{\boldsymbol{u_{nd}}\} \qquad [\boldsymbol{C}]=[\boldsymbol{X_{nd}}]^{-1}
\end{equation}

よって，

\begin{equation}
\{\boldsymbol{u}\}=[\boldsymbol{X}]\{\boldsymbol{\alpha}\}=[\boldsymbol{X}][\boldsymbol{C}]\{\boldsymbol{u_{nd}}\}
\end{equation}

であり，内挿関数 $[\boldsymbol{N}]$ は，以下のように求めることができる．

\begin{equation}
[\boldsymbol{N}]=[\boldsymbol{X}][\boldsymbol{C}]
\end{equation}

剛性マトリックスの組み立て

有限要素法における全体剛性方程式は，一般に以下の形をしており，その次元は $N=\text{(１節点自由度)}\times\text{(総節点数)}$ である．

\begin{equation}
\begin{Bmatrix}
F_1 \\
\vdots \\
F_i \\
\vdots \\
F_j \\
\vdots \\
F_k \\
\vdots \\
F_N
\end{Bmatrix}
=\begin{bmatrix}
K_{11} & \ldots & K_{1i} & \ldots & K_{1j} & \ldots & K_{1k} & \ldots & K_{1N} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \ldots & \ldots & \vdots \\
K_{i1} & \ldots & K_{ii} & \ldots & K_{ij} & \ldots & K_{ik} & \ldots & K_{iN} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \ldots & \ldots & \vdots \\
K_{j1} & \ldots & K_{ji} & \ldots & K_{jj} & \ldots & K_{jk} & \ldots & K_{jN} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \ldots & \ldots & \vdots \\
K_{k1} & \ldots & K_{ki} & \ldots & K_{kj} & \ldots & K_{kk} & \ldots & K_{kN} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \ldots & \ldots & \vdots \\
K_{N1} & \ldots & K_{Ni} & \ldots & K_{Nj} & \ldots & K_{Nk} & \ldots & K_{NN}
\end{bmatrix}
\begin{Bmatrix}
U_1 \\
\vdots \\
U_i \\
\vdots \\
U_j \\
\vdots \\
U_k \\
\vdots \\
U_N
\end{Bmatrix}
\end{equation}

これに対し，要素剛性方程式は，$n=\text{(１節点自由度)}\times\text{(１要素節点数)}$ の次元を有しており，これが要素数分存在するわけであるが，これらの共通の節点を持つインデックスの配列要素どうしを重ね合わせて，全体剛性方程式を作成する．

２次元応力解析のための三角形一次要素を例に，要素剛性方程式の具体的な形を説明する．
２次元応力解析用の三角形一次要素は，１要素３節点であり，平面問題のため，１節点２自由度（x方向変位およびy方向変位）を有する．
よって要素剛性方程式の次元は６（１節点２自由度ｘ１要素３節点）であり，具体的には以下の形をとる．ここに，三角形要素を構成する節点番号を反時計回りに，$i$，$j$，$k$ とする．

\begin{equation}
\begin{Bmatrix} f_{x,i} \\ f_{y,i} \\ f_{x,j} \\ f_{y,j} \\ f_{x,k} \\ f_{y,k} \end{Bmatrix}
=\begin{bmatrix}
k_{11} & k_{12} & k_{13} & k_{14} & k_{15} & k_{16} \\
k_{21} & k_{22} & k_{23} & k_{24} & k_{25} & k_{26} \\
k_{31} & k_{32} & k_{33} & k_{34} & k_{35} & k_{36} \\
k_{41} & k_{42} & k_{43} & k_{44} & k_{45} & k_{46} \\
k_{51} & k_{52} & k_{53} & k_{54} & k_{55} & k_{56} \\
k_{61} & k_{62} & k_{63} & k_{64} & k_{65} & k_{66} 
\end{bmatrix}
\begin{Bmatrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \end{Bmatrix}
\end{equation}

\begin{align}
&f_{x,i} & &\text{：節点$i$の$x$方向荷重} & & u_i & &\text{：節点$i$の$x$方向変位}\\
&f_{y,i} & &\text{：節点$i$の$y$方向荷重} & & v_i & &\text{：節点$i$の$y$方向変位}\\
&f_{x,j} & &\text{：節点$j$の$x$方向荷重} & & u_i & &\text{：節点$j$の$x$方向変位}\\
&f_{y,j} & &\text{：節点$j$の$y$方向荷重} & & v_i & &\text{：節点$j$の$y$方向変位}\\
&f_{x,k} & &\text{：節点$k$の$x$方向荷重} & & u_i & &\text{：節点$k$の$x$方向変位}\\
&f_{y,k} & &\text{：節点$k$の$y$方向荷重} & & v_i & &\text{：節点$k$の$y$方向変位}
\end{align}

上記 6 x 6 の次元の要素剛性マトリックスを NｘN の次元の全体剛性マトリックスにどのように配置するかを示したのが下図である．碁盤目に区切られた大きな正方形は全体剛性マトリックス（N x N）を示し，$i$，$j$，$k$ は要素を構成する節点番号，$k_{11}\sim k_{66}$ は要素剛性マトリックス（6 x 6）の成分を示す．

この操作を要素数分繰り返し足しあわせていくことにより，全体剛性マトリックスが完成する．

荷重ベクトル（列ベクトル）についても同様に，節点 $i$ に関連する行は，$2i-1$ 行と $2i$ 行になるので，$2+i-1$ 行に $f_{x,i}$ を，$2*i$ 行に $f_{y,i}$ を足し込む．節点 $j$ および節点 $k$ に関しても同様である．

なお，Pythonプログラムでは，配列要素をインデックス０から使用しているので，実装においては，インデックスが１ずれることに注意する．

連立一次方程式における既知量と未知量の取扱

一般に，有限要素法では，以下の形で表される連立一次方程式を解いている．

\begin{equation}
\{\boldsymbol{f}\}=[\boldsymbol{k}]\{\boldsymbol{u}\}
\end{equation}

ここに，応力解析の例を考えれば，$[\boldsymbol{k}]$ は全体剛性マトリックス，${\boldsymbol{u}}$ は節点変位ベクトル，${\boldsymbol{f}}$ は節点荷重ベクトルに相当する．
通常は，上式において，荷重ベクトル ${\boldsymbol{f}}$ を入力し，連立一次方程式を解くことにより，変位ベクトル ${\boldsymbol{u}}$ を求める作業を行うが，全体系で組み立てた剛性マトリックスは，そのままでは特異である．すなわち，力学的には変位拘束されていないため，構造系は反力をとれず，荷重により動いていってしまい，静止した釣り合い状態を保てない．そこで，変位拘束を行う必要があるわけである．

この変位拘束については，一般性を持たせるためには，必ずしもゼロでない変位（強制変位）を扱えるようにしたほうが，プログラムとしての使い勝手が良い．そこで，ここに述べるような処理が必要となってくる．

このような関係は，浸透流解析における全水頭と流量の関係，熱伝導解析における温度と熱流束の関係においても同様であり，プログラム上，同一の処理を行う必要がある．

ここで紹介しているFEM プログラムでは，全体剛性マトリックスの大きさを変えずに境界条件（変位指定）を導入し，連立方程式を解いている．この方法を用いる場合，連立一次方程式において，未知数と既知数の入れ替えが必要であるため，その考え方を以下に示す．

簡単な事例として，以下の連立一次方程式を考える．

\begin{align}
& k_{11} u_1 + k_{12} u_2 + k_{13} u_3 = f_1 \\
& k_{21} u_1 + k_{22} u_2 + k_{23} u_3 = f_2 \\
& k_{31} u_1 + k_{32} u_2 + k_{33} u_3 = f_3
\end{align}

ここで，未知数が $u_1, u_3, f_2$，既知数が $f_1, f_3, u_2$だとすれば，行列演算を行うためには，以下のように，未知数を左辺に，既知数を右辺に集めるような，書き換えが必要となる．

\begin{align}
& k_{11} u_1 + 0   + k_{13} u_3 = f_1 - k_{12} u_2 \\
& k_{21} u_1 - f_2 + k_{23} u_3 = 0   - k_{22} u_2 \\
& k_{31} u_1 + 0   + k_{33} u_3 = f_3 - k_{32} u_2
\end{align}

上記関係を，拡張して考えると，以下の通りとなる．

元の剛性方程式

\begin{equation}
\begin{bmatrix}
k_{11} & \ldots & k_{1i} & \ldots & k_{1j} & \ldots & k_{1n} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{i1} & \ldots & k_{ii} & \ldots & k_{ij} & \ldots & k_{in} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{j1} & \ldots & k_{ji} & \ldots & k_{jj} & \ldots & k_{jn} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{n1} & \ldots & k_{ni} & \ldots & k_{nj} & \ldots & k_{nn}
\end{bmatrix}
\begin{Bmatrix}
u_1 \\
\vdots \\
u_i \\
\vdots \\
u_j \\
\vdots \\
u_n
\end{Bmatrix}
=\begin{Bmatrix}
f_1 \\
\vdots \\
f_i \\
\vdots \\
f_j \\
\vdots \\
f_n
\end{Bmatrix}
\end{equation}

境界条件導入後の剛性方程式

$u_i$ および $u_j$ が既知とすれば，剛性マトリックスの $k_{ii}$ および $k_{jj}$ を $1$ とし，それ以外の $i$ 列および $j$ 列の要素をゼロとする処理を行う．また剛性マトリックスにおける $i$ 列および $j$ 列の効果を右辺に移項する必要がある．

\begin{equation}
\begin{bmatrix}
k_{11} & \ldots & 0      & \ldots & 0 & \ldots & k_{1n} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{i1} & \ldots & 1      & \ldots & 0      & \ldots & k_{in} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{j1} & \ldots & 0      & \ldots & 1      & \ldots & k_{jn} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{n1} & \ldots & 0      & \ldots & 0      & \ldots & k_{nn}
\end{bmatrix}
\begin{Bmatrix}
u_1 \\
\vdots \\
-f_i \\
\vdots \\
-f_j \\
\vdots \\
u_n
\end{Bmatrix}
=\begin{Bmatrix}
f_1 \\
\vdots \\
0 \\
\vdots \\
0 \\
\vdots \\
f_n
\end{Bmatrix}
-u_i
\begin{Bmatrix}
k_{1i} \\
\vdots \\
k_{ii} \\
\vdots \\
k_{ji} \\
\vdots \\
k_{ni}
\end{Bmatrix}
-u_j
\begin{Bmatrix}
k_{1j} \\
\vdots \\
k_{ij} \\
\vdots \\
k_{jj} \\
\vdots \\
k_{nj}
\end{Bmatrix}
\end{equation}

疎行列処理導入による連立１次方程式処理の高速化

前述の操作により，剛性マトリックスは非対称となるが，1000自由度程度の連立一次方程式であれば，numpy.linalg.solve(k, f) で大きなストレスなく解くことができる．ただし，自由度が1000を超えてくると，numpy.linalg.solve では処理時間が気になりだすため，疎行列処理を施した上で scipy.sparse.linalg.spsolv(k,f) を用いれば，計算時間短縮が図れる．プログラミング上の処理は以下のように行えば良い．

# 以下のインポートを追加
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix

# 連立方程式を解く部分で，剛性行列を csr_matrix に変換し，ソルバーに spsolve を採用
# solution of simultaneous linear equations
#disg = np.linalg.solve(gk, fp)
gk = csr_matrix(gk)
disg = spsolve(gk, fp, use_umfpack=True)

三角形一次要素を用いた２次元応力解析において，自由度10000のモデルで荷重条件を変えて３ケースの計算を行った計算時間実測例を以下に示す．疎行列処理を行ったケースのほうが圧倒的に早くなっている．

-----------------------------------------------------------
 自由度    numpy.linalg.solve   scipy.sparse.linalg.spsolve
             (疎行列処理なし)       (疎行列処理あり)
-----------------------------------------------------------
 n=10000      time=45.600 sec         time=5.154 sec
 n=10000      time=47.962 sec         time=5.265 sec
 n=10000      time=42.549 sec         time=5.680 sec
-----------------------------------------------------------

（参考リンク）Python：FEMでの疎行列計算利用による高速化（３次元骨組構造解析での事例）

三角形一次要素を用いた２次元応力解析の理論

比較的汎用性が高く基本的なケースである三角形一次要素を用いた２次元応力解析プログラムについて，解説する．

要素剛性方程式の基本形

直交座標系$(x-y)$平面上の三角形一次要素を考える．
２次元応力解析のための三角形一次要素は１要素３節点であり，平面問題のため，１節点２自由度（x方向変位およびy方向変位）を有する．
よって基本的な要素剛性方程式は，以下の形をとる．ここに，三角形要素を構成する節点番号を反時計回りに，$i$，$j$，$k$ とする．

\begin{equation}
\{\boldsymbol{f}\}=[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}
\end{equation}

\begin{equation}
\begin{Bmatrix} f_{x,i} \\ f_{y,i} \\ f_{x,j} \\ f_{y,j} \\ f_{x,k} \\ f_{y,k} \end{Bmatrix}
=\begin{bmatrix}
k_{11} & k_{12} & k_{13} & k_{14} & k_{15} & k_{16} \\
k_{21} & k_{22} & k_{23} & k_{24} & k_{25} & k_{26} \\
k_{31} & k_{32} & k_{33} & k_{34} & k_{35} & k_{36} \\
k_{41} & k_{42} & k_{43} & k_{44} & k_{45} & k_{46} \\
k_{51} & k_{52} & k_{53} & k_{54} & k_{55} & k_{56} \\
k_{61} & k_{62} & k_{63} & k_{64} & k_{65} & k_{66} 
\end{bmatrix}
\begin{Bmatrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \end{Bmatrix}
\end{equation}

\begin{align}
&f_{x,i} & &\text{：節点$i$の$x$方向荷重} & & u_i & &\text{：節点$i$の$x$方向変位}\\
&f_{y,i} & &\text{：節点$i$の$y$方向荷重} & & v_i & &\text{：節点$i$の$y$方向変位}\\
&f_{x,j} & &\text{：節点$j$の$x$方向荷重} & & u_i & &\text{：節点$j$の$x$方向変位}\\
&f_{y,j} & &\text{：節点$j$の$y$方向荷重} & & v_i & &\text{：節点$j$の$y$方向変位}\\
&f_{x,k} & &\text{：節点$k$の$x$方向荷重} & & u_i & &\text{：節点$k$の$x$方向変位}\\
&f_{y,k} & &\text{：節点$k$の$y$方向荷重} & & v_i & &\text{：節点$k$の$y$方向変位}
\end{align}

線形弾性論による基本式

平面応力状態での応力ーひずみ関係式

\begin{equation}
\begin{Bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{Bmatrix}
=\frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix}
\begin{Bmatrix} \epsilon_x-\alpha T \\ \epsilon_y-\alpha T \\ \gamma_{xy} \end{Bmatrix}
\end{equation}

平面ひずみ状態での応力ーひずみ関係式

\begin{equation}
\begin{Bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{Bmatrix}
=\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \cfrac{1-2\nu}{2} \end{bmatrix}
\begin{Bmatrix} \epsilon_x-(1+\nu)\alpha T \\ \epsilon_y-(1+\nu)\alpha T \\ \gamma_{xy} \end{Bmatrix}
\end{equation}

ひずみー変位関係式

\begin{equation}
\{\boldsymbol{\epsilon}\}
=\begin{Bmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{Bmatrix}
=\begin{Bmatrix} \cfrac{\partial u}{\partial x} \\ \cfrac{\partial v}{\partial y} \\ \cfrac{\partial u}{\partial y}+\cfrac{\partial v}{\partial x} \end{Bmatrix}
\end{equation}

記号の定義

\begin{align}
&\sigma_x     & &\text{：x方向直応力}    & & \epsilon_x  & &\text{：x方向直ひずみ}\\
&\sigma_y     & &\text{：y方向直応力}    & & \epsilon_y  & &\text{：y方向直ひずみ}\\
&\tau_{xy}    & &\text{：せん断応力}     & & \gamma_{xy} & &\text{：せん断ひずみ}\\
&E            & &\text{：要素弾性係数}   & & \nu         & &\text{：要素ポアソン比}\\
&\alpha       & &\text{：要素線膨張係数} & & T           & &\text{：要素温度上昇量}\\
&u            & &\text{：x方向変位}      & & v           & &\text{：y方向変位}
\end{align}

以下に，プログラムを作成するために必要な，具体的な要素を求めていく．

内挿関数

２次元応力解析に用いる三角形一次要素では，１節点２自由度３節点であるため，１要素あたり６自由度を有する．
このため，６個の未定係数$\alpha_{1\sim6}$を用い，以下の変位を仮定する．

\begin{align}
u=&\alpha_1+\alpha_2\cdot x+\alpha_3\cdot y \\
v=&\alpha_4+\alpha_5\cdot x+\alpha_6\cdot y
\end{align}

すなわち

\begin{equation}
\begin{Bmatrix} u \\ v \end{Bmatrix}
=[\boldsymbol{X}]\{\boldsymbol{\alpha}\}
=\begin{bmatrix}
1 & x & y & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & x & y
\end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \alpha_5 \\ \alpha_6 \end{Bmatrix}
\end{equation}

ここで，三角形の３つの頂点 $i$，$j$，$k$ を反時計回りに設定し，それぞれの節点変位 $(u,v)$ を節点座標 $(x,y)$ と $\alpha$ で表せば以下の通りとなる．

\begin{equation}
\begin{Bmatrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \end{Bmatrix}
=\begin{bmatrix}
1 & x_i & y_i & 0 & 0   & 0   \\
0 & 0   & 0   & 1 & x_i & y_i \\
1 & x_j & y_j & 0 & 0   & 0   \\
0 & 0   & 0   & 1 & x_j & y_j \\
1 & x_k & y_k & 0 & 0   & 0   \\
0 & 0   & 0   & 1 & x_k & y_k
\end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \alpha_5 \\ \alpha_6 \end{Bmatrix}
\end{equation}

これを ${\boldsymbol{\alpha}}$ について解けば

\begin{equation}
\{\boldsymbol{\alpha}\}=[\boldsymbol{C}]\{\boldsymbol{u_{nd}}\}
\end{equation}

\begin{equation}
[\boldsymbol{C}]=\frac{1}{2\Delta}
\begin{bmatrix}
a_i & 0   & a_j & 0   & a_k & 0   \\
b_i & 0   & b_j & 0   & b_k & 0   \\
c_i & 0   & c_j & 0   & c_k & 0   \\
0   & a_i & 0   & a_j & 0   & a_k \\
0   & b_i & 0   & b_j & 0   & b_k \\
0   & c_i & 0   & c_j & 0   & c_k
\end{bmatrix}
\qquad
\qquad
2\Delta=
\begin{vmatrix}
1 & x_i & y_i \\
1 & x_j & y_j \\
1 & x_k & y_k \\
\end{vmatrix}
\end{equation}

\begin{equation}
\begin{matrix}
a_i=x_j y_k-x_k y_j & \qquad & b_i=y_j-y_k & \qquad & c_i=x_k-x_j \\
a_j=x_k y_i-x_i y_k & \qquad & b_j=y_k-y_i & \qquad & c_j=x_i-x_k \\
a_k=x_i y_j-x_j y_i & \qquad & b_k=y_i-y_j & \qquad & c_k=x_j-x_i
\end{matrix}
\end{equation}

ここに，$\Delta$ は三角形要素の面積，${\boldsymbol{u_{nd}}}$は節点変位ベクトル，${\boldsymbol{\alpha}}$は一般化変位ベクトルである．$2\Delta$ を直接の計算式で示せば，以下の通りとなる．

\begin{equation}
2\Delta=(x_k-x_j)y_i+(x_i-x_k)y_j+(x_j-x_i)y_k
\end{equation}

以上より，２次元応力解析用三角形一次要素の内挿関数は以下の通り定まる．

\begin{equation}
[\boldsymbol{N}]=[\boldsymbol{X}][\boldsymbol{C}]
=\begin{bmatrix}
N_i & 0   & N_j & 0   & N_k & 0   \\
0   & N_i & 0   & N_j & 0   & N_k
\end{bmatrix}
\end{equation}

\begin{equation}
N_i=\frac{a_i+b_i x+c_i y}{2\Delta} \qquad
N_i=\frac{a_j+b_j x+c_j y}{2\Delta} \qquad
N_i=\frac{a_k+b_k x+c_k y}{2\Delta}
\end{equation}

ひずみ-変位関係行列

要素ひずみ ${\boldsymbol{\epsilon}}$ は変位 $(u,v)$ を微分することにより以下の通り表される．

\begin{equation}
\{\boldsymbol{\epsilon}\}=
\begin{Bmatrix}
\cfrac{\partial u}{\partial x} \\
\cfrac{\partial v}{\partial y} \\
\cfrac{\partial u}{\partial y}+\cfrac{\partial v}{\partial x}
\end{Bmatrix}
=\begin{bmatrix}
\cfrac{\partial N_i}{\partial x} & 0 & \cfrac{\partial N_j}{\partial x} & 0 & \cfrac{\partial N_k}{\partial x} & 0 \\
0 & \cfrac{\partial N_i}{\partial y} & 0 & \cfrac{\partial N_j}{\partial y} & 0 & \cfrac{\partial N_k}{\partial y} \\
\cfrac{\partial N_i}{\partial y} & \cfrac{\partial N_i}{\partial x} &
\cfrac{\partial N_j}{\partial y} & \cfrac{\partial N_j}{\partial x} &
\cfrac{\partial N_k}{\partial y} & \cfrac{\partial N_k}{\partial x}
\end{bmatrix}
\{\boldsymbol{u_{nd}}\}
=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\}
\end{equation}

よってひずみ-変位関係は以下の通りとなる．

\begin{align}
[\boldsymbol{B}]
=&\frac{1}{2\Delta}
\begin{bmatrix}
b_i & 0   & b_j & 0   & b_k & 0   \\
0   & c_i & 0   & c_j & 0   & c_k \\
c_i & b_i & c_j & b_j & c_k & b_k \\
\end{bmatrix}\notag \\
=&\frac{1}{2\Delta}
\begin{bmatrix}
y_j-y_k & 0       & y_k-y_i & 0       & y_i-y_j & 0       \\
0       & x_k-x_j & 0       & x_i-x_k & 0       & x_j-x_i \\
x_k-x_j & y_j-y_k & x_i-x_k & y_k-y_i & x_j-x_i & y_i-y_j
\end{bmatrix} \\
\Delta=&1/2\cdot[(x_k-x_j)y_i+(x_i-x_k)y_j+(x_j-x_i)y_k] \qquad \text{（要素面積）}
\end{align}

応力-ひずみ関係式

応力ーひずみ関係マトリックスは，以下の通り．

\begin{align}
&\text{平面応力：}& &[\boldsymbol{D_e}]=\frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix} \qquad
\begin{matrix} E & \text{：弾性係数} \\ \nu & \text{：ポアソン比} \end{matrix} \\
&\text{平面ひずみ：}& &[\boldsymbol{D_e}]=\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \cfrac{1-2\nu}{2} \end{bmatrix}
\end{align}

要素剛性マトリックス

要素剛性マトリックス$[\boldsymbol{k}]$は，ひずみ-変位関係行列$[\boldsymbol{B}]$および応力-ひずみ関係行列$[\boldsymbol{D_e}]$を用い，以下の通りとなる．一般的に要素剛性マトリックスは体積積分により求めるが，２次元弾性問題において，要素厚さを一定値 $t$ とすれば，要素体積 $t\cdot \Delta$ （$\Delta$ は三角形要素面積）に定数マトリックスである $[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]$ を乗じれば良い．

\begin{equation}
[\boldsymbol{k}]=t\cdot\int_A[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dA=t\cdot\Delta\cdot[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]
\end{equation}

節点温度荷重ベクトル

温度変化による節点温度荷重ベクトルは，以下の通りとなる．

\begin{equation}
\{\boldsymbol{f_t}\}=t\cdot\int_A[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dA=t\cdot\Delta\cdot[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}
\end{equation}

ここに，温度にずみ（初期ひずみ） ${\boldsymbol{\epsilon_0}}$ は以下のように求められる．平面応力状態と平面ひずみ状態では，違いがあることに注意する．

\begin{align}
&\text{要素平均温度変化量：}& & T_m=\cfrac{\phi_i+\phi_j+\phi_k}{3} \\
&\text{平面応力：}& & \{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_m \\ \alpha\cdot T_m \\ 0 \end{Bmatrix}
\\
&\text{平面ひずみ：}& & \{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_m \\ (1+\nu)\cdot \alpha\cdot T_m \\ 0 \end{Bmatrix}
\end{align}

ここに，$\alpha$は線膨張係数，$\nu$はポアソン比，$T_m$は要素を取り囲む節点の平均温度変化量（要素平均温度変化量とみなす），${\phi_i,\phi_j,\phi_k}^T$は節点温度変化量である．また，温度変化量の符号は，温度上昇を正とする．

節点慣性力ベクトル

面積座標で表された内挿関数の積分

慣性力など物体力による節点力ベクトルを求める場合，内挿関数の積分が必要となる．
このため，内挿関数の書き換えを行う．

\begin{equation}
[\boldsymbol{N}]=
\begin{bmatrix}
N_i & 0   & N_j & 0   & N_k & 0   \\
0   & N_i & 0   & N_j & 0   & N_k
\end{bmatrix}
=\begin{bmatrix}
\zeta_i & 0       & \zeta_j & 0       & \zeta_k & 0       \\
0       & \zeta_i & 0       & \zeta_j & 0       & \zeta_k
\end{bmatrix}
\end{equation}

\begin{align}
\zeta_i=\frac{\Delta_i}{\Delta} \qquad \qquad
\zeta_j=\frac{\Delta_j}{\Delta} \qquad \qquad
\zeta_k=\frac{\Delta_k}{\Delta} \\
2\Delta_i=a_i+b_i x+c_i y=
\begin{vmatrix}
1 & x   & y   \\
1 & x_j & y_j \\
1 & x_k & y_k \\
\end{vmatrix}\\
2\Delta_j=a_j+b_j x+c_j y=
\begin{vmatrix}
1 & x   & y   \\
1 & x_k & y_k \\
1 & x_i & y_i \\
\end{vmatrix}\\
2\Delta_k=a_k+b_k x+c_k y=
\begin{vmatrix}
1 & x   & y   \\
1 & x_i & y_i \\
1 & x_j & y_j \\
\end{vmatrix}
\end{align}

この座標 $(\zeta_i,\zeta_j,\zeta_k)$ は面積座標を表し，明らかに

\begin{equation}
\zeta_i+\zeta_j+\zeta_k=1 \qquad \qquad 0\leqq\zeta_r\leqq 1 \quad (r=i,j,k)
\end{equation}

である．

物体力を求める場合に必要となる，$\int_A[N]^T[N]dA$ の形の積分を行うため，まず $[N]^T[N]$ を求める．

\begin{equation}
[\boldsymbol{N}]^T[\boldsymbol{N}]=
\begin{bmatrix}
\zeta_i & 0       \\
0       & \zeta_i \\
\zeta_j & 0       \\
0       & \zeta_j \\
\zeta_k & 0       \\
0       & \zeta_k \\
\end{bmatrix}
\begin{bmatrix}
\zeta_i & 0       & \zeta_j & 0       & \zeta_k & 0       \\
0       & \zeta_i & 0       & \zeta_j & 0       & \zeta_k
\end{bmatrix}=
\begin{bmatrix}
\zeta_i{}^2    & 0              & \zeta_i\zeta_j & 0              & \zeta_k\zeta_i & 0              \\
0              & \zeta_i{}^2    & 0              & \zeta_i\zeta_j & 0              & \zeta_k\zeta_i \\
\zeta_i\zeta_j & 0              & \zeta_j{}^2    & 0              & \zeta_j\zeta_k & 0              \\
0              & \zeta_i\zeta_j & 0              & \zeta_j{}^2    & 0              & \zeta_j\zeta_k \\
\zeta_k\zeta_i & 0              & \zeta_j\zeta_k & 0              & \zeta_k{}^2    & 0              \\
0              & \zeta_k\zeta_i & 0              & \zeta_j\zeta_k & 0              & \zeta_k{}^2
\end{bmatrix}
\end{equation}

積分変数の変換公式より

\begin{equation}
\iint_A f(x,y) dx dy = \int_0^1\int_0^1 g(\zeta_i,\zeta_j,\zeta_k)\frac{\partial(x,y)}{\partial(\zeta_i,\zeta_j)} d\zeta_i d\zeta_j
\end{equation}

要素内任意点を内挿関数で表し，$\zeta_i+\zeta_j+\zeta_k=1$の関係を用いると

\begin{align}
x=&\zeta_i x_i+\zeta_j x_j+\zeta_k x_k \qquad &\rightarrow\qquad &x=\zeta_i(x_i-x_k)+\zeta_j(x_j-x_k)+x_k \notag \\
y=&\zeta_i y_i+\zeta_j y_j+\zeta_k y_k \qquad &\rightarrow\qquad &y=\zeta_i(y_i-y_k)+\zeta_j(y_j-y_k)+y_k
\end{align}

ヤコビアンは

\begin{align}
\frac{\partial(x,y)}{\partial(\zeta_i,\zeta_j)}=&
\begin{vmatrix}
\cfrac{\partial x}{\partial\zeta_i} & \cfrac{\partial x}{\partial\zeta_j} \\
\cfrac{\partial y}{\partial\zeta_i} & \cfrac{\partial y}{\partial\zeta_j}
\end{vmatrix}=
\begin{vmatrix}
x_i-x_k & x_j-x_k \\
y_i-y_k & y_j-y_k
\end{vmatrix} \notag \\
=&(x_k-x_j)y_i+(x_i-x_k)y_j+(x_j-x_i)y_k=2\Delta
\end{align}

よって

\begin{equation}
\iint_A f(x,y) dx dy = 2\Delta\int_0^1\int_0^1 g(\zeta_i,\zeta_j,\zeta_k) d\zeta_i d\zeta_j
\end{equation}

また面積座標 $(\zeta_i,\zeta_j,\zeta_k)$ に対し，以下の公式が成り立つ．

\begin{equation}
\int_0^1\int_0^1 \zeta_i{}^{\ell}\cdot \zeta_j{}^m\cdot \zeta_k{}^n \cdot d\zeta_i d\zeta_j
=\frac{\ell !\cdot m! \cdot n!}{(\ell+m+n+2)!} \qquad \text{（$\ell,m,n$は負でない整数）}
\end{equation}

たとえば

\begin{align}
&\ell=2, m=n=0     & \rightarrow \int_0^1\int_0^1 \zeta_i{}^{\ell}\cdot \zeta_j{}^m\cdot \zeta_k{}^n \cdot d\zeta_i d\zeta_j=1/12 \notag \\
&\ell=1, m=1, n=0  & \rightarrow \int_0^1\int_0^1 \zeta_i{}^{\ell}\cdot \zeta_j{}^m\cdot \zeta_k{}^n \cdot d\zeta_i d\zeta_j=1/24 \notag
\end{align}

したがって（$\zeta$での積分ではなく$A$での積分であることに注意）

\begin{align}
&\int_A \zeta_i{}^2 dA=\int_A \zeta_j{}^2 dA=\int_A \zeta_k{}^2 dA=1/12\cdot 2\Delta=1/6\cdot\Delta \notag \\
&\int_A \zeta_i\zeta_j dA=\int_A \zeta_j\zeta_k dA=\int_A \zeta_k\zeta_i dA=1/24\cdot 2\Delta=1/12\cdot\Delta \notag
\end{align}

以上により

\begin{equation}
\int_A [\boldsymbol{N}]^T[\boldsymbol{N}] dA=\Delta
\begin{bmatrix}
1/6  & 0    & 1/12 & 0    & 1/12 & 0    \\
0    & 1/6  & 0    & 1/12 & 0    & 1/12 \\
1/12 & 0    & 1/6  & 0    & 1/12 & 0    \\
0    & 1/12 & 0    & 1/6  & 0    & 1/12 \\
1/12 & 0    & 1/12 & 0    & 1/6  & 0    \\\
0    & 1/12 & 0    & 1/12 & 0    & 1/6
\end{bmatrix}
\end{equation}

節点慣性力ベクトル

物体力（慣性力）による節点力ベクトル ${\boldsymbol{f_{b}}}$ は，

\begin{equation}
\{\boldsymbol{f_{b}}\}
=\gamma\cdot t\cdot\int_{A}[\boldsymbol{N}]^T[\boldsymbol{N}]dA\cdot\{\boldsymbol{w_{nd}}\}
\end{equation}

と表せる．ここに，$\gamma$ は要素の単位体積重量（質量$\times$重力加速度），$t$ は要素板厚である．

ここで，ベクトル${\boldsymbol{w_ {nd}}}$を以下の通りおく．

\begin{equation}
\{\boldsymbol{w_{nd}}\}=\begin{Bmatrix} k_h & k_v & k_h & k_v & k_h & k_v \end{Bmatrix}^T
\end{equation}

ここに $k_h$ は水平方向加速度の重力加速度に対する比率（右向き：正），$k_v$ は鉛直方向加速度の重力加速度に対する比率（上向き：正）．これにより物体力（慣性力）による節点力ベクトル ${\boldsymbol{f_{b}}}$ は以下の通り定義される．

\begin{equation}
\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot t\cdot \Delta}{3}
\begin{bmatrix}
1/2 & 0   & 1/4 & 0   & 1/4 & 0   \\
0   & 1/2 & 0   & 1/4 & 0   & 1/4 \\
1/4 & 0   & 1/2 & 0   & 1/4 & 0   \\
0   & 1/4 & 0   & 1/2 & 0   & 1/4 \\
1/4 & 0   & 1/4 & 0   & 1/2 & 0   \\
0   & 1/4 & 0   & 1/4 & 0   & 1/2
\end{bmatrix}
\begin{Bmatrix} k_h \\ k_v \\ k_h \\ k_v \\ k_h \\\ k_v \end{Bmatrix}
\end{equation}

ここで，$\Delta$ は要素面積，$t$ は板厚，$\gamma$ は要素の単位体積重量である．

節点外力ベクトル

節点外力ベクトルについては，ここで紹介するプログラムでは，直接，等価節点外力を入力データファイルから入力することとしており，特別な計算処理は行っていない．

要素応力ベクトル

要素応力ベクトル ${\boldsymbol{\sigma}}$ は，以下のとおり求めることができる．

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} & &\text{（節点変位によるひずみ）}\\
&T_m=\cfrac{\phi_i+\phi_j+\phi_k}{3} & & \text{（要素平均温度変化量）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_m \\ (1+\nu)\cdot \alpha\cdot T_m \\ 0 \end{Bmatrix}
& & \text{（平面ひずみ）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_m \\ \alpha\cdot T_m \\ 0 \end{Bmatrix}
& & \text{（平面応力）}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
& & \text{（要素応力ベクトル）}
\end{align}

三角形一次要素を用いた２次元応力解析プログラム

プログラムの概要

プログラムの概要を以下に示す．特徴としては，強制変位（ゼロ以外の節点変位指定）を扱えること，連立一次方程式の解法に疎行列計算ライブラリの scipy.sparse.linalg.spsolve を用いていることであろうか．

• プログラミング言語はPythonとする
• 等方弾性体の２次元応力解析を行う
• 平面ひずみおよび平面応力問題を扱える
• 要素は，１要素３節点，１節点２自由度の三角形一次要素とする
• 座標システムは，平面上の直行座標系右方向を x 軸正方向，上方向を y 軸正方向とする．
• 荷重として，節点外力，節点強制変位（ゼロ変位含む），節点温度変化，要素への慣性力を扱える
• 連立一次方程式は，その大きさ（１節点自由度ｘ総節点数）を変えずに，節点強制変位導入による処理を行い解いている．
• 連立一次方程式は，大きな自由度での計算効率を考慮し，疎行列計算ライブラリの scipy.sparse.linalg.spsolve(A,b) で解く
• 出力は，節点変位および要素平均応力とする
• 入力データファイル，出力データファイルの文字区切りは空白とする
• 入力データファイル名，出力データファイル名は，コマンドライン引数として入力する

扱える境界条件の説明

項目	説明
節点外力	載荷節点数，積荷節点番号，荷重値を指定
節点温度変化	全節点の温度変化量を指定．温度上昇を正値とする
要素慣性力	材料特性データの中で水平方向および鉛直方向の加速度を重力加速度に対する比で指定
節点強制変位	強制変位を与える節点数，節点番号，変位値を指定

プログラムの構成（プログラム名：py_fem_pl3x.py）

1. モジュール読み込み
2. 関数：データ入力 (def inpdata_pl3)
   (1) 基本数値読み込み
   (2) 配列宣言
   (3) 材料特性データ読み込み
   (4) 要素構成節点番号および材料番号読み込み
   (5) 節点座標および節点温度変化量読み込み
   (6) 境界条件読み込み
   (7) 節点荷重データ読み込み
3. 関数：入力データ書き出し (def prinp_pl3)
4. 関数：計算結果書き出し (def prout_pl3)
5. 関数：Dマトリックス作成 (def dmat_pl3)
6. 関数：Bマトリックス作成 (def bmat_pl3)
7. 関数：要素剛性マトリックス作成 (def sm_pl3)
8. 関数：要素応力計算 (def calst_pl3)
9. 関数：主応力計算 (def pst_pl3)
10. 関数：節点温度荷重ベクトル作成 (def tfvec_pl3)
11. 関数：節点慣性力ベクトル作成 (def bfvec_pl3)
12. 関数：メインルーチン
    (1) コマンドラインより入出力ファイル名取得
    (2) データ入力
    (3) 入力データ書き出し
    (4) 剛性マトリックス組み立て
    (5) 境界条件処理（変位既知自由度の処理）
    (6) 剛性方程式を解く（変位計算）
    (7) 境界条件処理（既知変位の組み込み）
    (8) 要素応力計算
    (9) 計算結果書き出し
    (10) 計算時間等情報表示
13. メインルーチン実行

実行コマンドと入出力データ

解析実行コマンド

python3 py_fem_pl3x.py inp.txt out.txt

変数	説明
py_fem_pl3x.py	Python による FEM プログラム名
inp.txt	入力データファイル名（空白区切りデータ）
out.txt	出力データファイル名（空白区切りデータ）

入力データファイルフォーマット

npoin  nele  nsec  npfix  nlod  nstr  # Basic values for analysis
t  E  po  alpha  gamma  gkh  gkv      # Material properties
    ..... (1 to nsec) .....           #
node1  node2  node3   isec            # Element connectivity, material set number
    ..... (1 to nele) .....           #   (counter-clockwise order of node number)
x  y  deltaT                          # Coordinate, Temperature change of node
    ..... (1 to npoin) .....          #
node  kox  koy  rdisx  rdisy          # Restricted node and restrict conditions
    ..... (1 to npfix) .....          #
node  fx  fy                          # Loaded node and loading conditions
    ..... (1 to nlod) .....           #   (omit data input if nlod=0)

変数	説明
npoin, nele, nsec	節点数，要素数，材料特性数
npfix, nlod	拘束節点数，載荷節点数
nstr	応力状態（平面歪: 0，平面応力: 1）
t, E, po, alpha	板厚，弾性係数，ポアソン比，線膨張係数
gamma, gkh, gkv	単位体積重量，水平及び鉛直方向加速度（gの比）
node1, node2, node3, isec	節点1，節点2，節点3，材料特性番号
x, y, deltaT	節点x座標，節点y座標，節点温度変化
node, kox, koy	拘束節点番号，x及びy方向拘束の有無（拘束: 1，自由: 0）
rdisx, rdisy	x及びy方向変位（無拘束でも0を入力）
node, fx, fy	載荷重節点番号，x方向荷重，y方向荷重

出力データファイルフォーマット

npoin  nele  nsec npfix  nlod   nstr
    (Each value of above)
sec  t  E  po  alpha  gamma  gkh  gkv
    sec   : Material number
    t     : Element thickness
    E     : Elastic modulus
    po    : Poisson's ratio
    alpha : Thermal expansion coefficient
    gamma : Unit weight
    gkh   : Acceleration in x-direction (ratio to g)
    gkv   : Acceleration in y-direction (ratio to g)
    ..... (1 to nsec) .....
node  x  y  fx  fy  deltaT  kox  koy
    node   : Node number
    x      : x-coordinate
    y      : y-coordinate
    fx     : Load in x-direction
    fy     : Load in y-direction
    deltaT : Temperature change of node
    kox    : Index of restriction in x-direction (0: free, 1: fixed)
    koy    : Index of restriction in x-direction (0: free, 1: fixed)
    ..... (1 to npoin) .....
node  kox  koy  rdis_x  rdis_y
    node    : Node number
    kox     : Index of restriction in x-direction (0: free, 1: fixed)
    koy     : Index of restriction in y-direction (0: free, 1: fixed)
    rdis_x  : Given displacement in x-direction
    tdis_y  : Given displacement in y-direction
    ..... (1 to npfix) .....
elem  i  j  k  sec
    elem : Element number
    i    : Node number of start point
    j    : Node number of second point
    k    : Node number of third point
    sec  : Material number
    ..... (1 to nele) .....
node  dis-x  dis-y
    node  : Node number
    dis-x : Displacement in x-direction
    dis-y : Displacement in y-direction
    ..... (1 to npoin) .....
elem  sig_x  sig_y  tau_xy  p1  p2  ang
    elem   : Element number
    sig_x  : Normal stress in x-direction
    sig_y  : Normal stress in y-direction
    tau_xy : Shear stress in x-y plane
    p1     : First principal stress
    p2     : Second principal stress
    ang    : Angle of the first principal stress
    ..... (1 to nele) .....
n=(total degrees of freedom)  time=(calculation time)

変数	説明
dis-x, dis-y	x方向変位，y方向変位
sig_x, sig_y, tau_xy	x方向直応力，y方向直応力，せん断応力
p1, p2, ang	第一主応力，第二主応力，第一主応力の方向

プログラミング

注意事項

節点番号と配列のインデックス

• 入力データとしての節点番号は，慣例に習って１から番号付けを行っている．
• Pythonプログラム内での配列インデックスは0から始まっており，データの格納もインデックス0から行っている．
• このため，プログラム内での節点番号を示すインデックスは，入力データの節点番号から1を減じた値となっていることに注意する．

応力の符号

• 応力の符号は，材料力学の慣例に従い，引張を正，圧縮を負とする．

自重

• 座標系は，水平右向きをx軸正方向，鉛直上向きをy軸正方向と設定している．
• このため，自重を扱い場合は，材料特性にて，y方向の加速度（重力加速度との比）に gkv＝−１を入力する．

平面ひずみと平面応力

• 平面ひずみ解析と平面応力解析の切り替えは，入力データ１行目の変数 nstr で行う．
• 平面ひずみ：nstr＝０，平面応力：nstr＝１ である．
• 平面ひずみ解析の場合でも，要素板厚は単位（＝１）を入力しておく．

1. モジュール読み込み

• numpy：配列演算用
• scipy.sparse.linalg.spsolve：疎行列連立一次方程式解法
• scipy.sparse import csr_matrix：疎行列圧縮格納
• sys：コマンドライン引数使用
• time：計算時間計測用

#===========================
# 2DFEM Stress Analysis
#===========================
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time

2. 関数：データ入力 (def inpdata_pl3)

(1) 基本数値読み込み（npoin. nele, nsec, npfix, nlod, nstr）
(2) 配列宣言（ae, node, x, deltaT, mpfix, rdis, fp）
(3) 材料特性データ読み込み（配列 ae）
(4) 要素構成節点番号および材料番号読み込み（配列 node）
(5) 節点座標および節点温度変化量読み込み（配列 x, deltaT）
(6) 境界条件読み込み（配列 mpfix, rdis）
(7) 節点荷重データ読み込み（配列 fp）

def inpdata_pl3(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin =int(text[0]) # Number of nodes
    nele  =int(text[1]) # Number of elements
    nsec  =int(text[2]) # Number of sections
    npfix =int(text[3]) # Number of restricted nodes
    nlod  =int(text[4]) # Number of loaded nodes
    nstr  =int(text[5]) # 0: plane strain, 1: plane stress
    # array declaration
    ae    =np.zeros((7,nsec),dtype=np.float64)      # Section characteristics
    node  =np.zeros((nod+1,nele),dtype=np.int)      # Node-element relationship
    x     =np.zeros((2,npoin),dtype=np.float64)     # Coordinates of nodes
    deltaT=np.zeros(npoin,dtype=np.float64)         # Temperature change of node
    mpfix =np.zeros((nfree,npoin),dtype=np.int)     # Ristrict conditions
    rdis  =np.zeros((nfree,npoin),dtype=np.float64) # Ristricted displacement
    fp    =np.zeros(nfree*npoin,dtype=np.float64)   # External force vector
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i]=float(text[0]) #t    : :Plate thickness
        ae[1,i]=float(text[1]) #E    : Elastic modulus
        ae[2,i]=float(text[2]) #po   : Poisson's ratio
        ae[3,i]=float(text[3]) #alpha: Thermal expansion coefficient
        ae[4,i]=float(text[4]) #gamma: Unit weight of material
        ae[5,i]=float(text[5]) #gkh  : Acceleration in x-direction
        ae[6,i]=float(text[6]) #gkv  : Acceleration in y-direction
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #node_3
        node[3,i]=int(text[3]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]   =float(text[0]) # x-coordinate
        x[1,i]   =float(text[1]) # y-coordinate
        deltaT[i]=float(text[2]) # Temperature change of node
    # boundary conditions (0:free, 1:restricted)
    for i in range(0,npfix):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])             #fixed node
        mpfix[0,lp-1]=int(text[1])  #fixed in x-direction
        mpfix[1,lp-1]=int(text[2])  #fixed in y-direction
        rdis[0,lp-1]=float(text[3]) #fixed displacement in x-direction
        rdis[1,lp-1]=float(text[4]) #fixed displacement in y-direction
    # load
    if 0<nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[2*lp-2]=float(text[1]) #load in x-direction
            fp[2*lp-1]=float(text[2]) #load in y-direction
    f.close()
    return npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp

3. 関数：入力データ書き出し (def prinp_pl3)

def prinp_pl3(fnameW,npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>6s}'
    .format('npoin','nele','nsec','npfix','nlod','nstr'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:6d}'
    .format(npoin,nele,nsec,npfix,nlod,nstr),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>10s} {7:>10s}'
    .format('sec','t','E','po','alpha','gamma','gkh','gkv'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:10.3f} {7:10.3f}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i],ae[6,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s}'
    .format('node','x','y','fx','fy','deltaT','kox','koy'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d}'
        .format(lp,x[0,i],x[1,i],fp[2*i],fp[2*i+1],deltaT[i],mpfix[0,i],mpfix[1,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>15s} {4:>15s}'
    .format('node','kox','koy','rdis_x','rdis_y'),file=fout)
    for i in range(0,npoin):
        if 0<mpfix[0,i]+mpfix[1,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:15.7e} {4:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],rdis[0,i],rdis[1,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'
    .format('elem','i','j','k','l','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'
        .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne]),file=fout)
    fout.close()

4. 関数：計算結果書き出し (def prout_pl3)

def prout_pl3(fnameW,npoin,nele,disg,strs):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s}'.format('node','dis-x','dis-y'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e}'.format(lp,disg[2*lp-2],disg[2*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('elem','sig_x','sig_y','tau_xy','p1','p2','ang'),file=fout)
    for ne in range(0,nele):
        sigx =strs[0,ne]
        sigy =strs[1,ne]
        tauxy=strs[2,ne]
        ps1,ps2,ang=pst_pl3(sigx,sigy,tauxy)
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(ne+1,sigx,sigy,tauxy,ps1,ps2,ang),file=fout)
    fout.close()

5. 関数：Dマトリックス作成 (def dmat_pl3)

\begin{align}
&\text{平面ひずみ（nstr=0）：}& &[\boldsymbol{D_e}]=\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \cfrac{1-2\nu}{2} \end{bmatrix}\\
&\text{平面応力（nstr=1）：}& &[\boldsymbol{D_e}]=\frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix} \qquad
\begin{matrix} E & \text{：弾性係数} \\ \nu & \text{：ポアソン比} \end{matrix}
\end{align}

def dmat_pl3(nstr,E,po):
    d=np.zeros((3,3),dtype=np.float64)
    if nstr==0: #plane strain
        d[0,0]=1-po; d[0,1]=po
        d[1,0]=po  ; d[1,1]=1-po
        d[2,2]=0.5*(1-2*po)
        d=E/(1+po)/(1-2*po)*d
    else: # plane stress
        d[0,0]=1  ; d[0,1]=po
        d[1,0]=po ; d[1,1]=1
        d[2,2]=0.5*(1-po)
        d=E/(1-po**2)*d
    return d

6. 関数：Bマトリックス作成 (def bmat_pl3)

\begin{align}
&\Delta=1/2\cdot[(x_k-x_j)y_i+(x_i-x_k)y_j+(x_j-x_i)y_k] \qquad \text{（要素面積）}\\
&[\boldsymbol{B}]
=\frac{1}{2\Delta}
\begin{bmatrix}
y_j-y_k & 0       & y_k-y_i & 0       & y_i-y_j & 0       \\
0       & x_k-x_j & 0       & x_i-x_k & 0       & x_j-x_i \\
x_k-x_j & y_j-y_k & x_i-x_k & y_k-y_i & x_j-x_i & y_i-y_j
\end{bmatrix}
\end{align}

def bmat_pl3(x1,y1,x2,y2,x3,y3):
    bm=np.zeros((3,6),dtype=np.float64)
    area=0.5*((x3-x2)*y1+(x1-x3)*y2+(x2-x1)*y3)
    bm[0,0]=y2-y3; bm[0,1]=0    ; bm[0,2]=y3-y1; bm[0,3]=0    ; bm[0,4]=y1-y2; bm[0,5]=0
    bm[1,0]=0    ; bm[1,1]=x3-x2; bm[1,2]=0    ; bm[1,3]=x1-x3; bm[1,4]=0    ; bm[1,5]=x2-x1
    bm[2,0]=x3-x2; bm[2,1]=y2-y3; bm[2,2]=x1-x3; bm[2,3]=y3-y1; bm[2,4]=x2-x1; bm[2,5]=y1-y2
    bm=bm/2/area
    return bm,area

7. 関数：要素剛性マトリックス作成 (def sm_pl3)

\begin{equation}
[\boldsymbol{k}]=t\cdot\int_A[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dA=t\cdot\Delta\cdot[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]
\end{equation}

def sm_pl3(nstr,t,E,po,x1,y1,x2,y2,x3,y3):
    #Stiffness matrix [sm]=[B]T[D][B]*t*det(J)
    d=dmat_pl3(nstr,E,po)
    bm,area=bmat_pl3(x1,y1,x2,y2,x3,y3)
    sm=np.dot(bm.T,np.dot(d,bm))*t*area
    return sm

8. 関数：要素応力計算 (def calst_pl3)

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} & &\text{（節点変位によるひずみ）}\\
&T_m=\cfrac{\phi_i+\phi_j+\phi_k}{3} & & \text{（要素平均温度変化量）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_m \\ (1+\nu)\cdot \alpha\cdot T_m \\ 0 \end{Bmatrix}
& & \text{（平面ひずみ：nstr=0）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_m \\ \alpha\cdot T_m \\ 0 \end{Bmatrix}
& & \text{（平面応力：nstr=1）}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
& & \text{（要素応力ベクトル）}
\end{align}

def calst_pl3(nstr,E,po,alpha,tem,wd,x1,y1,x2,y2,x3,y3):
    eps0=np.zeros(3,dtype=np.float64)
    #stress vector {stress}=[D][B]{u}
    d=dmat_pl3(nstr,E,po)
    bm,area=bmat_pl3(x1,y1,x2,y2,x3,y3)
    eps=np.dot(bm,wd)
    #Thermal strain
    if nstr==0: # plane strain
        eps0[0]=tem*(1.0+po)*alpha
        eps0[1]=eps0[0]
        eps0[2]=0.0
    else: # plane stress
        eps0[0]=tem*alpha
        eps0[1]=eps0[0]
        eps0[2]=0.0
    stress=np.dot(d,(eps-eps0))
    return stress

9. 関数：主応力計算 (def pst_pl3)

主応力は次式で計算する．なお応力の符号は引張が正である．

\begin{align}
&\text{第一主応力} & &\sigma_1=\cfrac{\sigma_x+\sigma_y}{2}+\sqrt{\left(\cfrac{\sigma_x-\sigma_y}{2}\right)^2+\tau_{xy}{}^2} \\
&\text{第二主応力} & &\sigma_2=\cfrac{\sigma_x-\sigma_y}{2}-\sqrt{\left(\cfrac{\sigma_x-\sigma_y}{2}\right)^2+\tau_{xy}{}^2}
\end{align}

第一主応力の方向（x軸から反時計回りの角度）は，一般には次式で計算する．場合分けについてはプログラムを参照．

\begin{equation}
\theta=\cfrac{1}{2}\tan^{-1}\left(\cfrac{2 \tau_{xy}}{\sigma_x-\sigma_y}\right)
\end{equation}

def pst_pl3(sigx,sigy,tauxy):
    ps1=0.5*(sigx+sigy)+np.sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
    ps2=0.5*(sigx+sigy)-np.sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
    if sigx==sigy:
        if tauxy >0.0: ang= 45.0
        if tauxy <0.0: ang=135.0
        if tauxy==0.0: ang=  0.0
    else:
        ang=0.5*np.arctan(2.0*tauxy/(sigx-sigy))
        ang=180.0/np.pi*ang
        if sigx>sigy and tauxy>=0.0: ang=ang
        if sigx>sigy and tauxy <0.0: ang=ang+180.0
        if sigx<sigy:                ang=ang+90.0
    return ps1,ps2,ang

10. 関数：節点温度荷重ベクトル作成 (def tfvec_pl3)

\begin{align}
&T_m=\cfrac{\phi_i+\phi_j+\phi_k}{3} & & \text{（要素平均温度変化量）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_m \\ (1+\nu)\cdot \alpha\cdot T_m \\ 0 \end{Bmatrix}
& & \text{（平面ひずみ：nstr=0）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_m \\ \alpha\cdot T_m \\ 0 \end{Bmatrix}
& & \text{（平面応力：nstr=1）}\\
&\{\boldsymbol{f_t}\}=t\cdot\int_A[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dA=t\cdot\Delta\cdot[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}
& & \text{（節点温度荷重ベクトル）}
\end{align}

def tfvec_pl3(nstr,t,E,po,alpha,tem,x1,y1,x2,y2,x3,y3):
    eps0=np.zeros(3,dtype=np.float64)
    # {tfe=[B]T[D]{eps0}
    d=dmat_pl3(nstr,E,po)
    bm,area=bmat_pl3(x1,y1,x2,y2,x3,y3)
    #Thermal strain
    if nstr==0: # plane strain
        eps0[0]=tem*(1.0+po)*alpha
        eps0[1]=eps0[0]
        eps0[2]=0.0
    else: # plane stress
        eps0[0]=tem*alpha
        eps0[1]=eps0[0]
        eps0[2]=0.0
    tfe=np.dot(bm.T,np.dot(d,eps0))*t*area
    return tfe

11. 関数：節点慣性力ベクトル作成 (def bfvec_pl3)

\begin{equation}
\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot t\cdot \Delta}{3}
\begin{bmatrix}
1/2 & 0   & 1/4 & 0   & 1/4 & 0   \\
0   & 1/2 & 0   & 1/4 & 0   & 1/4 \\
1/4 & 0   & 1/2 & 0   & 1/4 & 0   \\
0   & 1/4 & 0   & 1/2 & 0   & 1/4 \\
1/4 & 0   & 1/4 & 0   & 1/2 & 0   \\
0   & 1/4 & 0   & 1/4 & 0   & 1/2
\end{bmatrix}
\begin{Bmatrix} k_h \\ k_v \\ k_h \\ k_v \\ k_h \\\ k_v \end{Bmatrix}
\end{equation}

def bfvec_pl3(t,gamma,gkh,gkv,x1,y1,x2,y2,x3,y3):
    mat=np.zeros((6,6),dtype=np.float64)
    area=0.5*((x3-x2)*y1+(x1-x3)*y2+(x2-x1)*y3)
    mat[0,0]=0.50;mat[0,1]=0.00;mat[0,2]=0.25;mat[0,3]=0.00;mat[0,4]=0.25;mat[0,5]=0.00
    mat[1,0]=0.00;mat[1,1]=0.50;mat[1,2]=0.00;mat[1,3]=0.25;mat[1,4]=0.00;mat[1,5]=0.25
    mat[2,0]=0.25;mat[2,1]=0.00;mat[2,2]=0.50;mat[2,3]=0.00;mat[2,4]=0.25;mat[2,5]=0.00
    mat[3,0]=0.00;mat[3,1]=0.25;mat[3,2]=0.00;mat[3,3]=0.50;mat[3,4]=0.00;mat[3,5]=0.25
    mat[4,0]=0.25;mat[4,1]=0.00;mat[4,2]=0.25;mat[4,3]=0.00;mat[4,4]=0.50;mat[4,5]=0.00
    mat[5,0]=0.00;mat[5,1]=0.25;mat[5,2]=0.00;mat[5,3]=0.25;mat[5,4]=0.00;mat[5,5]=0.50
    mat=mat*gamma*t*area/3
    w=np.array([gkh,gkv,gkh,gkv,gkh,gkv],dtype=np.float64)
    bfe=np.dot(mat,w)
    return bfe

12. 関数：メインルーチン

メインルーチン処理手順

(1) コマンドラインより入出力ファイル名取得

コマンドライン引数として入出力データファイル名を取得．
入力データファイル名：fnameR，出力データファイル名：fnameW

(2) データ入力

関数 inpdata_pl3 により入力データ取得

(3) 入力データ書き出し

関数 prinp_pl3 により入力データ書き出し

(4) 剛性マトリックス組み立て

各要素毎に，剛性マトリックスsm，温度荷重ベクトル tfe，慣性力ベクトル bfe を計算し，全体剛性方程式に組み込んでいく．ここで剛性マトリックス，温度荷重ベクトル，慣性力ベクトル計算の関数の引数は，全てスカラーとした．要素番号が指定されれば，要素を構成する節点座標や要素物性値は特定できるため，他の要素の値を含んで収納する配列よりは，スカラーを送り込んだほうがプログラム的にスッキリするように思う．
全体剛性マトリックスは変数 gk に，荷重ベクトルは変数 fp に加算されていく．
要素剛性マトリックスを全体剛性マトリックスに組み込むための位置情報は，配列 ir[ ] に格納されている．

(5) 境界条件処理（変位既知自由度の処理）

前述「連立一次方程式における既知量と未知量の取扱」で述べた処理が行われる．

変数	説明
npoin	総節点数
nfree	１節点自由度（ここでは 2）
mpfix	各節点毎の拘束条件を格納した配列（=1：拘束，=0：自由）
fp	節点荷重ベクトル
rdis	変位拘束する節点の変位量
gk	全体剛性マトリックス
disg	全体剛性方程式の解（変位）

(6) 剛性方程式を解く（変位計算）

csr_matrix により剛性マトリックスを疎行列として圧縮格納した後，spsolve にて解いている．

(7) 境界条件処理（既知変位の組み込み）

連立一次方程式を解いた後，既知変位に相当する自由度の項に，既知変位を入れ直すことを忘れないようにする．

(8) 要素応力計算

関数 calst_pl3 による要素毎の応力計算実施

(9) 計算結果書き出し

関数 prout_pl3 による計算結果出力

(10) 計算時間等情報表示

以下により，計算完了時刻とスタート時刻の差分を処理時間として出力する．

import time

start=time.time()
......
dtime=time.time()-start

メインルーチンプログラム

def main_pl3():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=3          # Number of nodes per element
    nfree=2        # Degree of freedom per node
    # data input
    npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_pl3(fnameR,nod,nfree)
    # print out of input data
    prinp_pl3(fnameW,npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp)
    # array declaration
    ir=np.zeros(nod*nfree,dtype=np.int)         # Work vector for matrix assembly
    gk=np.zeros((nfree*npoin,nfree*npoin),dtype=np.float64)   # Global stiffness matrix
    # assembly of stiffness matrix & load vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        m=node[3,ne]-1
        x1=x[0,i]; y1=x[1,i]
        x2=x[0,j]; y2=x[1,j]
        x3=x[0,k]; y3=x[1,k]
        t    =ae[0,m]
        E    =ae[1,m]
        po   =ae[2,m]
        alpha=ae[3,m]
        gamma=ae[4,m]
        gkh  =ae[5,m]
        gkv  =ae[6,m]
        tem=(deltaT[i]+deltaT[j]+deltaT[k])/3   # average temperature change
        sm=sm_pl3(nstr,t,E,po,x1,y1,x2,y2,x3,y3)  # element stiffness matrix
        tfe=tfvec_pl3(nstr,t,E,po,alpha,tem,x1,y1,x2,y2,x3,y3) # thermal load vector
        bfe=bfvec_pl3(t,gamma,gkh,gkv,x1,y1,x2,y2,x3,y3) # body force vector
        ir[5]=2*k+1; ir[4]=ir[5]-1
        ir[3]=2*j+1; ir[2]=ir[3]-1
        ir[1]=2*i+1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+tfe[i]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+sm[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    #disg = np.linalg.solve(gk, fp)
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                disg[iz]=rdis[j,i]
    # calculation of element stress
    strs=np.zeros((3,nele),dtype=np.float64) # Section force vector
    wd  =np.zeros(6,dtype=np.float64)        # element displacement
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        m=node[3,ne]-1
        x1=x[0,i]; y1=x[1,i]
        x2=x[0,j]; y2=x[1,j]
        x3=x[0,k]; y3=x[1,k]
        wd[0]=disg[2*i]; wd[1]=disg[2*i+1]
        wd[2]=disg[2*j]; wd[3]=disg[2*j+1]
        wd[4]=disg[2*k]; wd[5]=disg[2*k+1]
        E    =ae[1,m]
        po   =ae[2,m]
        alpha=ae[3,m]
        tem=(deltaT[i]+deltaT[j]+deltaT[k])/3
        strs[:,ne]=calst_pl3(nstr,E,po,alpha,tem,wd,x1,y1,x2,y2,x3,y3)
    prout_pl3(fnameW,npoin,nele,disg,strs)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()

13. メインルーチン実行

#==============
# Execution
#==============
if __name__ == '__main__': main_pl3()

（付録）Sympyによる行列演算の確認

Sympyを用いて，$2 \Delta$ および $2 \Delta [\boldsymbol{B}]$ の値を記号演算により確認してみよう．

まず，$2\Delta$ の値を確認してみる．

from sympy import *
init_printing()

x_i, y_i, x_j, y_j, x_k, y_k = symbols('x_i y_i x_j y_j x_k y_k')

D= Matrix([
[1,x_i,y_i],
[1,x_j,y_j],
[1,x_k,y_k],
])
D.det()

上記コードを実行すると，以下の結果が得られるはず．

\begin{equation}
x_{i} y_{j} - x_{i} y_{k} - x_{j} y_{i} + x_{j} y_{k} + x_{k} y_{i} - x_{k} y_{j}
\end{equation}

また，つづいて以下のコードを実行する．これは $2\Delta [\boldsymbol{B}]$ の値を計算するものである．

Xnd= Matrix([
[1,x_i,y_i,0,0,0],
[0,0,0,1,x_i,y_i],
[1,x_j,y_j,0,0,0],
[0,0,0,1,x_j,y_j],
[1,x_k,y_k,0,0,0],
[0,0,0,1,x_k,y_k]
])
C = Xnd.inv()*D.det()
simplify(C)   

上記コードにより，以下の結果が得られる．

\begin{equation}
\left[\begin{matrix}x_{j} y_{k} - x_{k} y_{j} & 0 & - x_{i} y_{k} + x_{k} y_{i} & 0 & x_{i} y_{j} - x_{j} y_{i} & 0\\y_{j} - y_{k} & 0 & - y_{i} + y_{k} & 0 & y_{i} - y_{j} & 0\\- x_{j} + x_{k} & 0 & x_{i} - x_{k} & 0 & - x_{i} + x_{j} & 0\\0 & x_{j} y_{k} - x_{k} y_{j} & 0 & - x_{i} y_{k} + x_{k} y_{i} & 0 & x_{i} y_{j} - x_{j} y_{i}\\0 & y_{j} - y_{k} & 0 & - y_{i} + y_{k} & 0 & y_{i} - y_{j}\\0 & - x_{j} + x_{k} & 0 & x_{i} - x_{k} & 0 & - x_{i} + x_{j}\end{matrix}\right]
\end{equation}

すなわち，$\Delta$ および $[\boldsymbol{B}]$ の値は確認されたことになる．

（余談）筆者のFEM歴など

• 学生時代，橋梁研究室に所属し，立体骨組の弾塑性有限変位解析で修論を書いた．
• 新入社員時代，現場に配属になり，ダムの２次元応力解析をやっていた．プログラム本体は本社のコンピュータに入っており，大きなグラフ用紙にダムの絵を書き自分でメッシュを切り，節点座標，要素ー節点関係を書き込み，それを紙テープ穿孔して本社のコンピュータに送り，戻ってきたデータをラインプリンタ出力するようなことをやっていた．骨組解析のプログラムはなかったので，学生時代に使っていた参考書などをネタに，Basic で自作して，トンネルの覆工設計などに使っていた．
• その後，転勤などしながら，２次元応力解析，熱伝導解析，浸透流解析などのプログラムを勉強と実益を兼ねて作ってきた．
• 言語も N88-Basic から Visual Basic，C++Builder，GNU-C, GNU-gfortan などを経由して，現在はたまーに Fortran 90 を使うが，ほぼ Python １本となっている．Python 便利！
• 現在は祖国日本より 6000km 離れたマレーシアの地で，発電所建屋の構造設計を実施中．でも残念ながら今年（2018年）10月には日本に戻らなければならない．

以　上