馬 勇， 王 飛， 杜秀忠

(廣東省水利水電科學(xué)研究院，廣東省巖土工程技術(shù)研究中心，廣東 廣州 510635)

基于MATLAB的土體三軸試驗(yàn)抗剪強(qiáng)度指標(biāo)求解計(jì)算

馬 勇， 王 飛， 杜秀忠

(廣東省水利水電科學(xué)研究院，廣東省巖土工程技術(shù)研究中心，廣東 廣州 510635)

莫爾應(yīng)力圓通常是利用作圖的方法獲得，進(jìn)而求得巖土體的強(qiáng)度指標(biāo)，但該方法易受人為影響，精度較低。通過MATLAB軟件編制相應(yīng)的計(jì)算程序求得巖土體的強(qiáng)度指標(biāo)，該軟件是基于最小二乘法原理，對最大主應(yīng)力和最小主應(yīng)力的線性回歸方程進(jìn)行求值，進(jìn)而通過抗剪強(qiáng)度指標(biāo)與最大主應(yīng)力和最小主應(yīng)力的關(guān)系反算求得c、φ值。通過相關(guān)文獻(xiàn)數(shù)據(jù)，驗(yàn)證本文方法的準(zhǔn)確性和可靠性，為巖土體抗剪強(qiáng)度指標(biāo)值的計(jì)算提供了一條新的解決途徑。

三軸試驗(yàn)；抗剪強(qiáng)度指標(biāo)；MATLAB；最小二乘法

計(jì)算巖土體抗剪強(qiáng)度有許多種方法[1-2]，三軸壓縮試驗(yàn)是測定巖土體抗剪強(qiáng)度的一種常用的試驗(yàn)方法，試驗(yàn)時(shí)，通常利用多個(gè)一定規(guī)格的圓柱形試驗(yàn)樣品，在保持一定的圍壓下分別進(jìn)行軸向壓力的施作，對試樣進(jìn)行剪切試驗(yàn)直至試樣破壞，獲得試樣在不同圍壓狀態(tài)下的極限軸向壓力，不同的軸向壓力和相對應(yīng)的圍壓也分別稱為大、小主應(yīng)力。然后，對獲得的數(shù)據(jù)根據(jù)強(qiáng)度準(zhǔn)則[3]進(jìn)行擬合處理，得出試樣的破壞包絡(luò)線，進(jìn)行求得試樣的巖土體強(qiáng)度參數(shù)。試驗(yàn)數(shù)據(jù)處理時(shí)，由于計(jì)算繪圖工作量大，需要做大量的重復(fù)性工作，且作圖法亦會引入一些人為誤差，隨意性較大。為了得到較為準(zhǔn)確的抗剪強(qiáng)度值，目前抗剪強(qiáng)度指標(biāo)的求解主要是在試驗(yàn)數(shù)據(jù)的基礎(chǔ)上采用較嚴(yán)謹(jǐn)?shù)膹?qiáng)度理論進(jìn)行推導(dǎo)或利用軟件分析，以減少手工作圖產(chǎn)生的隨意誤差[4]。阮波等[5]基于非線性規(guī)劃方法，利用EXCEL軟件對三軸剪切試驗(yàn)的破壞包絡(luò)線直接進(jìn)行規(guī)劃求解，進(jìn)而獲取了試樣的抗剪強(qiáng)度參數(shù)值。李小磊等[6]通過軟件對Matrix VB的功能進(jìn)行調(diào)用，利用合理的程序設(shè)計(jì)，基于最小二乘法原理對摩爾應(yīng)力圓的切點(diǎn)進(jìn)行線性回歸，進(jìn)而求解得到抗剪強(qiáng)度指標(biāo)值。王星輝[7]基于嚴(yán)格的理論推理，采用數(shù)據(jù)直接計(jì)算的方法來替代以前手工作圖或計(jì)算機(jī)交互圖解，取得了較好的效果。三軸剪切試驗(yàn)的破壞包絡(luò)線也可以采用CAD作圖的方法來繪制，但其繪圖的效率和精度均不理想，而且也無法直接獲得強(qiáng)度參數(shù)值，更由于個(gè)人的經(jīng)驗(yàn)和誤差的原因，得到的強(qiáng)度參數(shù)值極易出現(xiàn)多解的現(xiàn)象。

