汪 振 王 剛 辜峙钘 柏云清 龍鵬程

1（中國科學(xué)院核能安全技術(shù)研究所 中子輸運(yùn)理論與輻射安全重點(diǎn)實(shí)驗(yàn)室 合肥 230031）2（中國科學(xué)技術(shù)大學(xué) 合肥 230026）

加速器驅(qū)動鉛鉍冷卻自然循環(huán)次臨界堆束流瞬變事故研究

汪 振1,2王 剛1辜峙钘1,2柏云清1龍鵬程1

1（中國科學(xué)院核能安全技術(shù)研究所 中子輸運(yùn)理論與輻射安全重點(diǎn)實(shí)驗(yàn)室 合肥 230031）2（中國科學(xué)技術(shù)大學(xué) 合肥 230026）

利用FDS團(tuán)隊(duì)(Fission & Fusion Design Study)開發(fā)的中子學(xué)與熱工水力學(xué)耦合安全分析軟件，對一種加速器驅(qū)動鉛鉍自然循環(huán)次臨界反應(yīng)堆的束流中斷及束流超功率事故進(jìn)行了模擬分析。計(jì)算結(jié)果表明：加速器驅(qū)動次臨界潔凈核能系統(tǒng)(Accelerator Driven Sub-critical System, ADS)次臨界堆的功率對束流瞬變的響應(yīng)幾乎是瞬時的；事故工況下，自然循環(huán)會根據(jù)堆芯功率自動調(diào)整至重新達(dá)到穩(wěn)定；失束時間越長，材料溫度降得越低，功率瞬間恢復(fù)值越低，束流恢復(fù)后，材料溫度回升的速度越快；束流200%超功率事故發(fā)生后堆芯功率最終穩(wěn)定在初始功率值的192.2%，燃料溫度增幅最大，為286K，燃料和包殼不會發(fā)生損壞和熔化，冷卻劑不會發(fā)生沸騰。

加速器驅(qū)動次臨界潔凈核能系統(tǒng)，失束，束流超功率，中子學(xué)與熱工水力學(xué)耦合

“分離-嬗變”是當(dāng)今世界公認(rèn)的處理核廢物的有效策略。加速器驅(qū)動次臨界潔凈核能系統(tǒng)(Accelerator Driven Sub-critical System, ADS)被認(rèn)為是嬗變核廢料最有效的核裝置之一，目前國內(nèi)外均開展了大量的研究[1–3]。2011年，中國科學(xué)院正式啟動了戰(zhàn)略性先導(dǎo)科技專項(xiàng)“未來先進(jìn)核裂變能——ADS嬗變系統(tǒng)”研究計(jì)劃，核能安全技術(shù)研究所FDS團(tuán)隊(duì)承擔(dān)了該專項(xiàng)鉛鉍反應(yīng)堆項(xiàng)目的設(shè)計(jì)研究工作，積累了鉛鉍相關(guān)技術(shù)經(jīng)驗(yàn)[4–7]。

ADS是以加速器產(chǎn)生的高能強(qiáng)流質(zhì)子束轟擊靶核產(chǎn)生散裂中子作為外源中子驅(qū)動和維持次臨界堆運(yùn)行。因此，質(zhì)子束流的不穩(wěn)定性將對次臨界堆的功率水平產(chǎn)生影響，進(jìn)而對ADS的安全造成威脅。高能質(zhì)子加速器頻繁發(fā)生失束事故和可能存在束流超功率是兩種典型的ADS束流瞬變事故[8–9]。失束會引發(fā)外源的突然中斷，由此引起的功率驟變，直接影響反應(yīng)堆材料的溫度變化，引起熱應(yīng)力變化，導(dǎo)致材料疲勞與蠕變進(jìn)而影響反應(yīng)堆的壽命。加速器的束流超功率，將會導(dǎo)致次臨界堆芯功率急劇上升，進(jìn)而使堆芯溫度發(fā)生劇增，可能會超過材料的熔點(diǎn)或沸點(diǎn)，威脅反應(yīng)堆的完整性。另外，在這兩種事故下，功率的變化也將對堆芯的自然循環(huán)產(chǎn)生影響，堆芯能否重新達(dá)到穩(wěn)定有待考查。

