層狀彈性體系力學及其應用

王凱

(東南大學交通學院，南京210096)

層狀彈性體系力學是多層柔性路面、機場道面以及多層地基設計與計算的理論基礎。該理論及在其基礎上所編制的各種實用程序在世界各國的上述工程結構的科學研究和設計與計算中得到了廣泛的應用。

層狀彈性體系力學是彈性力學的一個分支。層狀彈性體系力學將所研究的物體看作是自上而下由若干彈性層和半無限彈性體組成的彈性體系。

盡管層狀彈性體系力學屬于彈性力學的范疇，但與傳統的彈性力學相比，層狀彈性體系力學有很大的不同。傳統的彈性力學只能解決半無限彈性體的力學分析和計算問題，而層狀彈性體系力學可以解決半無限層狀彈性體系的力學分析和計算問題。在層狀彈性體系力學中，經常采用漢克爾(Hankel)積分變換公式并頻繁使用各種特殊函數(伽馬(gamma)函數、貝塔(beta)函數、普西(psi)函數、橢圓積分、超幾何函數、貝塞爾(Bessel)函數、勒讓德(Legendre)函數、拉姆達(lambda)函數)的計算公式，而這兩種計算公式在傳統的彈性力學的力學分析與計算中很少出現或幾乎不出現。因此層狀彈性體系力學大大豐富了彈性力學的理論寶庫。

目前在我國，層狀彈性體系力學在交通工程界已得到廣泛的應用，而在力學界則鮮為人知。本文是迄今國內第一次以論文的形式全面、系統地向國內廣大讀者宣傳、介紹層狀彈性體系力學的發展歷史、基礎理論、計算公式和計算方法及其在科學研究和工程設計中的應用。盡管層狀彈性體系力學誕生已經幾十年了，文中的內容對于國內很多讀者來說，仍然是十分新鮮并具有理論價值和實用價值的知識，值得一讀。

1 層狀彈性體系力學的發展歷史

層狀彈性體系力學是在半無限彈性體力學的基礎上發展起來的。1885年，Boussinesq對半無限彈性體表面在單個垂直集中力作用下的應力和位移作出了理論解，它在近代土力學中獲得了廣泛的應用。1882年，Cerruti對半無限彈性體表面在單個水平集中力作用下的應力和位移作出了理論解。

由于數學和彈性力學的發展，到20世紀40～60年代，層狀彈性體系力學取得了長足的進步。1945年，Burmister[1]利用Love位移函數得到了在軸對稱垂直載荷作用下雙層和多層彈性體系應力和位移的理論解。1962年，Schiffman[2]又將其進一步推廣到多層彈性體系非軸對稱問題的求解。

20世紀60～70年代，電子計算機及計算方法發展很快，而工程實際也要求解決多層(層數N>3)彈性體系的力學計算，于是多層彈性體系力學計算程序在各國學者的努力下應運而生。

1973年，De Jong等[3]合作編制成功BISAR(bitumen structures analysis in roads)計算機程序，該程序可以計算在多圓均布復合載荷(包括垂直、單向水平載荷)作用下N層彈性連續-半結合體系內任一點的應力和位移分量、主應力和主應變分量以及主方向等。

除此之外，世界各國還有不少計算N層彈性體系應力和位移分量的計算機程序，但迄今為止世界上公認比較有名、實用且影響比較深遠的仍然是BISAR程序。

改革開放以來，隨著交通事業的發展，高等級公路和城市道路瀝青路面的大量設計與修建，我國迫切需要解決多層彈性體系的力學計算問題。但國外在這一方面對我國搞專利封鎖，著名的BISAR程序專利費高達100萬美元。

