基于ANSYS平臺的非巖性地基條件下核島廠房結構三維地震響應分析

尹訓強 趙文燕 王桂萱

（大連大學土木工程技術研究與開發中心，大連116622）

0 引 言

核電是高效型、安全型能源，然而由于國內不同地區地質條件差異較大，以及近些年來我國核電的飛速發展導致核電廠址選址問題日益突出，因此，內陸非巖性廠址已成為目前核電站建設的必然選擇［1-2］。其中，無限地基和近場地基非線性特征對核電結構地震響應的影響成為核電廠選址的關鍵技術問題。

考慮土-結構相互作用效應（SSI）的方法主要有比例邊界有限元法（SBFEM）［3］、阻尼抽取法（DSEM）［4］等子結構法以及人工邊界法［5］等直接法。與其他方法相比，人工邊界法通過在半無限地基的外邊界設置反射與折射等能量傳播方式，能夠有效吸收反射波動能量，同時可方便考慮樁土相互作用和近場地基非線性的影響，其中，由于黏性邊界［6］易于操作，概念清晰，在工程中應用較廣泛，同時也是核電抗震規范推薦方法。

近年來，非巖性地基［7］場地條件下核島廠房地震響應分析也備受國內外學者關注，并取得了很多有價值的成果。王桂萱，盛超等［8］建立了考慮近場非線性和土-結構相互作用的核島廠房三維地震響應分析模型，研究非巖性地基對結構響應的影響。鄒德高，隋翊等［9］在考慮地基土的非線性動力特性，采用人工邊界和等效結點荷載的地震動輸入方法的基礎上，考慮核島結構與深厚覆蓋土層地基的相互作用，研究結構動力響應的變化規律。尹訓強，薛志強等［10］采用等效線性和黏性人工邊界模擬近場地基特征，對土質地基條件下的核島廠房結構進行了擬三維地震分析；然而，擬三維分析需在兩垂直平面內簡化結構模型，也不能實現三方向上地震波的輸入，不能真實反映土體運動對變形特性參數的影響。

針對以上問題，本文基于等價線性法和黏性邊界的相關理論，開發了三維等價線性單元和黏性人工邊界單元，并結合UPFs二次開發特點嵌入至通用有限元軟件ANSYS中，建立了非巖性地基-核島廠房結構的動力相互作用分析模型，進而通過數值算例及實際工程對該計算模型進行合理性及實用性的驗證分析，并探討無限地基與近場非線性特征對核電工程結構地震響應的影響。

1 非巖性地基條件下核島廠房結構地震響應分析模型

如圖1所示，核島廠房結構-非巖性地基相互作用分析模型由近場地基與模擬無限地基的黏性人工邊界模型所構成。如前所述，近場非巖性地基的非線性動力特性和不均勻性，以及無限地基輻射阻尼是需要考慮的影響因素。

圖1 非巖性地基條件下的相互作用系統分析模型Fig.1 Structure-foundation interaction model of nuclear island factory building under the condition of non-lithology foundation

1.1 近場地基非線性模擬

在強震作用下，非巖性地基條件下土體的非線性動力特征顯著，直接非線性時域方法由于本構模型的復雜性以及數值計算的穩定性在核電工程中并未得到廣泛應用。美國土木工程師協會（ASCE）建議的等價線性法［11］概念明確且優勢突出，是不可替代的主流方法。

等價線性法［12-13］通過線性迭代求解非線性問題，關鍵思想是若動力分析迭代前后單元應變變化顯著，則利用每個單元的等效剪應變γeff=0.65γmax，參照所賦予材料的D-γ與G-γ關系曲線重新得出修正剪切模量G和阻尼比D，進行下一次動力分析，直到應變變化滿足收斂要求。具體計算過程如圖2所示。

圖2 等價線性解析流程圖Fig.2 Equivalent linear flow chart

1.2 黏性人工邊界模型及地震動輸入

由于近場地基由等價線性法來描述，每次迭代的動力分析中，土體則認為是線性材料。因此，無限地基輻射阻尼的影響可由黏性人工邊界來模擬。

