［バーガーズ方程式_2］

［バーガーズ方程式_2］

 

　前回は、準バーガーズ方程式（？）、

(1)

 

に対して座標変換、

(2)

 

を行い、

 

(3)

 

の形を得ました。ここでcは適当な定数係数、νは動粘性率です。そういう訳で(3)を、有限要素法へ持ち込みます。個人的には陽解法大好き星人なので(1)を直接ルンゲクッタなんかで解きたいところなのですが、「有限要素法にしろぉ～！」というお達しがネコ先生からあったので(^^;)。

 

４．解析領域に関する考察

　(1)のtとxは、時間と1次元の拡がりを表す空間変数です。つまり(1)は、一次元問題です。そうなんですが、時間tまで含めて有限要素法に持ち込むためには、まるで相対性理論の四次元時空のように(t，x)時空間を考えます(^^;)。

 

 

 

　(1)は時間を含んだ偏微分方程式なので、初期条件が必要です。図-1に示すように時刻t＝0におけるx方向へのuの分布u(0，x)が初期条件になります。x方向への解析領域の拡がりはLとしますが、図-1でu(0，x)はLの半分に分布しています。これが時刻t＝0 → t＝τ → t＝2τ ・・・ t＝nτと時間が進むに従い、どういう分布に変わるのかを（どのようにx方向へ伝わるかを）、時間ステップ幅τで追跡するのが目的です。

　偏微分方程式の解を一意に定めるためには、空間方向の端で境界条件を与える必要もあります。ここでは簡単のため、任意の時刻tでu(t，0)＝u(t，L)＝0とします。

　いま図-1を(2)で変換した場合、(t，x)系のt＝nτとx＝mhのラインは、

 

(4)

 

になります。ここでhは図-1に示すように、Lの分割幅です。(4)から図-1は図-2に変換されます。

 

 

 

 

　初期条件はu(0，η)。境界条件は、η＝－cξとη＝－cξ＋Lで、u(ξ，η)＝0です。

 

 

５．陽解法と陰解法の違い

　時間を含む微分方程式には因果律が成り立つので、陽解法は、図-2に赤で示したある時刻ξ＝nτ（t＝nτ）での格子点(ξ，η)でのu(ξ，η)の値を初期値として、後ろをいっさいふりかえる事無く、ξ＝(n＋1)τでの格子点（青）の値を差分計算などで求めます。この操作を本当の初期条件ξ＝0でのu(ξ，η)から始め、必要な時間tまで繰り返します。陽解法は直前のデータしか参照しないので、一般に計算手続きが簡単になり計算も速いです。しかし数値誤差が蓄積しやすいという欠点を持ちます。後ろをいっさいふり返らないからです。でも4次のルンゲクッタ法などは非常に高精度ですから、実用的には十分なケースもままあります。

　陰解法は最低でも1ステップ前までふり返り、ふり返った時間領域にあるすべての格子点で平均的に微分方程式が満たされるように調整します。その結果、数値誤差が抑えられ解は長持ちするのですが、調整計算が複雑で計算速度も遅くなりがちです。現実には問題に応じて使い分けます。ここではせっかくなので、ξ＝0～tの全ての格子点を計算対象とする陰解法をめざします。

 

６．重み付き残差法による弱形式

　ここで、関数解析などをやるわけにはいかないので、次の定理は認めてください。

 

［定理-1］

　{w_n(x)}， n＝1，2，・・・を閉区間[a，b]で少なくとも連続な関数列の完全系とする。任意のnで、

 

(5)

 

が成り立つなら閉区間[a，b]でf＝0。

 

　閉区間[a，b]上の関数列の完全系とは、[a，b]上の任意の関数f(x)を、

 

の形で展開可能な関数列の事を言います。ここにa_n， n＝1，2，・・・はxに無関係な定数係数です。{w_n(x)}は少なくとも連続としましたから、任意の関数とは任意の連続関数の事だと思ってもらってけっこうです。{w_n(x)}の実例はすぐ与えます。w_n(x)の事を重み関数と言います。

　［定理-1］により、次の重み付き残差法が可能になります。(3)に対して、

 

(6)

 

が任意のn＝1，2，・・・で成り立ちます。積分領域Ｒは閉領域です。

　(5)では1次元の関数だったw_n(x)が、(6)ではいきなり2次元のw_n(ξ，η)になっちゃいましたが、(6)の2次元の重積分は累次積分で、1次元のξとηの積分に直せるので［定理-1］が使えます。つまり任意のnに対して(6)が成り立つようにu(ξ，η)を決められたら、それは(3)の解だという事です。

　一般に微分方程式の微分の階数が高いほど、その数値解は暴れる君になります。そこで(6)をガウスの発散定理で部分積分し、微分の階数を減らしておきます。

(7)

 

　Ｓは積分領域Ｒの境界、dsはＳの外法線線素ベクトル、・は内積です。(7)の形を、重み付き残差法における弱形式といいます。あんまり使いませんけど、(6)は強形式という事になります。

 

７．要素の設計

　現在最も使い勝手の良い有限要素法は、重み付き残差法における弱形式に基づく定式化です。数値解法一般の発想は基本的には同じです。「どんな関数も局所的には多項式表現できる」です。有限要素法でもそう考えますが、有限要素法における局所の範囲は、図-2に示したような要素です。一個だけ取り出すと図-3になります。

 

 