為了打破國外的專利封鎖，中國學者們決心攻克這一國內難題。1980年筆者編制成功了在圓形均布垂直載荷作用下N層彈性連續體系的力學計算程序，并從此開始將筆者本人編寫的N層彈性體系力學計算程序統稱為NESCP(N-layer elastic system computer program)程序[4]。1981年筆者又編制成功了在雙圓均布復合載荷(垂直和單向水平載荷)作用下N層彈性連續體系的力學計算程序[5]。在此基礎上筆者經過數年努力，于1984年編制了功能更全面的多層彈性體系力學計算程序。該程序的功能已達到BISAR程序的功能[6-9]。1990年，由筆者完成的“層狀彈性體系應力與位移的一般分析與計算”課題獲陜西省科技進步一等獎。

除了筆者之外，我國還有不少學者按照不同思路編制了若干多層彈性體系力學計算程序，較有影響的有朱照宏的“高階矩陣代數法”程序[10]，郭文復的“分層求逆子陣的傳遞矩陣法”程序[11]，郭大智的“系數遞推法”程序[12]以及吳晉偉的“反力遞推法”程序[13]等。

應當指出，由于國外技術封鎖，我國多層彈性體系力學分析與計算的研究工作，基本上是在自力更生的基礎上進行的，但這也促使我們的研究在很多方面具有中國特色。

由多層彈性體系力學計算可知，如何從線性代數方程組快速地求算出積分常數，是保證快速計算出應力與位移分量的關鍵。在這一方面，中外的學者們提出了不同的方法，一種是筆者提出的“遞推回代法”，另一種是“傳遞矩陣法”(由于保密，國外文獻中從不介紹該計算方法的具體細節)。計算表明，“遞推回代法”的計算時間少，不會出現“病態矩陣”，對層間完全連續、層間完全光滑和層間半結合時積分常數的求解問題都能解決，比“傳遞矩陣法”(該方法不能解決層間完全光滑時積分常數的求解問題)明顯優越。

積分計算公式是多層彈性體系力學計算的另一個關鍵計算公式。由于保密，國外文獻中從不介紹積分計算公式的內容。筆者在前人工作的基礎上，通過數學推導和誤差分析，得到了論文中所述的積分計算公式并成功地應用于多層彈性體系的力學計算。

由于本文篇幅有限，有關內容的詳細介紹，請參看參考文獻[14]。

以上僅對層狀彈性體系力學的發展歷史作了簡單的回顧，下文將進一步介紹N層彈性體系的力學分析與計算。

2 N層彈性體系的力學分析與計算[14]

層狀彈性體系力學包含了豐富的內容，由于篇幅有限，本文只著重介紹在工程實踐中應用最廣泛的在圓形均布垂直載荷或單向水平載荷作用下N層彈性體系的力學分析與計算并在介紹過程中對NESCP程序和BISAR程序的力學計算公式和計算方法進行比較。所謂N層彈性體系通常指自上而下由N?1層有限厚的水平彈性層和最下層半無限彈性體組成的體系。當N=1時，層狀彈性體系即變為半無限彈性體。N層彈性體系在單圓均布垂直載荷或單向水平載荷作用下的力學計算簡圖如圖1所示。圖中q和p分別為垂直載荷和水平載荷集度，δ為載荷圓半徑，hi，Ei和μi分別為彈性體系各層的厚度、彈性模量和泊松比。

2.1 基本假定

為了進行層狀彈性體系的力學分析，在層狀彈性體系力學中關于層狀彈性體系，引入如下幾個基本假定：

(1)各層是連續的、完全彈性的、均勻的和各向同性的理想彈性體；

(2)最下層半無限彈性體在水平方向和垂直向下方向均為無限大，其上各層垂直方向具有有限的厚度，水平方向為無限大；

(3)各層水平無限遠處和最下層無限深處，應力和位移分量為零；

(4)層間結合狀況可以是完全連續的，或者是完全光滑的，或者是介于二者之間的半結合狀況，但層間不出現脫空現象；

(5)體系受力后位移和變形是微小的，在進行力學分析時體力可忽略不計。

2.2 應力和位移分量的表達式

采用位移函數法和漢克爾積分變換的公式[15]求解柱坐標系中彈性力學三維問題的基本方程，即可得到N層彈性體系在單圓均布垂直載荷或單向水平載荷作用下應力和位移分量的表達式：

