一般逆行列理論の流体問題への応用―Lagrange 記述と Euler 記述による手法―
横須賀洋平*1
松澤照男*2
梗

概本報では，不安定次数の高い構造物の挙動を示す形態解析手法に有用な一般逆行列による手法を流体問題へ応用し，Lagrange 記述と Euler 記述による解析事例を紹介する．Lagrange 記述による手法では，膜のポンディング現象を想定する．流体を理想流体としてモデル化し，膜面上で自重による形状変化を起こす流体と釣合う膜形状を求める．Euler 記述による手法では，膜と風の連成問題など流体の解析に注目がおかれる中，膜構造の解析に親和性のある一般逆行列による非圧縮性流体の非定常粘性流れにおける流れ場と圧力場の計算手法を提案する．

1．はじめに
一般逆行列による形態解析手法の特徴は，線形連立方程式の解が不定あるいは不能となる長方形行列を係数とする方程式に対して，一般逆行列を用いて最小二乗解やノルム最小解といった意味のある解を与える[1,2]．また，方程式の一般解における補解を用いることで，方程式の右辺ベクトルで表される量を 0とする解を提供することが可能となる．既往の研究では，ひずみや伸び、伸び速度あるいは外力ベクトルなどを 0 とする形態解析手法が多く提案されてきたが、本報では流体解析に用いられ，非圧縮性を示す連続の式を満たす速度を求める際に，一般逆行列を用いた手法を提案する．本報では，Lagrange 記述とEuler 記述によるそれぞれの手法を提案する．Lagrange 記述による手法では，膜のポンディング現象を想定した膜と流体の釣合形状を求める数値解析手法を提案する．Euler 記述による手法では，非圧縮性粘性流体の非定常流れ場，圧力場を求める数値解析手法を提案する．
2．Lagrange 記述による手法
1) ポンディング現象について
張力膜構造では，雨水や積雪，融雪水が貯水となるポンディング現象が起き，大きな荷重に進展する可能性が指摘されている．その発生メカニズムについては，文献[3]等を参照されたい．
ポンディング現象の興味深い点は，貯水による膜面にかかる荷重分布，荷重分布と釣合う膜形状が密接な関係にある．この関係性は下に凸の膜形状のライズスパン比が大きくなるように変位すれば，貯水領域が狭まる，すなわち荷重分布が変化し，また荷重分布に応じて膜形状も再度変化するといった膜と流体*1*2が不可分な関係にある．さらに，膜構造は面外の曲げ剛性が低いため，荷重分布の変化に応じて剛体変位を含む不安定挙動も生じると考えられ，弾性変形，応力弛緩，クリープ，材料非線形性も考慮すると極めて複雑な現象となる．
本報ではポンディング現象を流体と膜の連成問題と捉え，その基本性状を把握することを目的とする．そこで貯水領域を解析モデル化し，不安定挙動が生じる膜との連成解析手法を提案する．既往の研究では，貯水体積から膜面にかかる荷重分布を求める必要があるが，連成問題として解く場合は，自然と流体形状が膜に適合し，膜は流体の自重により形状が決定されると考える．複雑な膜曲面においても本手法は有効と思われる．本手法では，基本性状を把握するため，解析を簡略化し，以下のような仮定を与える．
a) 解析空間を二次元平面とする．動的影響を考慮しない準静的な問題と考える．
b) 外力ポテンシャルを停留する過程で，逐次変化過程を求める．釣合式から最終的な釣合形状の判定を行う．
c) 膜はトラス要素とし，弾性変形はせず不安定挙動のみを表す．膜全長が一定の仮定と節点滑動を許容する．
d) 膜の端部に膜の滑動を許容する支持点を設ける．これにより，支点間の膜長さが変化するため，面外変位が生じる．
e) 膜の端点に一定の水平力を与え，膜を引っ張る．流体と膜形状が逐次変化することで，最終的に釣合う水平力となる．
f) 流体は，完全流体の Lagrange 記述とし流体の移流速度

昭和女子大学環境デザイン学科 助教/北陸先端科学技術大学院大学情報科学研究科 博士後期課程
北陸先端科学技術大学院大学 シミュレーション科学研究センター 教授

