Mathematica在Winkler地基梁分析中的應用

樊軍偉 楊仕教 鄧波 孫冰 彭成

摘? 要：Winkler地基梁撓曲方程的分析涉及高階微分方程的求解，從而使得手工進行該項工作既耗時又枯燥。為提高本科生Winkler地基梁的分析水平，在土木工程本科教學中引入Mathematica軟件。將作用在Winkler地基梁上的均布荷載、集中力和集中力偶用不同類型的奇異函數(shù)在軟件中統(tǒng)一表達，借助軟件強大的計算能力，對兩個算例進行分析計算及內(nèi)力曲線繪制。兩個算例證明Mathematica的引入有利于學生對自己感興趣的具體邊界條件下的探索性問題運用Winkler理論進行分析。

關鍵詞：基礎工程；Winkler地基梁；撓曲微分方程；Mathematica軟件；奇異函數(shù)

中圖分類號：TP39；G434? ? 文獻標識碼：A 文章編號：2096-4706（2023）01-0186-05

Application of Mathematica in the Analysis of Winkler Foundation Beam

FAN Junwei， YANG Shijiao， DENG Bo， SUN Bing， PENG Cheng

（University of South China， Hengyang? 421001， China）

Abstract： The analysis of the deflection equation of Winkler foundation beam involves the solving of a higher-order differential equation， which makes this work by hand time-consuming and boring. In order to improve the undergraduates' analysis level of Winkler foundation beam， Mathematica software is introduced in the undergraduate teaching of civil engineering. Uniform loads， concentrated loads and concentrated moments acting on the beam are uniformly expressed in the software by different types of singular functions. With the help of the powerful calculating capacity of the software， two examples are analyzed and calculated， and the internal force curves are drawn. Two examples show that the introduction of Mathematica is helpful for students to analyze the exploratory problems by using Winkler theory under specific boundary conditions that they are interested in.

Keywords： foundation engineering; Winkler foundation beam;deflection of differential equation;Mathematica software; singularity function

0? 引? 言

基礎工程是普通高等學校土木工程及相關專業(yè)本科階段的核心課程。其中連續(xù)基礎設計時通常假定其下方作為支撐材料的地基處于彈性工作狀態(tài)，那么彈性地基與其上放置的連續(xù)基礎即構成了彈性地基梁系統(tǒng)。

所謂彈性地基梁，是指放置在具有一定彈性地基上，各點與地基表面緊密相貼的梁，如鐵路枕木、鋼筋混凝條形基礎梁等[1]。在分析彈性地基梁時，最常用的三種線彈性地基模型分別為：Winkler地基模型、彈性半無限空間地基模型和有限壓縮層地基模型[2]。其中 Winkler地基模型所含參數(shù)較少、便于工程應用，因此是目前彈性地基上梁分析時最常用的地基模型之一。

當不同類型的多個荷載橫向作用在Winkler地基上梁時，撓曲微分方程的手工計算求解極其復雜。即使計算特定位置的某一特征變量值都需要將不同類型的多個荷載在該特定位置的結果進行疊加，而同一地基梁對不同作用位置的荷載通常又屬于不同的計算模式（無限長梁、半無限梁、有限長梁等），這也更加劇了Winkler地基梁分析計算的復雜性。鑒于手工分析計算Winkler地基梁問題既耗時又枯燥，有學者將常用的數(shù)學軟件MATLAB[3]和Maple[4]引入Winkler地基梁的分析中，并取得了不錯的結果。然而MATLAB和Maple要求用戶具有一定的編程水平才能熟練分析具體邊界條件下的Winkler地基梁，這對沒有或極少編程經(jīng)驗的本科生來說極不友好。

然而作為“3M”之一的Mathematica語法規(guī)則簡單，其語法更接近數(shù)學運算的思維和表達方式，具有“所見及所得”的特點。同時Mathematica通過較少的語句和簡單的命令就能完成復雜微分方程（組）的求解及特征變量曲線繪制。Winkler地基梁的撓曲方程本質(zhì)上是一個四階常系數(shù)線性非齊次微分方程，因此本文引入Mathematica軟件將Winkler地基梁的撓曲微分方程作為純數(shù)學問題進行微分方程（組）求解并繪制出相應的特征值曲線具有可行性。

1? Winkler地基梁模型及撓曲方程

1.1? Winkler地基模型

1867年前后，捷克工程師Winkler提出了如下假設[5]：地基上任意一點所受的壓力p與該點的地基沉降s成正比，即：

p=ks? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （1）

