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

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

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

中圖分類號(hào)：TP39；G434? ? 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：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? 引? 言

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

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

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

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

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

1.1? Winkler地基模型

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

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

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

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

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

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

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

（b）梁的微段

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

（2）

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

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

（3）

（4）

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

將式（2）對(duì)坐標(biāo)x連續(xù)兩次求導(dǎo)可得：

（5）

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

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

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

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

（7）

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

（8）

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

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

（9）

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

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

（10）

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

（11）

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

2? 奇異函數(shù)

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

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

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

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

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

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

當(dāng)n≥0時(shí)

（13）

當(dāng)n＜0時(shí)

（14）

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

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

（15）

（16）

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

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

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

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

（17）

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

（18）

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

（19）

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

3? 兩個(gè)算例

3.1? 算例1

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

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

（20）

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

（21）

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

（22）

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

（23）

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

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

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

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

為獲得地基梁上具體計(jì)算點(diǎn)的特征變量值，以本算例獲取C點(diǎn)的撓度、彎矩和基底反力為例，在系統(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有限長梁采用疊加法手工計(jì)算得到C點(diǎn)的撓度、彎矩和基底反力分別為0.038 m、-1 126 kN·m和144.4 kPa。經(jīng)比較，兩種途徑計(jì)算所得C點(diǎn)三個(gè)特征變量值除彎矩值相差較大外，撓度和基底反力幾乎相等。

3.2? 算例2

以文獻(xiàn)[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所示，地基梁上作用有一個(gè)P=40 kN的集中力和一個(gè)q0=20 kN/m的均布荷載。則地基梁上總的荷載函數(shù)如下：

（24）

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

（25）

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

（26）

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

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

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

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

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

3.3? 算例討論

以上兩個(gè)算例均為梁端自由的邊界條件。當(dāng)梁端邊界為其他形式的約束時(shí)，只需參照上述兩個(gè)算例的命令進(jìn)行修改，并給出具體的邊界條件即可進(jìn)行其他形式的Winkler地基梁結(jié)構(gòu)的分析計(jì)算。

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

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

4? 結(jié)? 論

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

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

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

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

參考文獻(xiàn)：

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

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

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

[4] 丁洲祥，李濤，白冰，等.MAPLE在土力學(xué)與基礎(chǔ)工程研究型教學(xué)中的應(yīng)用 [J].力學(xué)與實(shí)踐，2013，35（6）：87-89.

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

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

[7] 富海鷹.基礎(chǔ)工程：第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在材料力學(xué)彎曲問題中的應(yīng)用 [J].現(xiàn)代信息科技，2022，6（11）：117-121+125.

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

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

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

收稿日期：2022-08-24

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