国立防災科学披術七ソター研究疑告第18号

1977年11月

556．14：556，332．342

不飽和浸透を考慮した降雨浸透の有限要素法による一解法
定

常

1充

倉

博＊

大

国立防災科学技術セ1／ター

An

Ana1ysis

Methodoffmm

Saturated−Unsat㎜rated

Rainfa11by

FiniteSteadyF1ow

E1ements

By

HiroshiNαOhkura・舳げ㈹ακ乃C・〃θけ・・眺α・伽〃舳・〃・・，1ψ伽AbstractAissaturatcd−uI］saturated

formu1ated

e1lipticasaboundary

type．

On

the

boundaries

which

steadyva1uerainingcann〇七beHowofprob1emofsurface

arr乱ngedandsoilwateral〕（1groun，vater

non1inear

seepageboundary

frompartia1diHerentialout

surface，there

conditionsarainfal1equationareofannoying

priori．Overcoming

the…1i…1i亡y・・dth・mi・・db・md・・y…diti…，th・b・㎜d・・y・・1・・p・・bl・mi・transformed
a11d

an

into

non1incar

alogorithmfor

sin1ultaneous

solvingtheequationsnon1inearbyGalcrkin丘nite

simuhaneous

equationse1ementismcthod，

proposcd．

By

th・・1g・・ithm・…士it・七・d・■t…ti・・m・th・d，・…id…b1y・t…g…1i・…p・・b1・m・canbc

simu1atesolved．

This

nonhncarsuggests

phenomem

that

which

there

can

are

be

nlany

applicationstransformed

intoof土hismn1inearalgorithm

tosill／ultaneous

equa七iOnS．

1．希者

堤防やロックフィルダムなどの土構造物内の浸透流の解析や地下水の解析に際し，不飽和浸透を考慮すると複雑になるため，飽和浸透だけを考慮した報告が多い（駒田・金沢，1975：

Kono，1974；Neuman，1971；山上・小田1975）．
不飽和浸透では毛管現象により負圧ヵミ生じ，この負圧が土構造物内の浸透流や地下水の流れに影響を与えるため，不飽和浸透を無視し得ない場合が生じる．
＊第4研究部計測研究室一51

国立防災科学技術セソター研究報告第18号

1977年11月

自然斜面の崩壊は斜面土中の含水率の分布に密接に関係し，斜面を構成する土の大部分が不飽和でも崩壊が発生することがある．このため，自然斜面の崩壊機構の解明を念頭に置き
含水率の分布を浸透流の解析より求めるには，不飽和浸透を考慮したけれぱならない．

植物の根に接する土壌は不飽和の場合が多く，農地の灌激と排水の問題も不飽和浸透を無視できない．

不飽和浸透が雨水の地盤への浸透の大部分を担う場合があり，水文学の流出解析における流出率の決定にも不飽和浸透の解析は有為な情報を与える．

不飽和浸透を考慮すると，透水係数が飽和度により異なるため，流れを支配する運動方程式は非線形偏微分方程式になる．さらに，自由水面が存在する不圧流れでは湧出点などの位置が事前に定まらず，事前に境界条件を既定できない窺界が存在する．このため，流れを支配する方程式の解析解を得ることが困難になる．解析的手法に代わる非線形偏徴分方程式を基礎にする数値近似計算法として差分法と有隈要素法が挙げられる．有限要素法は，差分法に比べ境界条件の処理が容易なため，境界条件の変更が生じるときに有利になる．さらに，材料の異方性，部分的不均一性も容易に扱える．このため，有限要素法は，降雨による土中の飽和・不飽和浸透流への有力な解析手法となる．

本報告で述べる解析法は，降雨による土中水と地下水の不飽和・飽和定常浸透流をラプラ

ス型非線形偏微分方程式の寛界値問題に定式化し，これを有隈要素法を用いて非線形連立方
程式に離散化した後，得られた非線形連立方程式を解いている．

境界値問題への定式化において，地下水・土中水の湧出する境界の範囲と雨水の浸透する境界の境界条件が事前に定まらたい．すなわち，ラプラス形偏微分方程式の境界条件が第1種の境界条件か，第2種の境界条件かが既定されない．この障害は非線形連立方程式を反復法で解く過程で処理されている．非線形連立方程式を解くにあたり，非線形性が強く一般の反復解法（赤坂（1967））で解が得られたい場合でも，新しく提案された反復解法を用いて収束解を得た．

2．

浸透のモデル

図1に示す定常浸透問題を例にとる．抑はデカルト座標をあらわし，κ。を鉛直方向座標とする（以下，式の表現の簡便化を計るため，特にことわらないかぎり，右脚添字を添えた

量は，脚添字が1個のものはベクトル，2個のものはテソソルをあらわす．さらに，式中の同一項に同じ添字をもつ量が繰返し現われるときは，その添字について総和を取るものとする）．9は土中の流れの領域（不飽和域も含む）をあらわす．「と■は浸透流域の境界をあらわし，1が添字されると全水頭んが既定されている境界（楕円型偏微分方程式の第1種
境：界条件）を，2が添字されると浸透水の領域内への流入量が既定されている境界（第2種境界条件）をあらわす．「は事前に「全水頭が既定されている」か「流入量が既定されて

