基于MATLAB的重力式擋土墻可靠度計(jì)算程序開(kāi)發(fā)及其應(yīng)用

宋云連 付 宇 李 強(qiáng),2

1(內(nèi)蒙古工業(yè)大學(xué)內(nèi)蒙古自治區(qū)土木工程結(jié)構(gòu)與力學(xué)重點(diǎn)實(shí)驗(yàn)室 內(nèi)蒙古 呼和浩特 010051) 2(鄭州市市政工程總公司 河南 鄭州 450001)

0 引 言

重力式擋土墻因其結(jié)構(gòu)簡(jiǎn)單、施工簡(jiǎn)便、可就地取材等特點(diǎn)被廣泛應(yīng)用。當(dāng)前的擋土墻穩(wěn)定性分析和驗(yàn)算軟件以安全系數(shù)法為主流，盡管該方法是一種相對(duì)有效的工程實(shí)用法，因其未考慮參數(shù)的不確定性，在實(shí)際工程中會(huì)出現(xiàn)安全系數(shù)的大小和工程結(jié)構(gòu)的安全程度不匹配的現(xiàn)象[1]。為此，將以概率論為基礎(chǔ)的可靠度理論引入到擋土墻穩(wěn)定性分析，國(guó)內(nèi)外已有一定的研究成果。杜永峰等[2]以一重力式擋土墻為例，采用JC法對(duì)擋土墻的兩種主要失效模式進(jìn)行計(jì)算，并用串聯(lián)理論進(jìn)行系統(tǒng)靠度分析。王海祥等[3]采用蒙特卡羅法對(duì)重力式擋土墻的抗傾覆和抗滑移穩(wěn)定性進(jìn)行分析計(jì)算并與定值分析法進(jìn)行比較，得出某些擋土墻安全系數(shù)足夠，而投入使用卻發(fā)生破壞的結(jié)論。Low[4]采用一次二階矩法在考慮擋土墻參數(shù)正態(tài)和非正態(tài)分布基礎(chǔ)上對(duì)擋土墻的計(jì)算進(jìn)行程序設(shè)計(jì)，并將其與Excel相結(jié)合使計(jì)算過(guò)程更加簡(jiǎn)化明了。Castillo等[5]將傳統(tǒng)安全系數(shù)法與概率理論方法相結(jié)合計(jì)算出重力式擋土墻的失效概率，并對(duì)其進(jìn)行敏感性分析和計(jì)算過(guò)程的優(yōu)化設(shè)計(jì)。Sivakumar等[6]采用可靠度理論對(duì)影響懸臂式擋土墻穩(wěn)定性的不確定參數(shù)隨機(jī)模擬，確定了抗滑移、抗傾覆等可靠指標(biāo)和失效概率的取值范圍。

綜上所述，以定值理論為基礎(chǔ)的擋土墻分析軟件已經(jīng)不能滿足實(shí)際工程需求；將可靠度引入擋土墻穩(wěn)定性分析恰好可以彌補(bǔ)定值計(jì)算的缺陷，且該方法已經(jīng)具備一定的理論支撐。此外，因擋土墻可靠度分析涉及復(fù)雜方程且運(yùn)算量極大，需輔以計(jì)算機(jī)。查閱計(jì)算機(jī)軟件類文獻(xiàn)發(fā)現(xiàn)，MATLAB的高效算法及其完善的GUI窗口開(kāi)發(fā)功能與該法的需求不謀而合[7-12]。因此，本文以蒙特卡羅法為基礎(chǔ)，結(jié)合MATLAB的GUI功能，開(kāi)發(fā)出一款面向?qū)ο蟮闹亓κ綋跬翂煽慷扔?jì)算和敏感性分析的通用計(jì)算程序軟件。該軟件具備在重力式擋土墻尺寸、墻后邊坡參數(shù)及墻后填土隨機(jī)參數(shù)已知的條件下，可直接計(jì)算出重力式擋土墻抗傾覆和抗滑移的失效概率與可靠指標(biāo)，并繪制隨機(jī)變量參數(shù)的敏感性分析曲線圖。此外，該軟件也具備定值計(jì)算的功能，可以求解抗傾覆和抗滑移安全系數(shù)。

1 重力式擋土墻模型及其功能函數(shù)的建立

1.1 擋土墻通用計(jì)算模型建立

對(duì)擋土墻的穩(wěn)定性進(jìn)行分析計(jì)算的前提是建立擋墻通用模型，如圖1所示。

圖1 重力式擋土墻通用模型

