求解動態子結構的大質量界面虛擬激勵法

孫萬泉，秦浩瑋，馬震岳

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

引言

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

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

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

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

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

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

將上式按第二項展開有：

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

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

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

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

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

式（6）整理后可得：

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

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

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

由此可得，在結構邊界處添加大質量M后，對邊界施加M?Peiωt的加速度激勵，即可得到邊界處加速度為Peiωt的結構內部動力方程。

2 大質量界面虛擬激勵法的提出

以最簡單的兩個子結構為例進行該方法的說明，具體實現步驟如下：

第一步，將結構劃分為子結構1 和2，其中結構外荷載位于子結構1 的內部節點上，記為。此時對于子結構1 和2 分別有：

式中 下標“1”和“2”為子結構的編號；角標“m”和“s”分別為子結構的內部自由度和界面自由度；和分別為子結構1 和2 的交界面上的界面力。

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

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

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

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

第六步，在求出能夠保證兩個子結構的界面加速度和界面內力都協調平衡的虛擬激勵F*(ω)后，就可以根據公式（15）和（16）求出各自子結構的真實內部響應。

第七步，采用以上方法對全頻段進行遍歷求解，即可得到結構在隨機激勵下的全部頻域響應，根據頻域解再應用Fourier 逆變換，就可得到時域解。

同理，如果系統包含n個子結構交界面，利用以上方法可以聯立方程組求解n個未知數(ω)。由于篇幅所限，不再贅述。

3 工程算例

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

3.1 集中質量模型

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

圖1 集中質量模型Fig.1 Lumped-mass model

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

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

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

3.2 地震和外荷載同時作用下的橋梁結構分析

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

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

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

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

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

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

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

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

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

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

3.3 地震作用下的大型復雜結構分析

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

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

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

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

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

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

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

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

4 結論

在直接求解的虛擬激勵法的原理上提出了一種求解動態子結構的大質量界面虛擬激勵法，通過理論推導與計算分析，可得到以下結論：

借助大質量點可以保證相鄰子結構的界面節點在相同的虛擬激勵下擁有相同的加速度。在保證界面加速度相同，界面內力數值相等，且方向相反的情況下，可以推導出保證各個子結構界面運動協調的虛擬激勵荷載。從而求解各個子結構的運動方程，實現大型結構的簡化計算。

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