汪 全 趙東拂,2,3,4 楊昕光
(1.北京建筑大學土木與交通工程學院,北京 100044;2.北京建筑大學大型多功能振動臺陣實驗室,北京 100044;3.北京建筑大學工程結構與新材料北京市高等學校工程研究中心,北京 100044;4.北京建筑大學北京節(jié)能減排與城鄉(xiāng)可持續(xù)發(fā)展省部共建協(xié)同創(chuàng)新中心,北京 100044;5.中冶建筑研究總院有限公司,北京 100088)
核電廠發(fā)生核事故時,安全殼作為阻止裂變產(chǎn)物向環(huán)境泄漏的最后一道屏障,確保其完整性尤為重要[1]。為了驗證安全殼的整體密封性、結構性能和強度是否滿足要求,在機組首次換料大修及后續(xù)每十年必須進行安全殼結構整體性試驗。
在該試驗中,預應力混凝土安全殼結構主要承受試驗壓力作用,并對其變形、應變、溫度、預應力鋼束力值等進行監(jiān)測,是一種典型的低荷載、小量程的精密試驗,測量精度要求較高[2]。例如,在最大試驗內壓作用下,穹頂豎向位移僅為安全殼高度的萬分之一左右,而變形測量精度需達幾十萬分之一。因此,需考慮諸如環(huán)境溫度、太陽輻射、風速等綜合環(huán)境因素對安全殼結構整體性試驗的影響,對安全殼變形進行合理的修正,以得到符合精度要求的測量結果。
目前,國內外關于環(huán)境因素作用下結構溫度場及作用效應的研究主要集中在橋梁、水工結構等工程領域[3-5],對于安全殼結構的相關研究尚不夠深入。林志伸等通過理論推導和編程計算,分析了外圍環(huán)境因素作用下安全殼的溫度場,計算結果與實測結果吻合較好[6]。劉世豪通過有限元方法分析了施工初期安全殼的溫度場及溫度應力,更準確地預測了安全殼整體和局部開裂情況[7]。黎鵬飛等采用大型有限元軟件對安全殼分別進行了正常運行、嚴重事故工況下的穩(wěn)態(tài)和瞬態(tài)分析,模擬了安全殼的溫度場變化,確定了最不利溫度效應[8]。張冰基于熱傳導和熱對流模擬計算了安全殼在整體密封試驗時內部溫度場分布情況,驗證了安全殼34 m以上空間溫度出現(xiàn)分層分布現(xiàn)象[9]。
綜上,為提高安全殼整體性試驗的變形測量精度,并對其進行合理的修正,研究環(huán)境溫度、太陽輻射等綜合環(huán)境因素作用下預應力混凝土安全殼的變形規(guī)律是十分必要的。通過有限元模擬,研究環(huán)境溫度、太陽輻射、輻射換熱以及不同熱邊界條件施加方法對安全殼變形的影響規(guī)律,探明預應力混凝土安全殼結構在綜合環(huán)境因素下的變形規(guī)律,研究可為安全殼結構整體性試驗與核電廠安全運行提供理論與技術支持。
基于傅里葉熱傳導理論[10],忽略混凝土內部產(chǎn)生的熱量,可得熱傳導方程如下:
(1)
式中:λ為混凝土的導熱系數(shù),W·m-1·K-1;T為結構表面溫度;n為表面法線方向;I、Qc、Qr分別為太陽總輻射、對流換熱和輻射換熱的熱流密度,W·m-2。
環(huán)境因素作用下安全殼的熱邊界可分為三類:太陽輻射、對流換熱和輻射換熱。太陽輻射包括太陽直射輻射、散射輻射以及太陽輻射到達地面后的反射輻射。對流換熱為安全殼與內外空氣的熱交換。輻射換熱包括結構向環(huán)境發(fā)射的輻射和環(huán)境向結構發(fā)射的輻射。
上述三類熱邊界施加主要有兩種方法:1)分別施加法[11],即將太陽輻射、對流換熱和輻射換熱邊界條件通過不同的方式同時施加于結構上;2)綜合系數(shù)法[12],即將太陽輻射、對流換熱和輻射換熱轉化為綜合溫度和綜合換熱系數(shù)進行施加。
1.2.1熱邊界分別施加
安全殼熱邊界條件復雜,為了簡化計算,將上述不同熱邊界條件使用Robin邊界進行求解[13]:
(2)
式中:α為混凝土太陽輻射吸收系數(shù);I為太陽總輻射強度,W·m-2;hc、hr分別為對流換熱系數(shù)和輻射換熱系數(shù),W·m-2·K-1;Ta、T分別為環(huán)境溫度和結構表面溫度,℃。
對流換熱系數(shù)與結構表面粗糙度、風向、風速等有關,其中風速影響最大,聶玉東對眾多對流換熱系數(shù)進行回歸分析[14],并按式(3)計算對流換熱系數(shù):
hc=4v+5.4
(3)
式中:v為風速,m/s。
輻射換熱主要與輻射發(fā)射者和輻射接受者的溫度以及物體有關,輻射換熱系數(shù)可由式(4)計算[15]:
hr=εC0(546+Ta+T)[(273+Ta)2+(273+T)2]
(4)
式中:ε為物體輻射發(fā)射率,混凝土結構取值為0.88,空氣取值為0.82;C0為斯蒂芬-玻爾茲曼常數(shù),其值為5.67×10-8W·m-2·K-4。
1.2.2綜合換熱系數(shù)和綜合溫度
將式(2)化簡成下式:
(5)
(6)
式中:ha為綜合換熱系數(shù);T′a為綜合溫度,式(6)表示環(huán)境溫度疊加太陽輻射、對流換熱和輻射換熱的一個綜合效果。
安全殼接收到的太陽輻射主要有三部分:太陽直射輻射、太陽散射輻射、太陽直射輻射和散射輻射到達地面后的反射輻射。
1.3.1太陽直射輻射強度
太陽輻射穿越大氣層時,不改變原來傳播方向的為太陽輻射直射強度[16]。一般認為太陽直射輻射強度在穿越大氣層時按指數(shù)規(guī)律衰減,Kehlbeck基于Bouguer-Lambert定律提出一種冪指數(shù)模型用于計算光線垂直的平面太陽直射輻射強度[17],即:
Id=I0Pm
(7)
其中P=0.9tu·ta
m=1/sinh
式中:I0為太陽常數(shù),取1 353 W·m-2;P為大氣透明系數(shù)[18];tu為林可氏渾濁度系數(shù);ta為相對壓強;m為大氣質量;h為太陽高度角。
1.3.2太陽散射輻射強度
太陽輻射穿越大氣層時被大氣層散射部分中的一部分從各個方向散射到結構表面的為太陽散射輻射強度[16]。太陽散射輻射強度通常與結構表面的方位角無關與表面傾角有關,水平面散射強度可按下式計算[19]:
Is=(0.271I0-0.942Id)sinh
(8)
其他傾斜面散射強度:
(9)
式中:θ為表面傾角,即斜面與水平面之間的夾角,筒身θ=90°。
1.3.3地面反射輻射強度
太陽直射輻射和散射輻射到達地面后,一部分被地面吸收,一部分被地面反射回去,反射回去的為地面反射輻射強度[20]。反射輻射強度可按下式計算:
(10)
式中:r為地面太陽輻射反射系數(shù),取0.2。
綜上所述,結構物接收到的太陽總輻射強度可按下式計算:
I=Idcosφ+Is+Ir
(11)
式中:φ為太陽入射角,穹頂可認為cosφ=sinh。
劉誠通過太陽輻射模型計算得到北京地區(qū)晴天水平地面的太陽輻射強度與實測值吻合良好[11],說明了該模型的適用性和正確性。
某核電站安全殼結構形式為帶密封鋼襯里的預應力混凝土結構,外形為帶圓穹頂?shù)膱A柱形筒體,剖面見圖1。該安全殼具有4個扶壁柱,沿筒均勻分布,且設備閘門孔豎向中線與206.666 gr(gr為角度單位,1 gr=0.9°)扶壁柱中線重合,筒體內徑37 m,壁厚0.9 m。穹頂內半徑為24 m,壁厚0.8 m,扁球殼。鋼襯里厚6 mm,底板厚5.5 m。安全殼混凝土強度等級PS40。由廊道頂部至筒體頂部布置144根豎向預應力鋼束導管;環(huán)向布置223根鋼束導管,每根鋼束呈360°包角,均勻錨固在4個扶壁柱上。穹頂設置174根鋼束導管,鋼束錨固于環(huán)梁。鋼絞線穿入鋼管,張拉后注入水泥漿,避免鋼絞線的腐蝕。
a—安全殼結構豎向剖面(高程單位為m,其他單位為mm);b—安全殼結構水平剖面(+22.9 m)。圖1 安全殼結構剖面Fig.1 Structural cross-section of containment
采用大型通用有限元軟件ABAQUS進行計算分析,預應力混凝土安全殼有限元模型如圖2所示??紤]到計算精度和計算成本,筒身和穹頂選用溫度位移耦合的C3D8T和C3D6T三維實體單元。預應力鋼束選用溫度位移耦合的T3D2T桿單元,預應力鋼束采用Embedded約束到混凝土中。鋼內襯選用溫度位移耦合的S4T和S3T殼單元進行模擬。鋼內襯與混凝土采用Tie綁定約束。
a—混凝土殼體;b—鋼內襯;c—預應力鋼束。圖2 預應力混凝土安全殼有限元模型Fig.2 Finite element model of pre-stressed concrete containment
在ABAQUS中可將太陽輻射定義為荷載,通過Thermal中的surface heat flux進行定義,屬于第二類熱邊界條件。對流換熱可在ABAQUS Interaction模塊,通過surface film condition定義,屬于第三類熱邊界條件。輻射換熱根據(jù)式(2)近似通過第三類熱邊界條件在Interaction模塊中surface film condition定義。如采用綜合換熱系數(shù)和綜合溫度,將三種熱邊界轉化成對流換熱邊界,可在ABAQUS Interaction模塊,通過surface film condition進行定義,綜合換熱系數(shù)即film coefficient,綜合溫度即sink temperature。
混凝土太陽輻射吸收系數(shù)α為0.6[21]。其他材料參數(shù)見表1。
表1 材料參數(shù)Table 1 Materials parameters
計算模擬時間選用5月下旬某連續(xù)6日,起算時間為上午10時,計算總時長為144 h,天氣晴,氣溫18~30 ℃,風力3級。安全殼結構初始溫度為24 ℃。
假設白天環(huán)境溫度按余弦變化,夜間環(huán)境溫度按線性變化,環(huán)境溫度可按下式計算[11]:
Ta(t)=
(12)
式中:Tmin和Tmax分別為日最低和最高氣溫,分別為18 ℃和30 ℃;tr和ts分別為日出和日落時間,以小時為單位,按24 h制計算,分別為6時和18時;td和tn分別為白天和夜間時長,td=ts-tr,tn=24-td。假設安全殼殼內空氣溫度保持恒定為24 ℃。
假設計算周期內風力為3級,故安全殼外表面風速取4 m/s,由式(3)計算得到外表面對流換熱系數(shù)為21.4 W·m-2·K-1。安全殼內部風速為0,內部對流換熱系數(shù)為5.4 W·m-2·K-1。
由式(4)可知輻射換熱系數(shù)是一個關于結構溫度和環(huán)境溫度的函數(shù)。文獻[14]表明,當輻射換熱系數(shù)取一個固定近似值時對計算結果影響較小,滿足精度要求。因此取輻射換熱系數(shù)近似值為5.4 W·m-2·K-1。
計算分析了4種工況下預應力混凝土安全殼的變形響應。工況1:將太陽輻射、對流換熱和輻射換熱三種熱邊界以不同施加方式同時施加在結構表面。工況2:將三類熱邊界轉化成綜合換熱系數(shù)和綜合溫度施加在結構表面。工況3:忽略太陽輻射,其他設置同工況1。工況4:忽略輻射換熱,其他設置同工況1。本文根據(jù)實際情況作了如下簡化與假設:1)穹頂不考慮太陽輻射到達地面后的反射輻射;2)筒身由于受到周圍核電廠房的遮擋,假設筒身不受太陽直射輻射影響。具體工況見表2。
表2 安全殼熱流邊界詳情Table 2 Details of containment heat flow boundary
通過ABAQUS進行了4種工況的模擬計算,各工況的安全殼位移等值線云圖如圖3~6所示。在4種工況下,隨著環(huán)境溫度的升高,安全殼變形整體外擴,隨著環(huán)境溫度的降低,安全殼變形整體回縮。在工況1、2和4下穹頂豎向位移最大值位于穹頂中心附近,而在工況3下穹頂豎向位移最大值位于穹頂邊緣靠近閘門一側。在4種工況作用下,安全殼筒身水平徑向位移和豎向位移最大值均位于設備閘門附近。
a—第104小時水平徑向位移;b—第104小時豎向位移。圖3 工況1下安全殼位移等值線云圖 mFig.3 Contour cloud of containment displacement under the working condition one
a—第104小時水平徑向位移;b—第104小時豎向位移。圖4 工況2下安全殼位移等值線云圖 mFig.4 Contour cloud of containment displacement under the working condition two
a—第104小時水平徑向位移;b—第104小時豎向位移。圖5 工況3下安全殼位移等值線云圖 mFig.5 Contour cloud of containment displacement under the working condition three
a—第104小時水平徑向位移;b—第104小時豎向位移。圖6 工況4下安全殼位移等值線云圖 mFig.6 Contour cloud of containment displacement under the working condition four
為研究上述4種工況下安全殼位移變化情況,分別選取了穹頂不同高程的2個點和筒身同一高程、不同方位的3個點,并輸出其豎向位移和水平徑向位移的時程曲線。圖7為點位示意,A點位于穹頂中心外表面,B點位于穹頂360.5 gr方向,高程54.56 m處外表面,C、D、E點位于筒身高程26 m,方向分別為352.8,255.4,119.5 gr處外表面。
a—穹頂點位示意;b—筒身點位示意。圖7 點位示意Fig.7 Point location schematic
圖8為4種工況下穹頂A點豎向位移時程曲線。由圖可知,3 d(上午10時至次日上午10時為1 d)后安全殼變形基本趨于平穩(wěn),初始溫度場對結構變形的影響逐漸可忽略,故后續(xù)分析取第4~6天的數(shù)據(jù)進行各工況對比分析,以消除初始溫度場對結構變形的影響[22]。
圖8 4種工況下穹頂A點豎向位移時程曲線Fig.8 Time course curve of vertical displacement at point A of the containment dome under four working conditions
圖9為4種工況下穹頂豎向位移時程曲線。由此可知,各工況下穹頂豎向位移隨環(huán)境溫度近似呈正弦曲線變化,且較環(huán)境溫度變化滯后約2 h。工況1下穹頂A、B兩點豎向位移平均幅值(波峰值與波谷值之差的一半)分別為0.459 mm和0.399 mm;工況2下穹頂A、B兩點豎向位移平均幅值分別為0.454 mm和0.395 mm。上述兩種工況下穹頂A、B兩點豎向位移平均幅值相差不足1%。這說明不同熱邊界條件施加方法對穹頂位移的影響可忽略。有學者通過熱邊界分別施加法[11]或綜合系數(shù)法[13-14]進行結構溫度場研究,模擬結果均與實測值吻合良好。
a—穹頂A點;b—穹頂B點。圖9 4種工況下安全殼穹頂豎向位移時程曲線Fig.9 Time course curve of vertical displacement of containment dome under four working conditions
工況3下穹頂A、B兩點豎向位移平均幅值分別為0.235 mm和0.208 mm,較工況1分別減小了49%和48%。這是由于工況3忽略了太陽輻射,導致穹頂豎向位移平均幅值較工況1大幅度減小。由此可知,太陽輻射是影響安全殼變形的主要因素,在計算時不可忽略。
工況4下穹頂A、B兩點豎向位移平均幅值分別為0.444 mm和0.387 mm,與工況1相差約4%。這說明輻射換熱是影響安全殼變形的次要因素,為了保證模擬的準確性,應在分析中考慮輻射換熱對安全殼結構變形的影響[11-13]。
圖10為4種工況下筒身水平徑向位移時程曲線。同穹頂豎向位移變化規(guī)律相似,各工況下筒身水平徑向位移隨環(huán)境溫度近似呈正弦曲線變化,較環(huán)境溫度變化滯后約2 h。工況1下筒身 C、D和E點水平徑向位移平均幅值分別為0.062,0.058,0.076 mm;工況2下筒身C、D和E點水平徑向位移平均幅值分別為0.062,0.057,0.075 mm。因此,上述兩種工況下筒身水平徑向位移平均幅值幾乎相同,與穹頂豎向變形規(guī)律一致,說明兩種熱邊界施加方法對安全殼結構變形的影響可忽略。
a—筒身C點;b—筒身D點;c—筒身E點。圖10 4種工況下安全殼筒身水平徑向位移時程曲線Fig.10 Time course curve of horizontal radial displacement of containment barrel under four working conditions
工況3下筒身C、D和E點水平徑向位移平均幅值分別為0.049,0.046,0.060 mm,較工況1分別小21%、20%和21%。由于筒身忽略了太陽直射輻射,工況1和工況3的筒身水平徑向位移平均幅值較穹頂豎向幅值有所減小,但太陽散射輻射和反射輻射對安全殼結構水平徑向變形的影響仍較大,應在分析時予以考慮。
工況4下筒身水平徑向位移平均幅值分別為0.058,0.054,0.069 mm,與工況1相差約9%,同樣說明輻射換熱對安全殼結構變形存在一定的影響。
綜上,通過各工況的對比分析,說明安全殼位移隨環(huán)境溫度近似呈正弦曲線變化,環(huán)境溫度和太陽輻射是影響安全殼變形的主要環(huán)境因素,輻射換熱對安全殼變形的影響有限,不同熱邊界條件施加方法對安全殼變形的影響可忽略。
通過有限元數(shù)值模擬,研究了環(huán)境溫度、太陽輻射、輻射換熱以及不同熱邊界條件施加方法對安全殼變形的影響,闡明了預應力混凝土安全殼結構在綜合環(huán)境因素下的變形規(guī)律,主要結論如下:
1)各工況下安全殼位移隨環(huán)境溫度近似呈正弦曲線變化,且較環(huán)境溫度變化滯后約2 h。
2)太陽輻射、環(huán)境溫度是影響安全殼變形的主要環(huán)境因素,應在安全殼結構整體性試驗時予以考慮,以對變形測量結果進行合理的修正。當忽略太陽輻射時,穹頂豎向位移平均幅值降低了約49%,筒身水平徑向位移平均幅值降低了約21%。
3)輻射換熱對安全殼結構變形影響有限,是影響安全殼結構變形的次要環(huán)境因素,為提高模擬精度,建議考慮輻射換熱的影響。
4)兩種熱邊界條件施加方法對安全殼變形的影響可忽略,為了計算簡便,可采用綜合系數(shù)法進行施加。