本文利用中子學(xué)與熱工水力學(xué)耦合安全分析軟件NTC程序，對FDS團(tuán)隊(duì)設(shè)計(jì)的一種熱功率為10MW依靠自然循環(huán)驅(qū)動鉛鉍冷卻ADS的失束事故及束流超功率事故進(jìn)行模擬，關(guān)注在這兩種事故下反應(yīng)堆的瞬態(tài)特性。

1 鉛鉍冷卻ADS計(jì)算模型

1.1 堆芯主要設(shè)計(jì)參數(shù)

本文模擬所用的鉛鉍冷卻ADS主要設(shè)計(jì)參數(shù)如表1所示，其堆芯采用一體化池式結(jié)構(gòu)，依靠自然循環(huán)的液態(tài)鉛鉍進(jìn)行冷卻；靶區(qū)采用有窗靶，依靠強(qiáng)迫循環(huán)的液態(tài)鉛鉍進(jìn)行冷卻。燃料包殼采用316Ti不銹鋼，組件壁采用316L不銹鋼。堆芯依據(jù)功能劃分由內(nèi)到外依次為：靶區(qū)、活性區(qū)、反射層和屏蔽層。燃料組件內(nèi)含有61根燃料棒，燃料棒呈三角形點(diǎn)陣排列，燃料棒間的中心間距為16.74mm。燃料棒外直徑為15mm，其包殼厚度為0.7mm，有效長度為800mm。屏蔽層組件由7根外直徑為23 mm的含天然碳化硼的不銹鋼棒組成，有效長度為1100mm。

表1 ADS的主要設(shè)計(jì)參數(shù)Table 1 Main data for ADS.

1.2 計(jì)算模型

FDS團(tuán)隊(duì)自主開發(fā)了中子學(xué)與熱工水力學(xué)耦合的安全分析程序NTC，并應(yīng)用于聚變驅(qū)動次臨界系統(tǒng)及聚變包層的瞬態(tài)分析[10–19]。其二維版本NTC-2D為R-Z兩維，中子學(xué)部分采用11群能量劃分、離散縱標(biāo)的準(zhǔn)靜態(tài)方法求解帶外源的中子輸運(yùn)方程，將中子通量密度函數(shù)分解為一個形狀函數(shù)和一個幅函數(shù)的乘積，可以模擬中子通量密度空間分布隨時間的變化，因此適用于ADS次臨界反應(yīng)堆事故模擬。

建模時需把三維堆芯實(shí)體幾何模型進(jìn)行簡化，等效成二維，即在徑向上把實(shí)際的堆芯組件布置按照面積相等等效成環(huán)形，將每一環(huán)的寬度作為每一個徑向流體力學(xué)網(wǎng)格的寬度，如圖1所示，共14環(huán)。軸向劃分了39個流體力學(xué)網(wǎng)格。在流體力學(xué)網(wǎng)格下再細(xì)分中子學(xué)網(wǎng)格，一般為0.5個擴(kuò)散長度即可。建模采用的換熱器模型為理想換熱器，其材料密度設(shè)置為無窮大，使任何溫度的流體經(jīng)過換熱器都能降到入口溫度，以便有效地模擬自然循環(huán)；靶區(qū)為強(qiáng)迫循環(huán)冷卻，不模擬質(zhì)子與鉛鉍的散裂反應(yīng)，而是由FLUKA軟件給出散裂反應(yīng)產(chǎn)生的外源中子源強(qiáng)的空間分布與能譜分布。

圖1 用于NTC計(jì)算的鉛鉍冷卻的ADS二維模型示意圖Fig.1 Two-dimensional model of lead-bismuth cooled ADS.

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

2.1 穩(wěn)態(tài)工況計(jì)算結(jié)果及分析

在進(jìn)行束流瞬變事故模擬之前，首先需要對鉛鉍冷卻ADS進(jìn)行穩(wěn)態(tài)工況的計(jì)算。反應(yīng)堆內(nèi)自持鏈?zhǔn)搅炎兎磻?yīng)的條件可以很方便地用有效增殖系數(shù)keff來表示，從中子平衡關(guān)系上keff定義為系統(tǒng)內(nèi)中子的產(chǎn)生率與系統(tǒng)內(nèi)中子的總消失（吸收+泄露）率的比值。ADS的有效增殖系數(shù)keff會隨著燃耗的加深產(chǎn)生波動，同時反應(yīng)堆的重要安全參數(shù)，例如多普勒系數(shù)、冷卻劑溫度系數(shù)等也會發(fā)生改變[20]。本次模擬選在壽期初，即keff最大時。表2給出了反應(yīng)堆在穩(wěn)態(tài)運(yùn)行時的計(jì)算結(jié)果。由表2可知，每項(xiàng)參數(shù)的相對偏差都在1.06%以內(nèi)，與設(shè)計(jì)值符合得很好。

