09:45 〜 10:00
[HDS09-04] 非線形長波方程式の有限差分法による津波計算における計算不安定の原因とその対処
キーワード：津波数値計算、計算不安定、非線形長波方程式
津波数値計算で最もよく用いられる非線形長波方程式を有限差分法で解く際に、CFL条件を満たしているにも関わらず、計算が不安定となり、その結果が発散することがある。リアルタイムシミュレーションを用いた津波警報の発出、更新などを考えた時、そのような計算の発散による計算の失敗は致命的な状況となりうる。そこで、それらの計算不安定を無くすために、それらの差分式の評価と計算条件について検討し、どのような方法であれば、計算不安定になることを防ぐことが出来るか検証した。
まず、筆者はこれまでに非線形長波方程式の摩擦項の時間差分の方法による違いについて検証を行っており(MINAMI 2020JpGUなど)、その中で非線形長波方程式での摩擦項の時間差分について、現在使われている差分法である、(1)Explicit (2)Semi-implicit(combined) (3)Semi-implicit(simple) (4)Full-implicit の4つについて、その離散化誤差と計算安定性について検討し、その中でも、(3)Semi-implicit(simple) と(4)Full-implicitが計算安定性で優れていることを示した。しかしながら現在の津波モデルでは、計算安定性のために様々な条件をつけ、(2)Semi-implicit(combined)を用いて計算を行っているものが多い。そこで本稿では、津波数値計算で頻繁に利用、参照されるTUNAMI-N2モデル(Imamura et al. 2006など)での計算条件を用いて同様の計算安定性の評価を行った。
TUNAMI-N2モデルでは、摩擦項の差分にSemi-implicit(combined)を用いており、さらに以下の4条件を加えて計算している。1.流束計算前に該当メッシュの水深が0.1cm未満の場合、その時間ステップの流束を0にする 2.摩擦項計算時の全水深が1cm未満の場合、それを1cmに切り上げする 3.移流項を計算する際に、全水深が0.1cm未満の場合、移流項を計算しない(0.1cm以上の時も、風上側の全水深が0.1cm未満の場合、風上の項のみ計算しない) 4.計算後の流束が10^-10mより小さい場合、その流束を0にする これらの条件を付し、様々なデータを用いて計算を行った結果、上記の条件があった場合でも、Fig.1のようにSemi-implicit(combined)では計算に大きな誤差が生まれる事例があった。これは、連続的に解いた場合には流束が0近くに収束する事例であり、この時に、摩擦項の差分化方法の違いにより、流束が0に収束せず、プラスとマイナスが逆転してしまう場合があった。この流束方向の逆転が起こると、移流項を計算する際に風上差分の風上方向が入れ替わってしまい、本来参照すべき風上の値が入れ替わってしまう。それが時間発展し、誤差が拡大することが原因であることが分かった。またそのような場合でも、(3)Semi-implicit(simple) もしくは(4)Full-implicitであれば、もともとの流束計算が正しく収束するため、風上の方向は入れ替わらず、計算が安定することが分かった。これらを上記の条件の数値を変えることによって回避することは可能であるが、それは即ち水深が小さい時に移流項と摩擦項がなくなることであり、それによっても誤差が大きくなってしまう。つまりこれらの計算安定性を回避するには、やはり(3)Semi-implicit(simple) もしくは(4)Full-implicitを使うことが適切と考えられる。
さらに計算領域の端に当たる部分での計算が不安定になる場合がある点についても検証した。津波計算において、有限差分法で放射境界を利用する場合、進行性長波の特性曲線を用いた方法(1982 後藤・小川)が採用されることが多い。この方法は線形長波を仮定しており、計算領域の端に水深の浅いメッシュがあった場合、摩擦項、移流項の影響が大きくなり、誤差が大きくなる。そこで放射境界での計算においても摩擦力を考慮して計算したところ誤差が小さくなることが分かった(Fig.2)。なお、実験の水深は一様で、周期境界条件を用いた計算との比較を行っている。この時、摩擦項をそのまま適用した場合は誤差が大きく、通常の1/2の摩擦項を適用すると、周期境界条件で計算した場合との誤差が少なかった。
以上の検証から、非線形長波方程式を有限差分法で解く際に、摩擦項の差分式をより陰的なものにする、計算領域の端に当たる部分の処理をより適切に行うことで、これまでよりも安定して計算することが出来ることが確認出来た。
津波数値計算の安定性は実務では非常に重要な要素であり、これからもより精緻にその条件を検証する必要がある。
まず、筆者はこれまでに非線形長波方程式の摩擦項の時間差分の方法による違いについて検証を行っており(MINAMI 2020JpGUなど)、その中で非線形長波方程式での摩擦項の時間差分について、現在使われている差分法である、(1)Explicit (2)Semi-implicit(combined) (3)Semi-implicit(simple) (4)Full-implicit の4つについて、その離散化誤差と計算安定性について検討し、その中でも、(3)Semi-implicit(simple) と(4)Full-implicitが計算安定性で優れていることを示した。しかしながら現在の津波モデルでは、計算安定性のために様々な条件をつけ、(2)Semi-implicit(combined)を用いて計算を行っているものが多い。そこで本稿では、津波数値計算で頻繁に利用、参照されるTUNAMI-N2モデル(Imamura et al. 2006など)での計算条件を用いて同様の計算安定性の評価を行った。
TUNAMI-N2モデルでは、摩擦項の差分にSemi-implicit(combined)を用いており、さらに以下の4条件を加えて計算している。1.流束計算前に該当メッシュの水深が0.1cm未満の場合、その時間ステップの流束を0にする 2.摩擦項計算時の全水深が1cm未満の場合、それを1cmに切り上げする 3.移流項を計算する際に、全水深が0.1cm未満の場合、移流項を計算しない(0.1cm以上の時も、風上側の全水深が0.1cm未満の場合、風上の項のみ計算しない) 4.計算後の流束が10^-10mより小さい場合、その流束を0にする これらの条件を付し、様々なデータを用いて計算を行った結果、上記の条件があった場合でも、Fig.1のようにSemi-implicit(combined)では計算に大きな誤差が生まれる事例があった。これは、連続的に解いた場合には流束が0近くに収束する事例であり、この時に、摩擦項の差分化方法の違いにより、流束が0に収束せず、プラスとマイナスが逆転してしまう場合があった。この流束方向の逆転が起こると、移流項を計算する際に風上差分の風上方向が入れ替わってしまい、本来参照すべき風上の値が入れ替わってしまう。それが時間発展し、誤差が拡大することが原因であることが分かった。またそのような場合でも、(3)Semi-implicit(simple) もしくは(4)Full-implicitであれば、もともとの流束計算が正しく収束するため、風上の方向は入れ替わらず、計算が安定することが分かった。これらを上記の条件の数値を変えることによって回避することは可能であるが、それは即ち水深が小さい時に移流項と摩擦項がなくなることであり、それによっても誤差が大きくなってしまう。つまりこれらの計算安定性を回避するには、やはり(3)Semi-implicit(simple) もしくは(4)Full-implicitを使うことが適切と考えられる。
さらに計算領域の端に当たる部分での計算が不安定になる場合がある点についても検証した。津波計算において、有限差分法で放射境界を利用する場合、進行性長波の特性曲線を用いた方法(1982 後藤・小川)が採用されることが多い。この方法は線形長波を仮定しており、計算領域の端に水深の浅いメッシュがあった場合、摩擦項、移流項の影響が大きくなり、誤差が大きくなる。そこで放射境界での計算においても摩擦力を考慮して計算したところ誤差が小さくなることが分かった(Fig.2)。なお、実験の水深は一様で、周期境界条件を用いた計算との比較を行っている。この時、摩擦項をそのまま適用した場合は誤差が大きく、通常の1/2の摩擦項を適用すると、周期境界条件で計算した場合との誤差が少なかった。
以上の検証から、非線形長波方程式を有限差分法で解く際に、摩擦項の差分式をより陰的なものにする、計算領域の端に当たる部分の処理をより適切に行うことで、これまでよりも安定して計算することが出来ることが確認出来た。
津波数値計算の安定性は実務では非常に重要な要素であり、これからもより精緻にその条件を検証する必要がある。