csv csv

2009年1月9日

2

次元自動要素分割プログラム

—

三角形要素・四角形要素

—

目次

1. はじめに 1 2. 三角形要素の自動分割 2 2.1 プログラムの特徴 . . . 2 2.2 データ入力書式（csv形式）. . . 3 2.3 出力ファイル書式（csv形式） . . . 3 3. 四角形要素の自動分割 4 3.1 プログラムの特徴 . . . 4 3.2 データ入力書式（csv形式）. . . 5 3.3 出力ファイル書式（csv形式） . . . 5 4. 単純な要素分割事例 6 5. 多くの内部境界を含む要素分割事例 8 6. 自動分割要素を用いた解析事例 11

1.

はじめに

本プログラムは，谷口健男著，「FEMのための要素自動分割 デローニー三角分割法の利用」に掲載の

FORTRANプログラムを，Microsoft Visual Basic 2008 Express Edition に書き直したものである．

機能は，円孔などの内部境界も含め，任意の形状の2次元領域を三角形要素あるいは四角形要素に自動分割 するものである．三角形要素と四角形要素の混在には対応していない． 三角形要素では，節点数，要素数，要素-節点関係，要素-領域関係，節点座標が出力され，物性値の異なる 領域にも自動対応できる． 四角形要素では，入力時は領域毎にデータ入力を行うが，出力時は要素-領域関係が出力されないため，手 動で対応する必要がある．なお，原典では四角形要素作成は3種類のプログラムにより行うこととなっている が，ここではVBのクラスの考え方を活用し， 1「三角形要素の生成」， 2「三角形要素から四角形要素への 変換」， 3「生成された四角形要素の形状修正」という過程を，1プログラムで順次実行し，結果をテキスト ファイル出力すると共に要素分割図をpng画像として出力している． 三角形要素・四角形要素ともに要素分割を行うのみで境界条件や荷重設定を行う機能は含んでいないため， これらは別途設定する必要がある．

参考文献

[1] 谷口健男著：FEMのための要素自動分割 デローニー三角分割法の利用，森北出版（株），2006年8月

2.

三角形要素の自動分割

2.1 プログラムの特徴 指定した総節点数となるよう三角形要素群が生成される． 外部境界，内部境界，境界内部節点として入力した節点は，出力時も同一番号・同一座標が保証される． 入力時に指定した以外の節点については法則性はなく自動的に生成・番号付けされる． 要素番号については外部境界数で指定した順に要素-領域関係の番号が付されるが，要素番号は法則性 はなく自動的に生成・番号付けされる． 出力は節点数，要素数，要素-節点関係，要素-領域関係，節点座標であり，物性値の異なる領域にも自 動対応できる． プログラムでは，Form1で「入出力ファイル指定と作図」，Form2で「要素分割計算」を行っている． 処理の流れは Form11 で入出力ファイルを指定， Form22 で要素分割実行と出力ファイルへの書き出 し， Form13 でForm2で書き出した出力ファイルを読み込み作図を実行，という流れになっている．

入出力ファイル名のForm間でのやりとりは，Form1で指定したファイル名をForm1のLabelコン

トロールに書き込み，これをForm2で「fname1 = My.Forms.Form1.Label1.Text」としてForm1の LabelコントロールのTextを読み込むことで行っている． 入力データ不良等により計算が途中で止まった場合は「Error終了メッセージ」を表示しプログラムを 終了する． 入力データ不良がなく「Error終了メッセージ」が表示された場合，配列を大きくとることにより解決 できる場合がある． 最終的に要素面積が0以下のものが存在する場合に「要素分割不良」のメッセージを出力する． 要素分割が細かい場合，節点番号・要素番号が見にくくなるため，Form1の画面で「節点図・要素図 拡大比率」を入力するようにしている．デフォルト値は１であるが，節点番号・数値番号が見にくい場 合，この数値を２あるいは３にすることにより数値が小さく表示されるようになるので，画像を拡大し て数値確認を行えるようにした．