表2 穩(wěn)態(tài)計(jì)算結(jié)果Table 2 Steady state calculation results.

圖2給出了堆芯最熱通道內(nèi)冷卻劑溫度、燃料溫度、包殼溫度的軸向分布。結(jié)果顯示，燃料溫度的最大值出現(xiàn)在軸向中心偏上，這與實(shí)際情況是相符的。因?yàn)槿剂蠝囟戎饕芏研井a(chǎn)熱和冷卻劑溫度兩方面影響，冷卻劑在軸向上一直被加熱，所以冷卻劑溫度隨軸向流動方向升高，在出口處達(dá)到最大。從中子學(xué)角度講，在軸向上堆芯的中平面處中子通量密度最大，核熱最多，受兩方面影響，燃料溫度的最大值出現(xiàn)在中平面和堆芯出口之間。

圖2 堆芯最熱通道內(nèi)燃料、包殼和冷卻劑溫度的軸向分布Fig.2 Axial temperature distributions of the fuel, cladding and coolant in the hottest channel.

圖3 給出了4個燃料組件通道軸向的歸一化功率分布曲線。由圖3可以看出，在軸向上，所有通道都呈現(xiàn)出中間大兩邊小的類似余弦分布，其中軸向網(wǎng)格14–23為活性區(qū)。另外，還可以看到4個燃料組件越靠近堆芯中心，功率越大，這與中子通量密度徑向分布有關(guān)。

圖3 堆芯4個燃料通道軸向功率分布Fig.3 Axial power distributions of the four fuel channels.

綜上所述，反應(yīng)堆穩(wěn)態(tài)工況計(jì)算得到的各參數(shù)都與設(shè)計(jì)值或?qū)嶋H情況符合得較好，可以作為束流瞬變事故的初始狀態(tài)。

2.2 失束事故計(jì)算結(jié)果及分析

在穩(wěn)態(tài)的基礎(chǔ)上進(jìn)行失束事故的模擬，由于現(xiàn)有加速器離子源及加速結(jié)構(gòu)的射頻系統(tǒng)事故會導(dǎo)致短時間的失束，參考經(jīng)濟(jì)合作與發(fā)展組織(Organization for Economic Co-operation and Development, OECD)/核能署(Nuclear Energy Agency, NEA)的研究成果[21]，束流中斷時間分別選為1 s、3 s、6 s、12 s。考查在4種情況下堆芯的功率和溫度響應(yīng)情況。

圖4給出了不同失束時間下堆芯功率隨時間的變化曲線。反應(yīng)堆在0 s時開始失束，堆芯功率瞬間降到初始功率的19.7%，失束后反應(yīng)堆處于次臨界狀態(tài)，功率繼續(xù)衰減。當(dāng)束流恢復(fù)時，功率同樣瞬間升到接近初始功率的水平，可見堆芯功率對外源的響應(yīng)是瞬時的，此時主要是瞬發(fā)中子在起作用。束流恢復(fù)后，由于緩發(fā)中子的慢響應(yīng)，使功率緩慢增至初始水平。

圖4 不同失束時間(1 s、3 s、6 s、12 s)下功率的變化Fig.4 Core power at different beam interruption time (1 s, 3 s, 6 s, 12 s).

比較4種不同失束時間的計(jì)算結(jié)果，可以看出，失束時間越長，功率瞬間恢復(fù)值越低，這是因?yàn)槭鴷r間越長，緩發(fā)中子先驅(qū)核積累的越多，束流恢復(fù)時，緩發(fā)中子參與鏈?zhǔn)搅炎兎磻?yīng)的相對份額就會增加，因此由瞬發(fā)中子導(dǎo)致的功率驟增就會減小。