該模型的主要參數(shù)有：擋土墻自身尺寸參數(shù)(高度H、頂寬dk、底寬Dk、墻背傾角α、墻趾高度hj、墻趾寬度bj、墻底傾角a0)；填料參數(shù)(填料容重γ、均值γmu、標(biāo)準(zhǔn)差γsig；填土內(nèi)摩擦角均值φmu、標(biāo)準(zhǔn)差φsig；填土與墻背間內(nèi)摩擦角δ、均值δmu、標(biāo)準(zhǔn)差δsig)；墻后邊坡參數(shù)(墻頂與路面高程差a、內(nèi)邊坡斜率α*、路肩寬b、行車(chē)道寬度L0、外邊坡斜率α1)。

1.2 擋土墻自身截面面積及重心計(jì)算

為簡(jiǎn)化擋土墻截面面積及其重心力矩的公式推導(dǎo)，現(xiàn)將圖1中的擋土墻分為固定的四個(gè)三角形模塊，經(jīng)過(guò)推導(dǎo)每一模塊的面積計(jì)算公式如下：

S1=0.5×Dk×hj

(1)

S2=0.5×[(hj+Dk×tana0)×tana+Dk]×

(hj+Dk×tana0)

(2)

S3=0.5×[(hj+Dk×tana0)×tana+Dk-bj]×

(H-hj-Dk×tana0)

(3)

S4=0.5×dk(H-hj-Dk×tana0)

(4)

經(jīng)推導(dǎo)各模塊重心距離墻趾的距離為：

l1=Dk/3

(5)

l2=[(hj+Dk×tana0)×tana+2×Dk]/3

(6)

l3=[(hj+Dk×tana0)×tana+Dk+2×bj+

(H-hj-Dk×tana0)×tana]/3

(7)

l4=[(hj+Dk×tana0)×tana+2×

(H-hj-Dk×tana0)×tana+Dk+dk]/3

(8)

式中：墻背俯斜時(shí)α為負(fù)值，仰斜時(shí)α為正值。經(jīng)計(jì)算可得重力式擋土墻的截面面積為：

S=S1+S2+S3+S4

(9)

其自重重心距墻趾的距離為：

l=(l1+l2+l3+l4)/2

(10)

1.3 最大主動(dòng)土壓力及其作用點(diǎn)計(jì)算

在程序設(shè)計(jì)過(guò)程中，擋土墻墻后主動(dòng)土壓力計(jì)算采用庫(kù)侖土壓力理論。假設(shè)墻后填土破裂面為直線型，取單位長(zhǎng)度的擋土墻墻后破裂棱體，利用力多邊形和正弦定理求得墻后主動(dòng)土壓力Ea的表達(dá)式如下[13]：

(11)

式中：Θ=a+δ+φ，α和δ角度的含義見(jiàn)圖1；G=γS為破裂棱體單位長(zhǎng)度上的自重。破裂棱體橫斷面的面積可以表示為：

S=A0tanθ-B0

(12)

式中：A0和B0是用來(lái)表示不同破裂棱體橫斷面的參數(shù)。由于破裂角θ交與位置不同時(shí)(如交于內(nèi)邊坡，路面上荷載內(nèi)部、中部及外部)，只有單位長(zhǎng)破裂棱體自重G的表達(dá)式不同，當(dāng)破裂棱體采用同一種土質(zhì)修筑時(shí)，G的不同之處表現(xiàn)在破裂棱體橫斷面面積S的不同。利用軟件計(jì)算S時(shí)，其破裂面的搜索范圍為0到90度，搜索步長(zhǎng)為0.01度。根據(jù)破裂面交于位置不同，各部位S中的主要參數(shù)及其作用點(diǎn)的計(jì)算公式如表1，各邊界條件的判定則是依據(jù)圖1中尺寸參數(shù)求出臨界破裂角θ的正切值。

表1 不同邊界條件對(duì)應(yīng)主動(dòng)土壓力主要參數(shù)計(jì)算公式

表1中，h1=(b-a×tanθ)/(tanθ+tanα)，h2=d/(tanθ+tanα)，h3=L0/(tanθ+tanα)，h4=H-h1-h2-h3。其中，Zy為各部位對(duì)應(yīng)的垂直作用點(diǎn)距離，水平力矩Zx=Dk-Zy×tanα。

在計(jì)算出破裂面交于各部位的所有主動(dòng)土壓力后，通過(guò)語(yǔ)句Eamax=max(Ea)找到最大主動(dòng)土壓力及其對(duì)應(yīng)破裂面，進(jìn)而通過(guò)查找語(yǔ)句追溯Eamax對(duì)應(yīng)的破裂角及相應(yīng)的作用點(diǎn)力矩大小。

