土力学是土建类专业的一门本科主干专业课程，通过土力学的学习，学生应牢固掌握土力学的基本原理、计算分析方法和基本实验技能，初步具备分析和解决相关工程问题的能力。然而，不同于其他力学类课程，土力学缺乏系统性，知识点多且分散，导致学生难以把握其内在知识脉络，学习效果大打折扣，解决实际问题的能力严重不足。

一、土力学教学现状及存在问题

当前的土力学教学，主要存在以下几个问题亟待解决。

(1) 教学重点把握不合理，实用性欠缺，学生解决实际问题的能力得不到锻炼。培养学生解决实际工程问题的能力是土力学教学重要目标之一，实际教学过程中，由于公式繁琐、计算量大，一些重要的工程问题在课堂教学中不得不被简化乃至放弃。比如，粘性土边坡稳定是工程实践中的常见问题，一般采用条分法计算，但是其计算过程极为繁琐，学生很难通过实例求解来掌握。由于缺少真正意义上的工程训练，学生通常在学习土力学过程中感到枯燥、空洞，在工作中遇到岩土工程问题时无从下手。

(2) 教学方法陈旧，教学效果欠佳，学生学习积极性不高。土力学是一门古老的学科，教材知识点相对陈旧，传统的教学模式完全忠于教材，学生只会机械地套用公式，缺乏创新性，不能灵活运用各种手段解决实际问题。比如，沉降计算是土力学教学中的重要内容，工程中经常使用分层总和法进行沉降计算，其计算原理是先将地基土分为若干层，计算每层土的压缩量，然后累积起来即为总的地基沉降量。该方法在计算过程中需反复查表，计算多点应力和多层土压缩量，手工计算效率非常低，而在当前课堂教学中，学生仍然采取手工计算方式求解。限于课时，又不得不将实际工程问题一再简化。长此以往，学生容易养成机械的学习习惯，对土力学的基本概念和原理一知半解。

(3) 工程求解手段落后，不利于学生创新能力培养。土力学教学内容中有大量对经验的成果总结，这些知识对于缺乏理论分析能力和工程实践经验的本科生来说，难以真正理解。比如一些参数无法由理论推导得出，需要根据作图法确定，先期固结压力就是其一。教材中一般采用Cassagrande作图法求取，该方法虽然简单、易行，但其是一种经验作图法，有很大局限性。该方法求得的结果很大程度上取决于作图精度，准确性难以得到保证。类似现象还有许多，使学生在土力学的学习中产生困惑，难以理解土力学作为一门专业基础课的科学性、严谨性，更不利于学生创新能力的培养。

近年来，计算机辅助教学技术的普及和数学计算软件的快速发展，为土力学教学中改善以上问题提供了新思路。MATLAB即是众多数学计算软件中的杰出代表，是美国Mathworks公司开发的主要面对科学计算、可视化以及交互式程序设计的高科技计算软件，其将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为众多领域的科学研究和工程设计提供了一种全面的解决方案。由于其强大的功能，MATLAB在土力学理论研究及实际工程中均得到了广泛应用^[1]，受到了科研工作者的青睐，部分高校教师也开始尝试将其引入土力学课堂教学，解决一些复杂的计算问题，取得了明显的教学效果^[2-3]。笔者将在前人研究的基础上，进一步总结MATLAB在土力学教学中的应用，为土力学教学改革提供参考。

二、MATLAB在土力学中的应用举例 (一) 在试验数据处理中的应用 1 在液塑限试验中的应用

规范^[4]规定，液塑限联合测定试验时，以含水率为横坐标，圆锥锥入深度为纵坐标绘制双对数坐标图。当三点共线时，分别取圆锥锥入深度为2 mm、17 mm时对应的含水率为塑限和液限；当三点不共线时，过高含水率的点连接其余两点得到两条直线，分别在两条直线上求取圆锥锥入深度为2 mm时的含水率，若两个含水率差值小于2%，则取两点的平均值与高含水率的点作直线，该直线上锥入深度为2 mm、17 mm的点即为塑限和液限，当所得两个含水率差值大于2%时，则应重做试验。按该方法作图时，不仅步骤繁琐、精度难以保证，而且无法预先判断试验数据是否合格。若编制MATLAB小程序求解该问题，则可达到事半功倍的效果，该程序如下：

$ \rm{x=}\left[\rm{w1, w2, w3} \right]\rm{; y=}\left[\rm{h1, h2, h3} \right]\rm{;} $