圖5和圖6分別給出不同失束時間下燃料最熱處的燃料芯部溫度和包殼溫度隨時間的變化。從圖5中可以看出，由于失束引起的功率驟降，燃料芯部溫度和包殼溫度會迅速降低，失束時間越長，溫度降得越低；而當(dāng)束流恢復(fù)后，燃料溫度和包殼溫度會迅速回升，失束時間越長，束流恢復(fù)后溫度回升的速率越快。這是因?yàn)椋剂系臏囟戎饕扇剂系漠a(chǎn)熱（功率）和燃料的散熱（與冷卻劑等的傳熱）決定，束流恢復(fù)后，功率的恢復(fù)值基本一樣，也就是說對于不同失束時間束流恢復(fù)后的燃料產(chǎn)熱是一樣的。圖7給出不同失束時間下燃料與冷卻劑溫差的變化曲線，可以看出失束時間越長，燃料與冷卻劑的溫差會越小，根據(jù)牛頓冷卻公式粗略的估計(jì)，燃料向外散熱的熱流密度就會越小。因此與產(chǎn)熱一樣，散熱慢的溫度回升得快。對于包殼溫度也是相同的道理。

圖5 不同失束時間(1 s、3 s、6 s和12 s)下燃料最熱處燃料芯部最大溫度的變化Fig.5 Max fuel temperature at different beam interruption time (1 s, 3 s, 6 s, 12 s).

圖6 不同失束時間(1s、3s、6s、12s)包殼溫度變化Fig.6 Cladding temperature at different beam interruption time (1s, 3s, 6s, 12s).

圖7 不同失束時間(1 s、3 s、6 s、12 s)燃料與冷卻劑溫差變化Fig.7 Temperature difference between fuel interior and coolant at different beam interruption time (1s, 3s, 6s, 12s).

2.3 束流超功率事故計(jì)算結(jié)果及分析

在穩(wěn)態(tài)計(jì)算的基礎(chǔ)上，0 s時事故開始，束流超功率瞬間使外源中子增加一倍，反應(yīng)堆不停堆且無其他動作。

圖8給出了事故過程堆芯功率和反應(yīng)性的變化曲線。由于束流功率加倍，堆芯功率瞬間上升至約188%初始功率值，這主要是因?yàn)樗舶l(fā)中子的快響應(yīng)。隨著時間的進(jìn)行，緩發(fā)中子逐漸積累，開始參與鏈?zhǔn)搅炎兎磻?yīng)，導(dǎo)致功率在突然增長之后又有一個緩慢增加的過程，最終堆芯功率在約100 s時重新穩(wěn)定，此時的堆芯功率凈增值為92.2%的穩(wěn)態(tài)功率值。由于負(fù)反饋，堆芯總反應(yīng)性隨溫度的上升而減小。

圖8 束流超功率過程功率及反應(yīng)性隨時間的變化Fig.8 Core power and reactivity change with time under the beam overpower.

圖9 給出束流超功率過程堆芯最熱通道處燃料、包殼和出口冷卻劑溫度的變化曲線。由于束流超功率導(dǎo)致堆芯功率迅速增加，引起相關(guān)的材料溫度上升。其中燃料芯部溫度漲幅最大，增加了約286K，穩(wěn)定值約為1160K，但低于燃料的溫度限值2573 K，因此燃料不會發(fā)生熔化或損壞。包殼溫度增加了約110 K，穩(wěn)定值約為825 K，也低于包殼的溫度限值973 K，因此包殼也不會發(fā)生損壞。冷卻劑溫度增加近108 K，達(dá)到約817 K，但仍遠(yuǎn)低于液態(tài)鉛鉍的沸點(diǎn)1 943 K，冷卻劑不會發(fā)生沸騰。

圖9 最熱通道出口處冷卻劑、包殼、燃料表面、燃料芯部溫度隨時間變化Fig.9 Temperature change with time of fuel, cladding and coolant in the hottest channel.

圖10給出了束流超功率過程中堆芯流量的變化。由于反應(yīng)堆采用的是自然循環(huán)，當(dāng)堆芯功率增加時，堆芯冷卻劑溫度增加，致使自然循環(huán)的驅(qū)動壓頭增大，冷卻劑流量增大?？梢詫ι厦嫒剂蠝囟仍龇淮笞龀鼋忉?。

圖10 堆芯的質(zhì)量流率隨時間變化Fig.10 Mass flow rate change with time in core.

