董曉輝 唐貞云 李振寶 杜修力

摘要 提出了基于Python和圖形處理器（GPU）數(shù)值求解的實(shí)時(shí)混合試驗(yàn)系統(tǒng)。將土?結(jié)相互作用系統(tǒng)作為試驗(yàn)?zāi)Ｐ?，使用Python?GPU代替CPU數(shù)值求解，對(duì)提出的實(shí)時(shí)混合試驗(yàn)系統(tǒng)進(jìn)行了仿真與試驗(yàn)驗(yàn)證。研究結(jié)果表明，使用Python?GPU對(duì)無條件穩(wěn)定算法求解，積分步長20 ms時(shí)GPU實(shí)時(shí)求解模型自由度超過24000，是同一臺(tái)計(jì)算機(jī)CPU求解規(guī)模的7倍左右，提升了實(shí)時(shí)混合試驗(yàn)的試驗(yàn)?zāi)芰Α?/p>

關(guān)鍵詞 實(shí)時(shí)混合試驗(yàn); 圖形處理器; Python; 數(shù)值積分算法; 求解效率

引 言

實(shí)時(shí)混合試驗(yàn)［1］是一種將數(shù)值仿真與物理試驗(yàn)相結(jié)合的試驗(yàn)方法，將試驗(yàn)結(jié)構(gòu)分為數(shù)值子結(jié)構(gòu)與物理子結(jié)構(gòu)，兩部分?jǐn)?shù)據(jù)實(shí)時(shí)交互傳輸，使得大型復(fù)雜結(jié)構(gòu)大比例試驗(yàn)成為可能，現(xiàn)已得到廣泛應(yīng)用。在實(shí)時(shí)混合試驗(yàn)中，數(shù)值子結(jié)構(gòu)動(dòng)力求解需在一個(gè)時(shí)間步長內(nèi)完成，數(shù)值子結(jié)構(gòu)模型自由度越高，求解耗時(shí)越大。過大的時(shí)間步長影響數(shù)值積分算法的精度及物理子結(jié)構(gòu)加載的準(zhǔn)確性，無法滿足實(shí)時(shí)混合試驗(yàn)的需要。在數(shù)值子結(jié)構(gòu)動(dòng)力分析過程中通常采用計(jì)算效率更高的顯式積分算法進(jìn)行求解，顯式積分算法（例如中心差分法）通常為條件穩(wěn)定算法，積分步長越小精度越高［2］。過小的步長限制了數(shù)值模型實(shí)時(shí)求解規(guī)模，難以實(shí)現(xiàn)大規(guī)模數(shù)值模型的實(shí)時(shí)求解［3］。為解決此問題眾多學(xué)者在改進(jìn)數(shù)值積分算法［3］及提升數(shù)值子結(jié)構(gòu)計(jì)算效率［4］等方面做出研究。

在改進(jìn)數(shù)值積分算法的研究中，Nakashima等［5］及Chang［6?7］分別提出擬動(dòng)力試驗(yàn)的無條件穩(wěn)定顯式積分算法，但這些方法在實(shí)時(shí)混合試驗(yàn)中為隱式算法。吳斌等［8］基于上述方法進(jìn)行研究，利用位移向前差分修正建立實(shí)時(shí)混合試驗(yàn)中無條件穩(wěn)定的顯式積分Chang算法。Chen等［9］基于離散控制原理和極值的概念提出CR法，在線性結(jié)構(gòu)和部分非線性系統(tǒng)中具有無條件穩(wěn)定性，且速度和位移均顯式。Tang等［10］基于運(yùn)動(dòng)方程離散化提出TL算法，在線性結(jié)構(gòu)系統(tǒng)中具有無條件穩(wěn)定性，且速度和位移均顯式表達(dá)。基于現(xiàn)有的動(dòng)力分析算法，在提升實(shí)時(shí)混合實(shí)驗(yàn)數(shù)值子結(jié)構(gòu)計(jì)算效率方面也有眾多學(xué)者研究。Nakashima等［11］提出Digital Signal Processor（DSP）方法，在一個(gè)時(shí)間步長內(nèi)完成數(shù)值模型仿真和動(dòng)力響應(yīng)信號(hào)處理，實(shí)現(xiàn)了10自由度模型330 ms時(shí)間步長或12自由度模型500 ms時(shí)間步長的實(shí)時(shí)混合試驗(yàn)。Chae等［12］開發(fā)Hybrid FEM程序，實(shí)現(xiàn)514自由度框架模型10 ms時(shí)間步長的實(shí)時(shí)混合試驗(yàn)。Zhu等［13］提出使用兩臺(tái)目標(biāo)計(jì)算機(jī)的實(shí)時(shí)混合試驗(yàn)D?RTHS系統(tǒng)提升求解規(guī)模，實(shí)現(xiàn)超過1240自由度時(shí)間步長為20 ms的實(shí)時(shí)混合試驗(yàn)。Lu等［14］在D?RTHS系統(tǒng)的基礎(chǔ)上結(jié)合計(jì)算機(jī)并行計(jì)算提出RTHS?W系統(tǒng)，使用條件穩(wěn)定的中心差分法實(shí)現(xiàn)時(shí)間步長為20 ms條件不超過7000自由度（Rayleigh阻尼β≠0）或40602自由度（Rayleigh阻尼β=0）模型的實(shí)時(shí)混合試驗(yàn)。綜上所述，通過改進(jìn)實(shí)時(shí)混合試驗(yàn)數(shù)值求解算法及求解的方法，可以提升數(shù)值子結(jié)構(gòu)模型實(shí)時(shí)求解能力，但目前的求解規(guī)模仍在7000以內(nèi)（Rayleigh阻尼β≠0）。