黏性人工邊界模型［14］的核心是在地基的邊界處布置阻尼器吸收邊界處的反射波波動能量（圖1），并通過施加節點等效荷載模擬截斷外邊界處實際波場的應力邊界條件。在三維計算模型中，輸入地震動轉化為等效荷載的具體計算公式為

式中：F為施加在人工邊界節點處的等效荷載；下標n表示當前邊界面的外法線方向；下標tj表示對應的切線方向；Ct與Cn分別為黏性人工邊界節點處的切向阻尼系數和法向阻尼系數；u?f為自由波場條件下邊界處節點的速度；σf為自由條件下邊界處的正應力；τf為相應的切應力。

對于非均質近場地基的地震動輸入，可先進行自由場分析得到外邊界處的節點響應，而后利用式（1）進行各節點處的等效荷載計算，在具體數值實現中，可建立土柱有限元模型進行自由場分析［15］。

2 基于ANSYS平臺的地震響應分析模型開發

2.1 三維等價線性單元

傳統的有限元軟件中，數值迭代往往是在某一荷載步中完成，無法實現線性動力分析的循環迭代。為了在ANSYS中基于等價線性法描述非巖性地基的動力特性，必須使地基單元能夠描述D-γ與G-γ關系曲線且能實現線性迭代過程。本文提出一種地基側三維等價線性單元，該單元形式是對已有等參元材料模型的修改。

在用戶可編程特性（UPFs）中，用戶可修改的UserElem.F是創建用戶單元的接口子程序。通過該接口程序，可從標準ANSYS程序中獲取原地基單元的源數據，包括節點號、節點坐標、彈性模量、密度、積分點等，該子程序接口中可同時創建多個單元，不同的單元類型由KEYOPT關鍵選項實現轉換。等價線性單元與標準單元的不同關鍵在于每次線性動力迭代分析后的剛度陣和阻尼陣的變化，該處需要通過UPFs提供的自定義函數來實現阻尼比和剪切模量更新后實時數據的傳遞，之后便可計算單元的剛度陣以及阻尼陣等輸出數據返回給標準程序，保證下一次動力分析時為修正的剛度陣及阻尼陣。該用戶單元子程序的編制流程如圖3所示，程序的正文部分與標準等參元的編制大致相同，僅為剪切模量G和阻尼比D的差異及誤差分析的添加，需要注意的是，第2步中需要接口子程序User01.f實現實時數據的更新。

圖3 等價線性用戶單元子程序的編制流程Fig.3 The compiling sketch of interface subroutine

2.2 黏性人工邊界單元

以往的文獻中，人工邊界的開發往往是將物理元件參數乘以每個節點的控制面積，然后在每個節點施加相應的單元類型，這樣的開發形式導致前處理及控制面積的計算均較為繁瑣。本文提出在原地基網格的基礎上創建一種三維黏性人工邊界面單元，該單元形式是對外已有邊界節點的修正，網格劃分與原始網格一致，易于操作。

如圖4所示為所開發單元的創建過程，通過節點1-4所圍成的區域實現人工邊界的施加。首先，計算該單元的面積并得到每個節點的控制區域；然后，可由物理元件系數得到邊界節點在該單元區域內對阻尼陣的貢獻并將其團聚到相應節點，從而可以得到阻尼對角矩陣，并將其添加到整體的計算矩陣之中。依照圖4所示的坐標系，邊界單元阻尼矩陣的表達式如下：

圖4 三維黏性邊界單元示意圖Fig.4 Sketch of three-dimensional viscoelastic boundary element

式中：CBN為所施加阻尼器的切向系數；CBT為法向系數；Ai（i=4）為邊界節點控制區域的面積。

另外，在單元子程序編制中，可方便獲取節點的坐標以及每一荷載步的響應（位移、速度以及加速度），從而可依據式（1）實現節點等效荷載的施加。

2.3 時域計算模型的建立

