李建波，楊 波，李志遠，丁志新

(1.大連理工大學海岸及近海工程國家重點實驗室，大連 116024；2.大連理工大學建設工程學部水利工程學院工程抗震研究所，大連 116024；3.中廣核工程有限公司核電安全監(jiān)控技術與裝備國家重點實驗室，深圳 518172)

核電工程中存在多種儲液結構，如AP1000、華龍1 號頂部的非能動冷卻水箱，具有形狀復雜、儲液量大等特點。在地震激勵下，液體發(fā)生震蕩與容器壁之間發(fā)生復雜的相互作用，直接影響核島結構的動力特性和地震反應。因此液體和儲液器壁之間的流-固耦合效應是核電結構抗震中需重點關注的內容。

現有核電結構的抗震分析，主要有兩種手段：一是基于理論推導的等效力學模型；二是基于有限元的數值分析。核電結構的抗震分析中更關注的是地震激勵下液體晃動對結構產生的相互作用，液體自身的復雜運動在抗震設計中通常可以忽略。因此，為了避免進行復雜的流體晃動計算，核電工程抗震規(guī)范中推薦使用Housner 附加質量模型。附加質量法是將液體等效為質量點，并以質量點運動產生的慣性力等效液體對器壁的動水壓力，從而復雜的流-固耦合分析解耦為等效質量點的計算和質量點-結構的動力響應分析兩個問題，極大地提升了計算效率。這方面的開創(chuàng)性工作由GRAHAM 和RODRIGUEZ[1]完成。他們首先基于線性勢流理論首先得到了矩形水箱液體晃動等效模型的精確公式，但無窮級數表達式不便于工程應用。隨后，HOUSNER[2-3]基于物理直觀分析，對文獻[1]中的級數公式進行了簡化，得到了脈沖質量和晃動質量的近似簡化公式。Housner 模型公式簡單且能夠較準確地模擬液體晃動中占主要作用的低階模態(tài)，因此在土木、水利工程中得到了廣泛使用。GOYAL 和CHOPRA[4]給出了圓形和矩形等不同截面進水塔結構動水壓力附加質量的經驗公式，為進水塔結構抗震安全設計的初步階段提供了技術手段。FARID 和GENDELMAN[5]將圓柱形水箱中的液體等效為彈簧質量模型，研究了油罐臨界區(qū)運動自由度與機械應力之間的數學關系，并確定了儲罐的失效模式和關鍵位置。ZOU 和WANG[6]基于彈簧-質量模型提出了一種新的水平激勵下流-固耦合晃動力學模型，并建立了晃動系統的耦合動力學方程，將計算結果與已有文獻結果進行比較，驗證了提出模型的合理性。劉煥忠等[7]基于石川島播磨重工業(yè)株式會社對儲液罐的一系列動力試驗[8]結果，給出了儲液罐附加質量的經驗公式。與日本電氣協會火力發(fā)電所的耐震設計指南[9]中建議的附加質量經驗公式相比，文獻[7, 9]的經驗公式得到附加質量在分布和數值上存在很大差異。這表明，儲液容器的抗震特性還缺乏一定的共識，現有的附加質量模型難以滿足工程抗震要求。

以上附加質量模型大都是依據圓柱形容器和矩形容器假設進行的推導，直接應用到形狀復雜的核電儲液容器將對核電結構抗震設計帶來較大的不確定性。為更準確地模擬液體的晃動特性，寶鑫等[10]在Housner 附加質量模型的基礎上給出了適用于環(huán)形容器的附加質量經驗公式，并從結構振動頻率、變形和內力三個方面比較了該方法與工程界常用的 60% 和 100% 均布附加質量法的精度。李小軍等[11]和雷墉[12]給出了底部為錐形的圓柱形PCS 水箱附加質量的經驗公式，并基于ADINA 軟件討論了簡化公式的模擬效果。但這些經驗公式仍主要針對環(huán)形或錐形等單一形狀，對于核電工程中更為一般的復雜形狀儲液容器的適用性仍存在不確定性。

