小型自然循環(huán)鉛鉍快堆無保護(hù)超功率事故不確定性分析研究

余清遠(yuǎn) 漆靜雯 趙鵬程 趙亞楠 于 濤

（南華大學(xué)核科學(xué)技術(shù)學(xué)院 衡陽421001）

采用液態(tài)鉛鉍冷卻的鉛鉍堆在嬗變核廢料和增殖核燃料方面具有獨(dú)特優(yōu)勢，擁有良好的經(jīng)濟(jì)性和固有安全性，被第四代國際核能論壇預(yù)計(jì)有望成為首個(gè)實(shí)現(xiàn)工業(yè)示范化的第四代核能系統(tǒng)［1］。事故瞬態(tài)是鉛鉍快堆熱工安全分析的重要內(nèi)容，有必要探究鉛鉍堆瞬態(tài)過程中關(guān)鍵參數(shù)的變化規(guī)律，評(píng)估鉛鉍堆的固有安全性。由于設(shè)計(jì)參數(shù)、初始邊界條件等因素的影響，鉛鉍堆瞬態(tài)模擬計(jì)算結(jié)果不可避免地存在不確定性，將不確定性分析方法應(yīng)用鉛鉍堆瞬態(tài)事故分析，可量化不確定性分析結(jié)果，評(píng)估瞬態(tài)工況中輸入?yún)?shù)對(duì)鉛鉍堆安全的影響，對(duì)鉛鉍堆的設(shè)計(jì)分析、安全審評(píng)具有現(xiàn)實(shí)意義。

國際上，國際性合作研究項(xiàng)目OECD（Organization for Economic Co-operation and Development）/NEA LWR UAM（Nuclear Energy Agency，Light Water Reactor，Uncertainty Analysis in Modeling）正在全力開展有關(guān)多物理場、多尺度耦合的不確定性分析［2-5］；日本基于京都大學(xué)臨界裝置開展了一系列加速器驅(qū)動(dòng)系統(tǒng)實(shí)驗(yàn)，并采用隨機(jī)采樣技術(shù)對(duì)實(shí)驗(yàn)開展了不確定性量化分析［6-9］；瑞典皇家理工學(xué)院基于ROAAM+、MELCOR 等程序?qū)Ρ睔W沸水堆開展了不確定性量化和敏感性分析［10-13］；由美國阿貢實(shí)驗(yàn)室牽頭的先進(jìn)反應(yīng)堆E級(jí)高性能模擬中心開展了堆芯物理計(jì)算的不確定量化分析和敏感性分析［14-15］。近年來，國內(nèi)有關(guān)快堆的不確定性研究主要集中于不確定性分析方法：西安交通大學(xué)先后采用不確定性分析方法開展了鈉冷快堆事故SHRT-17 試驗(yàn)數(shù)據(jù)和BN-600 快堆有效增值因數(shù)的不確定性和敏感性分析［16-17］；哈爾濱工程大學(xué)基于抽樣統(tǒng)計(jì)法開展了中國實(shí)驗(yàn)快堆（China Experimental Fast Reactor，CEFR）的不確定性分析和敏感性分析［18］；中國原子能研究院采用多種不確定性分析方法開展了CEFR堆芯設(shè)計(jì)計(jì)算的不確定度定量評(píng)價(jià)［19］。中國科學(xué)技術(shù)大學(xué)針對(duì)鉛鉍快堆開發(fā)了多物理多尺度程序Fluent/KMC-sub/NDK［20］，并對(duì)SNCLFR-100開展了啟堆條件下的流動(dòng)不穩(wěn)定性分析［21］。而鉛鉍快堆作為具有發(fā)展?jié)摿Φ暮蜻x堆型之一，其瞬態(tài)安全特性對(duì)反應(yīng)堆運(yùn)行至關(guān)重要，開展鉛鉍快堆的瞬態(tài)事故的不確定性量化分析，研究瞬態(tài)事故的不確定性來源，對(duì)鉛鉍快堆的設(shè)計(jì)、安全特性評(píng)估有重要意義。

