核電工程異形水箱流固耦合分布式參數模型與效用研究

李建波, 劉 佳, 李志遠, 林 皋

(1. 大連理工大學 海岸與近海工程國家重點實驗室,遼寧 大連 116024;2. 大連理工大學 建設工程學院,遼寧 大連 116024;3. 中國水利水電科學研究院 工程抗震研究中心,北京 100048)

2011年日本福島核事故后,核電結構的抗震分析模型變得日益精細化和體系化。在土-結構相互作用框架下,液體和屏蔽廠房之間的流固耦合效應成為結構抗震分析的重要組成部分。許多學者基于精細模擬對異形水箱的流固耦合問題開展了深入研究,發現了結構形狀對液體的晃動特性有顯著影響。

在核電結構異形水箱的流固耦合仿真分析中,主要有兩種思路:一是基于聲學-結構耦合法(Acoustic-Structural coupling,CAS)、光滑粒子流體動力學(smoothed particle hydrodynamics, SPH)[1]、耦合歐拉-拉格朗日(coupled Eulerian-Lagrangian,CEL)法[2-3]等精細模擬液體和器壁動力耦合相互作用的耦合數值模型;二是以附加質量表征流體動力學特征的集總參數模型,如在有限域水體中的Housner等[4-5]模型和無限域水體中基于Westergarrd[6]公式的動水附加質量模型。

雖然前者更為精細,但是具有建模復雜、求解耗時長等缺陷,并且其分析過程與輸入地震動密切耦聯,不適用參與復雜條件的工程計算。此外,與CAS法相比,雖然SPH和CEL方法能夠反映液體的形態,但由于流固耦合的特殊性,它常常導致土結構相互作用在穩定性方面出現嚴重問題,在CEL法中泄漏問題是經常遇到的情況,SPH法的計算精度偏低,計算耗時。因此,在流固耦合中CAS法是基于有限元的較為常用數值模擬方法,此外還有邊界元法[7]、比例邊界有限元法[8]等邊界類數值方法,都在流固耦合分析中也取得了進展。需要指出的是,邊界元法和屬于比例邊界有限元法本質上仍然屬于全耦合求解過程。

在以工程結構設計為目標的抗震分析通常使用一種標準化的穩定附加質量模型,在核電領域中也應該建立相匹配的標準化模型方法。雖然核電抗震設計規范中推薦使用Housner模型,但傳統的Housner附加質量模型主要有兩方面不足:一方面,Housner模型只有在剛性器壁下截面規則的斷面才能得到穩定解析解;另一方面,集中式附加質量模型只能反映整體效應,導致在局部應力分析中存在明顯不足。目前通過將波函數方法直接應用于Housner模型的推導公式,最終得到分布式質量模型[9]。如黃文等[10]利用Housner模型和流固耦合方法計算了圓柱形容器中的地震響應。寶鑫等[11-12]在Housner附加質量模型的基礎上研究了環形水箱、底部為錐形和圓柱形的PCS水箱,并給出了其經驗公式。

除以上方法,還有基于精細模擬流固耦合模型直接提取動水壓力的計算結果[13]將其轉化為附加質量的方法,但這種方法面臨在不同的激振信號下結果具有不確定性的問題。綜上所述,雖然等效附加質量的求解方式繁多,但Housner模型的基本建模思路仍然占據主要地位。

目前在流固耦合領域剝離附加質量模型主要存在兩個問題:其一是解耦動水附加質量與激勵地震動變量之間的解耦問題;其二現有的等效附加質量模型中常常忽略動水附加質量中的對流質量,而對流質量會對計算結果產生不確定性的影響[14]。

本文基于模態綜合法[15]提出一種復雜水箱的流固耦合分析方法,并且解耦動水質量與加速度激勵之間的關系,獲得了穩定的動水分布質量模型,并將其應用于核電異形水箱。模型具有以下特點:第一,從異形水體動力特性的理論角度出發,將流固耦合界面的動水壓力描述成動水壓力與地震激勵相乘的形式,剝離附加質量與地震激勵的耦聯關系。第二,晃動質量不再在整體器壁上進行積分,而是在各個點的控制處進行積分。同時,沖動質量的定義也發生了變化,它使用該層水體和其差值來定義,而不再使用總體差值來進行分攤,使得模型更加精確和準確。本模型通過與精細有限元模擬模型進行驗證,證實了該模型的準確性和有效性;同時,針對不同水位情況和不同格柵工況下非能動水箱的響應和減震作用進行了研究,得出了非能動水箱的較佳水位和格柵分布,為水箱的抗震和減晃措施研究提供了一種有效手段。

