土壤應變軟化行為之數值模擬研究

IN GEOTECHNICAL ENGINEERING IN TAIWAN I-Lan

A39-1

土壤應變軟化行為之數值模擬研究

林佩君

¹

張世龍

¹

蔡佩勳

²

1

朝陽科技大學營建工程系碩士班研究生

2

朝陽科技大學營建工程系助理教授

摘 要

本研究以Desai 提出之「擾動狀態概念」（Disturbed State Concept, DSC）為土壤組 成模式，配合有限差分數值分析軟體FLAC 程式來模擬土壤受剪時應變軟化行為，討 論有關擾動狀態概念的材料參數決定方法並比較數值結果與實驗值之差異。本研究討 論渥太華砂在壓密排水試驗之力學行為，分別探討渥太華砂(ASTM C-109)在三種不 同的圍壓(69kPa、207 kPa、345 kPa)下以及不同的鬆緊程度(相對密度 Dr 為 20%、90%) 之受剪行為，以得到渥太華砂在不同鬆緊程度及圍壓下的 DSC 材料參數。由本研究 結果發現，本研究所採用的 DSC 參數決定方法所得到的材料參數，可預估砂土應變 軟化的行為，且預估的曲線與實驗結果大致相符。

關鍵詞：擾動狀態概念、數值分析、應變軟化。

Numerical Simulation of the Strain-Softening Behavior of Soil by DSC Model ABSTRACT

The primary objective of this study is to simulate the strain-softening behavior of Ottawa sand by Disturbed State Concept (DSC) model. The model was proposed to simulate the damaged behavior of material by Desai. In the study, the constitute law will be used to simulate the triaxial tests with the FLAC software. The triaxial test with the soil samples of Dr=20% and 90%, under 3 confining pressures were performed. A proposed method was developed to decide the DSC parameters from experiment results in this study. Then the material parameters of DSC model were draw into the numerical modeling to compare their differences. From the results show that the predict responses by the proposed method almost agree with the experiment ones.

Key Words: Disturbed State Concept, Numerical Simulation,Strain Softening

一、簡介

土壤在受外力的過程中，因承受外力使得土壤內 部產生擾動，材料的強度等力學性能便會降低，導致 材料和力學性能產生微觀變化。這種有關材料損傷的 研究，過去學者都曾在這領域進行過研究，土壤在受 剪時產生的力學反應如果以損傷的角度來探討是另 外一個研究土壤在工程應用上的方向，因此本研究擬 以DSC 模式來評估土壤受剪時之損傷行為，並建立 一套DSC 材料參數的決定方法，以供工程應用上作 為預測之參考。

Desai（2000）[1]發表 DSC 模式之理論及其相關 應用模式，他認為 DSC 模式能被用來模擬彈塑性問

題、熱黏塑性問題、與時間有關的疲勞和退化之問 題、及材料因溫度弱化之問題，還可模擬接合處和介 面行為的特性，在其研究也提出有關材料剩餘壽命和 微裂縫密度的計算，他利用水泥和砂混合製成的土壤 水泥試體和陶瓷二種材料，透過數值模擬與試驗結果 做比較，由其結果顯示出 DSC 模式可有效模擬材料 損傷後之行為。Varadarajan 與 Sharma（2003）、(2006) [2]、[3]對兩種土石材料進行三軸排水試驗，求出擾 動狀態概念(DSC)來決定材料參數，他們認為擾動狀 態概念應用於土石材料將可得到令人滿意的結果。

Desai 與 Chen（2006）[4]對擾動狀態概念材料參數的 進行敏感度分析，並使用最小平方差法對砂土的三軸 壓縮試驗結果來確定 DSC 模式的最佳化材料參數。

A39-2

2 黃秀雅（2008）[5]以 Desai 提出之「擾動狀態概念模