本工作是針對(duì)中國科學(xué)技術(shù)大學(xué)和南華大學(xué)聯(lián)合設(shè)計(jì)的小型自然循環(huán)鉛基鉛鉍快堆（Small Natural Circulation Lead-bismuth Cooled Fast Reactor 10 MWth，SNCLFR-10），首先，通過文獻(xiàn)資料以及工程經(jīng)驗(yàn)判斷，考慮三類不確定性輸入?yún)?shù)：熱工水力模型參數(shù)、中子物理模型參數(shù)以及燃料制造公差；然后，基于課題組自主開發(fā)的物理/熱工耦合不確定性量化分析程序CFD/PFS 程序開展無保護(hù)超功率事故的不確定性分析，采用拉丁超立方（Latin Hypercube Sampling，LHS）對(duì)輸入?yún)?shù)進(jìn)行抽樣，對(duì)瞬態(tài)安全參數(shù)進(jìn)行不確定性量化分析，并通過相關(guān)系數(shù)分析法評(píng)估輸入?yún)?shù)對(duì)瞬態(tài)安全參數(shù)的敏感度。

1 計(jì)算模型

1.1 分析模型

SNCLFR-10 堆是由中國科學(xué)技術(shù)大學(xué)和南華大學(xué)聯(lián)合設(shè)計(jì)的小型自然循環(huán)鉛鉍快堆，具有良好的現(xiàn)實(shí)可行性、安全可靠性、實(shí)驗(yàn)靈活性和技術(shù)延續(xù)性［22］。圖1 為SNCLFR-10 堆的結(jié)構(gòu)設(shè)計(jì)圖，表1、2給出了反應(yīng)堆設(shè)計(jì)參數(shù)及堆芯主要設(shè)計(jì)參數(shù)。SNCLFR-10堆的熱功率為10 MW，燃料采用富集度為19.75%的UO2，燃料包殼材料采用316Ti不銹鋼。堆芯活性區(qū)由94盒組件組成，其中86盒燃料組件，8盒控制棒組件。一回路冷卻劑系統(tǒng)采用自然循環(huán)驅(qū)動(dòng)液態(tài)鉛鉍冷卻，冷卻劑入口溫度為533.15 K，出口溫度為663.15 K；二回路冷卻劑采用4 MPa 的加壓液態(tài)水。

表1 SNCLFR-10主要設(shè)計(jì)參數(shù)Table 1 Main design parameters of SNCLFR-10

圖1 SNCLFR-10堆結(jié)構(gòu)示意圖Fig.1 Structure diagram of SNCLFR-10

無保護(hù)超功率事故（Unprotected Transient Overpower，UTOP）是指向反應(yīng)堆內(nèi)突然引入一個(gè)正反應(yīng)性，導(dǎo)致功率急劇升高的瞬態(tài)過程，此過程中反應(yīng)堆停堆保護(hù)系統(tǒng)失效。UTOP 事故的瞬態(tài)條件：反應(yīng)堆初始工況為穩(wěn)態(tài)運(yùn)行，假設(shè)最大價(jià)值控制棒失去控制而不可抽出，故在2 s 內(nèi)迅速引入1$正反應(yīng)性，整個(gè)過程無停堆保護(hù)動(dòng)作。

1.2 CFD/PFS程序

CFD/PFS 是基于CFD 程序FLUENT 的用戶自定義函數(shù)（User Defined Function，UDF）耦合點(diǎn)堆中子動(dòng)力學(xué)模型（Point Kinetics Model，PKM）、燃料棒熱傳導(dǎo)模型（Fuel Pin Heat Transfer Model，F(xiàn)PH）以及不確定性分析程序SIMLAB開發(fā)的一套物理熱工耦合不確定性分析程序，可用于分析小型鉛鉍快堆的瞬態(tài)事故。SIMLAB 是研究不確定性的開源程序，由歐洲聯(lián)合研究中心主持研究開發(fā)，配備可視化操作界面，其中含有統(tǒng)計(jì)預(yù)處理子系統(tǒng)、模型執(zhí)行子系統(tǒng)、統(tǒng)計(jì)后處理子系統(tǒng)，廣泛用于不確定性和整體敏感性研究［23］。采用與SIMMER 程序code-to-code的方法進(jìn)行CFD/PFS 的對(duì)比驗(yàn)證，CFD/PFS 程序開發(fā)與詳細(xì)驗(yàn)證見參考文獻(xiàn)［24］。