當載荷為垂直載荷時

其中

當載荷為單向水平載荷時

其中

以上各式中，J0(x)，J1(x)和J2(x)分別為零階、一階和二階第一類貝塞爾函數，下標i表示各力學分量計算點所在的層數；Ai，Bi，Ci，Di，Fi，Ii為積分常數，實際為積分變量x的函數。如下所述，它們將根據具體力學問題的定解條件來確定，其他各變量解釋詳見參考文獻[14]，下同。

2.3 定解條件

2.3.1 NESCP程序

當載荷為垂直載荷時

(1)表面應力邊界條件(z=0)

(2)層間結合條件(z=zi)

此外，定解條件還有：當z→∞時，所有應力和位移分量均應趨于零，因此必有AN=CN=0。

當載荷為單向水平載荷時

(1)表面應力邊界條件(z=0)

(2)層間結合條件(z=zi)

此外，定解條件還有AN=CN=FN=0，理由同前。

2.3.2 BISAR程序

BISAR程序進行力學計算所采用的定解條件與上述NESCP程序所采用的定解條件基本一致，僅式(4)和式(6)中水平位移與水平剪應力的關系式分別修改為

式(7)中AKi稱為界面剪切彈簧系數，0≤AKi<∞(i=1,2,···,N?1)。當AKi=0時，ui=ui+1vi=vi+1，層間完全連續；當AKi→∞時，τzri=τzri+1=0，τzθi=τzθi+1=0，層間完全光滑；當0

由上述可知：

(1)當層間處于半結合狀況時，NESCP程序中采用的水平位移與水平剪應力的關系式是嚴格按照哥德曼(Goodman)法則的表達式，即

得到的，而BISAR程序中采用的關系式仍按照哥德曼法則的表達式但做了一點變化。

(2)在具體力學計算時，界面剪切彈簧系數AKi與層間結合系數Ki實際上是互為倒數關系。

2.4 根據定解條件建立求解積分常數的線性代數方程組

由于本文篇幅有限，下文僅介紹垂直載荷作用下N層彈性體系力學計算的有關內容。單向水平載荷作用下N層彈性體系力學計算的有關內容與此基本相似，請參看參考文獻[14]，不再一一敘述。將式(1)中有關應力分量的表達式代入表面應力邊界條件式(3)并對每一個等式兩邊施加漢克爾積分變換[15]即可得到兩個方程式(寫成矩陣形式)

將式(1)中有關力學分量的表達式代入層間結合條件式(4)并對每一個等式兩邊施加漢克爾積分變換[15]即可得到(4N?4)個方程式(寫成矩陣形式)

上述4N?2個方程式再加上AN=CN=0兩個方程式，正好構成4N元線性代數方程組，用于求解4N個未知積分常數Ai，Bi，Ci，Di(i=1,2,···,N)。

2.5 由線性代數方程組求解積分常數

由多層彈性體系力學計算可知，如何從上述線性代數方程組快速地求算出積分常數，是保證快速計算出應力與位移分量數值的關鍵。在這一方面，中外學者們提出了不同的方法，一種是筆者提出的“遞推回代法”并用于NESCP程序，另一種是“傳遞矩陣法”，為BISAR程序所采用，現分別敘述如下。

2.5.1 用“遞推回代法”求解積分常數

由前述可知，在垂直載荷作用下N層彈性體系求解積分常數的4N元線性代數方程組可以分成若干小組，其中由表面應力邊界條件得到的兩個方程式構成一個小組，由每一個層間界面上的四個層間結合條件得到的四個方程式分別構成N?1個小組。如果能利用這一特性建立相鄰結構層積分常數的遞推關系式，則4N元線性代數方程組的求解問題就可以簡化成若干個四元乃至二元線性代數方程組的求解問題，從而大大簡化計算。筆者提出的“遞推回代法”正是基于這一思路，而“傳遞矩陣法”也是基于這一思路。“遞推回代法”的具體過程，簡單來說就是先以各個小組為單位，將式(8)和式(9)內的Bi，Di表示成以Ai，Ci(i=1,2,···,N)為自變量的多項式函數表達式

