一個壓縮試驗所引發的“血案”

1. 背景

前段時間在和一個國外學生Adem交流一組試驗數據的處理方法時出現了以下的一個分歧：試驗目的主要是為了測試等向壓縮過程中不同基質吸力狀態下非飽和膨脹土的體變模量(E)問題，最終將這一參數用于計算淺層地基的變形，所以對應的荷載其實并不大。其中一組等吸力狀態下的體變-凈平均應力的關系曲線試驗數據如圖1所示，其他試驗數據與此類似。

圖1 原始試驗數據

方法一：

從圖中可以發現，由于試驗中所施加的凈平均應力并不大，最大也僅僅為200kPa，另一方面試樣的干密度又較大，加之又處于非飽和狀態，導致先期固結壓力較大，所以在整個壓縮過程中曲線并沒有顯現出明顯的非線性，使用線性擬合發現擬合效果也相當不錯，R²=0.9852.

如果根據體變模量的定義（D. G. Fredlund et al. 1993）將幾個變形參數公式寫在一起；

；(1)

不難發現，體變模量可以直接根據應力應變曲線每一個加載步的斜率進行計算，即E = 3*(1-2μ) / (slope in each load step)，但是這時存在一個問題：由于應力應變曲線并不嚴格為直線，那么在一個常吸力狀態下的等向壓縮過程中每一個加載步都會對應得到一個體變模量值，這樣的結果無論是在之后的模型建立中還是在實際的數值模擬使用中都是極為不便的。所以做出了如下的簡化處理，近似線性的應力應變曲線則可以直接使用整個直線的斜率來代替每個加載步的斜率，結果就可以在一個常吸力下對應一個體變模量值，這樣也符合我們實驗的最終的預期成果之一：分析吸力對土體變形模量的影響，進而探討與非飽和膨脹土體變本構模型應用方面的相關問題，而暫時忽略了不同應力水平下體變模量的差異問題。

如果假設泊松比為0.33，那么很容易得到變形模量為1902.99kPa，使用這個變形模量值來進行應力應變關系的計算時，由于是逆運算的關系，得到的自然是圖1中的線性擬合線。

方法二：

但是Adem堅持使用增量法來進行計算，用到了式(1)中的最后一項，如果將0-200kPa作為一個大的加載步，那么(σ-u_a)_ave=100kPa，壓縮指數C_t取e-lg(σ-u_a)曲線的直線段值（圖2所示），使用式 (2) 的最后一項來進行計算。

；(2)

圖2 原始數據e - lgp曲線

那么得到的結果為：

C_t = ln10*(the slope of linear part of e-lg(p) curve)

E = 6.908*(1-2μ)(1+e₀)/C_t*(σ-u_a)_ave = 406.1327 kPa.

發現兩種方法計算的結果相差巨大，已經不僅僅是誤差范圍內的差異了，從原理上來講，應力應變曲線和e-lg(σ-u_a)曲線都是描述土體變形的一種方法，僅僅狀態變量是選用體應變還是孔隙比之間的差別，應該不會有如此大的差別了。

那么問題究竟出在哪個地方呢？

在于對數坐標的使用上！

2. 對數坐標的來源

在國內的幾個土力學教材版本中，都是在土的壓縮性部分介紹出了半對數坐標系e-lg(p)曲線的使用，但是都沒有解釋為什么采用半對數坐標系。此外，在粒徑分布曲線和非飽和土中的SWCC等曲線中也都使用了對數坐標系。

回顧這些應用，不免會思考下面這個問題：其實更加簡單的雙曲線也在通常的坐標系下也會有類似的形狀，為什么單單會有如此多的地方使用對數坐標系？會不會有更深層次的原因呢？

其實，更重要的是這些方面的應用都有一些共同的特征：

(1) 變量之一在所研究的范圍內發生了幾個數量級的變化；

(2) 在自變量由零開始逐漸增大的初始階段，當自變量的少許變化引起因變量極大變化時，此時采用半對數坐標，曲線最大變化范圍可伸長，使圖形輪廓清楚；

(3) 需要將某種函數變換為直線函數關系。

回顧對數的產生過程：

盡管對數的產生要比指數更早，前者產生于1614年納皮爾出版的《奇妙的對數定律說明書》（Joost Bürgi獨立發現了對數，但直到納皮爾之后4年才發表），之后布里格斯于1624年在出版的《對數算術》中對其加以改造，并使之廣泛流傳。而指數符號則直到近二十年后的1637年才有法國數學家笛卡爾首次采用，只到1770年，歐拉在出版的一部著作中使用了y=a^x 來定義x=log_a(y)，他指出：“對數源于指數”。(引自：百度百科)

那么在巖土工程里面為何要使用對數呢？

