Mathematica在材料力學彎曲問題中的應用

曹建華,郭東旭

（黃山學院，安徽 黃山 245041）

0 引 言

對于機械類、材料類等專業(yè)，材料力學是一門必修的專業(yè)基礎(chǔ)課，而彎曲變形在材料力學中所占篇幅非常大，其所研究的問題是學生學習的難點。彎曲變形所研究的問題包括剪力，彎矩、轉(zhuǎn)角′（）和撓度（）的求解，一般采用截面法或積分法方法。積分法是利用剪力、彎矩、轉(zhuǎn)角′（）和撓度（）與分布載荷（）的積分關(guān)系，在作用載荷連續(xù)的一段梁內(nèi)會產(chǎn)生4個常數(shù)，當梁上出現(xiàn)不連續(xù)載荷時，如梁上作用著多個集中外力、多個集中外力偶，或者作用在多個不連續(xù)局部的分布載荷，需要分段寫出剪力，彎矩、轉(zhuǎn)角′（）和撓度（）等方程。假設(shè)載荷不連續(xù)的梁分成段求解，需要4個常數(shù)需要確定，是一個繁重的任務(wù)，且容易出錯。文獻[1]應用奇異函數(shù)法求解了一根直梁，也沒有相應的程序代碼。本文在詳述奇異函數(shù)性質(zhì)的基礎(chǔ)上，應用奇異函數(shù)法求解了一根或多根桿件組成的直梁彎曲變形問題，與文獻對比進行驗證，最后應用Mathematica數(shù)學軟件實現(xiàn)奇異函數(shù)法求解彎曲變形問題的過程。

1 奇異函數(shù)性質(zhì)

1.1 δ函數(shù)、單位階躍函數(shù)和單位斜坡函數(shù)

(1)狄拉克函數(shù)(Dirac Delta function)，具有如下性質(zhì)：

(2)單位階躍函數(shù)（unit step function）（-）為：

(3)單位斜坡函數(shù)（unit ramp function）（-a）為：

狄拉克函數(shù)（-）與單位階躍函數(shù)（-）的關(guān)系為：

1.2 奇異函數(shù)的定義

奇異函數(shù)定義如下：

其中，為整數(shù)，且當≥0時

當＜0時

由以上定義可知如下關(guān)系：

由定義可知，奇異函數(shù)有如下性質(zhì)：

（2）當≥0時，由（1）（2）性質(zhì)積分可得性質(zhì)(3)。

式中，為積分常數(shù)。當＜0時，R（0）是沒有意義的，不參與計算。

2 應用奇異函數(shù)求解梁彎曲變形問題

2.1 梁問題的統(tǒng)一表達式

如圖所示，結(jié)構(gòu)由兩根梁件通過9 m處的鉸鏈鏈接在一起，承受有均勻分布載荷，集中外力，集中外力偶、線性變化的分布載荷以及支座反力、固定端約束力。在眾多教材中，梁彎曲問題的求解一般步驟如下：

(1)以梁的軸左端為原點，表示沿軸線的橫坐標。

(2)按約束（支座和固定端）性質(zhì)，畫出約束反力，并通過平衡方程求解約束反力。

(3)對不連續(xù)載荷進行分析，計算載荷不連續(xù)處總數(shù)，則需要分-1段。如圖1所示，載荷不連續(xù)處包括分布載荷的起始終止處、集中外力處、集中力偶處、鉸鏈連接處和支座反力、固定端處，如果不連續(xù)處重疊（如3 m處支座反力和分布載荷終止處重合），算作1處不連續(xù)處，總共有7段不連續(xù)處，需要分6段。

圖1 受多個載荷的梁

(4)分段進行求解出剪力，彎矩、轉(zhuǎn)角′（）和撓度（），通過邊界條件求出積分數(shù)。

梁和梁的坐標系如圖1所示。以梁的軸左端為原點，表示沿軸線的橫坐標，表示梁的軸線撓度，向下為負。圖2顯示彎矩和剪力的符號規(guī)定。本文中，規(guī)定作用在梁上的集中外力和分布外力向下為正，逆時針方向的力偶為正。

