基于GPU的車輛-軌道-地基土耦合系統3D隨機振動并行計算方法

朱志輝 夏禹濤 王力東 劉禹兵

摘要：針對軌道不平順隨機特征導致車輛-軌道-地基土耦合系統隨機分析計算效率低的問題，采用虛擬激勵法降低大樣本分析的計算量；針對耦合系統等效剛度矩陣的稀疏特性，采用行壓縮（Compressed Sparse Row，CSR）格式存儲大型稀疏矩陣，采用預處理共軛梯度法（Preconditioned Conjugate Gradient，PCG）求解對稱正定的等效靜力平衡方程，最后通過MAT-LAB-CUDA（Compute Unified Device Architecture）混合平臺開發基于GPU的并行計算程序.數值算例表明：基于MATLAB-CUDA混合平臺求解等效靜力平衡方程的效率是串行多點同步算法的86.13倍，大大縮短了隨機振動分析的總計算時間，且內存占用小、易于在個人計算機上實施；采用PCG法求解車輛-軌道-地基土耦合系統形成的大型稀疏線性方程組時，建議以加速度指標作為迭代收斂精度的控制指標；可通過選取適當的迭代收斂精度，以達到計算精度和計算效率的平衡.

關鍵詞：隨機振動；GPU并行計算；3D有限元法；虛擬激勵法；車輛-軌道-地基土耦合模型

中圖分類號：U211.3；TP399文獻標志碼：A

基金項目：國家自然科學基金資助項目（51678576，52078498），National Natural Science Foundation of China（51678576，52078498）；國家重點研發計劃項目（2017YFB1201204），National Key R&D Program of China（2017YFB1201204）

A Parallel Computing Method for Three-dimensional Random Vibration of Train-track-soil Dynamic Interaction Based on GPU

ZHU Zhihui1，2，XIA Yutao1，WANG Lidong3，LIU Yubing1

（1. School of Civil Engineering，Central South University，Changsha 410075，China；2. National Engineering Laboratory for High Speed Railway Construction，Central South University，Changsha 410075，China；3. School of Civil Engineering，Changsha University of Science & Technology，Changsha 410114，China）

Abstract：Aiming at the issue of low efficiency in the random analysis and calculation of train-track-soil cou-pled system due to the random characteristics of track irregularity，the pseudo-excitation method（PEM）is used to re-duce the computing cost of large samples. Considering the sparsity of equivalent stiffness matrix of train-track-soil coupling system，the large-scale sparse matrix was stored in Compressed Sparse Row（CSR）format. The symmetric positive definite equation of equivalent static equilibrium is solved by the preconditioned conjugate gradient（PCG）method，and the parallel program is developed by the MATLAB-CUDA（Compute Unified Device Architecture）hy-brid platform. The numerical example shows that the efficiency of solving equivalent static equilibrium equation based on MATLAB-CUDA hybrid platform is 86.13 times that of the multi-point synchronization algorithm，which greatly reduces the total calculation time of the random vibration analysis，and the memory usage is small and it is easy to im-plement on a personal computer. When using PCG method to solve the large sparse linear equations of vehicle-trackfoundation soil coupling system，it is suggested to take the acceleration as the control index of convergence accuracy of iteration. By choosing the appropriate convergence accuracy of iteration，the balance of calculation accuracy and effi-ciency is achieved.

Key words：random vibration；GPU parallel computing；three dimensional finite element method；pseudo-exci-tation method；train-track-soil couple model

隨著列車速度、軸重、路網密度和行車密度的提高，鐵路運輸引起的振動對周圍環境的影響問題日漸突出.計算機性能的不斷提升使得有限元法[1-2]、有限元/無限元法[3]、邊界元法[4]等數值計算方法能夠廣泛用于研究車輛引起的周圍環境振動問題；根據所建立的模型維度可將數值模型分為2D模型[5]、2.5D模型[3]和3D模型[2].與3D模型相比，2D和2.5D模型假設軌道-地基土模型的幾何形狀和材料特性沿車輛運動方向保持不變或周期性變化，計算中只需建立線路的橫截面模型，具有計算效率高的優點.但車致環境振動實質上是3D空間上的工程問題，3D模型可以考慮地基土的不連續性及與附近結構的耦合[6]，因此仍具有不可替代的普適性.

