基于半解析有限元的水泥路面荷載響應分析

沈凱仁，陳先華，周 杰，張含宇

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

車輛荷載是水泥混凝土路面結構的主要外部作用，荷載引起的力學響應是路面結構設計的重要指標[1]。我國現行規范[2]中，水泥混凝土路面荷載應力的計算，仍主要使用彈性地基上的薄板模型；規范中只給出了臨界荷位處的荷載應力計算公式，但未對路面結構荷載應力的分布與變化趨勢進行表述。有限元方法通過對路面結構進行單元劃分，并逐一定義單元屬性，計算各節點的位移和應力，可以有效地反映路面結構的力學響應分布與變化趨勢[3]。因此，有限元方法在路面荷載響應分析中不斷得到應用，現有研究[4－5]證明了有限元是分析路面力學響應的有效工具。主流的有限元程序(如ABAQUS、ANSYS)分析路面荷載響應時，一般采用三維模型進行模擬，其對專業知識和計算成本要求較高，使其局限于研究領域，難以推廣到普遍的工程應用中[6]。因此，有必要結合路面結構的特殊性，對有限元方法進行適當的簡化，開發出既能保證計算精度，又能節約計算成本的有限元程序，并進一步將該程序引入實際工程中，以輔佐現行水泥路面設計規范。

將空間問題簡化為平面問題是提升有限元運算效率的重要方法之一。例如，擋土墻、隧道、管道等結構，其橫截面的幾何尺寸、材料屬性和外力約束一般沿縱向不發生變化，原本的三維有限元問題可以簡化為平面應變問題[7－8]。但對于典型的路面結構，通常各結構層假定為各向同性的均質材料，其橫斷面的幾何尺寸和材料屬性沿行車方向不發生變化，但是外部荷載和邊界條件沿該方向卻有著明顯的變化，因此，路面結構問題不能直接簡化為平面應變問題[9]。軸對稱分析也是簡化有限元模型的重要手段之一，但路面荷載的分布情況通常難以滿足軸對稱模型的要求，因此路面結構問題的軸對稱分析方法也受到了諸多限制[10]。基于路面結構的特殊性，可以在橫斷面上進行二維有限元網格劃分，并通過傅里葉級數表示單元節點位移和路面荷載沿行車方向的分量，利用傅里葉級數表示連續函數的能力，節點位移和荷載沿行車方向的變化均得到了良好的表達。同時，路面結構只在橫截面進行了二維有限元數值分析，荷載和位移沿行車方向的分量采用傅里葉級數的解析表達，在形成整體平衡方程的過程中，由于傅里葉級數的正交性，整體剛度矩陣變化為對角矩陣，原本的三維有限元問題簡化為多個獨立的二維平面問題與解析問題的組合，通過疊加二維平面問題的結果，即可得到路面結構真實的受力狀態，該方法也被稱作半解析有限元方法[11]。而且二維有限元往往只能獲得平面內的應力狀態，但半解析有限元可以通過選取橫斷面縱向坐標，獲得任意橫斷面的應力狀態，因此它的計算功能與三維有限元更為接近。圖1簡要對比了半解析有限元與三維有限元的模型建立方法。

已有的工程應用成果[12－16]，均表明半解析有限元在力學分析中具有良好的準確性與計算效率。Liu等[12－13]通過半解析有限元模擬瀝青路面彈性層狀體系，對比了試驗路的路表彎沉值與半解析有限元的數值解，結果具有良好的一致性。Wu等[14]將超聲波導波理論與半解析有限元相結合，開發了一種計算列車荷載作用下高鐵軌道導波特性的分析程序。Bian等[15－16]通過把Euler-Bernoulli梁模擬的鋼軌和半解析有限元模擬的軌下基礎相結合，分析了荷載作用下高鐵軌下結構的振動響應，分析結果與 ABAQUS的模擬結果吻合良好。基于此，將半解析有限元引入水泥路面荷載響應分析中，開發出準確、高效的有限元分析程序，具有一定的技術可行性和良好的工程應用前景。