式」為主軸，配合以有限差分為基礎之數值分析軟體 FLAC 程式來進行構件損傷後之承載分析，探討不同 損傷程度與其承載力的關係，建立一套損傷構件的承 載力評估方法。利用前人結合超音波量測與鋼筋拉拔 試驗的試驗結果，配合基因演算法得到混凝土構件之 DSC 材料參數，再由此 DSC 材料參數進行損傷構件 的承載行為之研究。整理出三種徑高比分別為 1:2、

1:5 與 1:10 之混凝土試體抗壓強度與損傷函數 D 的迴 歸公式，建立一套運用超音波波速配合迴歸公式，即 可預測受損後混凝土構件之抗壓強度的評估方法。

本研究的主要目的分為以下兩點：

1. 建立可在 FLAC 軟體使用的—擾動狀態概念 的FISH 副程式，作為發展土壤損傷評估模式之用。

2.利用三軸試驗找出土壤 DSC 材料參數，包 括:γ、 β、 n、a1 、η1、A 、B 等七個參數。

本研究首先建立Desai 在 2000 年[1]所發表有關 評估材料損傷程度之「擾動狀態概念」—DSC 損傷模 式的Fish 副程式，因為此模式並未內建於 FLAC 軟 體供使用者使用，故需以Fish 語言(FLAC 可辨讀之 程式語言)將此模式撰寫成一外掛程式，以便在進行 FLAC 分析時，以呼叫方式代入主程式中進行分析。

然後由三軸排水試驗結果計算其DSC 材料參數。材 料參數來進行三軸試驗數值模擬，其分析結果並與實 驗數值比對。本研究之流程圖如圖1 所示。

圖1. 研究流程圖

二、DSC 模式理論背景

當材料在受到載重或環境外力時，會對材料的微 觀結構造成改變或擾動，部分的材料會從相對完好狀 態（Relative Intact State, RI state）漸漸調整到完全調 適狀態（Fully Adjusted State, FA state）以平衡此載重 或環境外力，圖2 為擾動狀態概念的示意圖。初始狀

態為相對完好狀態（RI state），而完全調適狀態（FA state）則為破壞狀態或極限狀態，材料的行為則是介 於此兩者之間。材料的行為端看材料中有多少比例處 於RI 狀態，多少比例處於 FA 狀態。而此比例稱為損 傷函數（Disturbance Function）D，其值介於 0～1 之 間，當D＝0 時表示材料完全處於 RI 狀態，材料行為 完全由RI 狀態的行為控制，當 D＝1 時表示材料完全 處於FA 狀態，其行為完全由 FA 的行為控制。

D 可以是塑性變形、初始條件、溫度和含水量的 函數並且可以與可量測的物理量建立關係，如應力、

應變及孔隙比等。

將應力寫成增量型態則為下式(1)：

( )

^{i c}

(

^{c i}

)

a D d D d dD

dσ = 1− _σ ⋅ σ + _σ ⋅ σ + _σ σ −σ 式(1) 以上標i、c、a 分別代表 RI 狀態、FA 狀態和觀 察到的狀態。

RI 狀態和 FA 狀態則分別可以套用不同的組合 律，以模擬材料不同的非線性行為。本研究以HISS-δ0 model（Desai, [1]）作為 RI 狀態，以臨界狀態來當作 FA 狀態。

圖2. 擾動狀態概念示意圖 ( 摘自 Desai, 2000 [1]) 當彈性階段進入塑性階段，必須滿足降伏準則

（yielding criterion）；材料在未降伏前仍滿足線彈性 之虎克定律且降伏函數小於零，若進入塑性，則降伏 函數等於零，且產生塑性變形與塑性應變。因此，總 應變增量可分為彈性應變增量與塑性應變增量之總 和，如下式：

e p

dε =dε +dε

式(2)

其中

dεe_{：彈性應變增量，}dε ：塑性應變增量 ^p HISS-δ0降伏函數之定義為：

土 壤 應 變 軟 化 行 為 之 數 值 模 擬 研 究

