核電結構土-結相互作用分析分區混合計算方法1)

陳少林 郭琪超 周國良

?(南京航空航天大學土木與機場工程系，南京 210016)

?(環境保護部核與輻射安全中心，北京 100101)

引言

規范規定，核電屏蔽廠房的地震反應分析需要考慮土-結相互作用的影響[1].在核電結構的土-結動力相互作用分析中，阻尼是影響結構反應的一個重要因素.另外，出于安全性考慮，核電結構一般不允許進入非線性；而土體在地震作用下，容易進入非線性，因此土體非線性是影響土-結系統反應的另一重要因素.如何合理考慮阻尼和土體非線性是土-結動力相互作用分析的關鍵問題.其中，輻射阻尼一般通過人工邊界條件來考慮，如透射邊界[2]，黏彈性人工邊界[3-5]等.這里主要討論結構材料阻尼[6-17]和土體非線性.

土-結動力相互作用分析可采用頻域方法，如軟件SASSI[18].頻域方法可直接采用具有試驗觀測基礎的滯回阻尼模型，對于線性問題而言，可得到準確的結果，當考慮土體非線性時，需通過等效線性化進行處理.大多數研究表明，等效線性化適合于土體弱非線性，對于強非線性，宜采用時域逐步積分方法[19].若采用時域分析之模態疊加法進行土-結動力相互作用分析，可以直接采用規范規定的模態阻尼比，當體系的反應主要由低階模態控制時，具有運算速度快，阻尼輸入準確等優點，但原則上不適合于非線性分析.

當采用時域逐步積分方法時，土體可采用非線性黏彈性本構，能較為合理地反應其非線性特性.對于結構而言，時步積分法常采用的阻尼模型為瑞利阻尼模型和Caughy 阻尼模型.瑞利阻尼與質量和剛度成正比，通常稱為比例阻尼，通過兩階模態的頻率和阻尼比來確定兩個比例常數.當對地震反應起主要貢獻的結構模態數為兩個時，采用這兩個模態的頻率和阻尼比確定的瑞利阻尼可以較準確地反應這兩個模態的阻尼比，由瑞利阻尼模型計算的反應較為準確.但當對地震反應有貢獻的模態數較多時，瑞利阻尼能較準確反應指定阻尼的兩階模態阻尼，其余模態阻尼與真實情況有誤差，使得多數模態反應失真，造成地震反應與真實解有較大差異.鄒德高等[16]、李小軍等[17]分別改進了瑞利阻尼模型中兩系數的確定方法，并分別對土石壩和核電廠房進行了地震反應分析，但本質上未改變瑞利阻尼模型的實質，仍是通過同樣的函數來近似阻尼.Caughy 阻尼模型通過級數形式描述阻尼，可以使得更多階的模態阻尼比滿足規定值，但在通過滿足多個模態阻尼比確定級數的系數時，有時會出現系數矩陣為奇異的情形，造成系數求解的困難.另外，Caughy 阻尼在高頻時可能出現負阻尼情形，影響計算的失穩.Luco等[9]提出了一種Caughy 阻尼系數優化方法，可以避免上述相關問題.

綜上所述，核電結構不允許進入非線性，采用模態疊加法可簡便合理地記入模態阻尼.非基巖場地在地震作用下容易進入非線性，宜通過非線性黏彈性本構描述非基巖場地特性，采用時步積分法進行分析.目前土-結相互作用分析常用的軟件SASSI[18]，采用頻域分析方法，通過等效線性化方法考慮土體非線性，在強震時不能很好地體現土體特性，且不能考慮土與基礎間的接觸非線性[20-21].采用ANSYS，ABAQUS，OpenSees 等軟件，結構和土體只能采用相同的分析方法(要么都采用時步積分法，或都采用模態疊加法)，且計算效率較低[22](23 961 個節點，18 200 個單元，8000 時步數，在Intel Core i7 2.93 GHz，8 GB 內存的微機上用時4 周).因此，有必要發展一種能合理考慮結構阻尼和土體非線性的高效時域土-結相互作用分析方法.