3 結(jié)語

從上述計(jì)算結(jié)果和分析可以得出：(1) ADS中功率對外源變化的響應(yīng)幾乎是瞬時的；(2) 失束時間越長，材料溫度降得越低，功率瞬間恢復(fù)值越低，束流恢復(fù)后，材料溫度回升的速度越快；(3) 200%束流超功率事故發(fā)生后，最終功率為初始功率的192.2%，在此過程中燃料、包殼和冷卻劑溫度均低于安全限值，表明燃料和包殼都不會發(fā)生損壞，冷卻劑不會發(fā)生沸騰；(4) 在兩種事故下，堆芯溫度重新達(dá)到穩(wěn)定值，表明自然循環(huán)能夠在功率變化時重新達(dá)到穩(wěn)定。

致謝 本工作得到中國科學(xué)院核能安全技術(shù)研究所FDS團(tuán)隊(duì)其他成員的指導(dǎo)和幫助，在此向他們表示衷心的感謝。

1 Sugawara T, Nishihara K, Tsujimoto K. Transient analysis for lead-bismuth cooled accelerator-driven system[J]. Annals of Nuclear Energy, 2013, 55: 238–247

2 Wu Y, Bai Y, Wang W, et al. Overview of China lead alloy cooled reactor development and ADS program in China[C]. NUTHOS-9, Kaohsiung, Taiwan, September 9–13, 2012

3 Mikityuk K, Coddington P, Bubelis E. Comparison of the transient behavior of LBE- and gas-cooled experimental accelerator driven systems[J]. Nuclear Engineering and Design, 2006, 236: 2452–2473

4 Wu Y C, Bai Y Q, Song Y, et al. Overview of lead-based reactor design and R&D status in China[C]. International Conference on Fast Reactors and Related Fuel Cycles: Safe Technologies and Sustainable Scenarios(FR13), Paris, France, March 4–7, 2013

5 吳宜燦, 黃群英, 柏云清, 等. 液態(tài)鉛鉍回路設(shè)計(jì)研制與材料腐蝕性實(shí)驗(yàn)初步研究[J]. 核科學(xué)與工程, 2010, 30(3): 238–243

WU Yican, HUANG Qunying, BAI Yunqing, et al. Preliminary experimental study on the corrosion of structural steels in liquid lead bismuth loop[J]. Nuclear Science and Engineering, 2010, 30(3): 238–243

6 吳宜燦, 柏云清, 宋勇, 等. 中國鉛基研究反應(yīng)堆概念設(shè)計(jì)研究[J]. 核科學(xué)與工程, 2014, 34(2): 201–208

WU Yican, BAI Yunqing, SONG Yong, et al. Study on the conceptual design for China lead-bismuth research reactor[J]. Nuclear Science and Engineering, 2014, 34(2): 201–208

7 王改英, 柏云清, 高勝, 等. 液態(tài)鉛鉍氧濃度測量技術(shù)初步研究[J]. 核科學(xué)與工程, 2012, 32(2): 165–169

WANG Gaiying, BAI Yunqing, GAO Sheng, et al. Preliminary study on the measurement technology of oxygen concentration in liquid lead bismuth[J]. Nuclear Science and Engineering, 2012, 32(2): 165–169

8 Liu P, Chen X N, Rineiski A, et al. Transient analysis of the 400 MWth-class EFIT accelerator driven transmuter with the multi-physics code: SIMMER-III[J]. Nuclear Engineering and Design, 2010, 240(10): 3481–3494

9 于濤, 李吉根, 凌球, 等. ADS加速器束流瞬變分析程序開發(fā)[J]. 核動力工程, 2007, 28(2): 124–127

YU Tao, LI Jigen, LING Qiu, et al. Development of a beam transient code for ADS[J]. Nuclear Power Engineering, 2007, 28(2): 124–127

10 Bai Y Q, Ke Y, Wu Y C. Preliminary analysis of typical transients in fusion driven subcritical system (FDS-I)[C]. 15thInternational Conference on Nuclear Engineering, Nagoya, Japan, April 22–26, 2007

11 Wu Y C, FDS Team. Design status and development strategy of China liquid lithium-lead blankets and related material technology[J]. Journal of Nuclear Materials, 2007, 367–370: 1410–1415

