應用廣義塑性模式進行鯉魚潭水壩數值分析之研究

(1)

應用廣義塑性模式進行鯉魚潭水壩數值分析之研究

蔡佩勳 ¹ 馮正一 ² 趙世傑 ¹ 梁正育 ¹

1 朝陽科技大學營建工程系 ² 中興大學水土保持學系 E-mail: [email protected]

摘要

本研究以廣義塑性模式為主軸，利用Fish語言將此模式編寫成FLAC可執行之副程式，配合數值分析軟體 FLAC程式進行數值分析。副程式完成後，模擬砂土進行三軸不排水單調壓縮試驗和反覆載重試驗，分析結果與試驗結果作比較，發現此模式能有效模擬砂土承受反覆載重之情形，而且也驗證Fish副程式之正確性。土石壩之模擬分析，本研究係以鯉魚潭水庫為例，主要是因為鯉魚潭水庫之監測資料齊全，本文從分層滾壓階段至蓄水後滲流分析，最後施以集集地震之加速度進行動態反應分析。由分析結果與監測資料或前人利用摩爾庫倫模式所得的結果比較，此模式可用於大地工程問題之分析上。

關鍵詞：廣義塑性模式、土石壩、數值模擬

一、前言

台灣位處於歐亞板塊與菲律賓板塊上，斷層帶繁多且地震頻繁，最近一次的災害性地震發生於民國88年9 月21日，因車籠埔斷層錯動而引發了規模7.2的強震，

而如此強烈的地震對全台各地造成嚴重的損失，對於公共工程也造成極大損害，例如斷層附近的石岡大壩，因地盤隆起使得混凝土壩壩身斷裂，由於土石壩之安全與否對於民生及經濟上都有重要的影響，因此如何評估其安全性是重要的課題。

本研究利用大地工程數值分析軟體FLAC進行土石壩二維反應分析，包括模擬築壩過程壩體的受力行為、

蓄水後壩體滲流的情況、蓄水後之動態反應分析。主要利用Pastor(1986)【1】等人提出之「廣義塑性模式」

（Generalized Plasticity model），配合FLAC分析軟體進行土石壩之數值模擬，因廣義塑性模式未內建於FLAC 內供使用者使用，故需另以Fish語言(FLAC可辨讀之程式語言)將此模式撰寫為一外掛程式，在進行分析時以呼叫方式代入主程式中進行分析，故必須對其正確性做檢核。本文將以三軸試驗來進行驗證，然後再將「廣義塑性模式」應用至土石壩之數值模擬分析。主要研究方法分述如下：

1. 資料蒐集與文獻回顧：蒐集與本研究相關之資料及與「廣義塑性模式」架構相關的研究，了解土石壩相關分析方法及其發展演進。

2. 撰寫土壤模式副程式：本研究所用之土壤模式為「廣義塑性模式」，由於並非FLAC內建之模式，因此需利用Fish語言先將此模式編寫成FLAC可執行之副程式，再由FLAC主程式以CALL的指令呼叫，執行分

析運算。

3. 程式驗證：使用「廣義塑性模式」進行三軸試驗模擬，模擬結果與試驗結果比對，以驗證程式之正確性。

4. 土石壩靜態反應模擬：建置土石壩網格，模擬分層填築之築壩過程，其後進行蓄水後之滲流與應力平衡分析。

5. 土石壩動態反應模擬：施以集集地震之地表加速度，模擬壩體受震過程，進行動態反應分析。

二、廣義塑性模式

廣義塑性模式由 Pastor 與 Zienkiewicz(1986)【1】

應用於大地工程材料中，其後Pastor 與 Zienkiewicz 等人延續先前研究，多次對模式加以改良，此模式之特點為不需特別定義出降伏函數與塑性潛能函數，並可模擬應變硬化或軟化之行為。

原始的 Pastor- Zienkiewicz model (P-Z model)