在實(shí)時(shí)混合試驗(yàn)中，數(shù)值子結(jié)構(gòu)求解通?；谟?jì)算機(jī)的中央處理器（CPU）運(yùn)算，然而在進(jìn)行大量數(shù)據(jù)運(yùn)算時(shí)，CPU的計(jì)算能力難以滿足模型實(shí)時(shí)求解的需要。圖形處理器（GPU）芯片與CPU的架構(gòu)不同，GPU中有數(shù)倍于CPU的計(jì)算單元，在進(jìn)行大規(guī)模數(shù)值計(jì)算時(shí)GPU具有明顯優(yōu)勢。在土木工程領(lǐng)域，GPU的高性能求解已得到應(yīng)用。劉曉強(qiáng)［15］應(yīng)用GPU提升大規(guī)模數(shù)據(jù)繪制和拓?fù)浞治鲂阅?，?shí)現(xiàn)相比CPU數(shù)十倍的加速效果，滿足土木工程領(lǐng)域防災(zāi)減災(zāi)數(shù)據(jù)可視化需求。解琳琳等［16］應(yīng)用GPU加速OpenSees有限元仿真，對(duì)城市地震災(zāi)害進(jìn)行模擬，實(shí)現(xiàn)相比CPU計(jì)算39倍的加速效果。在實(shí)時(shí)混合試驗(yàn)研究方面，董曉輝等［17］已應(yīng)用MATLAB軟件實(shí)現(xiàn)GPU求解數(shù)值子結(jié)構(gòu)，基于GPU求解實(shí)時(shí)混合試驗(yàn)架構(gòu)并驗(yàn)證其性能，證明了應(yīng)用GPU求解數(shù)值模型可以提升試驗(yàn)?zāi)Ｐ鸵?guī)模，降低時(shí)間步長。在MATLAB?GPU實(shí)時(shí)混合試驗(yàn)系統(tǒng)架構(gòu)中，使用GPU加速條件穩(wěn)定的顯式算法求解，并且子結(jié)構(gòu)之間通訊存在3 ms的延遲。為解決上述問題，本文提出使用Python編程建立基于GPU計(jì)算的實(shí)時(shí)混合試驗(yàn)系統(tǒng)，優(yōu)化基于GPU求解的實(shí)時(shí)混合試驗(yàn)系統(tǒng)，以土?結(jié)相互作用系統(tǒng)為模型進(jìn)行振動(dòng)臺(tái)實(shí)時(shí)混合試驗(yàn)，測試架構(gòu)的可行性及性能。

1 試驗(yàn)系統(tǒng)組成

在基于Python?GPU求解的實(shí)時(shí)混合試驗(yàn)系統(tǒng)中，使用GPU代替CPU作為數(shù)值模型求解的硬件。為了實(shí)現(xiàn)數(shù)值模型的高效仿真及數(shù)據(jù)實(shí)時(shí)交互，試驗(yàn)系統(tǒng)如圖1所示，分為三個(gè)部分：數(shù)值子結(jié)構(gòu)部分，信號(hào)處理部分和試驗(yàn)子結(jié)構(gòu)部分。

1.1 數(shù)值子結(jié)構(gòu)部分

數(shù)值子結(jié)構(gòu)部分負(fù)責(zé)數(shù)值模型實(shí)時(shí)動(dòng)力分析，在此部分使用Python作為編程語言編寫動(dòng)力求解算法，調(diào)用CuPy函數(shù)庫實(shí)現(xiàn)GPU運(yùn)算［18］。在模型仿真程序中添加UDP網(wǎng)絡(luò)通訊模塊［19］，實(shí)現(xiàn)基于GPU的模型動(dòng)力求解與試驗(yàn)子結(jié)構(gòu)實(shí)時(shí)通訊。

1.2 信號(hào)處理部分

信號(hào)處理部分負(fù)責(zé)固定數(shù)值子結(jié)構(gòu)求解的時(shí)間步長，信號(hào)插值以及振動(dòng)臺(tái)控制補(bǔ)償。在實(shí)時(shí)混合試驗(yàn)中，數(shù)值子結(jié)構(gòu)快速且準(zhǔn)確的求解是試驗(yàn)成功的關(guān)鍵，試驗(yàn)的時(shí)間步長需要保持固定不變。由于數(shù)值模型仿真在非實(shí)時(shí)系統(tǒng)Windows［20］上運(yùn)行，每一步數(shù)值模型求解所用的時(shí)間不固定，對(duì)此需要保證每一步模型時(shí)間步長固定且求解耗時(shí)均小于時(shí)間步長。控制時(shí)間步長可以通過LabVIEW等軟件的定時(shí)循環(huán)或在數(shù)值求解計(jì)算機(jī)外添加實(shí)時(shí)硬件。