常采用的數值方法包括：光滑粒子流動力學法(SPH)[13-18]、耦合歐拉-拉格朗日法(CEL)[19-21]和任意拉格朗日-歐拉算法(ALE)[22-24]等。寶鑫等[25]借助Fluid80 流體單元研究了土-結構相互作用對儲液結構動力反應的影響。由于問題的復雜性，模擬液體與儲液結構的動力相互作用需要精細的網格和極短的時間步長，這對計算機的存儲和計算性能提出非常高的要求。近年來聲學-結構耦合(CAS)方法被用于流-固耦合的分析。王丕光等[26]推導了橢圓柱體地震動水壓力的解析解并利用CAS方法進行了驗證。RAWAT 等[27]采用CAS 方法，研究了三維矩形水箱和圓柱形水箱在地震激勵下的動力響應，推導了液體運動控制方程，并從晃動波高、動水壓力分布等方面進行了綜合評價。但是針對核電項目，安全殼和地基的有限元模型要百萬自由度，直接采用CAS 法考慮流-固耦合效應，需要采用顯示積分方法，步長一般為10-5s 量級，對于幾十秒長的地震動，其計算量遠超工程能接受的范圍。因此，在保證能夠滿足工程需求的前提下，將流體和結構的動力分析進行解耦從而提高流-固耦合分析的計算效率，對于核電結構動力分析具有重要的工程實踐價值。

為解決以上問題，本文基于CAS 方法提出一種流-固耦合分析的分布質量模型，兼顧計算效率和計算精度，適于核電工程中異形水箱的抗震分析?；贑AS 方法從動水壓強的剝離出了脈沖動水壓強，進而得到等效的分布質量模型，在后續(xù)分析中將以附加質量產生的慣性作用等效液體對結構的動水壓力，實現流-固耦合分析。在上述兩步計算中，計算量都要顯著小于直接的流-固耦合分析。通過比較AP1000 安全殼的CAS 模型和附加質量模型的結構響應、計算時間，驗證了提出模型的合理性和有效性。本文研究內容可為核電工程中流-固耦合抗震分析提供有益參考。

1 計算理論

本文重點關注CAS 方法計算得到的動水壓強，從而進行流-固耦合分析的解耦，參考RAWAT[27]中有限元理論部分進行推導。

1.1 液體運動控制方程

假設儲液容器的液體不可壓縮、無粘、無旋，容器壁為剛性。對于理想液體，流體運動可以采用一個速度勢函數φ(x,y,z,t)來表示，從而簡化液體波動方程的求解。三維的液體晃動，其控制方程可以用拉普拉斯方程來表示：

式中， ?2為三維拉普拉斯算子。液體晃動的加速度為ai=?2φ/?i2(i=x,y,z)。若不考慮質量力，用液體壓力p表示液體運動的控制方程為：

式中：c為流體中的聲速，c=；K為體積模量； ρ為密度。需要指出，當p為聲壓，c為聲波時，式(2)即為聲場波動方程。對于小幅晃動問題，僅考慮液體線性范圍內的晃動，則液體動水壓力p可表示為：

式中， ρL為液體的密度。方程的邊界條件主要從2 個方面進行考慮，液體和器壁接觸面B1和自由液面B2。

液體與儲液容器接觸的邊界B1如圖1 所示。在B1處的邊界條件可以表示為：

式中：n為向外指向流體區(qū)域的法線；an(x,y,z,t)為指向液體區(qū)域沿外法線方向的加速度分量。液體自由液面為B2，假定流體相對于靜止所時液面產生一個小幅晃動波高h，這個條件被稱為線性化表面波條件，其中包含了重力波。則由伯努利方程可得z=HL處的動水壓力為：

式中：g為重力加速度；HL為液體深度。自由液面小幅晃動h引起的動水壓力可表示為p=ρLgh。將液體晃動簡化為等效力學模型時，液體運動可被分解為同器壁運動一致的脈沖分量和自由振動的對流分量。當不考慮液體自由振動的對流分量時，在水平激勵作用下液體與器壁共同運動，在自由液面不會產生豎向運動，則z=HL處的動水壓力為：

由此可知，在ABAQUS 軟件中通過設置不同的邊界條件可將液體動水壓力進行有效分解。取式(6)作為邊界條件時，即設置自由液面壓力為0，則計算得到脈沖分量產生的動水壓強；若自由液面不設置邊界條件式(6)，計算得到的則為包括脈沖分量和對流分量的總動水壓強。

