界限含水率程序設計

1、界限含水率試驗數據處理技術10界限含水率概述界限含水率試驗是土工試驗中最常見、應用最廣的試驗之一。土的液塑限指標是細粒土進行分類和定名的最基本指標，在巖土工程中，液塑限指標的準確性涉及土壤定名的正確性，同時影響土樣狀態(tài)的確定，進而影響到土的承載能力的確定；在公路鐵路路基工程中，還直接影響細粒土的填料分組，所以正確地確定土壤的液塑限指標對工程具有很重要的意義。21 界限含水率試驗數據處理原理 根據極限平衡理論，沿圓錐與土的接觸面的極限剪應力（見圖1）為：Prlh3 （1）P圓錐重力錐角度數A圓錐與土的接觸面積，其計算方法如下： r為圓錐入土的半徑h入土深度l圓錐入土的斜邊長1 界限含水率試驗數據

2、處理原理 4一般的，由于為定值，令 則 （2） 對式(2)兩邊取對數得： （3）其中l(wèi)gKP為常數。由于重塑土的無側限抗壓強度與含水率存在雙對數關系，即 （4）1 界限含水率試驗數據處理原理 5將上式繪成雙對數線，是一條直線。對多種土進行不同含水率和抗剪強度與對應圓錐入土深度試驗，結果表明，理論曲線與試驗曲線在特性上相一致。經整理得： （5） 其中a和b是由試驗擬合確定的常數，根據以上原理，界限含水率試驗應采用雙對數坐標系。1 界限含水率試驗數據處理原理 6目前界限含水率數據處理常用的方法是規(guī)范作圖法，此外還有直接線性最小二乘法和公式法。下面列舉某次試驗中的兩組數據(見表1)，并且用以上各種方

3、法計算界限含水率偏差。錐沉量h(mm)含水率(%)錐沉量h(mm)含水率(%)第一組數據第二組數據3.4022.595.6036.709.2030.279.4543.2815.9036.0315.6753.212 試驗數據處理方法及偏差計算7直接線性擬合函數簡單，在Excel表中處理方便，有些試驗人員在試驗時直接用線性擬合（采用了最小二乘法），且在Excel表中有時擬合度達到0.99以上。但根據以上原理，顯然直接線性擬合是不合理的，且不能估計塑限差距是否滿足要求，造成數據計算偏差較大。表2是直線擬合法的計算結果。2.1直接線性擬合法第一組數據第二組數據擬合方程y = 1.0697x + 19.

4、468擬合方程y = 1.636x + 27.644擬合度R2=0.9847擬合度0.9997p (%)21.61p (%)30.9210mmL(%)30.1710mmL(%)44.0017mmL(%)37.6517mmL(%)55.468目前界限含水率試驗數據處理一般采用TB10102-2004鐵路工程土工試驗規(guī)程法，以含水率為橫坐標，圓錐下沉深度為縱坐標，在雙對數坐標紙上繪制如圖2所示的關系曲線，三點應連成一條直線。當三點不在一條直線上，則通過高含水率這一點與其余兩點連成兩條直線，在圓錐下沉深度為2mm處，可查得相應的兩個含水率，當這兩個含水率的差值小于2%時，應以這兩點的含水率平均值與高

5、含水率的點連成一條直線。當這兩個含水率之差值大于或等于2%時，則應重做試驗。 2.2 規(guī)范作圖法第一組數據11010010100w(%)h(mm)9規(guī)范作圖法是土工試驗規(guī)范中所建議的方法，但該法有如下缺點：作圖在雙對數紙上進行，由于對數坐標不均勻分格，1、2之間和3、4之間的間距不相同，作圖時一般是按線性讀數，顯然會造成讀數偏差；讀數受人為因素影響較大；對2%的要求無法計算；當塑性指數處于10和17附近時，直接影響了土的定名。所以在試驗中使用此法時求解不確定性較大。表3 作圖法計算第一組數據第二組數據塑限(%)18.50塑限(%)25.0010mmL(%)32.0010mmL(%)45.001