（1986）【1】已能夠模擬砂土在固定圍壓下，排水與不排水三軸試驗之行為，除此之外也能夠模擬反覆載重下之動態行為，但此時之模式尚未將砂土之應變硬化行為考慮進去，直至Pastor(1993)【2】才將土壤在反覆載重下的硬化行為加入此模式中。

一般材料之應力增量與應變增量之關係可寫為下式：

:

d σ = D

ep

d ε

………..（1）

上式中之D^ep如下式所示：

/

/ /

: : : : :

e T

ep e gL U

T e

e

L U gL U

D n n D

D D

H n D n

= −

+

……（2）

其中

d σ

：應力增量

d ε

：應變增量

D

ep：彈塑性勁度矩陣

D

e：彈性勁度矩陣

/

n

gL U：加壓與解壓塑性流方向向量

n

：加載方向向量

/

H

L U：加壓與解壓塑性模數

D e

彈性勁度矩陣如下所示：

(2)

4 2 3 3 0 0

2 4

3 3 0 0

0 0

2 2

3 3 0 0

e

G G

K K

G G

K K

D

G

G G

K K

⎡ + − ⎤

⎢ ⎥

⎢ − +

⎢ ⎥

=⎢ ⎥

⎢ ⎥

⎢ − − ⎥

⎢ ⎥

⎣ ⎦

0

⎥

p

………...（3）

土壤應變增量又可分為彈性應變增量與塑性應變增量，如下式所示：

d ε = d ε

e

+ d ε

………...（4）

d ε

e：彈性應變增量

d ε

p：塑性應變增量 2.1 加載階段

砂土在加壓階段下膨脹係數d_g如下式：

(1 )( )

p v

g p g

s

d d d

ε α M η

= ε = + −

………..（5）

其中：

p

d ε

v ：塑性體積應變增量

p

d ε

s ：塑性偏差應變增量

M

g：臨界階段p-q＇圖之原點與殘餘應力路徑切線之

η

：應力比斜率 q/p＇

α

：常數，通常以0.45 代入當

'

^g

q M

η = p =

時，則膨脹係數dg=0，而有二個階段可能出現dg=0 的現象，分別為特徵階段 (characteristic state)和臨界階段(critical state)。特徵階段之行為即砂土受力後，體積由壓縮進入膨脹之階段，臨界階段為砂土發揮殘餘強度且體積無變化時。

塑性流方向向量ng定義如下式：

gv g

gs

n n n

⎡ ⎤

= ⎢ ⎣ ⎦ ⎥

……….（6）

其中：

1

2 g gv

g

n d

d

= +

……….………….（7）

2

1

gs

1

g

n

d

= +

………（8）

加載方向向量n 如下式：

v

s

n n n

= ⎢ ⎥ ⎡ ⎤

⎣ ⎦

..………...（9）

其中：

1

2 f v

f

n d

d

= +

………（10）

2

1

s

1

f

n

d

= +

………（11）

(1 )( )

d

f

= + α M

_f

− η

……….（12）

M

f：由不排水三軸試驗決定。

在此模式中主要影響加載時塑性行為之塑性模數HL，其定義如下：

{ }

0

L f v s

H = H pH H + H H

_DM ……...….（13）

其中：

4 f

1

f

H η

η

⎛ ⎞

= ⎜ ⎜ ⎝ − ⎟ ⎟ ⎠

且

₁ ¹

f

M

f

η ⁼ ^⎛ ^⎜ _⎝ ⁺ α ^⎞ ^⎟ _⎠

…（14）

v 1

g

H M

= −

η

………...（15）

( )

0 1

exp

0

H

s

= β β − β ξ

………...（16）

max

H

DM

ζ

γ

ζ

⎛ ⎞

=⎜ ⎟⎝ ⎠且 ₁ ¹

p

1

M α η

α

ζ α

−

⎡ ⎛ ⎞ ⎤

= ⎢⎣ − ⎜⎝ + ⎟⎠ ⎥⎦

…...（17）