は要素の移動速度と等価と仮定する．
g) 流体界面の粘性応力と表面張力は無視できると仮定する．
h) 流体要素は三角形一次要素とする．
i) 流体要素の節点に自重に相当する等価な節点外力を与える．
2) 主要変数について連続体力学を基礎とした理論構成の中で，Lagrange 記述において現時刻𝑡0 = 𝑡を基準配置(解析済みの既知)として時刻𝑡0 = 𝑡 + 𝛥𝑡の配置を求める増分解析を updated Lagrange 法と呼ぶ．以下，「Lagrange 記述による手法」において同じ立場をとる．この場合，一般に次の関係式が成り立つ[4 ]．
(1)𝑬̇ = 𝑫𝜕𝑢̇𝜕𝑢̇𝑖𝑗̇ =1
𝑬̇ ≡ 𝐸𝑖𝑗+(2)2 (𝜕𝑥𝑗 𝜕𝑥𝑖 )1 𝜕𝑣𝑖 𝜕𝑣𝑗
𝑫 ≡ 𝐷𝑖𝑗 =+(3)2 (𝜕𝑥𝑗 𝜕𝑥𝑖 )全体座標𝑥𝑖 ，変位速度𝑢̇𝑖 ，速度𝑣𝑖 はそれぞれ現時刻𝑡0 = 𝑡の値を表す． Green-Lagrange 歪速度𝑬̇ は速度形境界値問題で現れる物理量であり，慣例的に物質点を固定した Lagrange 記述に対する問題で扱われ，変形速度テンソル𝑫は空間を固定した Euler 記述に対する問題でよく扱われる.本報では，式(1),(2),(3)は，これらの問題を統一的に扱うひとつの架け橋として考える．
式(2),(3)から明らかなように，現時刻𝑡を基準配置とする場合，𝒖̇ = 𝒗(4)が成り立つが，本手法では，膜の主要変数として変位速度ベクトル𝒖̇を用い，流体の主要変数として速度ベクトル𝒗を用いて表記を区別する．位置ベクトル𝒙を持つ物体形状が𝑡 → 𝑡 + 𝛥𝑡へ形状更新する際に， 𝒖̇𝛥𝑡 = 𝛥𝒖および𝒙𝑡+𝛥𝑡 − 𝒙𝑡 = 𝛥𝒖の関係から(5)
𝒙𝑡+𝛥𝑡 = 𝒙𝑡 + 𝒖̇𝛥𝑡 = 𝒙𝑡 + 𝒗𝛥𝑡で行う．また，本報の「Lagrange 記述による手法」において，時刻𝑡は運動エネルギーの影響を無視できる緩慢に進む時間パラメータとする．
3) 膜要素の定式化
3-1) 基礎式と離散化について
二次元平面を解析空間とし，膜をトラス要素とする．歪速度成分は単軸のみとなり，トラス要素は二次元空間に配する必要がある．このような場合，一般に有限要素法では，局所座標を参照してから全体座標へ座標変換を行うが，本報では，直接全体座標を参照して座標変換を行わない．膜要素において，以下のとおり歪速度と変位速度の関係は示される．
𝜕𝑢̇
𝜀̇ = tr𝑬̇ ≡ 𝐸𝑖𝑖̇ = 𝑖(6)𝜕𝑥𝑖
一般に変形速度テンソル𝑫の対角成分と対数歪𝜀の関係の物理的意味づけはなされている[4]．本報では Green-Lagrange 歪速度𝑬̇ で表しており，対角成分の和をとる意味は，各空間方向に分解して定義される歪速度成分の和をとることを意味する(𝜀̇ = 𝜀̇𝑥1 + 𝜀̇𝑥2 )．
式(1)を基礎式とし，アイソパラメトリック要素を用いて離散化を行う．アイソパラメトリック要素を用いる理由として，一般逆行列を用いれば自然座標から任意の座標系の空間微分値を導出することが可能な点にある．形状関数𝑁 (𝑘) (𝑘 = 1,2) を用いて変位速度と要素内座標を補間する．
11
𝑁 (1) = (1 − 𝑟), 𝑁 (2) = (1 + 𝑟)22(7)𝑢̇𝑖 = 𝑁 (1) 𝑢̇𝑖 + 𝑁 (2) 𝑢̇𝑖(1)(2)(8)(1)(1)(9)𝑥𝑖 = 𝑁 (1) 𝑥𝑖 + 𝑁 (2) 𝑥𝑖