本文在顯-隱式時域土-結相互作用分區算法的研究基礎上[23-24]，改用模態疊加法進行結構分析，土體仍采用集中質量顯式有限元方法結合人工邊界條件進行模擬，在每一時步，通過大質量法[25-27]進行多點激勵輸入，實現了模態疊加和時步積分結合的土-結相互作用分區混合算法.該方法能合理地考慮輻射阻尼和結構的材料阻尼，也可考慮土體的非線性[28].另外，結構和土體可采用不同的時間步距，且方便采用并行計算技術，具有較高的效率.通過簡單算例對該方法進行了驗證，并對某復雜場地上CAP1400 核電模型進行了土-結相互作用分析，對比分析了結構阻尼模型對核電結構反應的影響.

1 基本理論

圖1 為結構-基礎-土體模型示意圖，對該體系進行有限元離散，并將節點類型分為結構節點、結構與基礎的界面點、土體和基礎節點，以及人工邊界點.則體系的運動方程可寫為

式中，下標s,b,i 和a 分別表示結構節點、結構與基礎的界面節點、土體節點和人工邊界點.上標s 和g 分別表示結構和基礎.其中Kaa和Caa分別為黏彈性邊界的彈簧和阻尼系數矩陣，fa為地震波輸入時施加在人工邊界節點上的等效載荷[3-5].若采用透射邊界，可通過多次透射公式在人工邊界上施加位移[2].對方程(1)通過時步積分方法直接進行求解，即為土-結相互作用的直接法或整體解法.若采用隱式解法，則需每時步求解大型方程組，計算量很大，十分耗時.若采用集中質量顯式積分方法，每一時步不需求解大型方程組，但結構波速較大，穩定性要求時間步距較小，計算時步數較多，效率受影響.

圖1 土-基礎-結構整體分析模型示意圖Fig.1 Soil-foundation-structure model

若將式(1)變換到頻域，得到頻域形式的運動方程

其中，動力剛度

式中，U和u,F和f分別為傅里葉變換對.

消去方程(2)中土體節點的自由度，可得[30]

注意，這里的下標s 包含結構和基礎的節點，下標b 為基礎與土體的界面點，為基礎的動力剛度，為開挖掉的土體動力剛度，為土體與基礎界面點的自由場位移.求得自由場和動力剛度后，即可由式(4)求得結構和基礎的頻域響應，這種方法稱為子結構方法.子結構法在頻域內進行分析，原則上只適合于線性情形.

若將整個系統進行分區，分為上部結構、下部基礎和土體，按此分區，將方程(1)分開寫成如下形式

其中，式(5)的右端項為基礎對結構的作用力，式(6)的右端項第一分量為結構給基礎的作用力，兩者為一對作用力和反作用力.

考慮到土體自由度數目較大，采用集中質量顯式積分方法效率更高.因此，對式(6)采用集中質量形式，并采用顯式積分格式，如單邊中心差分格式

則式(6)中每一節點k的位移可通過如下方程求解

其中，N為與節點k相鄰的節點總數.?t為時間步距,分別為節點k在t=p?t時刻的加速度向量、速度向量和位移向量.mk為集中于節點k的質量,Ckj和Kkj分別為節點k與相鄰節點j之間的阻尼陣和剛度陣,為p時刻作用在節點k上的載荷向量.若k屬于基礎與結構相連的界面點，則為結構施加在基礎上的載荷；若k屬于人工邊界點，當采用黏彈性邊界時，為地震輸入時的等效載荷，當采用透射邊界時，該點的位移直接由多次透射公式求得；若k為基礎和土體的其余節點，則為零.

