應(yīng)用MATLAB確定鄧肯-張雙曲線模型中的K-n參數(shù)

應(yīng)用MATLAB確定鄧肯－張雙曲線模型中的K,n參數(shù)簡(jiǎn)介：接合承德中密砂常規(guī)三軸試驗(yàn)數(shù)據(jù)，介紹應(yīng)用Matlab語言編寫計(jì)算及繪圖程序來處理試驗(yàn)數(shù)據(jù)的方法，可顯著提高試驗(yàn)研究的數(shù)據(jù)處理效率和結(jié)果的可視化程度。

關(guān)鍵字：Matlab三軸試驗(yàn)鄧肯－張模型1前言基于廣義胡克定律的線彈性理論形式簡(jiǎn)單，參數(shù)少，物理意義明確，而且在工程界有廣泛深厚的基礎(chǔ)，得以應(yīng)用于許多工程領(lǐng)域中。早期土力學(xué)中的變形計(jì)算主要是基于線彈性理論的，只有在計(jì)算機(jī)得到迅速發(fā)展之后，非線性理論模型才得到較廣泛的應(yīng)用。鄧肯－張模型是建立在增量廣義胡克定律基礎(chǔ)之上的變模量的彈性模型，可以反映土變形的非線性，并在一定程度上反映土變形的彈塑性，很容易為工程界所接受，加之所用參數(shù)和材料參數(shù)不多，物理意義明確，只需用常規(guī)三軸壓縮試驗(yàn)即可確定這些參數(shù)及材料常數(shù)適應(yīng)的土類比較廣，所以該模型為巖土工程界所熟知，并得到了廣泛的應(yīng)用，成為土的最為普及的本構(gòu)模型之一。本文主要是應(yīng)用MATLAB編寫計(jì)算及繪圖程序來處理承德中密砂常規(guī)三軸試驗(yàn)數(shù)據(jù)。2基于MATLAB的計(jì)算過程實(shí)現(xiàn)現(xiàn)場(chǎng)的觀測(cè)數(shù)據(jù)經(jīng)過采集和整理后，按照一定的格式把數(shù)據(jù)存儲(chǔ)在數(shù)據(jù)文件中，然后可以使用MATLAB豐富的數(shù)值運(yùn)算功能可以非常容易地編制出數(shù)據(jù)處理程序，先用函數(shù)fopen()打開數(shù)據(jù)文件，fid=fopen(‘filename’,’r’)再用fscanf函數(shù)依次從文件中讀取格式化數(shù)據(jù)來完成對(duì)各變量地賦值，其使用語法為：matrix=fscanf(fid,format)。本文由于數(shù)據(jù)不是太多，所以在計(jì)算過程中沒有采取調(diào)用存儲(chǔ)文件地形式。直接在計(jì)算過程中輸入試驗(yàn)數(shù)據(jù)計(jì)算。2.1數(shù)據(jù)的處理對(duì)第一組數(shù)據(jù)，通過編寫Matlab語言，由軸向應(yīng)變和應(yīng)力差的試驗(yàn)數(shù)據(jù)可以作出～（）和～雙曲線關(guān)系圖形，主要用到的MATLAB命令為：plot(x1,y1)；axis([00.0403])；holdon%（1）plot(x1,x1./y1)；a=polyfit(x1,x1./y1,1)；t1=0:0.001:0.07;plot(x1,x1./y1,'.',t1,a(1)*t1+a(2))%（2）其中x1代表第一組軸向應(yīng)變，x2代表第一組應(yīng)力差。（具體的Matlab語言見后附錄），其他兩組的數(shù)據(jù)可以采取同樣的方法進(jìn)行處理。三種不同圍壓下生成的曲線如圖1所示。

