WANtaroHP (F90 FEM)

 

Contents


  • 連立一次方程式は,バンドマトリックスを作成しコレスキー法により解いています.
  • 三角形要素として,線形3節点要素が用いられています.
  • 四角形要素として,4ガウス積分点を有する4節点アイソパラメトリック要素が用いられています.
  • ここに掲載のプログラムでは,三角形要素と四角形要素の混合問題はとくことはできません. 三角形要素か四角形要素のどちらかを,問題を解く前に選ぶ必要があります.
  • 節点番号最適化プログラムが,2次元熱伝導解析・2次元浸透流解析・2次元応力解析・軸対称応力解析プログラムに対して準備されています. これらのプログラムの目的は,剛性行列のバンド幅を縮小することにより計算時間を短縮することにあります. 節点番号最適化プログラムの主要サブルーチンは,こちらのページ(http://www.archi.hiro.kindai.ac.jp/laboratory/SAL/dfujii/Report/fem/fem_5.pdf)を参照して作成しました.


2次元要素生成プログラム

プログラム概要

  • このプログラムは,三角形あるいは四角形の2次元要素を生成するプログラムです.
  • このプログラムは,谷口健男著,「FEMのための要素自動分割 デローニー三角分割法の利用」に掲載の FORTRAN プログラムの一部(全部ではない)を,Fortran 90 に書き直したものです.
  • 解説文書 (TeX_automesh.pdf) 中の図は GMT により作成されています.このため,GMTによる作図補助プログラム (f90_GMT_MESH) も掲載しています.
  • 実行用バッチファイル事例 (bat_a1T および bat_a1Q) の中には,GMT による作図コマンドが含まれています.
  • 三角形要素生成は,プログラム f90_FEM_MESH3 を1回実行することにより可能です.ただし入力データによっては失敗する場合もあります.
  • 四角形要素生成は,f90_FEM_MESH4A,f90_FEM_MESH4B,f90_FEM_MESH4C の3つのプログラムを順次1ずつ実行することにより可能です.ただし入力データによっては失敗する場合もあります.
  • 加えて各要素の領域を定義するプログラム f90_FEM_DOMAIN も示しています. 三角形要素に対しては,プログラムf90_FEM_MESH3の中で領域番号(材料番号)が定義されます.しかし三角形要素において領域の再定義を行いたい場合があるかもしれません. そこで,三角形要素における領域の再定義や,四角形要素における領域定義を行うため,このプログラムを作成しました.
  • 各プログラムの実行用バッチファイルは,「ソースコードおよび関連文書」に含まれています.

ソースコードおよび関連文書

三角形要素生成

FilenameDescription
f90_FEM_MESH3.txt2次元三角形要素生成プログラムソースコード
inp_test1T.txt簡単なテストデータ(1)
inp_test2T.txt簡単なテストデータ(2)
inp_test3T.txt簡単なテストデータ(3)
inp_3nohole.txt少し複雑なテストデータ(1)
inp_3r5mm50h.txt少し複雑なテストデータ(2)
inp_3rand19h.txt少し複雑なテストデータ(3)
bat_a1T.txt実行用バッチファイル事例

四角形要素生成

FilenameDescription
f90_FEM_MESH4A.txt2次元四角形要素生成プログラムソースコード(1)
f90_FEM_MESH4B.txt2次元四角形要素生成プログラムソースコード(2)
f90_FEM_MESH4C.txt2次元四角形要素生成プログラムソースコード(3)
inp_test1Q.txt簡単なテストデータ(1)
inp_test2Q.txt簡単なテストデータ(2)
inp_test3Q.txt簡単なテストデータ(3)
inp_4nohole.txt少し複雑なテストデータ(1)
inp_4r5mm50h.txt少し複雑なテストデータ(2)
inp_4rand19h.txt少し複雑なテストデータ(3)
bat_a1Q.txt実行用バッチファイル事例

解説文書とGMT作図補助プログラム

FilenameDescription
TeX_automesh.pdfプログラム使用法解説書
f90_GMT_MESH.txtGMT作図補助プログラム (メッシュ図のみ対応)
f90_FEM_DOMAIN.txt要素領域定義プログラム
a_domain.txt要素領域定義プログラム実行用バッチファイル事例 (GMT作図コマンド含む)
out_test3T.txtf90_FEM_MESH3 出力事例 (三角形要素)
out_test3Q.txtf90_FEM_MESH4C 出力事例 (四角形要素)
inp_dom.txt要素領域定義ファイル事例
out_test3T.txtf90_FEM_DOMAIN 出力事例 (三角形要素)
out_test3Q.txtf90_FEM_DOMAIN 出力事例 (四角形要素)

a_domain.bat および f90_FEM_DOMAIN 実行結果例を以下に示します.a_domain.bat にはGMT実行コマンドが含まれています.

fig_gmt_mesh_col3T.png
fig_gmt_mesh_col3Q.png

矩形領域分割プログラム

テスト用データ作成などに便利な,矩形領域を四角形要素に分割するプログラムを示す.

ソースコード
program f90_RECTMESH
implicit none
integer::nn,mm
integer::NODT,NELT
integer,allocatable::kakom(:,:)
real(8),allocatable::x(:),y(:)
real(8)::aa,bb,x0,y0
integer::i,j,k,ne,i0
character(len=50)::dummy

call getarg(1,dummy);read(dummy,*) aa
call getarg(2,dummy);read(dummy,*) bb
call getarg(3,dummy);read(dummy,*) nn
call getarg(4,dummy);read(dummy,*) mm
call getarg(5,dummy);read(dummy,*) x0
call getarg(6,dummy);read(dummy,*) y0

NODT = (nn + 1) * (mm + 1)
NELT = nn * mm
allocate(kakom(1:NELT, 1:4))
allocate(x(1:NODT),y(1:NODT))

k = 0
do i = 0,mm
do j = 0,nn
k = k + 1
x(k) = aa / dble(nn) * dble(j) + x0
y(k) = bb / dble(mm) * dble(i) + y0
end do
end do
ne = 0
do k = 0,mm - 1
i0 = 1 + (nn + 1) * k
do j = 0,nn - 1
i = i0 + j
ne = ne + 1
kakom(ne, 1) = i
kakom(ne, 2) = i + 1
kakom(ne, 3) = kakom(ne, 2) + nn + 1
kakom(ne, 4) = kakom(ne, 3) - 1
end do
end do

write(6,'(2(i10))') NODT,NELT
do ne=1,NELT
write(6,'(5(i10))') (kakom(ne,j),j=1,4),ne
end do
do i=1,NODT
write(6,'(2(e15.7),i10)') x(i),y(i),i
end do

end program f90_RECTMESH
実行命令書式
f90_RECTMESH aa bb nn mm x0 y0 > out_mesh.txt
f90_RECTMESH実行プログラム
aa x方向領域長さ
bb y方向領域長さ
nn x方向分割数
mm y方向分割数
x0 領域左下のx座標
y0 領域左下のy座標
out_mesh.txt出力ファイル名


2次元非定常熱伝導解析

プログラム概要

  • このプログラムは,2次元非定常熱伝導解析を行うためのものです.(source code: f90_FEM_HEAT.f90)
  • 4節点アイソパラメトリック要素(4ガウス積分点)が用いられています.1節点の自由度は1です.
  • 要素厚さは一定(単位厚さ)であり変えられません.
  • 以下の条件を扱えます.
    項目説明
    断熱境界 何も指定しなければ断熱境界となる.
    温度指定境界指定した節点に温度履歴を与える.
    熱伝達境界 指定した要素の辺に熱伝達率および外部温度履歴を与える.
    発熱材料 セメントコンクリートのような発熱材料を扱える.
  • 発熱材料の式
              T=K*[1-exp(-α*t)]
       T: 温度上昇量, K: 最高上昇温度, t: 時間, α: 発熱速度を表す係数
  • 時間の離散化には,Crank-Nicolson法が用いられています.
  • 節点番号最適化プログラム (source code: f90_NUMheat.f90) を用いれば,マトリックスのバンド幅と計算時間の短縮が図れます.
  • Windowsコマンドプロンプトでのバッチコマンドは以下のとおりです.
001 gfortran -o f90_NUMheat.exe f90_NUMheat.f90
002 gfortran -o f90_FEM_HEAT.exe f90_FEM_HEAT.f90
003
004 f90_NUMheat inp_org_heat.csv inp_ren_heat.csv
005 f90_FEM_HEAT inp_ren_heat.csv out_ren_heat.csv
001
'f90_NUMheat' (節点番号最適化プログラム)のコンパイル
002
'f90_FEM_HEAT' (熱伝導解析FEMプログラム)のコンパイル
004
'f90_NUMheat' (節点番号最適化プログラム)の実行, 入力ファイル: inp_org_heat.csv, 出力ファイル: inp_ren_heat.csv
005
'f90_FEM_HEAT' (熱伝導解析FEMプログラム)の実行, 入力ファイル: inp_ren_heat.csv, 出力ファイル: out_ren_heat.csv

ソースコードおよび関連文書

FilenameDescription
f90_FEM_HEAT.txt熱伝導解析FEMプログラムソース
f90_NUMheat.txt節点番号最適化プログラムソース
en_TeX_data_io_heat.pdf入出力書式説明用文書
TeX_theory_heat.pdf解析理論説明文書
en_TeX_ex_heat.pdf解析事例


2次元定常浸透流解析 (飽和-不飽和浸透流解析)

プログラム概要

  • このプログラムは,2次元定常浸透流解析を行うためのものです. (source code: f90_FEM_SEEP.f90)
  • 線形三角形要素あるいは4節点アイソパラメトリック要素(4ガウス積分点)を使用できます.1節点の自由度は1です. このプログラムでは三角形要素と四角形要素の混合問題は扱えません.使用に当たり三角形要素か四角形要素のどちらかを選ぶ必要があります.
  • このプログラムは,鉛直断面内の飽和-不飽和浸透流解析および水平断面内の飽和浸透流解析を行うことができます.
  • 不飽和地盤モデルとして Von Genuchten モデル が使われています.
  • 下記の条件が扱えます.
    項目説明
    不透水境界 何も指定しなければ不透水境界となる.
    全水頭指定境界節点番号と全水頭を指定する.
    流量指定境界 節点番号と流量を指定する.正の流量は流入,負の流量は流出を表す.
    浸潤点境界 浸出点となる可能性のある節点を指定する.それらの節点の圧力水頭はゼロあるいは負値であり流量の符号は負(流出)となる.
    解析断面 鉛直断面か水平断面かを選択する.鉛直断面での圧力水頭は全水頭と位置水頭との差で与えられる.水平断面での圧力水頭は全水頭と等しい.
  • 節点番号最適化プログラム (source code: f90_NUMseep.f90) を用いれば,マトリックスのバンド幅と計算時間の短縮が図れます.
  • Windowsコマンドプロンプトでのバッチコマンドを以下に示します.
001 gfortran -o f90_NUMseep.exe f90_NUMseep.f90
002 gfortran -o f90_FEM_SEEP.exe f90_FEM_SEEP.f90
003
004 f90_NUMseep inp_org_seep.csv inp_ren_seep.csv
005 f90_FEM_HEAT inp_ren_seep.csv out_ren_seep.csv
001
'f90_NUMseep' (節点番号最適化プログラム)のコンパイル
002
'f90_FEM_SEEP' (浸透流解析FEMプログラム)のコンパイル
004
'f90_NUMseep' (節点番号最適化プログラム)の実行.入力ファイル: inp_org_seep.csv, 出力ファイル: inp_ren_seep.csv
005
'f90_FEM_SEEP' (浸透流解析FEMプログラム)の実行.入力ファイル: inp_ren_seep.csv, 出力ファイル: out_ren_seep.csv

ソースコードおよび関連文書

FilenameDescription
f90_FEM_SEEP.txt浸透流解析FEMプログラムソースコード
f90_NUMseep.txt節点番号最適化プログラムソースコード
en_TeX_data_io_seep.pdf入出力書式説明文書
TeX_theory_seep.pdf解析理論説明文書
en_TeX_ex_seep1.pdf解析事例(ダム)
en_TeX_ex_seep2.pdf解析事例(トンネル)


2次元応力解析 (no-tension 解析)

プログラム概要

  • このプログラムは,2次元応力解析を行うためのものです. (source code: f90_FEM_PLNT.f90)
  • 線形三角形要素あるいは4節点アイソパラメトリック要素(4ガウス積分点)を使用できます.1節点の自由度は2です. このプログラムでは三角形要素と四角形要素の混合問題は扱えません.使用に当たり三角形要素か四角形要素のどちらかを選ぶ必要があります.
  • このプログラムは,平面応力あるいは平面ひずみ解析を行うことができます.
  • このプログラムは,no-tension 解析を行うこともできます.
  • no-tension 状態を追跡するために,Zienkiewicz により提案された応力遷移法が用いられています. この方法では,釣り合い状態は,線形剛性行列と内力の座標変換を用いた繰り返し計算により得られます. no-tension状態でのポアソン比はゼロとしています.
  • 以下の条件を扱うことができます.
    項目説明
    荷重 載荷される節点および荷重値を指定する.
    温度 全節点の温度変化量を指定する.
    慣性力 全要素に対し,重力加速度に対する比として水平・鉛直方向の慣性力を指定する.
    変位 変位指定する節点およびその変位量を指定する.ゼロを含む任意の値を入力できる.
    No-tension 材料全要素に対し,引張強さを指定する. ある要素において最大主応力が引張強さを超えた場合,その要素はno-tension材料として扱う.
  • 節点番号最適化プログラム (source code: f90_NUMplnt.f90) を用いれば,マトリックスのバンド幅と計算時間を短縮することができます.
  • プログラム 'f90_NUMplnt.exe' は内部で 'f90_num.exe' を呼んでいるので,'f90_NUMplnt.exe'と同じディレクトリに実行ファイル 'f90_num.exe' が存在することが必要です.
  • Windowsコマンドプロンプトでのバッチコマンドは以下のとおりです.
001 gfortran -o f90_num.exe f90_num.f90
002 gfortran -o f90_NUMplnt.exe f90_NUMplnt.f90
003 gfortran -o f90_FEM_PLNT.exe f90_FEM_PLNT.f90
004
005 f90_NUMplnt inp_org_plnt.csv inp_ren_plnt.csv
006 f90_FEM_PLNT inp_ren_plnt.csv out_ren_plnt.csv
001
節点番号最適化プログラム主要部 'f90_num' のコンパイル
('f90_num.exe' が 'f90_NUMplnt.exe' と同じディレクトリに存在すれば,'f90_num.exe' は自動的に 'f90_NUMplnt.exe' から呼び出されます)
002
節点番号最適化プログラム制御部 'f90_NUMplnt' のコンパイル
003
2次元応力解析FEMプログラム 'f90_FEM_PLNT' のコンパイル
004
節点番号最適化プログラム 'f90_NUMplnt' の実行.入力ファイル: inp_org_plnt.csv, 出力ファイル: inp_ren_plnt.csv
006
2次元応力解析FEMプログラム 'f90_FEM_PLNT' の実行.入力ファイル: inp_ren_plnt.csv, 出力ファイル: out_ren_plnt.csv

ソースコードおよび関連文書

FilenameDescription
f90_FEM_PLNT.txt2次元応力解析FEMプログラムソースコード
f90_NUMplnt.txt節点番号最適化プログラム制御部ソースコード
f90_num.txt節点番号最適化プログラム主要部ソースコード
en_TeX_data_io_plnt.pdf入出力書式説明文書
TeX_theory_plnt.pdf解析理論説明文書


軸対称応力解析 (no-tension 解析)

プログラム概要

  • このプログラムは,軸対称応力解析をおこなうためのものです.(source code: f90_FEM_ASNT.f90)
  • 4節点アイソパラメトリック要素(4ガウス積分点)が使われています.1節点の自由度は2です.
  • このプログラムではno-tension解析が行えます.
  • no-tension状態を解析する方法は Zienkiewicz により提案された応力遷移法です.
  • 下記の条件が扱えます.慣性力を扱うことはできません.
    項目説明
    荷重 載荷する節点と荷重値を指定する.
    温度 全節点の温度変化量を指定する.
    慣性力 全要素に対し,重力加速度に対する比として回転軸方向の慣性力を指定する.半径方向慣性力は考慮できない.
    変位 変位指定する節点と変位量を指定する.変位はゼロを含む任意の値を指定できる.
    No-tension 材料全要素の引張強さを指定する.ある要素において最大主応力が引張強さを越えた時点で,その要素はno-tension材料になったとみなす.
  • 節点番号最適化プログラム (source code: f90_NUMasnt.f90) を用いれば,マトリックスのバンド幅と計算時間を短縮することができます.
  • プログラム 'f90_NUMasnt.exe' は内部で 'f90_num.exe' を呼んでいるので,'f90_NUMasnt.exe'と同じディレクトリに実行ファイル 'f90_num.exe' が存在することが必要です.
  • Windowsコマンドプロンプトでのバッチコマンドを以下に示します.
001 gfortran -o f90_num.exe f90_num.f90
002 gfortran -o f90_NUMasnt.exe f90_NUMasnt.f90
003 gfortran -o f90_FEM_ASNT.exe f90_FEM_ASNT.f90
004
005 f90_NUMasnt inp_org_asnt.csv inp_ren_asnt.csv
006 f90_FEM_ASNT inp_ren_asnt.csv out_ren_asnt.csv
001
節点番号最適化プログラム主要部 'f90_num' のコンパイル
('f90_num.exe' が 'f90_NUMasnt.exe' と同じディレクトリに存在すれば,'f90_num.exe' は自動的に 'f90_NUMasnt.exe' から呼び出されます)
(このプログラムは2次元応力解析で用いているものと同一です)
002
節点番号最適化プログラム制御部 'f90_NUMasnt' のコンパイル
003
2次元応力解析FEMプログラム 'f90_FEM_ASNT' のコンパイル
004
節点番号最適化プログラム 'f90_NUMasnt' の実行.入力ファイル: inp_org_asnt.csv, 出力ファイル: inp_ren_asnt.csv
006
軸対称応力解析FEMプログラム 'f90_FEM_ASNT' の実行.入力ファイル: inp_ren_asnt.csv, 出力ファイル: out_ren_asnt.csv

ソースコードおよび関連文書

FilenameDescription
f90_FEM_ASNT.txt軸対称応力解析FEMプログラムソースコード
f90_NUMasnt.txt節点番号最適化プログラム制御部ソースコード
f90_num.txt節点番号最適化プログラム主要部ソースコード
en_TeX_data_io_asnt.pdf入出力書式説明文書
TeX_theory_asnt.pdf解析理論説明文書
en_TeX_ex_asnt.pdf解析事例


平面骨組構造解析

プログラム概要

  • このプログラムは,平面骨組構造解析をおこなうためのものです.(Source code: f90_FEM_FRAME.f90)
  • このプログラムは,梁要素を用いた骨組み構造の線形弾性解析を行うことができます.梁要素は1要素2節点,1節点当たり3自由度を有します.
  • 下記の条件が考慮可能です.
    項目説明
    荷重 載荷する節点と荷重値を指定する
    温度 全節点の温度変化量を指定する.
    慣性力重力加速度に対する比率として,全要素の水平方向・鉛直方向の慣性力を指定する.
    変位 変位指定する節点および変位値を指定する.変位値としては,ゼロを含む任意の値を指定できる.
  • Windowsコマンドプロンプトでのバッチコマンドを以下に示します.
001 gfortran -o f90_FEM_FRAME.exe f90_FEM_FRAME.f90
002
003 f90_FEM_FRAME inp_org_frame.csv out_org_frame.csv
001
平面骨組構造解析FEMプログラム 'f90_FEM_FRAME' のコンパイル
003
平面骨組構造解析FEMプログラム 'f90_FEM_FRAME' の実行.入力ファイル: inp_org_frame.csv, 出力ファイル: out_org_frame.csv

ソースコードおよび関連文書

FilenameDescription
f90_FEM_FRAME.txt平面骨組構造解析FEMプログラムソースコード
en_TeX_data_io_frame.pdf入出力書式説明文書
TeX_theory_frame.pdf解析理論説明文書


平面骨組座屈解析

プログラム概要

  • このプログラムは,Fortran90による,平面骨組の座屈解析を行うプログラムです.
  • このプログラムにより,座屈荷重,座屈変形モードが得られます.
  • このプログラムでは,固有値および固有ベクトルを求めるために,一般化ヤコビ法が用いられています.
  • このプログラムでは,変形の前後で荷重の方向は変化しないことが仮定されていることに注意してください. もし,変形の前後で荷重の方向が変化する場合の挙動を知りたい場合は,違うプログラムを用いる必要があります.

ソースコードおよび関連文書

FilenameDescription
f90_EGV_BFRM.txt平面骨組座屈解析プログラムF90ソース
TeX_EGV_B.pdf解析事例


平面骨組幾何学的非線形解析

プログラム概要

  • このプログラムは,Fortran90による,平面骨組の幾何学的非線形解析を行うプログラムです.
  • 弾性骨組み部材の大きな変位を追跡することができます.
  • 増分形剛性方程式は,弧長増分法により解かれています.
  • 骨組み部材の内力は,剛体変位を除去した変位と線形剛性行列を用いて算定されています.
  • このプログラムでは,変形の前後で荷重の方向は変化しないことが仮定されていることに注意してください.

ソースコードおよび関連文書

FilenameDescription
f90_FEM_GFRMS.txt平面骨組幾何学的非線形解析プログラムF90ソース
TeX_FDA_frame.pdf解析事例


平面骨組幾何学的非線形解析

プログラム概要

  • このプログラムは,Fortran90による,平面骨組の幾何学的非線形解析を行うプログラムです.
  • 弾性骨組み部材の大きな変位を追跡することができます.
  • 増分形剛性方程式は,弧長増分法により解かれています.
  • 骨組み部材の内力は,剛体変位を除去した変位と線形剛性行列を用いて算定されています.
  • このプログラムでは,変形の前後で荷重の方向は変化しないことが仮定されていることに注意してください.

ソースコードおよび関連文書

FilenameDescription
f90_FEM_GFRMS.txt平面骨組幾何学的非線形解析プログラムF90ソース
TeX_FDA_frame.pdf解析事例


平面トラス構造解析

プログラム概要

  • このプログラムは,平面トラスの構造解析を行うプログラムです.
  • ここ何十年?もトラスの解析をするニーズはなかったのですが,近い将来必要になりそうだったので,平面骨組用プログラムを改造してにわかに作成したものです.
  • 断面力は全体剛性行列 TSM に格納して表示させるようにしていましたが,原因不明のバグ?により途中で値が変わり,無視できない大きさのせん断力が入ってしまう場合があるため, 断面力を格納する配列 sresul を新たに作りました.Frame解析では大丈夫だったのになぜだろう?
  • 下記の条件が考慮可能です.
    項目説明
    荷重 載荷する節点と荷重値を指定する
    温度 全節点の温度変化量を指定する.
    慣性力重力加速度に対する比率として,全要素の水平方向・鉛直方向の慣性力を指定する.
    変位 変位指定する節点および変位値を指定する.変位値としては,ゼロを含む任意の値を指定できる.

ソースコードおよび関連文書

FilenameDescription
f90_FEM_TRUSS.txt平面トラス構造解析プログラムF90ソース


pic
inserted by FC2 system