潘昕懌　　蘭　兵　　張春明　　靖劍平　　攸國順（環(huán)境保護(hù)部核與輻射安全中心　　北京　100082）

?

多群核數(shù)據(jù)不確定性對堆芯物理計算的影響

潘昕懌蘭兵張春明靖劍平攸國順
（環(huán)境保護(hù)部核與輻射安全中心北京100082）

摘要核數(shù)據(jù)不確定性是造成反應(yīng)堆物理計算結(jié)果不確定性的重要因素之一?；谒璩闃雍藬?shù)據(jù)的協(xié)方差矩陣開發(fā)了隨機(jī)抽樣模塊（StochasticSampling，SAMP），在此基礎(chǔ)上利用SCALE（StandardizedComputerAnalyses forLicensingEvaluation）軟件包實現(xiàn)了混合法和隨機(jī)抽樣法兩種不確定性分析方法，以研究多群核數(shù)據(jù)不確定性對堆芯物理計算的影響。以3×3假想堆芯為對象，對兩種方法進(jìn)行了驗證，然后應(yīng)用于國際原子能機(jī)構(gòu)（InternationalAtomicEnergyAgency，IAEA）燃料管理基準(zhǔn)題中的Almaraz核電廠首循環(huán)堆芯。分析結(jié)果表明，兩種方法結(jié)果符合良好，Almaraz核電廠堆芯keff不確定性約為0.5%，堆芯徑向和軸向功率的最大不確定性分別為1.9%和0.45%。

關(guān)鍵詞多群核數(shù)據(jù)，不確定性，SCALE，隨機(jī)抽樣，協(xié)方差矩陣

國家科技重大專項（No.2013ZX06002001）資助

第一作者：潘昕懌，男，1983年出生，2010年畢業(yè)于清華大學(xué)，主要從事反應(yīng)堆物理和核安全研究

SupportedbyNationalScienceandTechnologyMajorProject（No.2013ZX06002001）

Firstauthor:PANXinyi，male，bornin1983，graduatedfromTsinghuaUniversityin2010，engagedintheresearchofreactorphysicsandsafety Correspondingauthor:YOUGuoshun，E-mail:youguoshun@chinansc.cn

Influence of multigroup nuclear data uncertainties on the reactor core physics calculation

PANXinyiLANBingZHANGChunmingJINGJianpingYOUGuoshun
（Nuclear and Radiation Safety Center， Ministry of Environmental Protection， Beijing 100082， China）

Abstract Background:Theuncertaintyofnucleardataisoneofthekeyfactorsresultingintheuncertaintyof reactorphysicscalculation.Purpose:Theinfluenceofmultigroupnucleardatauncertaintiesonthereactorcore physicscalculationwasstudiedinthispaper.Methods:ThestochasticsamplingmodularSAMPbasedoncovariance matrixofnucleardatawasdeveloped，andthehybridmethodandstochasticsamplingmethodwererealizedusing SCALE（StandardizedComputerAnalysesforLicensingEvaluation）softwarepackage.Thetwomethodswere validatedusing3×3hypotheticalcoreandthenappliedtothefirstcycleofAlmarazpressurized-waterreactor（PWR）intheIAEA（InternationalAtomicEnergyAgency）fuelmanagementbenchmark.Results:Resultsofthetwo methodsareingoodagreement.Theuncertaintyofcoreeffectivemultiplicationfactorisabout0.5%，andthe maximumuncertaintiesoftheradialandaxialpowerareabout1.9%and0.45%respectivelyinAlmarazPWR. Conclusion:Thetwo-stepmethodandstochasticsamplingmethodcanbothbeusedfortheuncertaintyanalysisof reactorcorecalculation.

Key words Multigroupnucleardata，Uncertainty，SCALE，Stochasticsampling，Covariancematrix

利用最佳估算技術(shù)進(jìn)行反應(yīng)堆設(shè)計和安全分析時，為確保足夠的安全裕量，必須進(jìn)行不確定性分析。反應(yīng)堆物理計算的不確定性主要包括工程不確定性（如燃料制造公差）、計算不確定性（如計算模型簡化和核數(shù)據(jù)等輸入?yún)?shù)不確定性）和現(xiàn)象不確定性（如密實化、棒彎曲和氙振蕩等效應(yīng)）[1]。近年來研究表明，核數(shù)據(jù)的不確定性是影響反應(yīng)堆物理計算的重要因素之一[2]。