同時，軌道不平順的隨機性導致輪軌相互作用力具有隨機性，移動車輛引起的任何特定地面位置的振動均是一個非平穩的隨機過程[7].采用蒙特卡羅法等常規方法分析隨機響應時由于樣本過多導致計算效率較低.針對隨機振動計算效率低的問題，林家浩等提出了高效的虛擬激勵法，該法已被廣泛地用于分析平穩或非平穩隨機激勵下線性系統的隨機響應[8].虛擬激勵法將平穩隨機振動分析轉化為諧響應分析，將非平穩隨機振動分析轉化為可用數值積分方法求解的確定性瞬態分析，與蒙特卡羅方法相比，其計算效率提高了1～2個數量級[9-10]. Wang等[2]采用虛擬激勵法分析車輛-軌道-地基土的隨機振動，基于填充元優化與多波前法提出了一種求解大型稀疏線性方程組的多點同步算法，進一步節省了計算時間.這些方法大部分都在傳統CPU串行平臺上開展并實施.由于車輛-軌道-地基土耦合系統3D有限元模型的自由度數量龐大，且使用虛擬激勵法進行隨機振動分析時仍需求解一系列大型稀疏線性方程組，因此計算效率低仍然是基于3D模型開展隨機振動分析面臨的主要問題.

近年來，并行計算等高性能計算技術的發展為大規模的工程計算問題提供了高效的解決方案.在隨機振動領域，聶旭濤等[11]結合有限元EBE并行策略和虛擬激勵法實現了結構隨機振動的并行計算，其中通過消息傳遞接口（MPI）實現多CPU節點間的通訊，并證明了該算法的準確性和高效性.張健[12]利用虛擬激勵法中每個頻點的隨機響應計算相互獨立的特點，采用包含8個CPU的工作站實現虛擬激勵法的并行計算，每個CPU平均分配各頻點的計算任務，結果表明計算加速比幾乎與CPU數目成線性關系.這些方法通過CPU集群實現程序的并行化加速，但是CPU集群平臺搭建與維護的成本較高.隨著GPU通用計算的發展，CUDA等統一計算架構和編程環境的出現為GPU并行計算提供了易用的編程接口.相對于傳統的CPU串行計算，GPU的浮點計算能力更強、內存帶寬更大，具有低功耗和高性價比的特點，因此越來越廣泛地應用于科學研究和工程計算中.陳曦等[13]通過使用CPU-GPU混合計算構架加速了預處理Krylov子空間迭代法，從而提升了巖土工程有限元分析的效率；蔡勇等[14]基于GPU并行計算平臺CUDA，使用多重網格法加快Jacobi迭代的收斂速度，在采用GTX 460顯卡的個人計算機上有效提高了大規模殼結構有限元分析的效率.

綜上所述，針對軌道不平順隨機特征導致車輛-軌道-地基土耦合系統隨機分析計算效率低的問題，采用虛擬激勵法降低大樣本分析的計算量；針對車輛-軌道-地基土耦合系統等效剛度矩陣的稀疏特性，采用行壓縮格式存儲大型稀疏矩陣，采用預處理共軛梯度法并行求解等效靜力平衡方程，最后通過MATLAB-CUDA混合平臺開發并行計算程序.通過數值算例驗證了該并行計算方法的高效性，并給出了使用該法進行車致場地振動分析時迭代收斂精度選取的建議.

1車輛-軌道-地基土耦合系統的3D隨機振動分析模型

1.1車輛-軌道-地基土耦合振動模型

車輛-無砟軌道-路堤-地基土系統可分為車輛子系統和軌道-地基土子系統，兩個子系統通過輪軌相互作用耦合為一個整體，如圖1所示.

2基于GPU的耦合時變系統隨機振動并行求解方法

使用虛擬激勵法分析隨機振動時，仍需求解一系列形如式（17）的等效靜力平衡方程，該步驟巨大的計算耗時是3D模型難以廣泛用于隨機振動分析的重要原因.根據并行系統的安姆達爾加速定律[17]（Amdahl’s law），對該步驟進行并行化加速能獲得更大的加速比.

2.1基于GPU的統一計算架構