圖1 半解析有限元與三維有限元模型示意圖Fig.1 Diagram of the semi-analytical finite element and 3D finite element model

本文基于半解析有限元原理和MATLAB軟件，開發了水泥混凝土路面荷載響應分析程序CP-SAFEM。在下面的章節中，將具體介紹半解析有限元基本理論和CP-SAFEM模型建立的方法，并結合具體工況與現行規范和 ABAQUS進行荷載響應的對比驗證，據此證明CP-SAFEM具有良好的計算精度與運行效率；進一步運用CP-SAFEM確定水凝混凝土路面的臨界荷位，為CP-SAFEM在具體工程中的應用提供一定的參考。

1 水泥路面半解析有限元的力學公式推導與模型建立

1.1 路面結構的半解析有限元模型

水泥路面半解析有限元模型如圖2所示，假定路面結構的幾何尺寸和材料屬性沿z方向保持不變。在此基礎上沿z方向對位移函數和荷載函數進行傅里葉變換，利用傅里葉級數表示連續函數的能力表征位移函數和荷載函數沿z坐標的變化趨勢。

圖2 水泥路面半解析有限元模型Fig.2 Semi-analytical finite element model of cement pavement

本文采用四節點矩形單元進行單元劃分，如圖3所示，與二維有限元不同的是，每個節點需要定義三個方向的位移分量。

圖3 半解析有限元四節點矩形單元Fig.3 The four-node rectangular unit of semi-analytical finite element

節點位移u的形函數表示成z的傅里葉級數，單元內任一點的位移u定義如下：

式中：k為單元節點編號；uk為節點位移；Nk(x,y,z)為對應的形函數；l為第l項傅里葉級數；L為傅里葉級數的總項數；是xy平面內的節點形函數，與二維平面有限元中使用的形函數相同。同理，位移v和w的表達式與u的位移表達式相似。

荷載函數也可以沿z方向進行相應的傅里葉變化，如式(2)。即若原荷載函數為作用平面xz內的面密度，則變換后的荷載函數為假定z坐標下沿x方向的線密度。

對于路面的邊界條件，也可以通過傅里葉變換的奇拓展或偶拓展實現。例如，假定在z=0和z=a的xy平面內僅存在沿z方向的位移，則應對位移u和v進行奇拓展，對位移w進行偶拓展，這樣在邊界平面內僅存在沿z方向的位移，單元內任一點的位移函數可以表示成下式：

考慮到路面結構三維應力狀態下單元應變和位移有如下關系：

將式(3)代入式(4)可得到式(5)。其中，應變-位移矩陣如式(6)。

此外，單元應力與單元應變有如下關系，其中D是各向同性的線彈性材料的本構矩陣。

根據最小勢能原理確定單元對整體平衡方程的貢獻，單元剛度矩陣和荷載向量表示見式(8)~式(9)：

式中：Nl是表面力作用面的形函數矩陣的傅里葉級數的第l項。同時，式(8)將包含下列積分：

由于積分的正交性，可以得到如下關系：

因此，單元剛度矩陣可表達成式(12)：

式中，Ae為單元積分域。通過疊加法組裝整體剛度矩陣和整體荷載向量，根據最小勢能原理獲得路面結構的平衡方程：

式(13)表明，最終的平衡方程簡化為L個獨立的方程，可簡寫為式(14)，其中，K11、F1分別為整體剛度矩陣和荷載向量在傅里葉變換下的第l項。

如果荷載向量的傅里葉變換僅包含特定的某幾項，那么只需求解對應的獨立方程，這進一步節省了計算成本。在得到節點位移的傅里葉級數后，通過傅里葉逆變換可以得到對應三維空間中的節點位移，進一步根據應變-位移、應力-應變關系可以分析整個路面結構的力學響應。