β0、β1：由試驗結果之應力應變曲線決定。

ξ：累積之偏差塑性應變量

ξ = ∫ d ε

_s^p

2.2 解載與反覆載重階段 解載時的塑性模數HU定義如下：

當 ^g

1

u

M

η ^>

^時，則 ⁰

u

g

U u

u

H H M

γ

η

⎛ ⎞

= ⎜ ⎟

⎝ ⎠

當 ^g

1

u

M

η ^≤

^時，則

^H

^U

⁼ ^H

^u⁰

H

U：第一次解載時其應力應變曲線的初始斜率

η

u：解載時的應力比

而解載時的塑性流方向向量定義如下式所示：

guv gU

gus

n n

n

⎡ ⎤

= ⎢ ⎥

⎣ ⎦

……….（18）

其中：

n

guv

= − abs n ( )

gv ………（19）

gus gs

n = + n

2.3 廣義塑性模式 Fish 副程式之驗證

將完成之廣義塑性模式Fish 副程式套用至 FLAC 中進行三軸試驗之模擬以驗證其正確性，而模擬之試驗採用Castro(1969)【3】對於砂土液化行為之研究結果，此研究為Castro 於哈佛大學發表之博士論文，主要利用三軸之不排水單調壓縮試驗和反覆載重試驗來觀察砂土在各種情況下之行為，故本論文將以此試驗資料進行驗證。

2.3.1 單調壓縮不排水三軸試驗

Castro 的試驗共分三個試體， a 試體相對密度 Dr=27%，b試體Dr=44%，c試體Dr=47%，所用圍壓皆為

(3)

4kg/cm²進行試驗，FLAC之模擬條件均參照Castro實驗內容，以應變控制的方式逐步施加速度使試體產生應變，而廣義塑性模式部分所需利用之參數則參考 Pastor(1990)【4】等人模擬此試驗時，所使用之參數如表1 所示，數值模擬結果與試驗結果如圖 1 所示。由圖可看出，廣義塑性模式對本試驗的模擬結果令人滿意。

表1 廣義塑性模式參數表

a b C

K0(kPa) 87.5 87.5 87.5 G0(kPa) 131.25 131.25 131.5

Mf 0.40 0.545 0.570 Mg 1.5 1.32 1.12

β0 4.2 4.2 4.2 β1 0.2 0.2 0.2 γ 1 1 1 γu 2 2 2 H0 350 350 350 Hu0 60000 60000 60000

α 0.45 0.45 0.45

0 1 2 3 4 5 6 7 8

0 1 2 3 4 5 6 7 8 9 10

Axial Strain (%)

Deviator Stress (kg/sq cm)

試體a(數值模擬) 試體b(數值模擬) 試體c(數值模擬) 試體 a(試驗結果) 試體 b(試驗結果) 試體 c(試驗結果)

（a）軸差應力與軸向應變之關係曲線圖

0 1 1 2 2 3 3 4 4 5

0 1 2 3 4 5 6 7 8 9 10

Axial Strain (%) Induced Pore Pressure (kg/sq cm)

試體a(數值模擬) 試體b(數值模擬) 試體c(數值模擬) 試體a(試驗結果) 試體b(試驗結果) 試體c(試驗結果)

(b)孔隙水壓與軸向應變之關係曲線圖圖1 單調壓縮不排水三軸試驗 2.3.2 無拉伸應力之反覆載重三軸試驗

此一試驗試體亦為上節所述之a試體(Dr=27%)，進

行反覆載重試驗，試驗過程採應力控制方式，反覆施加軸差應力，而所施加之軸差應力為0~1.44kg/cm²，圍壓為4kg/cm²，試驗結果如圖2。

0 1 2

0 5

Axial Strain (%)

D evi at or S tr es s (kg/ sq c m )

10

0 1 2

0 1 2 3 4

Effective Normal Stress (kg/aq cm)

Sh ea r S tre ss (k g/ sq c m )