本文在兩個(gè)子結(jié)構(gòu)之間添加Simulink Realtime xPC（下文簡稱xPC）實(shí)時(shí)硬件，xPC實(shí)時(shí)硬件有較好的實(shí)時(shí)性［21］，可以實(shí)現(xiàn)50 kHz的實(shí)時(shí)仿真，即可以實(shí)現(xiàn)最小時(shí)間步長為2×10-5的實(shí)時(shí)仿真，能滿足實(shí)時(shí)混合試驗(yàn)最小時(shí)間步長為1×10-3的要求。在xPC中添加UDP網(wǎng)絡(luò)通訊模塊［22］，連接數(shù)值子結(jié)構(gòu)與試驗(yàn)子結(jié)構(gòu)；添加信號(hào)插值模塊，對(duì)數(shù)值子結(jié)構(gòu)動(dòng)力響應(yīng)的離散信號(hào)進(jìn)行插值，使得輸入試驗(yàn)子結(jié)構(gòu)的信號(hào)指令更平滑。為降低試驗(yàn)子結(jié)構(gòu)加載裝置的加載誤差，在xPC中添加FSCS補(bǔ)償控制器［23］。完成插值及補(bǔ)償后，xPC將信號(hào)通過UDP傳輸至試驗(yàn)子結(jié)構(gòu)。

1.3 試驗(yàn)子結(jié)構(gòu)部分

試驗(yàn)子結(jié)構(gòu)部分與基于CPU的實(shí)時(shí)混合試驗(yàn)方法相同，通過物理加載的方法對(duì)試件試驗(yàn)。加載裝置的控制器接收到數(shù)值子結(jié)構(gòu)響應(yīng)后，對(duì)試驗(yàn)試件加載。在試驗(yàn)試件上布置加速度及位移傳感器，將時(shí)程響應(yīng)傳輸至加載系統(tǒng)控制器采集，再通過加載系統(tǒng)控制器將試件響應(yīng)反饋至數(shù)值子結(jié)構(gòu)。

2 基于Python的GPU求解性能

Python作為一種解釋型編程語言具有豐富的拓展庫，適用于各種系統(tǒng)平臺(tái)，隨著版本的不斷更新和語言新功能的添加，被廣泛應(yīng)用于科研領(lǐng)域。CuPy［18］是一個(gè)借助CUDA實(shí)現(xiàn)GPU加速的Numpy格式數(shù)據(jù)運(yùn)算庫，可以利用GPU的眾多CUDA核心提升數(shù)值計(jì)算效率。本文選擇Python作為編程語言，實(shí)現(xiàn)基于CuPy函數(shù)庫的數(shù)值子結(jié)構(gòu)模型GPU實(shí)時(shí)求解。

2.1 實(shí)施方案

使用Python進(jìn)行數(shù)值模型動(dòng)力求解需要建立數(shù)值子結(jié)構(gòu)動(dòng)力方程，模型的質(zhì)量、剛度和阻尼矩陣需要借助有限元軟件前處理得到。本文選擇使用ABAQUS有限元軟件建立數(shù)值子結(jié)構(gòu)模型，調(diào)整劃分網(wǎng)格完成模型前處理，添加模型參數(shù)矩陣輸出腳本，得到數(shù)值子結(jié)構(gòu)模型的質(zhì)量和剛度矩陣，模型阻尼采用Rayleigh阻尼。將模型的質(zhì)量、剛度、阻尼參數(shù)矩陣導(dǎo)入Python中建立結(jié)構(gòu)動(dòng)力方程。

使用CuPy庫進(jìn)行GPU加速計(jì)算前，將CuPy庫添加至數(shù)值子結(jié)構(gòu)求解的程序中，并將模型及算法所需的參數(shù)通過cupy.asarray函數(shù)轉(zhuǎn)換為CuPy數(shù)據(jù)格式［18］。為實(shí)現(xiàn)數(shù)值子結(jié)構(gòu)與試驗(yàn)子結(jié)構(gòu)之間數(shù)據(jù)交互，完成模型動(dòng)力分析后使用cupy.asnumpy函數(shù)將界面響應(yīng)轉(zhuǎn)換為Numpy數(shù)據(jù)格式，在每一步長動(dòng)力分析前后添加UDP通訊腳本，并在xPC軟件環(huán)境和試驗(yàn)子結(jié)構(gòu)部分中配置好UDP端口實(shí)現(xiàn)網(wǎng)絡(luò)通訊。使用Python?GPU進(jìn)行實(shí)時(shí)混合試驗(yàn)時(shí)，數(shù)值子結(jié)構(gòu)部分流程如圖2所示。

2.2 仿真參數(shù)及配置

