1. No category

５．差分法と有限要素法による圧密解析
層厚換算 51.doc
５．１ 層厚換算法と多層地盤の有限要素一次元圧密解析
コンピュータが容易に使えないような場合,軟弱地盤対策指針では、多層地盤の一次元圧密解
析に層厚換算法の利用が規定されていると思われる.漸増載荷に対する圧密解析でもTerzaghi
理論の適用には, 慣用的解析には,かなりの工夫が必要であり,限界もある.しかし,最近では能
力の高いパソコンが誰でも簡単に利用できる。層厚換算法を必要とするような場合や盛土の段
階載荷の場合には,地盤条件や荷重条件に対しより厳密に対応できる有限要素法や差分法など
の数値解析によるべきと考える.パソコン能力の向上で計算結果の図化も容易である.ここでは,
２層地盤の一次元圧密問題に対し,有限要素法と差分法そして層厚換算法による慣用解析法と
による数値解析結果を比較している.また,解析に用いる土モデルの線形あるいは非線形かによ
る影響も検討する.
（１）地盤ならびに解析条件
片面排水条件で２層地盤の一次元圧密を検討対象とする.上層の上端排水面,下層下端非排水
面とする.２層の層厚をそれぞれ５ｍ、地盤を線形弾性体とする解析では,表１の土質定数を用
いた.載荷重増分DP=0.5 kgf/cm2とした場合,両層の圧密沈下量は,それぞれ50cm,20cm,合計沈
下量 70cmである.有限要素解析に用いたBrittoらのプログラムでは,弾性係数E,ポアソン比ν
と透水係数ｋが必要である.ν=0.3と仮定し,E=(1+ν)*(1-2ν)/(1-ν)/mv、k= mv * cv *γωか
ら計算した.
参考文献；Britto,A.M. and Gunn,M.J. ： Critical state soil mechanics via finite elements,
ELLIS HORWOOD LIMITED,1987.
A層
B層
表１ 土質定数
層厚 体積圧縮係数 圧密係数
弾性係数 E
0.1 cv 3.71=364.4kPa
500
0.2
mv
500
0.08
0.025
9.29=910.9
2
2
cm /min
kgf/cm2
cm
cm /kgf
透水係数 k
3.33*10-9
3.33*10-10
差分法による２層地盤の計算では、両層における時間増分を等しくすると両層の要素長間に、
次式の関係がある.A層の要素長 l A は、下層のそれ l B の２倍にする必要がある.
しかし,有限要素法では,そのような関係の必要性はない.
c
l
0 .1
A
vA
lB
c vB
0 . 025
2
ここでは,表１に示した条件と表２の要素数で有限要素法と差分法で計算し, 層厚換算法に
よるそれと比較考察する.表２の要素数（３）では,A、B層の要素長は,１mで共通であるが,要素
数（１）、（２）では、A、B層の要素長比 l A / l B 2 A層の要素長は, 0.5mと1mである.
解析では, 各層の要素長は均等である.要素数（３）では, 差分法による上記の要素長比を無
視している.
表２ 計算に用いた要素数
要素数（１）要素数（２） 要素数（３）
A層
10
5
5
B層
20
10
5
全層(A+B) 30
15
10
S2= 20
cm
20
40
60
One-D. Consolidation
calculated by FEM
Layer H cm
A
500
B
500
101
Number
of elements
10
15
30
E
k
-9
364 3.33*10
-10
911 3.33*10
kPa m/sec
102
S1= 50 cm
Settlement ( cm )
0
0
S2= 20
cm
40
103
104
Time ( day )
STC1M.smp
図1 圧密沈下量時間曲線１（FEM）
STC3MDE.smp
図３圧密沈下量時間曲線３
～ 層厚換算法との比較
Number
of elements
20
60
S= S1+S2
15
30
One-D. Consolidation
calculated by FDM
Layer H cm
A
500
B
500
mv
cv
0.2
0.1
0.08 0.025
cm 2/kgf
101
cm2/min
102
S1= 50 cm
S= S1+S2
103
104
Time ( day )
STC2D.smp
図２圧密沈下量時間曲線２（FDM）
Settlement ( cm )
Settlement ( cm )
（２）計算結果と考察
図１と図２に、それぞれ有限要素法(FEM)と差分法(FDM)による圧密沈下量時間曲線の計算結
果を示した。両図から明らかなようにA、B層の圧密沈下量が正確に計算されている。FEMにより
計算された図１では, 要素数によらず計算結果はほぼ同じである.また,地層が変化しても,各
層間の要素長比を任意に選択しても同じ計算結果が得られることが,注目される.FDMによる図
２の計算結果では,圧密初期の沈下量に小さなズレが認められる.FDMでより正確な計算結果を
得るには,ある程度以上の要素数が必要と思われる.また, 計算に用いたFEMとFDMプログラム
は,TerzaghiのU-Tv関係でその妥当性を検証している.
0
S2= 20
cm
20
Equivalent
thickness method
40
S1= 50
cm
One-D. Consolidation
calculated by
FEM
& FDM
60
S=
S1+S2
Number of elements = 30
101
102
103
104
Time ( day )
図３は要素数３０で実施したFEMとFDMの計算結果（黒実線と赤点線）に 印で層厚換算法に
よる計算結果を加えた.層厚換算法では, ２層地盤における層境界の水の流れの連続条件を考
慮していないため, 圧密後半の圧密沈下量時間曲線が大きく食い違うことがわかる.
ある圧密度における圧密層内の過剰間隙圧分布を図４と図５に示した.それぞれFEMとFDMに
よる計算結果である.FEMによる計算結果は、要素数の影響を受けないが, FDMによるそれは, 要
素数が少ないと圧密初期の間隙圧が小さく計算される.その結果が図２の圧密沈下量時間曲線
２（FDM）の計算結果に対応している.図６に示した要素数３０の両方法による計算結果は, ほ
ぼ一致しているが,、圧密度が小さい場合, 若干FDMで計算した過剰間隙圧分布が小さい.また,
図６には層厚換算法による換算層で単一層とした場合の圧密度５０％における過剰間隙圧分布
の計算結果を 印で記入している.
Pervious
0
Number of 30
elements 15
10
Depth ( m )
Depth ( m )
0
2
4
Layer
A
6
B
Ue=0.1
Ue=0.7
8
Ue=0.5
Ue=0.3
Pervious
Number of 15
elements 30
2
4
Layer
A
6
B
Ue=0.1
Ue=0.7
8
Ue=0.3
Ue=0.5
Impervious
Impervious
10
10
0
10
20
30
40
50
Excess pore pressure ( k Pa )
Uty1.smp
図４ 過剰間隙圧分布１（FEM）
0
20
30
40
50
Excess pore pressure ( k Pa )
uty2.smp
図５過剰間隙圧分布２（FDM）
Depth ( m )
0
Uty３.smp
図６ 過剰間隙圧分布３
～ 層厚換算法との比較
10
Pervious
Number of
elements = 30
5
Ue=0.1
Layer
A
B
Ue=0.7
Ue=0.3
Ue=0.5
H2
10
Equivalent
thickness
Impervious
Calculated by
FEM
FDM
15
0
10
20
30
40
50
Excess pore pressure ( k Pa )
層厚換算法では,２層地盤の層厚５ｍの下層の層厚は,１０ｍになる.層厚換算後の圧密層全
体の平均圧密度Ue＝５０％における過剰間隙圧分布を一例として示している. 印の過剰間隙
圧分布は,FEMとFDMによる黒実線や赤点線より大きく圧密が遅れていることが明らかである.
層厚換算法では,多層地盤の層境界の水の流れの連続性を満足していない結果を示すものと考
えられる.
これまでの有限要素法による一次元圧密連成解析は,棒要素を採用している.一次元圧密問題
でも多次元要素で解析する場合も考えられる.Brittoらの三角形要素をFE解析プログラム
CRISP.forによる計算結果との比較を試みた.両有限要素解析とも変位を二次,間隙圧を一次形
状関数を用いている.１０要素（ただし,三角形要素では,１０列２０要素）で各層の要素長を均
等とした有限要素解析結果の比較を図７と図８に示した.
解析に用いた要素図を参考のため,本節末尾に示す.
CRISP.forによるものである.
棒要素は、Tiny.for 三角形要素は
S2= 20
cm
20
Rod
Trangle
element
One-D. Consolidation
calculated by FEM
40
Layer H cm
A
500
B
500
60
101
E
k
-9
364 3.33*10
-10
911 3.33*10
kPa m/sec
102
S1= 50 cm
S= S1+S2
Excess pore pressure ( kPa )
Settlement ( cm )
Number
= 10
of elements
0
50
u (B bottom)
40
u (A bottom)
30
20
Calculated
by FDM
( 15 elements )
FEM ( 10 elements )
Rod
Triangle
element
10
0
103
104
Time ( day )
101
102
103
104
Time ( day )
STC4M.ｓｍｐ
Utc1m.smp
図７ ＦＥＭによる圧密沈下量時間曲線４
図８ A,B層底部における過剰間隙圧
（棒要素と三角形要素）
の経時変化
0
mv
Cc
S2= 20
cm
20
40
60
Number
= 15
of elements
One-D.
Consolidation
calculated by FDM
S1= 50 cm
Layer Cc
f0 mv
A 1.993 6 0.2
B
0.513 3 0.08
2
cm /kgf
102
103
Settlement ( cm )
Settlement ( cm )
棒要素と三角形要素を用いたFEMによる計算結果図７と図８によれば,要素形状の違いによる
圧密沈下量時間曲線と両層底部の過剰間隙圧の計時変化に,ほとんど差は認められない.また,
図８には赤丸印でFDMによる計算結果を加えてある.その結果は,FEMによるそれと同じである.
また、以上の検討には,表１に示した土質定数,すなわちｍvやEによる線形応力ひずみ関係を
用いでいる.圧縮指数Ccや弾塑性モデル,非線形応力ひずみ関係を用いた圧密沈下量時間曲線の
計算結果を,層厚換算法によるそれと比較する.
Number
= 10
of elements
0
S2= 20
cm
20
Elastic
Elastoplastic
soil model
40
S1= 50 cm
One-D. Consolidation
S= S1+S2
104
Time ( day )
STC5DCc.smp
図９ FDMによる圧密沈下量時間曲線5
60
calculated by FEM CRISP.for
101
102
S= S1+S2
103
104
Time ( day )
STCM6MCc.smp
図１０ FEMによる圧密沈下量時間曲線６
図９ FDMによる計算結果では,図中実線で示した圧縮指数Cc,非線形土モデルを用いた圧密
沈下量時間曲線が,体積圧縮係数mv,線形によるそれより圧密沈下速度が速く計算される.この
計算結果は,一般的有効応力に関する圧密度を用いた圧密沈下量時間曲線とひずみのそれによ
るものに対応する.
三角形要素を用いた多次元圧密FE解析プログラムCRISP.forによる一次元圧密解析結果を図
q ( ｋPa )
１０に示した.赤丸印,弾性（線形）そして実線は弾塑性（非線形：修正 Cam Clay）土モデルに
よる計算結果である.図９FDMによる計算結果と異なり,土モデルが線形あるいは非線形かによ
る違いが計算されない.非線形土モデルのプログラミングに問題があると考えられる.精緻な土
モデルを用いたプログラムでもその結果が反映されていないようである.また,修正 Cam Clay
モデルには,モデル自身にも欠陥があることがよく知られている.一次元圧密では,水平方向の
有効応力が実際と異なる点である.図１１は,CRISP.forで計算された一次元圧密前後の有効応
力経路である.
η=1
70
Pq1.smp
Calculated by CRISP.for
Elastic
Elastoplastic
soil model
60
50
η=0.75
=const.
40
Layer A
図１１ CRISP.for（FEM）による
有効応力経路
Layer B
30
Before consolidation
20
20
30
40
50
60
70
p´ ( k Pa )
80
一次元圧密では圧密中静止土圧係数あるいは応力比η(=q/p)一定であることが知られている.
図１１中にη一定線を実線で示しているが,点線と破線で示したそれぞれ線形そして非線形土
モデルにより計算された有効応力経路は,ηが一定でない.計算された圧密沈下量時間曲線が実
測値とよく一致した土モデル（応力ひずみ関係）でも、計算された応力が実際と異なるのでは,
その妥当性に問題がある.その問題点として修正 Cam Clayモデルの場合,塑性ポテンシャルの
不具合が指摘されている.
また, 上記計算は, “ 層圧換算.bas ならびに ODCmv.bas ”で実行された.
参考図；
上端
y=0
排水面
Tiny.for の棒要素
変位３節点
水圧２節点
CRISP.forの三角形要素図
側面
y 変位
自由
x変位
拘束
非排水面
下端
y=H
非排水面
変位固定