其次，通過式(10)和式(11)中i=1的方程式小組以及式(11)中相鄰層間界面的方程式小組Bi和Di相等的關系，消去Bi，Di(i=1,2,...,N?1)，即可進一步建立Ai，Ci與Ai+1，Ci+1之間的遞推關系式

利用該遞推關系式以及AN=CN=0，從i=N?1開始通過自下而上逐層遞推，求出所需的Ai和Ci，這就是遞推過程；再將已求得的Ai和Ci回代到Bi和Di的多項式函數表達式中，求出所需的Bi和Di，這就是回代過程。

2.5.2用“傳遞矩陣法”求解積分常數

現將“傳遞矩陣法”的求解過程(仍以垂直載荷作用下N層彈性體系為例)敘述如下：

由式(9)并記[Ti]4×4=[Ki]?14×4[Ki+1]4×4(式中矩陣[Ti]稱為傳遞矩陣)，則可得到

又記[T]=[T1][T2]···[TN?1]并利用AN=CN=0，即可得到第一層積分常數與第N層積分常數之間的關系式為

將式(14)代入式(8)，得

求解時，首先利用式(15)確定BN和DN，然后利用式(13)和AN=CN=0，從第N層起逐層向上計算，直至計算出所需的Ai，Bi，Ci和Di為止。

2.5.3 “遞推回代法”與“傳遞矩陣法”的比較

首先，由于“傳遞矩陣法”的求算過程仍基于矩陣代數方法，在積分計算過程中需要頻繁進行矩陣求逆或矩陣相乘運算，費時頗多；而“遞推回代法”的求算過程完全成了普通的代數運算，比“傳遞矩陣法”省去了很多中間計算過程，從而大大節省了計算時間。其次，“傳遞矩陣法”不能解決層間完全光滑時積分常數的求解問題。這是由于層間完全光滑時式(9)中的系數矩陣[Ki]為奇異矩陣(有一行元素全部為零)，其對應的逆矩陣[Ki]?1不存在。而“遞推回代法”對層間完全連續、層間完全光滑和層間半結合時積分常數的求解問題都能解決。第三，由BISAR程序的計算可知，采用“傳遞矩陣法”進行矩陣運算時，有時會出現“病態矩陣”并導致計算結果很不精確。而“遞推回代法”不會出現“病態矩陣”。

2.6 貝塞爾函數計算

由式(1)和式(2)可知，N層彈性體系應力與位移分量積分表達式的被積函數中包含貝塞爾函數和指數函數多項式兩個部分。通過上述積分常數的求解和計算，被積函數中指數函數多項式部分已完全可以計算。下面再進一步介紹被積函數中貝塞爾函數的計算。在NESCP程序和BISAR程序中一般采用計算方法中的逼近多項式和逼近公式計算貝塞爾函數，前者采用阿倫(Allen)逼近公式，后者采用切比雪夫(Chebyshev)逼近公式，二者的計算精度和計算時間基本一致，前者稍佳。由于本文篇幅有限，下面僅介紹阿倫逼近公式。在NESCP程序中計算第一類零階和一階貝塞爾函數J0(x)和J1(x)的計算公式，如式(16)和式(17)所示。

當x≤3時，令t=(x/3)2，J0(x)和J1(x)的逼近多項式為

當x>3時，令t=3/x，J0(x)和J1(x)的逼近公式為

其中A0，B0，A1，B1采用的逼近多項式為

2.7 積分計算公式和積分上限計算公式

由式(1)和式(2)可知，N層彈性體系的應力與位移分量表達式都是含貝塞爾函數和指數函數的復雜被積函數的無窮積分，一般只能通過數值積分并采用有限數值的積分上限進行計算。在前人工作的基礎上，筆者通過數學推導和誤差分析，得到了如下述的積分計算公式(為簡單明了起見，未寫出應力與位移分量表達式積分號前的系數和無窮積分式大方括號外面的三角函數)和積分上限xs的計算公式并用于NESCP程序。

