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軟件輔助教學法還可以推廣應用到土力學課程的其他教學內容中。