核數(shù)據(jù)不確定性主要源自其測量過程及評價模型[3]，通常以協(xié)方差的形式與評價核數(shù)據(jù)庫同時發(fā)布。美國橡樹嶺國家實驗室開發(fā)的SCALE （StandardizedComputerAnalysesforLicensing Evaluation）軟件包制作了核數(shù)據(jù)協(xié)方差庫文件44groupcov，即包含401種核素44群核數(shù)據(jù)的相對方差和協(xié)方差[4]。

反應(yīng)堆物理計算分為組件計算和堆芯計算兩步進(jìn)行。目前研究核數(shù)據(jù)不確定性對組件計算影響的方法主要分為基于伴隨計算的微擾法和隨機(jī)抽樣法兩類，前者包括SCALE/TSUNAMI[5]、SUS3D[6]等程序，后者如XSUSA[7]程序。SCALE軟件包TSUNAMI-3D序列可處理包括堆芯在內(nèi)的三維幾何，但僅能對有效增殖因數(shù)keff進(jìn)行不確定性分析，而無法得到堆芯功率等其他參數(shù)的不確定度。本文以上述兩種組件計算不確定性分析方法為基礎(chǔ)，分別實現(xiàn)了混合法和隨機(jī)抽樣法用于堆芯物理計算不確定性分析，并以3×3假想堆芯為對象，驗證兩種方法的有效性。然后針對國際原子能機(jī)構(gòu)（InternationalAtomicEnergyAgency，IAEA）壓水堆燃料管理基準(zhǔn)題[8]中的Almaraz核電廠首循環(huán)堆芯，研究多群核數(shù)據(jù)不確定性對堆芯keff及堆芯功率等參數(shù)的影響。

1　分析方法

1.1抽樣模塊(Stochastic Sampling, SAMP)

在兩種堆芯不確定性分析方法中，需利用多群微觀核數(shù)據(jù)或雙群宏觀截面的協(xié)方差矩陣進(jìn)行抽樣，從而獲得相關(guān)核數(shù)據(jù)的隨機(jī)抽樣值。本文專門開發(fā)了計算模塊SAMP實現(xiàn)上述抽樣。

假設(shè)由所需抽樣的I個核數(shù)據(jù)構(gòu)成的向量為x，其相對協(xié)方差矩陣rΣ已知且為對稱半正定，則先對其進(jìn)行特征值分解：

將矩陣rE的各對角元素（即協(xié)方差矩陣rΣ的特征值，為非負(fù)實數(shù)）進(jìn)行開方后可得到矩陣然后通過式（2）得到矩陣

假設(shè)x中所有核數(shù)據(jù)均為正態(tài)分布，則可通過如下隨機(jī)抽樣得到擾動因子向量P[9]：

式中：i是元素全為1的向量。（0，1）N為I個均值為0、標(biāo)準(zhǔn)差為1的獨立正態(tài)分布隨機(jī)變量組成的隨機(jī)向量，可利用蒙特卡羅法或拉丁超立方體法等方法對（0，1）N進(jìn)行抽樣，進(jìn)而可得到x的抽樣值：

式中：μ為x的參考值或均值向量。根據(jù)μ和rΣ，可以得到x的絕對協(xié)方差矩陣Σ中的各個元素：

式中：1/2（0，1）

ΣN為擾動增量。保證在γ置信度下有α（百分比）的計算結(jié)果處于容差間隔之內(nèi)所需的最小抽樣數(shù)N可由Wilks公式[10]給出：

由式（7）可知，N與需抽樣核數(shù)據(jù)的個數(shù)I及其分布類型無關(guān)，當(dāng)γ和α取95%時N為93。綜上所述，SAMP模塊的輸入為N、μ、rΣ或Σ，輸出為x的N組抽樣值。

1.2混合法

混合法堆芯不確定性分析流程如圖1所示。首先利用基于伴隨計算的微擾法進(jìn)行組件不確定分析，得到每種燃料組件的雙群均勻化宏觀截面的協(xié)方差矩陣，然后利用SAMP模塊隨機(jī)抽樣獲得N組不同的雙群宏觀截面，之后進(jìn)行N次堆芯計算及不確定性分析。