為驗(yàn)證本文提出的Python?GPU求解數(shù)值模型方法的可行性與性能，采用如圖3所示的土?結(jié)相互作用模型作為算例，其中上部結(jié)構(gòu)為試驗(yàn)子結(jié)構(gòu)，下部土體為數(shù)值子結(jié)構(gòu)。數(shù)值子結(jié)構(gòu)動(dòng)力分析分別使用條件穩(wěn)定的中心差分法和無條件穩(wěn)定的Newmark?β算法。試驗(yàn)子結(jié)構(gòu)為單自由度鋁框架，質(zhì)量ms=14.5 t，剛度ks=10.03 kN/mm，阻尼cs=3.74 kN/（m?s?1），與數(shù)值子結(jié)構(gòu)頂部中心錨固連接。數(shù)值子結(jié)構(gòu)模型尺寸為30 m×30 m×15 m，密度為2.15×103 kg/m3，彈性模量為7.2×106 Pa，泊松比為0.35，阻尼比為0.05，模型四周及底部節(jié)點(diǎn)添加彈簧阻尼器模擬遠(yuǎn)場土邊界條件，節(jié)點(diǎn)法向彈簧剛度為20000 N/m，切向?yàn)?0000 N/m，法向阻尼器阻尼為1.437×106 N/（m?s?1），切向?yàn)?.45×106 N/（m?s?1）。數(shù)值模型在ABAQUS有限元軟件中建模，每個(gè)節(jié)點(diǎn)有三個(gè)自由度，通過調(diào)整網(wǎng)格劃分的間距調(diào)整模型的求解規(guī)模，網(wǎng)格越密，動(dòng)力求解計(jì)算量越大。仿真中模型荷載為Kobe地震波，加速度幅值調(diào)整為0.5g。仿真中使用GPU服務(wù)器作為求解計(jì)算機(jī)，軟件與硬件環(huán)境如表1所示。本文研究選擇的土?結(jié)相互作用模型僅為了驗(yàn)證Python?GPU求解的可行性與計(jì)算效率。使用Python?GPU的數(shù)值模型求解方法不限于本研究中的線彈性模型與模型邊界條件，適用于任何形式的實(shí)時(shí)混合試驗(yàn)。

2.3 仿真結(jié)果

數(shù)值子結(jié)構(gòu)求解的GPU計(jì)算機(jī)與物理子結(jié)構(gòu)仿真xPC之間通過UDP網(wǎng)絡(luò)連接，數(shù)據(jù)傳輸需要時(shí)間。為了測試UDP網(wǎng)絡(luò)通訊耗時(shí)在實(shí)時(shí)混合試驗(yàn)一個(gè)時(shí)間步長內(nèi)的占比，使用GPU計(jì)算機(jī)通過UDP網(wǎng)絡(luò)向xPC硬件發(fā)送數(shù)據(jù)10000次，計(jì)算單次UDP通訊耗時(shí)。使用Python中的time函數(shù)統(tǒng)計(jì)總耗時(shí)為0.078 s，平均單次通訊耗時(shí)為0.078/10000=7.87×10-3 ms，當(dāng)實(shí)時(shí)混合試驗(yàn)時(shí)間步長為1 ms時(shí)，UDP通訊耗時(shí)的占比為0.79%。使用UDP網(wǎng)絡(luò)傳輸數(shù)據(jù)的耗時(shí)在一個(gè)時(shí)間步長內(nèi)可以忽略不計(jì)。

為了測試GPU求解數(shù)值模型的精度與性能，分別選取自由度為2592，5400，10830以及24375的模型，時(shí)間步長為1 ms作為算例，中心差分法與Newmark?β法運(yùn)用CPU與GPU求解的試驗(yàn)子結(jié)構(gòu)加速度時(shí)程對(duì)比如圖4所示，GPU與CPU求解結(jié)果相同，說明GPU可以替代CPU作為實(shí)時(shí)混合試驗(yàn)的數(shù)值求解硬件。

在求解算法中使用Python中的time函數(shù)統(tǒng)計(jì)分別使用GPU與CPU，求解算法為中心差分法與Newmark?β法的計(jì)算耗時(shí)，結(jié)果如表2所示。表中SR（Speedup Ratio）為加速比，計(jì)算公式為SR=TCPU/TGPU，其中，TCPU為使用CPU求解每一步長平均耗時(shí)，TGPU為使用GPU求解每一步長平均耗時(shí)。

將表2中的數(shù)據(jù)繪圖如圖5所示，圖5中縱坐標(biāo)為單次求解模型所需的耗時(shí)，橫坐標(biāo)為模型自由度數(shù)。隨著模型自由度的提升，兩種算法求解所需的時(shí)間均有增加，使用CPU求解5400自由度模型耗時(shí)超過了實(shí)時(shí)混合試驗(yàn)的實(shí)時(shí)性要求（最大時(shí)間步長20 ms），而此時(shí)使用GPU求解僅需不到1 ms。模型規(guī)模越大，GPU相比CPU的加速越明顯，當(dāng)模型自由度為24375時(shí)，GPU相比CPU的加速比可達(dá)50余倍。無論是使用CPU還是GPU，條件穩(wěn)定的中心差分法求解所需的時(shí)間均小于無條件穩(wěn)定的Newmark?β法。原因是中心差分法的計(jì)算量相比Newmark?β法小，因此求解耗時(shí)短。

3 試驗(yàn)驗(yàn)證