1.4 可靠度計(jì)算功能函數(shù)的建立

(1) 抗滑移穩(wěn)定性功能函數(shù)。重力式擋土墻抗滑移穩(wěn)定性設(shè)計(jì)要求基底摩阻力大于擋墻的滑動(dòng)力。借助抗滑移安全系數(shù)定義，建立其抗滑功能函數(shù)表達(dá)式：

Z1=Mr-Ms=(Gy+Ex)μ-(Ey-Gx)=

(G×cosa0+Ea×sina0)μ-(Ea×cosa0-G×sina0)

(13)

式中：G為擋墻自重；Ea為最大主動(dòng)土壓力；Gy為垂直于基底的擋墻自重分力；μ為基底摩擦系數(shù)。

(2) 抗傾覆穩(wěn)定性功能函數(shù)。借助擋土墻抗傾覆安全系數(shù)定義可建立重力式擋土墻抗傾覆功能函數(shù)為：

Z2=Fr-Fs=GZG+EyZx-ExZy=

GZG+Ea×cosa0×Zx-Ea×sina0×Zy

(14)

式中：Zx和Zy分別為最大主動(dòng)土壓力作用點(diǎn)距墻趾的橫向和縱向力矩。

2 軟件可視化界面設(shè)計(jì)及程序編寫(xiě)

2.1 用戶界面設(shè)計(jì)

在MATLAB的命令行窗口輸入guide即可啟動(dòng)GUIDE工具。在該環(huán)境下，運(yùn)用Edit Text、Static Text、Push Button、Pop-up menu、Panel等控件設(shè)計(jì)出擋土墻計(jì)算的主界面窗口及各子界面窗口上的參數(shù)輸入、模型圖輸出及條件選擇等對(duì)話窗口[14]。各界面的菜單欄則采用Menu Editor控件設(shè)計(jì)，主界面布局見(jiàn)圖2。

圖2 主界面架構(gòu)

圖2中左上部分是已知信息輸入界面；左下部分是結(jié)構(gòu)可靠度和失效概論以及安全系數(shù)的計(jì)算結(jié)果輸出界面；右上是供擋墻計(jì)算圖形輸出預(yù)留的界面；右下部分是擋土墻計(jì)算過(guò)程中的中間結(jié)果輸出界面。本程序軟件還可以根據(jù)需要對(duì)界面進(jìn)行任意的預(yù)留或者輸出界面設(shè)計(jì)。

主子界面之間的切換則采用MATLAB中的“set (gcf,′visible′,′on′) 子界面名稱”語(yǔ)句來(lái)實(shí)現(xiàn)，而不同界面上數(shù)據(jù)的傳遞處理方法則是以兩界面間所需所有變量轉(zhuǎn)化為全局變量“Global 變量名”的形式實(shí)現(xiàn)的。

2.2 程序代碼編寫(xiě)

在設(shè)計(jì)出界面對(duì)話窗口后，可運(yùn)行GUI界面生成各個(gè)界面主程序框架，實(shí)現(xiàn)功能的代碼即可嵌入在該框架下的主函數(shù)或者各控件的回調(diào)函數(shù)架構(gòu)中。然后通過(guò)類似str2num(get(handles.編輯框編號(hào),′string′))的get語(yǔ)句獲取控件的輸入信息。在讀取所需參數(shù)后，按照重力式擋土墻可靠度計(jì)算流程進(jìn)行相應(yīng)的程序代碼編輯。首先，由各隨機(jī)參數(shù)的均值和標(biāo)準(zhǔn)差用normrnd(均值,標(biāo)準(zhǔn)差,M,N)語(yǔ)句生成M行N列服從正態(tài)分布的隨機(jī)矩陣。然后依據(jù)表1中計(jì)算公式，依次編寫(xiě)出各邊界條件下主動(dòng)土壓力Ea計(jì)算程序，最后通過(guò)以下搜索方法找到最大主動(dòng)土壓力：

j=0;

