孫萬泉，秦浩瑋，馬震岳

（1.華北電力大學水利與水電工程學院，北京 102206；2.大連理工大學建設(shè)工程學部水利工程學院，遼寧 大連 116023）

引言

隨機振動中的虛擬激勵法是一種既能保證計算精度，又能提高計算效率的分析方法［1］。該方法得到了國內(nèi)外學者的廣泛應(yīng)用［2］，其中，李永華等［3］提出了利用大質(zhì)量來直接求解絕對位移的虛擬激勵法；Zhao 等［4］從多點輸入的動力學方程入手提出了一種直接求解的多點虛擬激勵法；Jia 等［5］通過在多維質(zhì)量矩陣上“置大數(shù)”的方法，從理論上將一維多點激勵拓展到多維多點激勵。上述研究克服了傳統(tǒng)虛擬激勵法在有限元軟件中求解困難的缺點，解決了傳統(tǒng)虛擬激勵法需要將絕對位移進行分解后再求解的弊端，在多點隨機振動中已經(jīng)有了廣泛的應(yīng)用。

在實際工程中，很多大型復(fù)雜結(jié)構(gòu)或者不同類型的組合結(jié)構(gòu)的動力計算有時會采用動態(tài)子結(jié)構(gòu)方法。目前動態(tài)子結(jié)構(gòu)的計算方法主要有模態(tài)綜合類和機械導(dǎo)納類［6］。其中，模態(tài)綜合類可分為自由界面模態(tài)綜合法、固定界面模態(tài)綜合法和混合界面模態(tài)綜合法三類［7］。上述幾種模態(tài)綜合方法都是建立在先分解后復(fù)原的思想上，對子結(jié)構(gòu)進行高階模態(tài)的省略，再合并建立整體模態(tài)方程求解。另外，以模態(tài)綜合法為基礎(chǔ)衍生出的新型動態(tài)子結(jié)構(gòu)方法在不同的領(lǐng)域也有較好的應(yīng)用，例如，基于波分析的動態(tài)子結(jié)構(gòu)方法［8］和重復(fù)超單元法［9］等。這類方法都需要對模態(tài)坐標和物理坐標進行相互轉(zhuǎn)化，較為繁瑣。而機械導(dǎo)納類是通過建立以子結(jié)構(gòu)間的節(jié)點位移和節(jié)點力為未知數(shù)的綜合方程進行聯(lián)立求解。但當界面的節(jié)點數(shù)較多時，該方法的未知數(shù)較多，方程規(guī)模較大，求解困難。此外，當同時考慮結(jié)構(gòu)外荷載激勵和地震下基礎(chǔ)加速度激勵時，機械導(dǎo)納法的理論推導(dǎo)和求解也比較繁瑣。

為進一步提高動態(tài)子結(jié)構(gòu)方法的求解效率和便捷性，本文在絕對位移直接求解的虛擬激勵法［3，5］的啟發(fā)下，提出一種求解動態(tài)子結(jié)構(gòu)的大質(zhì)量界面虛擬激勵法。通過對各個子結(jié)構(gòu)界面上添加大質(zhì)量點，并在大質(zhì)量點上施加虛擬激勵荷載，在滿足子結(jié)構(gòu)間加速度相等和內(nèi)力平衡的條件下建立界面平衡方程組，求解能夠保持所有子結(jié)構(gòu)界面協(xié)調(diào)平衡的虛擬激振力幅值，并最終求解出各個子結(jié)構(gòu)的響應(yīng)。該方法與傳統(tǒng)的模態(tài)綜合法相比，避免了模態(tài)坐標和物理坐標的反復(fù)轉(zhuǎn)換，并能直接考慮子結(jié)構(gòu)的高頻特性；與機械導(dǎo)納法相比，綜合方程里減少了系統(tǒng)響應(yīng)的未知數(shù)，從而減小了求解方程的規(guī)模，使求解更為簡便。

1 絕對位移直接求解的虛擬激勵法

對于一個有N個支座和n個自由度的離散結(jié)構(gòu)，其多點地震激勵運動方程可寫成如下的分塊矩陣形式［1］：

式中m維列向量Xb，和分別代表N個 支座的地面強迫位移、速度和加速度；n維列向量Xs，和分別代表結(jié)構(gòu)系統(tǒng)所有非支座節(jié)點位移、速度和加速度；m維列向量Pb代表地面作用于N個支座的力；M，C和K分別為結(jié)構(gòu)的質(zhì)量、阻尼和剛度矩陣；下標“s”和“b”分別對應(yīng)結(jié)構(gòu)的非支座自由度和支座自由度。

將上式按第二項展開有：