那得從土的壓硬性說起。在側限壓縮中，并不會出現剪切中的破壞問題，隨著荷載的增大，土體不斷被壓縮，密度不斷增大，變得越來越“硬”。這一現象很容易聯想到我們常見的復利計算方法，隨著時間的增加，利息不斷加入原有本金而用于計算下一個計息周期的利息，在土體中，一定空間內土顆粒的不斷增多用于抵抗下一個階段的荷載。Oh yeah，看，是不是有點跟上節奏了。最終本息和將成為計息周期的指數函數：P=A*(1+i)^n，而土體所能承擔的荷載是土體密度（注：由于土體的最終變形是處于荷載全部由土體骨架所承擔，所以此處代表干密度）的指數函數，如果使用孔隙比變化來表征土體的干密度變化，就有了p=1.0E[N-(1+e)/λ]，底數也可以換成自然對數等其他base形式，進行變換就得到了常用的NCL表達式：ν=1+e=N-λ*lg(p)。當然，使用自然對數還是10為底的對數區別僅僅在于曲線的斜率之間會有一個常數倍的關系，但是他們的起源是不同的，筆者將在下一篇文章里面更加詳細的討論土的壓硬性和復利計算類比這一問題。

總結一下，我們應該在下列情況下應用對數坐標系：

(1) 如果所研究的函數y和自變量x在數值上變化了幾個數量級；

(2) 需要將曲線開始部分劃分成展開的形式；

(3) 當需要變換某種非線性關系為線性關系時。

這時，一個同學勇敢的舉手發言：難道說我們學習了N年的土力學教材就被你這隨便一YY就失效了? 嚴謹一點，好不。。。Too young, too simple !

好好好，這位同學，坐下喝杯茶，先壓壓驚，且聽慢慢道來。

在《土力學》教材中通常講解了四個與側限壓縮相關的參數：壓縮系數、壓縮指數、壓縮模量和體積壓縮系數。

壓縮系數：孔隙比增量和荷載增量的比值，a=-△e/△p；

壓縮指數：孔隙比增量和荷載對數增量的比值，C_c=-△e/△(lgp)；

側限壓縮模量：應力增量和應變增量的比值，E_s=△σ_z/△ε_z=(1+e₀)/a；

體積壓縮系數：側限壓縮模量的導數，m_v=1/E_s=a/(1+e₀)；

從計算公式來看，都從增量的角度進行定義。那么對于實際的使用來講，由于大多的土體的壓縮曲線是非線性的，如果使用“壓縮系數”、“壓縮模量”和“體積壓縮系數”的話，勢必會造成在每一個應力狀態下這些參數都有一個不同的參數數值的情況，而且參數隨著荷載的變化還會比較明顯，再加上什么含水率、什么應力比、應力水平、加載速率等等因素，是不是已經暈了，搞得這么復雜，這不是要讓本來就苦逼的工程師們罵街的節奏么! So, 半對數坐標系應運而生，除了壓力比較小的范圍有些乏力以外，一個λ或者C_c就可以描述相當范圍內比體積或者孔隙比隨荷載的變化規律，簡直是爽爆了有木有。而且，通常情況下自重應力再加上基礎傳來的附加應力導致地基土體中的應力都大于100kPa，在大于100kPa范圍內的壓縮曲線使用e-lg(p)來表示是沒有問題的，見圖2。這也是在《土力學》教材(李廣信，2013)和《建筑地基基礎設計規范》(GB 50007 - 2011)中對土體的壓縮性(低、中、高)分類是以100-200kPa荷載段來進行的原因之一。而且，在《建筑地基基礎設計規范》(GB 50007 - 2011)和《巖土工程勘察規范》(GB 50021 - 2001)(2009版)中也都明確指出壓縮系數和壓縮模量會隨著荷載而改變，在使用各向同性均質線性變形體理論進行地基變形計算時也是使用“土的自重壓力至土的自重壓力和與附加壓力之和的壓力段”所對應的壓縮模量，而非直接使用λ或C_c來進行變形計算。

所以說在規范所考慮的通常情況下，一般是不會出現本文中所涉及的問題的。

但是，如果不是規范中所默認的通常情況，而是需要計算幾個kPa、幾十個kPa情況下的基礎變形問題怎么辦?比如說月球上的基礎設計呢?或者超固結比很大的超固結土，其應力狀態處于線彈性范圍內的情況呢(本文的圖1)?這時就必須考慮在壓力比較小的情況下半對數坐標系的適用性問題了。

寫到這里發現本文其實也沒有解決任何實際的問題，僅供思考，謝謝觀看！:)

【參考文獻】

[1]   Fredlund, D. G. & Rahardjo, H. 1993. Soil mechanics for un-saturated soils[M]. New York: John Wiley & Sons Inc.

[2]   李廣信, 張丙印, 余玉貞. 2013. 土力學[M]. 北京: 清華大學出版社.

[3]   中國建筑科學研究院. 建筑地基基礎設計規范[S]. 北京: 建筑工業出版社.