吳 斌， 寧西占， 許國山， 王 貞

(1. 哈爾濱工業(yè)大學(xué) 結(jié)構(gòu)工程災(zāi)變與控制教育部重點實驗室，哈爾濱 150090；2. 哈爾濱工業(yè)大學(xué) 土木工程智能防災(zāi)減災(zāi)工業(yè)和信息化部重點實驗室，哈爾濱 150090)

混合試驗是一種將物理實驗和數(shù)值模擬通過在線方式整合到一起的混合模擬方法，在土木工程領(lǐng)域，它起源于擬動力試驗方法[1]。對于大型結(jié)構(gòu)，混合試驗將整體結(jié)構(gòu)拆分為若干個子結(jié)構(gòu)，其中無法準(zhǔn)確數(shù)值模擬、具有復(fù)雜非線性的部分作為物理子結(jié)構(gòu)在實驗室中進(jìn)行加載，剩余部分則進(jìn)行數(shù)值模擬，這就是所謂子結(jié)構(gòu)試驗方法。子結(jié)構(gòu)試驗降低了混合試驗對試驗設(shè)備的需求，可完成大比例尺甚至足尺試驗?zāi)Ｐ偷幕旌显囼?，有效地減小了尺寸效應(yīng)的影響。

子結(jié)構(gòu)試驗的關(guān)鍵問題之一，是物理子結(jié)構(gòu)邊界條件的模擬。在土木工程子結(jié)構(gòu)試驗中，物理子結(jié)構(gòu)的邊界條件一般由作動器實現(xiàn)。嚴(yán)格地來講，作動器的數(shù)量應(yīng)該與邊界自由度的數(shù)量相同。對于簡單邊界，即自由度較少的邊界，如以單根構(gòu)件或減震控制裝置為物理子結(jié)構(gòu)(或物理子系統(tǒng))的情況，子結(jié)構(gòu)試驗取得了成功的應(yīng)用[2-6]。但是對于更多情況，物理子結(jié)構(gòu)的邊界自由度往往大大超過一般結(jié)構(gòu)實驗室的作動器數(shù)量。例如，對于空間框架結(jié)構(gòu)，一個梁柱單元節(jié)點有6個自由度，如果取結(jié)構(gòu)底部若干層作為物理子結(jié)構(gòu)，那么即使對于一個雙向單跨的框架，其邊界自由度也有24個。即使個別大型試驗室有這么多作動器和雙向反力墻，但是豎向作動器的加載能力和反力架也將是對試驗室硬件資源的一個考驗。因此，在很多子結(jié)構(gòu)試驗中，不得不對邊界條件進(jìn)行簡化處理，甚至放棄某些自由度的模擬[7-10]。顯然，這樣做將不同程度地降低試驗精度，甚至可能導(dǎo)致結(jié)構(gòu)破壞模式發(fā)生變化，試驗結(jié)果嚴(yán)重偏離實際反應(yīng)。

我們將物理子結(jié)構(gòu)邊界條件無法完全實現(xiàn)稱為不完整邊界條件問題。為了解決子結(jié)構(gòu)試驗的不完整邊界條件問題，本文提出基于模型更新的在線數(shù)值模擬方法，稱之為考慮不完整邊界條件的混合試驗方法。本文余下部分主要介紹這種新型混合試驗方法的基本原理和所采用的本構(gòu)參數(shù)識別方法，通過一個子結(jié)構(gòu)試驗的數(shù)值模擬算例驗證所提方法的可行性。

1 基本原理

