MATLAB三维可视化求解在土力学教改中的应用★

吴雪婷 陈 斌

(1.中国地质大学(武汉)工程学院，湖北 武汉 430074； 2.杭州师范大学理学院，浙江 杭州 311121)

1 土力学教学现状及教学改革

土力学是一门建立在数学、力学基础之上偏于计算同时又和工程实践紧密结合的学科，传统的土力学教材中，从模型的建立，到公式的推导，再到实际工程案例的解决已形成了完整的逻辑体系。为了克服公式繁琐、计算量大的问题，传统解法主要采用将计算参数表格化，通过查表的方式来简化计算。尽管如此，还是给学生留下土力学计算机械重复、枯燥无味的印象，很难调动学生的学习积极性，教学效果欠佳。

计算机辅助教学技术的普及和数学计算软件的快速发展，为土力学教学中改善以上问题提供了新思路[1]。MATLAB是一款功能强大的数学计算软件，不仅可以实现矩阵的加、减、乘、除、可逆，以及转秩的各项功能运算；而且可以运用MATLAB建立数学模型，实现模型的优化设计并以此解决众多实际问题。同时，其可视化功能可以实现二维、三维图形的绘制，以便于对需要解决的问题进行直观分析[2]。因此，可以将MATLAB应用到土力学计算、实验数据处理、建模等土力学的教学中[3-5]，来改革传统的土力学教学模式[6-8]。

基于此，本文以土力学中竖向集中荷载作用下地基中的附加应力为例，在对其传统一般解法进行梳理的基础上，对于繁琐的计算、可视化呈现等方面，借助MATLAB软件，将学生从繁琐的计算中解脱出来，让学生能够关注土力学问题本身，帮助学生理解计算结果的物理和工程意义。

2 地基附加应力的传统查表法求解

附加应力计算中的半无限空间受法向集中力作用的问题是由法国数学家、物理学家布辛奈斯克(Boussinesq)在1885年依据弹性理论求解得出。目前在计算地基中的附加应力时都是借助于此方法，把地基看成是均质的弹性半空间，应用弹性力学理论求解[9]。

如图1所示，在地基表面作用有竖向集中荷载p时，在地基内任意一点M(x，y，z)会产生相应的应力分量及位移分量。竖向(z轴方向)地基附加应力值σz为：

(1)

其中，p为作用在坐标原点O的竖向集中荷载；x，y，z均为地基内任意点M的坐标，z为M点的深度，且r2=x2+y2,R2=x2+y2+z2。

(2)

其中，a为集中荷载作用下的地基竖向附加应力系数，是r/z的函数。根据不同的r/z值可算出附加应力系数a值，并列成表格。计算时再根据具体的r/z值，查表得到a值，最后再计算σz。

例如，在p=200 kN，求地基中z=2 m的水平面上不同半径r处的地基竖向附加应力时，传统教材和教学中是通过查表得到相应的附加应力系数a，再利用式(2)计算得到附加应力，计算结果如图2所示。

同样，在p=200 kN，求地基中r=0 m的竖直线上不同深度z处的地基竖向附加应力时，传统教材和教学中也是通过查表求得，计算结果如图3所示。

从上述求解的过程和绘制的图表可以看出，附加应力系数表格起着非常重要的作用，但给学生的印象就是机械地查表和重复地计算。同时，学生也无法从简单的数据点和二维曲线中想象出地基中附加应力的三维空间分布。

3 基于MATLAB的地基附加应力三维可视化求解

针对上述问题，本文利用MATLAB软件，对竖向集中荷载作用下地基中的附加应力进行了三维可视化求解，很好地弥补了传统教材、传统教学的解析解、查表法的不足。

3.1 地基附加应力随半径r的变化规律

和上述查表法例题一致，同样假设p=200 kN，求地基中z=2 m的水平面上，不同半径r处的地基竖向附加应力值。MATLAB软件可直接利用式(1)计算任意深度z处的附加应力随r的变化数据，并通过“Plot”命令将其以图形的形式呈现出来，计算结果如图4所示，分别给出了z=1 m,1.5 m,2 m,4 m深度处的附加应力随半径r的变化曲线。和图2对比可见，图4中的数据点更多，曲线显示的变化规律更明显，但是计算却更快捷。

此外，利用MATLAB三维可视化功能还可帮助学生建立地基附加应力的三维空间分布概念。图5给出了不同深度z=1.0 m,1.5 m,2.0 m,4.0 m时，地基附加应力在某一固定深度z处的三维空间分布图。图6展示了深度z=1.5 m和2 m时，地基附加应力的三维空间综合图(可通过MATLAB的surf函数实现)。由图4～图6可知，在地面下同一深度处，沿竖向集中力作用线上的附加应力最大，向两侧逐渐减小；在集中力作用线上(r=0 m)，深度越大，σz越小。

3.2 地基附加应力随深度z的变化规律

仍然假设p=200 kN，求地基中r=0 m的竖直线上，不同深度z处的地基竖向附加应力。利用查表法只是计算r=0 m处的附加应力随深度的变化已经有非常大的计算量了，但是，利用MATLAB编程可以快捷地计算并绘制r=0 m和r≠0 m(r=0.25 m,0.5 m,1 m)的任意竖直线上的地基附加应力随深度z的变化曲线，MATLAB计算结果如图7所示。同理，和图3对比可见MATLAB程序计算的优势。由图7可知，在集中力作用线上(r=0 m时)，当z=0 m时，σz→∞，随着深度增加，σz逐渐减小；在r≠0 m时的任一竖直线上，随着深度增加，σz先增大后减小。

3.3 地基附加应力等值面

如上所述，分别给出了某一深度不同半径，以及相同半径不同深度时附加应力的变化规律。那么，相同应力点在空间的分布情况又是怎样的呢。如果采用传统的解析计算和查表法反算这些相等的附加应力值对应的点坐标，计算量巨大繁琐，不仅浪费时间，同时也打消了学生学习的积极性。但是对于采用MATLAB软件计算而言，就是求式(1)中附加应力给一定值时，空间坐标(x,y,z)构成的等值面问题。

分别取地基附加应力σz=22 kPa,12 kPa,7 kPa,5 kPa,2.5 kPa,1.5 kPa，采用MATLAB的“isosurface”函数就可以实现三维可视化求解，并得到在空间分布的附加应力等值面，MATLAB计算结果如图8所示。利用编制好的附加应力等值面计算程序，只需要改变程序中σz值的大小，即可以求解任意地基附加应力σz值对应的等值面。附加应力等值曲面的空间形态如泡状，所以可称为“应力泡”。

根据上述MATLAB的计算结果，即可以方便、形象、直观地得出集中力作用下地基土中附加应力的分布规律。

4 结语

基于MATLAB的三维可视化求解方法在土力学地基附加应力的繁琐、机械、重复计算中得到了很好地应用，一方面提高了学生的学习兴趣，加深了学生对土力学知识的理解，帮助学生去更好地理解地基附加应力的空间分布和变化规律；另一方面教会了学生利用现代化手段解决工程实际问题的方法，在提高学生物理建模和数学计算能力的同时，也为其将来利用计算机软件解决工程问题打下了坚实的基础。此外，MATLAB软件辅助教学法还可以推广应用到土力学课程的其他教学内容中。