　要素の角はたいてい節点と呼ばれ、左回りに順序付けるのが普通です。節点1～4におけるu(ξ，η)の関数値u₁～u₄で要素内のu(ξ，η)の値を補間します。補間関数には「どんな関数も局所的には多項式表現できる」から、多項式が用いられます。さらに要素が図-3のような平行四辺形の場合、計算に便利なように、アイソパラメトリック座標と呼ばれる要素辺に沿った局所座標系(α，β)が使われます。

　平行四辺形の場合、局所座標(α，β)を適当に定めると、要素は辺長1の正方形として表せます。明らかに、

(8)

 

とすればOKです。(α，β)∈[0，1]×[0，1]。次に補間関数を決めます。補間関数はu₁～u₄の値で決まるα，βの多項式なので、4個以上の未定定数を持つ多項式である必要があります。1次補間だと未定定数は3個なので、完全2次多項式を考えます。

(9)

 

　そうすると未定定数は6個なので、係数a₁～a₆が2個余ります。それで普通はa₄＝a₆＝0とします。a₁～a₃を残すのは、いわゆる剛体変位と剛体回転を表現できるようにするためです。それが可能でないと、例えばu＝一定などの自明解を表現できないからです。よってa₄～a₆の中から2個まびく事になります。

　もしa₄＝a₅＝0とするとβ²だけが残り、β方向を偏重した要素になりますが、物理的にαよりβを優遇する理由はありません。そのような要素は癖が悪いわけです(^^;)。それでα，βに関する対称性を保つために、α²とβ²の項をまびき、αβの項を残します。不完全2次多項式、

 

(10)

 

 

の補間関数を持つ要素を双一次要素と言います。α＝一定またはβ＝一定だと(10)は一次関数になるからです。

　未定定数a₁～a₄は、(α，β)系で節点1～4の位置が、(α，β)＝(0，0)，(1，0)，(1，1)，(0，1)である事から、それらの座標を(10)に代入し、

 

(11)

 

 

で与えられます。(11)は3行目と4行目を交換すれば、

 

 

(12)

 

となるので、a₁～a₄は容易に計算できます。これがアイソパラメトリック座標の効果です。(11)または(12)より、

 

 

(13)

 

 

　従って補間多項式(10)は、

(14)

 

 

と書けます。こうして要素内でu(ξ，η)を多項式表現するという当初の目的は達せられましたが、じつは要素ので設計はこれだけでは済みません。

 

８．適合要素と非適合要素

　今まで未知関数を一要素上で多項式表現する事だけを考えてきました。一要素の（局所の）補間多項式が定まると、弱形式(7)を利用して、要素節点の関数値が数値解法上満たすべき方程式を出せます。しかし図-3に示すように、考慮すべき要素は一個だけではありません。有限要素法で答えを得るためには、解析領域内の全ての要素に対する全体系を組み上げる必要があるわけです。そのとき問題になるのが、隣の要素との接続条件です。

　6.で述べたように、未知関数u(ξ，η)は少なくとも連続なのが重み付き残差法の前提です。よって図-4左側に示すように隣り合う要素は、その接続辺において同じu(ξ，η)の値を持たねばなりません。

 

 

　しかし要素ごとに補間関数を考えるという事は、図-4右側のように接続された要素を分解して考えてる事になります。結果、接続辺において関数値u(ξ，η)が二重化されます。u(ξ，η)は少なくとも連続なので、少なくとも接続辺上の節点関数値は等しくなければなりません。図-4右側の記号を用いれば、

 

(14)

 

です。

　このような条件は、全体系の方程式を作る際に取り込まれます。しかしu₂⁽¹⁾とu₃⁽¹⁾の間のuの関数値と、u₁⁽²⁾とu₄⁽²⁾ 間のuの関数値も等しいことは保証されるでしょうか？。それは補間関数の決め方によります。

　双一次要素では、α＝一定またはβ＝一定でu(α，β)は一次関数になるのでした。上記2つの要素の接続辺は、α＝1とα＝0の辺です。よって(14)が成り立てば、一次関数は両端点で値が決まるので、辺上でuの値が一致する事が保証されます。このような要素を適合要素と言います。

　多項式表現は適合要素を比較的作りやすく、いわゆる剛体変位と剛体回転を表せる事も保証されるので、そういう理由もあって愛用されるのですが、多項式表現は次数によって未定定数の数が決まっています。それは適合性とか、要素形状に合わせて要素上に配置する節点の数などとは無関係な条件なので、要素設計において全ての条件を満たすには難儀する場合もあります。それで余りにも面倒な場合は、非適合要素を採用する場合もあり得ます。非適合要素の場合、接続辺上の関数値の一致は要素が十分細かく分割される事を前提に、そうみなされます。辺長が短ければ辺上では、どうせ一次関数みたくなるよね？。だから端点の節点関数値の一致さえ要求しとけば大丈夫さ・・・、となります(^^;)。

 

　(14)を弱形式(7)に代入し、重み関数を具体的に定め積分を実行すれば、一要素の節点関数値uに関する要素方程式と要素マトリックスが得られるのですが、けっこう長い計算になる上に、一気にやるしかない計算でもあります。という訳で［バーガーズ方程式_3］に続きます。