基于ANSYS平臺(tái)的非巖性地基條件下核島廠房結(jié)構(gòu)三維地震響應(yīng)分析

尹訓(xùn)強(qiáng) 趙文燕 王桂萱

（大連大學(xué)土木工程技術(shù)研究與開發(fā)中心，大連116622）

0 引 言

核電是高效型、安全型能源，然而由于國內(nèi)不同地區(qū)地質(zhì)條件差異較大，以及近些年來我國核電的飛速發(fā)展導(dǎo)致核電廠址選址問題日益突出，因此，內(nèi)陸非巖性廠址已成為目前核電站建設(shè)的必然選擇［1-2］。其中，無限地基和近場地基非線性特征對核電結(jié)構(gòu)地震響應(yīng)的影響成為核電廠選址的關(guān)鍵技術(shù)問題。

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

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

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

1 非巖性地基條件下核島廠房結(jié)構(gòu)地震響應(yīng)分析模型

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

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

1.1 近場地基非線性模擬

在強(qiáng)震作用下，非巖性地基條件下土體的非線性動(dòng)力特征顯著，直接非線性時(shí)域方法由于本構(gòu)模型的復(fù)雜性以及數(shù)值計(jì)算的穩(wěn)定性在核電工程中并未得到廣泛應(yīng)用。美國土木工程師協(xié)會(huì)（ASCE）建議的等價(jià)線性法［11］概念明確且優(yōu)勢突出，是不可替代的主流方法。

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

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

1.2 黏性人工邊界模型及地震動(dòng)輸入

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

黏性人工邊界模型［14］的核心是在地基的邊界處布置阻尼器吸收邊界處的反射波波動(dòng)能量（圖1），并通過施加節(jié)點(diǎn)等效荷載模擬截?cái)嗤膺吔缣帉?shí)際波場的應(yīng)力邊界條件。在三維計(jì)算模型中，輸入地震動(dòng)轉(zhuǎn)化為等效荷載的具體計(jì)算公式為

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

對于非均質(zhì)近場地基的地震動(dòng)輸入，可先進(jìn)行自由場分析得到外邊界處的節(jié)點(diǎn)響應(yīng)，而后利用式（1）進(jìn)行各節(jié)點(diǎn)處的等效荷載計(jì)算，在具體數(shù)值實(shí)現(xiàn)中，可建立土柱有限元模型進(jìn)行自由場分析［15］。

2 基于ANSYS平臺(tái)的地震響應(yīng)分析模型開發(fā)

2.1 三維等價(jià)線性單元

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

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

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

2.2 黏性人工邊界單元

以往的文獻(xiàn)中，人工邊界的開發(fā)往往是將物理元件參數(shù)乘以每個(gè)節(jié)點(diǎn)的控制面積，然后在每個(gè)節(jié)點(diǎn)施加相應(yīng)的單元類型，這樣的開發(fā)形式導(dǎo)致前處理及控制面積的計(jì)算均較為繁瑣。本文提出在原地基網(wǎng)格的基礎(chǔ)上創(chuàng)建一種三維黏性人工邊界面單元，該單元形式是對外已有邊界節(jié)點(diǎn)的修正，網(wǎng)格劃分與原始網(wǎng)格一致，易于操作。

如圖4所示為所開發(fā)單元的創(chuàng)建過程，通過節(jié)點(diǎn)1-4所圍成的區(qū)域?qū)崿F(xiàn)人工邊界的施加。首先，計(jì)算該單元的面積并得到每個(gè)節(jié)點(diǎn)的控制區(qū)域；然后，可由物理元件系數(shù)得到邊界節(jié)點(diǎn)在該單元區(qū)域內(nèi)對阻尼陣的貢獻(xiàn)并將其團(tuán)聚到相應(yīng)節(jié)點(diǎn)，從而可以得到阻尼對角矩陣，并將其添加到整體的計(jì)算矩陣之中。依照圖4所示的坐標(biāo)系，邊界單元阻尼矩陣的表達(dá)式如下：

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

式中：CBN為所施加阻尼器的切向系數(shù)；CBT為法向系數(shù)；Ai（i=4）為邊界節(jié)點(diǎn)控制區(qū)域的面積。

另外，在單元子程序編制中，可方便獲取節(jié)點(diǎn)的坐標(biāo)以及每一荷載步的響應(yīng)（位移、速度以及加速度），從而可依據(jù)式（1）實(shí)現(xiàn)節(jié)點(diǎn)等效荷載的施加。

2.3 時(shí)域計(jì)算模型的建立

非巖性地基條件下的核電結(jié)構(gòu)地震響應(yīng)時(shí)域模型在ANSYS的建立，三維等價(jià)線性單元及黏性人工邊界單元用戶子程序的編制與連接是關(guān)鍵技術(shù)。子程序的編制是對標(biāo)準(zhǔn)ANSYS的功能添加，需編譯連接成功后才能實(shí)現(xiàn)新的單元類型模塊的嵌入，進(jìn)而結(jié)合成熟的有限元軟件前后處理及求解器，實(shí)現(xiàn)本文模型的建立，詳細(xì)操作過程如表1所示。