一52一

不飽和浸透を考慮した降雨浸透の有限要素法による一解法一大倉

いるか」が定まっている

境界をあらわし，■は事前にそれらがX3

定められない境界をあらわす．

〈2

H，V，Rはそれぞれ既定全水頭，単位時問・単位面積あたりHの既定流入量，降雨強度をあら

h

n

わす．

兄
〈1

全水頭乃を未知数とし，んを

■↑V

る．

1）水は非圧縮性とする．

解析領域の断面図

図12）全水頭は，ψを圧力水頭

Fig．1

〈、

X1

0000000000

求めるために次の仮定を設け

Cross

Sec七ionof

Analysis

Regionとして
ん＝／z（抑）：ψ十北。

（1）であらわされる．ψは，飽和領域でψ＞0，白由水面上でψ＝0，不飽和域ではψ＜Oとたる．

3）浸透流は，飽和・不飽和域を間わず，ダルシー則肋吻＝一瓦・1、、（2）に従う．ここで，伽は拘方向のダルシー流速を，K幻は透水係数をあらわす．透水係
数K幻は飽和浸透に対する透水係数珊と比透水係数Kτの積であらわされる．

〜＝珊Kτ（S）（3）

ここで，K『は飽和度∫の一価関数とする．
4）

飽和度Sは，日三力水頭ψの一価非減少関数で∫＝∫（ψ）＝S（トκ。）（4）ただし，ψ≧0のときS＝1とする．
5）境界面からの蒸発は考えない．

6）ム，ム境界で浸透不能な雨水と■ユ境界より湧出した水は瞬時に流れ去る。

7）ム境界では全水頭んがそのκ。座標に等しく，水の単位時問・単位面積あたりの流入量は，単位面積に降る降雨強度を越えない．

8）■、境界では単位時問・単位面積あたりの流入水量が単位面積あたりの降雨強度に等
Lく，全水頭んはその点の。座標値を越えない．

一53一

国立防災科学技術セソター研究報告第18号

1977年11月仮定7），8）について説明を加える．事前に「全水頭が既定されている」か「流入量が既定されている」かが定まらない■境界では，仮定6）より湛水が生じないので流入水量は単位
面積あたりの降雨強度以下になり，さらに，圧力水頭は負または零，すなわち全水頭ゐはその狛座標値を越えない．■境界におげる全水頭と流入水量との状態の組合わせを考えると，
仮定7），8）の他に，次の組合せが考えられる．「単位時間・単位面積あたり流入量は降雨強度より少なく，全水頭はκ。座標より小さい．」ところが，流入量が降雨強度よりも少ないときは■境界に浸透できない雨水が存在（瞬時に流れ去る）するから境＝界では飽和度∫＝1となり・仮定2）・4）よりん＝π・となり矛盾する．故に，■境界では仮定7），8）以外の状態の組合わせは生じない．

3．流れを支配する方程式と境界条件
定常流を考え・水の非圧縮性とダルシー則を考慮すると，連続の式から次式が求まる．

∂1。（榊（∫）昔）一・（・）

境界条件は，．境1界r。において

ん＝H・

（6）

境界「。において∂ん桝『（∫）可＝V・（・）

ここで，物は境界の外向き単位法線ベクトルのゴ方向成分をあらわす，
境界■。において，

ん＝κ・・

（8）

境＝界■。において∂ん
物K高κ『（∫）r一一＝〃3R・

（9）

0〜
となる．

4・非線形偏微分方程式の離散化
有限要素法は・変分原理に基づき流れを支配する徴分方程式と等価な汎関数を導き，この汎関数の極小化を基礎とする・しかし・非線形徴分方程式は等価な汎関数が存在しない場合

や，存在が不明な場合が多い．かかる場合にガレルキソ法を用いると汎函数を求めなくとも，有限要素法による非線形境界値問題の解析が可能になる．
（5）式の未知数んを変分法の比較関数とし，んを未知廉数ゐ帆と基底関数ψ、（κ、）の有限

級数

一54一

不飽和浸透を考慮した降雨浸透の有限要素法による一解法一大倉
伽。）＝ん曲（γ。）（10）で近似する．ここで，ρル｛）は以下の手続で求1

める．流域を図2に示すように3節点1次A

ヱ。／
n

要素（藤野（1972））に分割し各々に通し番号をつける．以下，、、方向に領域と境界条件が一定とし，γ、，i冗e。平面の垂直2次元間題を考える．番号〃の節点（以下節点〃とする）を頂m、点とする三角要素よりなる多角形（図2の斜線部分）の領域をク，苗とする．9肌（五乞）はπ、の連
Fig．2

図2要素分割Ne士workof

Triangular

E1emen七s続関数で節点〃上で1，多角形ρ肌の辺上と外側で0と定義する．このとき，節点〃上でん＝ん冊となり，ん柵が節点〃における全水頭をあらわす．

ρ。を9皿に含まれる要素番号θの要素（以下，要素θと呼ぶ）の領域とすると
9π＝U9。