while (j

j=j+1;

if Eamax==Ea(j)

break;

end

end

w1=(j-1)/100+aaf % aaf為最大Ea對(duì)應(yīng)破裂角

擋土墻二維模型圖繪制方法是將擋墻模型主要轉(zhuǎn)折點(diǎn)的二維坐標(biāo)通過(guò)“plot([X,Y]”語(yǔ)句連接，隨后通過(guò)“hold on”語(yǔ)句將所有關(guān)鍵點(diǎn)連線繪制到同一幅坐標(biāo)圖中。最后，運(yùn)用類似“set(handles.fig1)”的set語(yǔ)句可將模型圖以及運(yùn)算結(jié)果輸出到GUI界面的指定區(qū)域[15]。擋土墻結(jié)構(gòu)抗滑和抗傾覆可靠度計(jì)算的具體程序流程設(shè)計(jì)框圖如圖3所示。

圖3 程序計(jì)算框圖

2.3 計(jì)算結(jié)果的保存與查看

擋土墻界面設(shè)計(jì)和各界面所需所有程序編寫(xiě)完成后，界面上的輸入?yún)?shù)、程序中間計(jì)算過(guò)程數(shù)據(jù)(主動(dòng)土壓力力矩、擋墻重心力矩、換算土柱高度等)和計(jì)算結(jié)果的保存可結(jié)合Excel表格實(shí)現(xiàn)。設(shè)計(jì)思路為：首先在計(jì)算機(jī)指定文件夾(如D:1 etaining)中建立一個(gè)Excel空表格(如retaining.xlsx)，然后由：“xlswrite(′D:1 etaining etaining.xlsx′,*,′Sheet1′,′I19′)”語(yǔ)句，將參數(shù)保存到該指定目錄下retaining.xlsx文件中相應(yīng)的工作表(Sheet1)中的預(yù)設(shè)位置(I19)。

界面上的曲線分析圖與擋土墻模型圖可采用以下語(yǔ)句存為JPG圖片。

axes(handles.*);

newfig=Fig. ;

set(newfig,′visible′,′off′);

set(newfig,′color′,′w′);

newaxes=copyobj(handles.fig1,newfig);

set(newaxes,′Units′,′default′,′Position′,′default′);

[filename, pathname]=uiputfile({′*.jpg′, ′Fig.type(*.jpg)′},′儲(chǔ)存圖片′);

if isequal(filename,0)||isequal(pathname,0)

return

end

str=fullfile(pathname,filename);

f=getframe(newfig);

f=frame2im(f);

imwrite(f,str);

close(newfig)

若要查看詳細(xì)的計(jì)算數(shù)據(jù)只需點(diǎn)擊GUI界面菜單欄中的“查看計(jì)算書(shū)”按鈕即可打開(kāi)保存有計(jì)算過(guò)程及計(jì)算結(jié)果的Excel表格，該功能由MATLAB語(yǔ)句“winopen(′D:1 etaining etaining.xlsx′)”實(shí)現(xiàn)。

2.4 應(yīng)用程序構(gòu)建及使用

使用MATLAB的Application Compiler工具箱可將GUI界面設(shè)計(jì)生成的主程序文件及關(guān)聯(lián)文件編譯成獨(dú)立使用的安裝包，安裝后的.exe應(yīng)用程序可以在不安裝MATLAB的情況下啟動(dòng)主界面的計(jì)算窗口。該功能的實(shí)現(xiàn)是基于一套獨(dú)立共享庫(kù)——MATLAB Runtime來(lái)完成的。此外，主界面菜單欄(見(jiàn)圖2)中的程序使用說(shuō)明書(shū)以及輸入框的提示信息可以為不熟悉該程序操作步驟的用戶提供相應(yīng)的幫助。

3 程序可行性驗(yàn)證及敏感性分析

3.1 程序可行性驗(yàn)證

某重力式擋土墻如圖4所示，擋墻各尺寸參數(shù)為：H=5.178 m，dk=0.5 m，Dk=0.98 m，hj=0.5 m，bj=0.3 m，L0=5 m，α=11.3°，墻底斜率為0.2。擋土墻填料容重23 kN/m3。

圖4 重力式擋土墻計(jì)算簡(jiǎn)圖

隨機(jī)變量的各參數(shù)具體值見(jiàn)表2，假設(shè)表2中各參數(shù)均服從正態(tài)分布。

表2 各隨機(jī)變量參數(shù)統(tǒng)計(jì)

在不計(jì)荷載的情況下，將已知參數(shù)輸入計(jì)算界面，通過(guò)運(yùn)行計(jì)算可得到對(duì)應(yīng)的計(jì)算結(jié)果如圖2所示。

為進(jìn)一步驗(yàn)證計(jì)算程序的正確性與可行性，現(xiàn)將本文所開(kāi)發(fā)的軟件與工程上常用的理正擋墻軟件就該算例的定值計(jì)算結(jié)果作對(duì)比分析，如表3所示。