CFD/PFS程序分為物理熱工耦合計(jì)算模型與不確定性分析計(jì)算，其計(jì)算流程如圖2 所示。物理熱工耦合計(jì)算模型包括CFD模型、點(diǎn)堆中子動(dòng)力學(xué)模型、燃料棒熱傳導(dǎo)模型。PKM是三維時(shí)空動(dòng)力學(xué)的簡化模型，采用6 組緩發(fā)中子的點(diǎn)堆動(dòng)力學(xué)方程主要用于計(jì)算反應(yīng)性變化：

圖2 CFD/PFS數(shù)據(jù)傳遞流程Fig.2 The data transfer route of CFD/PFS

式中：n（t）是中子密度；ρ（t）是反應(yīng)性；β是緩發(fā)中子總份額；βi是i組緩發(fā)中子的份額；Λ是中子代時(shí)間；λi是第i組緩發(fā)中子的衰變常數(shù)；Ci是i組緩發(fā)中子先驅(qū)核的濃度。

再根據(jù)反應(yīng)性通過點(diǎn)堆動(dòng)力學(xué)方程求解相對(duì)中子密度n（t），則總功率為：

反應(yīng)性反饋源于反應(yīng)堆溫度、壓力或者流量的變化，而溫度對(duì)反應(yīng)性的影響是主要反饋，故所有反應(yīng)性反饋、堆內(nèi)溫度總反應(yīng)性公式可表示為：

式中：KD為多普勒系數(shù)；αc，a為堆芯軸向膨脹反饋系數(shù)；αc，γ為堆芯徑向膨脹反饋系數(shù)；αc為冷卻劑溫度反饋系數(shù)；Tf（t）和Tf（0）分別為t時(shí)刻的燃料、冷卻劑平均溫度。

FTM包括燃料棒穩(wěn)態(tài)和瞬態(tài)熱傳導(dǎo)模型，主要用于計(jì)算燃料棒的溫度分布，包括燃料芯塊、氣隙以及包殼的溫度，提供熱流密度給CFD 模型，反饋芯塊平均溫度給PKM；CFD模型主要用于計(jì)算堆芯熱工水力參數(shù)，包括冷卻劑的溫度、速度等，反饋冷卻劑平均溫度給PKM，提供冷卻劑溫度和速度給FTM。不確定性分析計(jì)算流程如下：1）確定重要輸入?yún)?shù)及其分布范圍；2）利用SimLab程序?qū)Σ淮_定性輸入?yún)?shù)抽樣，生成輸入工況；3）將輸入工況傳遞給物理熱工耦合計(jì)算模型進(jìn)行計(jì)算，生成輸出工況；4）將輸出工況傳遞給SimLab程序，對(duì)結(jié)果進(jìn)行不確 定性分析和敏感性分析。

表2 SNCLFR-10堆芯主要設(shè)計(jì)參數(shù)Table 2 Core design parameters of SNCLFR-10

2 不確定性方法及關(guān)鍵參數(shù)識(shí)別

2.1 不確定性分析方法

本工作基于拉丁超立方抽樣的非參數(shù)統(tǒng)計(jì)法進(jìn)行不確定性分析，拉丁超立方采取分層抽樣，LHS抽樣過程如下：1）根據(jù)輸入?yún)?shù)的數(shù)量及分布范圍等概率的分成n個(gè)區(qū)域；2）在每個(gè)區(qū)域內(nèi)進(jìn)行隨機(jī)抽樣；3）將抽樣產(chǎn)生的樣本進(jìn)行隨機(jī)組合；4）挑選滿足約束條件的組合樣本作為輸入?yún)?shù)樣本。與隨機(jī)抽樣的蒙特卡羅方法比較，拉丁超立方抽樣可用較少樣本完成不確定性量化，且在均值和方差的估計(jì)上有顯著改善。

非參數(shù)統(tǒng)計(jì)法（wilks 公式）根據(jù)分布區(qū)間的置信度β以及輸出參數(shù)的概率γ確定所需最小抽樣數(shù)目n。wilks公式如下：

其中，式（4）適用于單側(cè)限值分布，式（5）適用于雙側(cè)限值分布。

由wilks 公式計(jì)算出的滿足不同置信度和容許區(qū)間的最小樣本量n值如表3 所示。根據(jù)非參數(shù)統(tǒng)計(jì)理論，最小樣本量選取59 時(shí)，目標(biāo)參數(shù)的容忍上限值取最大值；最小樣本量選取93 時(shí)，目標(biāo)參數(shù)的容忍上限值取次大值［25］。

表3 不同置信度和容許區(qū)間下最小樣本信息分布Table 3 The minimum number of samples under different confidence and tolerance intervals