自然座標𝑟 による形状関数の微分値は𝑑𝑁 (1) /𝑑𝑟 = −0.5 ,𝑑𝑁 (2) /𝑑𝑟 = 0.5で与えられる．ここで，全体座標による形状関数の微分値は以下の過程で与えられる．
(𝑘)⎡𝜕𝑁 ⎤𝑑𝑁 (𝑖) (𝑖) ⎢ 𝜕𝑥1 ⎥𝑥∑ 𝑑𝑟 2 ] ⎢𝜕𝑁 (𝑘) ⎥⎢⎥𝑖=1⎣ 𝜕𝑥2 ⎦2𝑑𝑁 (𝑘)𝑑𝑁 (𝑖) (𝑖)=𝑥∑𝑑𝑟𝑑𝑟 1[2𝑖=1(10)(𝑘)⎡𝜕𝑁 ⎤⎢ 𝜕𝑥1 ⎥= 𝑱 ̃ ⎢ (𝑘)⎥⎢𝜕𝑁 ⎥⎣ 𝜕𝑥2 ⎦̃𝑱 は，本来ヤコビ行列に相当するが，ここでは正則行列ではないため ̃ 記号を付す．したがって，ヤコビアンの演算は定義できないが，一般逆行列を用いれば，全体座標による形状関数の微分値は求まる．
(𝑘)⎡𝜕𝑁 ⎤(𝑘)⎢ 𝜕𝑥1 ⎥+ 𝑑𝑁(11)⎢𝜕𝑁 (𝑘) ⎥ = 𝑱 ̃ 𝑑𝑟⎢⎥⎣ 𝜕𝑥2 ⎦𝑱 ̃+ は，𝑱 ̃のムーア・ペンローズ一般逆行列を示す．トラス要素を用いる場合，𝑱 ̃は行ベクトルになるため，𝑱 ̃+ = 𝑱 ̃T /𝑱 ̃𝑱 ̃T で求めることができる．本離散化手法を応用すれば，一次元要素を二次元平面で表現する場合だけでなく，二次元要素を三次元空間で表現する場合にも有効といえる．
(1)

(1)

(2)

(2)

ひとつの要素における変位速度𝒖̇T𝑒 = {𝑢̇𝒙1 , 𝑢̇𝒙2 , 𝑢̇𝒙1 , 𝑢̇𝒙2 }で表し，全体座標による形状関数の微分値を用いて，式(6)を表すと以下のようになる．
𝜕𝑢̇𝑖𝜕≈𝜕𝑥𝑖 [𝜕𝑥1=

𝜕𝑁 (1)[ 𝜕𝑥1𝜕
𝑁 (1)𝜕𝑥2 ] [ 0𝜕𝑁 (1)𝜕𝑥2

𝜕𝑁 (2)
𝜕𝑥1

0𝑁 (1)𝑁 (2)00 𝒖̇
𝑁 (2) ] 𝑒𝜕𝑁 (2)𝒖̇ ≡ 𝑠𝑩𝑒 𝒖̇𝑒 = 𝜀̇𝑒𝜕𝑥2 ] 𝑒(12)

ここで 𝑠𝑩𝑒 マトリクスの下添字𝑠は流体と区別し，構造を意味する．
式(12)を全要素に拡張すれば，下式を得る．
̇𝑠𝑩 𝒅 = 𝜺̇(13)𝑠𝑩 マトリクスは𝑚𝑠 × 𝑛(𝑚𝑠 < 𝑛)，𝑚𝑠 は構造の要素数，𝑛は全体

また，膜端点に作用する水平の外力ベクトルを 𝑠𝒇 ，膜のトラス(連成する場合は流体も含む)の自由度数の長方形行列となる．
̇全体の変位速度𝒅 とし，𝑛次元列ベクトルとする．
本報では，膜に全長一定の仮定と節点の滑動を可能とする条件を与える．そのためには，式(13)では不十分であり，歪速度ではなく，伸び速度 ∆𝐿̇ で定義する必要がある． 𝑑𝐿を現時刻𝑡0 = 𝑡における微小な要素とし，歪速度を要素内で積分する．
𝐿 𝜕𝑢̇𝑖𝑑𝐿 ≈ 𝑠𝑩𝑒 𝒖̇𝑒 𝐿 = 𝑠𝑩𝑒∗ 𝒖̇𝑒 = ∆𝐿̇𝑒
(14)
∫0 𝜕𝑥𝑖式(14)を全要素に拡張すれば，下式を得る．

要素に働く軸力を𝑁とすれば，以下の関係が成り立つ．

(15)

∗ ̇̇
𝑠𝑩 𝒅 = ∆𝑳𝑠𝑩

速度形の関係式において上式が成立する理由として，一般に定義されている𝑩マトリクスは速度形も変位形も変わらない点にある． 𝑠𝑩 ∗∗ マトリクスにおいて，膜全体をひとつの要素と捉えているため，軸力𝑁がスカラーで表されることになる．
ここで，全長一定の仮定を与えるために，式(18)を用いて，𝑠𝑩∑∆𝐿̇′𝑖𝑘∗(2)𝑠 𝑩1∗(1)+ 𝑠 𝑩2(1)(2)

0𝑁 (2)𝑁 (3)0(2)(3)

(3)

∗(2)

𝜕𝑣𝑖
≈
𝜕𝑥𝑖𝜕
[𝜕𝑥1=𝜕𝑁 (1)𝜕𝑥2 ] [ 0𝜕𝑁 (1)[ 𝜕𝑥1𝜕𝑁 (1)𝜕𝑥2

0𝑁 (1)𝑁 (2)0𝜕𝑁 (2)
𝜕𝑥1

𝜕𝑁 (2)
𝜕𝑥2

𝜕𝑁 (3)
𝜕𝑥10 𝒗
𝑁 (3) ] 𝑒𝜕𝑁 (3)𝒗
𝜕𝑥2 ] 𝑒(22)≡ 𝑓 𝑩𝑒 𝒗𝑒 = 0∆𝐿̇′𝑖 =𝑖=1

∗(1)

(1)

ひとつの要素の速度𝒗T𝑒 = {𝑣𝑥1 , 𝑣𝑥2 , 𝑣𝑥1 , 𝑣𝑥2 , 𝑣𝑥1 , 𝑣𝑥2 }で表し，式(21)を表すと，式(14)を用いて，ひとつの小領域𝑑𝐿′ の伸び速度を表すと，以下となる．小領域の長さ𝐿′ で積分する．
𝐿′ 𝜕𝑢̇𝑖(17)𝑑𝐿′ ≈ 𝑠𝑩𝑒 𝒖̇𝑒 𝐿′ = 𝑠𝑩𝑒∗ 𝒖̇𝑒 = ∆𝐿̇′𝑒∫𝜕𝑥0𝑖
式(16)，(17)を用いて，全体(要素としてはひとつ)に拡張すると，膜要素数𝑚𝑠 = 1となり，一本の方程式が成立する．

[ 𝑠 𝑩1(20)𝒅̇ = 04) 流体要素の定式化
本報では，非圧縮性の連続の式を満たす速度を与えることで，流体要素が移動することを仮定する．したがって，連続の式を基礎式として離散化方程式を与える．密度変化のない非圧縮性流体の連続の式は，以下のとおり表すことができる．
𝜕𝑣∇ ∙ 𝒗 ≡ tr𝑳 ≡ 𝑖 = tr𝑫 = 0(21)𝜕𝑥𝑖
ここで速度勾配テンソル𝑳と変形速度テンソル𝑫の対角成分の和は同じ値となることを用いている．上式は式(6)とよく似た形となる．三角形一次要素を用いて，形状関数により離散化を行う．
一般的に有限要素法で用いられる離散化と同様となる．

(16)

𝑖=1

∑∗∗とおく．この方程式は流体要素で定式化される方程式と連立して解かれる．

𝑘∆𝐿̇𝑒 =(19)𝑁 = 𝑠𝒇

外力 𝑠𝒇 は𝑛次元列ベクトルとする．

上記の定式化はいずれも三次元空間を解析空間とする場合も同様に定式化できる．
3-2) 全長一定の仮定と節点の滑動
前節では膜をトラス要素とし，伸び速度と変位速度の関係式について離散化方程式を導出した．膜と流体の連成を成立させるために，流体要素の形状が自由に変形するには，膜要素も流体の節点を追従する必要がある．この場合，膜要素の節点が自由に動いてよいと仮定すると同時に，膜要素全体の長さは一定であるという仮定を与える．前者を節点の滑動と呼び，後者を全長一定の仮定と呼ぶことにする．両者を同時に満たすために，以下のように考える．
膜全体をひとつの要素と捉える．この要素をさらに複数の小領域に分割できるとする．与えられた分割点を節点とし，自由度を与える．ひとつの要素(膜全体)の伸び速度∆𝐿̇𝑒 と分割した小領域の伸び速度∆𝐿̇′𝑒 の関係は以下のようになる．