本文在Mohr-Coulomb強(qiáng)度準(zhǔn)則的假設(shè)下，利用最小二乘法的線性回歸原理，通過MATLAB軟件編制相應(yīng)的計(jì)算程序，利用σ1-σ3回歸方法進(jìn)行抗剪強(qiáng)度的求解，為三軸試驗(yàn)求解抗剪強(qiáng)度指標(biāo)提供了一條新思路。

1 數(shù)學(xué)模型

1.1 σ1-σ3法[8]

在對巖土體三軸試驗(yàn)數(shù)據(jù)進(jìn)行分析時(shí)，假設(shè)巖土的剪切破壞強(qiáng)度滿足M-C準(zhǔn)則：

τ=σtanφ+C

(1)

圖1 極限平衡狀態(tài)的應(yīng)力圓

利用破壞應(yīng)力圓知(見圖1)：

(2)

整理后得

σ1=A+Bσ3

(3)

其中，

(4)

巖土三軸試驗(yàn)可以直接測得試樣的大小主應(yīng)力σ1i和σ3i，然后根據(jù)線性回歸原理直接對(σ1i，σ3i)數(shù)據(jù)進(jìn)行最小二乘法線性回歸，得到A、B值，再利用式(4)通過反算就可以計(jì)算出抗剪強(qiáng)度c、φ值。

1.2 最小二乘法的一元線性回歸

在工程實(shí)踐和科學(xué)實(shí)驗(yàn)中，需要把實(shí)驗(yàn)數(shù)據(jù)擬合成經(jīng)驗(yàn)公式，以反映變量之間的關(guān)系，這稱為回歸分析[9]。對于樣本數(shù)據(jù)(xi，yi)，若存在yi=a+bxi+εi，并假定：

1) 誤差εi是數(shù)學(xué)期望值為0、方差不變的隨機(jī)變量，即E(εi)=0，Var(εi)不變；

2)xi為非隨機(jī)變量。

則使

∑(yi-a-bxi)2=min

(5)

滿足式(5)的a，b值是回歸方程y=a+bx的回歸系數(shù)的最佳線性無偏估計(jì)量。

a，b的值可以通過下式求得

(6)

2 MATLAB算法實(shí)現(xiàn)

MATLAB[10]軟件的功能很強(qiáng)大，特別是在數(shù)值計(jì)算和可視化處理方面。它強(qiáng)大的內(nèi)核中集成了海量的數(shù)值計(jì)算函數(shù)，而且覆蓋面很全面，用戶可以將所需解決的問題概化成數(shù)學(xué)表達(dá)式，計(jì)算結(jié)果就可以通過MATLAB強(qiáng)大的計(jì)算和可視化功能完美地展現(xiàn)出來。MATLAB具有函數(shù)和圖形窗口兩種進(jìn)行曲線擬合的方法，考慮到本文擬合的強(qiáng)度包絡(luò)線為線性函數(shù)，為了對不同的三軸試驗(yàn)數(shù)據(jù)進(jìn)行多次重復(fù)計(jì)算應(yīng)用，現(xiàn)采用第一種方法通過建立M函數(shù)對三軸試驗(yàn)數(shù)據(jù)進(jìn)行擬合。

2.1 多項(xiàng)式擬合函數(shù)polyfit

在MATLAB中，進(jìn)行多項(xiàng)式曲線擬合函數(shù)的調(diào)用格式為：

p=polyfit(x，y，n)

式中x，y采用矩陣形式，表示參與曲線擬合的實(shí)驗(yàn)數(shù)據(jù)；n為擬合多項(xiàng)式的次數(shù)，本文中為線性擬合，故n=1。

2.2 建立M函數(shù)

在MATLAB中，一個(gè)M函數(shù)就相當(dāng)于一個(gè)獨(dú)立的程序，本文在進(jìn)行程序設(shè)計(jì)中，將最大主應(yīng)力σ1和最小主應(yīng)力σ3分別以矩陣的形式作為數(shù)據(jù)輸入變量。只要按照一定的格式輸入三軸壓縮試驗(yàn)中的最大主應(yīng)力和最小主應(yīng)力的值，就可以得出抗剪強(qiáng)度指標(biāo)c、φ的值，并可以顯示強(qiáng)度包絡(luò)線和莫爾應(yīng)力圓的圖形。

程序代碼如下：

function moco (a3，a1) %a3為最小主應(yīng)力，即圍壓，a1為最大主應(yīng)力;

syms c fai;%定義粘聚力和內(nèi)摩擦角