因連接地面處的支座地震力Pb等于支座質(zhì)量Mb與地面加速度的乘積，即Pb=Mb，此時有：

在式（3）左右兩邊同乘Mb-1可得：

由此可得，只要支撐處的質(zhì)量矩陣足夠大，即在結(jié)構(gòu)總體質(zhì)量矩陣中支撐質(zhì)量矩陣處置入一個足夠大的數(shù)，就能保證支撐處的地面加速度與支撐地震響應(yīng)加速度相等［10］，從而使得式（1）的求解容易很多。

將式（5）代入式（1）第一項后展開可得：

式中S(iω)為加速度功率譜密度；ω為角頻率；上標“T”和“*”分別為矩陣的轉(zhuǎn)置和伴隨；e 為自然對數(shù)的底數(shù)；i 為虛數(shù)單位；SX為忽略局部場地效應(yīng)下的地面節(jié)點位移；t為加速度功率譜分解后得到的時間分量；P為n×r（r≤n）實陣；Q為與行波效應(yīng)相關(guān)的n×r實陣。

式（6）整理后可得：

將式（9）～（11）代入式（12）整理可得：

根據(jù)虛擬激勵法理論［1］假定阻尼力不是與絕對速度而是與相對速度成正比，即將公式（12）中的表示為動態(tài)相對速度，并用0 矩陣代替此時相當于將公式（13）中等效荷載項中的阻尼項忽略，因此可以得到虛擬激勵動力方程：

式（14）可以看作在每一個頻率點ω處的簡諧運動方程。

由此可得，在結(jié)構(gòu)邊界處添加大質(zhì)量M后，對邊界施加M?Peiωt的加速度激勵，即可得到邊界處加速度為Peiωt的結(jié)構(gòu)內(nèi)部動力方程。

2 大質(zhì)量界面虛擬激勵法的提出

以最簡單的兩個子結(jié)構(gòu)為例進行該方法的說明，具體實現(xiàn)步驟如下：

第一步，將結(jié)構(gòu)劃分為子結(jié)構(gòu)1 和2，其中結(jié)構(gòu)外荷載位于子結(jié)構(gòu)1 的內(nèi)部節(jié)點上，記為。此時對于子結(jié)構(gòu)1 和2 分別有：

式中 下標“1”和“2”為子結(jié)構(gòu)的編號；角標“m”和“s”分別為子結(jié)構(gòu)的內(nèi)部自由度和界面自由度；和分別為子結(jié)構(gòu)1 和2 的交界面上的界面力。

第二步，分別在子結(jié)構(gòu)1 和2 的交界面節(jié)點處添加大質(zhì)量單元M0（結(jié)構(gòu)整體質(zhì)量的105～108倍）。并對外荷載做Fourier 變換，得到外荷載在不同頻率ω下的Fourier 譜(ω)。

第五步，為保證兩個子結(jié)構(gòu)之間的界面加速度和界面內(nèi)力同時協(xié)調(diào)平衡，也就是要求在哪個虛擬激勵幅值（F*(ω)）下，可以使得兩個子結(jié)構(gòu)的界面加速度都等于，界面內(nèi)力絕對值也都等于N*（但方向相反）。因此，可以根據(jù)線彈性關(guān)系聯(lián)立方程如下：

以上方程類似兩條直線求交點，比較容易理解。整理方程（17）后可求得F*(ω)如下式所示：

可見，當系統(tǒng)只有一個交界面一個節(jié)點時，該方法也只有一個待求未知數(shù)F*(ω)。

第六步，在求出能夠保證兩個子結(jié)構(gòu)的界面加速度和界面內(nèi)力都協(xié)調(diào)平衡的虛擬激勵F*(ω)后，就可以根據(jù)公式（15）和（16）求出各自子結(jié)構(gòu)的真實內(nèi)部響應(yīng)。

第七步，采用以上方法對全頻段進行遍歷求解，即可得到結(jié)構(gòu)在隨機激勵下的全部頻域響應(yīng)，根據(jù)頻域解再應(yīng)用Fourier 逆變換，就可得到時域解。

同理，如果系統(tǒng)包含n個子結(jié)構(gòu)交界面，利用以上方法可以聯(lián)立方程組求解n個未知數(shù)(ω)。由于篇幅所限，不再贅述。

3 工程算例

首先以相對簡單的多自由度集中質(zhì)量模型進行分析，驗證這一求解方法。然后再以一個相對復(fù)雜的地震和外荷載同時輸入的橋梁模型，以及地震作用下的大型復(fù)雜水電站廠房結(jié)構(gòu)模型進行進一步計算驗證。其中，從工程角度，子結(jié)構(gòu)的劃分應(yīng)遵循以下原則［11］：按照實際復(fù)雜結(jié)構(gòu)幾何形狀和裝配關(guān)系劃分；盡量割斷較少的聯(lián)系，即盡量用較少的修改就能取得化整為零的最大效果等等。