由式(7)～式(9)求得土體和基礎(p+1)時刻的反應后，則式(5)的右端項已知，可求得結構(p+1)時刻的反應，包括結構作用在基礎上(p+1)時刻的載荷.我們稱該方法為分區方法(partitioned method)，分區方法的優點是土體和結構可獨立建模，采用適合各自特點的分析方法，在每一時步獨立分析，且可采用不同的時間步距，便于獨立開發各自的分析程序，或應用已有的分析程序，因此具有較大的靈活性和較高的效率.分區方法的缺點是有可能失穩，其失穩機理和穩定性條件還有待進一步研究.我們已編制相應的并行計算程序，實現了土-結相互作用的分區分析，稱之為PASSI (partitioned analysis of soil-structure interaction)[23-24,28-29].在文獻[23-24]中，結構采用Newmark 積分方法，這里采用模態疊加法.

圖2 弱耦合示意圖Fig.2 Illustration of loose coupling

求得結構基底的反應后，即可得到式(5)右端的載荷項，可由模態疊加法求得結構的響應以及結構給基礎的反力.若結構模型較為復雜，則每一時步求得式(5)右端項較為麻煩.由式(7)～式(9)求得結構與基礎界面點的加速度等響應后，由于連續條件，將界面點的響應施加于結構底部，相當于是多點激勵下的結構分析，并將結構對基礎的作用力反饋給基礎，如圖2 所示，關于剛性基礎情形，見文獻[23-24].這里，我們考慮加速度連續，通過大質量法[27]將加速度施加于結構底部，因此，上部結構的運動方程如下

大質量法通過在大質量基礎點上施加力載荷模擬地震作用；在數學處理上比較巧妙地通過在質量矩陣上“置大數”實現近似于真實值的地震動輸入，因此中的元素一般取結構總質量的106倍.

結構采用模態疊加方法,時間步距的選取滿足精度要求即可，可較土體分析的時間步距大，即結構和土體可以采用不同的時間步距.已知p時步及以前各時步土、基礎和結構的位移，求解(p+1)時步各點的位移，土-結相互作用分析的基本步驟如下：

(1)根據式(7)～式(9)可以計算土體和基礎節點及人工邊界節點(p+1)時步的響應，得到結構底部節點的加速度響應；

(2)以結構底部節點(p+1)時步的加速度為多點激勵，通過大質量法得到施加在結構底部大質量點上的力(式(11))，對式(10)采用模態疊加法，計算結構響應，并得到(p+1)時步結構對基礎的反力；

(3)重復以上各步，即可得到土-基礎-結構體系各時刻的反應.

2 算例分析

根據上述原理，我們編制了相應的計算程序，實現了模態疊加和時步積分結合的三維土-結相互作用分析的分區并行計算方法.對于無限域土體和基礎的動力響應，采用自編的Fortran 程序進行分析.對于結構響應，其每一時步的計算獨立于土體的計算，因此可使用ANSYS 等商業軟件進行分析，結構和土體可分別采用不同的時間步距.通過耦合算法和Fortran 程序與ANSYS 之間的交互，實現土-結相互作用動力分析.由于采用分區計算方式，土體和結構可以獨立進行建模，且在每一時步，兩者獨立進行計算.土體采用MPI協議，編程實現并行.結構可采用ANSYS 中的并行計算方案.土體和結構之間的并行通過異步傳輸數據實現，具體見文獻[17].下面，通過一簡單模型的算例對該方法進行驗證，并對某非均勻場址上CAP1400 核電結構的反應進行分析.

2.1 簡單模型土-結相互作用分析

2.1.1 方法驗證

計算模型如圖3所示，上部結構尺寸為2 m×2 m×10 m，采用1 m×1 m×1 m 的六面體八節點實體單元進行離散.選取三層水平成層場地，土層材料參數及相應厚度見表1.選取土體計算區域的尺寸為40 m×40 m×18 m，將土體沿X方向劃分為3 個子區域：土體1，土體2 和土體3，采用3 個進程進行并行計算.土體離散為1 m×1 m×1 m 的六面體八節點實體單元，單元總數為39 600，節點總數為31 939.剛性埋置基礎尺寸為6 m×6 m×4 m.

