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

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

0 引 言

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

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

近年來(lái)，非巖性地基［7］場(chǎng)地條件下核島廠房地震響應(yīng)分析也備受?chē)?guó)內(nèi)外學(xué)者關(guān)注，并取得了很多有價(jià)值的成果。王桂萱，盛超等［8］建立了考慮近場(chǎng)非線性和土-結(jié)構(gòu)相互作用的核島廠房三維地震響應(yīng)分析模型，研究非巖性地基對(duì)結(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］采用等效線性和黏性人工邊界模擬近場(chǎng)地基特征，對(duì)土質(zhì)地基條件下的核島廠房結(jié)構(gòu)進(jìn)行了擬三維地震分析；然而，擬三維分析需在兩垂直平面內(nèi)簡(jiǎn)化結(jié)構(gòu)模型，也不能實(shí)現(xiàn)三方向上地震波的輸入，不能真實(shí)反映土體運(yùn)動(dòng)對(duì)變形特性參數(shù)的影響。

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

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

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

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

1.1 近場(chǎng)地基非線性模擬

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

等價(jià)線性法［12-13］通過(guò)線性迭代求解非線性問(wèn)題，關(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)變變化滿(mǎn)足收斂要求。具體計(jì)算過(guò)程如圖2所示。

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

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

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

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

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

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

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

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

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

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

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

2.2 黏性人工邊界單元

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

如圖4所示為所開(kāi)發(fā)單元的創(chuàng)建過(guò)程，通過(guò)節(jié)點(diǎn)1-4所圍成的區(qū)域?qū)崿F(xiàn)人工邊界的施加。首先，計(jì)算該單元的面積并得到每個(gè)節(jié)點(diǎn)的控制區(qū)域；然后，可由物理元件系數(shù)得到邊界節(jié)點(diǎn)在該單元區(qū)域內(nèi)對(duì)阻尼陣的貢獻(xiàn)并將其團(tuán)聚到相應(yīng)節(jié)點(diǎn)，從而可以得到阻尼對(duì)角矩陣，并將其添加到整體的計(jì)算矩陣之中。依照?qǐng)D4所示的坐標(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à)線性單元及黏性人工邊界單元用戶(hù)子程序的編制與連接是關(guān)鍵技術(shù)。子程序的編制是對(duì)標(biāo)準(zhǔn)ANSYS的功能添加，需編譯連接成功后才能實(shí)現(xiàn)新的單元類(lèi)型模塊的嵌入，進(jìn)而結(jié)合成熟的有限元軟件前后處理及求解器，實(shí)現(xiàn)本文模型的建立，詳細(xì)操作過(guò)程如表1所示。

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

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

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

3.1 材料參數(shù)

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

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

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

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

3.2 有限元模型

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

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

3.3 輸入地震動(dòng)

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

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

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

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

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

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

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

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

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

4 工程應(yīng)用

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

4.1 計(jì)算模型

核島廠房結(jié)構(gòu)為質(zhì)點(diǎn)-梁模型，具體參數(shù)可參見(jiàn)文獻(xiàn)［8］，地基參數(shù)取值如表2所示，而各土層對(duì)應(yīng)的G-γ和D-γ曲線如圖5所示。地基模擬范圍為從筏板邊緣向四周延伸筏板長(zhǎng)度的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ì)算工況

為更好地研究土體的非線性特征和無(wú)限地基輻射阻尼效應(yīng)，選取三種不同的計(jì)算工況進(jìn)行對(duì)比分析，具體如表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é)果具有很好的代表性。通過(guò)對(duì)比工況1和工況2可以看出：采用黏性人工邊界模擬無(wú)限地基輻射阻尼后，核島的加速度峰值明顯降低，盡管兩者的加速度反應(yīng)譜曲線趨勢(shì)大致相同，但是加速度的量值卻降低了，水平方向反應(yīng)譜加速度峰值最大降低率約為31.2%，豎直方向反應(yīng)譜加速度峰值最大降低率約為28.6%。

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

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

考慮土體非線性特征對(duì)核島結(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)譜對(duì)比Fig.11 Comparison of acceleration reflection spectra under different working conditions of node 48

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

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

5 結(jié) 論

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

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

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

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

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