曲率計算と体積補正

本研究の壁境界モデルで用いる壁カーネル関数は曲率 0 の平面を仮定しているため，壁 境界の近傍粒子に壁が与える影響の考慮には，近接境界面の有する曲率に見合った体積補 正が必要となる．原田ら ³⁴⁾の手法では，粒子近傍の境界が成す平均曲率を平面角と関連付 けて体積補正を行い，平均曲率は距離関数の空間テーブルより差分計算することで，解析 領域内の任意座標における値を計算することができる．本研究では距離関数は用いないた め，境界面の表面形状から頂点位置の平均曲率を計算する．なお，佐々木ら ³⁵⁾の手法では 壁境界の成す角度から直接的に体積補正係数を算出しているが詳細は不明である．

本研究では，平均曲率の計算手法として，Meyer ら ³⁶⁾による三角形パッチの曲率計算手 法を採用した．この方法では次式により頂点の平均曲率ÂÃを計算する．

ÂÃ m< =⁹₈‖Å m< ‖ (4.100) Å m< =_8Æ⁹

ÇkÈ–Z∑ dcot !c c+ cot •cedmc− m<e (4.101) ここで，Å m< は曲率の大きさを持つ頂点法線ベクトルであり，本研究では SPH 法の計算 に用いる頂点法線ベクトルとは分けて計算した．!cと•cは頂点’まわりの頂点“との間にそれ ぞれ定義される角度であり図－4.8に示す．»Ê<$•Eは着目頂点回りの構成三角形において頂 点が占める面積の総和である．

図－4.8 着目頂点回りの頂点間に定義される角度と頂点の支配面積の総和計算

»Ê<$•Eは，アスペクト比が良好な三角形の場合にVoronoi領域の面積となる．そうでない

場合（鈍角を有する三角形）については，着目頂点の成す角が鈍角である場合に三角形面

積の1/2，そうでない場合に1/4をもって代替することをMeyerらは提案しており，本研究 においても採用した．

なお，上式より計算される曲率は符号を持たない．本研究では曲率に符号を持たせるた めに頂点法線ベクトルを用いて次式に示す計算を行った．

·<= − ∑ d³c <∙ jce − < (4.102) ここで，³<は頂点’の単位法線ベクトルであり，jcは頂点’まわりの頂点“の座標であり，_< は 頂点’と原点間の距離である．·<の符号は曲率の符号を意味し，頂点’が凹面にある場合に正 値となり，凸面にある場合に負値となる．

頂点位置における曲率の計算時には頂点法線ベクトルも計算した．曲率計算は壁面数が 多い場合には計算時間が長くなるが，変形の速度が非常に大きい場合を除き通常大きな変 化はない．よって，本研究では，頂点法線ベクトルは毎ステップ更新したが，効率化のた め頂点曲率計算は境界の変形速度に応じた間隔で更新した．

前項で説明した壁境界の計算を行う際には，近傍粒子の投影点における補間された曲率 を用いて体積補正を行う．壁カーネル関数は曲率 0 の平面内に仮想粒子を配置して積分値 を計算しており，曲率を持つ壁面に対しては影響半径内の壁体積に誤差が生じている．そ こで，壁の曲率に応じた補正係数を壁の積分値に乗じることで生じる誤差の低減を図った．

平面上における重みを実際の曲面上における重みに補正する係数は，図－4.9に示される幾 何関係から計算した．図中のËはそれぞれ凹面と凸面における曲率半径であり，平均曲率ÂÃ

の逆数である．

図－4.9 凹面と凸面における体積補正係数の計算

すなわち，壁面内にその一部が存在している影響半径5•の球体積の内，壁面の内側に包含さ れる体積を基準体積として体積補正係数を計算した．凹面（曲率が正）の場合，半径5•の球 体積から壁面に接する曲率半径Ëの球体積に含まれない部分球体積（図中の赤い領域）を計 算し，基準体積で除算した．凸面（曲率が負）の場合は，半径5•の球体積と，壁面に接する

曲率半径Ëの球体積の共通部分の体積（図中の青い領域）を計算し，基準体積で除算した．

つまり，補正係数の関数テーブルは壁面距離と曲率の関数として定義した．具体的に 3 次 元の場合の補正係数の計算式を以下に示し，その計算例を図－4.10にグラフで示す．

最初に曲率0の場合の基準体積Ìは次式から計算できる．

Ì =^S_W 5•− 5<” 8 25•+ 5<” (4.103)

次に，積分半径5•の球と曲率半径Ëの球の交差円と壁面までの距離Íおよび積分球の端ま での距離‚を計算する．距離Íと‚は壁面が凹面の場合は次式から計算できる．

Í =_{8 Î)6}^6–0)6k±0_k± (4.104)

‚ = 5•+ 5<”− Í (4.105)

また，壁面が凸面の場合は次式から計算できる．

Í =_{8 Îy6}^6–0)6k±0_k± (4.106)

‚ = 5•− 5<”− Í (4.107)

計算した 2 つの距離を用いて図中の赤い部分もしくは青い部分の体積は次式から計算で きる．

Ì =^S_WÏ‚⁸ 35•− ‚ + Í⁸ 3Ë − Í Ð (4.108) ただし，凹面で2Ë ≤ 5•+ 5<”の場合，および凸面で2Ë ≤ 5•− 5<”の場合は次式で示す曲率半 径球の体積で置き換える．

Ì =^USÎ_W² (4.109)

以上の体積計算より，本研究では曲率による体積補正係数 _™ª«を次式で計算した．2 次元 解析の場合についても，Ìと Ìの計算式が異なること以外は3次元の場合と同様である．

™ª« =_Ñ^Ñ

u (4.110)

図－4.10 3次元解析における体積補正係数の計算例

0 1 2 3 4 5

-50 -40 -30 -20 -10 0 10 20 30 40 50

Volume correction factor

Curvature factor z/re=0.0

z/re=0.2 z/re=0.4 z/re=0.6 z/re=0.8 z/re=1.0