Feel＆Think

今回からは、地下水流れや熱伝導のような場の問題に陽解法を適用してみます。動的陽解法と同様に、連立方程式を解かず時間刻みでの繰り返し計算によって解を求める方法を説明します。ポイントは、動的陽解法でも出てきた対角化により節点ごとでの除算を可能とすることと、既知の値だけから必要な解を得る時間差分です。

今回は地下水流れの問題です。最初に、地下水流れの支配方程式と離散化手法を導いておきます。長くなりますので、詳細は「地下水流れに関する支配方程式と陽解法」に示すとして、ここでは要点を述べます。 まず、間隙水（地下水）の質量保存則から次式が得られます。

ここに、φは間隙水圧、n、S_wはそれぞれ地盤の空隙率と飽和度、vは間隙水の平均流速です。なお、添え字の i、j は座標系を表します。

–∂S_w/∂φは比水分容量、-∂n/∂φは比貯留係数と呼ばれる物性値です。不飽和領域では間隙水圧変化が空隙率に影響を及ぼさないと考えられることから、-∂n/∂φには飽和領域で1、不飽和領域で0となる係数αを乗じています。また、K_bは水の入出量差によって生じた体積変化に対する間隙水圧変化を定めることから、以後は空隙弾性定数と呼ぶこととします。

他方、間隙水の平均流速に関しては、Darcyの法則が成り立つものとします。

ここに、ρ_w、gはそれぞれ間隙水の密度と重力加速度、kは透水係数です。

陽解法とするために、まず式(3)にガラーキン法（内挿関数を重みとした重み付き残差法）を適用し離散化します。これが一つ目のポイントです。このとき、一次の内挿関数を用い、要素重心で間隙水圧を定義し、流速は節点で求められることとします。一次の内挿関数を用いた変形問題では、要素重心で応力が定義され荷重が節点で求められますが、これと同じ方法です。この方法では要素内で応力と間隙水圧が同じ点で定義され、後に紹介する変形と地下水流れの連成解析などで、物理的な意味がわかりやすくなるメリットがあります。

要素単位での合算が終わった状態では、節点値について次のような式が得られます。

ここに、{F_e^t}は各節点における等価水圧勾配（重みの付いた水圧勾配）、{F_s^t}は境界面の水圧と等価な節点水圧勾配、Wは重みです。また、tは現在の時刻を表しており、後に示す時間差分では既知量として取り扱うことを示しています。

動的陽解法でも出てきましたが、左辺第1項の係数[W]は対角マトリクスとしています。これが二つ目のポイントです。これにより、節点ごとの割り算により節点流速を求めることができます。

式(1)の方程式に関しては、要素内の流速の勾配は内挿関数を用いて節点値から求め、時間微分は差分で近似することとします。

[N]は内挿関数からなるマトリクス、[∂]は空間に関する偏導関数マトリクスです。あるいは、勾配マトリクス[B]を用いると次式となります。

式(5)と式(7)を組み合わせることで、既知の値から次の時間の間隙水圧を求めることができ、陽解法により地下水流れの解析が可能となります。

すなわち、節点で求められた流速より要素内の流速勾配を求め、これに空隙弾性定数を乗ずることで間隙水圧の変化を算定します。

連立方程式を解かない陽解法では、少ないメモリで計算が可能である反面、安定な解を得るための時間刻みΔtに制約があります。筆者は、隣の要素とΔφの間隙水圧の差が生じている地下水流れにより、Δtの間にΔφ以上の間隙水圧変動が起こらないことが安定の条件と考えています。このようなことが起これば、Δtの間に地下水流れの向きが反転してしまい、解の振動が発生します。式(8)を参考にすれば、この条件は次式となります。

ここに、L>sub>minは解析対象を構成する要素の中で最も小さな節点間距離です。したがって、Δtには次の制約が課せられます。

筆者の経験では、この値の半分以下で安定した解が得られるようです。

解析例を一つ紹介します。この例では、図-1に示したモデル底面から45mの位置に初期地下水位を設定した後、右側境界の水位を100mまで上昇させています。解析に用いた物性値や不飽和特性を表-1と図-2に示します。

解析結果は、以下の動画でご覧ください。モデル右側から地下水の流入があり、新しい地下水面が形成されていく様子がわかります。陽解法でもきちんと解けていることが理解いただけると思います。次回は、物質移行の問題を陽解法で解いてみます。