式中，k稱為基床系數(shù)，單位為kN/m3或MN/m3。

Winkler地基模型把連續(xù)的地基分割成了一系列側面無摩擦的相互獨立的土柱，每個土柱的沉降僅與作用在該土柱上的豎向壓力有關并與之成正比。

1.2? Winkler地基梁撓曲微分方程

如圖1所示W(wǎng)inkler地基梁[6-7]，在外荷載作用下梁的撓度曲線和地基反力曲線如圖1（a）所示，沿x軸的方向在梁上取長度為dx的一個微段進行分析，見圖1（b）。

（a）梁上荷載及梁的撓曲

（b）梁的微段

材料力學中，由梁的純彎曲得到撓曲微分方程如下：

（2）

式中，w為梁的撓度，M為梁的彎矩，E為梁體彈性模量，I為梁的截面慣性矩。

如圖1（b）所示由微段的靜力平衡條件 、 可得：

（3）

（4）

式中，V為剪力，q為梁上作用的分布荷載，p為地基反力，b為梁底寬度。

將式（2）對坐標x連續(xù)兩次求導可得：

（5）

式（5）是地基梁的撓曲微分方程，對任意地基模型均適用。

當采用Winkler地基模型時，根據(jù)接觸條件，梁下地基表面任一點地基沉降等于相應位置梁的撓度：即s=w，將s替換為w并代入式（1）可得：

p=kw? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（6）

再將式（6）代入式（5）整理后可得：

（7）

式（7）即為Winkler地基梁的撓曲微分方程。為求解方便，令：

（8）

式中，λ為梁的柔度特征值，量綱為[L-1]。

將式（8）代入式（7）得：

（9）

式（9）是一個四階常系數(shù)線性非齊次微分方程，它的通解由對應齊次方程的通解和非齊次方程的一個特解經(jīng)疊加而成。

對于梁的無荷載段（q=0），式（9）變?yōu)椋?/p>

（10）

式（10）即構成了式（9）的線性齊次微分方程，其通解表示為：

（11）

式中，C1，C2，C3和C4為齊次方程（10）通解的4個待定的積分常數(shù)，可根據(jù)邊界條件求出，然后將求出的齊次方程（10）的通解與非齊次方程（9）的一個特解（由梁上荷載類型和作用位置確定）疊加即可獲得式（9）的通解。

2? 奇異函數(shù)

2.1? 奇異函數(shù)的定義及形式

Mathematica求解Winkler地基梁撓曲微分方程的關鍵是如何將間斷地作用在梁上的不同類型的多個荷載用奇異函數(shù)表達成關于x的連續(xù)函數(shù)表達式。

奇異函數(shù)是指本身具有不連續(xù)點（跳躍點）或其導數(shù)或積分具有不連續(xù)點的一族函數(shù)，又稱脈沖函數(shù)或麥考利函數(shù)（Macaulay function）[8-9]。其具體表達式為：

fn（x）=? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（12）

式中， 稱為麥考利括號（Macaulays Bracket），a為一任意實數(shù)。

當麥考利函數(shù)的指數(shù)n取不同值時，奇異函數(shù)有不同的形式，具體表現(xiàn)如下：

當n≥0時

（13）

當n＜0時

（14）

奇異函數(shù)具有特殊的定義與形式，對其進行微分和積分運算可以避免積分常數(shù)的求解，從而大大減少計算的工作量。

奇異函數(shù)的微分和積分的表達形式分別如下：

（15）

（16）

2.2? 梁上荷載的統(tǒng)一表達式

奇異函數(shù)最重要的優(yōu)點就是對函數(shù)不連續(xù)性的處理，Winkler地基梁上不同位置作用有不同類型的多個荷載時可以用不同形式的奇異函數(shù)分別表達，然后疊加成一個連續(xù)函數(shù)的形式統(tǒng)一表達出來，從而避免對微分方程分段積分的弊端，便于計算機編程。

作用在Winkler地基梁上的荷載通常有三種常見類型：均布荷載、集中力和集中力偶，其對應于不同形式的奇異函數(shù)[10]。

當n=0時，奇異函數(shù)可以將作用在梁上[a，b]區(qū)段的均布荷載q轉化為荷載函數(shù)：

（17）

當n=-1時，奇異函數(shù)可以將作用在梁上c點的集中力F轉化為荷載函數(shù)：

（18）

當n=-2時，奇異函數(shù)可以將作用在梁上d點的集中力偶M轉化為荷載函數(shù)：

（19）

