基于顯式有限元的無砟軌道系力學行為分析與程序開發

甘 愿，黃鈺潤

(西南交通大學土木工程學院，四川成都 610031)

無砟軌道在我國應用廣泛，但隨著我國高速鐵路的快速發展，行車速度越來越快，各部件的變形失效明顯加快，輪軌系統的動力作用也越來越明顯，在列車長期循環荷載和混凝土橋梁收縮徐變等作用下產生的橋墩變形和梁體結構變形、砂漿損傷、板端上拱變形等典型病害引起無砟軌道系統出現局部脫空，嚴重影響了無砟軌道系統的動力特性。

由于脫空的存在，系統的力學行為因可能的碰撞而演變為高度非線性，處理此類問題時，顯式有限元相對于隱式有限元具有先天的優勢，其是否收斂只與時間步長有關，具有良好的穩定性，且可直接解耦，不形成大型矩陣，占用內存空間較小。目前常用的商業顯式有限元軟件有LS-DYNA、ABAQUS/Explicit等，但由于代碼未開源，計算過程不透明，也不利于進行后續的參數化設計、結構優化、集成計算等工作，故基于顯式有限元分析方法，自主開發了應用于有脫空無砟軌道應力分析的通用計算軟件。

1 有限元模型建立

1.1 無砟軌道有限元離散

目前無砟軌道的理論計算方法中運用最多的是疊合梁模型、梁板模型和梁體模型，因為本文研究重點是軌道系統在有脫空狀態下應力狀態及程序應用可行性，故采用梁體模型。用減縮六面體實體單元離散軌道板，砂漿填充層，混凝土底座；鋼軌采用彈性點支承梁模型，用大轉動空間梁單元離散；考慮扣件的尺寸效應，用彈簧單元連接鋼軌節點與軌道板扣件尺寸范圍內節點，在不考慮地基不均勻沉降時，同樣也用彈簧單元模擬地基的支撐作用。

1.2 接觸問題顯式有限元列式

顯式有限元處理接觸問題的有限元理論與基本求解方法可以參考文獻[1-5]，在這里僅給出不考慮摩擦接觸的基本公式。

假定存在2個物體A與B，或者一個物體的2個部分A和B，在t1時刻占據空間VA和VB，其邊界分別為ΓA與ΓB，當至t2時刻時，系統的空間位置或者構型發生改變，2個物體發生了接觸。其基本控制方程與不含接觸系統的控制方程是一致的，但在接觸界面上需增加額外的動力學與運動學條件，運動學約束為不可侵入條件，即對于ΓA中任意一點xA與ΓB中距離最近一點xB滿足:

gn=(xA-xB)·n≥0

(1)

式中:n為xB處ΓB的外法線。同時，法向作用力只能為壓力:

Pn≤0

(2)

根據虛功原理，系統的運動方程弱形式可以表述為:

(3)

pN=kNgN(gN<0)

(4)

對式(3)進行離散，可以得到系統的半離散:

Mü=f

(5)

式中:f=fext-fint+fc。fint為系統內力，fext為系統外力，fc為接觸力，M為質量矩陣，ü為加速度向量。

上式以中心差分法直接積分，由于本程序中采用集中質量矩陣，因此可以將式(5)解耦。如已求得tn時刻的解，時刻tn+1的位移為:

(6)

2 有限元程序設計與開發

有限元程序的設計內容大致可以分為前處理，分析計算，后處理3個部分。其中分析計算部分由于有限元分析方法其需要描述真實世界，需要對實例進行多個層級的抽象與離散，天然的與面向對象程序語言的形式邏輯吻合；同時從節點到單元和材料、對于各類荷載與邊界條件的處理，再到集總計算，這之間涉及到大量的數據交換與處理問題，使得程序不可避免的變得龐大易錯。考慮到上述問題，本文選擇面向對象的c++語言在Visual Studio 2019集成開發環境下進行設計開發。將各層級的問題抽象為類進行組織，從而通過對不同類的實例化來實現有限元分析，提高程序的可重用性、易維護性和可拓展性。程序框架如圖 1所示。

圖1 程序框架

2.1 前處理模塊

用戶在通過前處理模塊完成與有限元分析系統的交互，控制一系列參數完成幾何模型的建立，同時確定材料特性、單元劃分、荷載和邊界條件、計算終止時間，脫空位置等數據，再進行整理并傳遞給分析計算模塊。

2.2 分析計算模塊

