1. はじめに

久しぶりに直営で2次元応力解析をやってみたので，その手順を残しておきます．

ソースコードはGistにリンクさせて閲覧可能にしてありますが，その目的はソースを示すことよりも，ちゃんと自分でプログラムを持っているよ，ということを示すためです．ソースよりも，こんな画像を出力しているんだ，というようなことが，訪問者の参考になれば幸いです．

FortranとPythonをごちゃ混ぜで使っていますが，大きな意味での使い分けとしては以下の様な感じです．

汎用プログラムはFortran．基本的に機能を変更するとき以外は書き換えない．
モデル毎に書き換えるものはPython．

もし，プログラムの使い方に興味のある方は，以下のWebページをご覧になってみてください．

http://kbtkt.web.fc2.com//sub_fem_mesh.html

[http://kbtkt.web.fc2.com//sue_fem_plnt.html]
(http://kbtkt.web.fc2.com//sue_fem_plnt.html)

使用プログラム類はほとんど自作ですが，自動メッシュ作成および２次元有限要素応力解析プログラムは，以下の教科書をベースとしています．

(1) 自動メッシュプログラム

教科書：FEMのための要素自動分割　デローニー三角分割法の利用，谷口健男著，2006年8月10日第1版第7刷発行，森北出版株式会社

教科書中の四角形メッシュ分割のプログラムを3分割し，Fortran90に書き換えています．

(2) ２次元有限要素法応力解析プログラム

教科書：コンピュータによる構造工学講座　有限要素法による構造解析プログラム，三本木茂夫・吉村信敏共著，昭和54年5月10日初版第9刷発行，株式会社培風館

教科書中の2次元弾性問題のプログラム(三角形一次要素)に対し，三角形一次要素と四角形アイソパラメトリック一次要素が使えるよう書き換えています．両者の混合モデルは不可．また連立一次方程式は，バンドマトリックス用のコレスキー法で解いています．

2. 使用プログラム

使用プログラムは以下のとおり．

(1) モデル表示

py_gmt_shape.py
：モデル構造の形状をGMTにより画像を作成，表示します．
やっていることは，

GMTに読み込ませるための座標をファイル出力させる
作図用GMTコマンドをファイル出力
bashでGMTコマンドを実行しeps画像を出力させる

報告書用に作成．入力データは自動メッシュ作成用のインプットデータ．プログラムとしての汎用性はありませんが，同一目的の画像出力のため，モデルが変わるごとに書き換えながら使い回すプログラムです．

py_gmt_model.py
：メッシュ図に荷重や境界条件を書き込んだ出力画像を作ります．
やっていることは，

FEM解析用のインプットデータを読み込む
メッシュ図や境界条件・荷重を出力させるデータを作成・出力
作図用GMTコマンドをファイル出力
bashでGMTコマンドを実行しeps画像を出力させる

(2) 自動メッシュ作成

自動メッシュを作成するための元データは，事前にエクセルやテキストエディタなどで作成しておきます．

f90_fem_mesh4a.f90
：三角形要素生成

f90_fem_mesh4b.f90
：四角形要素生成

f90_fem_mesh4c.f90
：四角形要素整形

f90_num_mesh.f90
：節点番号最適化．出力されたメッシュデータを，コレスキー法用に，バンド幅を縮小すべく，節点番号を書き換える．

f90_gmt_mesh.f90
：要素図をGMTにより画像出力する．下のスクリプトでは，GMTにより，fig_gmt_mesh.eps という画像ファイルが作成され，要素分割状況を確認します．

上記プログラムは，以下のスクリプトで実行する．

gfortran -o f90_fem_mesh4a f90_fem_mesh4a.f90
gfortran -o f90_fem_mesh4b f90_fem_mesh4b.f90
gfortran -o f90_fem_mesh4c f90_fem_mesh4c.f90
gfortran -o f90_gmt_mesh f90_gmt_mesh.f90
gfortran -o f90_num_mesh f90_num_mesh.f90

range=-6/6/175/185
scale=12/10
gmt set FONT_ANNOT_PRIMARY 10
gmt set MAP_ANNOT_OFFSET_PRIMARY 0.3c
gmt set FONT_LABEL 10p
gmt set MAP_LABEL_OFFSET 0.1c
gmt set MAP_TICK_LENGTH_PRIMARY -0.2c

inp=mesh_arch4.csv
./f90_fem_mesh4a $inp out_work_A.txt
./f90_fem_mesh4b out_work_A.txt out_work_B.txt
./f90_fem_mesh4c out_work_B.txt out_work_C.txt
./f90_num_mesh 4 out_work_C.txt out_mesh.txt
./f90_gmt_mesh 4 out_mesh.txt _inp_gmt.txt

inp=_inp_gmt.txt
fig=fig_gmt_mesh.eps
gmt psxy $inp -R$range -JX$scale -W0.5,black -G255/255/191 -P > $fig

py_model_dat.py
：メッシュデータ（上記の成果物であるテキストファイル out_mesh.txt）に，要素物性・境界条件・荷重条件を加えて，解析用入力データを作成します．
プログラムに汎用性はありませんが，モデルが変わるごとに書き換えながら使い回すプログラムです．

(3) 応力解析と結果表示

f90_fem_plnt.f90
：FEM応力解析プログラム

py_fem_cont.py
：第一主応力・第二種応力・最大せん断応力のコンター図．コンターといっても連続的に表示するのではなく，応力値に応じて要素に着色したもの．設計で使うにはこちらのほうが便利です．

py_fem_vect.py
：主応力ベクトル図・変位モード図作成．

上記プログラムは以下のスクリプトで実行．計算を実施し，結果を読み込んでpng画像を出力します．

gfortran -o f90_fem_plnt f90_fem_plnt.f90

./f90_fem_plnt inp_arch1.csv out_arch1.csv

python3 py_fem_cont.py out_arch1.csv
python3 py_fem_vect.py out_arch1.csv

この応力コンターを表示するプログラムは，GMTでもトライしてみましたが，
「任意の形状の四角形要素を応力値に応じて塗りつぶす」という作業が，塗りつぶす値別にファイルを作る必要があるなど，GMTでは面倒だったため，Pythonにより配列に覚えさせた座標に囲まれた領域（メッシュ）を応力値に応じて塗りつぶすという方法を採用しました．

なお，GMTではモデル領域を細かい矩形に分割して内挿により各点の応力値を計算し，コンターを作成する方法は，この事例では使っていませんが，うまくいっています．

3. 成果図