非巖性地基條件下的核電結構地震響應時域模型在ANSYS的建立，三維等價線性單元及黏性人工邊界單元用戶子程序的編制與連接是關鍵技術。子程序的編制是對標準ANSYS的功能添加，需編譯連接成功后才能實現新的單元類型模塊的嵌入，進而結合成熟的有限元軟件前后處理及求解器，實現本文模型的建立，詳細操作過程如表1所示。

表1 所開發時域模型的建立流程詳解Table 1 Flowchart of the establishment of proposed model

3 數值驗證

本節以某擬建內陸核電場址的地基資料為基礎，進行自由場動力分析，并分別采用所開發的計算模型與商用軟件SuperFLUSH/2D進行對比，驗證本文提出模型的合理性及精確性。

3.1 材料參數

依據地質勘察資料，該廠址地基水平層狀分布較為顯著，為簡化計算，假定為水平成層地基。各分層的材料屬性如表2所示。

表2 地基計算參數Table 2 Calculation parameters of soil

圖5給出由試驗得到的各土層對應的G-γ和D-γ曲線。

圖5 各土層對應的G-γ和D-γ曲線Fig.5 G-γand D-γcurves corresponding to each layer

3.2 有限元模型

基于ANSYS的三維土柱自由場有限元模型如圖6（a）所示，底部為半無限地基，四周為自由邊界，選取截斷范圍為4.0 m×4.0 m×47.7 m，各土層的深度由表1所示每層土厚度可知。基于SuperFLUSH/2D的二維有限元模型如圖6（b）所示，利用該模型計算時，兩個水平方向的響應通過二次動力計算獲得。在計算分析中，在基巖處的網格最大高度為3.0 m，在軟弱土層局部網格最大高度為1.0 m，滿足影響頻率范圍內簡諧波傳播對網格尺寸的要求。

圖6 土柱有限元模型Fig.6 Finite element model of soil pile

3.3 輸入地震動

根據實際場地地震安評報告，計算中所施加的地震動加速度時程曲線如圖7所示。水平X向

圖7 輸入的加速度時程曲線Fig.7 Enter acceleration time curve

不同方向隨地基深度的加速度幅值，不難看出，如圖9所示，新開發模型的計算結果與SuperFLUSH/2D的計算結果在三個方向上總體趨勢基本一致，數值上也差別較小。可進一步驗證的加速度峰值為0.105g，水平Y向的加速度峰值為0.104g，豎直Z向的加速度峰值為0.124g，時間步長為0.01 s，總持時為40 s。

圖9 加速度峰值沿高程的分布圖Fig.9 The distribution of Peak acceleration along the height

3.4 計算結果分析

通過對不同高程處的加速度反應譜及加速度幅值隨深度的變化規律來考察所開發模型。為方便比較，將所開發模型命名為SANSVE。

圖8所示為兩種計算模型在不同高程（頂部、中部和底部）處不同方向的加速度反應譜對比，不難發現，新開發模型SANSVE的計算結果與商用軟件SuperFLUSH/2D的計算結果整體趨勢相同，基本完全吻合，從數值上來看也差別很小，在水平和豎直方向三個高程的反應譜峰值相差最大分別0.000 212g、0.000 189g和0.000 261g。

圖8 不同高程處加速度反應譜對比Fig.8 Comparison of acceleration response spectra at different elevations

所開發計算模型的合理有效性。

4 工程應用

以該擬建核電廠所采用AP1000堆型為研究對象，根據水平成層地基的分布特征，進行非巖性地基條件下的核島廠房結構三維地震響應分析。

4.1 計算模型

核島廠房結構為質點-梁模型，具體參數可參見文獻［8］，地基參數取值如表2所示，而各土層對應的G-γ和D-γ曲線如圖5所示。地基模擬范圍為從筏板邊緣向四周延伸筏板長度的1倍，深度方向則自地表面至玄武巖層頂面取47.7 m。如圖10所示所建立核島廠房結構-地基系統有限元模型和核島廠房質點-梁模量（帶筏板）的局部模型，整體有限元模型共有281 277個單元、253 650個節點。

4.2 計算工況

為更好地研究土體的非線性特征和無限地基輻射阻尼效應，選取三種不同的計算工況進行對比分析，具體如表3所示。

