01:有限要素法プログラム解説(pdf)
概要( pdf ファイル)
- ここに掲載の有限要素法プログラムは,「三本木茂夫・吉村信敏:コンピュータによる構造工学講座I-1-B 日本鋼構造協会編 有限要素法による構造解析プログラム 考え方と解説,培風館,昭和54年5月」に掲載されている FORTRAN による三角形定ひずみ要素の2次元 FEM 応力解析プログラムおよび平面骨組構造解析プログラムのフローを参考に,各種機能を追加し VB 化あるいは C# 化したものです.
- 四角形要素は 4 節点パラメトリック要素で Gauss 積分点は 4 点,三角形要素は 3 節点定ひずみ要素としています.
- いずれのプログラムも一度全自由度を含む全体剛性行列を作成していますが,これは熱伝導解析の温度固定や応力解析での強制変位導入処理をわかりやすくするために行っているものです.
- 連立一次方程式はバンド化したコレスキー法で解いています.バンド化コレスキー法プログラムは,「宮沢健二:パソコン構造力学--力学挙動の可視化-- 啓学出版株式会社 1984年6月」掲載のものを参考にしました.掃き出し法と比べ格段にスピードアップしました.
- 有限要素法の解説では共通となる部分も多いため,解説書のみを1項目としてまとめました.2次元弾性応力解析でも取り入れている「強制変位導入の考え方」は平面骨組構造解析の解説書に記載しています.
- 本HPの有限要素法プログラムでは,各種配列の最初の数値は全て 1 としています.例えば,VB において座標 x(i) を i=1〜NODT(節点数)まで繰り返すとき,VB の配列宣言は添え字の最大値ですから,配列の大きさは,Redim x(NODT) となります.これに対しC#では,配列の大きさは,確保する個数になるので,double[] x=new double[NODT+1]; としていることに注意して下さい(使用していませんが,x[0] 分の個数を +1 しています).
Form に配置しているコントロール
FEM のプログラムでは Form に配置するコントロールを,以下を基本に統一しています.
- ToolStrip1
- ToolStripButton1
- DisplayStyle=Text
- Text="[Start]"
- OpenFileDialog1
- SaveFileDialog1
- TableLayoutPanel1
- ColumnCount=2
- RowCount=5
- Rows(コレクション)
- Row1〜Row5:サイズの型=Absolute(絶対サイズ)25ピクセル
- Columns(コレクション)
- Column1:サイズの型=Absolute(絶対サイズ)75ピクセル
- Column2:サイズの型=AutoSize(自動調整)
- Dock=Top
- CellBorderStyle=InsetDouble
- Label1〜10:TableLayoutPanel内に配置
- Label11:TableLayoutPanel外に配置
解説書
02:FEM2次元非定常熱伝導解析(四角形要素)
概要(プログラム:vbFEMondoCH, vcsFEMondoCH, gccFEMondo)
- このプログラムは,「三本木茂夫・吉村信敏:コンピュータによる構造工学講座I-1-B 日本鋼構造協会編 有限要素法による構造解析プログラム 考え方と解説,培風館,昭和54年5月」に掲載されている FORTRAN による三角形定ひずみ要素の2次元 FEM 応力解析プログラムを参考に,4 節点パラメトリック要素を用いた非定常熱伝導解析プログラムとして作成したものです.
- 挙動の検証結果は pdf 解説に記載しています.扱える境界条件は断熱・温度指定・熱伝達境界の 3 種であり,コンクリートの発熱を模擬した要素発熱を見込むことができます.なお定常解析には対応していません.
- 扱える条件を以下に示します.
断熱境界 | 境界条件として何も入力しなければ断熱境界となる |
温度固定境界 | 温度を指定する節点番号を入力しその温度時刻歴を指定する |
熱伝達境界 | 熱伝達する要素の辺と熱伝達率を指定し境界外部の温度時刻歴を指定する |
要素発熱 | コンクリートの断熱温度上昇形式の要素発熱を扱える |
- 扱える要素発熱は以下の式で表されるものです.
- T = K (1 - exp(-α*t))
- T は断熱温度上昇量,K は最終温度上昇量,α は発熱速度を表すパラメータ,t は時刻
- StatusStripには処理開始時刻(入出力ファイル選定後から)と処理終了時刻(ファイル出力完了後まで)および計算時間を示すようにしました.
- 入力ファイルは,基本情報入力用と境界条件の時刻歴入力用の2種が必要ですが,時刻歴入力用ファイル名は基本情報入力用ファイルの中で指定し,画面操作は入力ファイル指定と出力ファイル指定の 2 回としています.
- プログラムを開始するとまず OpenFileDialog が開きますので入力データファイルを指定します.次に SaveFileDialog が開きますので,出力ファイル名を指定します.計算が終わると MessageBox が計算完了を通知します.
気分的なもの
プログラムの基本部分は,2次元応力解析のものを用い,非定常熱伝導解析用に書き換えています.
熱伝導と応力解析の双方ができれば温度変化が構造物の発生応力に与える影響を算定できるため,それを目的に作成したものです.
入出力事例
- 入出力事例は,機能確認を行うための単純な発熱体あるいは単純な境界条件を有するモデルとしています.
- 解説書には水圧鉄管内面で水温変動を受けた場合の周辺の温度分布計算例を載せました.
- 入出力書式は解説書の項の pdf ファイルを参照してください.
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜12,CheckBox1
プログラムリスト
入出力事例
03:FEM2次元応力解析(四角形要素:No-tension,強制変位考慮可)
概要(プログラム:vbFEM4nodCHTS, vcsFEM4nodCHTS, gccFEM4nodTS)
- このプログラムは,「三本木茂夫・吉村信敏:コンピュータによる構造工学講座I-1-B 日本鋼構造協会編 有限要素法による構造解析プログラム 考え方と解説,培風館,昭和54年5月」に掲載されている FORTRAN による三角形定ひずみ要素の2次元 FEM 応力解析プログラムを参考に,四角形要素での2次元 No-tension 解析プログラムとして Visual Basic 化したものです.
- 考慮できる荷重・境界条件は以下の通りです.
外力 | 節点外力 |
温度変化 | 節点温度変化量 |
慣性力 | 水平・鉛直方向の要素慣性力 |
変位 | 節点強制変位入力可能 |
No-tension | 材料の引張強さを入力し,要素の主応力が引張強さを超えた場合 No-tension 材料として挙動する |
- 挙動の検証結果は pdf 解説に記載しています.
- No-tension 材料を含むケースでは「多重円環理論による圧力水路トンネルの設計」に掲載したプログラムの解析結果との比較を行っています.
- 実行画面の ListBox には No-tension 解析における収束状況を示すようにしました.
- No-tension の判定は,入力データとして各材料の引張り強さを指定し,これを超えたらひび割れが発生するなどして No-tension 状態になるようにしています.十分大きな引張強さを入力すれば通常の弾性計算になります.
- 収束判定は,収束計算過程での変位増分が十分に小さくなった時としています.温度変化や強制変位がある時は,たとえNo-tension状態になっていないとしても最低2回は収束計算を行います.
- プログラムを開始するとまず OpenFileDialog が開きますので入力データファイルを指定します.次に SaveFileDialog が開きますので,出力ファイル名を指定します.計算が終わると MessageBox が計算完了を通知します.
気分的なもの
近年多くの有限要素法プログラムが市販されておりそれらを便利に使うことができるようになりました.
しかし簡単なプログラムを手元に置くことにより,「ただで解析できる」,「自由度が高い(例えば要素番号を複数指定しその平均応力を出力したり,物性値を乱数入力したりすることも補助プログラムを作成することにより対応できます)」などのメリットがあります.このため有限要素法の勉強を兼ねてプログラムを自作してみました.温度応力算定機能および No-tension 機能はやってみたいことがあり追加しました.No-tension といっても線形剛性のまま座標変換により引張応力成分を除去しながら不平衡力が小さくなるまで繰り返し計算をしているものです.
入出力事例
- 入出力事例は,長方形構造物の単純な事例4ケースと岩盤内埋設型水圧鉄管のモデルです.
- 岩盤内埋設型水圧鉄管モデルは,鉄管背面の充填コンクリートの引張強さを 0 として No-tension 状態になるようにしたものです.
- 入出力書式は解説書の項の pdf ファイルを参照してください.
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト
入出力事例
04:FEM軸対称応力解析(四角形要素:No-tension,強制変位考慮可)
概要(プログラム:vbFEM4nodASTS, vcsFEM4nodASTS, gccFEM4nodASTS)
- 四角形要素を用いた軸対称有限要素法の No-tension 解析プログラムです.
- 四角形要素による No-tension 解析プログラムと概ね同等の機能です.
- 剛性方程式の考え方,解析事例は解説書の項の pdf ファイル(pdfNO-TENSION.pdf)を参照してください.
- 入出力書式は解説書の項の pdf ファイル(pdfFEMdata.pdf)を参照してください.
- 考慮できる荷重・境界条件は以下の通りです.
外力 | 節点外力 |
温度変化 | 節点温度変化量 |
変位 | 節点強制変位入力可能 |
No-tension | 材料の引張強さを入力し,要素の主応力が引張強さを超えた場合 No-tension 材料として挙動する |
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト
入出力事例
05:FEM2次元応力解析(三角形要素:No-tension,強制変位考慮可)
概要(プログラム:vbFEM3nodCHTS, gccFEM3nodTS)
- 三角形定ひずみ要素を用いた2次元有限要素法の No-tension 解析プログラムです.
- 四角形要素による No-tension 解析プログラムと同等の機能です.
- 入出力書式は解説書の項の pdf ファイルを参照してください.
- 考慮できる荷重・境界条件は以下の通りです.
外力 | 節点外力 |
温度変化 | 節点温度変化量 |
慣性力 | 水平・鉛直方向の要素慣性力 |
変位 | 節点強制変位入力可能 |
No-tension | 材料の引張強さを入力し,要素の主応力が引張強さを超えた場合 No-tension 材料として挙動する |
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト
06:FEM2次元応力解析(四角形要素:弾性解析,強制変位考慮可)
概要(プログラム:vbFEM4nodCH, vcsFEM4nodCH, gccFEM4nod)
- 四角形要素を用いた2次元有限要素法のプログラムです.
- No-tension 解析機能のない基本的なものです.
- 入出力書式は解説書の項の pdf ファイルを参照してください.
- 強制変位導入の考え方は連立方程式の未知数と既知数の取り扱いであり,プログラムによる違いはないため,平面骨組構造解析の解説書の記載を参考にしてください.
- 考慮可能な荷重・境界条件は以下の通りです.
外力 | 節点外力 |
温度変化 | 節点温度変化量 |
慣性力 | 水平・鉛直方向の要素慣性力 |
変位 | 節点強制変位入力可能 |
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト
07:FEM2次元応力解析(三角形要素:弾性解析,強制変位考慮可)
概要(プログラム:vbFEM3nodCH, vcsFEM3nodCH, gccFEM3nod)
- 三角形定ひずみ要素を用いた2次元有限要素法のプログラムです.
- No-tension 解析機能のない基本的なものです.
- 入出力書式は解説書の項の pdf ファイルを参照してください.
- 強制変位導入の考え方は連立方程式の未知数と既知数の取り扱いであり,プログラムによる違いはないため,平面骨組構造解析の解説書の記載を参考にしてください.
- 考慮可能な荷重・境界条件は以下の通りです.
外力 | 節点外力 |
温度変化 | 節点温度変化量 |
慣性力 | 水平・鉛直方向の要素慣性力 |
変位 | 節点強制変位入力可能 |
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト
08:FEM線形弾性平面骨組構造解析(梁要素:強制変位考慮可)
概要(プログラム:vbFRAME, vcsFRAME, gccFRAME)
- 梁要素を用いた線形弾性の平面骨組構造解析のプログラムです.
- プログラムの基本的フローは2次元応力解析のものと同じです.
- 入出力書式は解説書の項の pdf ファイルを参照してください.
- 考慮可能な荷重・境界条件は以下の通りです.
外力 | 節点外力 |
温度変化 | 節点温度変化量 |
慣性力 | 水平・鉛直方向の要素慣性力 |
変位 | 節点強制変位入力可能 |
入出力事例について
- 「inpbeam1.csv」は片持梁先端に鉛直集中荷重を受けるモデルです.
- 「inpbeam2.csv」では「inpbeam1.csv」の解析結果による載荷点での鉛直変位を入力し解析を行っています.
- 「inpframe3.csv」は内水圧を受ける円形断面の 1/4 部分を取り出したモデルです.
- 「inpframe4.csv」では「inpframe3.csv」の解析結果による載荷点での水平・鉛直変位を入力し解析を行っています.
- いずれのケースでも荷重載荷時と強制変位入力時で誤差の範囲で同一の結果が得られています.
- 「inp-gtest.csv」は水平片持梁に水平慣性力を作用させ発生軸力を見る事例です.
- 「inptempa1.csv」は両端固定梁が一様な温度変化を受ける場合の発生軸力を見る事例です.
- 「pdf_spangler.pdf」はこのプログラムを用いたSpanglerの式の再現解析事例です.
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト
入出力事例
09:FEM平面骨組有限変位解析(梁要素:弾性)
概要(プログラム:vbFRAME_GAL, vcsFRAME_GAL, gccFRAME_GAL)
- 梁要素を用いた平面骨組有限変位解析のプログラムです.
- 幾何学的非線形剛性を含む増分形剛性方程式を解きます.
- 要素内力は剛体回転を除去した変位に線形弾性剛性を乗じて求めています.
- 収束計算は弧長増分法により行っています.解析スタート時はコレスキー法により連立一次方程式を解いていますが,荷重低下が発生するとコレスキー法での計算が不能になりますので掃き出し法に変更して計算を継続します.
- 考慮可能な荷重・境界条件は以下の通りです.
外力 | 節点外力(増分で与える) |
変位 | 節点変位拘束のみ入力可能 |
入出力事例について
- vb(inp_flist1.csvのデータ並び)およびCの実行時パラメータは,「入力ファイル名,出力ファイル名,kprog,coefAB,deltaS0,coefS,maxload」です.
- 上記,整数 kprog は,第一中間点設定方法の指定です.「1」が「方法1」,「2」が「方法2」です.方法1および2の説明は「有限要素法プログラム解説:pdfFRAME_GAL.pdf」に記載しています.
- 上記,係数 coefAB は,スケーリングパラメータを定めるための係数です.
- 上記,係数 deltaS0 は,弧長増分法における初期の弧長です.
- 上記,係数 coefS は,収束性を高めるために,弧長を変更するパラメータです.1〜5位がいいようです.
- 上記,整数 maxload は,荷重載荷ステップ数です.50,100などの整数で与えます.
- 上記パラメータはモデルケースのテストランであった方が良いと思ったため導入しました.より大きな変形を追跡できるようにするため,解析対象によって変更したほうが良いようです.
使用コントロール
ToolStrip1(Button1),OpenFileDialog1
プログラムリスト
入出力事例
10:FEM平面骨組座屈固有値解析(梁要素:弾性)
概要(プログラム:vbFRAME_EV, vcsFRAME_EV, gccFRAME_EV)
- 梁要素を用いた平面骨組座屈固有値解析のプログラムです.
- 座屈荷重は,入力時に指定した初期荷重{F}に解析で求めた固有値 λ を乗じることで定まります.
- 固有値解析は一般化ヤコビ法にて行っています.
- 全自由度に対する固有値および固有ベクトルを求めていますが,出力は3次までの固有値と固有ベクトルとしています.
- 考慮可能な荷重・境界条件は以下の通りです.
外力 | 節点外力(荷重モードを指定するものであり原則小さな数字) |
変位 | 節点変位拘束のみ入力可能 |
入出力事例について
- vb・c#(inp_flist.csvのデータ並び)およびCの実行時パラメータは,「入力ファイル名,出力ファイル名」です.
- 軸力を受ける単純梁の座屈荷重を求める計算での所要時間です(単位は秒).なぜvbで非実用的な時間を有するかは不明.vb・c#はいずれもリリースビルドしたexeによるものです.
要素数 | 2 | 4 | 10 | 20 | 40 |
VB | 0.03 | 0.02 | 0.61 | 12.97 | 504.13 |
C# | 0.03 | 0.00 | 0.05 | 0.80 | 28.66 |
MinGW C | 0.02 | 0.00 | 0.09 | 1.05 | 18.69 |
使用コントロール
ToolStrip1(Button1),OpenFileDialog1
プログラムリスト
入出力事例
11:FEM2次元定常飽和浸透流解析(四角形要素)
概要(プログラム:vbSEEP4nod, vcsSEEP4nod, gccSEEP4nod)
- このプログラムは,「三本木茂夫・吉村信敏:コンピュータによる構造工学講座I-1-B 日本鋼構造協会編 有限要素法による構造解析プログラム 考え方と解説,培風館,昭和54年5月」に掲載されている FORTRAN による三角形定ひずみ要素の2次元 FEM 応力解析プログラムを参考に,4 節点パラメトリック要素を用いた定常飽和浸透流解析プログラムとして作成したものです.
- 挙動の検証結果は pdf 解説に記載しています.
- 扱える境界条件は不透水境界・全水頭指定・流量指定の 3 種です.透水係数の全体座標系での直交異方性を考慮できます.
- 連立一次方程式を解く上での水頭指定節点の処理方法は,非定常熱伝導解析の温度指定節点処理あるいは応力解析における変位指定節点処理と同一です.
- 非定常解析・不飽和解析には対応していません.
- 扱える条件を以下に示します.
不透水境界 | 解析領域境界で境界条件として何も入力しなければ不透水境界となる |
全水頭指定 | 解析領域境界・領域内部で全水頭を指定する節点番号を入力しその全水頭値を指定する |
流量指定 | 解析領域境界・領域内部で節点流量を指定する節点番号を入力しその流量値を指定する.領域への流入を正,流出を負とする |
解析断面 | 解析断面が水平断面か鉛直断面かを指定する.これにより解として得られた全水頭から位置水頭を減じた圧力水頭を出力するか判定している |
- StatusStrip には処理開始時刻と処理終了時刻および入出力時間を除いた計算時間を示すようにしました.
- プログラムを開始するとまず OpenFileDialog が開きますので入力データファイルを指定します.次に SaveFileDialog が開きますので,出力ファイル名を指定します.計算が終わると MessageBox が計算完了を通知します.
気分的なもの
定常飽和浸透流問題の基礎微分方程式は,非定常熱伝導解析のものを定常にしたものと基本的に同じですので,このプログラムは非定常熱伝導解析のプログラムを修正して作りました.ただし機能がシンプルな分コードは分かり易いと思います.熱伝導のプログラムがあったのでいつか浸透流も作ろうと思っていたものです.
入出力事例
- 入出力事例は,機能確認を行うための単純な境界条件を有するモデルとしています.
- 解説書には,圧力勾配を有する水平断面領域にドレーンを配置した場合の事例と,設置した鋼矢板の上下流で水頭差を有する場合の事例を載せました.
- 入出力書式は解説書の項のpdfファイルを参照してください.
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト
入出力事例
12:FEM2次元定常飽和浸透流解析(三角形要素)
概要(プログラム:vbSEEP3nod, vcsSEEP3nod, gccSEEP3nod)
- 四角形要素のものを三角形要素に書き換えたものです.
- 非定常解析・不飽和解析には対応していません.
- 扱える条件を以下に示します.
不透水境界 | 解析領域境界で境界条件として何も入力しなければ不透水境界となる |
全水頭指定 | 解析領域境界・領域内部で全水頭を指定する節点番号を入力しその全水頭値を指定する |
流量指定 | 解析領域境界・領域内部で節点流量を指定する節点番号を入力しその流量値を指定する.領域への流入を正,流出を負とする |
解析断面 | 解析断面が水平断面か鉛直断面かを指定する.これにより解として得られた全水頭から位置水頭を減じた圧力水頭を出力するか判定している |
- StatusStrip には処理開始時刻と処理終了時刻および入出力時間を除いた計算時間を示すようにしました.
- プログラムを開始するとまず OpenFileDialog が開きますので入力データファイルを指定します.次に SaveFileDialog が開きますので,出力ファイル名を指定します.計算が終わると MessageBox が計算完了を通知します.
使用コントロール
ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,TableLayoutPanel1,Label1〜11
プログラムリスト