（11）

9π＝U％。

（12）

％も同様に

であらわす．ここで，〜は9。内で一価連続，9色の外で零となる関数である．g肌、の具体的な関数形は付録2で示す．
（10）式ヵミ（5）〜（9）式の境界値問題を最良に近似する係数ん椛を決定するために，ガレルキ

ソ法を用いると∬、／、麦（榊（∫）㌶）一小一・（・・）となり，（10）式を代入して整理すると〜∬ム／、ζ（榊∫）峠）／ψ一・（・・）となる．∫。で要素θの重心の飽和度をあらわし，要素θ内でKτ（∫）を・一定

Kε㌧K『（∫侶）（15）とし，グリーソの第1定理を用いると（14）式は

早／い概肌努〜…い珊肌訂舳
一∬ム榊々い／一・
となる．

（・・）

ここで，中括弧内の第1項と第2項はそれぞれ節点πが「、，■、境界上にあるとき

一55一

国立防災科学技術セソター研究報告第18号

1977年11月

現われ，横分範囲几。とんはそれぞれ要素θの節点〃を含み「。，ム境界と一一致する辺である．％。の定義（付録2）により，（16）式を整理するとん。ん肌＝ρ。「・十ρ柵ム（17）㍍一那κ∬ム（訟繁一）〃1＝〜4∠κτ｛榊1・6㎜・十醐わ1l・㎜・・わ㎜1・伽）十蛉…榊｝

炸叱・〜〃

（18）

（19）

伽・一〜／ル・榊〃（20）

ここで，」，δ，、。，c，、。，δ、，。，c㎜。と伽は付録2を参照．「ユ，ノi。境界上にある既知水頭にダッシュをつけると，境界値問題（5）〜（9）式は，未知水頭ん、、の連立方程式

ん㎜ん肌：Bπ（21）

B冊＝Q、、＾十ρ，、！」ん〃（22）に離散化される．ここで添字1〃，〃は木知節点水頭だけを坂ることに注意されたい．

5．非線形連立方程式の解法
K『が∫の函数，∫がん一γ。の函数となるから，（21）式は，係数！、、、、がん、、の函数となり，非線型連立方程式になる．節点〃がノ11境界上にあると添字〃に対する（21）式が消失
し，λ，境界上にあると（22）式右辺第2項が現われる．また，節点一々がノi。境界にあると右辺第3項が現われる．よって，んムの範囲が（21）式の解んπに大きな影響を与える．故に，（21）式を解くには■i、，ノi。の範囲が楽前に定まらねばならない．しかし，2章7），8）より■1，ムの範胴を定めるにはん，すなわち（21）式の解が必要となり矛盾する．この矛盾を克服するため，始めにん冊、とノ1。，ムの範囲を仮定し，線形計算の繰返しにょり逐次ん㎜，ノ1。，

ムの範囲の修1H三を行ない（21）式の解を得る．

以下解法のアルゴリズムを述べるが，特に断わらないかぎり順を追って実行するものとする（図3参照）．

操作11一㎜のf直とノ1、，ノ1，の範閉の第0次近似ん肌o，」・o・ノ1・oを仮定する（以下，右肩添

字は反復回数をあらわす）．

操作2ゴ、を許容反復回数として
5：≦ら

（23）

を判定する．成立しないときは，収束解が得られないものとして計算を打切る．
操作3ん肌｛■1，ノi、｛一1，ノ1，｛■1を用いて（18），（19），（20），（22）式より（21）式の係数ん肌と定数項3柵を求め，それぞれ刈｛1，B凧｛■1とする．

一56一

不飽和浸透を考慮した降雨浸透の有限要素法による一解法一大倉START
ASSUME｝，パ，〃｛：：0

jl＝主十1

〜0

三≦三。

YE5
00舳τEハ幻8ケ1F尺0〃

1ノ；一〃一1

50L、！ε5〃、｛リしT＾川≡＝0〃S

し一，

1

EQリハ丁10〈

5ア0ρ
〈

δ

050Lリ丁0〜

ヱー1

ハ。加乃〃＝8η
JUD6E：一1わ守1／／1γεs〜0Jω6ε20
1≦ら

γE8

〈10

沽呂乃ξ1＋・（（μ一ゐξ1）

CON

PL1≡X

MOl）I

F

CAT

0〜

C0岬UTεノ1㍍〃

α、㌧ノヂ

〃0

）・（ノ三一ノξ）〈（JU06E＝7）
γES

5丁0P
00〜レER6EN
図3
Fig．3

T

60LUT10〜

解法アルゴリズムの流れ図
F1ow

Char七〇f

A1gorithm

to

So1ve

操作4 連立1次方程式
λ㌫1ん肌｛＝B椛｛一1

を解く．

操作5

εを許容誤差として，全ての刎に対して

一57一

（24）

国立防災科学技術セソター研究報告第18号

1977年11月

μ一んノ111≦ε（炸・）1 ㍑
1んLん肌｛一1j≦ε

（25）

（〃＝0）

（26）

を判定する，全ての舳に対して成立するときは，操作9へ飛越する．
操作6