表3 計算工況Table 3 Calculation conditions

限于篇幅，選取包絡性節點的計算結果作比較，所選節點位置如圖10（a）所示，分別為標高+333.13 m處的PCS水箱頂部（節點48）、標高+169.00 m處的內部結構頂部（節點210）。

圖10 核島廠房結構-地基系統有限元計算模型Fig.10 Structure-foundation system finite element calculation model of nuclear power plant site

4.3 計算結果分析

圖11、圖12所示為節點48和節點210的加速度反應譜曲線，由于節點位置的不同，兩節點的響應分析會有較明顯的差異，因此該分析的結果具有很好的代表性。通過對比工況1和工況2可以看出：采用黏性人工邊界模擬無限地基輻射阻尼后，核島的加速度峰值明顯降低，盡管兩者的加速度反應譜曲線趨勢大致相同，但是加速度的量值卻降低了，水平方向反應譜加速度峰值最大降低率約為31.2%，豎直方向反應譜加速度峰值最大降低率約為28.6%。

另外，由表4可知，工況1和工況2雖然出現峰值的頻率基本相同，但兩者的加速度峰值相差很大。這主要是因為剛性地基模型不能吸收散射波的波動能量，對模擬輻射阻尼效應具有一定的局限性。而黏性邊界模型可以很好地吸收散射波的逸散能量，從而使核島結構的響應減小，更加趨于安全。

表4 加速度反應譜的峰值加速度對比Table 4 Peak acceleration comparison of response spectrum

考慮土體非線性特征對核島結構響應的影響，從圖11、圖12可以看出，工況2和工況3之間存在差值，較之工況2不考慮土體的非線性特征的影響，工況3中考慮土體非線性效應影響后，加速度峰值顯著降低，其峰值頻率也移至較低頻率。在水平和豎直兩個方向的加速度幅值最大分別下降0.3g和0.4g，峰值頻率最大分別偏移了0.9 Hz和4.5 Hz。這主要是因為考慮土體非線性特征后，地基剛度變小，而阻尼比會增大，從而導致地基主頻減小，反應譜峰值也隨之減小。

圖11 不同工況條件下節點48處加速度反應譜對比Fig.11 Comparison of acceleration reflection spectra under different working conditions of node 48

圖12 不同工況條件下節點210的加速度反應譜對比Fig.12 Comparison of acceleration reflection spectra under different working conditions of node 210

由工況1和工況3的對比可知，考慮黏性邊界和近場非線性時，加速度反應譜有著非常明顯的變化，水平向加速度幅值最大下降了0.7g，峰值頻率最大偏移了0.9 Hz。豎直向加速度幅值最大下降了0.7g，峰值頻率最大偏移了4.5 Hz。反應譜峰值最大減小幅度約為46.7%，減小幅值約為0.7g，且加速度峰值頻率向較低頻偏移。進一步說明，在核島廠房結構軟土地基抗震安全性分析中，考慮無限地基輻射阻尼效應及土體非線性的必要性。

5 結 論

（1）基于UPFs的二次開發特點，在ANSYS進行非巖性地基條件下核島廠房結構三維地震響應分析計算模型的開發具有可行性，經算例驗證，該模型表現出良好的合理性及精確性，可操作性強，具有較高的工程實用價值。

（2）非巖性地基對核島廠房結構的地震響應有較明顯的影響，在抗震安全評價中考慮近場地基非線性和遠場無限地基輻射阻尼是十分必要的。

（3）基于通用有限元軟件的開發，可在滿足專業要求的前提下，最大限度地發揮軟件的優勢，利用其豐富的單元類型及強大的求解能力，解決更復雜的實際工程問題。

（4）所開發模型具有直接法的優點，可對不符合抗震設防要求的地基進行加固處理及優化數值模擬，考慮樁-土-結構的動力相互作用，進一步對地基處理后的方案進行安全性評價。

致謝 感謝同濟大學建筑設計研究院(集團)有限公司陳曦高級工程師在本項目中對作者的悉心指導和幫助。