本文使用Python作為編程語言，應(yīng)用CuPy庫進(jìn)行基于GPU求解的實(shí)時(shí)混合試驗(yàn)，設(shè)計(jì)了振動(dòng)臺(tái)實(shí)時(shí)混合試驗(yàn)系統(tǒng)，驗(yàn)證試驗(yàn)系統(tǒng)的可行性。數(shù)值子結(jié)構(gòu)部分GPU求解數(shù)值模型后將子結(jié)構(gòu)交界面動(dòng)力響應(yīng)發(fā)送至信號(hào)處理部分。在信號(hào)處理部分，配置好Simulink xPC軟件硬件環(huán)境，通過UDP網(wǎng)絡(luò)連接數(shù)值子結(jié)構(gòu)與物理子結(jié)構(gòu)部分。物理子結(jié)構(gòu)采用振動(dòng)臺(tái)作為加載裝置。本文試驗(yàn)?zāi)康脑谟隍?yàn)證基于Python?GPU求解的實(shí)時(shí)混合試驗(yàn)系統(tǒng)的性能及可行性，為避免物理子結(jié)構(gòu)試件建模誤差，試驗(yàn)中物理子結(jié)構(gòu)在xPC實(shí)時(shí)系統(tǒng)中仿真求解。當(dāng)振動(dòng)臺(tái)系統(tǒng)接收到數(shù)值子結(jié)構(gòu)的界面響應(yīng)后振動(dòng)臺(tái)空臺(tái)運(yùn)行，測量振動(dòng)臺(tái)臺(tái)面位移及加速度，傳輸至xPC中進(jìn)行物理子結(jié)構(gòu)仿真。物理子結(jié)構(gòu)仿真計(jì)算的動(dòng)力響應(yīng)通過信號(hào)處理部分傳輸至數(shù)值子結(jié)構(gòu)部分，由此構(gòu)成如圖6所示的基于Python?GPU求解的實(shí)時(shí)混合試驗(yàn)系統(tǒng)，圖中試驗(yàn)子結(jié)構(gòu)仿真部分為求解動(dòng)力平衡方程，M為試驗(yàn)子結(jié)構(gòu)質(zhì)量；C為試驗(yàn)子結(jié)構(gòu)阻尼；K為試驗(yàn)子結(jié)構(gòu)剛度；ü為加速度；u˙為速度；u為位移；üt為外荷載加速度激勵(lì)。試驗(yàn)中所用振動(dòng)臺(tái)照片如圖7所示。

3.1 試驗(yàn)參數(shù)

試驗(yàn)中所用模型為本文2.2節(jié)土?結(jié)相互作用模型的縮尺模型，采用Wang等［24］提出的縮尺方法，縮尺比為104。數(shù)值子結(jié)構(gòu)部分由ABAQUS有限元軟件前處理得到，通過添加腳本提取模型參數(shù)并使用MATLAB將參數(shù)轉(zhuǎn)換為矩陣形式。試驗(yàn)子結(jié)構(gòu)部分為單自由度結(jié)構(gòu)，采用Newmark?β算法仿真，質(zhì)量ms=1.45 kg，剛度ks=1003 N/m，阻尼cs=0.374 N/（m?s?1）。試驗(yàn)在振動(dòng)臺(tái)臺(tái)面布置位移及加速度傳感器，對(duì)振動(dòng)臺(tái)輸入白噪聲信號(hào)并進(jìn)行振動(dòng)臺(tái)系統(tǒng)辨識(shí)，使用3階傳遞函數(shù)辨識(shí)振動(dòng)臺(tái)特性，振動(dòng)臺(tái)傳遞函數(shù)如下式所示：

對(duì)振動(dòng)臺(tái)輸入白噪聲位移信號(hào)，振動(dòng)臺(tái)實(shí)際加載情況與期望的響應(yīng)有幅值和相位誤差。為此，在振動(dòng)臺(tái)信號(hào)輸入前添加Tang等［25］提出的FSCS補(bǔ)償控制器，對(duì)振動(dòng)臺(tái)系統(tǒng)進(jìn)行控制補(bǔ)償。補(bǔ)償前后振動(dòng)臺(tái)位移時(shí)程對(duì)比如圖8所示，補(bǔ)償后的振動(dòng)臺(tái)幅值和相位誤差相比未補(bǔ)償時(shí)減小。

3.2 試驗(yàn)結(jié)果

實(shí)時(shí)混合試驗(yàn)中數(shù)值子結(jié)構(gòu)模型自由度和求解時(shí)間步長對(duì)試驗(yàn)精度有影響，本文試驗(yàn)中分別選擇使用GPU與CPU作為求解硬件，采用不同的時(shí)間步長、模型自由度對(duì)數(shù)值模型GPU求解性能進(jìn)行對(duì)比。表3中工況1～8為時(shí)間步長為1 ms，5 ms，10 ms及20 ms時(shí)分別使用GPU與CPU求解數(shù)值子結(jié)構(gòu)模型的最大自由度數(shù)；工況9與工況7對(duì)比求解相同自由度的子結(jié)構(gòu)時(shí)，分別使用GPU與CPU可使用的最小時(shí)間步長。

對(duì)土?結(jié)相互作用模型進(jìn)行整體仿真，使用時(shí)間步長為1 ms，GPU求解2592自由度模型實(shí)時(shí)混合試驗(yàn)與整體仿真試驗(yàn)子結(jié)構(gòu)位移時(shí)程對(duì)比如圖9所示。圖9中混合試驗(yàn)與整體仿真結(jié)果吻合，說明基于Python?GPU的實(shí)時(shí)混合試驗(yàn)系統(tǒng)能滿足試驗(yàn)精度的需要。圖10為時(shí)間步長為1 ms，5 ms，10 ms以及20 ms分別使用GPU與CPU求解條件下，實(shí)時(shí)混合試驗(yàn)與整體仿真的試驗(yàn)子結(jié)構(gòu)位移時(shí)程對(duì)比，8組工況條件試驗(yàn)與仿真結(jié)果一致，說明此試驗(yàn)系統(tǒng)精度較好，除求解硬件外其他試驗(yàn)條件相同時(shí)，使用GPU求解可以實(shí)現(xiàn)更大規(guī)模的實(shí)時(shí)混合試驗(yàn)。