圖3 計算模型示意圖Fig.3 Numerical model

表1 下臥土體參數Table 1 Soil parameters

選取脈沖寬度0.15 s 的單位脈沖波(如圖4 所示)作為SV 波，垂直入射.土體采用顯式中心差分格式，時間步距?t1=2.0×10?4，計算步數為8192 步，結構分別采用Newmark 隱式積分格式和模態疊加法進行分析，不考慮結構阻尼，用于驗證模態疊加與時步積分結合的土-結相互作用分區計算方法的有效性.圖5～圖7 分別為土體A點、結構B,C,D點及基礎的位移，點的位置如圖3 所示.圖5～圖7 中實線為結構采用時步積分法的結果，虛線為結構采用模態疊加法的結果.從圖5～圖7 可以看出，兩種結果完全重合，這驗證了模態疊加方法結合時步積分法進行土-結相互作用分析的有效性.

圖4 輸入單位脈沖Fig.4 Pulse input

圖5 場地A 點的位移反應Fig.5 Displacements of point A

圖6 結構點的位移反應Fig.6 Displacements of structure

圖6 結構點的位移反應(續)Fig.6 Displacements of structure(continued)

圖6 結構點的位移反應(續)Fig.6 Displacements of structure(continued)

圖7 基礎位移反應Fig.7 Displacements of foundation

圖7 基礎位移反應(續)Fig.7 Displacements of foundation(continued)

2.1.2 阻尼模型影響分析

結構的阻尼分別采用Rayleigh 阻尼和模態阻尼，計算模型及其余參數與3.1 中相同.模態阻尼比取0.05，選取X方向質量影響系數最大的兩階自振頻率(分別為0.15 Hz 和0.92 Hz)，其對應的瑞利阻尼比同樣取為0.05，對應的瑞利阻尼系數α=0.081，β=0.014.由此得到的阻尼如圖8 所示.

圖8 阻尼曲線(右邊為局部放大圖)Fig.8 Damping curve(right:zoom in detail)

圖9 為剛性基礎的位移時程，其中實線為時步積分方法的結果，采用的是瑞利阻尼，虛線為模態疊加方法的結果，采用的是模態阻尼，后面的圖例與此相同，不再說明.可以看出，結構阻尼模型對基礎反應影響很小，可能是由于結構較小，土-結相互作用較小的原因.圖10 為B點和D點(位置如圖3 所示)的X方向位移時程.由圖中可以看出，結構阻尼模型對土體反應影響較小，但對結構反應影響較大.圖11 為D點X方向的位移頻譜，以及相對于基礎的傳遞函數HX(f)=Ux(f)/UFX(f)，其中UFX(f)為基礎X方向的位移頻譜.由圖中結果可看出，頻率大于2 Hz時，瑞利阻尼的結果要遠小于模態阻尼的結果，這與圖8 相一致，即大于2 Hz，瑞利阻尼遠大于模態阻尼.

圖9 基礎位移反應Fig.9 Displacements of foundation

圖10 土體B 點與結構D 點位移反應Fig.10 Displacements of point B and D

圖11 結構D 點的位移頻譜和傳遞函數Fig.11 Spectrum and transfer function of point D

2.2 核電結構土-結相互作用分析

2.2.1 場地模型

根據地脈動測試和鉆孔資料，獲得了某核電場地的剪切波速剖面，如圖12 所示.土體參數如表2 所示.選取土體計算區域的尺寸為640 m ×360 m×194 m，邊界采用黏彈性邊界.土體離散為2 m×2 m×2 m 的六面體八節點實體單元，單元總數為 5 587 200，節點總數為5 693 898.采用集中質量顯式有限元方法進行分析，時間步距?t1=1.0×10?4s.

2.2.2 結構模型