%三个点的含水率和锥入深度，按含水率由高到低排序

$ \rm{k12 = ( y( 1)-y( 2) ) /( x( 1) -x( 2) );} $

%过点1、点2直线的斜率

$ \rm{k13 = ( y( 1)-y( 3) ) /( x( 1) -x( 3) );} $

%过点1、点3直线的斜率

if k12 ~= k13

if abs((2-y(1)) /k12-(2-y(1)) /k13) < 2

wp=(2-y(1)) /(2* k12) +(2-y(1)) /(2* k13) +x(1)

k14 = (y(1)-2) /(x(1)-wp);

wl = (17-y(1)) /k14+x(1)

else

title ('试验数据不满足要求，请重做试验! ');

end

else

wp = (2-y(1)) /k12+x(1)

wl = (17-y(1)) /k12+x(1)

end

2 先期固结压力的确定

规范^[4]规定，原状土的先期固结压力，可按Cassagrande法确定：先在e~logp曲线上找出曲率半径最小的一点Q，过该点作水平线QA和切线QB及角AQB的角平分线QD，QD与曲线下半段的直线段延长线交于点E，E点对应的压力即为该土样的先期固结压力。该方法存在绘图繁琐、最小曲率半径点不易确定等缺点，作图误差较大。采用MATLAB程序求解，可以很好地解决这一问题。其中，压缩曲线的曲线段可采用指数形式的三次多项式拟合，具体推导可见文献[5]。该问题主要程序段如下：

a = polyfit(x, y, 3);   %寻找曲率半径最小点Q

xi = 25:5:1000;

yi = polyval(a, log10(xi));

plot(xi, yi, 'k-');

syms e Q;

e = a(1)*Q^3+a(2)*Q^2+a(3)*Q+a(4);

syms de de2 R;

de = diff(e, Q);

de2 = diff(e, Q, 2);

R = abs(de2)/(1+de^2)^(3/2);

xi = 25:1:800;

Ri = subs(R, Q, log10(xi));

R_min = min(Ri);

pos_R_min = find(Ri == R_min);

Q_min = xi(pos_R_min);

e_pos = subs(e, Q, log10(Q_min));

……

x_paral = [Q_min, Q_min + 150];   %过Q点作水平线QA

y_paral = [e_pos, e_pos];

plot(x_paral, y_paral, 'k:');

……

k1 = subs(de, Q, log10(Q_min));   %过Q点做曲线的切线QB

x_tan = [Q_min, Q_min + 150];

y_tan = [e_pos, k1*log10((Q_min+150)/Q_min)+e_pos];

plot(x_tan, y_tan, 'k:');

……

k = -tan(atan(-k1)/2);   %绘制角平分线QD

x_halve = [Q_min, Q_min + 300];

y_halve = [e_pos, k*log10((Q_min+300)/Q_min)+e_pos];

plot(x_halve, y_halve, 'k-');

……

xx = x(end-2:end);   %最小二乘法拟合最后三个点，绘制压缩曲线直线段

yy = y(end-2:end);

aa = polyfit(xx, yy, 1);

logpc = (e_pos-(aa(2)+k*log10(Q_min)))/(aa(1)-k);

xxi = [10.^logpc-100，1000];

yyi = polyval(aa, log10(xxi));

plot(xxi, yyi, 'k-');

……

plot(10.^logpc, k*(logpc-log10(Q_min))+e_pos, 'k.');   %求QD与拟合直线的交点

gtext(['(', num2str(Q_min), ', ', num2str(e_pos)，')']);

gtext(['(', num2str(10.^logpc), ', ', num2str(k*(logpc-log10(Q_min))+e_pos)，')']);

pc = 10^(logpc);

……

(二) 在工程计算中的应用 1 地基附加应力的计算

附加应力计算是地基沉降计算中重要的步骤，教材中通常根据荷载类型选用相应的表格，查取附加应力系数进行计算。在实际查表中，常常需要进行多次插值，耗费大量时间。若采用MATLAB编制程序求解，则非常方便快捷。以矩形均布荷载中心以下一定深度处的附加应力求解为例，给出其求解程序如下：

p=200; l=4; b=2; z=-10:0.1:0;

m=l/b; n=-2*z/b;

sigmaZ=2*p/pi().*(m.*n.*(m^2+2.*n.^2+1)./(m^2+n.^2)./(1+n.^2)./sqrt(m^2+n.^2+1)+asin(m./sqrt((m^2+n.^2).*(1+n.^2))));