ゴ、を許容単純反復回数（ゴ、≦ゴ召）として

ゴ≦ゴ岳

（27）

を判定する．成立しないときは操作8へ飛越する．

操作7αを適当な係数として，すべての舳に対して〃：＝〃

1＋α（㍑一〃

）

（28）

START
肋：＝1

〃0
加≦η4

γE8
切乃、1・片←1＋8加・4・・伽（乃ノー乃㌃I）

C0岬リ丁ε朋〃，曲F尺0州沽，ノ1
50LVε
舳・附1

5I〜uLTA〃E00S，ル1

ε0∪Aτ10W5

ηA〃 沁：切β。
○ノパ＝（切庇一んノ）／4−8加

D／1・カノーL〜

50しレε51MuLτANE0〃S

f

εQμAτ

0＾

5

C伽加・0ノ沽二力4一，十β〃閉
RETσ舳

ηど

図4

COmp1ex

Fig．4

F1ow

〜u〃βε尺0戸〃00εMOdi丘ca七ionの流れ図（解法アルゴリズムの流れ図の一部）

Chart

Algorithm

丁0τハL

to

of

Comp1ex

So1ve

Modi丘cation．（A

）

一58一

part

of

F1ow

Chart

of

不飽和浸透を考慮した降雨浸透の有限要素法による一解法一大倉を計算した後，操作9へ飛越する．以下，「：＝」は右辺を計算した後に左辺に代入することを意味する．

操作8a）ベクトル肌ゐの第〃成分肌伽を次式より求める（以下，操作8は図4参照）㎜々椛：＝乃旭｛一1＋δ肌肌・6・sign（ん、、乞一々皿ガ1）（29）

ここで，δ㎜冊はクロネッカーのデルタ，∂は定数，sign（Z）はZの符号をあらわす．
b）、、乃、、と■。｛一1，〃■を用いて，操作3と同様に（21）式のん1とB肌を求め肌んユ，

肌B冊とする．

c）肌伽を未知数とする連立1次方程式

帆ん＾肌伽＝肌B冊（30）を解く．ここで舳について総和は坂らない．

d）C三㎜とD1を次式より求める．
C1、、：＝（肌伽一伽｛）μ一δ〜，皿（31）

Dl：＝〃一」〃（32）

C1肌ル＝Dヱ（33）e）ルについての連立1次方程式を解く．

f）β｛■ユを適当なベクトルとして

ゐ肌包：＝〃11＋β壱．ソ帆（34）を計算する．

操作g乃肌也を用いて■1也，〃を求める（図3参照）．
操作10■1｛＝〃一

かつ■、｛＝∠、｛一1が成立し，かつ操作5が満されているならば，ゐ肌｛を連立方程式の収束解として計算を終了する．満されないときは操作2へ飛越し次回の反復を
繰返す．

以上がアルゴリズムの概要である．以下，説明の不十分な操作に説明を加える．
操作4と操作8の連立1次方程式，（24）式と（30）式はラプラス型偏徴分方程式（5）に対応するので，水頭既定の節点（「。，■。境界上の節点）が存在しないと不定になる．「。境界が存在しないとき，収束解に■。境界の存在が予想されるにもかかわらず反復計算中に■。｛一1境界が消減すると，（24）式と（30）式は不定となり，上述のアルゴリズムでは解が得られなくなる．この対策を述べる．操作3の実行にあたり「。，■1の存在を確認する．存在しない

ならぼ，事前に指定した1個の節点，たとえぼ節点冶がr。上にあるものとして操作3を実行する．ただし，ん㌃1，凪｛一1も計算する．操作4の実行の後，〃を次式より求める．

乃走ε：＝（即一1一λ㌫1乃肌）μ友1

一59一

（35）

国立防災科学技術セソター研究報告第18号

1977年11月（ここで，幻こついての総和は取らない，また，右辺分子第2項の舳の総和において刎＝冶
項は累加しない．）操作8に対しては節点后が「。境界上にあるものとして通常の計算を行
なう．

以上の対策を講じても，解が振動しながら収束する過程で比透水係数K、『がK、γ＝0となるとき，（24）式と（30）式が不定になることがある．K、『＝0の要素が連なって不透水層が生じ，流域が分離された形に（計算上）なってしまう．このとき分離された流域に「1または■エが存在しないと（24）式と（30）式は不定になる・この問題は，前もって，結果として予想される最少のK！よりも小さい楓i・（＞0）を定め，操作2または操作8のb）でK、τ＜楓i、ならぼ瓦『＝楓i。としてλが，3冊｛■1または肌んユ，刎B冊を求めて計算を続け，得られた収束解より全要素についてK！＞K二i、の成立を確認すれぼ回避できる．

操作7の修正に比べ，操作8の修正は非線型性のかなり強い連立方程式も確実に解けるが・1回の修正反復に対して連立1次方程式をその元数に等しい回数だけ解かなげれぼならず・演算時問が非常に長くなる・演算時問短縮のために，操作6が単純な操作7の修正をゴ、回反復しても収束解が得られないときに初めて，操作8の修正を行なうようにしている．