CUDA（Compute Unified Device Architecture）平臺是NVIDIA公司于2006年提出的GPU通用計算編程接口，它大大簡化了通用計算向GPU上移植的難度，縮短了程序開發的周期. CUDA稱CPU為主機（Host），稱GPU為設備（Device），使用GPU設備進行并行計算主要包括3個步驟：①在主機端與設備端分別分配計算數據所需的內存，將數據從主機端復制到設備端；②由設備端執行并行計算程序；③將計算結果數據從設備端復制到主機端.

基于CUDA的并行計算是建立在CPU和GPU協同工作的異構并行計算.如圖3所示，其中CPU主要負責程序中邏輯流程的控制及執行部分串行代碼，GPU負責執行大規模數據的并行計算部分.每個GPU設備由若干個流處理器簇（SM）組成，每個SM配備多個流處理器（SP），在GPU上執行的并行程序稱為核函數（Kernel），其中規定了線程（Thread）的組織和內存的管理.線程是并行計算的基本構建塊，CUDA的編程模型將線程組合形成了線程束（Warp）、線程塊（Block）以及線程網格（Grid）. GPU多線程的計算構架適合于大型3D有限元模型分析中涉及的大規模矩陣運算等操作.

2.2等效靜力平衡方程的并行求解

耦合系統隨機振動分析中的等效靜力平衡方程為大型稀疏線性方程組.線性方程組的解法一般分為直接法和迭代法，直接法經過有限步計算可得到方程組的精確解，但矩陣分解產生的大量填充元導致直接法的存儲量較大.迭代法可利用矩陣稀疏的特性，因此所需的內存占用小且易于編程實現，被廣泛用于大型稀疏線性方程組的求解，但需考慮計算精度和收斂速度的問題，一般來說，迭代循環的次數和計算時間隨收斂精度取值的減小而增大.考慮到GPU內存的限制和并行編程的難度，本文選用迭代法求解等效靜力平衡方程.

3D有限元方法形成的大型稀疏矩陣中含有大量不參與矩陣運算的零元素，圖2所示矩陣的稀疏度達到0.107‰，采用壓縮存儲的方式存儲稀疏矩陣可大大減少迭代求解時的內存占用、避免額外的計算開銷.為了滿足3D有限元法生成的無規則的系數矩陣，本文采用對矩陣結構無要求且適合并行編程的CSR格式[18]存儲等效剛度矩陣K，該法對稀疏矩陣進行逐行壓縮處理，只存儲非零元素及其位置索引.對于一個包含nz個非零元素的n×n維稀疏矩陣，分別使用三個一維數組存儲矩陣信息，如圖4所示.其中，Val數組含nz個元素，用于存儲所有非零元素的值；RowPtr數組含（n+1）個元素，用于存儲矩陣每行首個非零元素在數組Val中的位置；ColInd數組含nz個元素，用于存儲每個非零元素所在列位置.程序中，Val數組的數據類型設為雙精度浮點型（double），RowPtr和ColInd數組的數據類型設為整型（int）.

針對圖2所示的系數矩陣對稱正定的特征，采用Krylov子空間迭代法中的共軛梯度法進行求解.為了改善迭代求解的收斂性能，采用矩陣預處理技術減少系數矩陣條件數、改善其病態特征從而加快收斂速度.常用的預處理技術包括Jacobi預處理、不完全Cholesky分解預處理（IC）、多項式預處理和對稱超松弛預處理（SSOR）[19]等.本文采用的Jacobi預處理是一種簡單易行的預處理方法，能在不增加額外內存需求、通信開銷小的前提下達到較好的收斂性能，從而提高計算效率.對于線性方程組Ax = b，取系數矩陣的對角項形成對角預處理矩陣，如（18）式所示：

M = diag（A）（18）

采用預處理技術的共軛梯度法稱為預處理共軛梯度法（PreconditionedConjugateGradientMethod，PCG）.

PCG法將等效靜力平衡方程的求解轉化為矩陣和向量間的運算，其中主要包括稀疏矩陣向量乘、向量更新和向量內積運算，它們都適合于利用GPU進行的細粒度并行計算.基于CUDA的并行PCG法計算流程如表1所示，其中Jacobi預處理通過自編核函數實現，其余的矩陣向量運算通過調用CUBLAS和CUSPARSE函數庫實現.

2.3耦合時變系統隨機振動并行計算流程

