[境界要素法プログラムを設計する(アウトライン)] [ネコ騙し数学]
[境界要素法プログラムを設計する(アウトライン)]
1.開発範囲の決定
これは業務アプリの開発で言えば、顧客がどこまでの成果品を望むか?です。顧客の要求は常にドンダケ~なんですが、開発側としては出来るだけ簡単なものにしたい訳です(仕事したくないから(^^;))。今回の顧客は我々です。つまり開発者です。なので今回は、開発側の事情最優先です(^^)。
境界要素法プログラムの目的は、離散化された境界方程式、
を作る事にほぼ尽きます。しかし(1)はそのままでは解けないのでした。以下の条件を考慮する必要がありました。
[1]境界条件
これは、未知ベクトルψ,qの一部に値を与える補助方程式として実装するのが簡単でした。
[2]角点の処理
これは角点を表す節点jのqjが、qj(1),qj(2)と2つになり、qの成分数が増える事を意味します。それに伴い補助条件の追加も必要でした。しかも追加の仕方は、境界条件の与えられ方により変化しました。
最初から[2]を扱うと面倒そうです。例えば(1)+[1]だけなら、組み立てる連立方程式の形は、
になります。(D1 D2)ブロックは、境界条件で値を与えられる節点番号の列成分が1で、他は全て0の行列です。(d1 d2)tブロックは与える値のベクトルです。
ψとqの成分は、節点jのψjとqjを表しますが、z=(ψ q)tと一つの未知ベクトルとしてまとめると、qjのjとz上での成分番号が、既にこの時点でずれます。幸いψとqの成分数は同じなので、それをnとして、
節点j→ψj→z(j),節点j→qj→z(n+j) (3)
ではあります。z(j)などは、配列zのj成分の意味で使ってます。
これに[2]が加わると、qの成分数が角点が現れるたびに増え、(3)程度の対応ではすまなくなります。当然HとD2の列数も増え、どう増えるかなどを制御しなければなりません。さらにここに新たに補助方程式が加わりますが、加え方は一定ではないのでした。
やっぱり最初は(1)+[1]ですよね?。(1)+[1]のプログラムをテストして動作確認し、自信を持って複雑そうな[2]に進むのが妥当です。(1)+[1]の基本が出来れば、[2]への対処のポイントも絞れるはずです。要するに適切な実践に勝る教師なし!です。
もう一つ考慮すべきは、テストの容易さです。作ったプログラムは必ずテストしなければなりません。例えば全部作ってから動かなかったら、バグの特定にも難儀します。その時に「(1)+[1]まではOKよ」という事実があると、格段に楽になります。
業務アプリの開発でも、開発サイクルを出来るだけ短く区切って頻繁にテストし、より複雑な段階へ進む事が推奨されアジャイル開発と言います。その際、個々の開発サイクルでの成果は一つの製品であり、次の段階はその拡張仕様という見方をします。
2.大枠決め
「実践に勝る教師なし!」なんて言うものだから、すぐにプログラムをボトムアップで作成しようします。気持ちはわかります。ボトムアップ方式なら具体的な数式なんかを、どんどん書いて行けるので。でも「適切な」実践ですからね。開発の基本は常にトップダウンです。
プログラムの目的から始まって、何をやりたいか?何をやるべきか?のアウトラインをトップダウンで決めます。それを詳細化して行くと、局所的なボトムアップ作業へと自然につながります。こうすると常に全体を見渡しながら局所的な具体化作業に集中できます。こうすると「後で考えれば明らかだったよなぁ~」みたいな不用意で不必要なバグは激減します。
もう一つ大事なのは、トップダウン中の方針決定には、必ず明確な理由がある事。プログラム作成は、何をやろうと一長一短。すると局所的な便利さ不便さから、大域的な方針まで変更したくなります。その時、変更するのか我慢すべきか判断が必要です。全体系の整合性を保つ判断根拠になるのが、明確な決定理由です。
自分は一時期アプリ開発で飯を食ってたので、プロの言葉は信じて良いですよぉ~(←ほんとか?(^^;))。
まず実用的なプログラムの全体構成は一個しかありません。図-3です。これ以外ないですよね?。ここで問題になるのは[メモ欄]だと思いますが、これは必須です。ここにプログラムの仕様を書いて行きます。[メモ欄]を読めば、全体を見渡せるように書きます。これを怠ると、たいてい後でこっぴどい目に合います。
トップダウン思考でキーになるのが、「実行者をさがせ!」です。特定された実行者が、プログラムを動かすと考えて下さい。このプログラムの目的である式(2)を見てみます。
(D1 D2)と(d1 d2)tの部分は、どうやら値を並べるだけなのでどうとでもなりそうです。対して(B H)と(w)tの部分は、今までやってきた積分式を使うべきところです。
BやHは境界積分の項です。境界積分は個々の境界要素の境界節点を特定し、主に節点情報に基づいて積分値を得た後、それを全ての要素について足すのでした。その足し算を、行列という省略記法で書いたものがB,Hです。
wは領域積分項です。同様に、領域要素の節点を特定し節点情報に基づいて積分値を得た後、それを全ての要素について足します。
自分は、このプログラムの実行者は「要素」だと思います。
境界要素や領域要素は具体的には節点で定義されます。しかし図-1,2に示すように、同じ節点が隣の要素(一つとは限らない)に共有されたりします。従って要素-節点対応表が必要です。こんな具合です。
※ 領域要素の番号付けは、具体的な要素発生が不明なので例です。
ここでB_Countは境界要素数,BN_Countは境界節点数。R_Countは領域要素数,RN_Countは領域節点数。これらは[入力部]で具体的に要素データと節点データが与えられた時に、セットされるべきものです。表中のjの下の1,2と1,2,3,4は、境界要素では1:始端,2終端を表し、領域要素ではこの順番で左回りに節点を順序付ける事を意味します。
要素-節点対応表はもちろん配列で定義しますが、十進BASICでは可変長配列が使えないので、最初に配列寸法の上限を決めておく必要があります。ここでは、配列寸法上限:1000、としておきます。
境界要素に対する要素-節点対応表は、EN_B(k,2)の形にします。
kが境界要素番号,2の並びが先の1,2に対応します。EN_B(k,2)の内容は、(k,1),(k,2)で指定される境界要素kが持つ節点番号です。
同様に領域要素ではEN_R(k,4)の形にします。
kが領域要素番号,4の並びが先の1~4に対応し、(k,j),j=1~4で指定される領域要素kが持つ節点番号を格納します。
節点の具体的情報が必要です。節点データは当然、座標データで入力されます。これも配列で定義すべきです。境界節点座標はBN(j,2),領域節点座標はRN(j,2)とします。jが節点番号,2の並びがx座標,y座標です。
次に未知ベクトルψとqです。式(2)の計算方針では、ψとqを一つの未知ベクトルz=(ψ q)tとして扱います。それが連立方程式の組み方として、最も単純だったからです。
理想的には節点jにおいてψj,qjですが、(ψ q)tとまとめた時点でもうその扱いはできません。ψの成分数は明らかに境界節点数BN_Countと同じです。j→ψj→z(j)です。従ってj→qj→z(BN_Count+j)です。確かにこの規則があれば、各未知数のzの中での位置の特定は簡単です。しかしいずれ[2]を考慮するのでした。その時には、もっと複雑な対応付けに迫られます。節点-未知数対応表も用意しておきましょう。後で手を加えたとしても、基本の形だけは、ここで整えます。
節点-未知数対応表はBN_Z(j,2)とします。
jが節点番号で,(j,1)がψjに対応するz上の成分番号,(j,2)はqjに対応するz上の成分番号です。
式(1)の(D1 D2)と(d1 d2)tブロックは、境界条件が与えられる場所と、その値を示すものでした。これらの情報も[入力部]でセットされるべきものです。境界条件が与えられる場所の情報は、BN_BP(j,2)とします。
jは節点番号で、(j,1)にはψjが与えられる/与えられないで1または0をセットします。(j,2)にはqjが与えられる/与えられないで1または0をセットします。(D1 D2)ブロックで境界条件を与えるべき列番号は、節点番号jからBN_Z(j,1または2)の値でわかります。(D1 D2)の成分に与える値は1です。
境界条件の値の情報は、BN_BV(j,2)とします。BN_BP(j,2)と同様です。BN_BP(j,2)の指定に対応する形で、具体的な境界条件の値を格納します。(d1 d2)tブロックで境界条件の境界値を与えるべき行番号は、jに従って、zの中で(ψ q)tの後に順番に並べるだけです。
とりあえず最外殻のアウトライン決定作業としては、こんなものです。しかしここまででも、かなり沢山の事を決定し、わかった事もあります。だから[メモ欄]に決めた事を書いておくんですよ。平たく言えば、変数定義の説明です。また宣言部では必要とわかった変数(配列)を、この時点でちゃんと宣言します。後で勝手に変なものを作らせないためにです。[入力部]~[出力部]にも、わかった事を出来るだけ具体的に書いておきます。特に実行者が「要素」である事を明確にします。
※後で勝手に変なものを作るのは、じつは本人だったりするんですよね(^^;)。
3.プログラミング言語の選択
本当は、1.と2.を勘案しながら、システムを記述する言語を選択すべきなんですが、途中でちょろっと漏らしたように、まぁ~十進BASICで行けると思います。アマチュアにとって「無料」は大事です(^^)。
REM [メモ欄]
REM B_Count:境界要素数
REM EN_B(k,2):境界要素-節点対応表,k=1~B_Count
REM (k,1):要素kの始端の節点No.,(k,2):終端の節点No.
REM R_Count:領域要素数,Region Elements Countの略
REM EN_R(k,4):領域要素-節点対応表,k=1~R_Count
REM (k,j),j=1~4:jに従って要素kの要素節点No.を左回りに格納
REM BN_Count:境界節点数
REM BN(j,2):境界節点座標,j=1~BN_Count,(j,1):節点jのx座標,(j,2):y座標
REM RN_Count:領域節点数
REM RN(j,2):領域節点座標,j=1~RN_Count,(j,1):節点jのx座標,(j,2):y座標
REM BN_BP(j,2):節点の境界条件の有無,j=1~BN_Count
REM (j,1):節点jのψjが既知なら1、そうでないなら0
REM (j,2):節点jのqjが既知なら1、そうでないなら0
REM BN_BV(j,2):境界条件の値,j=1~BN_Count
REM (j,1):節点jのψjに与える値
REM (j,2):節点jのqjに与える値
REM z(i):未知ベクトル,i=1~2*BN_Count
REM BN_Z(j,2):節点-未知数対応表,j=1~BN_Count
REM (j,1):ψjのzの中の位置,(j,2):qjのzの中の位置
REM 節点j→ψj→z(j),節点j→qj→z(BN_Count+j)
REM [宣言部]
REM 配列寸法の上限は1000
REM B_Count:境界要素数
REM R_Count:領域要素数
Dim EN_B(1000,2) REM 境界要素-節点対応表,k=1~B_Count
Dim EN_R(1000,4) REM 領域要素-節点対応表,k=1~R_Count
REM BN_Count:境界節点数
REM RN_Count:領域節点数
Dim BN(1000,2) REM 境界節点座標,j=1~BN_Count
Dim RN(1000,2) REM 領域節点座標,j=1~RN_Count
Dim BN_BP(1000,2) REM 節点の境界条件の有無,j=1~BN_Count
Dim BN_BV(1000,2) REM 境界条件の値,j=1~BN_Count
Dim BN_Z(1000,2) REM 節点-未知数対応表,j=1~BN_Count
Dim z(1000) REM 未知ベクトル,i=1~2*BN_Count
REM [入力部]
REM境界要素数B_Countと領域要素数R_Countは、事前に与えられると仮定する。
REM 境界要素と領域要素データは、要素ごとに節点番号が与えられると仮定する。
REM 与えられる節点番号は、左回りに順序付けられたものが与えられると仮定する。
REM 境界要素を要素番号kで読み込むLoop
REM Loop内で要素ごとの節点番号を読み、境界要素-節点対応表EN_B(k,2)をセット
REM 領域要素を要素番号kで読み込むLoop
REM Loop内で要素ごとの節点番号を読み、領域要素-節点対応表EN_R(k,4)をセット
REM境界節点数BN_Countと領域節点数RN_Countは、事前に与えられると仮定する。
REM 境界節点を節点番号jで読み込むLoop
REM Loop内で境界節点座標を読み、BN(j,2)をセット
REM 領域節点を節点番号jで読み込むLoop
REM Loop内で領域節点座標を読み、RN(j,2)をセット
REM 境界条件の場所を節点番号jで読み込むLoop
REM Loop内で節点jの境界条件の有無を読み、BN_BP(j,2)をセット
REM 境界条件の値を節点番号jで読み込むLoop
REM Loop内で節点jの境界条件の値を読み、BN_BV(j,2)をセット
REM 節点番号jのLoopで、節点-未知数対応表BN_Z(j,2)をセット
REM (j,1):ψjのzの中の位置,(j,2):qjのzの中の位置
REM 対応は、節点j→ψj→z(j),節点j→qj→z(BN_Count+j)
REM [計算部]
REM 境界要素番号kで境界積分を行うLoop
REM 要素番号kからEN_B(k,2)に従い節点を特定し、要素ごとに境界積分を実行
REM 上記結果に基づき、係数マトリックスBとHを作成
REM 領域要素番号kで領域積分を行うLoop
REM 要素番号kからEN_R(k,4)に従い節点を特定し、要素ごとに領域積分を実行
REM 上記結果に基づき、既知ベクトルwを作成
REM 境界節点番号jで境界条件位置を指定するLoop
REM 節点番号jからBN_BP(j,2)とBN_Z(j,2) に従い、(D1 D2)ブロックを作成
REM 節点番号jに従いBN_BV(j,2)の値をwブロックの後に並べ、(d1 d2)ブロックを作成
REM 全体係数行列が出来たので、Gaussの掃き出し法にかけ、未知ベクトルzを求める
REM [出力部]
REM BN_Z(j,2) に従いz(j)から値を読み出し、出力
(執筆 ddt³さん)
第19回 ひずみテンソルと応力テンソルの関係 [ネコ騙し数学]
第19回 ひずみテンソルと応力テンソルの関係
§1 ひずみテンソルと応力テンソルの関係
弾性体が変形し、その相対的な変位がvであるとする。このとき、ひずみテンソルΦの成分は
となる。
等方的な弾性体では、応力テンソルΨとひずみテンソルΦの主軸は同じと考えられるので、ΦとΨの共通な主軸の方向の長さ1のベクトルをa、b、cとする。ひずみテンソルの主値と応力テンソルの主値をそれぞれε₁、ε₂、ε₃、σ₁、σ₂、σ₃としたとき、
という関係が成り立つものとする。
等方的なので、
と考えられので、(1)式を
と変形することができる。
ここで、
とすれば、(2)式は
また、
は不変量で、体積膨張率である。
したがって、
である。
同様に、
つまり、
である。
このλ、μをLamé(ラメ)の定数と呼ぶ。
また、応力テンソルΨは
これに(3)式を代入すると、
単位テンソルをIとすると、
また、
だから、
ひずみテンソルの成分を、応力テンソルの成分をとすると
となる。
テンソル積と単位テンソルを行列で表すと、
だから、
また、
§2 弾性体の釣り合い
弾性体内の各点の密度ρに比例した力ρKが単位体積あたりに働いているとする。さらに、弾性体の変形にともなう応力テンソルをΨとする。
弾性体内の任意の領域をVとし、その表面をS、Sの外法線ベクトルをnとする。
このとき、力の釣り合いは
(テンソルの)ガウスの発散定理から
Vは任意に選べるので、
となる。
これが弾性体の釣り合いの方程式である。
Kの成分を、Ψの成分をとすると、
特に、K=0のとき、
つまり、
である。
弾性体が等方的であるとき、
となるので、よって、(6)式は
となる。