當?shù)鼗荷喜煌恢米饔貌煌愋偷娜舾蓚€荷載時，只需將上述各種類型的多個荷載函數(shù)累加在一起即可獲得梁上總的荷載函數(shù) 。

3? 兩個算例

3.1? 算例1

以文獻[11]例3-1（P75）為例，如圖2所示的條形基礎抗彎剛度EI=4.3×103 MPa·m4，長l=17 m，底面寬b=2.5 m，基床系數(shù)k=3.8 MN/m3。試計算基礎中心C處的撓度、彎矩和基底反力。

作用在地基梁上有F1=1 200 kN，F(xiàn)2=2 000 kN，F(xiàn)3=2 000 kN及F4=1 200 kN四個集中力和M1=50 kNm及M4=-50 kNm兩個集中力偶共6個集中荷載。地基梁上總的荷載函數(shù)如下：

（20）

當n=-1時，奇異函數(shù)稱為狄拉克函數(shù)（Dirac delta function），表示為δ（x-a），在Mathematica中由內(nèi)置的DiracDelta函數(shù)實現(xiàn)；當n=-2時，奇異函數(shù)稱為二階脈沖函數(shù)（double-pulse function），表示為δ（x-a），在Mathematica中由內(nèi)置DiracDelta函數(shù)的導數(shù)實現(xiàn)。故該Winkler地基梁的微分撓曲方程表示為：

（21）

由材料力學理論及Winkler地基梁理論，Winkler地基梁的微分撓曲方程也可以用下列微分方程組的形式等效表示為：

（22）

由于地基梁兩端自由，地基梁的邊界條件表示為：

（23）

由邊界條件式（23），在系統(tǒng)命令行里輸入撓曲微分方程（21）或微分方程組（22）均可對微分方程求解，本文采用微分方程組的形式進行求解，輸入命令如下：

按下“Shift+Enter”鍵，系統(tǒng)自動計算。然后自動輸出轉角θ[x]、撓度w[x]、彎矩M[x]、剪力V[x]的解析解，限于篇幅，本算例四個解析解不在文中展示，感興趣的讀者可以參考上述命令自行運算以獲得各特征變量解的表達式。

為繪制地基梁的特征變量曲線，例如在系統(tǒng)中輸入如下命令：

按下“Shift+Enter”鍵，系統(tǒng)自動輸出基礎梁彎矩沿梁長x的曲線，如圖3所示。類似的命令在系統(tǒng)中輸入并執(zhí)行后，即可獲得基礎梁的轉角、撓度、剪力及地基反力沿梁長x的曲線圖，如圖4所示剪力曲線圖。

為獲得地基梁上具體計算點的特征變量值，以本算例獲取C點的撓度、彎矩和基底反力為例，在系統(tǒng)中輸入以下命令：

，按下“Shift+Enter”鍵，系統(tǒng)輸出“0.037 886 1” m；

，按下“Shift+Enter”鍵，系統(tǒng)輸出“-1 189.22” kN·m；

，按下“Shift+Enter”鍵，系統(tǒng)輸出“143.967” kPa。

算例1有限長梁采用疊加法手工計算得到C點的撓度、彎矩和基底反力分別為0.038 m、-1 126 kN·m和144.4 kPa。經(jīng)比較，兩種途徑計算所得C點三個特征變量值除彎矩值相差較大外，撓度和基底反力幾乎相等。

3.2? 算例2

以文獻[1]例題3-1（P58）為例，如圖5所示地基梁，長度l=4 m，寬度b=0.2 m。EI=1 333 kN·m2。地基的彈性壓縮系數(shù)k=40 000 kN/m3，梁的兩端自由。求截面1和截面2的彎矩。

如圖5所示，地基梁上作用有一個P=40 kN的集中力和一個q0=20 kN/m的均布荷載。則地基梁上總的荷載函數(shù)如下：

（24）

當n=0時，奇異函數(shù)稱為階躍函數(shù)（Heaviside step function），表示為H（x-a），在Mathematica中由內(nèi)置的HeavisideTheta函數(shù)實現(xiàn)，故該彈性地基梁的微分撓曲方程表示為：

（25）

由于地基梁兩端自由，本算例地基梁邊界條件表示為：

（26）

由邊界條件式（26），在系統(tǒng)命令行里輸入微分方程式（25）進行求解，命令如下：

按下“Shift+Enter”鍵，即可獲得該微分方程中w[x]的解析解。

然后輸入相應的命令，系統(tǒng)自動繪制該地基梁特征變量曲線，如圖6所示地基梁彎矩圖和圖7所示地基梁剪力曲線圖。