操作7の修正をこのまま行なうと，非線形性が強い場合に，αが大きく（α＞O．2）さらに乃ノの誤差が大きいと解が振動または発散し，操作8の初期値の誤差が大きくなる．このため演算時問の長い操作8の反復回数が多くなってしまう．しかし，αが小さい（α＜0．2）と非線形性が弱い場合でも操作7による解の収束性が悪くなる．そこで，次の対策を講じた．
すなわち，αを大きな値（α＞0．5）にし，非線形性の大きな（iゴK『μψ1の大きな）範囲内にψの値を適当な問隔で5・6個坂り，これを上昇順に並べ，づ番目の値をψ毒・とする．おのおのの節点に対し・操作7の結果・圧力水頭がψ！を1個以上飛び越えて変化するときは，ψ・を最初に飛び越されるψ！，κ肌ヨを節点舳のκ、座標値として

舳一κ伽・十ψ（36）と修正するわげである．

操作8が新しく提案された修正法である・理論は付録3を参照されたい．操作8e）の（33）式が付録3の（A皿・7）式に対応する・演算時問短縮のため，操作8では上述のa）〜f）の操作の他に次の操作を行った．「9ゴに含まれる要素の重心の飽和度∫色をそれぞれψ工とん、Hよ

り求め比較する・もし・巧内の全ての∫εがそれぞれ等しいならば，刎＝ノなる舳に対してb）・c）・d）を省略し・e）の連立1次方程式をルを除外して解き，f）で，閉＝プのとき

ゐ肌｛：一〃（37）とする．」このような省略が可能なのは次の理由による．巧に含まれない要素の∫、はψ互から求めても〃一1から求めても等しくなるから，巧に含まれる全ての要素の∫、がそれぞ
れこの2つの求める方に対して等しくなるならぱ，流域の全要素の∫、も等しくなる．この

一60一

不飽和浸透を考慮した降雨浸透の有限要素法による一解法一大倉

とき，要素の重心の比透水係数K、τがS、の1価関数であるから・全ての〃・1に対して（30）式と（24）式で！肌工一俳1，ゴB冊一B肌｛11（38）となり，c）のゴ石ヱと操作4の〃が等しくなる．よって，操作8e）の連立1次方程式の係数行列のノ列はC〃（＝一1）を除いて零となる．これはf）において・〃をゐ！一1の近傍で任意に修正したときに，次回の反復，操作4で得られる舳キノなる〃十1の値が変化しないことを意味する．一般に，これは，9ゴ内の全要素が自由水面下にあるときに一起きる．なぜたらぼ，自由水面下にある節点は，自由水面付近を除いて要素の和飽度は白由水面付近を除いて乃肌｛をん肌｛11一γ舳》0となるから，9ゴ内の1の近傍で任意に修正しても，S。＝1のままである．このとき，次回の操作4で（23）式の係数行列のノ行とプ列の要素Aふとλみ
さらに定数ベクトルのノ要素B！は乃肌也の修正に依存せず，（37）式の修正は次回の反復操
作4の解ん肌也十1に影響しない．

この操作が行われると，操作8c）の連立1次方程式を解く回数が少なくなり，さらに，e）の連立1次方程式の元数が小さくなるため，計算機演算時問を短縮出来る，特に浸透流域の大半が飽和している場合に，短縮が著しい．
f）における洋■1は

ザー讐

店、。
β｛一1＝1（鳥＞H、）（・・）（H、＞ツ篇㌧＞0）（41）より求めた．ここで，Hαは適当な定数，ツ鳳はツ焦為＝1刎olκ（ツ㎜）1とする．また，ツ肌α：Oのときは，（34）式をん肌｛：＝ん肌乞

（42）

とした．

操作gでは2章の7），8）の条件∂ん

境界■・上で境界■。上でん≦κ。

也蝋1、、≦

・R

（43）
（44）

を判定し，成立しない境界の範囲は，それぞれ■。境界を■。境界に，■。境界を■。境界に変更している．

6．数値計算例
図5に示す斜面の定常浸透を，境界ABを除き境界条件と物理定数を一定にし，ABの境界条件を自由排水（計算例1）または不透水層（計算例2）に変えた2例について解析した．

一61一

国立防災科学技術セソター研究報告第18号

1977年11月↓↓↓↓↓↓

500

↓F

↓

ω400
缶

↓

毒300

↓

V＝映ψ。

D

C

↓

↓／・

二、。。

…

↓

↓

R＝50・ψ。・・

缶

﹁2

↓

E

〈

巴
↓
。。lOO ↓

1↓
吊

1↓

V＝0剛蝋

↓

O

A〔。、r．
h＝Oom1・．

O

lO0

o0R

8
V＝Oom／s㏄

200

300

400

500

600

700

800

FiniteElement

Network

DlS1ANCE．CENTlMETERS
図5
Fig．5

解析例の境界条件と領域の要素分割
Cross

Sectionofan（1Boun（1arv

Ana1ysis

Regionwith

Conditions

境界ABは・計抑：例1ではr、境界になりん＝0cm

H，0．計算例2ではr，境界になり

