考慮數值積分算法的實時子結構試驗穩定性研究

2021-11-17 07:43:02魯旭杰周子豪李忠獻

工程力學 2021年11期

李寧，魯旭杰，周子豪，李忠獻

(1. 天津大學建筑工程學院，天津 300350；2. 濱海土木工程結構與安全教育部重點實驗室(天津大學)，天津 300350；3. 中國地震局地震工程綜合模擬與城鄉抗震韌性重點實驗室(天津大學)，天津 300350)

近年來，隨著隔震、減震/振結構中，采用動力特性與變形速率相關的新裝置大量涌現，傳統擬/靜力試驗方法無法適應測試和研究需要，快速或實時混合模擬(real-time hybrid simulation，RTHS)方法[1-2]正推廣用于結構動力性能研究。RTHS試驗常將原結構劃分為2 部分，其中本構關系復雜、無法準確建立模型或使用數值模擬的結果不理想的部分一般作為物理子結構，用作動器或振動臺等加載裝置進行加載；其余模擬方法成熟、較為準確的部分作為數值子結構，使用模擬軟件模擬。兩部分之間進行同步或高速數據交互。當反饋與命令達到實時交互時，也稱為實時子結構試驗(real-time substructure test，RTST)。

以圖 1 所示單層框架為例，RTST 在試驗開始時，首先輸入地震動，由數值子結構計算出命令位移，作動器或振動臺根據命令施加動邊界給物理子結構，再通過物理子結構上的傳感器測得其反應，通過DAQ 將反應量傳輸給數值子結構，兩者進行數據交互。如此循環往復，直至試驗結束。

圖1 實時混合模擬試驗Fig. 1 Real time hybrid simulation

近兩年，國內外學者通過實時混合試驗對多種結構的特性進行了試驗研究。例如，李騰飛等[3]通過基于OpenFresco 試驗平臺的混合試驗系統對高強鋼組合K 形偏心支撐鋼框架的抗震性能進行了研究；Mukai 等[4]應用實時混合試驗技術通過雙作動器對混凝土框架結構在主動調諧質量阻尼器控制下的結構性能進行了研究。Najafi 等[5]用實時混合試驗對配有非線性能量阱結構的動力特性進行了研究。

RTST 目前的研究熱點有：時滯補償、復雜邊界或加載方法，以及穩定性分析等。RTST 成功與否的關鍵是保證系統的穩定[6]。在RTST 中，由于數值子結構的求解屬于離散時間系統，物理子結構的響應以及作動器或振動臺等加載元件的運動是連續時間系統，RTST 應屬于連續-離散混合系統。對此，建立對應的分析模型和理論還是國內外熱點問題之一。目前主要采用連續時間系統或離散時間系統兩類方法對RTST 進行穩定性探究。

采用連續系統模型，Wallace 等[7]使用時滯微分方程導出了單自由度系統臨界時滯表達式，給出其精確解。Kyrychko 等[8]將單自由度系統臨界時滯計算方法推廣到了多自由度系統。為了全面考慮結構劃分對RTST 的影響，Maghareh 等[9]采用歐拉變換將時滯微分方程化簡，得到不同結構劃分方式(不同質量比、剛度比、阻尼)對系統穩定性的影響規律，認為結構劃分在設計RTST 時必須考慮。