2.2 關(guān)鍵參數(shù)的選取

反應(yīng)堆瞬態(tài)計(jì)算中的不確定性輸入?yún)?shù)主要考慮三部分：1）物理模型輸入?yún)?shù)；2）熱工水力模型輸入?yún)?shù)；3）燃料制造公差輸入?yún)?shù)。針對(duì)SNCLFR-10堆的UTOP事故，根據(jù)始發(fā)事件選取清單、文獻(xiàn)資料以及工程經(jīng)驗(yàn)判斷等方式選取不確定性輸入?yún)?shù)，并確定其范圍和概率密度分布。選取的不確定輸入?yún)?shù)分為三個(gè)部分：首先是幾何輸入?yún)?shù)，一些參數(shù)測量有一定的不安全性故缺乏，只能由專家假設(shè)估算；其次是中子學(xué)參數(shù)，因?yàn)閷?shí)驗(yàn)測量系統(tǒng)及統(tǒng)計(jì)誤差，故其引入的不確定性需考慮；最后是熱工水力參數(shù)，程序、模型及經(jīng)驗(yàn)關(guān)系式等都會(huì)引入不確定性，故選取其中最典型5 個(gè)參數(shù)。表4 列出了25 個(gè)不確定輸入?yún)?shù)信息，輸入?yún)?shù)服從均勻分布或正態(tài)分布，其中不確定性度表示參數(shù)名義值（設(shè)計(jì)值）和參數(shù)上下限值的相對(duì)偏差。

表4 不確定性輸入?yún)?shù)信息Table 4 Uncertainty distribution of input parameter

由于測量儀器精準(zhǔn)、技術(shù)先進(jìn)，反應(yīng)堆部件制造精度高，因此燃料制造公差輸入?yún)?shù)的不確定度設(shè)為±0.5%。中子代時(shí)間、緩發(fā)中子份額、衰變常數(shù)以及初始功率的不確定度參考其他類似反應(yīng)堆的不確定度；對(duì)于反應(yīng)性反饋系數(shù)，考慮它包含了一系列不確定性（材料屬性、幾何形狀、截面數(shù)據(jù)以及蒙特卡羅計(jì)算），為保證其保守性，將反應(yīng)性反饋系數(shù)的不確定度設(shè)為10%［26］。熱工水力現(xiàn)象極其復(fù)雜，熱工物性參數(shù)大多是根據(jù)實(shí)驗(yàn)數(shù)據(jù)擬合以及經(jīng)驗(yàn)關(guān)系式所得，因此熱工水力模型輸入?yún)?shù)通常存在較大的不確定度［27］。

2.3 關(guān)鍵參數(shù)的抽樣

基于SIMLAB程序采用拉丁超立方方法實(shí)現(xiàn)不確定性輸入?yún)?shù)抽樣，種子數(shù)設(shè)定為98 765。根據(jù)wilks 非參數(shù)統(tǒng)計(jì)法，在95%置信水平、95%概率水平下，所需最小計(jì)算次數(shù)為93次。考慮到抽樣可能存在壞點(diǎn)、程序可能運(yùn)行失效等原因，抽樣數(shù)目設(shè)定為100 組，其中每一組工況包括25 個(gè)輸入?yún)?shù)。為直觀地觀測輸入?yún)?shù)的概率分布特征、檢測其抽樣的合理性，對(duì)不確定性輸入?yún)?shù)進(jìn)行歸一化處理。圖3為不確定性輸入?yún)?shù)的拉丁超立方抽樣樣本歸一化結(jié)果。

圖3 不確定性輸入?yún)?shù)抽樣歸一化結(jié)果Fig.3 Normalized result of input uncertainty parameter sampling

當(dāng)輸入?yún)?shù)分布類型為正態(tài)分布，抽樣樣本在分布范圍內(nèi)呈現(xiàn)為中間密集，兩邊稀疏；當(dāng)輸入?yún)?shù)分布類型為均勻分布，抽樣樣本在分布范圍內(nèi)呈均勻態(tài)勢。由圖3對(duì)比表4發(fā)現(xiàn)，選定各輸入?yún)?shù)的概率分布類型準(zhǔn)確無誤。

3 不確定性與敏感性分析

3.1 不確定性量化分析

