4.4.1 時間増分の設定
数値解析では問題に応じて適切な時間増分を設定する必要があり，SPH 法の数値解析に おいても他の解析手法と同様の時間増分に関する制約条件がある．SPH 法による流体解析 で考慮する必要のある数値安定性上の制約条件を表－4.1にまとめて示す．
表－4.1 流体解析における時間増分に関する制約条件
時間増分に関する制約事項 制約事項の回避条件
縦波の伝搬速度による制約 圧力計算を陰的に行うことで回避可能 Courant 数の制約 非定常計算である限り回避不可能
拡散数の制約 粘性項の計算を完全陰解法で行うことで回避可能
縦波の伝搬速度による制約条件は圧力を構成式から陽的に求めた場合に発生する．例え ば水の場合，縦波の速度は約1500m/secであり，時間増分としては非常に小さい値を設定す る必要がある．陽解法を用いた一般の計算においては構成式のパラメータを調整し，実際 より縦波の速度を小さく設定することで時間増分を粗くしている．
Courant 数Dの制約条件 28)は，連続体解析において縦波の伝播速度による制約と同義であ
るが，圧力を陰的に計算することで縦波の伝播速度の制約を回避した場合においても，実 際の粒子移動速度によってこの制約が発生する．流体解析の時間増分を一般に支配する条 件であり，この条件が満たされないと粒子の異常貫入や貫通が生じ，計算が不安定化する．
D =∆$∆iŸ^ $ < 1 (4.70)
飛翔する孤立粒子の速度は大部分の流体速度に比べて非常に高速であり，時間増分をクー ラン数の条件から可変計算した場合，最大速度Ÿ^ $が極一部の飛翔粒子によって決定する ことで，時間増分が小さくなり過ぎて計算の続行が困難となったり，圧力振動が抑制でき なくなり妥当な解が得られなくなったりする．本研究においては，時間増分を可変計算す る場合には最小時間増分を設定することでこの問題の回避を行っている．なお，粒子法で 安定した計算を行うにはCourant数は0.2以下が好ましい4)．
拡散数¡の制約条件28)は，少なくとも土木建築分野における解析ではあまり問題となるこ とがない．それは計算体系が大きいため，拡散数があまり大きくならないためであり，非 常に高い分解能の計算を行う場合にはCourant数より厳しい条件となる．また半固体状の極 めて粘度が高い流体の解析においても拡散数は大きくなり，必要に応じて粘性項は陰解法 で計算する方が賢明である．
¡ = λ∆$∆i0< 0.5 (4.71)
ここでλは拡散係数であり，流体解析においては動粘性係数が相当する．拡散数についても 0.2以下程度の基準を設けるのが好ましいが，本研究においては粘性項のラプラシアンの計
算にMorris らによる近似式を用いた場合，極座標のラプラシアンを用いるより数値発散が
生じる閾値が厳しいことを確認している．
4.4.2 人工粘性の導入
時間増分の制約条件に関わらず時間増分は一般に細かく設定した方が，非定常現象の再 現性は高くなる．特に粒子法による非圧縮性流体の解析では，粒子配置に局所的に生じた 粗密が原因となって圧力解が不安定となり，結果として適切な解析結果を得られなくなる ことも多く，このような問題に対して時間増分を細かくすることは有効である．しかしな がら，計算コストの観点からは時間増分を粗くとる必要があり，本研究ではMonaghan12)に よってSPH法に導入された人工粘性を用いた．この人工粘性は従来差分法で用いられてい た密度変化に対する人工粘性応力と比較して，高速衝突などを含む密度変化の激しい問題 に対して有効であり次式で与えられる．
s<c 6i=66–d£kl∙mkle
kl0y¤6–0 (4.72)
q<c 6i= s<c 6id!w + •s<c 6ie (4.73)
−9h∇q<c 6i = − ∑ h^l
khlq<c 6i∇ d5<c, ℎe
c (4.74)
ここで，s<c 6i
は速度の次元をもち，q<c 6i
は圧力の次元をもつ人工粘性圧力である．wは縦波 の速度であり !，•，•は計算パラメータである．!と•は大きいほど人工粘性は大きくなり，
•は反対に大きいほど人工粘性は小さくなる．•は粒子が異常接近した際の安定化パラメー タであるが，本研究では常に 0 を用いている．本研究では人工粘性圧力を通常の圧力勾配
項の一部として流体粒子に加えており，密度の2乗で除算した次元を持つ点でMonaghanが 導入した式と異なる．これはプログラム内で圧力勾配項や粘性項と同時に処理しているた めである．他の研究にみられる人工粘性の多くは，上式と同じ定式に基づくが密度の 1 乗 で除算したものであり，比較すると流体密度の 1 乗ほど次元が合っていない．例えば流体 が水である場合，他の研究例に示されるパラメータに対して人工粘性は1/1000となる．人 工粘性の大きさはパラメータで調整できるため，計算上の問題はない．
4.4.3 表面張力計算
土木構造物のような大きなスケールを対象とした解析の場合はあまり問題となることは ないが，実験規模の小さなスケールのモデルについては表面張力の影響は無視できるもの ではない．特に液体金属などは表面張力係数が大きく，表面張力の効果を解析に組み込ま なければ定性的にも現象を再現することはできない．本研究においても自由表面の安定化 を意図して表面張力モデルを導入した．
表面張力モデルはFVMなどの格子法においてはCSFモデルが採用されており，粒子法に も野村ら4)によって導入されている．しかし，CSFモデルの計算においては界面の法線と曲 率を評価する必要があり，粒子法の計算においてはこれらを精度良く計算することが容易 ではない．本研究においても当初はこのモデルの導入を行ったが，後により安定して計算 可能な分子間力モデルを導入した．分子間力モデルは近藤ら29)によってMPS法に導入され た粒子法に適した表面張力モデルであり，界面の法線および曲率を計算する必要がない他，
壁面との接触角を計算できることが特徴である．
近藤らは粒子間に働くポテンシャルとして次式のポテンシャル関数q 5 を提案した．
q 5 =YWJ5 −W85 +985•L 5 − 5• 8 (4.75)
ここで，Dは表面張力の大きさを調整するポテンシャル係数であり，5は初期粒子間距離で ある．この関数を5で微分すると次式のポテンシャル勾配が得られる．
"‘ 6
"6 = −D 5 − 5 5 − 5• (4.76)
近接粒子間ではポテンシャル勾配による作用力が相互に働いており，対称性のない部位で はポテンシャル力による移動が発生して全体のポテンシャルを最小化する．D = 1として正 規化距離で計算したポテンシャル関数とポテンシャル勾配を図－4.4に示す．
図－4.4 ポテンシャル関数
図よりポテンシャル勾配は影響半径の内側で負の値を示しており，これは接近する粒子に 対して斥力を生じることを示している．MPS 法においてはあまり問題とならないことであ るが，SPH 法では初期粒子間距離以内に接近する粒子に対して圧力勾配項があまり大きく ならないため，局所的に生じる粗密が解消せずに問題となることがあり，このような斥力 を生じるモデルは有効である．一方で，MPS 法では重み関数の特性から異常近接粒子には 過大な圧力勾配項が発生するため，MPS 法の商用コードでは近藤らの提案したポテンシャ ル関数と異なり，全域で引力が生じるモデルが採用されている．
ポテンシャル関数の大きさを調整するポテンシャル係数Dは，表面張力係数と対応付けて 設定する．近藤らは図－4.5に示すように平面上に接地面積58の1列の粒子を配置し，これ らの粒子群（’）を平面の粒子群（“）から引き抜くのに必要なエネルギーを計算し，次式で 表面張力係数¥と対応づけた．
2¥58= ∑ ∑ qd5< c <ce (4.77) -0.35
-0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ポテンシャルポテンシャルポテンシャルポテンシャル
正規化距離 正規化距離 正規化距離 正規化距離
ポテンシャル関数 ポテンシャル勾配
図－4.5 ポテンシャル係数の計算
また，この分子間力モデルでは壁面との接触角の考慮が可能である．表面張力係数¥と接 触角Bの関係式であるYoungの式とポテンシャル係数を関連付け，近藤らは次式を提案した．
Y¦§
Y¦¦=98 1 + cos B (4.78)
ここで，D#;は流体-固体間のポテンシャル係数であり，D##は流体間のポテンシャル係数で ある．つまり，ポテンシャル勾配の影響半径内積分を実施する際に，壁面が相手粒子であ るときのみ固体-流体間のポテンシャル係数を選択することで接触角の考慮が可能である．
4.4.4 人工斥力モデル
粒子法の計算においては圧力勾配項以外の斥力モデルの導入が必要である．最大の理由 は自由表面粒子には圧力がなく，そのままでは反発力が得られないためであり，次に計算 の安定化の都合，異常接近粒子を引き離す必要があるためである．
越塚らによるMPS法の解析コード4)には次式で示めされる衝突計算モデルが含まれてお り，MPS法の商用コードにも同じ機能が組み込まれている．
j<©ª«= j<~−hkhyhll 9y¬ dj6l{)jk{e∙£kl
kl0 (4.79)
ここで，j<©ª«
は衝突解消後の速度であり，-は跳ね返り係数である．この計算は非弾性衝突 のモデル式であり，粒子間距離が初期粒子間距離より小さい閾値未満となった場合に，ダ ッシュポットにより跳ね返り計算を行う．
また，ダッシュポットによる跳ね返りモデルに加えてばねによる斥力モデルも追加する ことが可能であり，これはそのまま DEM の作用力伝達系モデルと同じものと考えられる．
岩本ら30)はSPH法において生じる局所的な粗密を解消するためにSPH粒子間にDEMの作 用力伝達系を挿入し，粒子間距離が閾値以下の場合に有効化した．ばね定数は時間増分か ら設定可能であり本研究では，κを1.0以下の値として次式から計算している．
3<c©ª«=∆i¯0^^k^l
ky^l (4.80)
このばねによって生じる力は∆€後に貫入を解消するばね力であり，ダッシュポットの反発 力と合わせて計算する．ダッシュポットの減衰定数は臨界減衰を基準とする条件や跳ね返 り係数から設定する．岩本らは地盤の解析において人工斥力モデルの発動する距離の閾値 を初期粒子間距離の0.8倍程度と設定したが，本研究においてはこの距離がスロッシング問 題などの波高応答に影響を及ぼすことを確認しており，自由表面粒子間以外では0.65倍の 値を基準として設定した．なお，自由表面については安定化のため0.9倍程度の値を基準距 離に設定しており，この設定は粒子数密度の小さい自由表面粒子が表面で凝縮することに よる計算の不安定性を解消する．