p=polyfit(a3，a1，1); %一次線性擬合

fai=(asin((p(1)-1)/(p(1)+1)))/pi*180; %求內(nèi)摩擦角

c=p(2)*(1-sin(fai/180*pi))/(2*cos(fai/180*pi)); %求粘聚力

c=vpa(c，4) %控制計(jì)算精度

fai=vpa(fai，4) %控制計(jì)算精度

[m，n]=size(a3); %獲取三軸試驗(yàn)組數(shù)

for i=1∶1∶n %畫出所有的應(yīng)力圓

alpha=0∶pi/50∶pi; %半圓角度[0，pi]

R=(a1(i)-a3(i))/2; %摩爾圓的半徑

x=R*cos(alpha)+(a1(i)+a3(i))/2; %x的坐標(biāo)

y=R*sin(alpha); %y的坐標(biāo)

plot(x，y，‘-’) %畫摩爾圓

hold on %多張圖畫在同一個(gè)坐標(biāo)上

axis equal

end

x=0∶1∶1.2*a1(n);

y=x*tan(fai/180*pi)+c;

plot(x，y，‘r-’) %畫強(qiáng)度包絡(luò)線

xlabel(‘主應(yīng)力σ/KPa’，‘fontsize’，14，‘fontname’，

‘Times new roman’);%x坐標(biāo)控制

ylabel(‘剪應(yīng)力τ/KPa’，‘fontsize’，14，‘fontname’，

‘Times new roman’);%y坐標(biāo)控制

axis([0 1.2*a1(n) 0 0.7*a1(n)]); %坐標(biāo)范圍控制

title(‘應(yīng)力圓和強(qiáng)度包絡(luò)線’，‘fontsize’，14);%標(biāo)題

axis on;axis xy;

3 算法驗(yàn)證與分析

為了驗(yàn)證MATLAB編程的正確性，本文參考選取了多組試驗(yàn)數(shù)據(jù)進(jìn)行對比分析。算例1：參考文獻(xiàn)[6]中的固結(jié)不排水剪切試驗(yàn)數(shù)據(jù)；算例2：參考文獻(xiàn)[11]中的泥質(zhì)粉砂巖三軸抗剪試驗(yàn)數(shù)據(jù)；算例3：參考文獻(xiàn)[12]中的不固結(jié)不排水三軸剪切試驗(yàn)數(shù)據(jù)；算例4：參考文獻(xiàn)[13]中砂巖的三軸試驗(yàn)數(shù)據(jù)。

根據(jù)不同參考文獻(xiàn)提供的試驗(yàn)數(shù)據(jù)(即巖土體達(dá)到破壞狀態(tài)下不同圍壓對應(yīng)的極限軸向壓應(yīng)力)，參考文獻(xiàn)中利用常規(guī)方法獲得了每一組試驗(yàn)數(shù)據(jù)下的巖土體強(qiáng)度參數(shù)值，本文在此基礎(chǔ)上利用筆者所編的程序進(jìn)行重新計(jì)算，并與參考文獻(xiàn)中得到的抗剪強(qiáng)度參數(shù)進(jìn)行對比，見表1～4所示。

表1 某固結(jié)不排水剪切試驗(yàn)結(jié)果對比

表2 某泥質(zhì)粉砂巖剪切試驗(yàn)結(jié)果對比

表3 某土樣不固結(jié)不排水剪切試驗(yàn)結(jié)果對比

表4 某砂巖三軸剪切試驗(yàn)結(jié)果對比

從表1～4的數(shù)據(jù)中可以看出:各參考文獻(xiàn)中提供的抗剪強(qiáng)度參數(shù)值與本文的計(jì)算值基本一致;φ值的整數(shù)部分基本完全一樣，只是小數(shù)部分有些差異，c值差別相對大一些，這可能是由于直線的截距隨著直線與橫坐標(biāo)夾角的微小變化而被放大的結(jié)果。不同參考文獻(xiàn)對比結(jié)果的誤差率詳見表5。

表5 本文計(jì)算值與參考文獻(xiàn)結(jié)果的誤差率

由表5可以看出：本文的計(jì)算方法與參考文獻(xiàn)實(shí)驗(yàn)室結(jié)果的平均誤差率c值為1.97%、φ值為0.28%，表明本文的方法是可靠的，而且本文的計(jì)算成果在精度上有了很大的提高，大大提高了數(shù)據(jù)處理效率，可以應(yīng)用于工程實(shí)踐中。