V＝0cm／sec．計算例1，2で不変な境界BCDEは不透水屑，すなわち，r，境界であり，V＝0cm3／secとする・境界EFAは雨に打たれる面で降雨強度50mm／h，すなわち，、4境界であり，R＝50mm／hourとなる．物理定数は全浸透流域において，飽和透水係数K高＝κ嘉＝1．6×10−2cm／sec
K鳶＝κ嘉＝0cm／sec

（45）
（46）

O，8

一30

ぶ

⁝⁝

比透水係数

州一（㌃

3

（47）

E
U−20
斗

O．6

ぎ

o．4

飽和度

O

生

山一10

S＝1

ψ≧0

O．2
（48）

∫＝0．3056arctan（0．2668ψ十5．78）
O．6

十0・58

ψ＜0

（49）

とした．図6に比透水係数と飽和庇，圧力水0．8

図6ψとκ『の飽和度による変化
Fig・6V・・i・ti…fψ…1κ
of，Vater Saturation

頭と飽和度との関係を図示した．

侵透流域の要素分割を図5に示す．全要素

一62一

1．O

w＾IERs＾τUR川㎝、S・ithD・g・㏄

不飽和浸透を考慮Lた降雨浸透の有隈要素法による一解法一大倉500
ζmト120

400

7

ω400
缶
缶

300

Σ

F300
竃

二

200

2←z200

■

＜
＞
1」」

／

I

／

100

」

山100

■■

・

0

、cAL、、。、、L．wvEL．C、丁首

■

！／

100

200

300

400

DISTANCE，C

500

ENTIM

600

700

800

ETERS

図7全水頭の分布と流速（矢印を省略）．＜計算例1＞
Fig．7

Contoursof（Arrow

Total

Heads

Water

are

Head

and

Vectorof

e1imina七ed）．＜Ana1ysis

F1ow

Velocity．

FxampIe1〉

O．6

50◎

0．4
0．6

0，8
1．O

ω400
匝
山
←
1」」
Σ

；300
Z

O．8

u」

ア

○

オ

◎200
胃

1．O

06

、〃

〈
＞
u」

一

Lu1OO
○島

100

200

300

400

DlSTANC

図8
Fig．8

E，C

500

ENT−ME

600

700

800

TERS

飽和度の分布．＜計算例1〉
ContoursofDegreeof，Va七er

Saturation．

一63一＜Ana1ysis

Examplel〉

国立防災科学技術セソター研究報告第18号

1977年11月

500
cmト120

400

ω400

匝
山
←
u」

300

Σ

／／

巨300
Z

u」

／

200

ゴ

o200

／

／／

i＝

／
■

事
巴
LL1100

／■

○

／

／

■
ノ

■
lOO

／

O
図9
Fig．9

■

／

■

／

■

■

／
■
■
■
／ ／■
／

lO0

1

SCAL、、。、FL．wv，L．C，T11竺

／

200

300

400

DlSTANCE．C

500

ENTlMETERS

600

700

800

全水頭の分布と流速（矢印を省略）．＜計算例2＞
Con七〇urs（Arrowof

Heads

Tota1Water

are

Head

and

Vectorof

e1iminated）．＜Analysis

F1ow

Ve1ocitv．

ExampIe2〉

500

O．6
O．4
0．6
0．8

．40．4

κ

ω

缶400

缶
…
←

話300
二

○暗

看

疹、．。

←200
≦
巴
1」」

100

0

100

200

300

400

500

600

700

8∞

DlSTANCE，CENTlMETERS
図10飽和度の分布．＜計算例2〉
Fig・10C・・士・・…fD・g・…fW・t・・S・t…七i・・．＜A・・1y・i・E・・mpl・2＞

一64一

不飽和浸透を考慮した降雨浸透の有限要素法による一解法一大倉

数はユ60個，全節点数は100個となる．5章のアルゴリズムに必要なバラメーターは
ゴ＝0．5

Hα＝20c皿H20

ゴ、＝150

ゴ、＝10

α＝D．6

ε＝2×10■4ψ。

＝＿40cエnH．O

ψ3＝＿20c蛆H，0ψ、

ψ4＝＿30c蛆Hρ

＝＿17cm−H．O

ψ51＝＿15c血H20として計算した．ここで，人為的に決定する．これらのバラメーターと要素分割は物理的には解に影響を与えない．計算例1の解析結果を図7，図8に示し，計算例2は図9，図10に示す．図7と図9において，流域中の線分は流速のベクトルをあらわすが，矢印を省略した．

計算例1と2はともに，湧出点付近を除き，飽和度0．8以上の領域では地下水が斜面に平行に流れ，飽和度0．8未満の領域では40度の急勾配斜面にもかかわらず鉛直方向に土中水が流れる．白由水面（図8，図10で∫＝1．0の等飽和度線）の形態が計算例1，2で異なる・
解析的考察（Kono，1974）によると，湧出点において白由水面と境界のなす角度は，流域外から見て境界がオーバーハソグしている場合は，90度，他の場合はO度となる．故に，計算例1は境界ABと白由水面が垂直に交わり，計算例2は湧出点で境界AFと自由水面の傾斜が等しくならねぼならない．しかし，言十算例では有限要素法の離散化誤f乍と要素分割の適否により，湧出点の境界と自由水面のなす角度が理論からずれている．