1 流固耦合分析的分布質量模型

1.1 液體對器壁的動水壓力

假定容器器壁為剛性,液體不可壓縮、無黏、無旋。以速度勢函數Ф(x,y,z,t)為變量的控制方程及結構交界面、自由液面的邊界條件為

(1)

(2)

(3)

圖1 水平激勵下的任意形狀液體容器Fig.1 Liquid container of arbitrary shape under horizontal excitation

(4)

在液固耦合分析中,通過設置不同的邊界條件可將液體動水壓力分解為與器壁共同運動的脈沖分量和液體自由振動的對流分量。對于包括脈沖分量和對流分量的總動水壓力,自由液面處的邊界條件由伯努利方程確定,即z=HL處的動水壓力為

p=ρLgh

(5)

式中:HL為液體深度;h為液面離開自由液面靜止位置的垂直位移。當不考慮液體自由振動的對流分量時,在水平激勵作用下液體與器壁共同運動,在自由液面不會產生豎向運動,則z=HL處的動水壓力為

p=0

(6)

求解式(1)～式(6)可以得到儲液容器內任意位置處的動水壓強p。將勢函數Φ分解成與液體晃動相關的運動勢ΦS和與液體脈動相關的沖動勢ΦU。液體對側壁總的動水壓力可表示為

(7)

基于剛性器壁假定,沖動勢ΦU函數可表示為器壁速度的線性函數。因此,求解動水壓力的問題關鍵為求解運動勢函數ΦS。

在分析域Ω上對ΦS的邊界條件和運動方程進行變分并使用格林函數進行化簡,可得控制方程

(8)

式中,φ*為權函數。液體的對流分量體現液體自由振動運動,因此ΦS可表示為液體自由振動模態的疊加

(9)

(10)

(11)

(12)

求解式(10)可以得到振型的廣義坐標Yn(t),代入式(9)中可得對流勢函數,代入式(7)可得液體對器壁的動水壓力。

1.2 流固耦合的分布質量模型

式(7)、式(11)和式(12)均為積分方程,對于簡單規則形狀的容器可直接求解,但對于異形水箱無法直接計算。針對異形水箱,采用對容器側壁進行有限元離散,對液面下的每個側壁單元上的動水壓強積分,可以得到這個單元上承擔的動水壓力值

(13)

式中,s為離散單元的面積。根據等參變化,式(13)的積分可在母單元的局部坐標下計算,通過雅可比矩陣將參變量轉換到實際單元中。

為了更加直觀地表示對流質量和脈動質量,將對流質量對應的廣義坐標采用液體晃動相對運動來表示。定義變量

(14)

水體相對于水箱晃動的相對位移表示為

(15)

則對應的液體對側壁的動水壓力表示為

(16)

(17)

(18)

每個節點處的等效脈沖質量和對流質量可通過,節點控制面積與單元面積的比值進行求解。至此對于任意形狀儲液容器液體的動水壓力分析,可以等效為分布式的對流質量和脈動質量,從而對復雜的流固耦合問題進行解耦并有效的提高計算效率。

1.3 基于建立模型流固耦合分析的計算流程

本文提出的基于分布式等效附加質量的流固耦合解耦模型,自主開發了相應程序,實現了液體脈沖質量、對流質量以及彈簧常數的自動計算和添加,及動力求解等功能。計算流程圖如圖2所示,基于提出模型的流固耦合分析主要包括兩個步驟:①對液體進行有限元建模,利用自編程序計算附加質量;②將附加質量以彈簧質量單元施加到結構整體模型中,進行地震動力計算。

圖2 分布式質量簡化模型計算流程圖Fig.2 Calculation flow chart of distributed mass simplified model

2 分布式質量模型的驗證與討論