∗∗T

∗(1)

⋯ 𝑠𝑩𝑘−1 + 𝑠𝑩𝑘∗(2)

̇𝑠 𝑩𝑘 ] 𝒅(18)≡ 𝑠𝑩 ∗∗ 𝒅 ̇ = ∆𝐿̇𝑒
上式は構造系が開いた形状を想定しているが，閉じた形状の
∗(1)
∗(2)
場合は， 𝑠𝑩1 + 𝑠𝑩𝑘 となる．また，分割した小領域は自由に動いてよいことになるため，小領域の長さ𝐿′ が 0 とならないよ
うに注意を要する．式(18)を膜の離散化方程式と考える．

ここで 𝑓 𝑩𝑒 マトリクスの下添字𝑓 は流体を意味する．前節と同様にアイソパラメトリック要素を用いて離散化することも可能である．
その場合，𝑱 ̃は正則行列となる．
全要素に拡張すれば，下式を得る．
̇𝑓 𝑩𝒅 = 𝟎(23)

膜要素と流体要素の節点が共有されている節点も存在するため，
̇速度𝒅 は構造と統一した記述としている．流体自重
流体要素
膜要素

固定支持点

水平力

膜の滑動を許容する支持点

Fig.1 解析条件

𝑓 𝑩 マトリクスは𝑚𝑓 × 𝑛(𝑚𝑓 < 𝑛)，𝑚𝑓 は流体要素数，𝑛は全体̇(連成する場合は構造を含む)の自由度数の長方形行列，𝒅 は膜の場合と同様，𝑛次元列ベクトルとする．
流体には，重力による物体力が働くとして，鉛直下向きの外力ベクトルが働く．時々刻々と変化する各流体要素の体積比に応じて，節点に等価な外力ベクトル 𝑓 𝒇 を与える．流体要素に一定6) 数値解析事例
これまでの定式化を基に Fig.1 に示された解析条件により，数値解析を行った． Fig.2 に解析結果を示す．Step=5000 でほぼ釣合状態となっている．水平力の大きさを変えることによって，膜の垂れ下がり具合は調整できる．Fig.3 のグラフは釣合式の収束判定を行っている．式(28)で表される量のノルムを計算している．これを釣合誤差ノルムと呼ぶことにする．の圧力𝑝が働くと想定すると，以下の関係が成り立つ．
𝑓𝑩T