組件不確定性分析使用SCALE軟件包TSUNAMI-2D序列完成。TSUNAMI-2D首先調(diào)用BONAMI和NITAWL模塊進(jìn)行共振計算，然后調(diào)用NEWT模塊進(jìn)行輸運計算和伴隨輸運計算，最后SAMS模塊利用44groupcov庫中44群核數(shù)據(jù)協(xié)方差矩陣iΣ計算雙群均勻化宏觀截面的不確定度，同時生成靈敏度系數(shù)矩陣S。其中靈敏度系數(shù)定義為雙群宏觀截面相對各個核數(shù)據(jù)的偏導(dǎo)數(shù)。SCALE軟件包的TSUNAMI-IP序列可通過式（8）計算得到雙群截面的相對協(xié)方差矩陣rΣ：

式中：矩陣iΣ大小為M×M，M等于目標(biāo)組件涉及到的核素-核反應(yīng)對數(shù)目與能群數(shù)（44群）的乘積；為I×I的對稱矩陣，I為堆芯不確定性分析中所需考慮的雙群均勻化宏觀截面的個數(shù)。本文后續(xù)計算中將考慮每種燃料組件的總截面（Σtot，1、Σtot，2）、吸收截面（Σabs，1、Σabs，2）、裂變相關(guān)截面（vΣfis，1、vΣfis，2、）和散射截面等11個宏觀截面，其中腳標(biāo)1和2分別代表快群和熱群。

圖1　混合法堆芯物理計算不確定性分析流程Fig.1　Flowchart　of　hybrid　uncertainty　analysis　method　for core　physics　calculation.

雙群宏觀截面參考值μ由SCALE軟件包TRITON/NEWT序列通過組件計算獲得，計算時采用SCALE內(nèi)置44群AMPX核數(shù)據(jù)庫。之后SAMP模塊抽樣生成N個宏觀截面文件供堆芯計算程序PARCS使用。

混合法的主要優(yōu)點是能夠一次性生成后續(xù)堆芯不確定性分析所需的雙群宏觀截面的協(xié)方差矩陣，避免重復(fù)進(jìn)行耗時較長的組件計算，同時可利用靈敏度系數(shù)衡量不同核數(shù)據(jù)的貢獻(xiàn)。其缺點是分析參數(shù)和范圍受到組件不確定性分析程序的限制，目前TSUNAMI-2D無法直接給出不連續(xù)因子、形狀因子和中子動力學(xué)參數(shù)等核數(shù)據(jù)的不確定度及靈敏度系數(shù)，同時實現(xiàn)組件燃耗計算和分支計算的不確定性分析較為困難。

1.3隨機(jī)抽樣法

隨機(jī)抽樣法堆芯不確定性分析流程如圖2所示。首先利用SAMP模塊和44groupcov庫對44群AMPX核數(shù)據(jù)庫中相關(guān)核數(shù)據(jù)進(jìn)行隨機(jī)抽樣，得到N組44群核數(shù)據(jù)庫后分別利用TRITON/NEWT序列和PARCS程序依次進(jìn)行N次組件計算和堆芯計算，最后對堆芯計算結(jié)果進(jìn)行不確定性分析。

圖2　隨機(jī)抽樣法堆芯不確定性分析流程Fig.2　Flowchart　of　stochastic　sampling　uncertainty　analysis method　for　core　physics　calculation.

在進(jìn)行隨機(jī)抽樣前，需先從44groupcov庫中讀取相關(guān)核素-核反應(yīng)對的相對方差和協(xié)方差數(shù)據(jù)，包括同一核素不同核反應(yīng)之間以及不同核素之間的協(xié)方差，構(gòu)成大小為M×M的相對協(xié)方差矩陣，其中M等于所考慮的核素-核反應(yīng)對數(shù)目與能群數(shù)（44群）的乘積。44groupcov庫中同一核素不同核反應(yīng)之間的協(xié)方差矩陣可能不對稱，在進(jìn)行式（1）特征值分解前需進(jìn)行對稱化處理，即對稱位置兩個元素取其中絕對值較大的一個。協(xié)方差矩陣還可能出現(xiàn)極少量負(fù)特征值，此時需在進(jìn)行式（2）計算之前將負(fù)特征值改為0。因為對結(jié)果不確定度起主要作用的是同一核素單個核反應(yīng)的方差，且不對稱元素和負(fù)特征值占比極小，上述對稱半正定化處理對不確定性分析結(jié)果的影響不大。