核電結構模型如圖13，結構單元數為597 686，相應的節點數為700 194.剛性基礎尺寸為92 m×60 m×16 m.結構分別采用Newmark 隱式時步積分方法和模態疊加方法進行分析，對應的阻尼分別采用瑞利阻尼和模態阻尼.模態阻尼比取為0.05，選取X方向質量影響系數最大的兩階自振頻率(分別為3.75 Hz 和6.88 Hz)，其對應的瑞利阻尼比同樣取為0.05，對應的瑞利阻尼系數α=1.509，β=0.001 5，由此得到的阻尼如圖14 所示.采用模態疊加法，選取了六階剛體模態以及300階非剛體模態.結構分析的時間步距?t2=25?t1，即?t2=2.5×10?3，本算例使 用DELL-Optiplex 小型工作站進行計算，CPU 是Intel@Xeon(R)CPU@E5-2667v4@3.2 GHz×32，主存大小為64 G，操作系統為Ubuntu16.04LTS.土體采用5 進程，結構采用2 進程，并行計算.

圖12 場地剪切波速剖面圖及模型Fig.12 Shear velocity profil and soil model

表2 土體參數Table 2 Soil parameters

圖13 計算模型示意圖Fig.13 Numerical model

圖14 阻尼曲線Fig.14 Damping curve

2.2.3 脈沖波輸入情形

考慮SV 波垂直入射，輸入的脈沖位移時程圖和頻譜圖如圖15 所示.

圖15 脈沖波輸入Fig.15 Pulse input

圖16 場地X 方向位移時程及頻譜Fig.16 X-direction displacements on site

圖16 場地X 方向位移時程及頻譜(續)Fig.16 X-direction displacements on site(continued)

圖16 為C,D,E三點(位置見圖14 所示)的位移時程及頻譜圖.對比圖中3 點的位移時程及頻譜圖，可以看出，離基礎越遠，結構阻尼模型對土體反應的影響越小.從C點位移頻譜圖可知，小于3.75 Hz 時，模態阻尼的結果要大于瑞利阻尼的結果，這與圖14 的阻尼曲線一致.總體而言，結構阻尼模型對土體反應影響較小.圖17 為基礎的位移時程及頻譜圖，可以看出，結果阻尼模型對基礎反應有較明顯的影響.由基礎X方向的位移頻譜圖可以看出，在結構自振頻率處(如3.75 Hz 和6.88 Hz)，基礎的穩態反應接近零，符合土-結相互作用理論.結合圖14，在3.75 Hz 至10 Hz，瑞利阻尼與模態阻尼較為接近，所以兩種阻尼模型對應的基礎X方向穩態反應差別不大，但小于3.75 Hz 時，瑞利阻尼遠大于模態阻尼，因此在2 Hz 左右，其反應明顯小于模態阻尼.

圖17 基礎位移時程及頻譜Fig.17 Displacements of foundation

圖17 基礎位移時程及頻譜(續)Fig.17 Displacements of foundation(continued)

圖18 為結構11 m 高度處的截面圖，給出其中2107 號節點的反應.圖19 為核島結構屏蔽廠房剖面及參考點位置圖，給出其中136 340，136 367，185 547，184 783，64 139 號節點的反應.由于SV 波垂直入射，主要產生X方向的位移，因此只給出上述節點X方向的位移及其頻譜圖，見圖20.

圖18 核島11 m高度截面圖及參考點位置Fig.18 11 m height section of nuclear island and reference point

圖19 核島結構屏蔽廠房剖面及參考點位置Fig.19 Section of nuclear island and reference points

圖20 結構點的位移時程及其頻譜Fig.20 Displacements of structure

圖20 結構點的位移時程及其頻譜(續)Fig.20 Displacements of structure(continued)

圖20 結構點的位移時程及其頻譜(續)Fig.20 Displacements of structure(continued)