在RTST 中數值子結構需采用積分算法逐步求解，應考慮積分步長和積分算法假定對RTST穩定性的影響。許多學者對此問題進行了研究，例如：基于離散系統模型，Chen 等[10]結合作動器時滯假設和離散系統控制原理，得到了考慮數值積分算法后的系統穩定邊界；王倩穎等[11]使用譜半徑法研究了使用Operator-Split 積分算法的RTST試驗穩定性隨結構參數的變化；鄧麗霞等[12]引入速度和加速度的附加假定，認為穩定性必須考慮試驗和物理子結構質量之比的影響；Zhu 等[13]使用離散根軌跡方法研究了不同數值積分算法對多自由度RTST 系統時滯穩定性影響，認為時滯假設在低頻時對加載系統的動力特性描述更精確；唐貞云等[14]通過對電液伺服加載系統進行建模，討論了試件與加載系統對RTST 穩定性的影響；郭珺等[15]結合振型疊加法和增益裕度，認為時滯補償對RTST 系統穩定性的影響不能確定。但是，既有的離散穩定性分析方法也存在一些不足，例如：Chen 等[10]提出的離散方法引用的時滯模型使其只適用于線彈性剛度和小時滯的情況；Zhu 等[13]提出的離散方法需要將要研究模型劃分比例寫成系統傳遞函數增益的形式進而通過根軌跡方法進行穩定性分析，使其局限為只能對質量劃分比例進行研究的方法。這也表明離散系統穩定性分析方法的探索有待深入，也需要一種較通用的分析方法。

本文基于離散控制理論[16]，建立了引入積分算法的RTST 時滯系統傳遞函數，進而對系統傳遞函數的極點分布和系統的穩定性進行分析。在此基礎上：一方面對RTST 中數值和物理子結構劃分導致的系統穩定性的變化規律進行了詳細的討論，給出了其穩定區域以及對應的臨界時滯，為RTST 中數值和物理子結構的劃分及時滯補償的可行性提供了依據；另一方面，本文也使用了其他兩種連續系統假定的穩定性分析方法與本文所提方法得到的穩定性分析結果進行了對比，說明了本文方法更為安全可靠。

1 穩定性分析理論概要

1.1 穩定性分析概要

實時子結構系統可表述為一個閉環系統，如圖 2 所示。

圖2 實時子結構閉環控制系統的框圖Fig. 2 Block diagram representation of real-time hybrid testing of an SDOF system

閉環控制系統的性能是由閉環傳遞函數Gcl(z)決定的：

1.2 引入數值積分算法的穩定性分析方法

使用基于極點分布的穩定性分析方法對SDOF的RTST 穩定性進行分析，其運動方程如下：

表1 傳遞函數系數Table 1 Transfer function coefficients

1.3 基于連續時間假設的臨界時滯精確解

式中：j=1, 2, ···；k=0, 1, 2, ···。式(14)中正的臨界時滯解才具有物理意義，且最小的正數解對RTST系統的穩定性最為關鍵，因為其代表了RTST 系統初次從穩定狀態進入到不穩定狀態。

1.4 基于連續時間假設的臨界時滯近似解

2 考慮積分算法的RTST 穩定性分析

當系統時滯為0 時，RTST 的傳遞函數的穩定性由積分算法自身的穩定性控制(CR 積分算法此時無條件穩定)。當時滯不為0 時，考慮到RTST中時滯的真實情況，本文研究的時滯范圍為0 s～1 s，并將結構在時滯范圍為0 s～1 s 內都穩定的現象定義為“無條件穩定”。

某 SDOF 的質量M=52.56 kg，剛度K=87 600 N/m，阻尼比ζ=0.03，物理子結構的加載步長與數值子結構的積分步長都為0.01 s。本節研究物理子結構與數值子結構的質量比、剛度比，以及物理子結構與原結構的阻尼比作為自變量時，子結構劃分對RTST 系統穩定性的影響。

2.1 質量之比作為自變量的穩定性探究

當以物理子結構、數值子結構的質量比作為變量時，令物理子結構的剛度及阻尼為0，此時其結構簡圖如圖3(a)所示，使用勞斯-赫爾維茨穩定性判據可得該系統近似臨界時滯的表達式為：

式(18)始終為負，故對于該工況，近似方法失效。該工況下的時滯穩定域圖像如圖3(b)所示。

圖3 實時子結構的時滯穩定域(Ke=0，Ce=0)Fig. 3 Delay stability region of RTST (Ke=0, Ce=0)