實際應(yīng)用時通常僅對具有協(xié)方差數(shù)據(jù)且靈敏度較高的核數(shù)據(jù)進(jìn)行抽樣，為保證AMPX庫中所有核數(shù)據(jù)之間的一致性，在SAMP模塊隨機(jī)抽樣之后需做如下修正：（1）總反應(yīng)截面、吸收截面等集總核數(shù)據(jù)可表示為其他若干截面之和，若對其中某個分量截面進(jìn)行了抽樣，則與其相關(guān)的所有集總核數(shù)據(jù)均需加上相同的擾動增量；（2）某能群的散射截面可表示為從該能群向所有能群的散射截面分量之和，這些分量保存為二維散射表，若對44群散射截面進(jìn)行了抽樣，則需將表中所有相關(guān)分量乘以相同的擾動因子。同時需對抽樣后的各裂變核素裂變能譜進(jìn)行歸一化。最后將原始AMPX庫中所有相關(guān)核數(shù)據(jù)替換為抽樣值或修正值后生成N組44群核數(shù)據(jù)庫文件。

SAMP模塊能夠處理問題中所有涉及到的核素和核數(shù)據(jù)類型。對于本文中的常規(guī)壓水堆新燃料組件，可僅考慮H、10B、16O、235U、238U和天然Zr等6種核素（或元素）的彈性散射截面、非彈性散射截面、裂變截面、輻射俘獲截面、裂變能譜和每次裂變釋放平均中子數(shù)。對于非零燃耗燃料組件，可考慮其他錒系核素和裂變產(chǎn)物。

隨機(jī)抽樣法的主要優(yōu)勢是可考慮任意核數(shù)據(jù)的不確定性在反應(yīng)堆物理兩步法計算流程中的傳播，并可考慮所有組件計算結(jié)果不確定性對堆芯計算的影響。其缺點是組件不確定性分析耗時長，且作為一種“黑盒”方法，難以確定不同核數(shù)據(jù)的貢獻(xiàn)。

2　方法驗證

以一個由9盒相同燃料組件裝載的3×3假想反應(yīng)堆為對象進(jìn)行方法驗證。燃料富集度為3.1%，不含可燃毒物和控制棒，慢化劑和反射層均為水，具體組件設(shè)計參數(shù)見文獻(xiàn)[8]。分別使用混合法和隨機(jī)抽樣法對該假想堆芯進(jìn)行不確定分析，兩種方法中抽樣次數(shù)N均取150。組件計算時硼濃度取零，燃料溫度和慢化劑溫度均為291.4oC。堆芯計算工況為零燃耗、無可溶硼和熱態(tài)零功率（HotZeroPower，HZP）。

表1給出了混合法和隨機(jī)抽樣法得到的堆芯keff統(tǒng)計分析結(jié)果，同時與TSUNAMI-3D序列結(jié)果進(jìn)行比較。TSUNAMI-3D使用蒙特卡羅計算模塊KENO進(jìn)行堆芯三維建模和輸運計算，并使用微擾法進(jìn)行不確定性分析，可作為keff計算值和不確定度的參考標(biāo)準(zhǔn)。結(jié)果對比表明，本文兩種方法得到的keff平均值和不確定度與TSUNAMI-3D均基本一致，其中keff相對標(biāo)準(zhǔn)差約為0.5%。

表2比較了兩種方法得到的3.1%燃料組件11個宏觀截面的相對標(biāo)準(zhǔn)差。其中混合法的標(biāo)準(zhǔn)差結(jié)果由TSUNAMI-2D直接給出，隨機(jī)抽樣法則由150次組件計算結(jié)果統(tǒng)計得到。結(jié)果表明，每一種宏觀截面的相對標(biāo)準(zhǔn)差均符合良好。TSUNAMI-2D作為組件不確定性分析的標(biāo)準(zhǔn)，進(jìn)一步驗證了隨機(jī)抽樣方法中核數(shù)據(jù)抽樣的準(zhǔn)確性。

表1　3×3假想堆芯keff統(tǒng)計分析結(jié)果Table 1 Statistical analysis results for keffof 3×3 hypothetical core.

