基于GPU的車輛-軌道-地基土耦合系統(tǒng)3D隨機(jī)振動(dòng)并行計(jì)算方法

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

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

關(guān)鍵詞：隨機(jī)振動(dòng)；GPU并行計(jì)算；3D有限元法；虛擬激勵(lì)法；車輛-軌道-地基土耦合模型

中圖分類號(hào)：U211.3；TP399文獻(xiàn)標(biāo)志碼：A

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目（51678576，52078498），National Natural Science Foundation of China（51678576，52078498）；國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（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

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

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

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

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

1車輛-軌道-地基土耦合系統(tǒng)的3D隨機(jī)振動(dòng)分析模型

1.1車輛-軌道-地基土耦合振動(dòng)模型

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

2基于GPU的耦合時(shí)變系統(tǒng)隨機(jī)振動(dòng)并行求解方法

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

2.1基于GPU的統(tǒng)一計(jì)算架構(gòu)

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

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

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

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

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

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

M = diag（A）（18）

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

PCG法將等效靜力平衡方程的求解轉(zhuǎn)化為矩陣和向量間的運(yùn)算，其中主要包括稀疏矩陣向量乘、向量更新和向量?jī)?nèi)積運(yùn)算，它們都適合于利用GPU進(jìn)行的細(xì)粒度并行計(jì)算.基于CUDA的并行PCG法計(jì)算流程如表1所示，其中Jacobi預(yù)處理通過自編核函數(shù)實(shí)現(xiàn)，其余的矩陣向量運(yùn)算通過調(diào)用CUBLAS和CUSPARSE函數(shù)庫(kù)實(shí)現(xiàn).

2.3耦合時(shí)變系統(tǒng)隨機(jī)振動(dòng)并行計(jì)算流程

通過MATLAB-CUDA混合平臺(tái)開發(fā)了計(jì)算程序，其中CPU端的程序在擅長(zhǎng)矩陣數(shù)值運(yùn)算的MATLAB平臺(tái)上編制，通過調(diào)用CUDA程序在GPU端實(shí)現(xiàn)并行PCG法求解等效靜力平衡方程.其中，采用cudaMemcpy函數(shù)實(shí)現(xiàn)CPU和GPU間的數(shù)據(jù)傳輸，通過NVIDIA Visual Profiler軟件進(jìn)行CUDA程序的耗時(shí)分析.結(jié)果表明：對(duì)于本文研究對(duì)象，每一時(shí)間步內(nèi)CPU和GPU間的數(shù)據(jù)傳輸時(shí)間小于0.5 s，占該步總計(jì)算時(shí)間的1%以內(nèi)，因此數(shù)據(jù)傳輸對(duì)計(jì)算效率的影響較小，該并行計(jì)算流程適合本文耦合時(shí)變系統(tǒng)的隨機(jī)動(dòng)力分析.

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

3數(shù)值算例

3.1有限元模型和計(jì)算參數(shù)

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

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

3.2計(jì)算精度分析

為了驗(yàn)證本文方法的計(jì)算精度，考慮到PCG法中迭代收斂精度ε對(duì)車輛-軌道-地基土耦合系統(tǒng)3D隨機(jī)振動(dòng)計(jì)算精度的影響，分別取ε為10-7和10-6進(jìn)行計(jì)算，并與文獻(xiàn)[2]中串行的ANSYS求解所得結(jié)果進(jìn)行對(duì)比.選取了如圖8所示的A、B、C和D四個(gè)采樣點(diǎn)，它們與軌道中心線的距離分別是0 m、7.5 m、20 m和40 m.本文所采用的計(jì)算平臺(tái)參數(shù)如表2所示.圖9和圖10分別給出了四個(gè)采樣點(diǎn)豎向速度和豎向加速度標(biāo)準(zhǔn)差的ANSYS求解和不同精度下本文方法計(jì)算的結(jié)果.

結(jié)果表明：

1）當(dāng)收斂精度ε取10-7時(shí)，本文方法所得四個(gè)采樣點(diǎn)的計(jì)算結(jié)果與ANSYS所得結(jié)果均吻合良好，具有較高的精度.