圖 3(b)中綠色區域代表通過離散方法得到的穩定域，其他區域代表系統不穩定。可以看出，當質量比小于0.062 時，系統為無條件穩定，不受時滯條件影響。隨著物理、數值子結構的質量比Me/Mn增加，該系統的穩定性顯著下降，且當0.062＜Me/Mn＜0.427 時，系統出現穩定性切換現象，即隨著時滯增加，系統在穩定狀態與不穩定狀態之間變換。精確方法繪制出的臨界時滯曲線的上界全部收斂至Me/Mn=1 處，且穩定區域的上界也為1，說明了該工況下當且僅當Mn＞Me時，系統才可能穩定，反之，系統始終失穩。

雖然精確方法和離散方法總體上得到的穩定區域基本一致，但由局部放大圖可以看出，離散方法得到的穩定域更小，說明了離散方法較精確方法更保守、接近實際。

2.2 剛度之比作為自變量的穩定性探究

當以物理、數值子結構的剛度比作為變量時，研究了如下4 種工況：1)Me=0；2)Me/Mn=0.5 ；3)Me/Mn=1 ；4)Me/Mn=2。圖 4(a)為工況1 所對應的結構簡圖，圖 4(b)為該工況的時滯穩定域。

如圖 4 所示，隨著物理、數值子結構的剛度比Ke/Kn增加，該系統的穩定域越來越小，穩定性越來越差，說明物理剛度的引入會降低RTST穩定性，而當Ke/Kn＜0.062時，該系統為無條件穩定，當 0.062 ＜Ke/Kn＜0.625時，該系統同樣會出現穩定性切換現象。當Ke/Kn=1.5時，只要時滯足夠小，通過連續方法也可繪制出臨界時滯曲線，證明此時系統也是穩定的。說明了剛度比并不存在絕對的穩定上界。這一點與質量比顯然不同。當時滯遠小于1 ms 時，近似方法與精確方法繪制出的臨界時滯邊界十分接近；時滯較大時，近似方法繪制的臨界時滯邊界會迅速衰減，說明近似方法僅適用于小時滯情況；當時滯較大時，近似方法無法準確判斷穩定性。此外，相對于連續時間系統穩定性分析方法，離散方法更加保守安全。

圖4 實時子結構的時滯穩定域(Me=0，Ce=0)Fig. 4 Delay stability region of RTST (Me=0, Ce=0)

圖5 分別為工況2～工況4 對應的穩定域圖像。對比三幅圖所繪制出穩定區域的大小可以看出，Me/Mn=0.5 時穩定域最大，Me/Mn=1得到的穩定域要小于Me/Mn=0.5時的穩定域，而當Me/Mn=2時，系統為無條件不穩定。一方面，證明了前述得到的“RTST 要保持穩定必須要滿足質量比小于1”的結論(Me/Mn＜1)；另一方面，也說明物理子結構比重的增加會破壞RTST系統的穩定性，Me/Mn越大，系統越容易失穩。同時，Me/Mn=0.5和1 的工況下，結構無條件穩定域分別出現在-0.078 ＜Me/Mn-Ke/Kn＜0.075和Ke/Kn≈1 ，因此，可以得出結論，在Me/Mn＜1的條件下，當且僅當Me/Mn-Ke/Kn≈0時，系統有無條件穩定特性。與只考慮物理子結構的質量或只考慮物理子結構的剛度比較，可以發現同時考慮物理子結構的剛度和質量后，系統不再是簡單的隨著質量比或者剛度比增大而穩定性降低。說明了該系統受質量比與剛度比這兩個參數的耦合作用影響。

圖5 實時子結構的時滯穩定域(Ce=0)Fig. 5 Delay stability region of RTST (Ce=0)

2.3 阻尼之比作為自變量的穩定性探究

當以物理子結構與原結構的阻尼之比作為變量時，同樣也研究了四種工況：1)Ke=0 ；2)Ke/Kn=0.5 ；3)Ke/Kn=1 ；4)Ke/Kn=2。