2.2 データ入力書式（csv形式） 入力項目 備考 NODT 総節点数（指定値） NEX,MIN 外部境界数，内部境界数 ibex(i) 個々の外部境界上の節点数 ・・・外部境界数(NEX行)繰り返し入力・・・ ibin(i) 個々の内部境界上の節点数 ・・・内部境界数(MIN行)繰り返し入力・・・ 内部境界数が0なら省略 ibno(i,1),・・・ibno(i,ibex(i)) 個々の外部境界上の節点番号を1行に入力 ・・・外部境界数分(NEX行)繰り返し入力・・・ NOB,NIB 境界上の節点総数，境界内部の節点総数 px(i),py(i) 節点のx座標およびy座標 ・・・節点数分(NOB+NIB行)繰り返し入力・・・ ■注意事項 外部境界は時計回りに，内部境界は反時計回りに節点番号を入力する． 外部境界上の節点番号ibno(i,j)は循環させず1個手前の節点までの入力とする． 節点番号は， 1外部境界， 2内部境界， 3境界内部節点の順につけ，節点座標は節点番号順に入力する． 2.3 出力ファイル書式（csv形式） 出力内容 NODE,NELM mtj(i,1),mtj(i,2),mtj(i,3),idm(i) mtj(i,j)：要素-節点関係 idm(i)：要素-領域関係 ・・・要素数分(NELM行)繰り返し出力・・・ px(i),py(i) 節点の(x,y)座標 ・・・節点数分(NODE行)繰り返し出力・・・

3.

四角形要素の自動分割

3.1 プログラムの特徴 総節点数は指定はできず，自動設定される． 外部境界・内部境界上では指定した節点の中間に1節点追加されて四角形要素群が生成される． 要素の大きさは境界上の節点間隔と領域毎入力する格子間隔に支配される． 外部境界，内部境界，境界内部節点として入力した節点は，出力時も同一番号・同一座標が保証される． 入力時に指定した以外の節点については法則性はなく自動的に生成・番号付けされる． 出力は節点数，要素数，要素-節点関係，節点座標である．入力時は領域毎にデータ入力を行うが，出 力時は要素-領域関係が出力されないため，材料種別の指定は手動で対応する必要がある．

プログラムでは，Form1で「入出力ファイルの指定と作図」，Form2で「三角形要素の生成」，Form3

で「三角形要素から四角形要素への変換」，Form4で「生成された四角形要素の形状修正」という過程

を順次実行しており，それぞれのデータの受け渡しは，ファイル出入力で行っている．このため入力

ファイル格納フォルダに，Form2での計算結果を格納するworkf2.csvとForm3での計算結果を格納

するworkf3.csvという作業用ファイルが作成される．

処理の流れは Form11 で入出力ファイルを指定， Form2∼42 で要素分割実行と出力ファイルへの書

き出し， Form13 でForm4で書き出した出力ファイルを読み込み作図を実行，という流れになって

いる．

入出力ファイル名のForm間でのやりとりは，Form1で指定したファイル名をForm1のLabelコント

ロールに書き込み，これをForm2以降で「fname1 = My.Forms.Form1.Label1.Text」としてForm1 のLabelコントロールのTextを読み込むことで行っている． 入力データ不良等により計算が途中で止まった場合は「Error終了メッセージ」を表示しプログラムを 終了する． 入力データ不良がなく「Error終了メッセージ」が表示された場合，配列を大きくとることにより解決 できる場合がある． 最終的に要素面積が0以下のものが存在する場合に「要素分割不良」のメッセージを出力する． 要素分割が細かい場合，節点番号・要素番号が見にくくなるため，Form1の画面で「節点図・要素図 拡大比率」を入力するようにしている．デフォルト値は１であるが，節点番号・数値番号が見にくい場 合，この数値を２あるいは３にすることにより数値が小さく表示されるようになるので，画像を拡大し て数値確認を行えるようにした．

3.2 データ入力書式（csv形式） 入力項目 備考 NEX,MIN 外部境界数，内部境界数 ibex(i),delx(i) ibex(i)：個々の外部境界上の節点数 ・・・外部境界数(NEX行)繰り返し入力・・・ delx(i)：各部分領域の格子間隔 ibin(i) 個々の内部境界上の節点数 ・・・内部境界数(MIN行)繰り返し入力・・・ 内部境界数が0なら省略 ibno(i,1),・・・ibno(i,ibex(i)) 個々の外部境界上の節点番号を1行に入力 ・・・外部境界数分(NEX行)繰り返し入力・・・ NOB,NIB 境界上の節点総数，境界内部の節点総数 px(i),py(i) 節点のx座標およびy座標 ・・・節点数分(NOB+NIB行)繰り返し入力・・・ ■注意事項 外部境界は時計回りに，内部境界は反時計回りに節点番号を入力する． 各部分領域の格子間隔delx(i)は，外部境界で指定した座標間隔と同程度が良いようである． 外部境界上の節点番号ibno(i,j)は循環させず1個手前の節点までの入力とする． データ入力における三角形要素生成プログラムとの違いは，総節点番号(NODT)の入力がないこと， 各部分領域の格子間隔(delx)の入力を必要とすることの2点のみである． 3.3 出力ファイル書式（csv形式） 出力内容 NODE,NELM mtj(i,1),mtj(i,2),mtj(i,3),mtj(i,4) mtj(i,j)：要素-節点関係 ・・・要素数分(NELM行)繰り返し出力・・・ px(i),py(i) 節点の(x,y)座標 ・・・節点数分(NODE行)繰り返し出力・・・