図6より，∫＝0．6付近でψが変化するとK『が大きく変化し，非線形性が強いことがわ
かる．不飽和浸透を考慮した有限要素法によるロックフィルダムや河川堤防の解析例の多く

は降雨浸透を考えず，図1のr。境界のような水位を有する水源からの浸透を対象にしている．この場合，不飽和域を通過する浸透流量の全体の浸透量に対する割合が少なく，たとえ，不飽和域で非線形性が強くても，非線形性の流れ全体に及ぽす影響はさほど大きくない．しかし，水位をもつ貯水池や河川が存在せず，降雨のみを水源とする本報告の計算例
は，斜面土中を流れる浸透流の大半が非線形性の強い不飽和域を通過しており，不飽和浸透の流れ全体への影響ヵミ大きくなる．このため，計算例は，前述のロックフィルダムや堤体の解析例に比べはるかに不安定になっている．にもかかわらず，5章のアルゴリズムを用いる
ことにより，解が安定して得られた．

7．結語

降雨による土中水と地下水の定常浸透流を不飽和浸透を考慮して境界値問題に定式化した．この境界値問題はラプラス型の非線形偏微分方程式になるが，雨水の浸透する境界と土

中水・地下水の湧出する境界で第1種の境界条件（水頭が既定）か，第2種の境界条件（浸一65一国立防災科学技術センター研究報告第1・号1…年11月

透水量が既定）かが事前に定まらない場合があ／，解析解を得るのが困難になる．このた二ガレルキソ有隈要素法を用いて境界f直問題を非線形連立方程式に離散化し，これを線形計算を繰返す反復法で解いた・ここで提案された反復法は，不飽和浸透域で非線形性のかな

り強い場合や・やや不安定嫡界条件をもつ場合でも安定して収束解が得られた．これは，非線形性が強くかな／不安定な他の問題への有隈要素法の応用を示唆している．

参考文献
ユ）赤坂隆（

2）鰯勉（

3）駒田

967）：数値計算・応用数学講座第・巻，コロナ杜，1・・。，。。1、、．

972）：熱伝導と熟応力・一ンピー一タによる構造工学講座皿一・．。，培風館，、5、、．

金沢（975）：フ／ルダムの貯水池水位急降下時の非定常浸透流解析，および安定解析，土木学会論文報告集第240号．
4）Kono・L（974）：Fi・it・・1・…t・ml・・i…1・t・・・・・・・…1・・i。・r。、。wat，rF、。W，。。。。．。。

1∫CE，No・228，AUG．1974，ユ（）9一ユ16．
5）Neuman（1971）：A・・1・・i・・f・・・・・・・…1…it・・・・…。。f。。、。、i，g、、、F、、、t，E、、m、、、

Method・Wα伽R舳κθ・・舳・仏・・．・，・・．・，1。。。1・。1，。11．。。。．

6）山田．小田（975）：有限要素法1こよる非ダルシン浸透流の一解法，土木学会論文報告集第23、号．

ユ11pP．

（1977年7月21日

一66一

原稿受理）

不飽和浸透を考慮Lた降雨浸透の有限要素法による一解法一大倉

付録1記号
本文に用いられた記号を以下に示す．
ん㎜

境界値問題を離散化して得られる非線形連立方程式の係数行列．

月仰；；1

第ゴー1次反復解より求めた連立1次方程式の係数行列．

㎜ん一㎜κより得た連立1次方程式の係数行列．

B冊

境界値問題を離散化して得られる非線形連立方程式の定数ベクトル．

B冊｛■1

第｛一1次反復解より求めた連立1次方程式の定数ベクトル．

、B冊

抑、ゐより得た連立1次方程式の定数ベクトル．

6批巴％、の定数係数（付録2参照）。

C三肌

ゴー1次反復解と収束解との誤差の近似値を与える連立1次方程式の係数行列．

c椛e％巴の定数係数（付録2参照）．

D工

ゴー1次反復解と収束解との誤差の近似値を与える連立1次方程式の定数ベクト
ル．

6■1より肌ゐを作る任意定数．

H

既知水頭，既知全水頭．

ん水頭，全水頭．

Hαβ｛を求めるための任意定数．

1工1

既知水頭の節点々における水頭．

乃冊

節点〃の水頭．

乃．、｛節点〃の水頭のゴ次反復解．

、、κ〃■1と♂より作られるベクトル．

η、乃肌㎜ゐの第〃要素，肌々，。＝乃，、｛一

、、、石、

ゴ

十δ㎜冊・ゴ

Sign（

一ん，、｛

1）．

連立1次方程式，肌λ，、召mん、＝肌凪の未知数．

非線形連立方程式の反復解法の反復回数・

ゴε

非線形連立方程式の反復解法の許容反復回数．

ゴ品

非線形連立方程式の反復解法の単純な修正に対する許容反復回数．

K。ゴ

透水係数．

珊

飽和透水係数．

K『

比透水係数．

Kε「要素θの重心における比透水係数．