表3 本文軟件和已有軟件計(jì)算結(jié)果對(duì)比

由表3可知，本文開(kāi)發(fā)軟件的定值計(jì)算結(jié)果與理正軟件計(jì)算結(jié)果相仿。此外，本文軟件還可以計(jì)算出結(jié)構(gòu)在兩種模式下的具體失效概率和可靠指標(biāo)，為實(shí)際工程擋土墻結(jié)構(gòu)的安全可靠性評(píng)價(jià)提供可靠理論依據(jù)。

3.2 敏感性分析

通過(guò)點(diǎn)擊穩(wěn)定性分析按鈕可以切換到對(duì)該算例的參數(shù)敏感性分析窗口，如圖5所示。

圖5 參數(shù)敏感性分析計(jì)算窗口

假設(shè)各參數(shù)均值不變，變異系數(shù)以設(shè)定步長(zhǎng)變化，可得擋土墻抗傾覆與抗滑移可靠指標(biāo)隨各參數(shù)變異系數(shù)變化的曲線如圖6-圖9所示。

圖6 可靠指標(biāo)與容重變異系數(shù)關(guān)系

圖7 可靠指標(biāo)與填土內(nèi)摩擦角變異系數(shù)關(guān)系

圖8 可靠指標(biāo)與基底摩擦系數(shù)變異系數(shù)關(guān)系

圖9 可靠指標(biāo)與填土和墻背間內(nèi)摩角變異系數(shù)關(guān)系

在對(duì)四個(gè)主要影響因素進(jìn)行敏感性分析后，由圖6-圖9可知：填土容重變異系數(shù)、基底摩擦系數(shù)變異系數(shù)、填土內(nèi)摩擦角變異系數(shù)三者在0.4之前變化時(shí)對(duì)抗滑移穩(wěn)定性的影響比較大，填土容重變異系數(shù)、填土內(nèi)摩擦角變異系數(shù)兩者在0.4之前變化時(shí)對(duì)抗傾覆穩(wěn)定性的影響比較大，而填土與墻背間的內(nèi)摩擦角變異系數(shù)則剛好相反。因此，在實(shí)際設(shè)計(jì)時(shí)，若相應(yīng)參數(shù)變異系數(shù)低于0.4時(shí)，應(yīng)準(zhǔn)確統(tǒng)計(jì)其變異系數(shù)，方能更客觀地分析結(jié)構(gòu)的可靠度。

4 結(jié) 語(yǔ)

(1) 基于蒙特卡羅法，結(jié)合MATLAB軟件矩陣算法、GUI界面設(shè)計(jì)等功能，開(kāi)發(fā)一款重力式擋土墻結(jié)構(gòu)可靠度計(jì)算及敏感性分析的計(jì)算機(jī)程序軟件，彌補(bǔ)了重力式擋土墻結(jié)構(gòu)傳統(tǒng)定值分析法的不足，為進(jìn)行擋墻結(jié)構(gòu)可靠度計(jì)算提供了適用計(jì)算工具。

(2) 本文開(kāi)發(fā)的程序可以計(jì)算不同類型(仰斜、俯斜、垂直、有無(wú)擴(kuò)展墻趾、是否計(jì)算路面荷載)重力式擋土墻的結(jié)構(gòu)可靠度。穩(wěn)定性計(jì)算模塊不僅可以計(jì)算擋土墻結(jié)構(gòu)在兩種失效模式下的失效概率和可靠指標(biāo)，還可以計(jì)算出基于定值理論下抗滑移和抗傾覆安全系數(shù)。

(3) 本文軟件開(kāi)發(fā)生成了以計(jì)算過(guò)程和計(jì)算結(jié)果為主要內(nèi)容的擋土墻結(jié)構(gòu)可靠度計(jì)算書(shū)。查看計(jì)算書(shū)，可判斷擋土墻設(shè)計(jì)是否符合規(guī)范要求。此外，軟件將擋土墻模型圖及分析曲線圖也一并存儲(chǔ)在相應(yīng)的文件夾下，可供隨時(shí)查看。

(4) 通過(guò)一個(gè)典型擋土墻算例的結(jié)構(gòu)可靠度計(jì)算和敏感度分析，獲得了隨機(jī)參數(shù)取值及其變異系數(shù)范圍對(duì)結(jié)構(gòu)可靠性和安全性的影響規(guī)律，驗(yàn)證了本文開(kāi)發(fā)程序軟件的正確性和實(shí)用性。