圖2 梁的內(nèi)力正方向

設(shè)（）代表梁上某一連續(xù)段（，）的載荷集度，符號規(guī)定如下：向上的載荷為正，順時針的力偶為正，反之為負。

在＜＜梁內(nèi)，剪力和分布載荷函數(shù)的積分關(guān)系表示如下：

其中為常數(shù)。對上式再次積分，可得彎矩表達式為：

由 = 可知，通過積分可得轉(zhuǎn)角′（）和撓度（）表達式：

其中，、、、為積分常數(shù)，（）是梁的抗彎剛度。

由以上可知，如果將梁上的所有不連續(xù)載荷都能包含在載荷函數(shù)（）表達式中，就能夠無需分段，并求出剪力，彎矩、轉(zhuǎn)角′（）和撓度（）。奇異函數(shù)具有將不連續(xù)載荷包含在一個表達式并進行計算的能力，將不連續(xù)的問題變成了連續(xù)問題處理，簡化了很多計算。下面敘述采用奇異函數(shù)法來表示（）。

(1)作用于=處的集中外力，其載荷函數(shù)）表達式為：

(2)對于在區(qū)間[，]內(nèi)線性變化的分布載荷：

(3)若是在區(qū)間[，]的分布載荷函數(shù)F（），載荷函數(shù)（）表達式為：

(4)作用于=處的集中外力偶，其載荷函數(shù)（）表達式為：

(5)=（0＜＜）處鉸鏈連接處或剪力連接，在此處截面轉(zhuǎn)角發(fā)生突變，其載荷函數(shù)（）表達式為：

(6)=（0＜＜）處彎矩連接，在此處截面處位移發(fā)生突變，其載荷函數(shù)（）表達式為：

2.2 梁的邊界條件

(2)剛性支撐：=0，如圖1中3 m處的支撐。

(3)左端簡支邊界：=0，+M=0（M為作用于左端的集中外力偶），如圖3中左端邊界。

圖3 彎矩連接的梁

(4)右端簡支邊界：=0，-M=0（M為作用于右端的集中外力偶）。

(5)左端自由邊界：+V=0，+M=0（M、V為作用于左端的集中外力偶和集中外力），如圖1中左端邊界。

(6)右端自由邊界：-V=0，-M=0（M、V為作用于右端的集中外力偶和集中外力）。

(7)左端引導邊界：+V=0，=0（M、V為作用于左端的集中外力偶和集中外力）。

(8)右端引導邊界：-V=0，=0（M、V為作用于右端的集中外力偶和集中外力），如圖4所示。

圖4 右端是引導邊界的直梁

(9)鉸鏈連接（a hinge or moment release）：=0，如圖1中9 m處連接，只能傳遞剪力，不能傳遞彎矩。

(10)彎矩連接（shear release）：=0，如圖3所示，只能傳遞彎矩，不能傳遞剪力。

2.3 梁問題的求解步驟

用奇異函數(shù)法求解剪力，彎矩、轉(zhuǎn)角′（）和撓度（），步驟如下：

(1)建立坐標系，默認以左端點為原點，平行于軸線為軸。

(2)奇異函數(shù)寫出沿梁長度的載荷集度（）。

(3)代入邊界條件，求出常 數(shù)、、、以及Δ、Δ。

(4)作出剪力，彎矩、轉(zhuǎn)角′（）和撓度（）圖形。

3 求解實例及Mathematica求解

3.1 求解實例的演算過程

如圖5所示，承受兩根桿件組成的直梁，通過鉸鏈連接，承受2個集中外力。本實例來自劉鴻文《材料力學.I》6版中習題4.6(b)。

圖5 梁的幾何和受力

解：第一步：寫出梁的載荷函數(shù)：

第二步：根據(jù)載荷函數(shù)，通過積分寫出剪力，彎矩、轉(zhuǎn)角′（）和撓度（）的方程分別為：

第三步：寫出邊界條件，并代入第二步中相應的四個方程中，求出積分常數(shù)。本題是兩端簡支梁，其邊界條件為：