通過MATLAB-CUDA混合平臺開發了計算程序，其中CPU端的程序在擅長矩陣數值運算的MATLAB平臺上編制，通過調用CUDA程序在GPU端實現并行PCG法求解等效靜力平衡方程.其中，采用cudaMemcpy函數實現CPU和GPU間的數據傳輸，通過NVIDIA Visual Profiler軟件進行CUDA程序的耗時分析.結果表明：對于本文研究對象，每一時間步內CPU和GPU間的數據傳輸時間小于0.5 s，占該步總計算時間的1%以內，因此數據傳輸對計算效率的影響較小，該并行計算流程適合本文耦合時變系統的隨機動力分析.

MATLAB中調用CUDA并行程序的步驟如下：利用nvcc編譯器將CUDA代碼（.cu）編譯成目標文件（.obj），使用mex命令編譯C++代碼（.cpp）并將其鏈接到CUDA目標文件（.obj），最終得到可執行的動態鏈接庫文件（.mexw64）.其中，CUDA代碼（.cu）完成了主機端和設備端間的數據傳輸和PCG法的并行求解流程. C++代碼（.cpp）中定義了被MATLAB調用的外部子程序的入口地址、MATLAB系統和子程序間傳遞的參數.

3數值算例

3.1有限元模型和計算參數

以京滬高鐵某路基段試驗場地為例，建立了如圖6所示的有限元模型，模型參數通過現場實測獲得[20].為控制有限元模型的尺寸，除地表外的其他五個邊界面采用3D黏彈性人工邊界[21]模擬半無限域場地.地基土參數、有限元模型橫截面幾何參數、有限元模型單元尺寸及車輛參數，參考文獻[2]中數據.整個3D有限元模型由210 686個單元、193 135個節點和542 002個自由度組成，采用CSR格式存儲后的等效剛度矩陣所占內存僅為296.7 MB.

本文算例僅考慮單線行車，車輛選取8車編組的CRH3型高速列車組，車速設為316 km/h.選取德國低干擾高低不平順軌道譜作為隨機激勵，如圖7所示取其空間頻域范圍為0.0125×2π～1×2π，并在該范圍內等間距離散為100個頻點（即m = 100）.軌道-地基土系統有限元模型采用瑞利阻尼，阻尼比取0.05.時間積分步長為0.002 s，當考慮車輛的運行時間為4 s時，總時間步為2 000.

3.2計算精度分析

為了驗證本文方法的計算精度，考慮到PCG法中迭代收斂精度ε對車輛-軌道-地基土耦合系統3D隨機振動計算精度的影響，分別取ε為10-7和10-6進行計算，并與文獻[2]中串行的ANSYS求解所得結果進行對比.選取了如圖8所示的A、B、C和D四個采樣點，它們與軌道中心線的距離分別是0 m、7.5 m、20 m和40 m.本文所采用的計算平臺參數如表2所示.圖9和圖10分別給出了四個采樣點豎向速度和豎向加速度標準差的ANSYS求解和不同精度下本文方法計算的結果.

結果表明：

1）當收斂精度ε取10-7時，本文方法所得四個采樣點的計算結果與ANSYS所得結果均吻合良好，具有較高的精度.

2）當收斂精度ε取10-6時，A點的計算結果吻合良好；3 s后離軌道中心線較遠的B、C和D點的豎向速度標準差輕微偏離，豎向加速度標準差嚴重偏離.經過分析，造成該現象的主要原因包括：①每一步的計算精度下降導致PCG算法的截斷誤差隨著時間步的積累效應逐漸顯著.②3 s之后振動響應的數值減小，同時距離軌道中心線越遠的振動響應數值越小，此時計算機字長限制導致的舍入誤差愈發明顯.

3）加速度指標對收斂精度ε的敏感性更高，采用PCG法求解車輛-軌道-地基土耦合系統形成的大型稀疏線性方程組時，建議以加速度指標作為迭代收斂精度的控制指標；可通過選取適當的迭代收斂精度，以達到計算精度和計算效率的平衡，如當關注的振動范圍較小時，取較大的ε可減少迭代次數、縮短求解時間.

3.3計算效率分析