4.

単純な要素分割事例

下図に示す3個の領域（I，II，III）と領域IIの内部に1個の内部境界を有するモデルを要素分割する．境 界上の総節点数は10とし，内部節点1個を領域I内に配置する． • • • • • • • • • • • 1 2 3 4 5 6 7 8 9 10 11 領域I 領域II 領域 III 表1 三角形要素分割入力データ 入力データ 変数記号 備考 20 NODT 総節点数 3,1 NEX,MIN 外部境界数，内部境界数 4 ibex(1) 1番目(領域I)の外部境界節点数 4 ibex(2) 2番目(領域II)の外部境界節点数 3 ibex(3) 3番目(領域III)の外部境界節点数 3 ibin(1) 1番目の内部境界節点数 1,2,3,4 ibno(1,1),ibno(1,2),ibno(1,3),ibno(1,4) 1番目(領域I)の外部境界上の節点番号 4,3,5,6 ibno(2,1),ibno(2,2),ibno(2,3),ibno(2,4) 2番目(領域II)の外部境界上の節点番号 6,5,7 ibno(3,1),ibno(3,2),ibno(3,3) 3番目(領域III)の外部境界上の節点番号 10,1 NOB,NIB 境界上の節点総数，境界内部の節点総数 0,0 px(1),py(1) 外部境界上の節点座標（節点1） 0,10 px(2),py(2) 外部境界上の節点座標（節点2） 5,10 px(3),py(3) 外部境界上の節点座標（節点3） 5,0 px(4),py(4) 外部境界上の節点座標（節点4） 25,10 px(5),py(5) 外部境界上の節点座標（節点5） 25,0 px(6),py(6) 外部境界上の節点座標（節点6） 30,5 px(7),py(7) 外部境界上の節点座標（節点7） 10,3 px(8),py(8) 内部境界上の節点座標（節点8） 15,3 px(9),py(9) 内部境界上の節点座標（節点9） 12,7 px(10),py(10) 内部境界上の節点座標（節点10） 2.5,5 px(11),py(11) 境界内部節点座標（節点11） 表2 四角形要素分割入力データ 入力データ 変数記号 備考 3,1 NEX,MIN 外部境界数，内部境界数 4,5 ibex(1),delx(1) 1番目(領域I)の外部境界節点数と格子間隔 4,5 ibex(2),delx(2) 2番目(領域II)の外部境界節点数と格子間隔 3,5 ibex(3),delx(3) 3番目(領域III)の外部境界節点数と格子間隔 3 ibin(1) 1番目の内部境界節点数 1,2,3,4 ibno(1,1),ibno(1,2),ibno(1,3),ibno(1,4) 1番目(領域I)の外部境界上の節点番号 4,3,5,6 ibno(2,1),ibno(2,2),ibno(2,3),ibno(2,4) 2番目(領域II)の外部境界上の節点番号 6,5,7 ibno(3,1),ibno(3,2),ibno(3,3) 3番目(領域III)の外部境界上の節点番号 10,1 NOB,NIB 境界上の節点総数，境界内部の節点総数 0,0 px(1),py(1) 外部境界上の節点座標（節点1） 0,10 px(2),py(2) 外部境界上の節点座標（節点2） 5,10 px(3),py(3) 外部境界上の節点座標（節点3） 5,0 px(4),py(4) 外部境界上の節点座標（節点4） 25,10 px(5),py(5) 外部境界上の節点座標（節点5） 25,0 px(6),py(6) 外部境界上の節点座標（節点6） 30,5 px(7),py(7) 外部境界上の節点座標（節点7） 10,3 px(8),py(8) 内部境界上の節点座標（節点8） 15,3 px(9),py(9) 内部境界上の節点座標（節点9） 12,7 px(10),py(10) 内部境界上の節点座標（節点10）