(1)積分計算公式

當z≥h1時

當z

以上二式中，z為計算點的z坐標，h1為多層彈性體系第一層的厚度，J(x)表示各力學分量積分表達式中被積函數的貝塞爾函數部分，E(x)表示被積函數的指數函數多項式部分，e(x)是當彈性體為半無限彈性體時E(x)的表達式，積分和采用高斯(Gauss)求積公式，利用數值積分計算出，而積分則采用列普司切茲?漢克爾(Lipschitz-Hankel)積分公式計算出[16]。

(2)積分上限計算公式

其中δ為載荷圓半徑，z，h1的意義同式(19)。

試算表明，當積分上限xs取式(20)計算且在0～xs積分區間內設置足夠多的數值積分結點時，即可滿足計算結果至少精確到小數點后四位數。

由于保密，國外文獻中從不介紹積分計算公式的內容，但根據對有關資料和對比算例的綜合分析可知，BISAR程序的積分計算公式與上述NESCP程序基本一致，主要不同之處在于NESCP程序的積分上限xs是由式(20)通過計算確定的，而BISAR程序的積分上限xs是數值積分過程中根據計算精度分析確定的，見下文2.8節數值積分中所述。

2.8 數值積分

在NESCP程序中筆者采用高斯?勒讓德求積公式對各應力與位移分量表達式進行數值積分，在利用式(20)計算出積分上限xs后，首先將積分區間[0，xs]劃分成長度相等的m個子區間，每個子區間上采用34結點高斯?勒讓德求積公式計算其積分值，然后再將m個子區間的積分值相加，得到整個積分區間上的積分值。試算表明，在通常情況下當使用由式(20)得到的積分上限xs且m取3～5時，應力與位移分量的計算結果與國內外已公開發表的計算結果或者完全一致，或者十分接近并且滿足計算結果至少精確到小數點后四位數。

在BISAR程序中各應力與位移分量表達式的數值積分過程為：

由式(1)和式(2)可知，各力學分量表達式都含有貝塞爾函數，貝塞爾函數是一個波動衰減函數，在積分區間[0，xs]內有很多零點。BISAR程序中的數值積分是在原點(x=0)到第一個貝塞爾函數乘積的零點(x=z1)以及在后續的一個貝塞爾函數乘積的零點(x=zn)到下一個貝塞爾函數乘積的零點(x=zn+1)(n=1,2,···)之間的區間上進行的。在第一個區間[0，z1]上，采用8結點高斯?勒讓德求積公式計算區間的積分值。從第二個區間開始，采用4～15結點高斯?洛巴多(Lobatto)求積公式計算每一個區間的積分值。設INT為累計積分值，在積分計算過程中，如果連續兩個貝塞爾函數乘積的零點之間的區間上積分的絕對值小于10?5(當|INT|<0.01)或區間上的積分值與INT之比的絕對值小于10?4(當|INT|≥0.01)，則積分計算結束。

當計算積分常數的矩陣運算過程中出現病態矩陣，數值積分會過早地停止并顯示中止的原因。

由上述可知，NESCP程序數值積分過程簡單明了，而BISAR程序則細致繁復，這大大增加了數值積分計算的時間。除此之外，在BISAR程序積分過程中，當由于出現病態矩陣使得數值積分過早停止時，將導致數值積分的結果很不精確。

2.9 后續的力學計算

以上內容介紹了在單圓載荷作用下N層彈性體系的應力與位移計算，在BISAR程序和NESCP程序中，后續的力學計算還有如下幾項。

(1)將各單圓載荷在計算點處產生的柱坐標中的應力與位移分量變換成公共直角坐標中的應力與位移分量，再通過求和得到多圓載荷作用下在計算點處產生的公共直角坐標中的總應力與總位移分量并利用彈性力學公式進一步計算出公共直角坐標中的總應變分量；