表2　3.1%燃料組件雙群宏觀截面相對標(biāo)準(zhǔn)差Table 2 Relative standard deviation results for tow-group macroscopic cross-section of 3.1% fuel assembly (%).

圖3比較了堆芯徑向功率分布的統(tǒng)計分析結(jié)果。其中相對功率參考值由TRITON/NEWT和PARCS依次進(jìn)行組件計算和堆芯計算得到，組件計算采用SCALE內(nèi)置44群AMPX核數(shù)據(jù)庫。結(jié)果表明，混合法和隨機(jī)抽樣法計算得到的各組件相對功率平均值與參考值符合良好，二者得到的功率相對標(biāo)準(zhǔn)差也基本相符。上述結(jié)果驗證了混合法和隨機(jī)抽樣法的準(zhǔn)確性，二者均可用于堆芯keff和功率分布等參數(shù)的不確定性分析。

圖3　3×3假想堆芯徑向功率分布統(tǒng)計分析結(jié)果Fig.3　Statistical　analysis　results　for　radial　power　distribution of　3×3　hypothetical　reactor　core.

3　方法應(yīng)用

圖4顯示了IAEA壓水堆燃料管理基準(zhǔn)題[8]中Almaraz核電廠的首循環(huán)堆芯裝載情況。

該反應(yīng)堆包含157盒燃料組件，根據(jù)燃料富集度和可燃毒物棒（硼玻璃）數(shù)量的不同分為7種新組件。分別使用混合法和隨機(jī)抽樣法對該堆芯進(jìn)行不確定分析，抽樣次數(shù)N取150?；旌戏ㄖ懈鶕?jù)式（8）得到的協(xié)方差矩陣大小為77×77。組件計算時硼濃度取零，燃料溫度和慢化劑溫度均為291.4oC，堆芯計算工況為壽期初（BeginningofLife，BOL）、控制棒全提出堆外（AllRodsOut，ARO）、無可溶硼和HZP。

圖4　Almaraz核電廠首循環(huán)堆芯裝載圖（1/4堆芯）Fig.4　Loading　arrangement　for　the　first　cycle　of　Almaraz PWR　plant?。?/4　core）.

下述分析結(jié)果中堆芯keff和功率參考值由TRITON/NEWT和PARCS依次進(jìn)行組件和堆芯計算得到。表3給出了堆芯keff統(tǒng)計分析結(jié)果，混合法和隨機(jī)抽樣法給出的keff平均值與參考值非常接近，且二者得到的keff標(biāo)準(zhǔn)差基本相符，均約為0.5%。

表3　Almaraz核電廠堆芯（BOL，HZP，ARO，無硼）keff統(tǒng)計分析結(jié)果Table 3 Statistical analysis results for keffof Almaraz PWR plant (BOL, HZP, ARO, no boron).

圖5顯示了堆芯徑向功率分布的計算結(jié)果，混合法和隨機(jī)抽樣法得到的相對功率平均值與參考值偏差很小。圖6顯示了堆芯徑向功率分布的不確定性分析結(jié)果。兩種方法得到的組件功率相對標(biāo)準(zhǔn)差基本一致。由于堆芯布置方案的徑向不均勻性（包括燃料富集度、可燃毒物棒根數(shù)和反射層等的布置），導(dǎo)致堆芯中心和外圍的組件功率不確定度相對較大，最大值為1.9%左右。

圖7顯示了堆芯軸向功率分布的統(tǒng)計分析結(jié)果。混合法和隨機(jī)抽樣法得到的軸向相對功率平均值與參考值符合良好，二者的相對標(biāo)準(zhǔn)差計算結(jié)果也基本一致。相對于徑向功率，多群核數(shù)據(jù)不確定性對堆芯軸向功率的影響相對較小，由于軸向反射層造成的軸向不均勻性，最大相對標(biāo)準(zhǔn)差出現(xiàn)在堆芯活性段的頂部和底部，為0.45%左右。

圖5　Almaraz核電廠堆芯徑向功率分布平均值（BOL，HZP，ARO，無硼）Fig.5　Mean　values　for　radial　power　distribution　of　Almaraz　PWR　plant?。˙OL，　HZP，　ARO，　no　boron）.

圖6　Almaraz核電廠徑向功率分布不確定性分析結(jié)果（BOL，HZP，ARO，無硼）Fig.6　Uncertainty　analysis　results　for　radial　power　distribution　of　Almaraz　PWR　plant?。˙OL，　HZP，　ARO，　no　boron）.