針對(duì)SNCLFR-10 堆的UTOP 事故不確定性分析，通過不確定性輸入?yún)?shù)的組合量化反應(yīng)性、堆芯功率、燃料中心溫度以及包殼溫度的不確定性。圖6~9 為不同目標(biāo)參數(shù)的變化情況，目標(biāo)參數(shù)的名義值均處于95/95不確定包絡(luò)帶內(nèi)，不確定性分析結(jié)果具有合理性，且不確定度上限沒有出現(xiàn)過于保守的情況；對(duì)比不同目標(biāo)參數(shù)的不確定性包絡(luò)帶，可知不同時(shí)刻輸入?yún)?shù)的不確定性傳播到目標(biāo)參數(shù)的不確定性不同。

圖4、5 為100 組工況下反應(yīng)性與功率的變化情況，發(fā)生UTOP事故的2 s內(nèi)，總反應(yīng)性、總功率急劇增強(qiáng)。此階段不確定帶較小，隨后由于燃料和冷卻劑溫度上升，負(fù)反應(yīng)性反饋增強(qiáng)，且反應(yīng)性反饋系數(shù)不確定度較大，促使反應(yīng)性、功率的不確定帶相對(duì)變寬，隨著負(fù)反應(yīng)性反饋?zhàn)饔?，反?yīng)性與功率逐漸進(jìn)入一個(gè)新的穩(wěn)態(tài)。圖6、7為100組工況下包殼溫度、燃料溫度的變化情況，UTOP事故發(fā)生后，隨著功率的增大，燃料與包殼的溫度迅速上升，之后由于負(fù)反饋效應(yīng)增強(qiáng)，且反應(yīng)性反饋系數(shù)不確定度較大，燃料與包殼溫度的不確定帶愈加明顯，隨著功率下降，燃料 與包殼溫度逐漸進(jìn)入穩(wěn)態(tài)。

圖4 反應(yīng)性變化情況（100組工況）(a)與不確定性包絡(luò)帶(b)Fig.4 Changes of reactivity(100 sets of outputs)(a)and uncertainty bands(b)

圖5 功率變化情況（100組工況）(a)與不確定性包絡(luò)帶(b)Fig.5 Changes of power(100 sets of outputs)(a)and uncertainty bands(b)

圖6 包殼溫度變化情況（100組工況）(a)與不確定性包絡(luò)帶(b)Fig.6 Changes of cladding temperature(100 sets of outputs)(a)and uncertainty bands(b)

圖7 燃料溫度變化情況（100組工況）(a)與不確定性包絡(luò)帶(b)Fig.7 Changes of fuel temperature(100 sets of outputs)(a)and uncertainty bands(b)

基于非參數(shù)統(tǒng)計(jì)理論，在SNCLFR-10 堆的UTOP 事故中取100 組工況的次大值作為目標(biāo)參數(shù)容忍上限值。計(jì)算可得，包殼峰值溫度的95/95容忍上限值為1 208.66 K，燃料峰值溫度的95/95 容忍上限值為2 757.25 K。

3.2 敏感性分析

針對(duì)SNCLFR-10 堆UTOP 事故的敏感性分析，本工作基于SIMLAB 程序采用Spearman 秩相關(guān)分析方法研究輸入?yún)?shù)對(duì)燃料最高溫度以及包殼最高溫度的敏感度。Spearman 秩相關(guān)分析法用于全局敏感性分析，可分析多個(gè)輸入?yún)?shù)對(duì)目標(biāo)參數(shù)的影響。相關(guān)系數(shù)取值范圍為-1~+1，正號(hào)表示正相關(guān)，負(fù)號(hào)表示負(fù)相關(guān)，零表示不相關(guān)，相關(guān)系數(shù)絕對(duì)值越靠近1，敏感性越強(qiáng)。Spearman 秩相關(guān)系數(shù)表達(dá)式如下：

式中：ρs為Spearman秩相關(guān)系數(shù)；RXi為Xi在X中的大小排序；RYi為Yi在Y中的大小排序；n為樣本數(shù)。