要素分割結果は下図のとおりである． 三角形要素では，節点数50を指定しているため，節点数50，要素数74の分割が得られている． 四角形要素では，節点数52，要素数42の分割が得られている． いずれも，入力として指定した節点については節点番号・座標とも入力と同一のものが得られており，その 他節点は自動生成されたものとなる． 三角形要素 節点番号 三角形要素 要素番号 四角形要素 節点番号 四角形要素 要素番号 図1 三角形および四角形要素分割図

5.

多くの内部境界を含む要素分割事例

寸法が幅150mm，高さ300mmの矩形板に孔を開けた試験体の要素分割を行った．供試体の上部30mmは 載荷盤をモデル化している． 基本的な境界の与え方は以下の通り． 外部境界は供試体の周に沿って15mmピッチで指定した． 内部境界（孔の部分）の設定条件は以下の通り． • 供試体の周辺境界から孔中心位置の離隔は半径の2倍以上とした． • 内部境界どおしの中心位置の離隔は半径の2.5倍以上とした． • 楕円の場合の離隔判定のための半径はそれぞれの長軸半径とした． • 内部境界は円形もしくは楕円形とし，周の分割数は12とした．なお円孔では中心角を12等分の 分割とし，楕円孔では短軸長を6等分することにより境界座標を設定している． • 内部境界の周りにメッシングの補助として境界内部節点を配置することも考えられるが，ここでは 境界内部節点は設置しない（NIB=0） 入力境界条件に対し，三角形要素では与えた節点数となるよう分割が進められ，四角形要素では与えた境界 節点の中間に1節点新たな節点を追加して分割が進められる． 下図に外部境界・内部境界の与え方の事例を示す． 外部境界は節点を時計回りに入力することから，供試体部では1∼21∼31∼51∼60の順に節点番号を，載 荷盤部では61∼71∼31∼21の順に節点番号を入力する． 内部境界は節点を反時計回りに入力することから，72∼76，77∼81の順に節点番号を入力する． 入力は，外部境界，内部境界，境界内部節点（設定する場合）の順に行う．節点座標は節点番号順に入力 する． - 150mm(15mm@10) 6 ? 300mm(15mm@20) 6 ?30mm 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 51 50 49 48 47 46 45 44 43 42 41 40 39 38 37 36 35 34 33 32 31 22 23 24 25 26 27 28 29 30 60 59 58 57 56 55 54 53 52 61 62 63 64 65 66 67 68 69 70 71 • • • • • • • • • • • • 72 73 74 75 76 77 78 79 80 81 図2 外部境界・内部境界の与え方

■内部境界座標設定のプログラム・テクニック

• 乱数による楕円半径・回転角指定

Dim rando as New Random ’ 乱数は Random クラスの Next メソッドを使用 For i = 1 To MIN ’MIN は内部境界数

aa(i) = CDbl(rando.Next(2, 10)) ’ 横半径を 2∼10mm の乱数で指定 bb(i) = CDbl(rando.Next(2, 10)) ’ 縦半径を 2∼10mm の乱数で指定

alpha(i) = CDbl(rando.Next(180)) / 180.0 * pi ’ 傾き（時計回り）を 0∼180 度の乱数で指定 Next i

• 内部境界の節点指定

Private Sub HOLEXY  _

(ByRef nd As Integer, ByVal n As Integer, ByVal xx As Double, ByVal yy As Double, _ ByVal a As Double, ByVal b As Double, ByVal alp As Double, _

ByRef px() As Double, ByRef py() As Double, ByVal idiv As Integer) Dim i As Integer : Dim m As Integer : Dim dy As Double

Dim xw As Double : Dim yw As Double : Dim r As Double Dim x1 As Double : Dim y1 As Double : Dim theta() As Double ReDim theta(n)

Select Case idiv ’ 中心角の指定 Case 0 ’ 中心角等間隔指定 For i = 1 To n theta(i) = 2.0 * pi / CDbl(n) * CDbl(i - 1) Next i Case 1 ’y 方向等分割指定，分割数 n は偶数に限る m = CInt(n / 2) dy = 2.0 * b / CDbl(m) For i = 1 To m yw = b - dy * CDbl(i - 1) xw = a / b * (-Math.Sqrt(b ^ 2 - yw ^ 2)) r = Math.Sqrt(xw ^ 2 + yw ^ 2) theta(i) = Math.Acos(Math.Abs(yw) / r) If yw < 0.0 Then theta(i) = pi - theta(i) Next i For i = 1 To m yw = -b + dy * CDbl(i - 1) xw = a / b * (Math.Sqrt(b ^ 2 - yw ^ 2)) r = Math.Sqrt(xw ^ 2 + yw ^ 2) theta(m + i) = Math.Acos(Math.Abs(yw) / r) If yw < 0.0 Then theta(m + i) = pi + theta(m + i) Else theta(m + i) = 2.0 * pi - theta(m + i) End If Next i End Select For i = 1 To n nd = nd + 1 ’nd は節点番号

