基于低围压土体的三轴试验数据拟合及其软件实现

董俊言，薛大文，毛明锋，蔡足铭

(浙江海洋大学港航与交通运输工程学院，浙江舟山 316022)

岩土工程中，岩土的三轴试验是是研究岩土体材料破坏规律的基础实验，其主要是通过拟合求解三轴试验应力摩尔圆包络线以获取土样的粘聚力c和内摩擦角φ。低围压情况下(σ3＜10 MPa)，岩土体破坏的摩尔库伦线性模型是在实际应用中最为常用的模型[1]。作为低围压土体的典型代表，软土在我国分布广泛，沿海、内陆湖盆、洼地及河流两岸地域随处可见，如长江、珠江三角洲地带，黄淮冲击平原，洞庭湖、洪泽湖等内陆湖盆等均富含软土分布。随着我国社会经济发展，人们与海洋软土资源的接触日益频繁。因此，对包含海洋软土在内的低围压土体的三轴试验数据处理具有及其重要的现实意义。

在常用的岩土三轴试验数据处理中，通常做法是，对同一种岩土样的多个试件分别在不同的周围压力σ3下进行试验[2]。土样在周围压力σ3的作用下，通过增大轴向压力σ1，直至破坏。通过该方式可以做出一个摩尔圆。在此基础上，改变围压σ3，重复上一步骤，得到多个应力摩尔圆。此后通过求解这些应力圆的公切线，根据包络线的截距和斜率，以此获得该土样粘聚力c和内摩擦角φ的数值。

目前来说，已有众多学者研究了对三轴试验摩尔圆的不同的处理方法。邹宗兴等[3]采用先对相邻的应力摩尔圆求得c、φ值，然后再通过求均值来获得土样的抗剪强度参数。该方法计算较为简单，但其结果较为粗糙，精度较低，误差较大。《水利水电岩石试验规程》[4]中的做法是，对同组的试验中的σ1、σ3作最小二乘线性拟合，在直线上等间距取5～6组σ＇1、σ＇3，再做摩尔圆，最后再做这些圆的公切线得到c、φ值。该方法精度较文献[3]的方法高，但操作复杂，计算繁琐。王星辉[5]采用三轴强度指标回归的方法，求出每两圆公切线的切点，然后拟合所有的切点作求得参数c、φ值。此外，还有杨同等[6]、王艳芳等[7]、吴力平等[8]、李驰等[9]均提出了各自不同的方案对三轴试验数据处理。总体来说，此类方法都操作可行，但对试验数据的可靠性论证稍显不足。马惠彪等[10]、田兆阳等[11]、孙红等[12]、陈超斌等[13]均对海洋软土进行过三轴试验研究，但他们并未提出对三轴试验数据编程拟合及其软件实现。

综上所述，本文基于最小二乘法对低围压土体的三轴试验数据进行线性拟合，并对拟合程序进行软件实现。文章重点对低围压土体摩尔圆包络线的线性拟合进行详细的公式推导，采用牛顿迭代求解拟合过程中出现的非线性方程组以此获得抗剪强度参数粘聚力和内摩擦角。在此基础上，通过MFC类库搭建用户界面，将以上计算方法实现软件化，以此替代手工图解或计算机交互图解，促进并指导相应专业和领域中的应用。

1 线性拟合方法详解

1.1 图表法

如图1所示，图中摩尔圆代表了受力物体内某一点的应力情况，其中，横坐标的数值为正应力值，而纵坐标的数值为剪应力值。经过截面上的正应力值和剪应力值可以由摩尔圆圆周上所对应的坐标来表示。由于摩尔圆的圆周有无数点，其代表了无数个土样不同斜度的截面上的应力值，那么整个圆就代表了一个应力状态。半径长度代表了最大剪应力，而圆心到坐标原点的距离则代表了平均应力，差应力（σ3-σ1）的值为圆的直径。通过手工绘制在图上，观察摩尔圆的大小和距坐标原点的远近，则可以大致了解其应力状态。由图像分析得出，粘聚力与轴向破坏应力大致成正比关系。该方法清晰易懂，但操作繁琐，耗时较长，对三轴试验数值的求解精确度低，因而逐步形成以计算机求解为主的拟合方法（图中摩尔圆第2个角标表示为摩尔圆的序数）。

图1 摩尔应力圆线性拟合图表法Fig.1 Graphical method for linear fitting of Mohr stress cycles

1.2 算法推导

本文通过最小二乘法对摩尔圆包络线线性拟合，采用牛顿迭代求解拟合过程中出现的非线性方程组从而获得较为精确数值。具体计算思路和算法推导如下：

定义σ1和σ3为岩土正应力，假设拟合直线为（k为斜率，c为截距）：

为了求出上式极小值，分别对c和k求偏导得到：

为计算简便，令

式子中的σoi为对k求偏导的表示式子，由于上述对偏导的求解过程中出现非线性项，采用牛顿迭代的方法对该方程组迭代求解，以求得近似的土样粘聚力c和内摩擦角φ的数值。牛顿迭代法[14]具体如下：