x=0，U=0，S=0

x=6，U=0，M=0

x=4，M=0

代入第二步中相應的方程，可得：

第四步：將第三步所求出的常數(shù)代入第二步的方程中，可得彎曲問題的方程，同上，略去不表。

3.2 Mathematica實現(xiàn)奇異函數(shù)法

Mathematica具有強大的數(shù)值和符號計算能力、強大的圖形功能等，深受廣大從事科研和工程人員的歡迎。Mathematica有內(nèi)置函數(shù)DiracDelta實現(xiàn)δ函數(shù)，HeavisideTheta函數(shù)實現(xiàn)單位階躍函數(shù)。當n＜-1時，R（x-a）可以用DiracDelta導數(shù)表示，非常方便。下面是基于Mathematica軟件實現(xiàn)的奇異函數(shù)法程序代碼，只需要根據(jù)具體問題，修改粗體部分數(shù)據(jù)，即可用于求解任意直梁彎曲問題。代碼內(nèi)容如下：

Clear[“Global`*”]

SetDirectory[]

(*集中外力對應的載荷函數(shù)q(x) *)

Fc[{P_, a_}][x_] := -P DiracDelta[-a + x];

(*分布載荷對應的載荷函數(shù)q(x) *)

Fdist[{p_, a_, q_, b_}][x_] := -((b p - a q - p x + q x)(HeavisideTheta[-a + x] - HeavisideTheta[-b + x]))/(-a + b);

(*集中外力偶對應的載荷函數(shù)q(x) *)

Cp[{M_, a_}][x_] := -M DiracDelta’[x - a];

(*鉸鏈連接對應的載荷函數(shù)q(x) *)

Sc[{dtheta_, a_}][x_] := -dthetaDiracDelta’’[x - a];

(*彎矩連接對應的載荷函數(shù)q(x) *)

Mc[{dw_, a_}][x_] := -dwDiracDelta’’’[x - a];

(*梁的邊界條件和參數(shù)輸入 *)

bnd={“Fixed”,”Pinned”};(*邊界條件設(shè)置*)

Le=6; (*梁的總長*)

EI=1; (*梁的抗彎剛度，為了方便設(shè)為1*)

PointLoad={50,2,50,5};(*集中外力設(shè)置*)

LineLoad ={};(*分布載荷設(shè)置*)

PointCouple= {} ; (*集中外力偶設(shè)置*)

InSpan ={};(*剛性支撐設(shè)置*)

ShearConnection ={4};(*鉸鏈連接設(shè)置*)

MomentConnection = {};(*彎矩連接設(shè)置*)

(*生成由載荷函數(shù)各項組成的鏈表*)

px = Join[

Map[ Fc[#1][x] &, DeleteCases[Join[Partition[PointLoad,2],InSpan], {_, 0 | 0.0} | {_, Le}] ],

Flatten[Map[ Fdist[#1][x] &, Partition[LineLoad, 4] ]] /.HeavisideTheta[x - Le] :> 0,

Map[ Sc[#1][x] &, ShearConnection ],

Map[ Mc[#1][x] &, MomentConnection ],

Map[ Cp[#1][x] &, DeleteCases[ Partition[PointCouple, 2],{_, 0 | 0.0} | {_, Le} ]]

];

(*通過積分求解剪力表達式*)

V=Append[Integrate[px, x], c1];

(*通過積分求解彎矩表達式*)

M = Append[Integrate[V, x], c2];

(*通過積分求解轉(zhuǎn)角表達式*)

S = Append[Simplify[ (Integrate[ Flatten[Map[Times[1/EI, #1]&, M]], x ] /. urule[x, Le]) /. HeavisideTheta[x +a_] ->HeavisideTheta[x + a/Le] ] /. HeavisideTheta[x + a_]->HeavisideTheta[x + a *Le], c3];

(*通過積分求解撓度表達式*)

U = Append[Integrate[S, x], c4];

(*取出左端邊界上的集中外力和力偶*)

nhL = {Cases[Partition[PointLoad, 2], {_, 0 | 0.0}],Cases[Partition[PointCouple, 2], {_, 0 | 0.0}]} /.{{{P_, _}} :> P,{} :> 0};

(*根據(jù)左端邊界條件，將x=0代入到所對的條件中，若邊界條件中有彎矩和剪力，將左端的集中外力和力偶加進去*)

bndL = Switch[ First[bnd], “Fixed”, {U, S}, “Pinned”,{U, Append[M, Last[nhL]]},

“Free”, {Append[M, Last[nhL]], Append[V,First[nhL]]}, “Guided”, {S, Append[V, First[nhL]]}] /. x -> 0;

(*取出右端邊界上的集中外力和力偶*)

nhR = {Cases[Partition[PointLoad, 2], {_, Le}],Cases[Partition[PointCouple, 2], {_, Le}]} /.{{{P_, _}} :> P, {}:> 0};

(*根據(jù)右端邊界條件，將x=Le代入到所對的條件中，若右端邊界條件中有彎矩和剪力，將右端的集中外力和力偶加進去*)

bndR = Switch[ Last[bnd], “Fixed”, {U, S}, “Pinned”,{U, Append[M, -Last[nhR]]},

“Free”, {Append[M, -Last[nhR]], Append[V,-First[nhR]]}, “Guided”, {S, Append[V, -First[nhR]]}] /. x ->Le;

(*將中間剛性支撐，鉸鏈連接，彎矩連接的條件，與邊界條件合在一個鏈表中*)

intList = Join[bndR, bndL, evalIntCon[x, U, InSpan],evalIntCon[x, M, ShearConnection], evalIntCon[x, V,MomentConnection]];

(*將左端邊界條件的鏈表中各項加在一起，并等于0，對右端邊界條件，中間剛性支撐，鉸鏈連接，彎矩連接的條件進行同樣操作*)

eqns = Thread[ Apply[Plus, intList, 1] == 0 ] /. urule[x, Le] /.{DiracDelta->0};

(*求解上述命令得出的方程*)

soln = Flatten[ Solve[eqns, Join[ {c1, c2, c3, c4}, First /@Join[InSpan, ShearConnection, MomentConnection]]] ];

(*將求解結(jié)果代回到V, M, S, U的表達式中*)

res = Apply[Plus, {V, M, S, U} /. {DiracDelta->0} /. soln /.urule[x, Le], 1];

(*繪制剪力圖*)

Plot[res[[1]] /. {HeavisideTheta ->UnitStep}, {x, 0,Le},

PlotPoints -> 1500, Frame -> True,

FrameLabel -> {x, “Shear Force”}, LabelStyle -> {12,Bold},

FrameStyle ->Directive[Black, 14, FontFamily -> “Times New Roman”],

Filling -> Axis]

(*繪制彎矩圖*)

Plot[res[[2]] /. {HeavisideTheta ->UnitStep}, {x, 0, Le},

PlotPoints -> 1500, Frame -> True,

FrameLabel -> {x, “Bending Moment”}, LabelStyle -> {12,Bold},

FrameStyle ->Directive[Black, 14, FontFamily -> “Times New Roman”],

Filling -> Axis]

(*程序結(jié)束*)

運行程序，可得如圖6所示的剪力圖和彎矩，與文獻[5]對比，結(jié)果一致。

圖6 剪力圖和彎矩圖

4 結(jié) 論

文中在詳述奇異函數(shù)法理論的基礎(chǔ)上，應用奇異函數(shù)法，將受彎梁上不連續(xù)的載荷寫成統(tǒng)一的載荷函數(shù)，并通過多個實例的計算進行驗證。奇異函數(shù)法演算過程無需求解支座反力，無需受力分析，能夠求解各類彎曲問題，便于編制統(tǒng)一的計算程序。采用Mathematica數(shù)學軟件實現(xiàn)了奇異函數(shù)法算法，并繪制了彎曲問題的剪力、彎矩、轉(zhuǎn)角和撓度等圖形，統(tǒng)一了解決材料力學中各類問題的計算過程，對教與學都有著一定的參考意義。