圖 6(b)為工況1 僅考慮物理子結構的阻尼時結構的穩定域圖像，可以看出，當阻尼之比Ce/C＜0.5時，系統始終是無條件穩定的，當Ce/C＞0.5后，系統隨時滯增加呈現出周期性的穩定性切換現象，同時，穩定區域隨著時滯增加越來越小。對該工況，并不存在小時滯情況，故近似方法并不適用于該工況。其余兩種方法對比的結論與第2.2 節剛度之比為自變量時，工況1 的穩定性變化規律類似，不再贅述。

圖6 實時子結構的時滯穩定域(Me=0, Ke=0)Fig. 6 Delay stability region of real time substructure (Me=0, Ke=0)

圖7 分別為工況2～工況4 對應的穩定域圖像，對比三幅圖可以看出，隨著Ke/Kn不斷增大，系統的穩定域越來越小，即結構越來越容易失穩。由于研究中選取的積分步長為0.01 s，使用離散方法進行研究時，時滯最小只能取至與積分步長相同，故該積分步長下，離散方法無法直接判斷時滯小于10 ms 時結構是否穩定，而當時滯小于10 ms 時，三幅圖中都存在一條使用精確方法繪制出的臨界時滯曲線，由于該區間內時滯較小，故認為這條曲線左側都是穩定的，并用陰影區域表示，即該穩定區域是通過精確方法獲得的。同時，可以觀察到每幅圖中的臨界時滯曲線近似可以看做多條平行的豎向線段，即系統的穩定性受阻尼之比影響很小，反而對剛度比更加敏感。故可以認為在設計實時子結構試驗時，阻尼之比不起到決定性的影響作用。

圖7 實時子結構的時滯穩定域(Me=0)Fig. 7 Delay stability region of real time substructure (Me=0)

可知精確方法得到的最左側的臨界時滯曲線與近似方法繪制的近似臨界時滯曲線吻合程度較高。充分說明近似方法只能滿足小時滯情況，當時滯較大時，近似方法與其他兩種方法都無法準確判定穩定性。

2.4 阻尼比和積分時間步長對穩定性的影響探究

為說明阻尼比與積分步長對單自由度實時子結構穩定性的影響，取圖4(a)結構模型進行研究。圖 8 為阻尼比分別為0.03、0.05 和0.10 時的時滯穩定域圖像。可以看出不同阻尼比的結構穩定域輪廓基本一致，但系統的穩定域會隨著阻尼比的增大而增大，說明增大結構的阻尼可以提高系統的穩定性。

圖8 不同阻尼比條件下結構穩定域Fig. 8 Stability region of an SDOF system with different damping ratios

圖9 為積分步長分別為0.01 s、0.02 s 和0.03 s時結構的時滯穩定域圖像。從圖9 可以看出，當積分步長為0.01 s 時，精確法和離散方法繪制的穩定域基本一致，隨著積分步長增大到0.02 s 和0.03 s后，離散方法得到的穩定區域顯著降低，與連續方法得到的穩定區域差異也逐漸增大。說明了當積分步長小于或等于0.01 s 時，不需要考慮積分算法對穩定性的影響；反之，當積分步長大于0.01 s后，就不得不考慮積分算法對穩定性的改變，此時應該使用更為保守安全的離散穩定性分析方法。

圖9 不同積分步長下的結構穩定域Fig. 9 Stability region of an SDOF system with different with different integration steps

3 數值驗證

為了驗證本文得到的結果，本文搭建了單自由度實時子結構仿真模型進行仿真驗證，Simulink模型如圖 10 所示。

圖10 時滯系統Simulink 模型Fig. 10 Simulink model of delay system

取物理子結構僅考慮剛度的工況繪制的穩定域圖像上“●”標記點進行數值模擬，如圖 11 所示，即剛度比Ke/Kn=0.3，時滯為80 ms、130 ms和200 ms。積分步長取1/1024 s，激勵方式為單位脈沖，得到結構的脈沖響應如圖 12 所示。

