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

樊军伟 杨仕教 邓波 孙冰 彭成

摘  要：Winkler地基梁挠曲方程的分析涉及高阶微分方程的求解，从而使得手工进行该项工作既耗时又枯燥。為提高本科生Winkler地基梁的分析水平，在土木工程本科教学中引入Mathematica软件。将作用在Winkler地基梁上的均布荷载、集中力和集中力偶用不同类型的奇异函数在软件中统一表达，借助软件强大的计算能力，对两个算例进行分析计算及内力曲线绘制。两个算例证明Mathematica的引入有利于学生对自己感兴趣的具体边界条件下的探索性问题运用Winkler理论进行分析。

关键词：基础工程；Winkler地基梁；挠曲微分方程；Mathematica软件；奇异函数

中图分类号：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  引  言

基础工程是普通高等学校土木工程及相关专业本科阶段的核心课程。其中连续基础设计时通常假定其下方作为支撑材料的地基处于弹性工作状态，那么弹性地基与其上放置的连续基础即构成了弹性地基梁系统。

所谓弹性地基梁，是指放置在具有一定弹性地基上，各点与地基表面紧密相贴的梁，如铁路枕木、钢筋混凝条形基础梁等[1]。在分析弹性地基梁时，最常用的三种线弹性地基模型分别为：Winkler地基模型、弹性半无限空间地基模型和有限压缩层地基模型[2]。其中 Winkler地基模型所含参数较少、便于工程应用，因此是目前弹性地基上梁分析时最常用的地基模型之一。

当不同类型的多个荷载横向作用在Winkler地基上梁时，挠曲微分方程的手工计算求解极其复杂。即使计算特定位置的某一特征变量值都需要将不同类型的多个荷载在该特定位置的结果进行叠加，而同一地基梁对不同作用位置的荷载通常又属于不同的计算模式（无限长梁、半无限梁、有限长梁等），这也更加剧了Winkler地基梁分析计算的复杂性。鉴于手工分析计算Winkler地基梁问题既耗时又枯燥，有学者将常用的数学软件MATLAB[3]和Maple[4]引入Winkler地基梁的分析中，并取得了不错的结果。然而MATLAB和Maple要求用户具有一定的编程水平才能熟练分析具体边界条件下的Winkler地基梁，这对没有或极少编程经验的本科生来说极不友好。

然而作为“3M”之一的Mathematica语法规则简单，其语法更接近数学运算的思维和表达方式，具有“所见及所得”的特点。同时Mathematica通过较少的语句和简单的命令就能完成复杂微分方程（组）的求解及特征变量曲线绘制。Winkler地基梁的挠曲方程本质上是一个四阶常系数线性非齐次微分方程，因此本文引入Mathematica软件将Winkler地基梁的挠曲微分方程作为纯数学问题进行微分方程（组）求解并绘制出相应的特征值曲线具有可行性。

1  Winkler地基梁模型及挠曲方程

1.1  Winkler地基模型

1867年前后，捷克工程师Winkler提出了如下假设[5]：地基上任意一点所受的压力p与该点的地基沉降s成正比，即：

p=ks                                      （1）

式中，k称为基床系数，单位为kN/m3或MN/m3。

Winkler地基模型把连续的地基分割成了一系列侧面无摩擦的相互独立的土柱，每个土柱的沉降仅与作用在该土柱上的竖向压力有关并与之成正比。

1.2  Winkler地基梁挠曲微分方程