在實(shí)時(shí)混合試驗(yàn)中，數(shù)值子結(jié)構(gòu)求解采用數(shù)值積分算法，數(shù)值積分算法求解的精度與計(jì)算步長有關(guān)，求解算法的精度隨著時(shí)間步長的增大而下降，過大的積分步長會(huì)導(dǎo)致條件穩(wěn)定的數(shù)值積分算法計(jì)算發(fā)散。因此，算法的時(shí)間步長不能無限增大，在實(shí)時(shí)混合試驗(yàn)中時(shí)間步長最大取20 ms。為測試增大時(shí)間步長對(duì)試驗(yàn)精度的影響，試驗(yàn)中選擇相同自由度數(shù)值子結(jié)構(gòu)模型，不同時(shí)間步長進(jìn)行測試。本試驗(yàn)中使用表1配置的CPU求解3549自由度模型，需要時(shí)間步長為20 ms，使用GPU求解相同模型僅需1 ms。使用1 ms，20 ms時(shí)間步長進(jìn)行實(shí)時(shí)混合試驗(yàn)，試驗(yàn)結(jié)果與1 ms整體模型仿真對(duì)比如圖11所示。圖11中時(shí)間步長1 ms基于GPU求解的試驗(yàn)結(jié)果與整體仿真一致，時(shí)間步長20 ms基于GPU求解的試驗(yàn)結(jié)果與整體仿真峰值誤差為2.88%，能滿足實(shí)時(shí)混合試驗(yàn)的精度要求。董曉輝等［17］對(duì)實(shí)時(shí)混合試驗(yàn)中不同步長對(duì)中心差分法精度的影響進(jìn)行測試，20 ms時(shí)間步長試驗(yàn)結(jié)果相比1 ms時(shí)間步長整體仿真誤差達(dá)到10.28%。時(shí)間步長對(duì)無條件穩(wěn)定的Newmark?β算法精度的影響相比條件穩(wěn)定的中心差分法小，使用Newmark?β算法與大時(shí)間步長（20 ms）求解大規(guī)模數(shù)值模型仍具有較高的精度。

3.3 試驗(yàn)結(jié)果分析

由圖10（a）～（b）可見，當(dāng)時(shí)間步長為1 ms時(shí)，使用CPU求解的實(shí)時(shí)混合試驗(yàn)數(shù)值子結(jié)構(gòu)模型為1215自由度，使用GPU求解能達(dá)到5400自由度，使用GPU能較大幅度提升數(shù)值子結(jié)構(gòu)實(shí)時(shí)求解規(guī)模。得益于UDP網(wǎng)絡(luò)的高速、高質(zhì)量通訊，使得數(shù)值子結(jié)構(gòu)與試驗(yàn)子結(jié)構(gòu)實(shí)時(shí)交互界面響應(yīng)，盡可能降低因信號(hào)通訊造成的精度影響。相較于使用LabVIEW作為信號(hào)通訊軟件并通過模擬信號(hào)通訊的實(shí)時(shí)混合試驗(yàn)方法［17］，此試驗(yàn)系統(tǒng)架構(gòu)可實(shí)現(xiàn)更小的GPU實(shí)時(shí)求解時(shí)間步長，優(yōu)化試驗(yàn)系統(tǒng)架構(gòu)，使得大規(guī)模數(shù)值模型實(shí)時(shí)求解實(shí)時(shí)混合試驗(yàn)成為可能，更小的時(shí)間步長可以提升試驗(yàn)精度。

由圖10（g）～（h）可見，當(dāng)時(shí)間步長為20 ms時(shí)，使用CPU求解的實(shí)時(shí)混合試驗(yàn)數(shù)值子結(jié)構(gòu)模型為3549自由度，使用GPU求解時(shí)為24375自由度。此時(shí)GPU的求解性能相比CPU有了顯著的提升，實(shí)現(xiàn)了CPU無法實(shí)時(shí)求解的大規(guī)模模型實(shí)時(shí)混合試驗(yàn)。使用LabVIEW?MATLAB的試驗(yàn)架構(gòu)實(shí)現(xiàn)基于條件穩(wěn)定的中心差分法20 ms時(shí)間步長24000自由度的實(shí)時(shí)混合試驗(yàn)，由于算法的條件穩(wěn)定性，對(duì)數(shù)值模型具有一定限制［17］。本試驗(yàn)采用無條件穩(wěn)定的Newmark?β算法，對(duì)模型參數(shù)無條件限制，相比條件穩(wěn)定的算法進(jìn)行實(shí)時(shí)混合試驗(yàn)拓展了應(yīng)用范圍。

分別使用1 ms和20 ms的時(shí)間步長，Newmark?β算法求解3549自由度數(shù)值子結(jié)構(gòu)模型，記錄其試驗(yàn)子結(jié)構(gòu)位移時(shí)程結(jié)果。將1 ms和20 ms時(shí)間步長的混合試驗(yàn)結(jié)果與1 ms時(shí)間步長的整體仿真結(jié)果對(duì)比如圖11所示。20 ms時(shí)間步長條件下使用無條件穩(wěn)定的Newmark?β算法具有良好的精度，時(shí)間步長對(duì)Newmark?β算法的影響較小。因此，在基于GPU求解的實(shí)時(shí)混合試驗(yàn)中可以使用較大的時(shí)間步長（20 ms），實(shí)時(shí)求解更大規(guī)模的數(shù)值子結(jié)構(gòu)模型。