（b）剪應力與有效應力之關係曲線圖

0 1 2 3 4 5

0 5

Axial Strain (%) Induced Pore Pressure (kg/sq cm)

10

（c）孔隙水壓與軸向應變之關係曲線圖

圖2 無拉伸應力反覆載重三軸試驗(Castro 1969)

FLAC 模擬結果如圖 3，比較圖 2 與圖 3 發現，兩者的結果也近乎相等，所以廣義塑性模式可合理模擬砂土在反覆載重之力學行為。

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

0 1 2 3 4 5 6 7 8 9 10

(4)

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

0 1 2 3 4 5

Effective Normal Stress (kg/sq cm)

D evi at or S tr es s (kg /s q cm )

（b）剪應力與有效應力之關係曲線圖(模擬結果)

0 1 1 2 2 3 3 4 4 5

0 1 2 3 4 5 6 7 8 9 10

Axial Strain (%) Induced Pore Pressure (kg/sq cm)

（c）孔隙水壓與軸向應變之關係曲線圖

圖3 無拉伸應力反覆載重三軸試驗(數值模擬結果) 2.3.3 反覆載重三軸試驗

同樣以a試體(Dr=27%)為進行試驗之對象，應力加載的方式亦為應力控制，而軸差應力加載之範圍為 1.4kg/cm² ~ -1kg/cm²，圍壓採用相同的4kg/cm²，試驗之結果如圖4 所示。數值模擬仍採應力控制方式施加軸差應力，軸差應力加載之範圍為1.4kg/cm² ~ -1 kg/cm²，圍壓為4kg/cm²，數值模擬結果如圖 5。由圖 2 與圖 3 發現，兩者的結果近乎相等，所以廣義塑性模式可合理模擬砂土在反覆三軸載重之力學行為。

-1.5 -1 -0.5 0 0.5 1 1.5 2

-2 0 2 4 6 8 10

圖4 反覆載重三軸試驗(Castro 1969)

-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5

-1 0 1 2 3

Axial Strain (%) D e vi a tor S tres s (k g /sq c m )

4

-2 -1 0 1 2

0 1 2 3 4 5

Effective Normal Stress (kg/ sq cm)

S h e a r S tre ss ( k g /s q cm )

（b）軸差應力與軸向應變之關係曲線圖

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5

-1 0 1 2 3 4 5

Axial Strain (%) D e viator S tr e ss (k g/ sq c m )

（c）孔隙水壓與軸向應變之關係曲線圖圖5 反覆載重三軸試驗

三、分析結果與討論

由三軸試驗之模擬結果可以證明「廣義塑性模式」

套用於FLAC 程式中，可有效模擬土壤在單調壓縮載重或反覆載重後之行為反應，且有其準確性。因此，本文將再利用此模式於土石壩之數值模擬上。土石壩之分析步驟如下：

1. 靜態分析：為模擬土石壩之滾壓填築過程，網格將以分層方式依序向上堆疊以與現地施工狀況相符，

而分析之步驟為將網格採取彈性模式進行分析，計算其初始應力，其後將上下游殼層網格之分析模式更改為「廣義塑性模式」，心層則採莫爾庫倫模式，

再進行計算以模擬土壤之塑性行為。

2. 滲流分析：當壩體滾壓構築完成後，為模擬其蓄水後之滲流所產生之孔隙水壓與有效應力，將再進行滲流分析，而蓄水高度以水壩之滿水位來進行模擬。

3. 動態反應：選用地震加速度為 921 集集地震發生時於水里國小測站所監測之地震加速度資料，將此加速度施加水壩底部，以探討其在遭遇此等地震後壩體之動態行為。

3.1 廣義塑性模式之材料參數計算

鯉魚潭水庫之廣義塑性模式材料參數之訂定，係以鯉魚潭水庫施工中安全性評估計畫（1992）【5】中之壓密不排水三軸試驗(CU Test)相關結果為依據，並以 FLAC 模擬其試驗，調整相關參數使模擬結果與試驗結果相符，以得到廣義塑性模式之相關參數，模擬與試驗結果如圖6。

