用MATLAB做三轴试验强度包线最小二乘法拟合φ和c的完整套路做三轴剪切试验写报告或者做毕设的同学到了数据处理这一关十有八九会被“画莫尔圆、作公切线、读内摩擦角φ和粘聚力c”这套流程折磨过。早几年我自己也是拿直尺和三角板在坐标纸上一条条贴公切线贴完还要用两把尺子量角度费时不说不同人贴出来的切线位置差个一两度都是常事——而这一两度换算到φ和c上面结论可能就差出一大截。后来我彻底改成用MATLAB处理输入每级围压和破坏时主应力差用最小二乘法自动拟合强度包线φ和c直接出数值图也顺手画好省下来的时间用来复查数据质量比手工画图划算得多。这篇就把整套思路和代码完整写出来讲清楚每一行在干什么、为什么这么算以及哪些坑我踩过之后希望你直接绕开。先说适用对象手里有三轴试验数据、需要按摩尔-库仑准则求总应力或有效应力抗剪强度指标内摩擦角φ、粘聚力c的本科生、研究生和工程师。代码只用到MATLAB最基础的polyfit和plot函数不依赖任何额外工具箱R2016b之后的版本都能跑。即使你MATLAB还不熟按下面的步骤抄也能出图出数。1. 先搞懂原理三轴试验、莫尔圆和强度包线怎么串起来1.1 三轴试验到底能给我们哪些数常规三轴压缩试验的做法是把3到4个相同土样分别放在不同的恒定围压σ3下然后施加轴向压力让试样剪切破坏。围压就是压力室里的水压一般取100、200、300、400 kPa这一档。每个试样剪切过程中仪器记录轴向应变和主应力差σ1-σ3的关系曲线试验结束时从曲线上取破坏点对应的主应力差σ1-σ3f——如果曲线有峰值就取峰值没有峰值就按15%轴向应变对应的值取。有了围压σ3和破坏时主应力差σ1-σ3f就能算出试样破坏时的大主应力σ1f σ3 (σ1-σ3)f一组围压对应一个试样一个试样就能画出一个莫尔圆。比如我用过的典型数据试样编号围压σ3 (kPa)破坏主应力差 (σ1-σ3)f (kPa)破坏大主应力σ1f (kPa)11002853852200455655330062592544007901190这就是全部输入数据了。后面所有的圆、所有的拟合直线都是从这四个数推出来的。1.2 莫尔圆和强度包线究竟是什么关系一个试样破坏时破坏面上的正应力和剪应力并不是任意值而是落在莫尔圆上。莫尔圆的圆心在σ轴上横坐标是p (σ1f σ3) / 2半径是q (σ1f - σ3) / 2。在τ-σ坐标系里画这个圆圆上每一个点就代表该试样在某一个方向平面上的应力状态圆顶点的含义就是最大剪应力对应的面也就是与大主应力面夹角45°的平面。摩尔-库仑强度理论的核心假设是土体破坏时强度包线是τ c σ·tanφ这条直线其中c是粘聚力φ是内摩擦角。当某个莫尔圆恰好和这条包线相切时土体就处于极限平衡状态。所以数据处理的目标很明确找到一条直线使之与所有试验得到的莫尔圆都“尽可能相切”然后读出这条直线的截距和倾角就是c和φ。1.3 别再用直尺贴公切线了手画公切线的经典流程是把所有莫尔圆画在同一张坐标纸上用直尺在圆族外侧找一条能与所有圆相切的直线然后用量角器量倾角读纵截距。这个方法最大的问题是主观性太强。实际试验数据很少能完美落在同一条切线上总有一个圆稍微凸出一点另一个圆稍微缩进去一点直尺到底贴哪一个不同人贴法不同结果也就不同。尤其当某个圆的数据质量不佳时人眼容易被异常点带偏。最小二乘法做的事情本质上就是用“残差平方和最小”这个明确标准替代人眼的直觉判断。拟合结果不受操作者影响同一份数据谁跑都一样还能顺手给出拟合优度指标用来判断整组试验数据靠不靠谱。这就是我推荐用MATLAB跑最小二乘而不是手工贴线的原因。2. 拟合方法选型把“找公切线”变成“拟合一条直线”2.1 关键的坐标变换从莫尔圆到q-p平面直接对莫尔圆族求公切线在数学上是个非线性问题新手很容易卡在这里。工程上更常用的做法是绕一步不直接拟合包线而是对每个莫尔圆的圆心横坐标p和半径q做线性回归拟合这条线叫Kf线在q-p坐标系里画出来就是一条直线。这条线贯穿各个莫尔圆的顶点所以有些教材也叫它“顶点线”或“强度线”。为什么可以绕这一步回到相切条件。设强度包线为τ c σ·tanφ任意一个莫尔圆圆心在(p, 0)半径为q。圆心到包线的距离等于半径时莫尔圆与包线相切。把点到直线的距离公式代进去整理会得到一个极其干净的线性关系q c·cosφ p·sinφ也就是说在q-p平面上各莫尔圆的半径q和圆心横坐标p之间本身就呈线性关系斜率是sinφ截距是c·cosφ。所以我们对(p, q)数据点做一元线性回归再把斜率、截距反算回φ和c就等价于找到了一条与所有圆“相切程度最优”的强度包线。2.2 斜率、截距怎么换算成φ和c设一元回归得到的直线是q a k·p对照上面的公式k sinφa c·cosφ于是φ arcsin(k)c a / cosφ举我前面那组数据为例。算出的4个莫尔圆圆心横坐标和半径分别是试样p (σ1fσ3)/2 (kPa)q (σ1f-σ3)/2 (kPa)1192.5142.52327.5227.53462.5312.54595.0395.0对这4个点做最小二乘线性拟合斜率和截距大约是k≈0.662、a≈13.6。然后换算φ arcsin(0.662) ≈ 41.4°cos(41.4°) ≈ 0.75所以c 13.6 / 0.75 ≈ 18.1 kPa这是一组典型的“有一定粘聚力、内摩擦角较大”的土体参数。如果斜率偏小φ就小说明土的摩擦强度低如果截距偏大c就大说明土体本身具备较好的“黏性”强度。这里有个细节值得强调三角函数计算时MATLAB默认角度单位都是弧度。asin返回的是弧度必须乘以180/pi才是度数而算c a / cos(phi_rad)时必须用cos(弧度值)不要先转成角度再算cos那样会得到完全错误的结果。我第一次写这个脚本就在这个细节上翻过车。2.3 什么时候需要考虑更高级的拟合方法普通最小二乘对三轴试验数据基本够用但有一个前提各级围压下的数据误差相对均匀。如果你发现低围压的数据点明显更离散、高围压点更集中可以考虑给每个点分配不同权重做加权最小二乘如果有个别异常点严重偏离趋势也可以用稳健回归减少它的影响。不过这两种方法对常规课程设计和大多数科研项目来说属于“升级选项”不是必选项。我自己的习惯是先用普通最小二乘出结果然后看R²和残差如果R²掉到0.95以下再回头检查是不是有数据点本身有问题而不是急着换高档算法。3. MATLAB完整实现从试验数据到φ、c和包线图3.1 数据准备与起步代码先把输入数据准备好。最直接的方式是在代码开头用数组把围压和主应力差写死数据量不多的时候这样最清晰。如果数据量很大建议放到Excel里用readmatrix读入但核心计算逻辑完全一样。% 三轴试验原始数据单位统一为 kPa sigma3 [100; 200; 300; 400]; % 各级围压 dsigma [285; 455; 625; 790]; % 破坏时主应力差 % 计算破坏大主应力、莫尔圆圆心横坐标p与半径q sigma1 sigma3 dsigma; p (sigma1 sigma3) / 2; % 圆心横坐标 q (sigma1 - sigma3) / 2; % 莫尔圆半径注意单位必须统一。如果你实验室仪器读数是MPa就全部用MPa是kPa就全部用kPa中途尽量不要混着写否则算出来的c和φ数值会很离谱。3.2 最小二乘拟合与指标换算核心代码最小二乘线性拟合最省事的实现是polyfit。polyfit(x, y, 1)返回的是降幂排列的系数向量第一个元素是斜率第二个元素是截距coef polyfit(p, q, 1); k coef(1); % 斜率即 sin(phi) a coef(2); % 截距即 c * cos(phi) % 换算抗剪强度指标 phi_rad asin(k); % 内摩擦角弧度 phi_deg phi_rad * 180 / pi; % 内摩擦角角度 c a / cos(phi_rad); % 粘聚力 % 拟合优度 R^2 q_fit polyval(coef, p); res q - q_fit; SS_res sum(res.^2); SS_tot sum((q - mean(q)).^2); R2 1 - SS_res / SS_tot; fprintf(内摩擦角 phi %.2f°\n, phi_deg); fprintf(粘聚力 c %.2f kPa\n, c); fprintf(拟合优度 R^2 %.4f\n, R2);这里有几个点值得展开说一下polyfit本身不需要统计工具箱属于MATLAB基础函数所以哪怕你用的是学校机房精简版也没关系。R²的计算方法就是统计学里那个标准定义1减去残差平方和与总平方和的比值。R²越接近1说明(p, q)数据点越接近一条直线也就是各莫尔圆的尺寸和位置关系越符合摩尔-库仑直线包线的假设。如果数据量少比如只有2个围压点R²必然等于1但这不代表结果可靠只是两点决定一条直线的必然结果。至少3到4级围压的R²才有评价意义。3.3 绘制莫尔圆、强度包线与切点标注绘图是整个流程里比较能体现细节的部分。首先绘制莫尔圆推荐用参数方程一个圆心在(p, 0)、半径为q的圆参数方程可以写成σ p q·cos(2θ)τ q·sin(2θ)其中θ对应实际物理面的法线与大主应力方向的夹角。θ从0取到π画的就是上半圆想画完整圆把θ改成0到2π即可。教科书里强度包线讨论的通常都是上半平面所以我习惯只画上半圆图面也干净。figure(Color,w,Position,[100 100 800 600]); hold on; axis equal; grid on; xlabel(正应力 \sigma (kPa)); ylabel(剪应力 \tau (kPa)); title(sprintf(三轴试验莫尔圆与强度包线 \\phi%.1f°, c%.1f kPa, phi_deg, c)); colors lines(length(p)); % 为每个圆分配不同颜色 for i 1:length(p) theta linspace(0, pi, 200); sigma_circle p(i) q(i) * cos(2 * theta); tau_circle q(i) * sin(2 * theta); plot(sigma_circle, tau_circle, Color, colors(i,:), LineWidth, 1.5); % 标记莫尔圆顶点 plot(p(i), q(i), o, MarkerFaceColor, colors(i,:), MarkerEdgeColor, k); end % 绘制真实强度包线 tau c sigma * tan(phi) sigma_line linspace(0, max(p) max(q), 200); tau_line c sigma_line * tan(phi_rad); plot(sigma_line, tau_line, k-, LineWidth, 2.5);这段代码有3个细节容易踩坑第一axis equal必须加。不加的话图形窗口会根据σ和τ的取值范围自动拉伸坐标轴明明应该是圆形的莫尔圆会被画成椭圆。这个坑我见很多同学踩过圆一变形包线和圆是否相切就完全看不出来了。第二绘制包线用的是tan(phi_rad)角度单位必须是弧度。如果你存了phi_deg就直接写tan(phi_deg)算出来的斜率完全不对画出来的包线会飘到天上去。第三不要在莫尔圆图上直接把Kf线画出来。Kf线是在q-p坐标系里的直线坐标轴是p和qσ-τ坐标系里画莫尔圆时如果把Kf线原样叠上去它的位置并不对应强度包线很容易给读者或者审稿人造成误解。想展示拟合效果可以单独画一张q-p散点加拟合直线的图在主图上只画真实包线τ c σ·tanφ就够了。另外如果想标注强度包线与每个莫尔圆的切点公式也不复杂。切点位于上半圆对应中心角2θ 90° φ所以切点坐标是σ* p - q·sinφ τ* q·cosφ把这个点画在每个莫尔圆上可以很直观地看到包线是否真的与圆相切% 计算并标记切点位置 for i 1:length(p) sigma_t p(i) - q(i) * sin(phi_rad); tau_t q(i) * cos(phi_rad); plot(sigma_t, tau_t, rx, MarkerSize, 8, LineWidth, 1.5); end如果拟合效果好这些叉号会正好落在黑色包线上如果某个点偏离明显说明该级围压的试验数据可能有异常值得回去查原始记录。3.4 输出结果说明与拟合效果判断运行完整脚本后命令行会输出类似这样的结果内摩擦角 phi 41.42° 粘聚力 c 18.10 kPa 拟合优度 R^2 0.9993R²达到0.9993说明这组数据线性程度很好σ3取100到400 kPa范围内摩尔-库仑直线包线的假设是合适的。如果R²偏低要分情况看待一是数据本身波动大二是强度包线在这个应力范围内本身就明显弯曲。对于后者比如某些超固结黏土在低围压段的包线呈现明显曲率这时候强行用直线拟合虽然也能出数但要把适用应力范围在报告里写清楚不能外推到围压范围之外。4. 实际操作中的常见坑与排查建议4.1 莫尔圆被画成“椭圆”八成是坐标轴比例问题sigma轴和tau轴虽然单位都是kPa但取值范围不同sigma从0到1000多tau最大也就400多。MATLAB默认会按数据范围自动缩放坐标轴两个轴的单位长度不一致圆立刻变扁或者变瘦。解决方式就是前面代码里那句axis equal。加了之后两个轴按相同比例缩放圆才是真正的圆。这个坑非常好避免但几乎每次有人把图发给我看第一眼就能看到圆变形——所以无论如何画莫尔圆之前先把这句话写上。4.2 拟合斜率超过1怎么办理论上k sinφ打死也不可能大于1。但试验数据是有噪声的当某个围压级别数据偏差太大时最小二乘拟合出来的斜率完全可能变成1.02甚至更大。这时候asin(k)会返回复数或者直接报错。遇到这种情况不要慌先按以下顺序排查检查数据录入是否有误比如把395错录成935。检查该级围压试样的破坏模式。如果试样出现明显端部约束或者先沿着某个软弱面破坏数据点会系统性偏离趋势线。检查是不是不同试样之间土体不均匀。三轴试验要求同一组试样初始状态尽量一致如果孔隙比相差很大强度自然对不上。如果所有检查都找不到明确错误可以作一个“手工约束版”的拟合强制斜率不大于0.99同时输出警告。但说实话这种情况更可能的结论是这组试验不宜拿来求抗剪强度指标重做试样可能比强行拟合更有意义。4.3 总应力指标与有效应力指标千万别混用UU、CU、CD三类三轴试验得到的c、φ含义完全不同。UU试验得到的是不排水总应力指标适用于模拟饱和黏性土快速加载的短期稳定问题CU试验如果不同步测孔压得到的也是总应力指标通常记为ccu和φcuCD试验以及CU试验中同步测定孔压后换算出的结果是有效应力指标记为c和φ适用于排水长期工况。画强度包线时同一张图里要么全是总应力莫尔圆要么全是有效应力莫尔圆绝对不能混用。我见过一份课程报告低围压用有效应力圆、高围压用总应力圆拟合出来的φ和c完全不知道属于什么工况这种数据拿到答辩现场基本就是送分题给老师挑毛病。4.4 数据点偏少、围压范围偏窄时怎么处理如果只有两个围压点比如σ3 100和300 kPa那么两个点一定能连出一条直线φ和c看起来也有模有样但这只是被两个点“硬锁”出来的结果没有任何冗余来验证合理性。常规三轴试验都要求至少3到4级围压并且围压范围要覆盖你关心的应力区间。整理报告时如果条件确实只允许做2个有效试样建议在结论部分明确写明“本次仅由2级围压数据拟合强度指标供参考”避免给后续使用造成误导。4.5 数据录入与单位引发的怪异结果有一种很隐蔽的错误是把某组主应力差输成主应力。比如把σ1f 385输入到dsigma的位置相当于把莫尔圆半径翻了一倍拟合出来的φ会异常偏高c甚至会算出负值。碰到φ超过50°或者c出现负数的“可疑结果”第一时间核原始数据不要急着怀疑算法。我在处理同学的数据时至少有一半的“异常结果”最后都出在录入错误上。5. 批量处理和扩展方向5.1 把拟合过程封装成函数多组土样一次跑完如果你要对多组土样分别求c和φ最清晰的组织方式是把核心计算封装成一个函数输入sigma3和dsigma输出phi_deg、c和R2。这样主脚本只需要循环调用函数读不同的数据文件最后汇总结果到一个表格。function [phi_deg, c, R2] triaxial_fit(sigma3, dsigma) sigma1 sigma3 dsigma; p (sigma1 sigma3) / 2; q (sigma1 - sigma3) / 2; coef polyfit(p, q, 1); k coef(1); a coef(2); phi_rad asin(k); phi_deg phi_rad * 180 / pi; c a / cos(phi_rad); q_fit polyval(coef, p); SS_res sum((q - q_fit).^2); SS_tot sum((q - mean(q)).^2); R2 1 - SS_res / SS_tot; end调用的时候先读Excel再循环data readmatrix(试验数据.xlsx); % 每一列依次是围压和主应力差 for g 1:size(data, 3) % 按实际情况调整维度 [phi(g), c(g), R2(g)] triaxial_fit(data(:,1,g), data(:,2,g)); end这样几十组土样也能几分钟内出齐表格比手动画图快一个量级。5.2 从总应力指标扩展到有效应力指标很多三轴试验课程项目做到CU试验时会记录孔隙水压力u。这时候可以分别用总应力(σ1, σ3)和有效应力(σ1 σ1 - u, σ3 σ3 - u)各算一组c、φ。两者的物理含义完全不同在报告里都值得给出来。计算有效应力指标时只要替换sigma3和dsigma的输入逻辑把孔压扣掉就行核心拟合代码完全不用改。5.3 非线性强度包线的处理思路如果你的土样在较宽的围压范围内数据显示包线明显弯曲单纯直线拟合就不够了。一个常用的折中做法是把围压范围分段每段分别做直线拟合分别读取c、φ另一个做法是改用非线性准则描述这就超出本文最小二乘直线拟合的范畴了。对于大多数课程设计和常规岩土工程问题摩尔-库仑直线包线仍然是通行做法重点是把拟合的适用应力范围交代清楚。最后分享一个我自己的小习惯每次拟合完我都会顺手把每个莫尔圆圆心到包线的垂直距离打印出来和半径放在一起对比。垂直距离由公式d (p·tanφ c) / sqrt(tan²φ 1)算出理论上应该等于半径q。如果两者差值都在0.5 kPa以内说明这组拟合线的“相切程度”是可信的一旦某个圆明显偏离一定是那级试验数据本身有情况。这个检查动作花不了两分钟但能让你的结论在答辩或报告审查时经得起追问我自己已经把它列进每次数据处理的标准流程了。
MATLAB三轴试验强度包线拟合：最小二乘法求φ和c完整教程
用MATLAB做三轴试验强度包线最小二乘法拟合φ和c的完整套路做三轴剪切试验写报告或者做毕设的同学到了数据处理这一关十有八九会被“画莫尔圆、作公切线、读内摩擦角φ和粘聚力c”这套流程折磨过。早几年我自己也是拿直尺和三角板在坐标纸上一条条贴公切线贴完还要用两把尺子量角度费时不说不同人贴出来的切线位置差个一两度都是常事——而这一两度换算到φ和c上面结论可能就差出一大截。后来我彻底改成用MATLAB处理输入每级围压和破坏时主应力差用最小二乘法自动拟合强度包线φ和c直接出数值图也顺手画好省下来的时间用来复查数据质量比手工画图划算得多。这篇就把整套思路和代码完整写出来讲清楚每一行在干什么、为什么这么算以及哪些坑我踩过之后希望你直接绕开。先说适用对象手里有三轴试验数据、需要按摩尔-库仑准则求总应力或有效应力抗剪强度指标内摩擦角φ、粘聚力c的本科生、研究生和工程师。代码只用到MATLAB最基础的polyfit和plot函数不依赖任何额外工具箱R2016b之后的版本都能跑。即使你MATLAB还不熟按下面的步骤抄也能出图出数。1. 先搞懂原理三轴试验、莫尔圆和强度包线怎么串起来1.1 三轴试验到底能给我们哪些数常规三轴压缩试验的做法是把3到4个相同土样分别放在不同的恒定围压σ3下然后施加轴向压力让试样剪切破坏。围压就是压力室里的水压一般取100、200、300、400 kPa这一档。每个试样剪切过程中仪器记录轴向应变和主应力差σ1-σ3的关系曲线试验结束时从曲线上取破坏点对应的主应力差σ1-σ3f——如果曲线有峰值就取峰值没有峰值就按15%轴向应变对应的值取。有了围压σ3和破坏时主应力差σ1-σ3f就能算出试样破坏时的大主应力σ1f σ3 (σ1-σ3)f一组围压对应一个试样一个试样就能画出一个莫尔圆。比如我用过的典型数据试样编号围压σ3 (kPa)破坏主应力差 (σ1-σ3)f (kPa)破坏大主应力σ1f (kPa)11002853852200455655330062592544007901190这就是全部输入数据了。后面所有的圆、所有的拟合直线都是从这四个数推出来的。1.2 莫尔圆和强度包线究竟是什么关系一个试样破坏时破坏面上的正应力和剪应力并不是任意值而是落在莫尔圆上。莫尔圆的圆心在σ轴上横坐标是p (σ1f σ3) / 2半径是q (σ1f - σ3) / 2。在τ-σ坐标系里画这个圆圆上每一个点就代表该试样在某一个方向平面上的应力状态圆顶点的含义就是最大剪应力对应的面也就是与大主应力面夹角45°的平面。摩尔-库仑强度理论的核心假设是土体破坏时强度包线是τ c σ·tanφ这条直线其中c是粘聚力φ是内摩擦角。当某个莫尔圆恰好和这条包线相切时土体就处于极限平衡状态。所以数据处理的目标很明确找到一条直线使之与所有试验得到的莫尔圆都“尽可能相切”然后读出这条直线的截距和倾角就是c和φ。1.3 别再用直尺贴公切线了手画公切线的经典流程是把所有莫尔圆画在同一张坐标纸上用直尺在圆族外侧找一条能与所有圆相切的直线然后用量角器量倾角读纵截距。这个方法最大的问题是主观性太强。实际试验数据很少能完美落在同一条切线上总有一个圆稍微凸出一点另一个圆稍微缩进去一点直尺到底贴哪一个不同人贴法不同结果也就不同。尤其当某个圆的数据质量不佳时人眼容易被异常点带偏。最小二乘法做的事情本质上就是用“残差平方和最小”这个明确标准替代人眼的直觉判断。拟合结果不受操作者影响同一份数据谁跑都一样还能顺手给出拟合优度指标用来判断整组试验数据靠不靠谱。这就是我推荐用MATLAB跑最小二乘而不是手工贴线的原因。2. 拟合方法选型把“找公切线”变成“拟合一条直线”2.1 关键的坐标变换从莫尔圆到q-p平面直接对莫尔圆族求公切线在数学上是个非线性问题新手很容易卡在这里。工程上更常用的做法是绕一步不直接拟合包线而是对每个莫尔圆的圆心横坐标p和半径q做线性回归拟合这条线叫Kf线在q-p坐标系里画出来就是一条直线。这条线贯穿各个莫尔圆的顶点所以有些教材也叫它“顶点线”或“强度线”。为什么可以绕这一步回到相切条件。设强度包线为τ c σ·tanφ任意一个莫尔圆圆心在(p, 0)半径为q。圆心到包线的距离等于半径时莫尔圆与包线相切。把点到直线的距离公式代进去整理会得到一个极其干净的线性关系q c·cosφ p·sinφ也就是说在q-p平面上各莫尔圆的半径q和圆心横坐标p之间本身就呈线性关系斜率是sinφ截距是c·cosφ。所以我们对(p, q)数据点做一元线性回归再把斜率、截距反算回φ和c就等价于找到了一条与所有圆“相切程度最优”的强度包线。2.2 斜率、截距怎么换算成φ和c设一元回归得到的直线是q a k·p对照上面的公式k sinφa c·cosφ于是φ arcsin(k)c a / cosφ举我前面那组数据为例。算出的4个莫尔圆圆心横坐标和半径分别是试样p (σ1fσ3)/2 (kPa)q (σ1f-σ3)/2 (kPa)1192.5142.52327.5227.53462.5312.54595.0395.0对这4个点做最小二乘线性拟合斜率和截距大约是k≈0.662、a≈13.6。然后换算φ arcsin(0.662) ≈ 41.4°cos(41.4°) ≈ 0.75所以c 13.6 / 0.75 ≈ 18.1 kPa这是一组典型的“有一定粘聚力、内摩擦角较大”的土体参数。如果斜率偏小φ就小说明土的摩擦强度低如果截距偏大c就大说明土体本身具备较好的“黏性”强度。这里有个细节值得强调三角函数计算时MATLAB默认角度单位都是弧度。asin返回的是弧度必须乘以180/pi才是度数而算c a / cos(phi_rad)时必须用cos(弧度值)不要先转成角度再算cos那样会得到完全错误的结果。我第一次写这个脚本就在这个细节上翻过车。2.3 什么时候需要考虑更高级的拟合方法普通最小二乘对三轴试验数据基本够用但有一个前提各级围压下的数据误差相对均匀。如果你发现低围压的数据点明显更离散、高围压点更集中可以考虑给每个点分配不同权重做加权最小二乘如果有个别异常点严重偏离趋势也可以用稳健回归减少它的影响。不过这两种方法对常规课程设计和大多数科研项目来说属于“升级选项”不是必选项。我自己的习惯是先用普通最小二乘出结果然后看R²和残差如果R²掉到0.95以下再回头检查是不是有数据点本身有问题而不是急着换高档算法。3. MATLAB完整实现从试验数据到φ、c和包线图3.1 数据准备与起步代码先把输入数据准备好。最直接的方式是在代码开头用数组把围压和主应力差写死数据量不多的时候这样最清晰。如果数据量很大建议放到Excel里用readmatrix读入但核心计算逻辑完全一样。% 三轴试验原始数据单位统一为 kPa sigma3 [100; 200; 300; 400]; % 各级围压 dsigma [285; 455; 625; 790]; % 破坏时主应力差 % 计算破坏大主应力、莫尔圆圆心横坐标p与半径q sigma1 sigma3 dsigma; p (sigma1 sigma3) / 2; % 圆心横坐标 q (sigma1 - sigma3) / 2; % 莫尔圆半径注意单位必须统一。如果你实验室仪器读数是MPa就全部用MPa是kPa就全部用kPa中途尽量不要混着写否则算出来的c和φ数值会很离谱。3.2 最小二乘拟合与指标换算核心代码最小二乘线性拟合最省事的实现是polyfit。polyfit(x, y, 1)返回的是降幂排列的系数向量第一个元素是斜率第二个元素是截距coef polyfit(p, q, 1); k coef(1); % 斜率即 sin(phi) a coef(2); % 截距即 c * cos(phi) % 换算抗剪强度指标 phi_rad asin(k); % 内摩擦角弧度 phi_deg phi_rad * 180 / pi; % 内摩擦角角度 c a / cos(phi_rad); % 粘聚力 % 拟合优度 R^2 q_fit polyval(coef, p); res q - q_fit; SS_res sum(res.^2); SS_tot sum((q - mean(q)).^2); R2 1 - SS_res / SS_tot; fprintf(内摩擦角 phi %.2f°\n, phi_deg); fprintf(粘聚力 c %.2f kPa\n, c); fprintf(拟合优度 R^2 %.4f\n, R2);这里有几个点值得展开说一下polyfit本身不需要统计工具箱属于MATLAB基础函数所以哪怕你用的是学校机房精简版也没关系。R²的计算方法就是统计学里那个标准定义1减去残差平方和与总平方和的比值。R²越接近1说明(p, q)数据点越接近一条直线也就是各莫尔圆的尺寸和位置关系越符合摩尔-库仑直线包线的假设。如果数据量少比如只有2个围压点R²必然等于1但这不代表结果可靠只是两点决定一条直线的必然结果。至少3到4级围压的R²才有评价意义。3.3 绘制莫尔圆、强度包线与切点标注绘图是整个流程里比较能体现细节的部分。首先绘制莫尔圆推荐用参数方程一个圆心在(p, 0)、半径为q的圆参数方程可以写成σ p q·cos(2θ)τ q·sin(2θ)其中θ对应实际物理面的法线与大主应力方向的夹角。θ从0取到π画的就是上半圆想画完整圆把θ改成0到2π即可。教科书里强度包线讨论的通常都是上半平面所以我习惯只画上半圆图面也干净。figure(Color,w,Position,[100 100 800 600]); hold on; axis equal; grid on; xlabel(正应力 \sigma (kPa)); ylabel(剪应力 \tau (kPa)); title(sprintf(三轴试验莫尔圆与强度包线 \\phi%.1f°, c%.1f kPa, phi_deg, c)); colors lines(length(p)); % 为每个圆分配不同颜色 for i 1:length(p) theta linspace(0, pi, 200); sigma_circle p(i) q(i) * cos(2 * theta); tau_circle q(i) * sin(2 * theta); plot(sigma_circle, tau_circle, Color, colors(i,:), LineWidth, 1.5); % 标记莫尔圆顶点 plot(p(i), q(i), o, MarkerFaceColor, colors(i,:), MarkerEdgeColor, k); end % 绘制真实强度包线 tau c sigma * tan(phi) sigma_line linspace(0, max(p) max(q), 200); tau_line c sigma_line * tan(phi_rad); plot(sigma_line, tau_line, k-, LineWidth, 2.5);这段代码有3个细节容易踩坑第一axis equal必须加。不加的话图形窗口会根据σ和τ的取值范围自动拉伸坐标轴明明应该是圆形的莫尔圆会被画成椭圆。这个坑我见很多同学踩过圆一变形包线和圆是否相切就完全看不出来了。第二绘制包线用的是tan(phi_rad)角度单位必须是弧度。如果你存了phi_deg就直接写tan(phi_deg)算出来的斜率完全不对画出来的包线会飘到天上去。第三不要在莫尔圆图上直接把Kf线画出来。Kf线是在q-p坐标系里的直线坐标轴是p和qσ-τ坐标系里画莫尔圆时如果把Kf线原样叠上去它的位置并不对应强度包线很容易给读者或者审稿人造成误解。想展示拟合效果可以单独画一张q-p散点加拟合直线的图在主图上只画真实包线τ c σ·tanφ就够了。另外如果想标注强度包线与每个莫尔圆的切点公式也不复杂。切点位于上半圆对应中心角2θ 90° φ所以切点坐标是σ* p - q·sinφ τ* q·cosφ把这个点画在每个莫尔圆上可以很直观地看到包线是否真的与圆相切% 计算并标记切点位置 for i 1:length(p) sigma_t p(i) - q(i) * sin(phi_rad); tau_t q(i) * cos(phi_rad); plot(sigma_t, tau_t, rx, MarkerSize, 8, LineWidth, 1.5); end如果拟合效果好这些叉号会正好落在黑色包线上如果某个点偏离明显说明该级围压的试验数据可能有异常值得回去查原始记录。3.4 输出结果说明与拟合效果判断运行完整脚本后命令行会输出类似这样的结果内摩擦角 phi 41.42° 粘聚力 c 18.10 kPa 拟合优度 R^2 0.9993R²达到0.9993说明这组数据线性程度很好σ3取100到400 kPa范围内摩尔-库仑直线包线的假设是合适的。如果R²偏低要分情况看待一是数据本身波动大二是强度包线在这个应力范围内本身就明显弯曲。对于后者比如某些超固结黏土在低围压段的包线呈现明显曲率这时候强行用直线拟合虽然也能出数但要把适用应力范围在报告里写清楚不能外推到围压范围之外。4. 实际操作中的常见坑与排查建议4.1 莫尔圆被画成“椭圆”八成是坐标轴比例问题sigma轴和tau轴虽然单位都是kPa但取值范围不同sigma从0到1000多tau最大也就400多。MATLAB默认会按数据范围自动缩放坐标轴两个轴的单位长度不一致圆立刻变扁或者变瘦。解决方式就是前面代码里那句axis equal。加了之后两个轴按相同比例缩放圆才是真正的圆。这个坑非常好避免但几乎每次有人把图发给我看第一眼就能看到圆变形——所以无论如何画莫尔圆之前先把这句话写上。4.2 拟合斜率超过1怎么办理论上k sinφ打死也不可能大于1。但试验数据是有噪声的当某个围压级别数据偏差太大时最小二乘拟合出来的斜率完全可能变成1.02甚至更大。这时候asin(k)会返回复数或者直接报错。遇到这种情况不要慌先按以下顺序排查检查数据录入是否有误比如把395错录成935。检查该级围压试样的破坏模式。如果试样出现明显端部约束或者先沿着某个软弱面破坏数据点会系统性偏离趋势线。检查是不是不同试样之间土体不均匀。三轴试验要求同一组试样初始状态尽量一致如果孔隙比相差很大强度自然对不上。如果所有检查都找不到明确错误可以作一个“手工约束版”的拟合强制斜率不大于0.99同时输出警告。但说实话这种情况更可能的结论是这组试验不宜拿来求抗剪强度指标重做试样可能比强行拟合更有意义。4.3 总应力指标与有效应力指标千万别混用UU、CU、CD三类三轴试验得到的c、φ含义完全不同。UU试验得到的是不排水总应力指标适用于模拟饱和黏性土快速加载的短期稳定问题CU试验如果不同步测孔压得到的也是总应力指标通常记为ccu和φcuCD试验以及CU试验中同步测定孔压后换算出的结果是有效应力指标记为c和φ适用于排水长期工况。画强度包线时同一张图里要么全是总应力莫尔圆要么全是有效应力莫尔圆绝对不能混用。我见过一份课程报告低围压用有效应力圆、高围压用总应力圆拟合出来的φ和c完全不知道属于什么工况这种数据拿到答辩现场基本就是送分题给老师挑毛病。4.4 数据点偏少、围压范围偏窄时怎么处理如果只有两个围压点比如σ3 100和300 kPa那么两个点一定能连出一条直线φ和c看起来也有模有样但这只是被两个点“硬锁”出来的结果没有任何冗余来验证合理性。常规三轴试验都要求至少3到4级围压并且围压范围要覆盖你关心的应力区间。整理报告时如果条件确实只允许做2个有效试样建议在结论部分明确写明“本次仅由2级围压数据拟合强度指标供参考”避免给后续使用造成误导。4.5 数据录入与单位引发的怪异结果有一种很隐蔽的错误是把某组主应力差输成主应力。比如把σ1f 385输入到dsigma的位置相当于把莫尔圆半径翻了一倍拟合出来的φ会异常偏高c甚至会算出负值。碰到φ超过50°或者c出现负数的“可疑结果”第一时间核原始数据不要急着怀疑算法。我在处理同学的数据时至少有一半的“异常结果”最后都出在录入错误上。5. 批量处理和扩展方向5.1 把拟合过程封装成函数多组土样一次跑完如果你要对多组土样分别求c和φ最清晰的组织方式是把核心计算封装成一个函数输入sigma3和dsigma输出phi_deg、c和R2。这样主脚本只需要循环调用函数读不同的数据文件最后汇总结果到一个表格。function [phi_deg, c, R2] triaxial_fit(sigma3, dsigma) sigma1 sigma3 dsigma; p (sigma1 sigma3) / 2; q (sigma1 - sigma3) / 2; coef polyfit(p, q, 1); k coef(1); a coef(2); phi_rad asin(k); phi_deg phi_rad * 180 / pi; c a / cos(phi_rad); q_fit polyval(coef, p); SS_res sum((q - q_fit).^2); SS_tot sum((q - mean(q)).^2); R2 1 - SS_res / SS_tot; end调用的时候先读Excel再循环data readmatrix(试验数据.xlsx); % 每一列依次是围压和主应力差 for g 1:size(data, 3) % 按实际情况调整维度 [phi(g), c(g), R2(g)] triaxial_fit(data(:,1,g), data(:,2,g)); end这样几十组土样也能几分钟内出齐表格比手动画图快一个量级。5.2 从总应力指标扩展到有效应力指标很多三轴试验课程项目做到CU试验时会记录孔隙水压力u。这时候可以分别用总应力(σ1, σ3)和有效应力(σ1 σ1 - u, σ3 σ3 - u)各算一组c、φ。两者的物理含义完全不同在报告里都值得给出来。计算有效应力指标时只要替换sigma3和dsigma的输入逻辑把孔压扣掉就行核心拟合代码完全不用改。5.3 非线性强度包线的处理思路如果你的土样在较宽的围压范围内数据显示包线明显弯曲单纯直线拟合就不够了。一个常用的折中做法是把围压范围分段每段分别做直线拟合分别读取c、φ另一个做法是改用非线性准则描述这就超出本文最小二乘直线拟合的范畴了。对于大多数课程设计和常规岩土工程问题摩尔-库仑直线包线仍然是通行做法重点是把拟合的适用应力范围交代清楚。最后分享一个我自己的小习惯每次拟合完我都会顺手把每个莫尔圆圆心到包线的垂直距离打印出来和半径放在一起对比。垂直距离由公式d (p·tanφ c) / sqrt(tan²φ 1)算出理论上应该等于半径q。如果两者差值都在0.5 kPa以内说明这组拟合线的“相切程度”是可信的一旦某个圆明显偏离一定是那级试验数据本身有情况。这个检查动作花不了两分钟但能让你的结论在答辩或报告审查时经得起追问我自己已经把它列进每次数据处理的标准流程了。