1.2 層間連接和邊界條件的半解析表達

對于普通水泥混凝土路面結構，如圖4所示，選取一塊水泥混凝土面板作為分析對象，考慮橫縫構造和傳力桿的傳荷作用，傳力桿的設置不應妨礙相鄰混凝土板的自由伸縮[17]，因此對z=0和z=a的截面僅固定xy方向的位移，z方向自由伸縮。可以通過式(15)的傅里葉變換滿足上述邊界條件：

圖4 普通水泥混凝土路面結構示意圖Fig.4 Diagram of cement concrete pavement structure

縱縫邊界設置為自由端，同時基層和路基選擇與路面板相同的邊界條件。此外，在路基底部設置豎向位移為零的法向約束。層間接觸面假定為完全連續，即接觸面上、下單元共用節點。

1.3 車輛荷載的半解析表達

依照我國路面結構厚度設計規范，本文采用單軸雙輪 100 kN作為分析荷載，輪胎接地形狀簡化為矩形，荷載位置分為板中、縱縫中部和角隅三種情況，如圖5所示。對荷載項應進行特定的傅里葉變換以確保在Z域中正確模擬荷載情況，荷載函數的傅里葉變換[18]可以寫成：

式中：Pt為荷載壓應力；t為輪載個數；Zt1、Zt2分別為輪載起始和結束的z軸坐標。

圖5 不同荷載位置示意圖Fig.5 Diagram of different load positions

1.4 基于MATLAB的半解析有限元分析程序

本文采用 MATLAB編寫力學響應分析程序。通過發揮 MATLAB矩陣運算的優勢，快速獲得路面力學響應，用于指導層位結構和厚度設計。需要指出，為了避免積分的不可積性，對剛度矩陣和荷載向量分別采用高斯插值求解。同時，在組裝整體剛度矩陣的過程中，由于單元形函數需要從局部坐標系過渡到全局坐標系，因此需要借助雅可比矩陣進行坐標變換。

式中：J為雅可比矩陣；|J|為對應的行列式；Wi、Wj為高斯點對應的權系數。

2 CP-SAFEM各模塊功能介紹

水泥路面半解析有限元程序由前處理模塊、分析模塊和后處理模塊三個部分組成，軟件總體結構如圖6所示，各模塊又由對應的子程序組成，具體功能如下。

2.1 前處理模塊

包括建立路面結構的半解析有限元模型和表達車輛荷載的半解析有限元格式的功能。通過輸入模型尺寸、材料參數和荷載參數，定義單元類型和尺寸，實現模型的網格自動劃分。進一步自動形成節點編號、單元編號和單元對應節點。不同結構層的層間接觸通過相鄰單元共用節點模擬為層間完全連續。根據輸入的荷載形狀與大小，形成在模型空間中的荷載函數。

2.2 分析模塊

包括計算分析路面結構的力學響應的功能。依照有限元單元剛度矩陣和單元荷載向量的形成方法，并沿z坐標對位移形函數和荷載形函數進行傅里葉變換，形成單元剛度矩陣和單元荷載向量在局部坐標系下的傅里葉表達。并根據雅各比矩陣和單元在整體中的位置，實現單元剛度矩陣和單元荷載向量的傅里葉表達從局部坐標系下到全局坐標的轉換。進一步根據直接剛度法和疊加法逐一裝配整體模型中各單元的剛度矩陣和荷載向量。依據最小勢能原理求解傅里葉級數各維度下的節點位移，并通過傅里葉逆變換求解出原三維空間中的節點位移。最后根據位移-應變-應力關系求解出整體模型各項力學響應。

圖6 CP-SAFEM模塊流程圖Fig.6 Flow diagram of CP-SAFEM

2.3 后處理模塊