圖8、9 分別為輸入?yún)?shù)與包殼最高溫度、燃料最高溫度的Spearman秩相關(guān)系數(shù)，深灰色柱狀條表示輸入?yún)?shù)與目標(biāo)參數(shù)呈正相關(guān)，淺灰色柱狀條表示輸入?yún)?shù)與目標(biāo)參數(shù)呈負(fù)相關(guān)。由圖8、9 可得，燃料制造公差輸入?yún)?shù)中，包殼內(nèi)直徑與燃料最高溫度、包殼最高溫度的相關(guān)系數(shù)分別為0.345 1、-0.352 2，呈中等相關(guān)；熱工水力輸入?yún)?shù)中，燃料平均溫度與燃料最高溫度、包殼最高溫度的相關(guān)系數(shù)分別為0.912、0.908 6，呈極強(qiáng)相關(guān)，這是由于燃料平均溫度自身存在較大的不確定度，且對(duì)反應(yīng)性反饋有直接影響；中子物理輸入?yún)?shù)中，第1組衰變常數(shù)與目標(biāo)參數(shù)的相關(guān)系數(shù)分別為-0.619 4、0.565 2，功率與目標(biāo)參數(shù)的相關(guān)系數(shù)分別為0.386、-0.241 1，均呈中度相關(guān)；燃料多普勒常數(shù)與目標(biāo)參數(shù)的的相關(guān)系數(shù)分別為0.934 5、0.945 6，徑向膨脹系數(shù)與目標(biāo)參數(shù)的的相關(guān)系數(shù)分別為0.703、0.742 2，軸向膨脹系數(shù)與目標(biāo)參數(shù)的的相關(guān)系數(shù)分別為0.477 99、0.554，說明反應(yīng)性反饋系數(shù)對(duì)目標(biāo)參數(shù)具有較高的靈敏度，這是由于其自身存在較大的不確定度，且反應(yīng)性反饋系數(shù)是引發(fā)UTOP工況的主要參數(shù)。

圖8 不確定性輸入?yún)?shù)與包殼最高溫度的Spearman秩相關(guān)系數(shù)Fig.8 Correlation coefficient of Spearman rank between the input parameters and peak cladding temperature

圖9 不確定性輸入?yún)?shù)與燃料最高溫度的Spearman秩相關(guān)系數(shù)Fig.9 Correlation coefficient of Spearman rank between the input parameters and peak fuel temperature

綜上所述，在SNCLFR-10 堆的UTOP 事故中，顯著影響燃料、包殼峰值溫度的輸入?yún)?shù)包括包殼直徑、穩(wěn)態(tài)燃料平均溫度、初始功率、第1 組衰變常數(shù)以及反應(yīng)性反饋系數(shù)；其中燃料平均溫度以及反應(yīng)性反饋系數(shù)對(duì)瞬態(tài)安全參數(shù)的影響最大。

4 結(jié)語

本工作基于CFD/PFS 程序開展了小型自然循環(huán)鉛鉍快堆SNCLFR-10 的無保護(hù)超功率瞬態(tài)事故不確定性分析，并研究了輸入?yún)?shù)對(duì)瞬態(tài)安全參數(shù)的影響程度。結(jié)論如下：

1）基于非參數(shù)統(tǒng)計(jì)法的不確定性計(jì)算，得到的目標(biāo)參數(shù)95/95不確定帶能夠較好地包絡(luò)名義值；對(duì)比目標(biāo)參數(shù)不確定性包絡(luò)帶，可得不確定性輸入?yún)?shù)傳播到輸出參數(shù)的不確定性隨時(shí)間而變化；包殼峰值溫度95/95 容忍上限值為1 208.66 K，燃料峰值溫度的95/95容忍上限值為2 757.25 K。

2）通過Spearman 相關(guān)系數(shù)法對(duì)SNCLFR-10 的UTOP 事故進(jìn)行敏感性分析。結(jié)果表明：包殼內(nèi)直徑、穩(wěn)態(tài)燃料平均溫度、初始功率、第1 組衰變常數(shù)以及反應(yīng)性反饋系數(shù)顯著影響燃料和包殼峰值溫度，其中燃料平均溫度和反應(yīng)性反饋系數(shù)對(duì)瞬態(tài)安全參數(shù)的敏感度最強(qiáng)。

作者貢獻(xiàn)聲明余清遠(yuǎn)：實(shí)施研究，文章撰寫；漆靜雯：采集數(shù)據(jù)，分析/解釋數(shù)據(jù)；趙鵬程：統(tǒng)計(jì)分析，獲取研究經(jīng)費(fèi)；趙亞楠：對(duì)文章的知識(shí)性內(nèi)容作批評(píng)性審閱；于濤：行政、技術(shù)支持，指導(dǎo)，支持性貢獻(xiàn)。