2.1 圓柱形水箱的自振頻率及動水壓力

表1 前3階振型的固有頻率

為了驗證聲學-結構耦合法和分布式附加質量法的合理性,在圓柱形平底水箱的情況下提取了位于水箱中部位置的動水壓力,并將其與Housner模型[16]的結果進行對比,如圖3所示。由圖3可知,本文提出的分布式附加質量模型和CAS模型在分析規則圓柱形水箱時產生的等效壓力結果與Housner模型的結果是一致的,這進一步驗證了這兩種方法的合理性。

圖3 距自由液面2 m處動水壓力時程曲線Fig.3 Time history curve of dynamic water pressure 2 m away from the free liquid surface

2.2 PCS水箱的自振頻率及動水壓力

李建波等建立的AP1000模型如圖4所示,重力水箱在核電廠屏蔽廠房的頂部,質量約3 000 t,正常運行水深為9.8 m。安全殼材料采用C40標號混凝土,質量密度為2 400 kg/m3,泊松比0.17,彈性模量3.25×104MPa;液體部分密度為1 000 kg/m3,剪切模量1.96 GPa。PCS水箱和結構廠房的尺寸如表2所示。

表2 材料參數表Tab.2 Material parameter table

圖4 AP1000 PCS有限元模型剖面圖和關鍵點示意圖Fig.4 AP1000 PCS finite element model and profile Schematic diagram of key points

結構的模態和頻率是結構的固有特性,為驗證此分布式附加質量模型的正確性以及其精度,將本模型與流固耦合分析中常用的CAS數值模擬的前3階振型的固有頻率(Hz)相對比如表3所示。

表3 前3階振型的固有頻率

根據法規標準譜RG1.60[17-18]合成地震波加速度時程輸入,持續時間為27.99 s,時間步長為0.01 s,加速度峰值為0.3g。由1.2節的計算得到的脈沖質量和對流質量,以彈簧阻尼器和附加質量的形式施加在AP1000頂部的水箱側壁和斜底面,進行動力求解計算。

提取位于水箱中部位置的動水壓力,如圖5所示。不難看出,本文分布式附加質量模型與CAS模型所產生的等效壓力具有較好的一致性。

圖5 距自由液面3 m處動水壓力時程曲線Fig.5 Time history curve of dynamic water pressure 3 m away from the free liquid level

2.2.1 計算效率

為了比較兩種方法的計算效率并確保計算時間的可比性,以上的計算分析都在同一臺計算機上進行。該計算機的配置如下:內存容量為32 G,CPU主頻為2.61 GHz,擁有6核12線程。在相同的安全殼網格劃分下,不同計算所需的時間如表4所示。

表4 兩種方法計算時間比較

CAS模型的計算時間相較于附加質量模型高出一倍。分析其主要原因有兩點:首先, CAS模型在計算過程中需要同時處理液體網格,導致計算量顯著增加;其次,CAS模型需要同時考慮液體和器壁之間的接觸關系,在判斷邊界條件時需要耗費大量時間。因此,本文提出的分布式附加質量模型在計算效率上明顯優于CAS模型,顯著提升了求解效率。

2.2.2 結構響應的加速度時程

對核島結構進行動力計算并提取兩種計算方法的3個觀測點在RG1.60地震動下的加速度時程響應結果如圖6～圖8,在地震動下加速度響應峰值對比如表5所示,并分析兩個模型計算結果的相對誤差。

表5 RG1.60地震下加速度峰值對比

圖6 觀測點N16加速度時程Fig.6 Observation point N16 acceleration time history

圖7 觀測點N30加速度時程Fig.7 Observation point N30 acceleration time history

圖8 觀測點N585加速度時程Fig.8 Observation point N19 acceleration time history

通常在箱線圖中樣本數據超出上邊緣線和下邊緣線的概率低于0.7%,故使用箱線圖可以極好地反映各點加速度響應相對誤差的離群程度,并且此種統計方法將誤差分析從峰值點拓展到全時程。

如圖9所示3個觀測點中各點加速度響應的相對誤差最大值都小于4%,本分布式附加質量模型在保證理想計算精度的條件下,節省了計算時間,極大地提高了計算效率,證明了本分布式附加質量模型與CAS計算模型在模擬結果中的高度一致性。