1.2 流-固耦合解耦分析

抗震設計中的核電結構一般包括地基、安全殼、PCS 水箱、液體四部分。本文提出的流-固耦合分析方法包括兩個步驟：第一，針對PCS 水箱和液體，采用CAS 方法計算簡諧荷載激勵下的脈沖動水壓力，從而得到分布附加質量；第二，將附加質量施加到地基、安全殼、PCS 水箱的模型中，進行抗震分析即可。計算流程圖見圖2。下文針對以上步驟進行詳細說明。

圖2 解耦思路示意圖Fig.2 Schematic diagram of decoupling method

在步驟一中，可以先將PCS 水箱以外的單元刪除，只保留PCS 相關單元，并建立三維聲學單元，在水箱底部輸入任意正弦荷載，得到容器各節(jié)點脈沖動水壓強P(t)。假設單元節(jié)點控制的面積上脈沖動水壓力均勻分布，則動水壓力需要乘以該節(jié)點的控制面積，可表示為：

式中：Fp為節(jié)點處的動水壓力；t為時刻；Anode為該節(jié)點的等效控制面積。如圖3 所示，灰色陰影部分表示節(jié)點E的等效控制面積，可表示為：

圖3 節(jié)點等效面積示意圖Fig.3 Schematic diagram of node equivalent area

式中：Ai為與該節(jié)點關聯單元的面積；di為關聯單元節(jié)點的個數。

核電結構中，儲液容器多為具有鋼襯的混凝土結構，滿足剛性假定，器壁節(jié)點的加速度與激勵相同。與液體的脈沖分量比較，動水壓力中的對流分量占比通常很小，并且對流分量的振動相位與器壁的運動相位不同，實際上可能起到減震的作用。因此，在計算等效質量時僅考慮脈沖動水壓力。當液體運動只考慮脈沖分量時，液體和側壁加速度保持一致。此時動水壓強均由脈沖質量產生，則節(jié)點處的動水壓力又可以表示為：

式中，a(t)為外部加速度激勵。由于只考慮脈沖分量，Fp(t)與a(t)為同相位。為了減小計算誤差，脈沖質量取二者的最大值進行相除得到：

基于CAS 方法計算的分布質量不同于簡單的將水體總質量取一定的百分比進行平均分布[28]。由CAS 法得到的脈沖動水壓力在離散意義上更加準確。計算等效質量時，僅增加了單元節(jié)點控制的面積上脈沖動水壓力均勻分布的假定，沒有引入更多不確定性，因此得到的等效質量在數值和分布上都具有更高的精度。

2 分布質量分析

2.1 CAS 計算模型驗證

為驗證基于CAS 方法的分布質量模型的精度和效率，采用提出的模型分析矩形儲液容器在正弦激勵下的動水壓力時程、脈沖動水壓力產生的作用力以及水體對器壁產生的總作用力，并與Housner 分布質量模型[2]的結果進行比較。

二維矩形水箱，長L和高h均為2 m，液體高HL為1.5 m，單元尺寸為0.1 m。底部加速度激勵(t)=sin(4t)，激勵持時為5 s，步長 Δt為0.01 s。提取距離自由液面0.1 m 和0.9 m 兩個位置動水壓強時程，如圖4 所示。從圖4 中可以看出，提出模型的計算結果與Housner 分布質量模型計算結果吻合性非常好。

圖4 距自由液面不同深度處動水壓力時程曲線Fig.4 Time-history curve of hydrodynamic pressure at different depths from free liquid level

為了進一步驗證本文采用的基于CAS 方法計算的模型的正確性，還提取了容器底部脈沖動水壓力產生的作用力以及水體對器壁產生的總作用力進行比較。計算結果如圖5 所示，液體晃動產生的脈沖反力時程同外部激勵相位保持一致；同時可以看出，Housner 分布質量模型計算結果與本文提出模型計算得到的結果幾乎完全吻合，證明了本文方法的合理性和準確性。

圖5 基底反力時程曲線比較Fig.5 Comparison of substrate reaction time history curve

2.2 激勵對分布質量的影響