𝒑 = 𝑓𝒇(24)

Step=0

圧力は𝑚𝑓 次元列ベクトル，外力は𝑛次元列ベクトルである．上式の左辺は圧力の勾配を意味する．したがって，上式は静止流体の釣合式を意味する．

Step=250

Step=7505)

膜と流体の連成
流体と膜の連成を定式化する．「Lagrange 記述による手法」
で定義している流体要素は，連続の式を満たす速度を導出できるが，運動方程式を満たすとは限らない．本報では，膜と流体の釣合形状が得られればよく，移動中，運動方程式を満たす必要はないと考える．したがって，最終的に式(19)，(24)で表される釣合式を満足することで釣合形状が得られると考える．
上記の定式化から式(20)，(23)を連立して連成式を与える．
𝑓𝑩𝒅 ̇ ≡ 𝑚𝑖𝑥𝑩 𝒅 ̇ = 𝟎[ 𝑠𝑩 ∗∗ ]Step=1000

Step=1250

Step=2000(25)
Step=5000

ここで，一般逆行列を用いて，上式を満たす解は，以下のように表される．
𝒅 ̇ = [𝑰𝑛 − 𝑚𝑖𝑥𝑩 + 𝑚𝑖𝑥𝑩 ]𝜶 ≡ 𝑯𝜶(26)

Fig.2 数値解析結果

𝑰𝑛 は𝑛 × 𝑛，𝑛は自由度数の単位行列，𝑚𝑖𝑥𝑩 + はムーア・ペンロ++𝒑T
T[𝑁 ] = ( 𝑚𝑖𝑥𝑩 ) [ 𝑓 𝒇 + 𝑠𝒇 ] ≡ ( 𝑚𝑖𝑥𝑩 ) 𝑚𝑖𝑥𝒇