圖9 誤差箱線圖Fig.9 Error box diagram

2.2.3 結構響應的樓層譜

根據加速度反應的時程曲線計算得出樓層反應譜,圖10～圖12列出了AP1000模型在地震動激勵下3個觀測點的樓層反應譜情況,在水平方向激勵下兩種計算模型的反應譜曲線在各個頻段基本吻合。為了更清楚地反映其差異大小,圖13繪制了3個觀測點樓層反應譜的誤差箱線圖。由箱線圖可以得出在100 Hz以內所有頻率點中絕大部分的誤差都小于5%,說明兩種計算模型結果具有一致性。

圖10 觀測點N16反應譜對比Fig.10 Observation point N16 acceleration time history

圖11 觀測點N30反應譜對比Fig.11 Observation point N30 acceleration time history

圖12 觀測點N585反應譜對比Fig.12 Observation point N19 acceleration time history

圖13 誤差箱線圖Fig.13 Error box diagram

2.2.4 安全殼應力包絡圖

圖14為本分布式附加質量模型和CAS計算模型所計算得到的安全殼最大Mises應力包絡分布云圖。由圖14可知,兩種方法計算得到的包絡云圖分布基本一致,兩種計算方法所得到的最大Mises應力值均出現在安全殼側壁轉角的位置,即為安全殼應力包絡圖紅色部分所示。其中聲學-結構耦合(CAS)法所得到的最大值為55.11 MPa,本分布式附加質量模型的應力最大值為53.73 MPa,相對誤差為2.50%。在考慮流固耦合效應的實際工程應用中,應在安全殼側壁轉角處以及安全殼底部充分考慮安全裕度。

圖14 本分布式附加質量模型(左)與聲學-結構耦合(CAS)計算模型(右)最大Mises應力分布圖Fig.14 Distributed additional mass model (left) and Acoustic-structural coupling (CAS) calculation model (right) maximum Mises stress distribution

2.3 分布式附加質量的討論

通常在動水壓力中對流質量的占比較小,由于核電結構精細模擬的必要性,在核電結構設計中仍然考慮對流質量。

圖15和圖16給出了脈沖質量和對流質量隨液體高程分布圖。圖15中曲線①表示水箱內壁的脈沖質量隨高程的分布,曲線②表示水箱外壁的脈沖質量分布。

圖15 脈沖質量隨高程分布圖Fig.15 Pulse mass distribution with elevation

圖16 對流質量隨高程分布圖Fig.16 Distribution of convective mass with elevation

本分布式附加質量模型克服了經典Housner模型的局限性,將規則水箱的求解拓展到了不規則水箱。

為得出附加質量沿容器圓周分布的規律,獲取從X軸正向逆時針偏轉30°、60°、90°的脈沖質量,如圖17所示。可以得出最大質量值截面在X軸0°方向上,并且隨著偏轉角度的增加,各截面的沿高程分布的附加質量值逐漸減少。當截面旋轉到圖17中圓形圖例截面位置時,該位置上附加質量接近于0。

圖17 附加質量隨高程分布圖Fig.17 Additional mass distribution map with elevation

3 AP1000安全殼動力響應的規律討論

異形水箱的不規則性主要體現在以下兩點:截面的不規則性和網格劃分的不規則性。在面臨液體內有夾雜的復雜工況時,例如容器存在減震格柵、容器內部存在異種液體分層、容器邊角倒角等情況,本文以格柵為代表進行探討,計算不同數目格柵情況下的加速度響應以及傾覆力矩等。

網格劃分的不規則性主要出現在對容器內部液體的網格中,一般在規則水箱中的液體網格劃分以六面體為主,很少出現不規則的網格。但是在不規則容器中,液體內部可能出現四面體單元。在參考文獻[15]中只對圓柱形平底水箱進行了數值分析驗證,本文將其推廣到了異形水箱,并使用有限元的方法建立附加質量模型可以基于四面體單元進行離散,解決網格劃分不規則性的問題。

3.1 水位對AP1000安全殼結構響應的影響