设z=f(x,y)在点的(x0,y0)某一邻域内连续且直到2阶的连续偏导数，(x0+k，y0+k)为此邻域内任一点，则有：

则方程可表示为z=0可近似表示为：

同理z=g(x，y)可得到如下方程式：

方程z=g(x，y)=0可表示为

由此得到了一组方程组

求解出如下的方程组为：

记符号为：

求得的迭代公式为：

2 软件实现及简介

对上述算法在计算机中通过编写C++程序，程序主要分为两部分(主函数和子函数)。主函数部分主要包括三轴试验数据输入，摩尔圆及拟合直线输出。子函数部分包含牛顿迭代法的循环迭代过程。主函数通过调用子函数处理所输入的三轴试验大、小主应力数据，从而求得土样的抗剪强度指标，即粘聚力和内摩擦角。在此基础上，再通过MFC类库搭建用户界面，将以上计算方法实现软件化。本程序软件界面及运行的效果如图2所示。图中窗口左边为主应力输入框，目前最多可分别输入大、小主应力各6组的数据。右边上半部的静态文本框用于输出结果，分别输出拟合直线的斜率、内摩擦角（弧度和角度）以及粘聚力（单位为kPa）。右边下半部分的图框用于显示摩尔应力圆以及所拟合的直线，通过该图框图像显示拟合的过程与结果。

本软件程序除了能实现各摩尔圆的线性拟合之外，还能检验三轴试验数据的可靠性。如图2b和2c所示为 3 组主应力值和 4 组主应力值测试结果，主应力值分别为(20.5，10.45)，(50.65，30.98)，(90.15，60.63)和(20，10)，(50，30)，(86，53)，(62，38)。通过图中的图像拟合结果看出，各摩尔圆与拟合直线相切程度较高，表明这2组三轴试验数据可靠性较高，真实有效。反观图2d中，对于该五组数值拟合结果，有部分摩尔圆与拟合直线相切程度不高，这是由于程序算法的本质正是为了使拟合直线与所有圆的相切程度最高所致。从图中可以看出，其中第五组数据（77，43）较其余4组数据可靠性低，误差较大，系测量粗大误差，在试验数据处理时应剔除。

图2 摩尔圆线性拟合软件运行效果图Fig.2 Results of the linear fitting of Mohr circles with the s software

3 实例验证

3.1 实例1

淮河入江水道软土数据见表1。

表1 淮河入江水道的软土数据[15]Tab.1 Data of soft soil in Huaihe River

根据软件拟合结果(图3)，得到该软土的粘聚力c为6.370 6 KPa，内摩擦角φ为12.261 6°，实验测得的总应力抗剪强度指标c为6.4 KPa，内摩擦角φ为12.3°。软件拟合数值相对误差分别为0.45%和0.31%。

3.2 实例2

天津港某软黏土三轴试验数据见表2。

图3 淮河入江水道的软土拟合结果Fig.3 Fitting results of soft soil in Huaihe River

表2 某软黏土土样最大及最小主应力[8]Tab.2 Maximum and minimum principal stress

软件拟合结果(图4)显示该软黏土的粘聚力c为10.704 KPa，内摩擦角φ为19.513°，实验测得的总应力抗剪强度指标c为11.1 KPa，内摩擦角φ为20.1°。软件拟合数值相对误差分别为3.56%和2.92%。

3.3 实例3

某泥质砂岩三轴试验数据见表3。

软件拟合结果(图5)显示该软黏土的粘聚力c为47.570 KPa，内摩擦角φ为25.956°，实验测得的总应力抗剪强度指标c为45.33 KPa，内摩擦角φ为26.9°。软件拟合数值相对误差分别为4.94%和3.51%。

图4 软件拟合效果图Fig.4 The fitting results of the software

图5 软件拟合效果图Fig.5 The fitting results of the software

表3 泥质砂岩数据[2]Tab.3 Data of argillaceous sandstone

4 结论

低围压条件下的河口、海洋软土通常均符合摩尔库伦准则，此时土体轴向破坏应力与围压呈线性关系。为了对该类土体获得较准确的抗剪强度指标，本文对三轴试验数据进行最小二乘拟合，通过C++程序对摩尔应力圆线性拟合进行程序编写，采用牛顿迭代法处理非线性方程组，并通过MFC类库搭建用户界面将程序过程软件实现。最后对不同河域、海域软土试验数据进行实例验证，以对比验证软件的准确度和精确性。本文总体结论如下：

(1)本文软件程序除了能实现对摩尔应力圆的线性拟合外，还能检验三轴试验数据的可靠性，可为实验数据处理中判别测量误差，删选有效实验数据提供帮助。

(2)本文采用的计算方法精度高，迭代时间短，收敛速度快。通过实例对比验证，软件拟合数值相对误差在5%以下。

(3)本软件界面清晰，操作简单，使用方便。与传统的人工拟合方法或计算机交互图解方法相比，本软件程序优势明显。本软件广泛适用于各低围压软土土体，可直接促进并指导相应专业和领域中的应用。