(5)

0 1 2 3 4 5 6 7 8

0 2 4 6 8 10 12 14 16 18

模擬結果(圍壓1kg/sq cm) 試驗結果(圍壓1kg/sq cm) 模擬結果(圍壓4kg/sq cm) 試驗結果(圍壓4kg/sq cm)

-2.00E-01 0.00E+00 2.00E-01 4.00E-01 6.00E-01 8.00E-01

0.00E+00 5.00E+00 1.00E+01 1.50E+01 2.00E+01

In d u ced P o re P re ss u re (k g /s q cm )

試驗結果(1kg/sq cm) 試驗結果(1kg/sq cm)

模擬結果(4kg/sq cm) 模擬結果(4kg/sq cm)

圖6 鯉魚潭水庫壓密不排水三軸試驗模擬與試驗結果 FLAC 分析時之網格與材料區分如圖 7，所算出的材料參數如表 2，因李俊男(2005) 【6】曾針對鯉魚潭進行相關研究，故本案例心層之莫爾庫倫之材料參數採用李俊男使用之數值，本文之數值模擬結果也將與李俊男(2005)之相關分析成果進行比較。

圖7 網格與材料區分圖（鯉魚潭水庫）

表2 鯉魚潭材料參數表

3.2 鯉魚潭水庫靜態反應分析

壩體經靜態分析所得之沈陷量如圖8，最大之沈陷

量分別發生在上游與下游殼層約 15cm，水平位移量如圖9，其最大位移量約為 15cm，與李俊男(2005)相比較如圖10、圖 11，最大沈陷量部位皆為上、下游殼層部分，李俊男(2005)所分析之結果最大沈陷量發生處較接近於心層部位約 30cm，而本研究最大沈陷量發生位置較靠近殼層的中央位置，沈陷量也較小，就水平位移量部分來觀察趨勢大致一致。

圖8 靜態平衡後 Y 方向變位量（鯉魚潭水庫）

圖9 靜態平衡後 X 方向變位量（鯉魚潭水庫）

圖10 靜態平衡後 Y 方向變位量（鯉魚潭水庫）

【李俊男(2005)】

圖11 靜態平衡後 X 方向變位量（鯉魚潭水庫）

【李俊男(2005)】

3.3 鯉魚潭水庫滲流分析

分析之結果如圖12，與鯉魚潭管理中心所監測之孔隙水壓分佈情形相比較，如圖13所示，其孔隙水壓分部之大小與趨勢皆為相近似。

圖12 孔隙水壓分佈圖（鯉魚潭水庫）

莫爾庫倫模式(Mohr-Coulomb Model) 區

域 γ (kg/m³)

E (Pa)

υ c (Pa)

φ (°)

permeability (m²/(Pa sec))‧ 1 2100 1.08E9 0.34 14500 41 5.0E-12 2 2100 1.65E9 0.45 83500 27 1.02E-13 3 2200 1.04E9 0.3 30500 35 5.0E-12

廣義塑性模式(Generalized Plasticity Model) 區

域 K0 G0 Mg Mf αg，αf β0 β1 H0 γ，

γu

1 140 190 1.4 0.6 0.45 2.5 0.02 25000 1 3 140 190 1.4 0.6 0.45 2.5 0.02 25000 1

Finn Model 區

域 γ

(kg/m³) E

(Pa) υ c

(Pa) φ

(°) ff_c1 ff_c2 ff_s wit ch 1 2100 1.08E9 0.34 14500 41 0.0185 2.157 1 3 2200 1.04E9 0.3 30500 35 0.0185 2.157 1

1 2 3

(6)

圖13 孔隙水壓力分佈圖(鯉魚潭水庫定期監測分析,2000)【7】

3.4 鯉魚潭水庫動態分析