(27)

本報では，釣合式を用いて収束判定を行う．以下の式を満たすことで流体と膜の釣合形状が得られることを意味する．
𝒑T
T +𝑚𝑖𝑥𝑩 [𝑁 ] − 𝑚𝑖𝑥𝒇 = 𝑚𝑖𝑥𝑩 ( 𝑚𝑖𝑥𝑩 ) 𝑚𝑖𝑥𝒇 − 𝑚𝑖𝑥𝒇 ≈ 𝟎
T(28)

釣合誤差ノルムーズ一般逆行列，𝜶は任意の𝑛次元列ベクトルを表す．𝜶は外力ポテンシャルを停留する過程において，膜と流体の形状を逐次変化させながら，求めることができる．詳細は[1]等に詳しい．
流体の圧力と膜を置換したトラス要素における軸力は以下のように求めることができる．ただし，移動中誤差を含むと考える．

1.510.501000

20003000
Step4000

Fig.3 釣合誤差ノルム

50006000

3．Euler 記述による手法
1) 非圧縮性流体の非定常粘性流れについて
非圧縮性流体の非定常粘性流れの解析においては Navier－Stokes 方程式と連続の式をそれぞれ満たして解く必要がある．これらは速度や圧力が未知変数として記述される．しかし，速度計算に用いる式をNavier－Stokes方程式，圧力計算に用いる式を連続の式とすると，連続の式には圧力が陽に含まれていないため，困難を要する．この困難を克服するため，非圧縮性流体の解析に様々な工夫が施されたアルゴリズムが提案されており，ペナルティ関数法などは連続の式に圧力を陽に含める工夫を施したものといえる．一般に，よく用いられている手法といえば，分離解法である MAC 法の流れを汲むものといえる．
分離解法は速度と圧力を分離し，順番に解くことで行列の規模を小さくするという考えがある．この場合，方程式をうまく組み合わせて圧力に関する Poisson 方程式を導き，圧力を求め，速度を計算するアルゴリズムとなる．
上記のような流体解析手法が数多く提案されている中で，一般逆行列理論を流体問題に応用するにあたり，非圧縮性の連続の式に注目する．Navier－Stokes 方程式と連続の式をそれぞれ満足する必要があり，分離解法では圧力に関する Poisson 方程式を解くことで連続の式を満たしている．一般逆行列理論によれば，連続の式のように方程式の右辺が 0 となる速度の解を提供する可能性を持つ．これは形態解析手法で言うならば，剛体変位モードと同様となる連続の式を満たす速度モードを抽出できる．また，圧力は連続の式の係数行列の一般逆行列を用いて表現することが可能であり，この速度モードと圧力を用いて，Navier－Stokes 方程式を表すことで，両式を満たす速度が得られることになる．
本手法は，毎ステップ Poisson 方程式を収束計算により解く分離解法と比較して，毎ステップ代入計算のみで求めることができるため，効率的といえる．分離解法の場合，Poisson 方程式を解く際，誤差判定による収束計算が必要となり，計算速度を優先すると，誤差が大きくなる．本手法は，厳密に連続の式を満たすため，精度が高いと思われる．一方で，本手法は速度モードを算出する段階で，一般逆行列演算が必要となる．計算機メモリの負担が大きくなるが，Euler 記述であれば，一回の逆行列演算で済み，この逆行列は，境界条件や初期条件とは独立した格子情報のみに依存するため，再利用することも可能である．
2) 支配方程式と空間方向離散化について
非圧縮性流体の非定常粘性流れの計算に用いられるNavier－Stokes 方程式と連続の式は以下のとおりである．
𝜕𝒗
1+ (𝒗 ∙ ∇)𝒗 = − ∇𝒑 + 𝜈∇2 𝒗 + 𝑿(29)𝜕𝑡𝜌∇∙𝒗=𝟎(30)
ここで，𝜈 = 𝜇/𝜌は動粘性係数である．本手法では，空間を固定した Euler 記述とし，スタガード格子を用いて差分法により離散化を行う．この場合，一般に有限要素法と比較して圧力振動が起きないとされている．また，差分法は要素毎に方程式を満たすため，全体行列を作成する必要がないが，本手法においては全体行列を作成する．
解析空間を二次元平面とし，𝑥1 = 𝑥, 𝑥2 = 𝑦，𝑣1 = 𝑢, 𝑣2 = 𝑣
1111とする．また，速度ベクトルを𝒗T𝑒 = {𝑢𝑖−2,𝑗 , 𝑢𝑖+2,𝑗 , 𝑣𝑖,𝑗−2 , 𝑣𝑖,𝑗+2 }11とする．差分近似として𝑥方向に𝜕𝑢/𝜕𝑥 ≈ (𝑢𝑖+2,𝑗 − 𝑢𝑖−2,𝑗 )/∆𝑥を与える．𝑦方向も同様となる．
連続の式に基づくひとつの格子の離散化式は以下となる．
𝜕𝑣
1111−∇ ∙ 𝒗 ≡ 𝑖 ≈ [−𝒗∆𝑦 ∆𝑦] 𝑒𝜕𝑥𝑖
∆𝑥 ∆𝑥(31)≡ 𝑑𝑚𝑩𝑒 𝒗𝑒∆𝑥, ∆𝑦はひとつの微小な格子幅を表す．𝑑𝑚𝑩𝑒 の添字𝑑𝑚は前節の𝑩マトリクスと区別して差分法を意味する．格子全体に拡張
すると，
𝑑𝑚