進 行 三 軸 C D 試 驗 從 渥 太 華 砂 中 以

兩 種 不 同 的 相 對 密 度 D r = 2 0 % 、 9 0 % 及 三 種 不 同 圍 壓 6 9 k P a 、 2 0 7 k P a 、 3 4 5 k P a

的 求 得 其 砂 土 的 相 關 參 數 。

分 析 結 果

計 算 完 D S C 所 需 參 數 之 後 套 入 三 軸 試 驗 結 果 進 行 比 對 ， 可 得 到 與 實 驗 相 似 的 良

好 的 比 對 結 果 。 計 算 土 壤 D S C 材 料 參 數 以 三 軸 試 驗 試 驗 結 果 ，

軸 差 應 力 - 軸 向 應 變 與 軸 向 應 變 - 體 積 應 變 及ٛ角 ， 來 計 算 砂 土 的 D S C 材 料 參 數 。

IN GEOTECHNICAL ENGINEERING IN TAIWAN I-Lan

A39-3

(

1

)

0

]

[ ^0.5

2 1 1

2

2 ⎟⎟⎠ − =

⎜⎜ ⎞

⎝ + ⎛

⎟⎟⎠

⎜⎜ ⎞

⎝

− ⎛

−

= _r ⁻

a n

a a

D S

P J P

J P

F J α γ β

式(3)

sγ

D

F J n J

1 1 2

2 1 2 ⎟⎟⎠

⎜⎜ ⎞

⎝

−⎛

= 式(4)

(

1−

)

^−0.5

= _r

s s

F β 式(5)

3/2 2

3

2 27

D D

r J

S = J

式(6)

1

η1

α = ξ ^a

式(7)

( )

∫ ^⋅

= d ε

^p

d ε

^{p 12}

ξ

式(8)

33 22 11

1

= σ

_ii

= σ + σ + σ

J

式(9)

ji

J σ

ij

σ 2 1

2

=

式(10)

mi km

J σikσ σ 3

1

3=

式(11)

6

2 1 2 2

J J

J _D= −

式(12)

3 1 2

1 3

3 27

2 3

2J J J

J

J _D = − +

式(13) 式中

J1、J2和J3為應力不變量 J2D、J3D：偏差應力不變量 Pa：大氣壓力

F：HISS-δ0模式之降伏函數 dε ：塑性應變增量 p

ξ

：累積塑性應變

a1、η1、n、 γ、 β 為 DSC 材料參數

一般彈塑性材料之應力增量與應變增量之關係 可寫為下式：

ε σ C d

d =

^ep

⋅

式(13) 如果材料的降伏函數等於塑性潛能函數即滿足 諧和流法則（Associative flow rule），此時 C^ep如下式 所示：

2 / 1

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

∂

⎟ ∂

⎠

⎜ ⎞

⎝

⎛

∂

∂

∂

−∂

∂

⎟ ∂

⎠

⎜ ⎞

⎝

⎛

∂

∂

⎟⎠

⎜ ⎞

⎝

⎛

∂

∂

∂

∂

−

=

σ σ ξ σ σ

σ σ

F F F C F

F

F C C F

C C

T e

T

e T e

e ep

式(14)

其中

d

σ

_{：應力增量}

dε ：應變增量

C

ep：彈塑性勁度矩陣

C

e：彈性勁度矩陣

式(15)中之

C

^e彈性勁度矩陣在平面應變的情形下 為：

⎥⎥

⎥⎥

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢

⎣

⎡

+

−

−

− +

−

−

− +

=

G G K G K G K

G K G K G K

G K G K G K

C^e

2 0 0

0

3 0 4 3

2 3

2

3 0 2 3

4 3

2

3 0 2 3

2 3

4

式(15)

式中

K：體積模數 G：剪力模數

FA 狀態之勁度矩陣

C

^c如式(16)所示，式(17)~式 (19)中的 Cij是指 RI 狀態的彈塑性矩陣

C

^i(ep)的第 i 行第j 列的值。

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

=

=

0 0 0 0

0 0 0

3 2 1

3 2 1

3 2 1

) (

C C C

C C C

C C C C

C^{c iep}

式(16)

3

31 21

1

C

11

C C

C = + +

式(17)

3

32 22

2

C

12

C C

C = + +

式(18)

3

33 23

3

C

13

C C

C = + +

式(19)

損傷函數 D（Disturbance function）隨著偏差塑 性應變ξ_D而改變[3]，如圖 3 所示。

圖3.偏差塑性應變-損傷函數之關係曲線(摘自 Varadarajan 與 Sharma, 2003 [2]) 損傷函數D 如下所示：

D D

B D A

ξ ξ

= +

式(20) 其中

A、B 為材料參數

A39-4

4 ξD：累積偏差塑性應變

( )

^1/2

∫

= _ij^p _ij^p

D dE dE

ξ

式(21)

ij ii p ij p

E

ij

ε ε δ 3

− 1

=

式(22) DSC 模式中的應力增量

d σ

^a可用HISS-δ0 的應 力增量

d σ

ⁱ與臨界狀態的應力增量

d σ

^c及本時階所 計算的損傷函數計算出來，即

( ) (

nⁱ

)

c n n c

n n i n n

a D d D d dD

dσ = 1− ⋅ σ ₊₁+ ⋅ σ ₊₁+ σ ₊₁−σ ₊₁ 式(23) 其中

n n

n

D D

dD =

₊₁

−

式(24) n：時階(timestep)

三、渥太華砂壓密排水試驗

本研究使用的砂樣為渥太華砂，為美國材料試驗 標準砂之一(ASTM C-109)，渥太華砂為均勻的石英 砂，其物理性質為Gs=2.66、D50=0.58mm、均勻係數

=1.38、曲率係數=0.97，屬均勻級配。渥太華砂的顆 粒形狀接近圓形，因此以其當實驗材料，以剔除顆粒 形狀之影響。本研究於取得渥太華砂烘乾後，砂土中 的雜質去除，以供試驗使用。

本研究中的砂土三軸試驗之主要目的在於求得 一個合理且通用的應力–應變模式以計算 DSC 的材料 參數。圍壓選擇69kPa、207kPa、345kPa，相對密度 Dr 分別為 20%、90%，在進行三軸試驗時，為避免砂 土於試驗過程中顆粒被壓碎，採用砂土皆不重覆使 用，資料的擷取及計讀由電腦自動記錄(圖 4)。

圖4. 三軸試驗儀器及計讀系統 本研究共進行6 組三軸試驗，包含了 3 組相對密 度為20%飽和砂、3 組相對密度為 90%飽和砂，各試 驗之結果與圖 5、6 所示。由這些試驗結果中發現鬆 砂的軸差應力隨著圍壓與軸向應變率增加而增加，當 軸差應力達到定值時會趨近穩定。

0 200 400 600 800 1000

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Deviator stress(kPa)

345kPa 207kPa 69kPa

圖5.不同圍壓下飽和砂相對密度 Dr=20%之三軸應力 -應變行為

-2.5 -1.5 -0.5 0.5 1.5 2.5

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Volume strain(%) 345kPa

207kPa 69kPa

圖6.不同圍壓下飽和砂相對密度 Dr=20%之三軸體積 -應變行為

由這些試驗結果中發現緊砂的軸差應力隨著圍 壓與軸向應變率增加而增加，當軸差應力達到定值時 會開始緩慢下降 (圖 7)，且試體之體積變化達到一個 定值後會開始上升如圖8 所示。

0 200 400 600 800 1000 1200 1400

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Deviator stress(kPa)

345kPa 207kPa 69kPa

圖7.不同圍壓下飽和砂相對密度 Dr=90%之三軸應力 -應變行為

-10 -5 0 5 10

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Volume strain(%)

345kPa 207kPa 69kPa

圖8.不同圍壓下飽和砂相對密度 Dr=90%之三軸體積 -應變行為

四、數值分析

本研究主要是以FLAC程式結合Fish語言所撰寫 之DSC模式副程式，代入本研究由三軸壓縮式驗所得 到DSC材料參數，再依此進行後續的研究。以下將分 別就本研究所採用的FLAC程式，DSC材料參數解決

IN GEOTECHNICAL ENGINEERING IN TAIWAN I-Lan

A39-5 方法分別說明。

FLAC 程式(Fast Lagrangian Analysis of Continua) 由美國Itasca Consulting Group, Inc. 發展以外顯式有 限差分法（Explicit finite difference method）為原則的 數值分析軟體，適用於深開挖穩定、隧道開挖、邊坡 滑動、基礎承載、土石壩分析等多項大地工程問題。

其運算過程是以「時階」方式(time-step)來求解每一 個節點的運動方程式，在每一個時階中，利用系統每 一節點的不平衡力，利用運動方程式，求出節點的新 速度和位移，然後決定節點之應變增量；再利用應力 應變組合律決定新的應力增量，藉此應力增量又可以 更新節點的不平衡力，如此程式將反覆執行此運算步 驟，直到各節點中最大的不平衡力趨近於零，達靜力 平衡為止。

本研究求DSC 參數之計算方法說明如下:

1. 由排水三軸試驗的實驗結果得到其基本參數，計 算出軸差應力-軸向應變及體積應變-軸向應變的關 係曲線並求其φ角等參數。

2. 利用式(25)與式(26)計算 P1 與 P2 再求出γ、β。

( )

[ ]

_⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛

= −

−

= ⁻

c c

P c

φ β φ

γ 3 sin

sin 3 1 2

1 ^1/4

式(25)

( )

[ ]

_⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛

= + +

= ⁻

E E

P E

φ β φ

γ 3 sin

sin 3 1 2

2 ^1/4

式(26)

其中

4

2 1 ⁻

⎟⎠

⎜ ⎞

⎝

=⎛

′ P

P P

式(27)

4

1 1 ⁻

⎟⎠

⎜ ⎞

⎝

⎛ + ′

− ′

= P

β P

式(28) 3. 計算應力不變量 J1、J2、J3以及S^r。

4. 將極限狀態的 J1、J2、J3及S^r以及步驟2 所得之γ、

β，代入降伏函數中得到 n 值。

5. 將 n、γ、β並令 F=0 求得各應力值所對應的α。

6. 將實驗所得之體積應變及軸向應變減去彈性應變 的部分得到塑性應變量再由其計算ξ以及ξD。 7. 由 D 與ξD的數值迴歸方法計算出A、B。

8. 由α、ξ的數值迴歸求出 a1、η1 值。

本研究所得的DSC 材料參數如表 1、2 所示:

表1. Dr=20%時之 DSC 所需參數

σ3(kPa) φ β γ n a1 η1 A B

345 32.7 0.6223 0.2 5.999 3.0E-7 0.612 0.418 -7E-15 307 33.4 0.6306 0.2 5.657 4.6E-6 0.633 0.428 -2E-14

69 34.9 0.6389 0.2 6.446 5.5E-5 0.719 0.401 1E-14

表2. Dr=90%時之 DSC 所需參數

然後利用FLAC 建立的 DSC 模式，將試驗數據 以及數值分析結果進行比對，以驗證分析模式的正確 性。

本試驗為三軸試驗，模擬的試驗材料為渥太華砂 其鬆砂單位重為 1700kg/m³、 緊砂單位重為 1810 kg/m³，本研究之數值模擬結果與三軸試驗結果的軸 差應力-軸向應變關係如圖7~圖 12 所示。

0 100 200 300 400 500 600 700 800 900

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Deviator stress(kPa)

實驗 DSC model

圖7. DSC 模式模擬曲線與三軸試驗結果比較 (Dr=20%,σ₃=345kPa)

0 100 200 300 400 500 600

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Deviator stress(kPa)

實驗 DSC model

圖8. DSC 模式模擬曲線與三軸試驗結果比較 (Dr=20%,σ₃=207kPa)

0 20 40 60 80 100 120 140 160 180 200

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Deviator stress(kPa)

實驗 DSC model

圖9. DSC 模式模擬曲線與三軸試驗結果比較 (Dr=20%,σ₃=69kPa)

0 200 400 600 800 1000 1200

0 5 10 15 20

Axial strain(%)

Deviator stress(kPa)

實驗 DSC model

圖10. DSC 模式模擬曲線與三軸試驗結果比較 (Dr=90%,σ₃=345kPa)

σ3(kPa) φ β γ n a1 η1 A B

345 38.1 0.6827 0.2 14.84 1.3E-21 2.309 0.353 3E-14 307 38.6 0.6883 0.2 11.94 2.8E-14 1.207 0.356 -4E-14

69 34.4 0.7149 0.2 12.95 1.9E-10 1.051 0.359 -3E-14

A39-6

6

0 100 200 300 400 500 600 700 800

0 2 4 6 8 10 12 14 16 18

Axial strain(%)

Deviator stress(kPa)

實驗 DSC model

圖11. DSC 模式模擬曲線與三軸試驗結果比較 (Dr=90%,σ₃=207kPa)

0 50 100 150 200 250 300

0 2 4 6 8 10 12 14 16 18

Αxial strain(%)

Deviator stress(kPa)

實驗 DSC model

圖12. DSC 模式模擬曲線與三軸試驗結果比較 (Dr=90%,σ₃=69kPa)

由分析結果顯示，本研究所採用的 DSC 參數計 算方法所求得的參數計算結果再套入 FLAC 數值模 擬後，其結果與三軸試驗結果大致吻合，可以得知本 研究所完成DSC 模式 Fish 副程式與 HISS-δ0 模式參 數計算結果有其正確性。

五、結論

本研究以 DSC 模式為主軸，配合以有限差分為 基礎之數值分析軟體 FLAC 程式來探討土壤之應變 軟化行為。以三軸試驗結果反算 DSC 材料參數最佳 值，並模擬三軸試驗。經由上述分析結果，歸納以下 之結論：

1. 在本研究的討論範圍中 Desai 發展的 DSC 模式可 模擬土壤之應變軟化的力學行為。

2. 本研究由三軸試驗結果推算所得的擾動狀態概 念「DSC 模式」材料參數，此計算方法所得的參 數有其正確性，因為後續數值模擬結果與試驗結 果大致相符。

參考文獻

[1]. Desai, Chandra S., “Mechanics of Materials and Interface the Disturbed State Concept”, CRC Press, Boca Raton, Florida (2000).

[2]. A. Varadarajan, K. G. Sharma, K. Venkatachalam, A.

K. Gupta.,“Testing and Modeling Two Rockfill Materials,” International Journal of Geomechanics, Vol. 129, No. 3 (2003).

[3]. A. Varadarajan, K. G. Sharma, S. M. Abbas and A.

K. Dhawan., “Constitutive Model for Rockfill Materials and Determination of Material Constants,” Vol. 6, No. 4 (2006).

[4]. Chandra S. Desai and Joseph Y. Chen.,“Parameter Optimization and Sensitivity Analysis for Disturbed

State Constitutive Model,” International Journal of Geomechanics, Vol. 6, No. 2 (2006).

[5]. 黃秀雅，「利用擾動概念評估損傷構件之承載行 為」，朝陽科技大學營建工程系碩士論文，2008。

[6]. 陳志瑋，「粉土質砂土中摩擦型群錨之垂直受力 行為」，朝陽科技大學營建工程系碩士論文，

2008。