subplot(1, 1, 1);

plot(sigmaZ, z);

grid on;

title('sigmaZ随z的关系曲线');

xlabel('sigmaZ(kPa)');

ylabel('z(m)');

以上程序求解的是长4 m、宽2 m，大小为200 kPa的矩形均布荷载中心点下0到10 m的附加应力分布，如图 1所示。

2 滑坡稳定性计算

滑坡稳定性计算是土力学中的重要问题，通常用条分法求解，即先假定一个滑动面，对滑动土体进行条分并求取安全系数，然后逐一搜索其他滑动面并计算安全系数，最后找出最危险的滑动面及其对应的安全系数。该过程计算量非常大，且容易出错，而采用MATLAB程序可方便求解。文章以Fellenius条分法为例，简述其实现过程。

如图 2所示，设均质土坡沿圆弧滑动面AC滑动，圆心为O，半径为R，将滑动土体ABC分为若干土条，忽略土条间的作用力，并假定各土条底部滑动面上的安全系数等于整个滑动面的安全系数，则该滑动面安全系数为

$ \mathit{K}=\frac{\sum{({{\mathit{c}}_{\mathit{i}}}{{\mathit{l}}_{\mathit{i}}}+{{\mathit{G}}_{\mathit{i}}}\rm{cos}{{\mathit{\theta }}_{\mathit{i}}}\rm{tan}{{\mathit{\varphi }}_{\mathit{i}}})}}{\sum{{{\mathit{G}}_{\mathit{i}}}\rm{sin}{{\mathit{\theta }}_{\mathit{i}}}}} $

(1)

式中，c_i和φ_i分别是第i个土条的粘聚力和内摩擦角，G_i为该土条的重力，I_i为该土条范围内的滑弧长，θ_i该土条底面中点的法线与竖直线的夹角。

Fellenius条分法的程序框图如图 3所示，文献[3]给出了一种解法，文章在其基础上对最危险滑弧的搜索方法进行了优化，其求解程序如下：

function [Fsmin, xb, yb, Rb]= Fellenius (b, h, gama, phi0, c)

 Fsmin=100.0;

 phi0=phi0*pi/180.0;

 alpha=atan(b);

 L=h/sin(alpha);

 m=L*cos(alpha);

 for theta=0:alpha/100:(alpha-phi0)

  L0=h/sin(alpha-theta);

  m0=L0*cos(alpha-theta);

  x0=m0/2;

  y0=h/2;

  d0=L0/2*tan(theta);

  for d=d0:(L/1000.0):(1.25*L)

   sum1=0;sum2=0;

   x=x0-cos(pi/2-alpha+theta)*d;

   y=y0+sin(pi/2-alpha+theta)*d;

   R=sqrt(x^2+y^2);

   d1=m0/1000;

   %for xd=0:d1:x1

   for xd=0:d1:m0

    yd=y-sqrt(R^2-(xd-x)^2);

    beta=atan((xd-x)/(y-yd));

    n=d1/cos(beta);

    if xd < =m

     y2=tan(alpha)*xd;

     h1=abs(y2-yd);

     W1=gama*h1*d1;

     sum1=sum1+W1*cos(beta)*tan(phi0)+c*n;

     sum2=sum2+W1*sin(beta);

    else

     h1=abs(h-yd);

     W1=gama*h1*d1;

     sum1=sum1+W1*cos(beta)*tan(phi0)+c*n;

     sum2=sum2+W1*sin(beta);

    end

   end

   Fs=sum1/sum2;

   if(Fs < Fsmin)

    Fsmin=Fs;

    xb=x;

    yb=y;

    Rb=R;

   end

  end

 end

end

以文献[6]中第266页例题10-2为算例，均质黏性土土坡坡高20 m，坡度1：2，填土黏聚力为10 kPa，内摩擦角20°，重度18 kN/m³，利用程序可求得该土坡安全系数为1.091，最危险滑弧圆心坐标为(7.054，40.963)。

三、结语

文章介绍了MATLAB在土力学教学中的应用情况，并给出了教学实践中的部分应用实例。实践表明，在土力学教学中引入MATLAB程序，不仅可以激发学生学习的兴趣，还可加深学生对土力学原理的理解，增强学生对工程实际问题的认识，并提高其力学建模和数学计算能力，为将来更好地解决工程实际问题打下坚实基础。