為驗證本文方法的高效性，采用表3所示的3種不同方法求解車輛-軌道-地基土耦合系統3D隨機動力分析中的等效靜力平衡方程.其中，Case 1與Case 2的計算方法相同，均采用文獻[2]中填充元優化的串行多點同步算法；Case 3采用本文提出的MATLAB-CUDA混合編程的并行PCG方法，迭代收斂精度ε取10-7. Case 1和Case 3均在如表2所示的計算平臺上實施；Case 2的計算平臺硬件環境[2]為Intel Xeon E5-2620 CPU，內存為64 G.

表4給出了不同方法的計算效率對比結果，其中計算加速比以Case 1一個時間步的計算耗時為基準計算得到.考慮到該算例完整的隨機振動計算包含2 000個時間步，而Case 1的總計算時間預計可達180 d左右，故未開展Case 1的全過程計算.

由表4數據可得以下結論：

1）對比Case 1和Case 3可知，相同計算平臺下計算一個時間步，本文方法的計算效率是多點同步算法的86.13倍，大大減少了計算耗時，并將總計算時間減少到可接受的3.5 d.

2）對比Case 1和Case 2可知，多點同步算法在不同的計算平臺實施時，Case 2的計算效率是Case 1的13.74倍.主要原因是多點同步法[2]屬于線性方程組解法中的直接法，對內存的需求大，在內存較小的計算機上使用直接法進行大規模3D有限元方程求解時，會導致計算效率大幅降低.

3）對比Case 1、Case 2和Case 3可知，本文方法采用的稀疏矩陣CSR壓縮存儲格式和預處理共軛梯度法有效降低了算法的內存需求、提高了計算效率，易于在一般個人計算機上實施計算分析.

4結論

本文結合虛擬激勵法與并行的預處理共軛梯度法兩種高效算法，提出了一種基于GPU求解車輛-軌道-地基土耦合系統3D隨機振動的并行計算方法.通過不同工況的對比分析，得出以下結論與建議：

1）本文方法將GPU并行計算技術應用于車輛-軌道-地基土耦合系統3D隨機振動分析，在保證精度的前提下大大提升了計算效率，且對計算機內存的需求小，在一般個人計算機的硬件配置下就能達到較好的計算效率提升，有利于3D模型在隨機振動分析中的推廣應用.

2）在工程設計中，本文方法有助于更快速地對鐵路運輸引起的環境振動問題進行評估.例如，使用本文所得結果同時結合3σ法則可高效地確定車輛-軌道-地基土耦合系統隨機響應的上下限值.

3）本文方法結合了不同計算平臺的優勢，采用商業數學軟件MATLAB串行處理隨機振動計算涉及的大部分運算及計算流程中的各種邏輯控制，通過CUDA調用GPU實現系統等效靜力方程的高效并行求解.

4）采用PCG法求解車輛-軌道-地基土耦合系統形成的大型稀疏線性方程組時，建議以加速度指標作為計算精度的控制指標.同時可通過選取適當的迭代收斂精度ε，以達到計算精度和計算效率的平衡.

參考文獻

[1]KOUROUSSIS G，VAN PARYS L，CONTI C，et al. Using three-di-mensional finite element analysis in time domain to model railwayinduced ground vibrations[J]. Advances in Engineering Software，2014，70：63—76.

[2]WANG L D，ZHU Z H，BAI Y，et al. A fast random method for three-dimensional analysis of train-track-soil dynamic interaction[J]. Soil Dynamics and Earthquake Engineering，2018，115：252—262.

[3]YANG Y B，HUNG H H. A 2.5D finite/infinite element approach for modelling visco-elastic bodies subjected to moving loads[J]. Inter-national Journal for Numerical Methods in Engineering，2001，51（11）：1317—1336.

[4]SHENG X，JONES C J C，THOMPSON D J. Prediction of ground vi-bration from trains using the wavenumber finite and boundary ele-ment methods[J].Journal of Sound and Vibration，2006，293（3/4/ 5）：575—586.

[5]GODINHO L，AMADO-MENDES P，PEREIRA A，et al. A coupled MFS-FEM model for 2-D dynamic soil-structure interaction in the frequency domain[J].Computers & Structures，2013，129：74—85.

[6]GALVIN P，ROMERO A，DOMINGUEZ J. Fully three-dimensional analysis of high-speed train-track-soil-structure dynamic interac-tion[J].Journal of Sound and Vibration，2010，329（24）：5147—5163.