圖7　Almaraz核電廠軸向功率分布統(tǒng)計分析結(jié)果（BOL，HZP，ARO，無硼）Fig.7　Statistical　analysis　results　for　axial　power　distribution　of　Almaraz　PWR　plant?。˙OL，　HZP，　ARO，　no　boron）.

4　結(jié)語

本文通過開發(fā)基于核數(shù)據(jù)協(xié)方差矩陣的隨機(jī)抽樣模塊SAMP，實現(xiàn)了混合法和隨機(jī)抽樣法兩種方法以研究多群核數(shù)據(jù)不確定性對堆芯物理計算的影響，并以3×3假想堆芯為對象對兩種方法的有效性進(jìn)行了驗證。Almaraz壓水堆核電廠首循環(huán)堆芯不確定性分析結(jié)果表明，兩種方法結(jié)果符合良好，其中堆芯keff不確定度約為0.5%，堆芯徑向和軸向功率分布的最大不確定度分別為1.9%和0.45%左右。目前兩種方法仍有其不足之處，混合法的應(yīng)用受組件不確定性分析程序分析能力的限制，而隨機(jī)抽樣法計算時間較長，且未考慮共振計算中引入的隱含不確定性，在后續(xù)工作中需進(jìn)一步研究改進(jìn)。

參考文獻(xiàn)

1 司勝義.堆芯核設(shè)計程序系統(tǒng)的不確定性綜述[C].第十一屆反應(yīng)堆數(shù)值計算和粒子輸運學(xué)術(shù)會議暨2006年反應(yīng)堆物理會議，哈爾濱，中國，2006

SIShengyi.Reviewontheuncertaintyofsoftwaresystem forreactorcorenucleardesign[C].11stAcademic ConferenceontheReactorNumericalCalculationand ParticleTransportandReactorPhysicsConference2006，Harbin，China，2006

2 IvanovK，AvramovaM，KamerowS.Overviewand discussionofphaseIoftheOECDLWRUAM benchmarkactivity[C].InternationalConferenceon Mathematics，ComputationalMethods＆ReactorPhysics （M＆C2009），NewYork，USA，2009

3 ZhouXM，LiuGM，LiD，et al.Usingactivationmethod tomeasureneutronspectruminanirradiationchamberof aresearchreactor[J].NuclearScienceandTechniques，2014，25（1）:010603.DOI:10.13538/j.1001-8042/nst.25. 010603

4 WilliamsML，WiardaD，ArbanasG，et al.SCALE nucleardatacovariancelibrary[R].ORNL/TM-2005/39，OakRidgeNationalLaboratory，2011

5 ReardenBT，JessseeMA，WilliamsML.TSUNAMI-1D: controlmoduleforonedimensionalcross-section sensitivityanduncertaintyanalysisforcriticality[R]. ORNL/TM-2005/39，OakRidgeNationalLaboratory，2005

6 KodeliI.Multidimensionaldeterministicnucleardata sensitivityanduncertaintycodesystem:methodand application[J].NuclearScienceandEngineering，2001，138（1）:45-66.DOI:10.13182/NSE00-43

7 ZwermannW，Krzykacz-HausmannB，GallnerL，et al. Uncertaintyanalysiswithnuclearcovariancedatain reactorcorecalculations[J].JournaloftheKorean PhysicalSociety，2011，59（2）:1256-1259.DOI: 10.3938/jkps.59.1256

8 AhnertC，AragonesJM，MerinoF，et al.In-corefuel managementcodepackagevalidationforPWRs[R]. IAEA-TECDOC-815，InternationalAtomicEnergy Agency，1995

9 BevingtonPR，RobinsonK.Datareductionanderror analysisforthephysicalsciences[M].Boston: McGraw-Hill，2003

10 WilksS.Determinationofsamplesizesforsetting tolerancelimits[J].TheAnnalsofMathematicalStatistics，1941，12（1）:91-96.DOI:10.1214/aoms/1177731788

收稿日期：2015-08-31，修回日期：2015-10-23

通信作者：攸國順，E-mail:youguoshun@chinansc.cn

DOI:10.11889/j.0253-3219.2016.hjs.39.010602

中圖分類號TL329