利用表1中編號583土樣的數(shù)據(jù)可以繪制該土樣的應(yīng)力圖及抗剪強(qiáng)度包絡(luò)線，見圖2。

圖2 應(yīng)力圓和抗剪強(qiáng)度包絡(luò)線

4 結(jié)語

本文利用MATLAB數(shù)值計(jì)算軟件，基于巖土體的抗剪強(qiáng)度準(zhǔn)則，編制相應(yīng)的計(jì)算程序處理三軸試驗(yàn)數(shù)據(jù)，計(jì)算方法比較合理，精度和效率較高，計(jì)算結(jié)果比較合理且具有唯一性，較傳統(tǒng)繪圖方法優(yōu)勢比較明顯，可供土工試驗(yàn)人員等參考。

[1] 孫敬.不同方法測定軟土參數(shù)的比較分析[J]. 廣東水利水電，2012(4)：20-22．

[2] 朱思軍， 楊光華， 陳富強(qiáng)，等. 十字板剪切試驗(yàn)在珠三角深厚軟土基坑工程中的應(yīng)用[J]. 廣東水利水電， 2016(6)：28-33．

[3] 南京水利科學(xué)研究院.土工試驗(yàn)規(guī)程：SL237—1999[S]. 北京:中國水利水電出版社，1999．

[4] 江穩(wěn)瓊. 淺談固結(jié)、三軸試驗(yàn)中存在的一些問題及改進(jìn)措施[J].廣東水利水電，2009(5): 29-33．

[5] 阮波， 張向京， 彭意. Excel規(guī)劃求解三軸試驗(yàn)抗剪強(qiáng)度指標(biāo)[J]. 鐵道科學(xué)與工程學(xué)報(bào)，2009，6(5)：57-60．

[6] 李小磊， 吳云剛.求解c、φ值的軟件開發(fā)及試驗(yàn)驗(yàn)證[J]. 安全與環(huán)境工程，2010，17(1)：103-106．

[7] 王星輝. 三軸剪切總應(yīng)力強(qiáng)度指標(biāo)的回歸求值[J]. 地球科學(xué)與環(huán)境學(xué)報(bào)，2004，26(1)：52-54．

[8] 陳立宏，陳祖煜，李廣信.三軸試驗(yàn)抗剪強(qiáng)度指標(biāo)線性回歸方法的討論[J]. 巖土力學(xué)，2005，26(11)：1 785-1 789．

[9] 鄭咸義，姚仰新，雷秀仁，等，應(yīng)用數(shù)值分析[M]. 廣州：華南理工大學(xué)出版社，2008．

[10] 劉正君. MATLAB科學(xué)計(jì)算寶典[M]. 北京：電子工業(yè)出版社， 2012．

[11] 李雷， 陳尚星. 巖石三軸抗剪強(qiáng)度的計(jì)算[J].華北水利水電學(xué)院學(xué)報(bào)， 2002， 23(4):46-48．

[12] 張書憲， 劉春平. 求解三軸抗剪強(qiáng)度指標(biāo)的一種方法[J]. 勘察科學(xué)技術(shù)， 1993(3):15-18．

[13] 楊同， 徐川. 巖土三軸試驗(yàn)中的粘聚力與內(nèi)摩擦角[J]. 中國礦業(yè)， 2007， 16(12):104-107．

(本文責(zé)任編輯 馬克俊)

The Solution Computation of Triaxial Tests Shear Strength Index Based On Matlab

MA Yong, WANG Fei, DU Xiuzhong

(Guangdong Research Institute of Water Resources and Hydropower, The Geotechnical EngineeringTechnology Center of Guangdong Province, Guangzhou 510635，China)

Mohr’s circle is widely used in soil mechanics; it is generally finite precision and affected by human factors due to be obtained by mapping method. This paper can obtain shear strength index of triaxial test by using Matlab software compiled corresponding program. The procedure can calculate the value of c、φ through the relation between the hear strength index and principal stress by using the principle of least square method. The accuracy and reliability of this method are verified by the relevant literature data, which provides a new way to calculate the shear strength of soil.

triaxial tests; shear strength index; matlab; least square method

2016-11-06;

2016-11-29

馬勇(1988)，男，碩士，主要從事巖土工程設(shè)計(jì)咨詢工作。

TU43

：B

：1008-0112(2016)012-0018-04