為了驗證附加質量與激勵的關系，取半徑為3 m，高為6 m，水深為4.5 m 的圓柱形水箱作為分析對象，分別作用5 種不同的正弦激勵。正弦激勵如表1 所示，激勵1～激勵3 幅值相同，圓頻率不同；激勵3～激勵5 幅值不同，圓頻率相同。基于CAS 分布質量模型，計算各激勵作用下的脈沖附加質量，沿高程分布如圖6 所示。從計算結果可以看出， 5 種不同的激勵作用下，沿容器高程的分布質量基本一致，表明分布質量是水箱中液體晃動的固有特征，僅與水箱形狀、液體密度、液面高度等因素相關，與外荷載無關。因此基于CAS 計算分布質量時，施加形式簡單且短持時的正弦荷載激勵即可，從而達到簡化計算，提高計算效率的效果。

圖6 不同激勵作用下附加質量沿高程分布圖Fig.6 Elevation distribution of the additional mass under different excitations

3 AP1000 模型分析

3.1 AP1000 數值模型建立

3.1.1 模型建立

以第三代先進非能動反應堆AP1000 安全殼為例，進一步表明提出模型的良好精度和廣泛適用性。建立了AP1000 安全殼整體模型，為了計算的效率和后處理的簡便，核電廠模型僅考慮頂部重力水箱和安全殼兩部分，通風口、擋板等構造物在建模時忽略。AP1000 安全殼幾何特征[29]見圖7，PCS 水箱內部直徑Rin、外部直徑Rout、高度Hw和厚度Ww分別為 10.668 m、27.13 m、11.8 m 和0.6 m，水箱底部傾斜角度為15°；屏蔽廠房直徑Rs、高度Hs和厚度Ws分別為 44.2 m、71 m 和0.92 m。如圖8 所示，AP1000 安全殼整體模型采用殼單元進行模擬，共計4766 個單元，最大單元尺寸為4.29 m，其中PCS 水箱模型共計2688 個單元，最大單元尺寸為1.33 m。

圖7 AP1000 安全殼幾何尺寸Fig.7 AP1000 containment geometry

圖8 AP1000 安全殼有限元模型Fig.8 AP1000 containment finite element model

作為對照，同時建立流-固耦合整體模型，液體深度HL為9.8 m，采用聲學單元進行模擬，共計6400 個單元，最大單元尺寸為1.33 m。

3.1.2 材料參數和結構阻尼

安全殼結構材料整體采用C40 混凝土建模，材料密度為2700 kg/m3，彈性模量為30 GPa，泊松比為0.3；液體部分采用水進行模擬，密度為1000 kg/m3，剪切模量為1.96 GPa。

在ABAQUS 中進行模態(tài)分析時是采用無阻尼系統，但在動力時程分析時，需要對結構設置一定阻尼。本文選取Rayleigh 阻尼，阻尼比為ξ=5%。

3.1.3 地震動的選擇

美國ASCE4-16 規(guī)范與GB 50267-2019 規(guī)范在選取設計地震動方法上相同，取實測地震動或通過設計譜生成人工波，本文選用根據RG 1.60設計譜進行人工合成加速度時程，人工地震波持續(xù)時間為28 s，步長為0.01 s，加速度峰值為0.3g，如圖9 所示。在計算中僅考慮X方向的地震波。

圖9 RG1.60 人工地震波加速度時程Fig.9 RG1.60 artificial acceleration time history

3.2 PCS 水箱分布質量計算

按照圖3 中給出的計算思路，完成對AP1000安全殼的抗震分析。首先計算得到PCS 水箱側壁的質量分布，圖10 是PCS 水箱的示意圖。圖11所示的是水箱內外側壁各水箱沿X方向切割開后展開的分布質量包絡云圖。

圖10 PCS 水箱示意圖Fig.10 View of the PCS tank

圖11 PCS 水箱側壁分布質量分布云圖 /kgFig.11 Mass distribution of the side wall of the PCS water tank