r = a * b / Math.Sqrt(b ^ 2 * Math.Sin(theta(i)) ^ 2 + a ^ 2 * Math.Cos(theta(i)) ^ 2) xw = -r * Math.Sin(theta(i)) ’ 内部境界上節点の x 座標

図3に要素分割結果を示す． ケースの設定は以下の通りであり，上段が三角形要素，下段が四角形要素である． Case 1 内部境界を有しない均一板モデル Case 2 半径5mmの円孔を領域内に50個ランダムに設置したもの Case 3 長軸半径・短軸半径・設置角・設置位置を乱数によりランダムに19個設定したもの 長軸半径・短軸半径は5∼15mmの範囲とした． 三角形要素では領域の指定が出力されるため，供試体部と載荷盤部を色分けしたが，四角形要素では領域指 定は出力されないため供試体部と載荷盤部の色分けは行っていない．しかし四角形要素は指定した外部境界 点の中間に1節点を追加してメッシングされるため，全体領域の上部2段が載荷盤であることが認識できる．

Case 1 Case 2 Case 3

均一板 半径5mm，孔数50 半径5∼15mm，孔数19 三 角 形 要 素 節点数=1000，要素数=1873 節点数=2000，要素数=3371 節点数=2000，要素数=3623 四 角 形 要 素 節点数=903，要素数=840 節点数=3215，要素数=2602 節点数=1708，要素数=1436

6.

自動分割要素を用いた解析事例

前節で作成した要素を用い，2次元弾性FEMによる平面応力状態での引張試験を行った．供試体はコンク リートの物性値を想定し，載荷盤は鉄鋼の物性値を想定した． 表3に解析条件を，表4に解析結果の概要を示す． また図4に三角形要素での解析結果を，図5に四角形要素での解析結果を示す． 主応力分布図については，三角形要素は１要素１応力状態として，四角形要素では１要素４積分点の応力が 算定されるため要素を４分割しそれぞれの積分点での応力で色分けしている． 表3 解析条件 供試体 寸法 幅150mm×高300mm×厚10mm 弾性係数 20600 MPa ポアソン比 0.2 載荷盤 寸法 幅150mm×高30mm×厚10mm 弾性係数 206000 MPa ポアソン比 0.3 応力状態 – 平面応力状態 荷重 位置・大きさ 載荷盤上部より一様引張応力5MPaを載荷 境界条件 x（横）方向 供試体下端中央部・載荷盤上端中央部変位拘束 y（縦）方向 供試体下端全面変位拘束 表4 解析結果

要素種別 項目 Case 1 Case 2 Case 3

均一モデル 円孔50個 ランダム楕円19個 （半径5mm） （半径5∼15mm） 三角形要素 節点数 1000 2000 2000 要素数 1873 3371 3623 計算時間（sec）*1) 29.06 254.36 254.22 供試体上端平均縦変位（mm） 0.07266 0.08768 0.1015 供試体中央横方向変位（mm） 0.00734 0.00937 0.00994 供試体最大引張主応力（MPa） 5.364 17.045 22.969 供試体弾性係数（MPa） 2.064e4 1.711e4 1.479e4

供試体ポアソン比（MPa） 0.202 0.214 0.196 四角形要素 節点数 903 3215 1708 要素数 840 2602 1436 計算時間（sec）*1) 23.35 1227.53 169.26 供試体上端平均縦変位（mm） 0.0727 0.08955 0.10253 供試体中央横方向変位（mm） 0.00734 0.01007 0.01010 供試体最大引張主応力（MPa） 5.431 18.481 24.448

供試体弾性係数（MPa） 2.064e4 1.675e4 1.463e4

Case 1 Case 2 Case 3 均一板 半径5mm，孔数50 半径5∼15mm，孔数19 要 素 分 割 第 一 主 応 力 分 布 第 二 主 応 力 分 布 図4 三角形要素分割および主応力分布（主応力は引張が正）

Case 1 Case 2 Case 3 均一板 半径5mm，孔数50 半径5∼15mm，孔数19 要 素 分 割 第 一 主 応 力 分 布 第 二 主 応 力 分 布