𝑩𝒗 = 𝟎(32)𝑑𝑚

𝑩 マトリクスは𝑚 × 𝑛(𝑚 < 𝑛)，𝑚は格子数，𝑛は自由度数の長方形行列となる．以下，添字𝑑𝑚は省略して表記する．ここで，一般逆行列を用いて，上式を満たす解は，以下のように表される．
𝒗 = [𝑰𝑛 − 𝑩 + 𝑩 ]𝜶 ≡ 𝑯𝜶(33)
𝑰𝑛 は𝑛 × 𝑛，𝑛は自由度数の単位行列，𝑩 + はムーア・ペンローズ
一般逆行列，𝑯は連続の式を満たす速度モード，𝜶は任意の列ベクトルを表す．𝑩マトリクスを用いて圧力の勾配は下式となる．
(34)∇𝒑 = −𝑩 T 𝒑3)

時間方向離散化について時間方向の離散化は SMAC 法にならって行う．ここで，速度𝒗𝑛+1 の右肩の添字は𝑡𝑛+1 時間ステップを表し，移流項，拡散項の係数行列を𝑨(𝒗𝑛 ), 𝑫 とする．式(34)を用いると Navier－Stokes 方程式は下式となる．
1𝒗𝑛+1 − 𝒗𝑛
(35)
+ 𝑨(𝒗𝑛 )𝒗𝑛 = −(−𝑩 T 𝒑𝑛+1 ) +𝑫𝒗𝑛
𝑅𝑒∆𝑡
𝑅𝑒はレイノルズ数である．𝑛ステップは既知量とすると，𝒑𝑛+1 が未知であるため，𝒑𝑛 と置き換え，
1𝒗̃ − 𝒗𝑛(36)+ 𝑨(𝒗𝑛 )𝒗𝑛 = −(−𝑩 T 𝒑𝑛 ) +𝑫𝒗𝑛
𝑅𝑒∆𝑡𝒗̃は𝒗𝑛+1 の予測値となる．式(35)，(36)の辺々を引くと，下式となる．
𝒗𝑛+1 − 𝒗̃(37)= 𝑩 T (𝒑𝑛+1 − 𝒑𝑛 )∆𝑡
本来 SMAC 法では，上式から𝒑𝑛+1 − 𝒑𝑛 = 𝝋と置いて，𝒗𝑛+1 について整理し，連続の式に代入することで圧力の修正量𝝋のPoisson 方程式を導く．本報では，𝒗𝑛+1 = 𝑯𝜶と置くことで，連続の式を満たす速度を抽出できることになる．したがって，式(33)，(35)を用いて，𝒑𝑛+1 を表すと，
1𝑯𝜶 − 𝒗𝑛
𝒑𝑛+1 = (𝑩 + )T (+ 𝑨(𝒗𝑛 )𝒗𝑛 −𝑫𝒗𝑛 )
𝑅𝑒
∆𝑡

(38)

+

ここで，(𝑩 T ) = (𝑩 + )T の関係を用いている．
式(33)，(38)を式(37)に代入して，
𝑯𝜶 − 𝒗̃1
𝑯𝜶 − 𝒗𝑛= 𝑩 T (𝑩 + )T (+ 𝑨(𝒗𝑛 )𝒗𝑛 −𝑫𝒗𝑛 )
∆𝑡
∆𝑡
𝑅𝑒

