基于半解析有限元的水泥路面荷載響應(yīng)分析

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

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

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

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

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

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

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

1 水泥路面半解析有限元的力學(xué)公式推導(dǎo)與模型建立

1.1 路面結(jié)構(gòu)的半解析有限元模型

水泥路面半解析有限元模型如圖2所示，假定路面結(jié)構(gòu)的幾何尺寸和材料屬性沿z方向保持不變。在此基礎(chǔ)上沿z方向?qū)ξ灰坪瘮?shù)和荷載函數(shù)進(jìn)行傅里葉變換，利用傅里葉級數(shù)表示連續(xù)函數(shù)的能力表征位移函數(shù)和荷載函數(shù)沿z坐標(biāo)的變化趨勢。

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

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

圖3 半解析有限元四節(jié)點(diǎn)矩形單元Fig.3 The four-node rectangular unit of semi-analytical finite element

節(jié)點(diǎn)位移u的形函數(shù)表示成z的傅里葉級數(shù)，單元內(nèi)任一點(diǎn)的位移u定義如下：

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

荷載函數(shù)也可以沿z方向進(jìn)行相應(yīng)的傅里葉變化，如式(2)。即若原荷載函數(shù)為作用平面xz內(nèi)的面密度，則變換后的荷載函數(shù)為假定z坐標(biāo)下沿x方向的線密度。

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

考慮到路面結(jié)構(gòu)三維應(yīng)力狀態(tài)下單元應(yīng)變和位移有如下關(guān)系：

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

此外，單元應(yīng)力與單元應(yīng)變有如下關(guān)系，其中D是各向同性的線彈性材料的本構(gòu)矩陣。

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

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

由于積分的正交性，可以得到如下關(guān)系：

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

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

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

如果荷載向量的傅里葉變換僅包含特定的某幾項(xiàng)，那么只需求解對應(yīng)的獨(dú)立方程，這進(jìn)一步節(jié)省了計(jì)算成本。在得到節(jié)點(diǎn)位移的傅里葉級數(shù)后，通過傅里葉逆變換可以得到對應(yīng)三維空間中的節(jié)點(diǎn)位移，進(jìn)一步根據(jù)應(yīng)變-位移、應(yīng)力-應(yīng)變關(guān)系可以分析整個路面結(jié)構(gòu)的力學(xué)響應(yīng)。

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

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

圖4 普通水泥混凝土路面結(jié)構(gòu)示意圖Fig.4 Diagram of cement concrete pavement structure

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

1.3 車輛荷載的半解析表達(dá)

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

式中：Pt為荷載壓應(yīng)力；t為輪載個數(shù)；Zt1、Zt2分別為輪載起始和結(jié)束的z軸坐標(biāo)。

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

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

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

式中：J為雅可比矩陣；|J|為對應(yīng)的行列式；Wi、Wj為高斯點(diǎn)對應(yīng)的權(quán)系數(shù)。

2 CP-SAFEM各模塊功能介紹

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

2.1 前處理模塊

包括建立路面結(jié)構(gòu)的半解析有限元模型和表達(dá)車輛荷載的半解析有限元格式的功能。通過輸入模型尺寸、材料參數(shù)和荷載參數(shù)，定義單元類型和尺寸，實(shí)現(xiàn)模型的網(wǎng)格自動劃分。進(jìn)一步自動形成節(jié)點(diǎn)編號、單元編號和單元對應(yīng)節(jié)點(diǎn)。不同結(jié)構(gòu)層的層間接觸通過相鄰單元共用節(jié)點(diǎn)模擬為層間完全連續(xù)。根據(jù)輸入的荷載形狀與大小，形成在模型空間中的荷載函數(shù)。

2.2 分析模塊

包括計(jì)算分析路面結(jié)構(gòu)的力學(xué)響應(yīng)的功能。依照有限元單元剛度矩陣和單元荷載向量的形成方法，并沿z坐標(biāo)對位移形函數(shù)和荷載形函數(shù)進(jìn)行傅里葉變換，形成單元剛度矩陣和單元荷載向量在局部坐標(biāo)系下的傅里葉表達(dá)。并根據(jù)雅各比矩陣和單元在整體中的位置，實(shí)現(xiàn)單元剛度矩陣和單元荷載向量的傅里葉表達(dá)從局部坐標(biāo)系下到全局坐標(biāo)的轉(zhuǎn)換。進(jìn)一步根據(jù)直接剛度法和疊加法逐一裝配整體模型中各單元的剛度矩陣和荷載向量。依據(jù)最小勢能原理求解傅里葉級數(shù)各維度下的節(jié)點(diǎn)位移，并通過傅里葉逆變換求解出原三維空間中的節(jié)點(diǎn)位移。最后根據(jù)位移-應(yīng)變-應(yīng)力關(guān)系求解出整體模型各項(xiàng)力學(xué)響應(yīng)。

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

2.3 后處理模塊

包括生成各項(xiàng)響應(yīng)指標(biāo)的數(shù)據(jù)表格、繪制力學(xué)響應(yīng)的結(jié)構(gòu)云圖和查詢指定節(jié)點(diǎn)的力學(xué)響應(yīng)的功能。后處理模塊主要是對力學(xué)響應(yīng)進(jìn)行可視化處理。通過指定橫斷面的z坐標(biāo)，可以生成各項(xiàng)力學(xué)響應(yīng)指標(biāo)與對應(yīng)單元節(jié)點(diǎn)的數(shù)據(jù)表格；同時(shí)也可以通過指定節(jié)點(diǎn)坐標(biāo)查詢該節(jié)點(diǎn)的各項(xiàng)力學(xué)指標(biāo)。并依據(jù)生成的數(shù)據(jù)表格繪制響應(yīng)指標(biāo)沿某一坐標(biāo)的變化趨勢或任意截面的力學(xué)響應(yīng)云圖。