4 結(jié) 論

本文使用Python作為編程語言，使用GPU建立了實(shí)時(shí)混合試驗(yàn)系統(tǒng)，通過仿真與試驗(yàn)驗(yàn)證了GPU求解數(shù)值子結(jié)構(gòu)模型的性能，仿真與試驗(yàn)結(jié)果如下：

（1）相同時(shí)間步長下，使用Python?GPU相比于傳統(tǒng)CPU可以實(shí)時(shí)求解更大規(guī)模的數(shù)值子結(jié)構(gòu)模型。使用本文選擇的GPU和CPU作為求解硬件，無條件穩(wěn)定的Newmark?β算法在20 ms時(shí)間步長條件下，可實(shí)現(xiàn)24375自由度模型的實(shí)時(shí)混合試驗(yàn)。

（2）相同數(shù)值子結(jié)構(gòu)自由度條件下，使用Python?GPU相比于CPU可以實(shí)現(xiàn)更小積分步長的實(shí)時(shí)混合試驗(yàn)。使用本文選擇的GPU和CPU作為求解硬件，數(shù)值積分算法為無條件穩(wěn)定的Newmark?β算法時(shí)，CPU求解3549自由度模型時(shí)需要時(shí)間步長為20 ms，而GPU求解僅需1 ms的時(shí)間步長，更小積分步長有利于提高實(shí)時(shí)混合試驗(yàn)精度。

（3）Python?GPU相比于LabVIEW?MATLAB的GPU試驗(yàn)架構(gòu)，提升了子結(jié)構(gòu)之間的數(shù)據(jù)通訊效率，使得基于GPU求解的實(shí)時(shí)混合試驗(yàn)最小步長由4 ms降低至1 ms，并且實(shí)現(xiàn)了無條件穩(wěn)定Newmark?β算法的GPU加速，避免了算法穩(wěn)定性對(duì)模型的限制。時(shí)間步長增加對(duì)Newmark?β算法的影響相比中心差分法小，Newmark?β算法求解大規(guī)模數(shù)值模型時(shí)使用大時(shí)間步長仍有較高精度，進(jìn)一步拓展了實(shí)時(shí)混合試驗(yàn)的應(yīng)用范圍。本文選擇的阻尼模型為Rayleigh阻尼，數(shù)值積分算法為中心差分法與Newmark?β算法，若采用更高效的阻尼模型和數(shù)值求解算法，與GPU并行計(jì)算相結(jié)合有望進(jìn)一步提升模型實(shí)時(shí)求解的效率。

參考文獻(xiàn)

1Nakashima M， Kato H， Takaoka E. Development of real‐time pseudo dynamic testing［J］. Earthquake Engineering & Structural Dynamics， 1992， 21（1）： 79-92.

2孟凡濤， 趙建鋒， 于廣明. 實(shí)時(shí)子結(jié)構(gòu)混合試驗(yàn)中的數(shù)值積分方法對(duì)比分析［J］. 地震工程與工程振動(dòng)， 2011， 31（5）： 60-67.

Meng Fantao， Zhao Jianfeng， Yu Guangming. Study on numerical integration methods in real-time hybrid testing experiment［J］. Earthquake Engineering and Engineer Dynamics， 2011， 31（5）： 60-67.

3王倩穎， 吳斌， 歐進(jìn)萍. 考慮作動(dòng)器時(shí)滯及其補(bǔ)償?shù)膶?shí)時(shí)子結(jié)構(gòu)實(shí)驗(yàn)穩(wěn)定性分析［J］. 工程力學(xué)， 2007， 24（2）： 9-14.

Wang Qianying， Wu Bin， Ou Jinping. Stability analysis of real-time substructure testing considering actuator delay and compensation［J］. Engineering Mechanics， 2007， 24（2）： 9-14.

4洪越， 唐貞云， 何濤， 等. 大尺寸非線性實(shí)時(shí)動(dòng)力子結(jié)構(gòu)試驗(yàn)實(shí)現(xiàn)［J］. 振動(dòng)工程學(xué)報(bào)， 2017， 30（6）： 913-920.

Hong Yue， Tang Zhenyun， He Tao， et al. The implementation of nonlinear real-time dynamics substructuring for large scale specimen［J］. Journal of Vibration Engineering， 2017， 30（6）：913-920.

5Nakashima M， Ishida M， Ando K. Integration techniques for substructure pseudo dynamic test： pseudo dynamic test using substructuring techniques［J］. Journal of Structural and Construction Engineering（Transactions of AIJ）， 1990， 417： 107-117.

6Chang S Y. Explicit pseudodynamic algorithm with unconditional stability［J］. Journal of Engineering Mechanics， 2002， 128（9）： 935-947.

7Chang S Y. A new family of explicit methods for linear structural dynamics［J］. Computers & Structures， 2010， 88（11-12）： 755-772.

8吳斌， 保海娥. 實(shí)時(shí)子結(jié)構(gòu)實(shí)驗(yàn)Chang算法的穩(wěn)定性和精度［J］. 地震工程與工程振動(dòng)， 2006， 26（2）： 41-48.

Wu Bin， Bao Haie. Stability and accuracy of Chang algorithm for real-time substructure testing［J］. Earthquake Engineering and Engineering Vibration， 2006， 26（2）： 41-48.