考慮不完整邊界條件的混合試驗方法的主要思想是：針對待模擬的結(jié)構(gòu)建立三套模型，其一是結(jié)構(gòu)的整體數(shù)值模型，其二是物理子結(jié)構(gòu)即試件，其三是物理子結(jié)構(gòu)的數(shù)值模型；將結(jié)構(gòu)整體數(shù)值模型計算出來的物理子結(jié)構(gòu)邊界上部分自由度的位移通過作動器實現(xiàn)；利用物理子結(jié)構(gòu)的位移和反力實測信息，以及物理子結(jié)構(gòu)的數(shù)值模型，采用適當(dāng)識別方法，識別物理子結(jié)構(gòu)的本構(gòu)參數(shù)；假定全結(jié)構(gòu)與物理子結(jié)構(gòu)具有相同的本構(gòu)關(guān)系，因此可用識別得到的本構(gòu)參數(shù)更新全結(jié)構(gòu)數(shù)值模型的本構(gòu)參數(shù)值，進(jìn)而進(jìn)行下一步的結(jié)構(gòu)反應(yīng)計算和參數(shù)識別更新，如此反復(fù)直至試驗結(jié)束。本文提出的新方法與傳統(tǒng)的混合試驗方法最重要的區(qū)別在于，前者需建立整體結(jié)構(gòu)的數(shù)值模型，而后者無整體結(jié)構(gòu)的數(shù)值模型。本文方法相對于傳統(tǒng)方法的最大優(yōu)點在于，可通過整體結(jié)構(gòu)的數(shù)值模型降低物理子結(jié)構(gòu)不完整邊界條件對模擬結(jié)果的影響。

關(guān)于軟件實現(xiàn)，如圖1所示，本文所提考慮不完整邊界條件的新型混合試驗方法由以下4個模塊組成：① 時間積分模塊，用于求解結(jié)構(gòu)運動方程，本文采用Matlab編程；② 整體結(jié)構(gòu)有限元分析模塊，本文利用OpenSees求解整體結(jié)構(gòu)的靜恢復(fù)力；③ 物理子結(jié)構(gòu)分析模塊，采用OpenSees建模，為識別模塊提供在部分邊界自由度位移輸入下物理子結(jié)構(gòu)的恢復(fù)力；④ 本構(gòu)參數(shù)識別模塊，根據(jù)作動器加載物理試件得到的位移及其對應(yīng)的恢復(fù)力和物理子結(jié)構(gòu)分析模塊計算得到的恢復(fù)力，在線識別本構(gòu)模型參數(shù)，采用Matlab編程。新型混合試驗方法的實施步驟為：

圖1 考慮不完整邊界條件的新型混合試驗原理

步驟1采用逐步積分方法求解結(jié)構(gòu)運動方程，獲取整體結(jié)構(gòu)動力自由度上的位移向量dk，并將其發(fā)送給整體結(jié)構(gòu)有限元模型中對應(yīng)的自由度。

步驟5以最新估計的本構(gòu)模型參數(shù)θk更新整體有限元模型中本構(gòu)模型參數(shù)θk-1，再次根據(jù)結(jié)構(gòu)的位移向量dk完成一步靜力分析，并提取相應(yīng)節(jié)點反力Rk返回時間積分模塊。

步驟6重復(fù)步驟1～5，直至試驗結(jié)束。

2 基于梯度的材料本構(gòu)參數(shù)估計

本文方法不限于特定的本構(gòu)模型和參數(shù)估計方法，僅以單軸混凝土本構(gòu)和基于梯度的參數(shù)識別方法為例，說明本文方法的有效性。

2.1 材料本構(gòu)模型

通過材料本構(gòu)關(guān)系直接反應(yīng)結(jié)構(gòu)或構(gòu)件力學(xué)行為的截面纖維模型常用于模擬鋼筋混凝土框架結(jié)構(gòu)，梅竹等[11]采用有限元軟件OpenSees中Concrete01對應(yīng)的Kent-Scott-Park模型[12]作為混凝土的本構(gòu)模型，提出了混凝土結(jié)構(gòu)混合試驗的模型更新方法。本文也采用此本構(gòu)模型，其骨架曲線(見圖2)的數(shù)學(xué)表達(dá)式為

當(dāng)εc≤0.002K

(1)

σc=Kfc[1-Zm(εc-0.002K)]

當(dāng)εc>0.002K

(2)

圖2 OpenSees中的單軸混凝土本構(gòu)模型

式中：fc為非約束區(qū)混凝土峰值應(yīng)力；σc為混凝土的軸向應(yīng)力；εc為混凝土的軸向應(yīng)變；K為約束引起的混凝土強度增加系數(shù)；Zm為應(yīng)變軟化段斜率?？紤]到約束區(qū)與非約束區(qū)混凝土參數(shù)之間關(guān)系，式(1)和式(2)所描述的本構(gòu)模型可表示為如下非線性函數(shù)

σc=φ(εc,fc,ε0,εu,Kuc,Ks1,Ks2,Ks3,ρs)

(3)

Ks1=fyh/fc

(4)