3.1 集中質(zhì)量模型

為了驗證本文提出的界面虛擬激勵法求解動態(tài)子結(jié)構(gòu)的準確性和計算精度，先以六個自由度的集中質(zhì)量桿件系統(tǒng)在水平簡諧激勵作用下的響應(yīng)計算為例進行分析。計算模型和子結(jié)構(gòu)劃分如圖1 所示，水平簡諧激勵F=sin(10π?t)kN，作用于第一層質(zhì)量單元。各單元間距為2 m，各單元質(zhì)量為M=10 kg。

圖1 集中質(zhì)量模型Fig.1 Lumped-mass model

利用動態(tài)子結(jié)構(gòu)模型的計算結(jié)果與整體結(jié)構(gòu)計算結(jié)果進行相互對比驗證，記結(jié)果誤差值P=×100%，其中，本文的子結(jié)構(gòu)求解方法的幅值解為S1，整體模型求解的幅值結(jié)果為S2。計算結(jié)果對比如表1 所示。

表1 利用界面虛擬激勵法的子結(jié)構(gòu)求解與整體結(jié)構(gòu)求解結(jié)果對比Tab.1 Comparison of the substructure calculating results and the global structure calculating results using the interface virtual excitation method

由表1 結(jié)果可知，利用本文方法的計算結(jié)果與整體直接計算結(jié)果幾乎相等，誤差最大為1.06%，因此可以有效證明本文方法的準確性和精度。

3.2 地震和外荷載同時作用下的橋梁結(jié)構(gòu)分析

以某混凝土橋梁為例，全橋長為360 m，跨徑布置 為110 m+160 m+100 m，寬 為13 m，厚度為3.5 m，兩墩分別為高65 和105 m、寬3.5m的實心矩形墩。該橋三維模型和子結(jié)構(gòu)劃分如圖2 所示，模型采用六面體實體網(wǎng)格剖分，共計376 個單元和805 個節(jié)點。子結(jié)構(gòu)界面上的節(jié)點全部施加大質(zhì)量單元和虛擬激勵求解。墩底處完全約束，橋臺處約束豎向位移。橋墩處受到垂直于橋向的地震加速度荷載，同時，橋面中心受到豎直向下，均方根值為10 kN 的白噪聲隨機激勵的外荷載F。

圖2 三維有限元模型Fig.2 Three-dimensional finite element model

地震加速度模型選用胡聿賢和周錫元的修正Kanai-Tajimi 加速度隨機模型：

式中S0為譜強度因子；ωg和ξg分別為場地特征頻率和阻尼比；ωc為控制低頻含量參數(shù)。根據(jù)文獻［12］進行參數(shù)取值，ωg=13.96 rad/s，ξg=0.8，ωc=0.6π，S0=8.6697 cm2/s3。

根據(jù)結(jié)構(gòu)的自振頻率，取頻率積分區(qū)間為ω∈[0.5，60] rad/s，步長為0.5rad/s。地震虛擬力荷載為地震加速度與大質(zhì)量的乘積，其中地震加速度是由地震譜轉(zhuǎn)換而來的，如公式（7）～（9）所示。為更清晰地驗證本文方法在不同頻段下都具有足夠的準確性，取本文動態(tài)子結(jié)構(gòu)方法的計算結(jié)果與整體結(jié)構(gòu)直接求解的結(jié)果進行比較。選取子結(jié)構(gòu)1 上的節(jié)點A，子結(jié)構(gòu)1 和2 連接界面上的節(jié)點B，以及子結(jié)構(gòu)3 上的節(jié)點C 進行頻域結(jié)果比較，并單獨列出高頻區(qū)的結(jié)果進行對比。在地震加速度與橋面外荷載共同激勵下的計算結(jié)果如圖3～5 所示。

圖3 節(jié)點A 的豎向位移幅值Fig.3 Vertical displacement amplitudes of node A

圖4 節(jié)點B 的豎向位移幅值Fig.4 Vertical displacement amplitudes of node B

圖5 節(jié)點C 的豎向位移幅值Fig.5 Vertical displacement amplitudes of node C

從圖3～5 中可以看出，求解動態(tài)子結(jié)構(gòu)的界面虛擬激勵法在全頻段內(nèi)，無論是子結(jié)構(gòu)界面還是子結(jié)構(gòu)內(nèi)部，在求解隨機激勵問題時與整體結(jié)構(gòu)直接求解結(jié)果都非常吻合，該方法的子結(jié)構(gòu)計算結(jié)果與整體求解的誤差在0.6%左右。極少數(shù)點的誤差雖能達到6%，但是從計算結(jié)果來看已經(jīng)達到10-6量級，對整體響應(yīng)結(jié)果的影響很小。