包括生成各項響應指標的數據表格、繪制力學響應的結構云圖和查詢指定節點的力學響應的功能。后處理模塊主要是對力學響應進行可視化處理。通過指定橫斷面的z坐標，可以生成各項力學響應指標與對應單元節點的數據表格；同時也可以通過指定節點坐標查詢該節點的各項力學指標。并依據生成的數據表格繪制響應指標沿某一坐標的變化趨勢或任意截面的力學響應云圖。

3 有限元程序 CP-SAFEM 的合理性驗證

結合單層板模型和雙層板模型兩種常見的水泥路面結構，對比了規范公式、CP-SAFEM 和ABAQUS的荷載應力求解結果；并針對一種典型路面結構，進行了板中荷位作用下 CP-SAFEM 與ABAQUS的各項應力指標對比，進一步比較了兩者的精算精度和計算效率。

3.1 單層板模型荷載應力對比結果

對于粒料基層上的單層板水泥路面應力分析，規范中采用彈性地基單層板模型[2]。面層板底面以下部分按彈性地基處理，采用地基當量回彈模量Et表征板下結構層的整體彈性特性，計算公式如式(20)~式(23)：

式中：E0/MPa為路基回彈模量；Ex/MPa為粒料層當量回彈模量；α為回歸系數；hx為粒料層總厚度；Ei和hx/m為第i結構層的回彈模量與厚度。

設計軸載在臨界荷位處產生的荷載應力σps(路面板層底最大彎拉應力)按式(24)~式(26)計算：

式中：Ps/kN為設計軸載軸重；Ec/MPa、hc/m和vc分別為混凝土面層板的彈性模量、厚度和泊松比；r/m為路面板的相對剛度半徑；Dc/(MN·m)為混凝土面層板的截面彎曲剛度。

對于典型的單層板路面結構，結構層材料參數如表1所示，路面單元寬度3.5 m，長度4.5 m，荷載選取100 kN，臨界荷位選擇路面板縱縫中部[2]，如圖7所示。

在CP-SAFEM和ABAQUS中，橫截面上采用相同的劃分形式，如表2，對于縱向行車方向，ABAQUS劃分為50個單元，CP-SAFEM選取40項傅里葉級數。荷載選取單軸雙輪100 kN，輪印采用 0.2 m×0.18 m 的矩形模擬，均布壓應力為0.694 MPa。通常，對于典型彈性層狀體系，傅里葉級數維度選擇 40即可保證良好的計算精度[19]，因而本文的傅里葉維度均選擇40。

表1 單層板水泥路面各結構層材料參數Table 1 Material parameters of single-layer cement pavement

圖7 臨界荷位示意圖Fig.7 Diagram of the critical load position

表2 路面橫斷面網格劃分規格Table 2 Meshing division of pavement cross section

三種方法的計算結果如圖8。對比板底最大彎拉應力，規范計算結果為 1.652 MPa，CP-SAFEM計算結果為 1.643 MPa，ABAQUS計算結果為1.577 MPa。CP-SAFEM 與規范計算結果相差0.54%。對比結果表明，CP-SAFEM 的計算結果與規范、ABAQUS具有良好的一致性。另外，CP-SAFEM不僅求解出最大荷載應力，同時也能繪制不同方向的荷載應力的變化趨勢，例如路面板板底的縱向荷載應力遠大于橫向荷載應力[20]，這對路面板結構和配筋設計具有重要意義。

圖8 單層板模型路面層層底拉應力Fig.8 The tensile stress of single-layer plate model pavement

3.2 雙層板模型荷載應力對比結果

對于設置了半剛性基層的水泥路面結構，規范中采用彈性地基上的雙層板模型，半剛性基層以下部分按彈性地基處理，采用地基當量回彈模量Et表征板下結構層的整體彈性特性，分別求解路面板層底和半剛性基層層底的荷載應力。

設計軸載Ps在臨界荷位處產生的路面板層底荷載應力σps(最大彎拉應力)按下式計算：