Ks3=0.9fyh/300εu

(5)

式中：ε0，εu分別為非約束區(qū)混凝土峰值應(yīng)變、極限應(yīng)變；Kuc為混凝土極限應(yīng)力和峰值應(yīng)力的比值；ρs為體積配箍率；fyh為箍筋的屈服強度，Ks2=(1～5)Ks1。

2.2 參數(shù)估計方法

參數(shù)識別有很多可用的方法，這里僅以基于梯度的優(yōu)化方法為例來說明本文所提新型混合試驗方法的效果。相應(yīng)的目標(biāo)函數(shù)可表達(dá)為

R(d1:k))

(6)

θ=[fc,εc,εu,Kuc,Ks1,Ks2,Ks3]T

(7)

式中：下標(biāo)k為當(dāng)前時刻；θ為由待識別的本構(gòu)模型參數(shù)組成的決策變量；d1:k為從初始時刻到時刻k的位移向量；R和F分別為實測和預(yù)測的物理子結(jié)構(gòu)部分邊界自由度上恢復(fù)力組成的列向量，其維數(shù)是nk×1，n為邊界自由度個數(shù)。

顯然對基于有限元的混凝土本構(gòu)參數(shù)優(yōu)化問題，恢復(fù)力是參數(shù)θ的非線性函數(shù)。以式(6)為目標(biāo)函數(shù)的參數(shù)識別問題可歸結(jié)為非線性最小二乘問題，然后通過迭代求解。求解式(6)非線性最小二乘問題的迭代方法較多，本文選用Levenberg-Marquardt法。

按Levenberg-Marquardt方法，參數(shù)的迭代表達(dá)式為

(8)

其中,

(9)

(10)

(11)

式中：i為迭代步；α的選取及計算方法參見文獻(xiàn)[13]，δ為一小參數(shù)，本文取10-9。由于梯度函數(shù)無法解析表達(dá)，因此在式(9)中我們采用差分方法近似計算梯度。

3 數(shù)值模擬

以鋼筋混凝土框架結(jié)構(gòu)為例，通過對比傳統(tǒng)混合試驗、模型更新混合試驗、考慮不完整邊界條件的新型混合試驗和OpenSees純數(shù)值分析結(jié)果來驗證本文所提方法的優(yōu)越性和有效性。結(jié)構(gòu)整體數(shù)值模型、虛擬物理子結(jié)構(gòu)模型與Matlab之間采用基于TCP/IP協(xié)議的socket通訊機制進(jìn)行數(shù)據(jù)交互[14]，參數(shù)識別程序和物理子結(jié)構(gòu)數(shù)值模型之間采用文本方式進(jìn)行數(shù)據(jù)交互。

3.1 虛擬混合試驗概況

本文算例取自一正在進(jìn)行的實際混合試驗，其全結(jié)構(gòu)為一榀2層2跨鋼筋混凝土框架結(jié)構(gòu)，如圖3所示，層高分別為3.6 m和3.0 m，跨度6 m，梁柱截面分別為500 mm×500 mm，300 mm×550 mm。為考慮樓板對梁的約束作用，取梁的有效翼緣為1 740 mm。建模時假定基礎(chǔ)與地基剛接，梁柱均采用基于力的非線性梁柱單元，每個單元取5個Gauss-Lobatto積分點，截面為纖維截面；每層質(zhì)量為8 150.4 kg?；炷晾w維和鋼筋纖維的本構(gòu)模型分別為單軸Kent-Scott-Park (Concrete01)和 G-M-P(Steel02)模型，鋼筋和混凝土材料特征點參數(shù)見表1，表中未列參數(shù)建模時采用OpenSees的推薦值。地震動輸入選為El-Centro波，時間間隔為0.01 s，峰值加速度為310 gal，持時為10 s。采用OpenSees作的結(jié)構(gòu)整體純數(shù)值分析結(jié)果作為參考準(zhǔn)確解來評估三種混合試驗的模擬效果。

(a) 框架原型

(b) 梁柱截面尺寸