6、7mmL(%)37.517mmL(%)55.002.2 規(guī)范作圖法10設試驗參數如圖3所示。為便于計算，我們將lg(h)作為自變量，lg()作為因變量。AB的直線方程為同理AC的方程為 2.3 公式法(6)(8)11分別把h2mm代入式(6)(8)可得AB2和AC2，則由兩條直線引起的塑限偏差值為：p=AB2-AC2 (10) 當p2%公式法的優(yōu)點是能明確p的大小，控制試驗是否需要重做。但該法以錐沉量最大者為準確值，調整B、C兩點顯然也不盡合理，且不能保證三點按規(guī)范所擬合的誤差最小，特別是塑限偏差值接近2%時。第一組數據第二組數據塑限(%)18.93塑限(%)24.1410mmL(%)31.1

7、210mmL(%)44.7817mmL(%)36.7917mmL(%)54.352.3 公式法13公式法計算Dim Kab As Double, Kac As Double, Wab2 As Double, Wac2 As Double, Wpcha As DoubleDim ha As Double, hb As Double, hc As Double, wa As Double, wb As Double, wc As Doubleha = Val(Text48)hb = Val(Text47)hc = Val(Text46)wa = Val(Text51) / 100wb = Val(

8、Text50) / 100wc = Val(Text49) / 100If ha = 0 Or hb = 0 Or hc = 0 Or wa = 0 Or wb = 0 Or wc = 0 Or ha = hb Or hb = hc Or wa = wb Or wb = wc Then MsgBox 原始數據輸入有誤，請重輸 Exit SubEnd If輸入ha,hb,hcWa,wb,wc判別輸入數據有沒有空值14Kab = (Lg(wa) - Lg(wb) / (Lg(ha) - Lg(hb)Kac = (Lg(wa) - Lg(wc) / (Lg(ha) - Lg(hc) 兩段直線的斜率W

9、ab2 = 10 (Kab * Lg(2) + Lg(wa) - Kab * Lg(ha)Wac2 = 10 (Kac * Lg(2) + Lg(wa) - Kac * Lg(ha) 計算兩段直線所對應的塑限Wpcha = (Wab2 - Wac2)塑限差值wText60(0) = Format(Wpcha * 100, #0.000)Dim tmp As Double, Wp As Double, Wl10 As Double, Wl17 As Doubletmp = (Lg(wa) - Lg(Wab2 + Wac2) / 2) / (Lg(ha) - Lg(2)tmp = tmp * (L

10、g(2) - Lg(ha) + Lg(wa)Wp = 10 tmpText60(1) = Format(Wp * 100, #0.000) 塑限求解Public Function Lg(ByVal x As Double) As DoubleLg = Log(x) / Log(10)End Function公式法的優(yōu)點，可以判定符合規(guī)范與否15tmp = (Lg(wa) - Lg(Wab2 + Wac2) / 2) / (Lg(ha) - Lg(2)tmp = tmp * (Lg(10) - Lg(ha) + Lg(wa)Wl10 = 10 tmpText60(2) = Format(Wl10

11、 * 100, #0.000) 10mm液限求解tmp = (Lg(wa) - Lg(Wab2 + Wac2) / 2) / (Lg(ha) - Lg(2)tmp = tmp * (Lg(17) - Lg(ha) + Lg(wa)Wl17 = 10 tmpText60(3) = Format(Wl17 * 100, #0.000) 17mm液限求解Text60(4) = Format(Wl10 - Wp) * 100, #0.000)Text60(5) = Format(Wl17 - Wp) * 100, #0.000) 塑性指數的求解特點：變量隨時可以定義，但在一個過程中不能兩次聲明一個變量

12、162.4先求對數再最小二乘擬合法 最小二乘法是德國數學家C.F.高斯在1794年解決行星軌道預測問題時首先提出的。線性擬合的最小二乘法基本原理為：設經驗方程是y=f(x)，方程中含有一些待定系數an，給出真實值(xi,yi)|i=1,2,.n，將這些x，y值代入方程然后作差，可以描述誤差：yi- f (xi)，為了考慮整體的誤差，可以取平方和，之所以要平方是考慮到誤差可正可負直接相加可以相互抵消，所以記誤差為： 172.4先求對數再最小二乘擬合法 如果經驗方程是線性的，形如y=ax+b，就是線性回歸。按上面的分析，誤差函數為： 于是得到關于a ，b的線性方程組： 要令e為最小，則e對a和b求

13、偏導數為零得：n為數據點數設 (13)(14)(15)(16)(17)(18)182.4先求對數再最小二乘擬合法 則方程化簡為： (19) 解出上述方程組得：根據以上原理，利用Visual Basic編寫的最小二乘法線性擬合的關鍵函數程序段如下：192.4先求對數再最小二乘擬合法 Private Sub Nihe(x() As Double, y() As Double, a As Double, b As Double, e As Double)Dim n As IntegerDim aa As Double, bb As Double, cc As Double, dd As Double

14、, delta As Doubleaa = 0: bb = 0: cc = 0: dd = 0n = UBound(x) + 1For i = 0 To UBound(x) aa = aa + x(i) * x(i) bb = bb + x(i) cc = cc + x(i) * y(i) dd = dd + y(i)Nextdelta = aa * n - bb * bbIf Abs(delta) 0.0000000001 Then MsgBox errorElse a = (cc * n - bb * dd) / delta b = (aa * dd - cc * bb) / deltaE

15、nd Ife = 0For i = 0 To UBound(x)求解誤差 e = e + (y(i) - a * x(i) - b) 2Next iEnd Sub根據原理編寫如何調用過程（函數）20對上述函數調用Private Sub Command18_Click()Dim x(2) As Double, Y(2) As DoubleDim a As Double, b As Double, e As DoubleDim Wl10 As Double, Wl17 As Double, Wp2 As Double, Ip As Doublex(0) = Lg(Val(Text48)x(1) =

16、 Lg(Val(Text47)x(2) = Lg(Val(Text46)Y(0) = Lg(Val(Text51)Y(1) = Lg(Val(Text50)Y(2) = Lg(Val(Text49) 先求對數再擬合不求對數就是直接線性擬合法21過程調用續(xù) Nihe x(), Y(), a, b, e或者call nihe(x(),y(),a,b,e)過程調用的兩種寫法22反求液塑限續(xù)Wp2 = 10 (a * Lg(2) + b)Wl10 = 10 (a * Lg(10) + b)Wl17 = 10 (a * Lg(17) + b) 反求液塑限Ercheng(1) = Format(Wp2,

17、#0.000)Ercheng(2) = Format(Wl10, #0.000)Ercheng(3) = Format(Wl17, #0.000)Ercheng(4) = Format(Wl10 - Wp2, #0.000)Ercheng(5) = Format(Wl17 - Wp2, #0.000)Ercheng(0) = Format(1 - e, 擬合度 & #0.0000)23任意沉降對應的含水量Dim Wren As DoubleIf Val(Text55) = 30 Thentext56.text= 出錯Exit SubElsetmp = 10 (a * Lg(Val(Text55

18、) + b)text56.text= Format(tmp, 0.#) 任一沉降對應的含水率求解End IfEnd Sub242.4先求對數再最小二乘擬合法 注：x()數組是A、B、C三點的橫坐標的對數，y()數組是相應的縱坐標的對數，a、b返回要求解的線性方程系數，e返回誤差的平方和。把A(lgA,lghA)，B(lgB,lghB)，C(lgC,lghC)三點坐標組成的數組x(0)=lgA，x(1)=lgB，x(2)=lgC和y(0)=lghA，y(1)=lghB，x(2)=lghC兩個數組調用上述函數，即可得出擬合后的線性方程參數a，b。把y=lg(2)及y=lg(10)代入擬合公式可得塑限和液限。計算結果見表5。但該法也有自身的缺點，即擬合度很高時，也無法保證p2%的規(guī)范要求。25第一組數據第二組數據塑限(%)19.21塑限(%)25.1210m