3 有限元程序 CP-SAFEM 的合理性驗(yàn)證

結(jié)合單層板模型和雙層板模型兩種常見的水泥路面結(jié)構(gòu)，對比了規(guī)范公式、CP-SAFEM 和ABAQUS的荷載應(yīng)力求解結(jié)果；并針對一種典型路面結(jié)構(gòu)，進(jìn)行了板中荷位作用下 CP-SAFEM 與ABAQUS的各項(xiàng)應(yīng)力指標(biāo)對比，進(jìn)一步比較了兩者的精算精度和計(jì)算效率。

3.1 單層板模型荷載應(yīng)力對比結(jié)果

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

式中：E0/MPa為路基回彈模量；Ex/MPa為粒料層當(dāng)量回彈模量；α為回歸系數(shù)；hx為粒料層總厚度；Ei和hx/m為第i結(jié)構(gòu)層的回彈模量與厚度。

設(shè)計(jì)軸載在臨界荷位處產(chǎn)生的荷載應(yīng)力σps(路面板層底最大彎拉應(yīng)力)按式(24)~式(26)計(jì)算：

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

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

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

表1 單層板水泥路面各結(jié)構(gòu)層材料參數(shù)Table 1 Material parameters of single-layer cement pavement

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

表2 路面橫斷面網(wǎng)格劃分規(guī)格Table 2 Meshing division of pavement cross section

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

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

3.2 雙層板模型荷載應(yīng)力對比結(jié)果

對于設(shè)置了半剛性基層的水泥路面結(jié)構(gòu)，規(guī)范中采用彈性地基上的雙層板模型，半剛性基層以下部分按彈性地基處理，采用地基當(dāng)量回彈模量Et表征板下結(jié)構(gòu)層的整體彈性特性，分別求解路面板層底和半剛性基層層底的荷載應(yīng)力。

設(shè)計(jì)軸載Ps在臨界荷位處產(chǎn)生的路面板層底荷載應(yīng)力σps(最大彎拉應(yīng)力)按下式計(jì)算：

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

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

表3 單層板水泥路面各結(jié)構(gòu)層材料參數(shù)Table 3 Material parameters of single-layer cement pavement

表4 路面橫斷面網(wǎng)格劃分規(guī)格Table 4 Meshing division of pavement cross section

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

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

3.3 板中荷位下的荷載應(yīng)力對比結(jié)果

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

表5 水泥路面各結(jié)構(gòu)層材料參數(shù)Table 5 Material parameters of each structural layer of the cement pavement

表6 路面橫斷面網(wǎng)格劃分規(guī)格Table 6 Meshing division of pavement cross section

圖10 路面結(jié)構(gòu)ABAQUS剖分模型Fig.10 Pavement structure model of ABAQUS

圖11 路面結(jié)構(gòu)CP-SAFEM剖分模型Fig.11 Pavement structure model of CP-SAFEM

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

3.4 CP-SAFEM計(jì)算效率的驗(yàn)證

圖12 各項(xiàng)力學(xué)響應(yīng)指標(biāo)對比圖Fig.12 Comparison of various mechanical indicators

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

表7 各項(xiàng)應(yīng)力指標(biāo)最大值對比Table 7 Maximum comparison of various stress indicators

表8 CP-SAFEM與ABAQUS單元、節(jié)點(diǎn)和運(yùn)行時(shí)間對比Table 8 Comparison of CP-SAFEM and ABAQUS units,nodes and operation time

通過本節(jié)分析可以發(fā)現(xiàn)，半解析有限元不僅在計(jì)算精度上基本達(dá)到 ABAQUS水準(zhǔn)，同時(shí)利用傅里葉變換和正交函數(shù)的優(yōu)勢，可以極大地減少計(jì)算成本。隨著分析對象結(jié)構(gòu)尺寸的增加和材料復(fù)雜程度的加深，這一優(yōu)勢將變得更加突出。

3.5 不同荷位下路面板彎拉應(yīng)力對比分析

在水泥混凝土路面路面荷載應(yīng)力分析中，通常選擇使面層板產(chǎn)生最大應(yīng)力的荷載位置作為應(yīng)力計(jì)算的臨界荷位，臨界荷位一般發(fā)生在水泥路面板板底?，F(xiàn)通過半解析有限元程序CP-SAFEM，采用3.4節(jié)水泥混凝土路面結(jié)構(gòu)模型，縱縫邊界不考慮拉桿，分別考慮板中、縱縫中部、角隅三種荷位，作路面板層底橫向和縱向拉應(yīng)力的對比分析，據(jù)此尋找臨界荷位。

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

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

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

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

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

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

4 結(jié)論

本文介紹了基于半解析有限元方法的水泥混凝土路面力學(xué)響應(yīng)計(jì)算理論，并根據(jù)相關(guān)理論和MATLAB開發(fā)了半解析有限元程序 CP-SAFEM，結(jié)合具體工況驗(yàn)證了 CP-SAFEM 具有良好的準(zhǔn)確性與計(jì)算效率。

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

(2)對于板中荷位下的水泥路面，CP-SAFEM和ABAQUS各項(xiàng)力學(xué)指標(biāo)誤差均在5%以內(nèi)，進(jìn)一步驗(yàn)證了CP-SAFEM的準(zhǔn)確性。

(3)針對同一路面結(jié)構(gòu)，采用相同的橫截面網(wǎng)格劃分形式，對比了CP-SAFEM和ABAQUS的模型復(fù)雜程度和計(jì)算時(shí)間，CP-SAFEM模型的單元數(shù)目僅為ABAQUS的1/50，運(yùn)行時(shí)長僅為1/6，證明了CP-SAFEM良好的計(jì)算效率。

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