虛擬混合試驗的物理子結(jié)構(gòu)取為底層邊柱2/3柱高以下部分，幾何相似比為1/2，并施加恒定軸力荷載P=135.8 kN，相應(yīng)的試驗裝置圖如圖4所示。這里假定此柱反彎點在2/3柱高處，邊界取在反彎點處可以避免轉(zhuǎn)角自由度的模擬。顯然，這里物理子結(jié)構(gòu)的邊界條件是不完整，僅水平位移可以完全實現(xiàn)。為模擬真實物理試驗過程，采用均值為零的隨機數(shù)作為物理子結(jié)構(gòu)恢復(fù)力的觀測噪聲，其幅值為恢復(fù)力峰值的1%。

3.2 模擬結(jié)果

三種虛擬混合試驗?zāi)M結(jié)果與參考解的結(jié)構(gòu)首層響應(yīng)如圖5所示，定量誤差見表2。從圖5和表2可以看出，傳統(tǒng)混合試驗與參考解的差異最大，模型更新混合試驗次之，本文新方法與參考解幾乎完全重合，顯示了本文方法的優(yōu)越性。傳統(tǒng)混合試驗的數(shù)值子結(jié)構(gòu)材料本構(gòu)參數(shù)采用初始值，與真實值存在較大差別，所以其結(jié)果也自然誤差最大。模型更新混合試驗方法，利用物理子結(jié)構(gòu)的試驗數(shù)據(jù)對數(shù)值子結(jié)構(gòu)的混凝土參數(shù)本構(gòu)進(jìn)行更新，很大程度上改善了模擬結(jié)果。但是模型更新混合試驗的邊界條件與傳統(tǒng)混合試驗一樣存在誤差：因為反彎點不一定在三分之二柱高處，亦即三分之二柱高處可能仍然存在彎矩；試件柱軸力恒定，也與實際受力不相符。而本文新方法的整體結(jié)構(gòu)數(shù)值模型無邊界條件問題，同時本構(gòu)參數(shù)也得到了更新，所以其精度高于另外兩種方法。

圖4 試驗裝置及試件

(a) 位移時程整體圖

(b) 位移時程局部放大圖

圖6給出了底層左柱2/3高處軸力和彎矩在各種方法下的對比結(jié)果，為了便于比較，這里的結(jié)果均指原型結(jié)構(gòu)(即相似比為1∶1)。結(jié)構(gòu)在動力荷載作用下軸力的變化幅值高達(dá)35%；彎矩且假定反彎點處的彎矩高達(dá)底部最大彎矩參考解的30%。由此可看出，簡化邊界條件將顯著改變結(jié)構(gòu)的受力狀態(tài)，影響結(jié)構(gòu)性能評估的準(zhǔn)確性。

表1 鋼筋和混凝土材料參數(shù)

圖7給出了模型更新混合試驗與本文所提方法混凝土峰值應(yīng)力和峰值應(yīng)變的識別結(jié)果，兩種情況下本構(gòu)參數(shù)在2.4 s內(nèi)歷經(jīng)一定波動后均收斂于真實值。

(a) 軸力時程(b) 彎矩時程

圖6 邊界軸力與彎矩時程對比

(12)

式中：xref，xsim分別為參考解和混合模擬結(jié)果。

4 結(jié) 論

為解決邊界條件不完整對混合試驗的不利影響，提高混合試驗的模擬精度，本文提出基于模型更新的

(a) fc(b) ε0

圖7 混凝土本構(gòu)模型部分參數(shù)識別結(jié)果

Fig.7 The identification results of constitutive parameters of concrete

在線數(shù)值模擬方法，稱之為考慮不完整邊界條件的新型混合試驗方法。與傳統(tǒng)混合試驗相比，本文所提方法物理子結(jié)構(gòu)的恢復(fù)力不再返回時間積分模塊參與運動方程的求解，而是返回到本構(gòu)參數(shù)識別模塊用于本構(gòu)模型在線識別。新方法中結(jié)構(gòu)反應(yīng)由整體結(jié)構(gòu)數(shù)值模型求得，邊界條件自然得到滿足；而且數(shù)值模型的本構(gòu)參數(shù)不斷根據(jù)物理試件的試驗結(jié)果在線修正，提高了數(shù)值計算的準(zhǔn)確性。對一邊界條件不完整的鋼筋混凝土框架結(jié)構(gòu)的混合試驗進(jìn)行了數(shù)值仿真分析，結(jié)果表明，新方法幾乎消除了邊界條件不足帶來的不利影響，極大地提高了混合試驗的模擬精度。