從圖11 中可以看出，在X方向激勵作用下，內外側壁的分布質量關于Y方向表現較強的對稱性，但并不是嚴格意義上的對稱，大小上仍具有一定的區(qū)別。對于水箱外壁，分布質量在沿X(-X)方向取得最大且從側壁頂部向底部逐漸增大，在底部取得最大；分布質量從X(-X)方向往Y(-Y)軸方向呈弧線逐漸減小，在Y(-Y)軸時取得最小值為0。對于水箱內壁，分布質量在沿X(-X)方向取得最大，且從側壁頂部向底部逐漸減小，在底部取得最??；分布質量仍表現為從X(-X)方向往Y(-Y)軸方向呈弧線逐漸減小，在Y(-Y)軸時取得最小值為0。表明在單向激勵作用下，內外側壁分布質量沿激勵方向向兩側逐漸減小，在垂直于激勵方向取得最小值。

3.3 AP1000 安全殼動力計算

將3.2 節(jié)中計算得到的PCS 水箱側壁各節(jié)點的分布質量值以附加質量點的形式施加到AP1000安全殼頂部的水箱側壁對應節(jié)點進行動力分析，從而實現流-固耦合的解耦分析。分別從計算時間、安全殼響應等多個方面對本模型(后稱附加質量模型)和直接采用CAS 方法進行流-固耦合分析(后稱流-固耦合模型)計算的結果進行了比較。

3.3.1 計算效率

流-固耦合模型只能采用顯式計算方法進行計算；而附加質量模型采用隱式計算方法進行計算。為了比較兩種方法的計算效率，充分保證計算時間的可比性，以上的計算分析均在同一臺計算機中進行。計算機的配置如下：內存為64 G，CPU 主頻為3.6 GHz，6 核12 線程。安全殼網格劃分相同，計算所需要的時間如表2 所示。

表2 兩種方法計算時間比較Table 2 Comparison of the calculation time between the two methods

流-固耦合模型的計算時間要比附加質量模型高出2 個數量級。分析其主要原因有以下3 點：流-固耦合模型在計算時液體網格同時參與計算，導致計算量大大增加；在處理液體和器壁之間接觸關系，判斷邊界條件上需花費大量時間；相關文獻[30 - 31]指出，采用ABAQUS 軟件進行顯示動力計算時，考慮結構瑞利阻尼會對計算時間和計算穩(wěn)定性造成影響。因此本文提出的附加質量模型在計算效率上遠優(yōu)于流-固耦合模型，大大提升了求解效率。

3.3.2 基底反力和傾覆力矩

表3 中給出了流-固耦合模型和附加質量模型在水平地震激勵下，安全殼底部的基底反力和傾覆力矩的最大值以及計算誤差。從表中可以看出，附加質量模型計算結果均大于流-固耦合模型計算結果，基底反力最值的誤差約為10.7%，傾覆力矩最值的誤差約為18.6%。證明本文基于CAS 的解耦的附加質量模型計算結果偏保守，對工程而言具有一定的安全裕度。

表3 AP1000 底部基底反力和傾覆力矩比較Table 3 Comparison of AP1000 base reaction and overturning moment

圖12 比較的是流-固耦合模型基底反力時程和對流分量產生的基底反力時程。對流分量產生的基底反力通過總的基底反力時程減去脈沖分量產生的基底反力時程得到。從圖中可以看出，對流分量產生的基底反力在流-固耦合及附加質量模型計算得到的基底反力中占比較小，對基底反力貢獻不大。

圖12 兩種模型基底反力時程和對流分量基底反力時程比較Fig.12 Comparison of base reaction time histories and convective component base reaction time histories between the two models

表3 以及圖12 的5 s～10 s 段，流-固耦合模型基底反力計算結果明顯小于附加質量模型。其原因是，流-固耦合模型考慮了對流分量，對流分量一般與脈沖分量不同步，在一定程度上可以抵消脈沖分量對側壁的沖擊，相當于提供了一種流體阻尼，從而使結構的基底反力小于附加質量模型的計算結果。從而可以證明本方法的假設，即在安全殼結構抗震分析時，僅考慮脈沖分量作用而忽略對流分量即可滿足工程設計需要。

3.3.3 安全殼不同高程加速度時程和反應譜

圖13 給出了安全殼外壁不同高程處的加速度時程曲線。沿安全殼側壁從頂部依次往下取了4 個觀測點，編號為N1～N4，如圖8 所示。從圖13 可以看出，流-固耦合模型和附加質量模型計算得到的不同高程處加速度時程基本吻合。加速度隨著高程的增加而增加，在安全殼頂部N1觀測點處取得最大值。