2）當(dāng)收斂精度ε取10-6時(shí)，A點(diǎn)的計(jì)算結(jié)果吻合良好；3 s后離軌道中心線較遠(yuǎn)的B、C和D點(diǎn)的豎向速度標(biāo)準(zhǔn)差輕微偏離，豎向加速度標(biāo)準(zhǔn)差嚴(yán)重偏離.經(jīng)過分析，造成該現(xiàn)象的主要原因包括：①每一步的計(jì)算精度下降導(dǎo)致PCG算法的截?cái)嗾`差隨著時(shí)間步的積累效應(yīng)逐漸顯著.②3 s之后振動(dòng)響應(yīng)的數(shù)值減小，同時(shí)距離軌道中心線越遠(yuǎn)的振動(dòng)響應(yīng)數(shù)值越小，此時(shí)計(jì)算機(jī)字長(zhǎng)限制導(dǎo)致的舍入誤差愈發(fā)明顯.

3）加速度指標(biāo)對(duì)收斂精度ε的敏感性更高，采用PCG法求解車輛-軌道-地基土耦合系統(tǒng)形成的大型稀疏線性方程組時(shí)，建議以加速度指標(biāo)作為迭代收斂精度的控制指標(biāo)；可通過選取適當(dāng)?shù)牡諗烤?，以達(dá)到計(jì)算精度和計(jì)算效率的平衡，如當(dāng)關(guān)注的振動(dòng)范圍較小時(shí)，取較大的ε可減少迭代次數(shù)、縮短求解時(shí)間.

3.3計(jì)算效率分析

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

表4給出了不同方法的計(jì)算效率對(duì)比結(jié)果，其中計(jì)算加速比以Case 1一個(gè)時(shí)間步的計(jì)算耗時(shí)為基準(zhǔn)計(jì)算得到.考慮到該算例完整的隨機(jī)振動(dòng)計(jì)算包含2 000個(gè)時(shí)間步，而Case 1的總計(jì)算時(shí)間預(yù)計(jì)可達(dá)180 d左右，故未開展Case 1的全過程計(jì)算.

由表4數(shù)據(jù)可得以下結(jié)論：

1）對(duì)比Case 1和Case 3可知，相同計(jì)算平臺(tái)下計(jì)算一個(gè)時(shí)間步，本文方法的計(jì)算效率是多點(diǎn)同步算法的86.13倍，大大減少了計(jì)算耗時(shí)，并將總計(jì)算時(shí)間減少到可接受的3.5 d.

2）對(duì)比Case 1和Case 2可知，多點(diǎn)同步算法在不同的計(jì)算平臺(tái)實(shí)施時(shí)，Case 2的計(jì)算效率是Case 1的13.74倍.主要原因是多點(diǎn)同步法[2]屬于線性方程組解法中的直接法，對(duì)內(nèi)存的需求大，在內(nèi)存較小的計(jì)算機(jī)上使用直接法進(jìn)行大規(guī)模3D有限元方程求解時(shí)，會(huì)導(dǎo)致計(jì)算效率大幅降低.

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

4結(jié)論

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

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

2）在工程設(shè)計(jì)中，本文方法有助于更快速地對(duì)鐵路運(yùn)輸引起的環(huán)境振動(dòng)問題進(jìn)行評(píng)估.例如，使用本文所得結(jié)果同時(shí)結(jié)合3σ法則可高效地確定車輛-軌道-地基土耦合系統(tǒng)隨機(jī)響應(yīng)的上下限值.

3）本文方法結(jié)合了不同計(jì)算平臺(tái)的優(yōu)勢(shì)，采用商業(yè)數(shù)學(xué)軟件MATLAB串行處理隨機(jī)振動(dòng)計(jì)算涉及的大部分運(yùn)算及計(jì)算流程中的各種邏輯控制，通過CUDA調(diào)用GPU實(shí)現(xiàn)系統(tǒng)等效靜力方程的高效并行求解.

4）采用PCG法求解車輛-軌道-地基土耦合系統(tǒng)形成的大型稀疏線性方程組時(shí)，建議以加速度指標(biāo)作為計(jì)算精度的控制指標(biāo).同時(shí)可通過選取適當(dāng)?shù)牡諗烤圈?，以達(dá)到計(jì)算精度和計(jì)算效率的平衡.

參考文獻(xiàn)

[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]林家浩，張亞輝，趙巖.虛擬激勵(lì)法在國(guó)內(nèi)外工程界的應(yīng)用回顧與展望[J].應(yīng)用數(shù)學(xué)和力學(xué)，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策略實(shí)現(xiàn)結(jié)構(gòu)動(dòng)力響應(yīng)的并行計(jì)算[J].振動(dòng)與沖擊，2007，26（10）：51—55. NIE X T，F(xiàn)AN 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]張健.車輛-軌道耦合系統(tǒng)動(dòng)力學(xué)的數(shù)值方法研究[D].大連：大連理工大學(xué)，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ì)算構(gòu)架在巖土工程有限元分析中的應(yīng)用[J].土木工程學(xué)報(bào)，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]蔡勇，李光耀，王琥.基于多重網(wǎng)格法和GPU并行計(jì)算的大規(guī)模殼結(jié)構(gòu)快速計(jì)算方法[J].工程力學(xué)，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]朱志輝，王力東，龔?fù)?，?多種垂向輪軌關(guān)系的對(duì)比及改進(jìn)的車-線-橋系統(tǒng)迭代模型的建立[J].中南大學(xué)學(xué)報(bào)（自然科學(xué)版），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]李佳佳，張秀霞，譚光明，等.選擇稀疏矩陣乘法最優(yōu)存儲(chǔ)格式的研究[J].計(jì)算機(jī)研究與發(fā)展，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.