圖11 數值模擬驗證工況Fig. 11 Numerical simulation verification case

圖12 反應的是結構在不同時滯情況下的單位脈沖響應。結構位移響應呈現發散的趨勢，則說明系統是不穩定的。從圖中可以看出該單自由度實時子結構系統在時滯為80 ms 和200 ms 時位移是發散的，在時滯為130 ms 時位移是收斂的。即系統在80 ms 是失穩的，當時滯增大至130 ms，系統又進入了穩定狀態，當時滯進一步增大至200 ms系統再次失穩。說明隨著時滯增大，系統出現了穩定性切換現象。該數值模擬得到的結論與理論分析得到的結論相同，證明了本文理論的可靠性。

圖12 結構的單位脈沖響應Fig. 12 Impulse response of structure

4 試驗驗證

4.1 SDOF 穩定性驗證試驗

本文進行了單時滯源的RTST 實測。試驗設備為濱海土木工程與安全教育部重點試驗室的MTS伺服液壓作動器，該作動器本身存在時滯，時滯大小約為21/1024 s。試驗中的數值計算通過控制器完成，積分步長取1/1024 s，積分方法為CR 積分算法。如圖 13，試驗模型為SDOF 結構，取結構的部分剛度作為物理子結構。考慮到失穩會造成構件和試驗設備造成不必要的損失，故物理子結構也通過數值建模，作動器加載端不放置實際的構件。而物理子結構為線彈性結構時，物理子結構的反饋力即彈性力，因此只需要將作動器的位移信號作為反饋信號反饋給數值子結構即可。

圖13 實時子結構試驗加載裝備Fig. 13 Facility for real-time hybrid testing

采用反補償方法[19]對作動器時滯進行補償，補償后該實時子結構系統的時滯減小至約3 ms。由于結構的體量較小，且物理子結構是通過數值建模的線彈性結構，因此可以在位移反饋信號中使用Simulink 提供的delay 模塊，進而達到改變時滯大小的效果。取結構參數如表 2，此時結構劃分情況為Ke/Kn=1/10。使用本文提出方法計算并繪制出不同時滯下控制極點分布圖，如圖 14 所示。

表2 結構參數Table 2 Structural parameters

圖14 不同時滯下控制極點模長圖Fig. 14 Norm of the control poles with different time delay

圖14 中縱坐標為控制極點的模長，橫坐標為系統時滯，控制極點模長大于1 時系統失穩，模長小于1 時，系統穩定。由圖14 可以看出在該劃分情況下，系統時滯在9 ms～63 ms 是失穩的，其余時滯區間內是穩定的，故取圖14 上標記的6 個不同時滯進行實測，試驗中結構參數和工況分別如表 2 和表 3，試驗激勵為周期為4 s 的正弦激勵。

表3 試驗工況Table 3 Test conditions for SDOF test

圖 15 為不同時滯下試驗結構的位移響應，圖中深實線為命令位移，淺虛線為測量位移。可以看出，當時滯為4 ms 時，結構的位移響應呈現出等幅震蕩的現象但沒有發散，說明該時滯下系統接近于臨界穩定；當時滯增大至20 ms、35 ms、50 ms 時，這三種時滯下的測量位移都呈現出高頻振蕩和發散現象；而增大時滯至90 ms 時，結構的測量位移收斂。說明當時滯為4 ms 和90 ms時，該單自由度實時子結構系統是穩定的，時滯為20 ms、35 ms 和50 ms 時，該系統失穩。因此，試驗結果與理論分析得到的穩定狀態是一致的。通過該試驗證明了單自由度實時子結構存在穩定性切換現象，隨著時滯增大，該單自由度實時子結構從穩定狀態進入不穩定狀態，當時滯進一步增大后，又從不穩定狀態回到了穩定狀態。圖 15(b)為時滯7 ms 的試驗結果，從圖上可以看出該時滯下結構的的測量位移出現了高頻震蕩的現象，最終在13 s 時失穩了。該試驗結果與理論分析相反，該時滯下結構應該是穩定的，其位移響應不應該呈現發散的狀態。該工況下理論與試驗結果不符可能是因為試驗中存在的噪聲等干擾會對結構的穩定性造成影響，另一方面，作動器的時滯是在小范圍內波動的而不是恒定不變的，故在該工況下，該實時子結構系統的時滯已接近或超過臨界時滯(9 ms)，故導致最終發散。