圖13 觀測點N1～N4 加速度時程曲線Fig.13 Time-history curve of acceleration at observation points N1-N4

圖14 給出了兩種模型在觀測點N1～N4的樓層加速度反應譜的比較曲線。從圖14 中可以看出，在安全殼頂部觀測點N1處反應譜峰值最大，且隨著高程的降低呈現減小趨勢；附加質量模型各觀測點的反應譜值稍大于流-固耦合模型計算得到的結果。觀測點N1～N4反應譜峰值誤差分別為4.82%、10.09%、8.48%和6.51%，說明本文通過流-固耦合解耦得到的附加質量模型計算得到的樓層加速度反應譜結果偏保守。對工程設計而言，其結果是有利的。

圖14 觀測點N1～N4 樓層加速度反應譜曲線Fig.14 Response spectrum of floor acceleration at observation point N1-N4

3.3.4 安全殼不同高程加速度峰值

提取了安全殼沿激勵方向所有高程位置節(jié)點的加速度最大值，提取點位置如圖15 所示。不同高程處最大加速度和高程H之間的關系如圖16所示，從圖16 可以看出，兩種方法計算得到的加速度最大值均隨著高程的增加逐漸增加，在頂部時最大。附加質量模型計算得到的結果整體來說均大于流-固耦合模型計算所得的結果，僅在高程為38.7 m～47.2 m 的范圍內略小于流-固耦合模型。因此總體而言，附加質量模型的計算結果偏保守。

圖15 加速度最大值提取節(jié)點示意圖Fig.15 Node path for extracting maximum acceleration

圖16 加速度峰值隨高程的變化Fig.16 Variation of the peak acceleration with elevation

3.3.5 安全殼應力包絡圖

圖17(a)、圖17(b)所示分別為流-固耦合模型和附加質量模型計算得到的安全殼最大Mises 應力包絡分布云圖。從圖中可以看出，兩種方法計算得到的包絡云圖分布基本一致。在地震激勵作用下，上部PCS 水箱上應力值整體較小，在PCS 水箱和殼體交接處較大且分布較為對稱；下部安全殼殼體應力值整體偏大，在高程約為60 m 的位置處較小，呈對稱分布。

圖17 RG1.60 地震波作用下安全殼最大應力包絡云圖Fig.17 Maximum Mises stress envelope diagram under RG1.60

兩種方法最大Mises 應力值均出現在沿激勵作用方向上安全殼側壁轉角位置，即圖中引線標注部分。其中流-固耦合模型計算值為11.42 MPa，附加質量模型計算結果為10.74 MPa，誤差約為5.95%。若采用本文提出的模型計算所得應力進行配筋設計，可以通過乘以裕度系數的方法來進行考慮，即可保證計算結果的保守性。

4 結論

本文提出的模型基于CAS 方法實現了流-固耦合問題的解耦分析，提高了核電工程儲液容器的抗震分析效率。利用CAS 方法對水箱形狀適應性強的特性，提高了基于Housner 簡化模型的流-固耦合分析的精度。分別采用CAS 模型和本文提出的基于CAS 的分布質量模型對正弦激勵下簡單矩形水箱和地震激勵作用AP1000 安全殼進行了動力分析，比較了兩種方法的計算效率和計算精度等內容。主要結論如下：

(1) 采用本模型進行核電結構流-固耦合分析時，首先基于CAS 方法施加簡諧荷載得到水箱附加質量，再將附加質量施加到核電結構整體模型中，進行抗震分析，大大提升了計算效率。

(2) CAS 模型適用于任意形狀水箱的流-固耦合分析。通過設置邊界條件可以計算得到各節(jié)點位置動水壓強的脈沖分量，計算結果與推導的Housner 分布模型的解析解吻合較好。

(3) 針對AP1000 安全殼的抗震分析，本文提出的基于CAS 的分布質量模型計算效率明顯高于直接進行流-固耦合分析的CAS 模型；計算得到的最大基底反力、傾覆力矩、不同高程樓層加速度反應譜均大于CAS 模型的計算結果，結果偏于保守。兩種模型計算得到的最大應力包絡云圖基本一致。