表1 所開發(fā)時(shí)域模型的建立流程詳解Table 1 Flowchart of the establishment of proposed model

3 數(shù)值驗(yàn)證

本節(jié)以某擬建內(nèi)陸核電場址的地基資料為基礎(chǔ)，進(jìn)行自由場動(dòng)力分析，并分別采用所開發(fā)的計(jì)算模型與商用軟件SuperFLUSH/2D進(jìn)行對比，驗(yàn)證本文提出模型的合理性及精確性。

3.1 材料參數(shù)

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

表2 地基計(jì)算參數(shù)Table 2 Calculation parameters of soil

圖5給出由試驗(yàn)得到的各土層對應(yīng)的G-γ和D-γ曲線。

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

3.2 有限元模型

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

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

3.3 輸入地震動(dòng)

根據(jù)實(shí)際場地地震安評報(bào)告，計(jì)算中所施加的地震動(dòng)加速度時(shí)程曲線如圖7所示。水平X向

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

不同方向隨地基深度的加速度幅值，不難看出，如圖9所示，新開發(fā)模型的計(jì)算結(jié)果與SuperFLUSH/2D的計(jì)算結(jié)果在三個(gè)方向上總體趨勢基本一致，數(shù)值上也差別較小?？蛇M(jìn)一步驗(yàn)證的加速度峰值為0.105g，水平Y(jié)向的加速度峰值為0.104g，豎直Z向的加速度峰值為0.124g，時(shí)間步長為0.01 s，總持時(shí)為40 s。

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

3.4 計(jì)算結(jié)果分析

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

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

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

所開發(fā)計(jì)算模型的合理有效性。

4 工程應(yīng)用

以該擬建核電廠所采用AP1000堆型為研究對象，根據(jù)水平成層地基的分布特征，進(jìn)行非巖性地基條件下的核島廠房結(jié)構(gòu)三維地震響應(yīng)分析。

4.1 計(jì)算模型

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

4.2 計(jì)算工況

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

表3 計(jì)算工況Table 3 Calculation conditions

限于篇幅，選取包絡(luò)性節(jié)點(diǎn)的計(jì)算結(jié)果作比較，所選節(jié)點(diǎn)位置如圖10（a）所示，分別為標(biāo)高+333.13 m處的PCS水箱頂部（節(jié)點(diǎn)48）、標(biāo)高+169.00 m處的內(nèi)部結(jié)構(gòu)頂部（節(jié)點(diǎn)210）。

圖10 核島廠房結(jié)構(gòu)-地基系統(tǒng)有限元計(jì)算模型Fig.10 Structure-foundation system finite element calculation model of nuclear power plant site

4.3 計(jì)算結(jié)果分析

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

另外，由表4可知，工況1和工況2雖然出現(xiàn)峰值的頻率基本相同，但兩者的加速度峰值相差很大。這主要是因?yàn)閯傂缘鼗Ｐ筒荒芪丈⑸洳ǖ牟▌?dòng)能量，對模擬輻射阻尼效應(yīng)具有一定的局限性。而黏性邊界模型可以很好地吸收散射波的逸散能量，從而使核島結(jié)構(gòu)的響應(yīng)減小，更加趨于安全。

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

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

圖11 不同工況條件下節(jié)點(diǎn)48處加速度反應(yīng)譜對比Fig.11 Comparison of acceleration reflection spectra under different working conditions of node 48

圖12 不同工況條件下節(jié)點(diǎn)210的加速度反應(yīng)譜對比Fig.12 Comparison of acceleration reflection spectra under different working conditions of node 210

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

5 結(jié) 論

（1）基于UPFs的二次開發(fā)特點(diǎn)，在ANSYS進(jìn)行非巖性地基條件下核島廠房結(jié)構(gòu)三維地震響應(yīng)分析計(jì)算模型的開發(fā)具有可行性，經(jīng)算例驗(yàn)證，該模型表現(xiàn)出良好的合理性及精確性，可操作性強(qiáng)，具有較高的工程實(shí)用價(jià)值。

（2）非巖性地基對核島廠房結(jié)構(gòu)的地震響應(yīng)有較明顯的影響，在抗震安全評價(jià)中考慮近場地基非線性和遠(yuǎn)場無限地基輻射阻尼是十分必要的。

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

（4）所開發(fā)模型具有直接法的優(yōu)點(diǎn)，可對不符合抗震設(shè)防要求的地基進(jìn)行加固處理及優(yōu)化數(shù)值模擬，考慮樁-土-結(jié)構(gòu)的動(dòng)力相互作用，進(jìn)一步對地基處理后的方案進(jìn)行安全性評價(jià)。

致謝 感謝同濟(jì)大學(xué)建筑設(shè)計(jì)研究院(集團(tuán))有限公司陳曦高級(jí)工程師在本項(xiàng)目中對作者的悉心指導(dǎo)和幫助。