圖18設定了AP1000水箱三種不同水位的示意圖,其中圖18(a)為高儲液率水位,取水位值為9.8 m;圖18(b)為中儲液率水位,取水位值為4.9 m;圖18(c)為低儲液率水位,取水位值為2.45 m。

圖18 AP1000模型水箱不同水位工況示意圖Fig.18 Schematic diagram of different water levels of AP1000 model water tank

圖19給出了三種儲液率工況下的附加質量分布圖。圖20給出了不同水位工況下AP1000安全殼不同高度點的峰值加速度在RG1.60地震作用下的變化情況,在AP1000頂點(N16)的加速度響應隨著水箱水位的降低表現為逐漸增大的變化規律,從高儲液率下降到中儲液率時增加幅度較不明顯,但中儲液率下降到低儲液率時增加幅度較大。在AP1000結構中部位置(N30和N33),中儲液率的峰值加速度高于高儲液率和低儲液率工況的峰值加速度,隨水箱水位高度的變化規律表現為先增大再減小。在AP1000結構中下部位置(N585)在各個水箱水位工況下的峰值加速度變化幅值不大,規律變化也表現為先增大再減小。水箱的液體確實在一定程度上起到了降低加速度響應的作用,如表6所示。

表6 RG1.60地震下加速度峰值和水位高度的關系

圖19 不同水位工況下的附加質量分布圖Fig.19 Additional mass distribution maps under different water level conditions

圖20 RG1.60下加速度峰值和水位高度的關系Fig.20 Relationship between peak acceleration and water level at RG1.60

3.2 格柵對AP1000安全殼結構響應的影響

當核電工程遭遇地震作用時,PCS重力水箱中的動水壓力對結構造成不利影響甚至導致水箱破壞,因此抑制液體晃動的格柵具有極高的實際應用價值[19-20]。四格柵和八格柵AP1000水箱示意圖,如圖21所示,在工況設計中將格柵高度考慮為從容器頂部連接到容器底部,同時水箱水位都假定為9.8 m。

圖21 AP1000模型水箱四、八格柵圖Fig.21 AP1000 model water tank four, eight grid diagram

圖22給出了三種不同格柵數目工況下的附加質量分布圖。圖23表示AP1000安全殼不同觀測點位置與不同格柵工況的峰值加速度之間的變化關系。在AP1000頂點(N16)的加速度響應隨著格柵數量的增加表現為減小-增加的變化規律。格柵數量從無格柵增加到四格柵時峰值加速度降低,格柵數量從四格柵增加到八格柵時峰值加速度增加。

圖22 不同格柵數目下的附加質量分布Fig.22 Additional mass distribution under different number of grilles

圖23 RG1.60下加速度峰值與格柵數量的關系Fig.23 Relationship between peak acceleration and number of grilles at RG1.60

在AP1000結構中部位置(N30和N33),中儲液率的峰值加速度高于高儲液率和底儲液率工況,隨水箱水位高度的變化規律表現為先增大再減小。在AP1000結構中下部位置(N585)在各個水箱水位工況下的峰值加速度變化幅值不大,規律變化也表現為先增大再減小。表7可得出在水箱內部布設一定數目的格柵能在一定程度上能降低加速度響應,在計算的三種格柵數目工況中四格柵表現最佳。

表7 RG1.60地震下加速度峰值與格柵數量的關系

4 結 論

本文提出了一種適用于異形水箱的液體晃蕩動力性質模擬的分布式質量簡化模型,這一模型能夠有效地克服傳統的Housner模型僅適用于規則水箱,且過于強調整體水體晃動效應的局限性。

從理論中基于勢流理論和模態綜合法進行分析,不僅適用于不規則水體形狀,而且有效分離了液體晃動特征與激振條件之間的耦合作用。在數值方面,通過與精細流固耦合模型進行比較,驗證了該模型具有良好的精度。

從工程實踐中圍繞核電異形水箱的液固耦合動力響應分析,在不同水位及采用不同數量的格柵減晃措施下,容器壁的動力響應反映出顯著的差異,說明本文模型在工程適用性方面具有較好的參數敏感性特征,并對工程的優化設計問題具有一定的參考作用。