圖15 試驗結果Fig. 15 Test results

4.2 2DOF 穩定性驗證試驗

在SDOF 試驗的基礎上進行2DOF 穩定性驗證試驗，將2DOF 上層結構作為數值子結構，下層作為數值子結構，如圖 16 所示。物理子結構為一單層框架結構，其參數為：質量Me=5 kg，剛度Ke=5500 N/m，底部應變與頂部剪力的對應關系為f(t)=68 000ε(f(t)為剪力，ε 為柱腳應變)。注意，物理子結構的慣性力應該扣除以得到物理子結構應當反饋給數值子結構的反力真實值。

圖16 2DOF 試驗Fig. 16 2DOF test

試驗工況見表 4，激勵為El Centro 地震動，試驗中積分步長取值為1/1024 s。試驗結果如圖 17 所示。

表4 試驗工況Table 4 Test conditions for 2DOF test

從圖 17 可以看出，當Me/Mn=Ke/Kn=10，2DOF在時滯28 ms 時的測量位移與目標位移幾乎完全吻合，且位移呈現收斂的趨勢，說明2DOF 在該劃分參數和時滯下處于穩定狀態，與理論穩定分析結果一致。當Me/Mn=6、Ke/Kn=3 而時滯為28 ms時，命令位移與測量位移在5 s 前吻合較好，5 s后測量位移呈現高頻振蕩并發散的趨勢，說明該劃分方式相應時滯情況下2DOF 系統是失穩的。當Me/Mn=3、Ke/Kn=3 而時滯為48 ms 時，測量位移與命令位移在2 s 后迅速呈現發散的現象，說明該工況下2DOF 失穩嚴重。三種工況下的試驗結果都與理論分析結果吻合，工況1 和工況2 結果對比可知，當時滯相同時，2DOF 的穩定性會隨著結構劃分情況變化而變化，因此在2DOF 試驗中，適當的調整數值子結構占原結構的比例可以保證試驗的穩定性。

圖17 2DOF 試驗結果Fig. 17 2DOF test results

5 結論

本文建立了RTHS 基于積分算法的離散穩定性分析方法，并使用提出方法與已有的連續穩定性分析方法研究了結構劃分參數、阻尼比、積分步長對結構穩定性的影響，模擬和實測對比了兩類方法的結果，證明了所提出方法的正確性和可靠性。

(2) 對每一種結構劃分方式，會存在多個時滯穩定區域，說明結構隨著時滯的增加，在穩定狀態與不穩定狀態進行變換。即使時滯很大時，結構也可能是穩定的。實際研究中，一般最關心第一次出現穩定性變換時(即第一次發生失穩現象)的臨界時滯。

(3) 對比三種方法可知，近似方法僅適用于小時滯情況(τ＜1 ms)。當積分步長小于或等于0.01 s，不需要考慮積分算法對穩定性的影響，可以直接使用精確法研究結構的穩定性；反之，當積分步長大于0.01 s，必須考慮積分算法對穩定性的影響，此時應該使用更為保守安全的離散穩定性分析方法，本文方法即是。

(4) 經RTHT 試驗檢驗，本文所提出離散化方法對多自由度問題也有正確的效果，但相關研究擬于后續開展。本文方法目前僅對線性系統進行了驗證，后續會對更為復雜的非線性系統展開研究。