(39)

T 𝑛−𝑩 𝒑

上式を整理すると，
T+ T(𝑯 − 𝑩 (𝑩 ) 𝑯 )𝜶T= 𝒗̃ + 𝑩 T (𝑩 + )(−𝒗𝑛 + ∆𝑡 (𝑨(𝒗𝑛 )𝒗𝑛 −1𝑫𝒗𝑛 ))
𝑅𝑒

(40)

− ∆𝑡𝑩 T 𝒑𝑛

ここで，左辺の𝜶の係数を見て，正射影行列𝑯 の性質から以下の関係が成立する．
𝑯 − 𝑩 T (𝑩 + )T 𝑯 = (𝑰𝑛 − 𝑩 T (𝑩 + )T )𝑯 = 𝑯𝑯 = 𝑯(41)

また，𝑯 は特異行列でもあり，逆行列を定義できない．ただし，ムーア・ペンローズ一般逆行列を用いると，𝑯 + = 𝑯 となる．したがって，任意の列ベクトル𝜶は下式で表せる．
𝜶=𝑯1
𝒗̃ + 𝑩 T (𝑩 + )T −𝒗𝑛 + ∆𝑡 (𝑨(𝒗𝑛 )𝒗𝑛 −𝑫𝒗𝑛 ))((𝑅𝑒− ∆𝑡𝑩 T 𝒑𝑛

Fig.4 cavity 流れ(Re=100)(42))

上式の右辺は，すべて既知量で表せる．𝜶を式(33)に代入することで，𝒗𝑛+1 を求めることができる．
また，定式の中では，移流項と拡散項の係数行列と速度は分けて表現しているが，格子毎に各項を計算してベクトルで表現
すると計算効率化が図れる． また，Euler 記述において解析中𝑩マトリクスは変化しない．毎ステップ収束計算の必要もなく，代入計算のみで求解できることは効率的であると考える．
4) 数値解析事例
これまでの定式化を基にcavity流れの数値解析を行った．非圧縮性非定常粘性流れのベンチマーク問題として扱われる解析事例である．レイノルズ数を 100 としたときの流れ場の様子をFig.4 に示す．定量的な検討は今後必要といえるが，流れ場の特徴をよく表しているといえる．

4．おわりに本報では，Lagrange 記述，Euler 記述の手法における流体問題において一般逆行列理論の有用性を示した．Lagrange 記
述においては，今後運動方程式を導入したモデルを検討した
い．Euler 記述においては，定量的な検討が必要といえる．
[参 考 文 献]1) 半谷裕彦, 川口健一, 「形態解析」, 培風館, 1991
2) 川口健一, 「一般逆行列と構造工学への応用」, コロナ社,
20113) 宮田真紀, 斎藤公男, 他, 「張力膜構造のポンディング現象に関する基礎的研究」, 日本建築学会関東支部研究報告集,
20024) 久田俊明, 野口裕久, 「非線形有限要素法の基礎と応用」,
丸善, 19965) 中山司, 「流れ解析のための有限要素法入門」, 東京大学
出版会, 2008

Application of generalized inverse to fluid analysis‐methods based on Lagrangian and Eulerian descriptions Yohei Yokosuka *1)
Teruo Matsuzawa*2)
SYNOPSIS
The feature of methods involving the shape-finding analysis by a generalized inverse is that a generalized inverse provides
meaningful values such as the least squares solution and the minimum norm solution with the analysis. In addition, it is capable by using acomplementary space wherein a generalized inverse provides the solution in which the right side of an equation is a zero vector. The
methods wherein the right side of an equation is a vector of strains, deformations, and external forces have been proposed in past studies.
In this paper, we propose methods of fluid analysis by using a generalized inverse. These are based on each expression of the Lagrangianand Eulerian descriptions. In the Lagrangian description, we assume a numerical analysis model that is a phenomenon of ponding on amembrane structure. In the Eulerian description, we propose the method of numerical analysis of an unsteady fluid field with anincompressible fluid.
*1） Assistant Prof. at Showa Women’s University, Dept. of Human Environment Science / Graduate Student, School of Information Science, JAIST*2）Prof., Research Center for Simulation Science, JAIST