3曲線擬合通過應(yīng)用最小二乘法，利用Matlab中polyfit函數(shù)進(jìn)行多項(xiàng)式擬合，具體命令為：Polyfit(X,Y,n),其中X、Y為擬合數(shù)據(jù)，n為擬合多項(xiàng)式的階數(shù)。對(duì)生成的曲線進(jìn)行數(shù)據(jù)擬合，生成的擬合曲線如圖2所示。圖中的三條線性曲線的截距為試驗(yàn)中的初始變形模量的倒數(shù)。繪出與的關(guān)系圖，兩者近似成直線關(guān)系，見圖3。所以可得下式：，其中，為大氣壓（=101.4kpa）,量綱與相同；K、n為試驗(yàn)常數(shù)，分別代表與直線的截距和斜率。通過對(duì)大量的數(shù)據(jù)信息轉(zhuǎn)化為直觀、易理解的圖象信息，便于對(duì)客觀規(guī)律的認(rèn)識(shí)。使用下面的語句繪制圖2的圖形：figure；plot（，）；xlabel（’’,）;ylabel(‘’）;至此，已經(jīng)求出鄧肯－張模型中的兩個(gè)重要參數(shù)K、n。

4結(jié)論就本文介紹的用Matlab方法處理承德中密砂常規(guī)三軸試驗(yàn)數(shù)據(jù)而言，雖然從計(jì)算過程并不能看出其明顯區(qū)別與其他程序語言強(qiáng)大的矩陣運(yùn)算功能。但是，它實(shí)現(xiàn)了從數(shù)據(jù)到圖形的無縫對(duì)接，體現(xiàn)了Matlab對(duì)計(jì)算結(jié)果可視化的無可比擬的優(yōu)勢(shì)，可以極大的提高諸如巖土工程等各個(gè)工程技術(shù)領(lǐng)域試驗(yàn)研究的工作效率和研究水平。附表承德中密砂常規(guī)三軸試驗(yàn)數(shù)據(jù)表(圍壓σ3=100kPa)應(yīng)力差軸向應(yīng)變體應(yīng)變?chǔ)?-σ3（100kPa)ε1εv0.508.00225.000741.002.00449.001301.463.00674.001761.849.00898.002232.149.01123.002232.331.01348.002142.477.01572.001762.587.01797.001582.665.02021.001112.730.02246.000562.777.02471.000092.812.02695.000032.845.03032‐.000122.878.03369‐.001952.889.03706‐.002882.894.04043‐.003812.889.04492‐.004742.879.04941‐.005672.869.05390‐.006502.853.05839‐.007522.836.06289‐.008272.809.06738‐.00891承德中密砂常規(guī)三軸試驗(yàn)數(shù)據(jù)表(圍壓σ3=300kPa)應(yīng)力差軸向應(yīng)變體應(yīng)變?chǔ)?-σ3（100kPa)ε1εv0.909.00125.000752.286.0035.001513.533.00575.002364.3410.008.003025.203.01025.003295.781.01251.003296.202.01476.003296.529.01701.003296.789.01926.003027.065.02151.002837.233.02376.002367.390.02601.001897.515.02939.001607.714.03276.000947.846.03614.000287.939.03952‐.000478.016.04402‐.001238.043.04852‐.002078.058.05302‐.002928.047.05753‐.003778.026.06203‐.004527.992.06653‐.005377.932.07216‐.006137.859.07778‐.00707承德中密砂常規(guī)三軸試驗(yàn)數(shù)據(jù)表(圍壓σ3=500kPa)應(yīng)力差軸向應(yīng)變體應(yīng)變?chǔ)?-σ3（100kPa)ε1εv0.756.00025.000572.454.0025.001144.587.00476.001996.349.00701.002757.689.00926.003038.751.01151.003419.573.01376.0035110.224.01602.0036010.759.01827.0036011.197.02052.0034111.551.02277.0031311.837.02502.0029412.164.02840.0024612.439.03178.0019912.652.03516.0015212.808.03854.0010412.953.04304.0000913.069.04754‐.0006613.148.05205‐.0017113.21.05655‐.0024613.239.06106‐.0032213.234.06556‐.0040813.211.0719‐.0042313.117.07682‐.00569附錄>>x1=[0.002250.004490.006740.008980.011230.013480.015720.017970.020210.022460.024710.026950.030320.033690.037060.040430.044920.049410.05390.058390.062890.06738];>>x1=x1';>>y1=[0.5081.0021.4631.8492.1492.3312.4772.5872.6652.732.7772.8122.8452.8782.8892.8942.8892.8792.8692.8532.8362.809];>>y1=y1';>>plot(x1,y1)>>axis([00.0403])>>x1./y1ans=Columns1through40.00440.00450.00460.0049Columns5through80.00520.00580.00630.0069Columns9through120.00760.00820.00890.0096Columns13through160.01070.01170.01280.0140Columns17through200.01550.01720.01880.0205Columns21through220.02220.0240>>plot(x1,x1./y1)>>a=polyfit(x1,x1./y1,1)a=0.31100.0019>>t1=0:0.001:0.07;>>plot(x1,x1./y1,'.',t1,a(1)*t1+a(2))>>x2=[0.001250.00350.005750.0080.010250.012510.014760.017010.019260.021510.023760.026010.029390.032760.036140.039520.044020.048520.053020.057530.062030.066530.072160.07778];>>y2=[0.9092.2863.5334.3415.2035.7816.2026.5296.7897.0657.2337.3907.5157.7147.8467.9398.0168.0438.0588.0478.0267.9927.9327.859];>>plot(x2,y2)>>axis([00.0509])>>x2./y2;>>plot(x2,x2./y2)>>b=polyfit(x2,x2./y2,1)b=0.11130.0008>>t2=0:0.001:0.07;>>plot(x2,x2./y2,'.',t2,b(1)*t2+b(2))>>x3=[0.000250.00250.004760.007010.009260.011510.013760.016020.018270.020520.022770.025020.028400.031780.035160.038540.043040.047540.052050.056550.061060.065560.07190.07682];>>y3=[0.7562.4544.5876.3497.6898.7519.57310.22410.75911.19711.55111.83712.16412.43912.65212.80812.95313.06913.14813.2113.23913.23413.21113.117];>>plot(x3,y3);>>axis([00.06014])>>c=polyfit(x3,x3./y3,1)c=0.06700.0005>>t3=0:0.001:0.07;>>plot(x3,x3./y3,'.',t3,c(1)*t3+c(2))>>E=[1/a(2)1/b(2)1/c(2)]>>o3=[135]>>E=[1/(a(2)*10^5)1/(b(2)*10^5)1/(c(2)*