如图1所示Winkler地基梁[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连续两次求导可得：

（5）

式（5）是地基梁的挠曲微分方程，对任意地基模型均适用。

当采用Winkler地基模型时，根据接触条件，梁下地基表面任一点地基沉降等于相应位置梁的挠度：即s=w，将s替换为w并代入式（1）可得：

p=kw                                     （6）

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

（7）

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

（8）

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

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

（9）

式（9）是一个四阶常系数线性非齐次微分方程，它的通解由对应齐次方程的通解和非齐次方程的一个特解经叠加而成。

对于梁的无荷载段（q=0），式（9）变为：

（10）

式（10）即构成了式（9）的线性齐次微分方程，其通解表示为：

（11）

式中，C1，C2，C3和C4为齐次方程（10）通解的4个待定的积分常数，可根据边界条件求出，然后将求出的齐次方程（10）的通解与非齐次方程（9）的一个特解（由梁上荷载类型和作用位置确定）叠加即可获得式（9）的通解。

2  奇异函数

2.1  奇异函数的定义及形式

Mathematica求解Winkler地基梁挠曲微分方程的关键是如何将间断地作用在梁上的不同类型的多个荷载用奇异函数表达成关于x的连续函数表达式。

奇异函数是指本身具有不连续点（跳跃点）或其导数或积分具有不连续点的一族函数，又称脉冲函数或麦考利函数（Macaulay function）[8-9]。其具体表达式为：

fn（x）=                             （12）

式中， 称为麦考利括号（Macaulays Bracket），a为一任意实数。

当麦考利函数的指数n取不同值时，奇异函数有不同的形式，具体表现如下：

当n≥0时

（13）

当n＜0时

（14）

奇异函数具有特殊的定义与形式，对其进行微分和积分运算可以避免积分常数的求解，从而大大减少计算的工作量。

奇异函数的微分和积分的表达形式分别如下：

（15）

（16）

2.2  梁上荷载的统一表达式

奇异函数最重要的优点就是对函数不连续性的处理，Winkler地基梁上不同位置作用有不同类型的多个荷载时可以用不同形式的奇异函数分别表达，然后叠加成一个连續函数的形式统一表达出来，从而避免对微分方程分段积分的弊端，便于计算机编程。

作用在Winkler地基梁上的荷载通常有三种常见类型：均布荷载、集中力和集中力偶，其对应于不同形式的奇异函数[10]。

当n=0时，奇异函数可以将作用在梁上[a，b]区段的均布荷载q转化为荷载函数：

（17）

当n=-1时，奇异函数可以将作用在梁上c点的集中力F转化为荷载函数：

（18）

当n=-2时，奇异函数可以将作用在梁上d点的集中力偶M转化为荷载函数：

（19）

当地基梁上不同位置作用不同类型的若干个荷载时，只需将上述各种类型的多个荷载函数累加在一起即可获得梁上总的荷载函数 。

3  两个算例

3.1  算例1

以文献[11]例3-1（P75）为例，如图2所示的条形基础抗弯刚度EI=4.3×103 MPa·m4，长l=17 m，底面宽b=2.5 m，基床系数k=3.8 MN/m3。试计算基础中心C处的挠度、弯矩和基底反力。

作用在地基梁上有F1=1 200 kN，F2=2 000 kN，F3=2 000 kN及F4=1 200 kN四个集中力和M1=50 kNm及M4=-50 kNm两个集中力偶共6个集中荷载。地基梁上总的荷载函数如下：

（20）

当n=-1时，奇异函数称为狄拉克函数（Dirac delta function），表示为δ（x-a），在Mathematica中由内置的DiracDelta函数实现；当n=-2时，奇异函数称为二阶脉冲函数（double-pulse function），表示为δ（x-a），在Mathematica中由内置DiracDelta函数的导数实现。故该Winkler地基梁的微分挠曲方程表示为：

（21）

由材料力学理论及Winkler地基梁理论，Winkler地基梁的微分挠曲方程也可以用下列微分方程组的形式等效表示为：

（22）

由于地基梁两端自由，地基梁的边界条件表示为：

（23）

由边界条件式（23），在系统命令行里输入挠曲微分方程（21）或微分方程组（22）均可对微分方程求解，本文采用微分方程组的形式进行求解，输入命令如下：

按下“Shift+Enter”键，系统自动计算。然后自动输出转角θ[x]、挠度w[x]、弯矩M[x]、剪力V[x]的解析解，限于篇幅，本算例四个解析解不在文中展示，感兴趣的读者可以参考上述命令自行运算以获得各特征变量解的表达式。

为绘制地基梁的特征变量曲线，例如在系统中输入如下命令：

按下“Shift+Enter”键，系统自动输出基础梁弯矩沿梁长x的曲线，如图3所示。类似的命令在系统中输入并执行后，即可获得基础梁的转角、挠度、剪力及地基反力沿梁长x的曲线图，如图4所示剪力曲线图。

为获得地基梁上具体计算点的特征变量值，以本算例获取C点的挠度、弯矩和基底反力为例，在系统中输入以下命令：

，按下“Shift+Enter”键，系统输出“0.037 886 1” m；

，按下“Shift+Enter”键，系统输出“-1 189.22” kN·m；

，按下“Shift+Enter”键，系统输出“143.967” kPa。

算例1有限长梁采用叠加法手工计算得到C点的挠度、弯矩和基底反力分别为0.038 m、-1 126 kN·m和144.4 kPa。经比较，两种途径计算所得C点三个特征变量值除弯矩值相差较大外，挠度和基底反力几乎相等。

3.2  算例2

以文献[1]例题3-1（P58）为例，如图5所示地基梁，长度l=4 m，宽度b=0.2 m。EI=1 333 kN·m2。地基的弹性压缩系数k=40 000 kN/m3，梁的两端自由。求截面1和截面2的弯矩。

如图5所示，地基梁上作用有一个P=40 kN的集中力和一个q0=20 kN/m的均布荷载。则地基梁上总的荷载函数如下：

（24）

当n=0时，奇异函数称为阶跃函数（Heaviside step function），表示为H（x-a），在Mathematica中由内置的HeavisideTheta函数实现，故该弹性地基梁的微分挠曲方程表示为：

（25）

由于地基梁两端自由，本算例地基梁边界条件表示为：

（26）

由边界条件式（26），在系统命令行里输入微分方程式（25）进行求解，命令如下：

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

然后输入相应的命令，系统自动绘制该地基梁特征变量曲线，如图6所示地基梁弯矩图和图7所示地基梁剪力曲线图。

梁截面1和截面2的弯矩值通过系统计算与文献[1]计算结果对比如表1所示。

文献[1]中梁截面1和截面2的弯矩值采用初参数法计算所得。由表1可知，地基梁上截面1和截面2弯矩手工计算值与系统计算值相差不大。

3.3  算例讨论

以上两个算例均为梁端自由的边界条件。当梁端边界为其他形式的约束时，只需参照上述两个算例的命令进行修改，并给出具体的边界条件即可进行其他形式的Winkler地基梁结构的分析计算。

另外，运用Mathematica系统对土木工程中常见的地基梁进行分析计算时，无须事先判断Winkler地基梁的类型（无限长梁、半无限梁、有限长梁、短梁），只需在系统中输入微分方程及具體的边界条件就能从纯数学的角度得到各特征变量的解析解，并绘制相应特征变量曲线。

事实上，并非所有的高阶微分方程均能求出精确解，当具体边界条件下Winkler地基梁的微分方程无法求得其精确解时，可以用NDSolve命令替代DSolve求微分方程的数值解。

4  结  论

本文以基础工程课程中的重点和难点-Winkler地基梁的分析计算为例，将Mathematica系统引入课堂教学环节通过两个算例进行Winkler地基梁挠曲微分方程（组）的求解。具有以下优点：

（1）系统强大的微分方程（组）求解及曲线绘制功能能够吸引学生课堂学习的注意力，可视化的曲线输出让学生直观地、沉浸式地感受地基梁各特征变量沿梁长的变化规律。

（2）系统的引入让学生通过简单的编程事先获得课后习题的计算结果，有利于提高学生对课后习题进行手工分析计算的积极性，并可以采用系统计算解来验证手工计算结果的正确性。

（3）在本文已有命令的基础上稍加修改就能进行具体边界下其他类型的Winkler地基梁结构的分析计算，为后续课程中桩基础和格构梁等探索性课题的分析计算打下良好的数理基础。

参考文献：

[1] 崔振东.地下结构设计 [M].北京：中国建筑工业出版社，2017.

[2] 周景星，李广信，张建红等.基础工程：第3版 [M].北京：清華大学出版社，2015.

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

[4] 丁洲祥，李涛，白冰，等.MAPLE在土力学与基础工程研究型教学中的应用 [J].力学与实践，2013，35（6）：87-89.

[5] 王协群，章宝华.基础工程 [M].北京：北京大学出版社，2006.

[6] 陈晓平.基础工程设计与分析 [M].北京：中国建筑工业出版社，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].现代信息科技，2022，6（11）：117-121+125.

[10] 王建午，楼京俊，李欣一，等.基于梁变形微分方程与奇异函数的轴系校中计算研究 [J].舰船科学技术，2019，41（11）：71-75.

[11] 华南理工大学，浙江大学，湖南大学.基础工程：第4版 [M].北京：中国建筑工业出版社，2019.

作者简介：樊军伟（1983—），男，汉族，河北邯郸人，讲师，博士研究生，研究方向：岩土工程方面的教学与科研；通讯作者：杨仕教（1964—），男，汉族，湖南浏阳人，教授，博士生导师，博士研究生，研究方向：矿山岩土工程方面的教学与科研。

收稿日期：2022-08-24

基金项目：湖南省普通高等学校教学改革研究项目（HNJG-2021-0622，HNJG-2022-0754）；湖南省教育厅科研项目（20C1608，21C0269）；南华大学校级教学改革项目（2019YB-XJG14，2021YB-XJG34）