(2)計算主應力及其方向余弦和主應變；

(3)計算極值剪應力、極值剪應變及其相應正應力以及它們的方向余弦；

(4)計算彈性體內每單位體積的總應變能和畸變應變能。

由于本文篇幅有限，上述后續力學計算的計算公式和計算過程請參看參考文獻[14,17-20]的有關論述，本文不再一一介紹。

2.10 計算算例

計算算例是一個如圖1所示的N等于5的五層彈性體系，表面承受單圓均布垂直載荷或單向水平載荷的作用。由于本文篇幅有限，下文僅列出該算例用NESCP程序和BISAR程序分別計算的在單圓均布垂直載荷或單向水平載荷作用下，各計算點處產生的柱坐標應力與位移分量。

(1)五層彈性體系的結構參數

五層彈性體系的結構參數如表1所示。

(2)載荷計算參數

載荷計算參數如表2所示。

表2 載荷計算參數

(3)數值計算結果

第一個載荷即垂直載荷單獨作用時數值計算結果如表3所示(表中從左至右依次為計算點所在的層數i，計算點的x，y，z坐標，柱坐標應力與位移分量。每個計算點的計算數據的第一行為NESCP程序的計算結果，第二行為BISAR程序的計算結果)

表3 垂直載荷作用下計算點處產生的柱坐標應力與位移分量

第二個載荷即單向水平載荷單獨作用時數值計算結果如表4所示(表中從左至右依次為計算點所在的層數i，計算點的x，y，z坐標，柱坐標應力與位移分量。每個計算點的計算數據的第一行為NESCP程序的計算結果，第二行為BISAR程序的計算結果)。

由兩個程序的對比計算結果可知，絕大多數對比數據完全一致，少數對比數據雖不完全一致，但十分接近(小數點后至少前4位數據能保持一致)。

3 層狀彈性體系力學在科學研究和工程設計中的應用

如本文開頭所述，層狀彈性體系力學是多層柔性路面、機場道面以及多層地基設計與計算的理論基礎。該理論及在其基礎上所編制的各種實用程序在世界各國的上述工程結構的科學研究以及設計與計算中得到了廣泛的應用。

以BISAR程序為例，英荷殼牌石油公司在最初編寫BISAR程序時，對BISAR程序的使用功能作了兩點定位：一是多層彈性體系通用力學計算程序，二是為殼牌路面設計方法SPDM(shell pavement design mothod)編制路面設計圖表。自1973年BISAR程序問世以來，該程序為世界各國的路面與機場道面的科學研究和工程設計與計算做出了重大的貢獻。

同樣筆者編寫的NESCP程序從20世紀80年代開始，為我國的路面與機場道面的科學研究也做出了很大的貢獻，先后計算了幾十萬個數據。筆者在NESCP程序的基礎上編制了多個公路和城市道路瀝青路面設計程序，服務于我國交通建設事業和高校教學。到目前為止，這些程序已在全國幾百個設計單位和高校得到了成功應用，取得了顯著的社會經濟效益。

4 結語

(1)本文是迄今國內第一次以論文的形式全面、系統地向國內廣大讀者宣傳、介紹層狀彈性體系力學及其應用，讓層狀彈性體系力學發揚光大并得到新的更大的發展。層狀彈性體系力學包含了豐富的內容而本文篇幅有限，為了比較深入和全面地了解層狀彈性體系力學，筆者建議有興趣的讀者可以進一步閱讀參考文獻[14]。

(2)BISAR程序是當今世界上最負盛名的多層彈性體系力學計算程序。本文第二節對筆者編寫的NESCP程序和BISAR程序的力學計算公式和計算方法進行了系統的對比。對比結果表明，NESCP程序的功能不僅可以與BISAR程序相媲美，而且在某些方面甚至更優越。本文用事實證明，外國人能做到的，中國人通過努力也一定能做到，而且還可能做得更好。