動力子結(jié)構(gòu)方法能夠大幅度地縮減動力分析的規(guī)模，實現(xiàn)大規(guī)模復(fù)雜模型的分塊并行求解，避免了整體高階動力學方程的求解，提高了計算效率。計算效率的定量分析與具體動力模型的規(guī)模和子結(jié)構(gòu)的劃分有關(guān)。本算例中，整體模型計算時間為9.3 s，運用子結(jié)構(gòu)計算時間為4 s，計算效率提高了57%。隨著模型規(guī)模和計算步數(shù)的增加，該子結(jié)構(gòu)法的計算優(yōu)勢也更加突出。

3.3 地震作用下的大型復(fù)雜結(jié)構(gòu)分析

水電站廠房是由大體積不規(guī)則混凝土結(jié)構(gòu)和上部框架板肋結(jié)構(gòu)組成的大型復(fù)雜結(jié)構(gòu)體系。以某水電站廠房為例，根據(jù)其結(jié)構(gòu)形式特點，劃分為上部和下部兩個子結(jié)構(gòu)，如圖6 所示。廠房總長為65 m、寬為35 m、高 為48.5 m，模型共 計83050 個單元 和53108個節(jié)點。邊界條件為基礎(chǔ)底面與尾水管上游側(cè)完全約束，采用與3.2節(jié)相同的加速度模型在基礎(chǔ)底面施加與廠房橫河向呈45°角的水平地震加速度荷載。

圖6 三維有限元模型Fig.6 Three-dimensional finite element model

根據(jù)結(jié)構(gòu)的自振頻率，取頻率區(qū)間為ω∈[0，35] Hz，取步長為1 Hz。為驗證本文方法在大型復(fù)雜結(jié)構(gòu)下也具有足夠的先進性與準確性，取本文動態(tài)子結(jié)構(gòu)方法的計算結(jié)果與整體結(jié)構(gòu)直接求解的結(jié)果進行比較。選取子結(jié)構(gòu)1 上的節(jié)點a 和b，子結(jié)構(gòu)2 上的節(jié)點c 和d 進行橫向與縱向的頻域結(jié)果比較，如圖7～10 所示。

圖7 節(jié)點a 在0～35 Hz 頻段的位移幅值Fig.7 Displacement amplitude of node a in 0～35 Hz frequency band

圖8 節(jié)點b 在0～35 Hz 頻段的位移幅值Fig.8 Displacement amplitude of node b in 0～35 Hz frequency band

圖9 節(jié)點c 在0～35 Hz 頻段的位移幅值Fig.9 Displacement amplitude of node c in 0～35 Hz frequency band

圖10 節(jié)點d 在0～35 Hz 頻段的位移幅值Fig.10 Displacement amplitude of node d in 0～35 Hz frequency band

從圖7～10 中可以看出，求解動態(tài)子結(jié)構(gòu)的大質(zhì)量界面虛擬激勵法在整體模型較大、較復(fù)雜時也具有足夠的精度。其中，各個點位的求解誤差均小于5%。該大型復(fù)雜結(jié)構(gòu)在整個頻率區(qū)間的整體模型計算時間為560 s，分為兩個子結(jié)構(gòu)的情況下計算時間為455 s，計算效率提高了19%。

4 結(jié)論

在直接求解的虛擬激勵法的原理上提出了一種求解動態(tài)子結(jié)構(gòu)的大質(zhì)量界面虛擬激勵法，通過理論推導(dǎo)與計算分析，可得到以下結(jié)論：

借助大質(zhì)量點可以保證相鄰子結(jié)構(gòu)的界面節(jié)點在相同的虛擬激勵下?lián)碛邢嗤募铀俣?。在保證界面加速度相同，界面內(nèi)力數(shù)值相等，且方向相反的情況下，可以推導(dǎo)出保證各個子結(jié)構(gòu)界面運動協(xié)調(diào)的虛擬激勵荷載。從而求解各個子結(jié)構(gòu)的運動方程，實現(xiàn)大型結(jié)構(gòu)的簡化計算。

本文提出的方法在界面單節(jié)點與界面多節(jié)點以及一維與多維的情況下都具有較高的準確性，并且大幅減小了模型整體計算求解的規(guī)模，實現(xiàn)了各個子結(jié)構(gòu)模型單獨求解，提高了模型的計算效率。相較于傳統(tǒng)動態(tài)子結(jié)構(gòu)方法，在處理混合荷載作用下的結(jié)構(gòu)時也能保證足夠的精度和簡便。