式中：rg/m為路面板的相對剛度半徑；Dc/(MN·m)為混凝土面層板的截面彎曲剛度；Eb/MPa、hb/m和vb分別表示半剛性基層的彈性模量、厚度和泊松比。半剛性基層層底的荷載應力σbps(最大彎拉應力)按式(30)計算：

對于典型的雙層板路面結構，結構層材料參數如表3所示，路面單元寬度3.5 m，長度4.5 m，在CP-SAFEM和ABAQUS中，定義相同的橫截面網格劃分形式，如表4，對于縱向行車方向，ABAQUS劃分為50個單元，CP-SAFEM選取40項傅里葉級數。荷載形式與單層板模型中相同。

表3 單層板水泥路面各結構層材料參數Table 3 Material parameters of single-layer cement pavement

表4 路面橫斷面網格劃分規格Table 4 Meshing division of pavement cross section

三種方法的計算結果如圖9。對比路面板板底最大彎拉應力，規范計算結果為 1.039 MPa，CP-SAFEM計算結果為0.968 MPa，ABAQUS計算結果為0.947 MPa，CP-SAFEM與規范計算結果相差6.83%。對比半剛性基層板底最大彎拉應力，規范計算結果為0.341 MPa，CP-SAFEM計算結果為0.335 MPa，ABAQUS計算結果為 0.316 MPa，CP-SAFEM與規范計算結果相差1.76%。

圖9 雙層板模型路面層層底拉應力Fig.9 The tensile stress of double-layer plate model pavement

3.3 板中荷位下的荷載應力對比結果

通過 CP-SAFEM 分析板中荷載作用下水泥混凝土路面力學響應，通過與 ABAQUS分別對比道路表面豎向位移、輪跡下豎向壓應力、水泥板底部橫向和縱向拉應力四項指標，以此驗證CP-SAFEM的正確性和可靠性。各結構層材料參數如表5所示，荷載選擇單輪-雙軸荷載，軸重100 kN，輪印寬度0.24 m，輪印高度0.2 m，接地壓力0.521 MPa。路面寬度取4.8 m，長度取5.0 m。在CP-SAFEM和ABAQUS中，定義相同的結構層，橫截面上采用相同的劃分形式，如表6，對于縱向行車方向，ABAQUS劃分為50個單元，CP-SAFEM選取40項傅里葉級數。同時在CP-SAFEM中選擇四節點矩形單元，ABAQUS中選擇八節點六面體單元，模型結構如圖10、圖11所示。

表5 水泥路面各結構層材料參數Table 5 Material parameters of each structural layer of the cement pavement

表6 路面橫斷面網格劃分規格Table 6 Meshing division of pavement cross section

圖10 路面結構ABAQUS剖分模型Fig.10 Pavement structure model of ABAQUS

圖11 路面結構CP-SAFEM剖分模型Fig.11 Pavement structure model of CP-SAFEM

對于典型的水泥路面結構，最大拉應力發生在輪載下水泥板板底，需對水泥板板底橫向和縱向拉應力進行研究；豎向壓應力是保證路基穩定性的重要指標。因此本文選擇上述四項指標作CP-SAFEM與ABAQUS的對比驗證，驗證結果如圖12所示。從圖12中可以發現，CP-SAFEM 的計算結果與ABAQUS均保持良好的一致性，關于各項指標的最大值的具體對比結果見表7。通過對比發現各項指標與 ABAQUS誤差均在 5%以內，說明自編程序CP-SAFEM的計算精度滿足要求。理論上提高傅里葉變換維度，計算精度會進一步提高。

3.4 CP-SAFEM計算效率的驗證

圖12 各項力學響應指標對比圖Fig.12 Comparison of various mechanical indicators

為了評估CP-SAFEM的計算效率，CP-SAFEM和ABAQUS均在intel core i5-8500 CPU、8G內存上運行，單元、節點個數和計算時間如表8所示。通過對比可知，CP-SAFEM 需要的單元個數僅為ABAQUS的1/50，計算時間僅為ABAQUS的1/6，如果開啟 MATLAB多線程優化，計算時間將會進一步縮減[21]。