[7]SI L T，ZHAO Y，ZHANG Y H，et al. Random vibration of an elastic half-space subjected to a moving stochastic load[J].Computers & Structures，2016，168：92—105.

[8]林家浩，張亞輝，趙巖.虛擬激勵法在國內外工程界的應用回顧與展望[J].應用數學和力學，2017，38（1）：1—32. LIN J H，ZHANG Y H，ZHAO Y. The pseudo-excitation method and its industrial applications in China and abroad[J].Applied Mathe-matics and Mechanics，2017，38（1）：1—32.（In Chinese）

[9]LU F，LIN J H，KENNEDY D，et al. An algorithm to study non-sta-tionary random vibrations of vehicle-bridge systems[J]. Computers& Structures，2009，87（3/4）：177—185.

[10]ZENG Z P，LIU F S，LOU P，et al. Formulation of three-dimensional equations of motion for train-slab track-bridge interaction system and its application to random vibration analysis[J]. Applied Math-ematical Modelling，2016，40（11/12）：5891—5929.

[11]聶旭濤，范大鵬.基于EBE策略實現結構動力響應的并行計算[J].振動與沖擊，2007，26（10）：51—55. NIE X T，FAN D P. Implementation of structural dynamic response parallel computation based on ebe policy[J]. Journal of Vibration and Shock，2007，26（10）：51—55.（In Chinese）

[12]張健.車輛-軌道耦合系統動力學的數值方法研究[D].大連：大連理工大學，2014. ZHANG J.Studies on the numerical methods for dynamics of coupled vehicle-track system[D].Dalian：Dalian University of Technology，2014.（In Chinese）

[13]陳曦，王冬勇，任俊，等. CPU-GPU混合計算構架在巖土工程有限元分析中的應用[J].土木工程學報，2016，49（6）：105—112. CHEN X，WANG D Y，REN J，et al. Application of hybrid CPUGPU computing platform in large-scale geotechnical finite element analysis[J].China Civil Engineering Journal，2016，49（6）：105—112.（In Chinese）

[14]蔡勇，李光耀，王琥.基于多重網格法和GPU并行計算的大規模殼結構快速計算方法[J].工程力學，2014，31（5）：20—26. CAI Y，LI G Y，WANG H. A fast calculation method for large-scale shell structure based on multigrid method and GPU parallel comput-ing[J]. Engineering Mechanics，2014，31（5）：20—26.（In Chi-nese）

[15]ZHU Z H，GONG W，WANG L D，et al. A hybrid solution for study-ing vibrations of coupled train-track-bridge system[J]. Advances in Structural Engineering，2017，20（11）：1699—1711.

[16]朱志輝，王力東，龔威，等.多種垂向輪軌關系的對比及改進的車-線-橋系統迭代模型的建立[J].中南大學學報（自然科學版），2017，48（6）：1585—1593. ZHU Z H，WANG L D，GONG W，et al. Comparative analysis of several types of vertical wheel/rail relationship and construction of an improved iteration model for train-track-bridge system[J]. Jour-nal of Central South University（Science and Technology），2017，48（6）：1585—1593.（In Chinese）

[17]WAIVIO N. Parallel test description and analysis of parallel test system speedup through Amdahl’s law[C]//2007 IEEE Autotestcon. Baltimore，MD，USA：IEEE，2007：735—740.

[18]李佳佳，張秀霞，譚光明，等.選擇稀疏矩陣乘法最優存儲格式的研究[J].計算機研究與發展，2014，51（4）：882—894. LI J J，ZHANG X X，TAN G M，et al. Study of choosing the optimal storage format of sparse matrix vector multiplication[J]. Journal of Computer Research and Development，2014，51（4）：882—894.（In Chinese）

[19]GEORGESCU S，CHOW P，OKUDA H. GPU acceleration for FEMbased structural analysis[J]. Archives of Computational Methods in Engineering，2013，20（2）：111—121.

[20]BIAN X C，JIANG H G，CHANG C，et al. Track and ground vibra-tions generated by high-speed train running on ballastless railway with excitation of vertical track irregularities[J]. Soil Dynamics and Earthquake Engineering，2015，76：29—43.

[21]LIU J B，DU Y X，DU X L，et al. 3D viscous-spring artificial bound-ary in time domain[J]. Earthquake Engineering and Engineering Vibration，2006，5（1）：93—102.