此模塊為有限元計算的核心，接收前處理模塊整理后的數據，實例化節點類，再以節點類為基礎，結合材料參數實例化單元類。其中節點類儲存了節點的節點號、位置、慣性、速度、加速度、節點力等信息，單元類儲存了此單元成員節點的索引、單元號、單元集中質量矩陣、單元體積(長度)、自由度等信息。再將各實例化后的節點與單元輸入系統集總類，形成整體集中質量矩陣，遍歷節點，將集總得到的慣性賦予節點相應的自由度，同時進行加載與邊界條件的處理。最后再引入運算類進行迭代計算，其中包含節點內力計算、接觸搜索、接觸力計算與沙漏力計算等過程。分析計算模塊程序執行流程如圖2所示。

圖2 分析計算模塊程序執行流程

2.3 后處理模塊

TECPLOT是一款功能強大的數據分析和可視化處理的軟件，它提供了豐富的繪圖格式，能后根據有限元分析輸出的節點信息生成網格，支持云圖繪制等有限元后處理過程中不可缺少的功能。

故本文使用TECPLOT軟件進行可視化后處理，后處理模塊整理計算結果，根據用戶需求提取不同時間步結果，輸出為結構化數據導入TECPLOT軟件。

3 程序驗證

3.1 算例建立

本文依據CRTS Ⅰ 型板式無砟軌道結構參數建立有限元模型作為算例，與通用顯式有限元軟件ANSYS/LS-DYNA進行對比驗證。

在ANSYS/LS-DYNA軟件中，用Beam161單元劃分軌道，對應本文程序中的大轉動空間梁單元；用Solid164單元劃分軌道板、砂漿與混凝土底座板，對應于本文程序中的減縮六面體實體單元；用Spring-damping165單元模擬地基與扣件，對應于本文程序中的彈簧單元。

軌道采用CHN60鋼軌；扣件間距為0.625 m，剛度取50 MN/m；軌道板長度為5 m、寬度為2.4 m、厚度0.19 m；軌道板選用C50混凝土，彈性模量取36 000 MPa；砂漿充填層的長寬與軌道板相同，厚度為0.03 m，彈性模量取300 MPa；混凝土底座板寬3 m，厚0.3 m，選用C40混凝土，彈性模量取32 500 MPa；路基基礎面剛度取190 MPa/m。

為了避免邊界條件對結果的影響，建立3塊軌道板長度的模型，取中間塊作為分析和研究的對象，模型如圖3所示。同時為了更好的模擬出CA砂漿離縫脫空狀態下軌道的力學行為，選擇最不利位置進行加載。據參考文獻[12]可知，當荷載加載于扣件附近的鋼軌之上時，結構的受力最明顯，故本文在模型縱向方向8 500 mm位置設置長1 500 mm，寬2 400 mm，高0.1 mm的離縫，再選擇在離縫范圍內中間位置的扣件上方作為最不利位置，按照單軸雙輪荷載施加集中力P，示意見圖4，其中P取150 kN，以隨時間線性增大的方式加載于處于最不利位置的軌道上，以達到擬靜態的效果。

圖3 軌道模型

圖4 荷載示意

3.2 結果分析

提取ANSYS/LS-DYNA與本文程序結果進行對比研究，發現：

(1) 考慮到有限元計算過程包含接觸分析，故提取由于離縫脫空而與砂漿發生接觸行為的軌道板部分計算結果進行研究分析。通過對比軌道板部分最大垂向擾度關于時間的發展曲線,如圖5所示，ANSYS/LS-DYNA與本文程序結果符合良好，發展趨勢一致，最大相對誤差為2.6%，發生在0.19 s時，參考圖6軌道板部分最大垂向加速度關于時間的發展曲線，可以看到本文程序與ANSYS/LS-DYNA計算結果的相對誤差同樣也是在0.19 s附近波動較大，而這正是發生接觸的時刻，故推測誤差可能是因為接觸算法部分的區別造成的。

(2) 提取1.0 s時刻模型整體的垂向擾度計算結果繪制云圖(圖7為ANSYS/LS-DYNA結果,圖8為本文程序結果)，可以看到垂向擾度分布規律一致，最大垂向擾度都出現在同一位置，位于加載位置扣件的下方。

圖5 軌道板垂向位移

圖6 軌道板垂向加速度

圖7 ANSYS/LS-DYNA垂向位移云圖

圖8 程序垂向位移云圖

4 結束語

本文基于顯式有限元的基本思想，針對有脫空軌道板力學分析時難以處理的高度非線性接觸問題，采用面向對象的c++語言設計開發了一套完整的有限元程序，并通過與通用商業有限元軟件ANSYS/LS-DYNA的對比驗證了程序的正確性與、實用性與高效性，為有脫空軌道板力學分析計算提供了參考，為未來更進一步的顯式有限元軟件開發工作打下基礎。