12 Wu Y C, FDS Team. Conceptual design and testing strategy of a dual functional lithium-lead test blanket module in ITER and EAST[J]. Nuclear Fusion, 2007, 47(11): 1533–1539

13 Wu Y C, FDS Team. Conceptual design activities of FDS series fusion power plants in China[J]. Fusion Engineering and Design, 2006, 81(23–24): 2713–2718

14 Wu Y C, FDS Team. Design analysis of the China Dual-functional Lithium Lead (DFLL) test blanket module in ITER[J]. Fusion Engineering and Design, 2007, 82: 1893–1903

15 Wu Y C. Progress in fusion-driven hybrid system studies in China[J]. Fusion Engineering and Design, 2002, 63–64: 73–80

16 Wu Y C, Zhu X X, Zheng S L, et al. Neutronics analysis of dual-cooled waste transmutation blanket for the FDS[J]. Fusion Engineering and Design, 2002, 63–64: 133–138

17 Wu Y C, FDS Team. Overview of liquid lithium lead breeder blanket program in China[J]. Fusion Engineering and Design, 2011, 86(9–11): 2343–2346

18 柏云清. 磁約束聚變堆包層瞬態(tài)安全特性分析研究[D].中國科學(xué)院合肥物質(zhì)科學(xué)研究院, 2007

BAI Yunqing. Transient safety characteristic for magnetic confinement fusion reactor blanket[D]. Hefei Institutes of Physical Science, Chinese Academy of Sciences, 2007

19 柯嚴(yán). 聚變驅(qū)動次臨界堆瞬態(tài)安全分析研究[D]. 中國科學(xué)院合肥物質(zhì)科學(xué)研究院, 2006

KE Yan. Transient safety analysis for fusion-driven subcritical reactor[D]. Hefei Institutes of Physical Science, Chinese Academy of Sciences, 2006

20 Schikorr W M. Assessments of the kinetic and dynamic transient behavior of sub-critical systems (ADS) in comparison to critical reactor systems[J]. Nuclear Engineering and Design, 2001, 210: 95–123

21 D’Angelo A, Gabrielli F. Benchmark on beam interruptions in an accelerator-driven system-final report on Phase I calculations[R]. OECD/NEA, 2003

CLC TL364

Beam transient accident for lead-bismuth cooled ADS by natural circulation

WANG Zhen1,2WANG Gang1GU Zhixing1,2BAI Yunqing1LONG Pengcheng1
1(Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei 230031, China) 2(University of Science and Technology of China, Hefei 230026, China)

Background: Beam interruption and beam overpower are two typical transient accidents for accelerator-driven system. Purpose: To investigate the safety characteristics of the lead-bismuth cooled Accelerator Driven Sub-critical System (ADS) under beam transient accident, the steady state, beam trip and beam overpower accident of lead-bismuth cooled accelerator-driven system were simulated by using NTC-2D. Methods: NTC-2D is a two-dimensional version of neutronics and thermo-hydraulics coupled simulation program NTC developed by FDS team. Results: As for the beam interruptions, the temperature variations of cladding and fuel pellet at different time (1s, 3s, 6s, 12s) were given. The longer the beam interruptions, the lower temperature of cladding and fuel. After beam overpower occurred, core power finally stabilized at 192.2% of the initial power value. The temperatures of the fuel, cladding and coolant were all smaller than the safety limits. The fuel and cladding would not melt and the coolant would not boil. Conclusion: The reactor was safe under two accidents. The transient response of the power for ADS under beam transient accidents is instantaneous. Delayed neutron has less effect on neutron generation time. Natural circulation can reach new steady under two accidents.

Accelerator Driven Sub-critical System (ADS), Beam trip, Beam overpower, Neutronics and thermal-hydraulics coupled

TL364

10.11889/j.0253-3219.2015.hjs.38.010604

中科院戰(zhàn)略性先導(dǎo)科技專項(xiàng)(No.XDA03040000)和國家自然科學(xué)基金項(xiàng)目(No.91026004)資助

汪振，男，1989年出生，2012年畢業(yè)于南華大學(xué)，現(xiàn)為博士研究生，研究領(lǐng)域?yàn)榉磻?yīng)堆事故安全分析

王剛，E-mail: gang.wang@fds.org.cn

2014-06-26，

2014-09-13