梁截面1和截面2的彎矩值通過系統(tǒng)計算與文獻[1]計算結果對比如表1所示。

文獻[1]中梁截面1和截面2的彎矩值采用初參數(shù)法計算所得。由表1可知，地基梁上截面1和截面2彎矩手工計算值與系統(tǒng)計算值相差不大。

3.3? 算例討論

以上兩個算例均為梁端自由的邊界條件。當梁端邊界為其他形式的約束時，只需參照上述兩個算例的命令進行修改，并給出具體的邊界條件即可進行其他形式的Winkler地基梁結構的分析計算。

另外，運用Mathematica系統(tǒng)對土木工程中常見的地基梁進行分析計算時，無須事先判斷Winkler地基梁的類型（無限長梁、半無限梁、有限長梁、短梁），只需在系統(tǒng)中輸入微分方程及具體的邊界條件就能從純數(shù)學的角度得到各特征變量的解析解，并繪制相應特征變量曲線。

事實上，并非所有的高階微分方程均能求出精確解，當具體邊界條件下Winkler地基梁的微分方程無法求得其精確解時，可以用NDSolve命令替代DSolve求微分方程的數(shù)值解。

4? 結? 論

本文以基礎工程課程中的重點和難點-Winkler地基梁的分析計算為例，將Mathematica系統(tǒng)引入課堂教學環(huán)節(jié)通過兩個算例進行Winkler地基梁撓曲微分方程（組）的求解。具有以下優(yōu)點：

（1）系統(tǒng)強大的微分方程（組）求解及曲線繪制功能能夠吸引學生課堂學習的注意力，可視化的曲線輸出讓學生直觀地、沉浸式地感受地基梁各特征變量沿梁長的變化規(guī)律。

（2）系統(tǒng)的引入讓學生通過簡單的編程事先獲得課后習題的計算結果，有利于提高學生對課后習題進行手工分析計算的積極性，并可以采用系統(tǒng)計算解來驗證手工計算結果的正確性。

（3）在本文已有命令的基礎上稍加修改就能進行具體邊界下其他類型的Winkler地基梁結構的分析計算，為后續(xù)課程中樁基礎和格構梁等探索性課題的分析計算打下良好的數(shù)理基礎。

參考文獻：

[1] 崔振東.地下結構設計 [M].北京：中國建筑工業(yè)出版社，2017.

[2] 周景星，李廣信，張建紅等.基礎工程：第3版 [M].北京：清華大學出版社，2015.

[3] 羅汀，姚仰平，胡賀祥.高等基礎工程 [M].北京：人民交通出版社，2013.

[4] 丁洲祥，李濤，白冰，等.MAPLE在土力學與基礎工程研究型教學中的應用 [J].力學與實踐，2013，35（6）：87-89.

[5] 王協(xié)群，章寶華.基礎工程 [M].北京：北京大學出版社，2006.

[6] 陳曉平.基礎工程設計與分析 [M].北京：中國建筑工業(yè)出版社，2005.

[7] 富海鷹.基礎工程：第3版 [M].北京：中國鐵道出版社，2019.

[8] BOEDO S. Singularity Functions Revisited：Clarifications and Extensions for Construction of Shear–Moment Diagrams in Beams [J].International Journal of Mechanical Engineering Education，2020，48（4）：351-370.

[9] 曹建華，郭東旭.Mathematica在材料力學彎曲問題中的應用 [J].現(xiàn)代信息科技，2022，6（11）：117-121+125.

[10] 王建午，樓京俊，李欣一，等.基于梁變形微分方程與奇異函數(shù)的軸系校中計算研究 [J].艦船科學技術，2019，41（11）：71-75.

[11] 華南理工大學，浙江大學，湖南大學.基礎工程：第4版 [M].北京：中國建筑工業(yè)出版社，2019.

作者簡介：樊軍偉（1983—），男，漢族，河北邯鄲人，講師，博士研究生，研究方向：巖土工程方面的教學與科研；通訊作者：楊仕教（1964—），男，漢族，湖南瀏陽人，教授，博士生導師，博士研究生，研究方向：礦山巖土工程方面的教學與科研。

收稿日期：2022-08-24

基金項目：湖南省普通高等學校教學改革研究項目（HNJG-2021-0622，HNJG-2022-0754）；湖南省教育廳科研項目（20C1608，21C0269）；南華大學校級教學改革項目（2019YB-XJG14，2021YB-XJG34）