9Chen C， Ricles J M. Development of direct integration algorithms for structural dynamics using discrete control theory［J］. Journal of Engineering Mechanics， 2008， 134（8）： 676-683.

10Tang Y， Lou M L. New unconditionally stable explicit integration algorithm for real-time hybrid testing［J］. Journal of Engineering Mechanics， 2017， 143（7）：04017029.

11Nakashima M， Masaoka N. Real-time on-line test for MDOF systems［J］. Earthquake Engineering & Structural Dynamics， 1999， 28（4）： 393-420.

12Chae Y， Kazemibidokhti K， Ricles J M. Adaptive time series compensator for delay compensation of servo‐hydraulic actuator systems for real‐time hybrid simulation［J］. Earthquake Engineering & Structural Dynamics， 2013， 42（11）： 1697-1715.

13Zhu F， Wang J T， Jin F， et al. Simulation of large-scale numerical substructure in real-time dynamic hybrid testing［J］. Earthquake Engineering and Engineering Dynamics， 2014， 13（4）： 599-609.

14Lu L Q， Wang J T， Zhu F. Improvement of real-time hybrid simulation using parallel finite-element program［J］. Journal of Earthquake Engineering， 2020， 24（10）： 1547-1565.

15劉曉強(qiáng). 面向土木工程數(shù)值模擬的可視化方法研究［D］. 北京： 清華大學(xué)， 2013.

Liu Xiaoqiang. Research on visualization method for civil engineering numerical simulation［D］. Beijing： Tsinghua University， 2013.

16解琳琳， 韓博， 許鎮(zhèn)， 等. 基于OpenSees的大型結(jié)構(gòu)分析GPU高性能計(jì)算方法［J］. 土木建筑工程信息技術(shù)， 2014， 6（5）： 22-25.

XIE Linlin， HAN Bo， XU Zhen， et al. GPU powered high-performance computing method for the analysis of large-scale structures based on OpenSees［J］. Journal of Information Technology in Civil Engineering and Architecture， 2014， 6（5）： 22-25.

17董曉輝， 唐貞云， 李振寶， 等. 應(yīng)用GPU求解的實(shí)時(shí)子結(jié)構(gòu)試驗(yàn)架構(gòu)與性能驗(yàn)證［J］. 振動(dòng)工程學(xué)報(bào)，2022，35（1）：64-71.

Dong Xiaohui， Tang Zhenyun， Li Zhenbao， et al. Performance verification of GPU-based framework for real-time hybrid testing［J］. Journal of Vibration Engineering，2022，35（1）：64-71.

18Okuta R， Unno Y， Nishino R， et al. CuPy： a NumPy-compatible library for NVIDIA GPU calculations［C］. 31st Confernce on Neural Information Processing Systems （NIPS 2017）. Long Beach， CA， 2017.

19Goerzen J，Bower T， Rhodes B. Foundations of Python Network Programming［M］. Berkeley： Apress， 2010.

20李彬. 基于Windows的計(jì)算機(jī)數(shù)字控制系統(tǒng)實(shí)時(shí)性的研究［D］. 哈爾濱： 哈爾濱工業(yè)大學(xué)， 2008.

Li Bin. Research on real-time performance of computer numerical control system based on Windows［D］. Harbin： Harbin Institute of Technology， 2008.

21Rosenquist C. Hard realtime rapid prototyping development platform［D］. Link?ping： Link?ping University， 2003.

22進(jìn)兵. 基于xPC Target的無人機(jī)飛行控制軟件快速原型設(shè)計(jì)［D］. 南京： 南京航空航天大學(xué)， 2008.

Jin Bing. Rapid prototype design for flight control software of UAV based on xPC Target［D］. Nanjing： Nanjing University of Aeronautics and Astronautics， 2008.

23Tang Z Y， Dietz M， Hong Y， et al. Performance extension of shaking table-based real-time dynamic hybrid testing through full state control via simulation［J］. Structural Control & Health Monitoring， 2020， 27（10）： e2611.

24Wang Qiang， Wang Jinting， Jin Feng， et al. Real-time dynamic hybrid testing for soil-structure interaction analysis［J］. Soil Dynamics and Earthquake Engineering， 2011， 31（12）： 1690-1702.

25Tang Z Y， Dietz M， Li Z B， et al. The performance of delay compensation in real-time dynamic substructuring［J］. Journal of Vibration and Control， 2018， 24（21）： 5019-5029.

Research on real-time hybrid simulation by Python-GPU computing

DONG Xiao-hui 1 ?TANG Zhen-yun 1，2 ?LI Zhen-bao 1DU Xiu-li 1

1. The Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education， Beijing University of Technology， Beijing 100124， China；

2. Hebei Key Laboratory of Earthquake Disaster Prevention and Risk Assessment， Sanhe 065201， China

Abstract This article establishes an RTHS framework based on Python and graphics processor unit （GPU） computation. Using a soil-structure interaction system as a testing model， the performance of the testing framework is verified by numerical simulation and experiment. The results show that using Python-GPU to speed the unconditionally stable algorithm allows a numerical substructure with 24000 DOFs to be solved in a time step of 20 ms. Because the scale of GPU calculation is almost 7 times greater than that of CPU calculation in the same machine， the capability of RTHS is greatly extended.

Keywords real-time hybrid simulation; graphics processing unit; Python; numerical integration algorithm; calculation efficiency