表7 各項應力指標最大值對比Table 7 Maximum comparison of various stress indicators

表8 CP-SAFEM與ABAQUS單元、節點和運行時間對比Table 8 Comparison of CP-SAFEM and ABAQUS units,nodes and operation time

通過本節分析可以發現，半解析有限元不僅在計算精度上基本達到 ABAQUS水準，同時利用傅里葉變換和正交函數的優勢，可以極大地減少計算成本。隨著分析對象結構尺寸的增加和材料復雜程度的加深，這一優勢將變得更加突出。

3.5 不同荷位下路面板彎拉應力對比分析

在水泥混凝土路面路面荷載應力分析中，通常選擇使面層板產生最大應力的荷載位置作為應力計算的臨界荷位，臨界荷位一般發生在水泥路面板板底。現通過半解析有限元程序CP-SAFEM，采用3.4節水泥混凝土路面結構模型，縱縫邊界不考慮拉桿，分別考慮板中、縱縫中部、角隅三種荷位，作路面板層底橫向和縱向拉應力的對比分析，據此尋找臨界荷位。

在板中、縱縫中部、角隅荷位下的橫向彎拉應力分別如圖13所示，每一單元大小為0.12 m×0.02 m，為了使云圖更為細致地表現應力狀態，豎向尺寸放大了 10倍。路面板上部表現為壓應力，最大彎拉應力出現在面層層底。其中，對于板中和板邊荷位，輸出應力截面為z=2.5 m，角隅輸出截面為z=4.9 m。類似地，縱向彎拉應力如圖14所示。

綜合考慮不同荷位下橫縱向彎拉應力最大值，結果如表9所示。最大彎拉應力出現在板邊荷位下的縱向拉應力，因此選擇板邊荷位作為臨界荷載。對于角隅荷位，由于傳力桿的約束作用，彎拉應力相較于其他荷位要小很多。

圖13 不同荷位下路面板橫向彎拉應力/MPaFig.13 Transverse bending stress of road panel under different loading positions

圖14 不同荷位下路面板縱向彎拉應力/MPaFig.14 Longitudinal bending stress of road panel under different loading positions

表9 不同荷位下的最大彎拉應力Table 9 Maximum bending stress at different load levels

依據上述分析，對于橫縫設置傳力桿的水泥混凝土路面，臨界荷位應選擇板邊荷位，這與文獻[22]中水泥路面荷載應力的臨界荷位分析結果一致。最大拉應力為左輪荷載下路面板底部縱向拉應力。

4 結論

本文介紹了基于半解析有限元方法的水泥混凝土路面力學響應計算理論，并根據相關理論和MATLAB開發了半解析有限元程序 CP-SAFEM，結合具體工況驗證了 CP-SAFEM 具有良好的準確性與計算效率。

(1)針對單層板和雙層板兩種典型的水泥路面結構，按規范要求選擇縱縫中部荷位，對比了規范、CP-SAFEM和ABAQUS的荷載應力計算結果，單層板和雙層板下層的計算誤差約為 1%，雙層板上層計算誤差約為6%。

(2)對于板中荷位下的水泥路面，CP-SAFEM和ABAQUS各項力學指標誤差均在5%以內，進一步驗證了CP-SAFEM的準確性。

(3)針對同一路面結構，采用相同的橫截面網格劃分形式，對比了CP-SAFEM和ABAQUS的模型復雜程度和計算時間，CP-SAFEM模型的單元數目僅為ABAQUS的1/50，運行時長僅為1/6，證明了CP-SAFEM良好的計算效率。

(4)通過CP-SAFEM對比了不同荷位下水泥混凝土路面的荷載應力，發現荷位為縱縫中部時出現最大荷載應力，證明了臨界荷位應選擇縱縫中部，這與現有研究結論相同。