本研究分析之結果如圖14~16，李俊男的分析結果如圖17~19所示。x方向變位最大值約為1.25cm，發生位置在下游殼層部分，整個壩體有向下游移動之趨勢。最大沈陷量發生於壩頂部分約1.25cm。

圖14 動態分析之 X 方向變位量

圖15 動態分析之 Y 方向變位量

圖16 動態分析之孔隙水壓分佈圖

圖17 動態分析之 X 方向變位量【李俊男(2005)】

圖18 動態分析之 Y 方向變位量【李俊男(2005)】

圖19 動態分析之孔隙水壓分佈圖【李俊男(2005)】

四、結論

本研究以廣義塑性模式對鯉魚潭水庫進行數值模擬，與前人所作之研究進行比對，以確定此模式之正確性，歸納出下列結論：

1. 使用「廣義塑性模式」來進行三軸試驗模擬時無論是單調壓縮載重或者反覆加載皆與試驗結果相符，具有不錯之模擬能力。

2. 本研究之案例二鯉魚潭水庫部分，以廣義塑性模式分析結果與李俊男(2005)【6】之研究比較，本研究以廣義塑性模式所得之結果，水平位移量與沈陷量趨勢一致。

五、參考文獻

【1】 Pastor, M., and Zienkiewicz, O. C., “A generalized plasticity hierarchical model for sand under monotonic and cyclic loading,” Proc. 2^nd Int. Conf.

on Numerical Models is Geomechanics, Ghent(Belgium), G. N. Pande and W. F. Van Impe (eds), M. Jackson and Son Pub., pp. 131-150, 1986.

【2】 Pastor, M., Zienkiewicz, O.C., Xu, G.D., and Peraire, J., “Modeling of sand behavior: cyclic loading, anisotropy and localization,” Modern Approaches to Plasticity, D. Kolymbas, Ed., Esevier, pp. 469-492, 1993.

【3】 Castor, G., “Liquefaction of Sands,” Havard Soil Mechanics Series, No. 81, Havard Unversity, 1969.

【4】 Pastor, M., Zienkiewicz, O.C., and Chan, AHC.,

“Generalized plasticity and the modeling of soil behavior,” International Journal for Numerical and Analyticl Method in Geomechanics, Vol. 14, pp.151-190, 1990.

【5】經濟部水利處中區水資源局，「鯉魚潭水庫施工中安全性評估計畫」，1992。

【6】李俊男，「土石壩地震反應二維數值分析」，碩士論文，朝陽科技大學營建工程系，台中，2005。

【7】經濟部水利處中區水資源局，「鯉魚潭水庫定期監測分析-總報告」，2000。

(7)

Implementation of generalized plasticity model in numerical analysis

of Li-yu-tan Earthdam

Tsai, P.H., Feng, Z., Chao, S.C., Liang, C. Y.

Department of Construction Engineering Chaoyang University of Technology

ABSTRACT

In 1986 Pastor and Zienkiewicz brought up the

“Generalized Plasticity Model” and also applied it into the geotechnical engineering issue. This research used the generalized plasticity model as soil constitutive law program of FLAC, a kind of numerical analysis software and based on the finite difference method to perform the analysis. Because this model has not been built in the program of FLAC, it is necessary to use the Fish language to compile this model into a subprogram that can be implemented by FLAC. After completing the subprogram, this research then verified its rationality by simulating the triaxial tests for comparison reason. We compared the numerical result of this study and the result of Castro’s (1969) experience study on sandy soil, the result is very close to that of Castro with the identical condition. Then, we used the user-written constitutive model to study earth-rockfill dam simulation analysis. For the Li-Yu-Tan Reservoir, the numerical result of this study will be compared with those of former researches in the construction stage and the dynamic stage. It shows that little differences between the result of this study and that of other research. It is further verified that this model is feasible in for the geotechnical engineering.

Keywords: generalized plasticity model, earthdam,

numerical simulation.

應用廣義塑性模式進行鯉魚潭水壩數值分析之研究