〃毒境界の外向き単位1法線ベクトノレ．

ρ皿八

非線形方程式の定数項B，、を構成する数．「。境界より派生する．

Q，、ム

非線形方程式の定数項B犯を構成する数．ム境界より派生する．

一67一

国立防災科学技術セソター研究報告

ダルシー流速の幻座漂方向の成分．

R

雨量強度．

s

飽和度．

∫召

の

第18号

！977年11月

要素θの重心の飽和度．

Signsign（Z）＝Z／lZ1

＝O

（Zキ0）

（Z＝0）

V

「。境界より浸透流域に単位時問・単位面積に流入する水量．

均

ゴ方向デカルト座標．

κ3

鉛直方向デカルト座標．

ツ肌

乃肌也■1と収束解伽との誤差の近似値．

｛一1
ツ㎜α皿

ルの絶対値の最大値．

α

反復解法の修正係数（5章操作7参照）．

β｛1

反復解法の修正係数（5章操作8f）参照）．

r

境界条件が既定される境界．

1□。

水頭伽の値が既知（乃、、＝H）の境界．

r．

単位時問・単位面積の流入水量が既知な境界．

11。。

要素θの「。境界となる辺．

∠

要素θの面積．

δ舳

クロネッーカーデルタ

δ㎜れ＝1

（舳＝〃）

＝0

（舳キ〃）

ε

非線形連立方程式反復解法の許容誤差．

∠

境界条件が既定されない境界．

■1

結果として水頭が既定された■境界．

ム

結果として浸透水量が既定された■境界．

∠。｛

ゴ次反復解より得られた■。境界の範囲．

■。｛

ゴ次反復解より得られたム境界の範囲．

％

んを線形結合，乃＝伽％で近似するガレルキソ法の基底関数．

％ε

要素θより派生して％を構成する関数．

ψ

圧力水頭．

ψ！

反復法による解の振動・発散を防ぐパラメーター（ψ！＜0）．上昇順に，ゴ：1，ゴ＝2，…，とする．

ρ

解析領域，解析流域．

ρ召

要素θの占める領域．

ム

節点〃を頂点に持つ要素の占める領域．

一68一

不飽和浸透を考慮した降雨浸透の有限要素法による一解法一大倉

付録2伽、の関数型
図A皿．1に示すように，要素θの占める領域をρ、，要素θの3節点を1，舳，〃とし，3節点の座標をそれぞれ（灼1，π1、），（∬㎜・，κ舳）・（κ舵・・κ冊・）とする・gπ侶は9色の外でO・9εの中でπ。，九ヨの1次関数とする．

伽（舳）一士（＾沽篶）（κ1，x3）∈9、

11

1

x王、

∠＝＿1

κ肌1

1

冊。

2

％、（

1，κ3）＝0

・（・・1…づ（A1I，1）

仙π棚（A皿．2）冊。
（κ1，κ3）守ρ召（AII．3）とする・ここで・」は要素1の面積にな工（・1l，・1・）る．％の値が節点〃で1・辺1舳で0にな

m（・・1。…）

図A皿．1

るから
α冊ε十わ閉κ11＋c冊日κ13＝0
α冊、十5帆、κ肌1＋c肌、κ舳＝0（AI［。4）

α閉十ろ冊επ肌1＋c閉κ榊＝2∠
α犯。，ゐ冊。，c冊召について解くと

α肌、＝κ王1π肌3＿γ伽凶3（A皿．5）

5珊、＝π工3一κ舳

（A］I・6）

c冊、＝κ肌1＿κ王1（A皿．7）

灼、，％、も同様に求めると，従属関係
灼、十％、十？冊、：1（五、，κ、）∈9、（A］I・8）が成立する．

付録3第f次反復解と真の解との誤差を求める修正法
係数月肌㎜が未知数12パノ＝1，2，…，舳，…，〃，…）の単調連続関数となる準線形連立1次方

程1式
λ冊肌（んゴ）ん肌；石、（んゴ）（A皿。1）

を考える．近似解をルとし，ルと真の解幻との誤差をεゴとする．
勾＝ツゴ十εゴ（A皿．2）（A皿．1）式のんバこツゴを代入し，2肌を未知数とする連立1次方程式
λ肌肌（ル）2肌＝B冊（ル）一69一（A皿・3）

固立防災科学技術セソター研究報告第18号

1977年11月

を解くと，2腕はルの関数となる．
2刷＝2腕（ル）（A皿．4）

上式のルに〜を代入すると2πの定義より
乃帆＝ル十επ＝z帆（ル十ε1）（A皿．5）

Tay1ar展開して，ε｛の二次以上の項を無視すると∂2

ル十ε㎜＝2π（ル）十Σユε｛1∂ツ｛（A皿．6）となる．故にε4はoπ灼＝δ肌（A皿．7）

・一・一1誇一1剛
あ肌＝カー2（ル）（・皿・）（A皿．9）を解いて近似的に求められる．ここで，δ㎜1はクロネッカーのデルタをあらわす．βを適当な定数として，ル十βε5を新たにルとして反復を繰返すと，ルの〜への収束が予想される．

一70一