結合圖15 與圖20 進行分析，總體而言，在3.75 Hz 至6.88 Hz，瑞利阻尼與模態阻尼較為接近，所以兩種阻尼模型對應的穩態反應差別不大，但頻率小于3.75 Hz 或大于6.88 Hz 時，瑞利阻尼遠大于模態阻尼，因此其反應明顯小于模態阻尼的反應，這種差異在結構中上部點越為明顯.對于136 367 點的反應，在5.5 Hz 左右，模態疊加反應明顯大于瑞利阻尼的結果，可能與局部模態有關.

2.2.4 地震波輸入情形

采用地震安全性評價得到的人工地震波.圖21 是露頭基巖處的X方向的加速度時程及其傅里葉幅值譜.這里假定為SV波垂直入射，在邊界區近似為水平成層場地，將露頭基巖處的地震波折減一半做為輸入，采用傳遞矩陣方法，計算得到邊界自由場，進而得到黏彈性邊界的等效載荷，做為土-結相互作用分析的輸入.

圖21 人工地震波加速度時程及其頻譜Fig.21 Artificia seismic wave

圖22 所示為C,D,E三點的位移、加速度和加速度反應譜(阻尼比均為5%).從圖22 可以看出，地震波輸入時，結構阻尼模型對土體的位移響應基本沒有影響.對離基礎較近(C點)的土體加速度存在影響，但對離基礎較遠的D和E點，則影響很小.

圖22 場地X 方向反應時程及加速度反應譜Fig.22 X-direction displacements on site

圖22 場地X 方向反應時程及加速度反應譜(續)Fig.22 X-direction displacements on site(continued)

圖23 為基礎的位移時程，由于SV 波垂直入射，主要產生X方向的響應，其余方向的響應較小.從圖中可以看出，阻尼模型對基礎位移的影響較小.

圖23 基礎位移時程Fig.23 Displacements of foundation

圖23 基礎位移時程(續)Fig.23 Displacements of foundation(continued)

圖24 為結構點的響應(點號及位置見圖18 和圖19)，從圖中可以看出，結構阻尼模型對結構的響應影響較大，采用瑞利阻尼時結構的響應要比采用模態阻尼的小.從加速度反應譜可以看出，瑞利阻尼的結果要小于模態阻尼的結果，越靠近結構頂部越明顯(如圖中185 547 號點).這與圖15 所示的核電結構各振型阻尼比一致，即盡在3.75 Hz 至6.88 Hz 的很小頻段范圍內，瑞利阻尼與模態阻尼較為接近，在其余頻率，瑞利阻尼遠大于模態阻尼，因此其反應明顯小于模態阻尼的反應.

圖24 結構X 方向位移、加速度及反應譜Fig.24 X-direction responses of structure

圖24 結構X 方向位移、加速度及反應譜(續)Fig.24 X-direction responses of structure(continued)

圖24 結構X 方向位移、加速度及反應譜(續)Fig.24 X-direction responses of structure(continued)

圖24 結構X 方向位移、加速度及反應譜(續)Fig.24 X-direction responses of structure(continued)

3 結論

本文提出并實現了模態疊加和時步積分相結合的土-結動力相互作用分析的分區算法，可較為合理地考慮核電結構土-結相互作用分析中的阻尼.該方法具有如下特點：(1)核電結構采用模態疊加法進行分析，可采用實測的模態阻尼；(2)土體采用顯式時步積分分析，可考慮土體的非線性和滯回阻尼；(3)半無限土體的輻射阻尼可通過黏彈性邊界或透射邊界計入；(4)結構和土體可分別采用滿足各自精度和穩定性要求的時間步距，且便于并行計算，效率較高.

通過算例驗證了該分區算法的有效性.對某復雜場地上CAP1400 核電模型的分析結果表明，結構分別采用瑞利阻尼和模態阻尼，對場地反應的影響不大，但對結構反應的影響較為明顯.因此，在核電結構的土-結相互作用分析中，需要合理選擇阻尼模型.