CEFR一臺一回路泵停運疊加失廠外電非對稱工況的三維瞬態熱工水力特性數值模擬

馬翔鳳，陸道綱，梁江濤，符精品，劉一哲，楊 軍，郭忠孝，張鈺浩,*

CEFR一臺一回路泵停運疊加失廠外電非對稱工況的三維瞬態熱工水力特性數值模擬

馬翔鳳1,2，陸道綱1,2，梁江濤1,2，符精品1,2，劉一哲3，楊 軍3，郭忠孝3，張鈺浩1,2,*

（1. 華北電力大學 核科學與工程學院，北京102206；2.非能動核能安全技術北京市重點實驗室，北京102206；3. 中國原子能科學研究院，北京 102413）

一臺一回路泵停運疊加失廠外電事故下堆容器及堆內非對稱三維熱工水力特性對于池式快堆的設計與安全運行十分重要。池式鈉冷快堆結構復雜，流動路徑多，現有系統程序難以準確獲得快堆非對稱事故下的三維熱工水力特征。本研究基于三維數值計算方法，建立CEFR冷熱鈉池全尺寸一體化模型，模擬了CEFR一臺一回路泵停運疊加失廠外電這一典型非對稱瞬態事故工況下的三維熱工水力行為，特別是主泵惰轉及返流的瞬態過程，揭示了鈉池內三維非對稱流場、溫度分布及IHX進出口溫度等關鍵熱工參數的瞬態變化特性。計算結果表明，在事故前1 500 s，冷卻劑自停運環路IHX出口向上返流至熱鈉池并通過正常環路IHX出口進入冷鈉池，正常環路IHX出口平均溫度在600 s左右出現極大值約491.9 ℃，而停運環路IHX出口溫度持續上升并逐漸與正常環路趨于一致。該計算結果可為該工況下反應堆安全評價及結構應力分析提供關鍵數值參考。

CEFR一體化模型；主泵停運；失廠外電；非對稱；三維數值模擬

中國實驗快堆（CEFR）在發生一臺一回路泵停運疊加失廠外電等設計基準事故時，堆內會呈現特殊的非對稱瞬態變化特點，該事故下的堆內瞬態溫度分布變化特性對于鈉冷快堆的安全運行具有十分重要的意義。由于池式鈉冷快堆結構布置復雜，堆內冷卻劑流動呈現多路徑的特點，且不同路徑之間空間尺寸差異較大，瞬態事故下影響流動和傳熱的因素較多，因此明確該事故工況下快堆鈉池內的三維熱工水力特性及其瞬態變化一直存在較大困難。通過對鈉冷快堆熱工水力特性的研究，確保在穩態和瞬態的條件下，冷卻劑能夠及時帶走堆芯熱量，維持結構完整性，這對快堆的安全評價具有非常重要的意義。因此對池式鈉冷快堆在關鍵事故工況下的熱工水力特性進行分析與安全評價是核電技術研發和安全評審中的關鍵因素，但鈉冷快堆開展實驗成本高昂、難度很大，因此對于各種事故工況下池式鈉冷快堆三維瞬態熱工水力研究一直存在嚴峻挑戰。

采用兩環路和多模塊蒸汽發生器設計方案的核電廠，在運行過程中，可能會由于某一環路的泵或蒸汽發生器模塊導致熱量排出功能受限，從而出現不對稱工況[1-2]。在CEFR運行期間，可能發生失廠外電事故，之后反應堆觸發緊急停堆信號。堆芯功率迅速降低，主循環泵和給水泵失去供電，進行自然惰轉。故障環路中的主泵惰轉，轉速逐漸降低至0，正常環路中的主泵依靠應急柴油機的支持，轉速降低并維持到一定水平。后續熱鈉池中的熱流體直接通過IHX進入冷鈉池，正常環路流量大部分進入堆芯，另一部分通過故障主泵返流，經故障環路IHX進入熱鈉池。因此在事故期間鈉池會經歷復雜的功率及流量變化，使得池式鈉冷快堆內的流動及溫度變化難以準確獲得，而快堆結構復雜性及事故進程的不確定性同樣使得實驗研究非常困難，以上對于池式快堆的安全分析與評價帶來巨大挑戰。

國內外已有多種數值模擬軟件和系統程序用于核工程領域，如美國SASSYS-1，法國OASIS等[3-5]。大部分系統軟件程序是建立一維或簡化的多維模型以分析核電站瞬態過程中的系統響應，對于反應堆整體三維參數進行的研究相對較少。

任麗霞等[6]采用FASYS程序分析CEFR40%功率下失去廠外電后單臺主泵停運的一回路瞬態特性分析；張熙司等[7]運用OASIS程序通過模擬鈉泵切除試驗；楊曉燕等[8]分析了非對稱工況下對鈉冷快堆堆芯入口溫度的影響。國內學者對于CEFR的一回路系統也進行了一定研究。朱桓軍等[9]運用CFD軟件對CEFR進行簡化建模，獲得了其在額定運行工況下的速度、溫度分布，并分析了堆內冷、熱鈉池的初步三維熱工分布；崔滿滿[10]利用Fortran語言開發了CEFR一回路熱工水力穩態計算程序，對CEFR的一回路系統在滿功率下的穩態熱工水力特性進行計算分析。

在正常運行工況下，鈉處于高溫狀態，冷熱鈉池中存在復雜的熱分層現象，堆內眾多的設備、池式結構的復雜性以及硬件條件的限制，使得同時建立冷熱鈉池三維系統分析方案存在較大困難。本文基于CEFR結構和運行特點，采用分區建模、一體化耦合計算的方法，建立詳細的冷、熱鈉池三維模型，運用商用CFD軟件FLUENT，實現一臺一回路泵停運疊加失廠外電典型非對稱事故工況下三維瞬態計算。通過對1 500 s長瞬態計算分析，獲得事故工況下一回路主熱傳輸系統的熱工水力現象及三維非對稱溫度、流場變化規律。此三維計算結果可為反應堆結構部件的力學分析提供輸入條件，以及為核電廠在該事故工況下的安全運行分析提供數值參考。

1　CEFR計算模型及簡化方案

本研究以CEFR為研究對象，采用1:1全尺寸計算模型，考慮主要部件的結構功能完整性和整體網格數量保守性對模型進行合理簡化[11]，具體結構介紹如下：

CEFR一回路系統由兩條并聯的環路組成，每條環路設置一臺鈉循環泵、兩臺IHX、一臺DHX及其相應的管道，共同組成反應堆主容器內一體化池式結構。一回路主冷卻劑泵從冷池吸入冷鈉，由葉輪加壓，通過壓力管傳遞到柵板聯箱，后流經堆芯，吸收熱量，從堆芯流出后進入熱鈉池。熱鈉經IHX冷卻，將熱量傳遞給二回路鈉，進入冷池后，然后重新被循環泵吸入完成一回路循環流動。

經過上述流動分析，將整體一回路系統分為冷鈉池、熱鈉池、堆芯、主容器冷卻系統等各部分，采取模塊化建模及網格劃分方式，最后進行耦合計算。一回路系統主要包括堆芯、主泵、IHX、DHX、堆內屏蔽、支承板等關鍵部件，采用多孔介質設置方法對于堆芯等部分進行相應的簡化，增加部件壁厚，利用FLUENT用戶自定義函數（UDF）實現堆芯、IHX、DHX功率及堆內流量的瞬態變化，更加精準地模擬CEFR堆內的三維熱工水力現象。

1.1　冷鈉池區域

在冷鈉池中，兩組環路之間由隔板分隔，此外還有中間隔板將冷池分隔成上下兩部分區域，相同環路中的上下兩部分可以連通。大部分設備集中分布在冷池的上半部區域，每臺泵分別對應貫穿冷池中部隔板的三通型壓力管，模擬實際泵加壓注入功能，通過入口邊界條件將流體注入壓力管，然后進入柵板聯箱。冷池中塊布置多塊肋板，部分肋板中間開孔，在起支撐作用的同時，將冷池空間連通或隔離等。簡化后的冷鈉池模型示意如圖1所示。

圖1　冷鈉池結構示意圖

1.2　熱鈉池區域

1.2.1堆芯及柵板聯箱

堆芯為系統提供熱源邊界條件；柵板聯箱主要對進入堆芯的流量進行分配。本方案簡化堆芯組件內部燃料棒的復雜結構，按照堆芯整體形狀將堆芯劃分為活性區和非活性區，設置不同功率密度的體積熱源，模擬堆芯功率分布，進而為冷、熱鈉池計算提供精確的邊界條件。柵板聯箱內部簡化為多孔介質模型，通過調節阻力系數以滿足流量分配的要求。

1.2.2一回路主循環泵及泵支承冷卻系統

泵結構是整體計算域的入口及出口邊界邊界。簡化泵外部的熱屏蔽、波紋補償器以及泵內部的電機、葉輪等結構。將泵吸入口設置為壓力出口邊界條件；將泵的加壓排出口設置為入口邊界條件，實現驅動堆內冷卻劑的循環流動的功能，按照設計參數及計算過程中監測數據設置泵及泵支承冷卻系統入口、出口參數。

1.2.3堆芯圍筒

簡化堆芯圍筒內部棒束結構，內部設置多孔介質。基于 1/6實際結構建模，開展數值實驗，得到阻力系數等參數。水平熱屏蔽分隔冷池和熱池，由三塊水平板組成，該區域流動較弱，簡化為固體區域，各換熱相關物性參數按照鈉和鋼的加權平均值設置。

1.2.4IHX、DHX

IHX和DHX簡化內部傳熱管結構，換熱區域設置為體積冷源，通過UDF實現功率變化，模擬二次側傳熱冷卻過程，主要為計算提供精確邊界條件。正常運行過程中主要依靠IHX將堆芯熱量傳遞給二次側回路；而事故工況下，特別是IHX功率完全喪失以后，主要依靠DHX建立自然循環，從而帶走堆芯衰變熱。在建模過程中，將IHX和DHX簡化為多層圓筒形結構，增加筒體和支承壁面厚度，內部簡化為多孔介質模擬一次側流動過程，基于工程設計壓降進行數值實驗以確定阻力系數。

1.2.5主容器冷卻系統及徑向熱屏蔽

主容器冷卻系統主要為減小高溫鈉池對主容器的熱沖擊，通過隔板將流道分隔成上升和下降通道，與熱池內三層徑向熱屏蔽實體采用結構化網格進行一體化建模。下封頭區域為橢球型，采用非結構化網格進行網格劃分。

簡化后的熱鈉池模型如圖2所示。

圖2　熱鈉池結構示意圖

1.3　非對稱工況建模處理方法

一臺一回路泵停運-疊加失廠外電事故是典型的非對稱工況，在建模過程中需考慮實現非對稱計算的建模處理方法。在正常運行或主泵惰轉期間，需設置兩個環路的泵入口及出口以及通過調節IHX的阻力參數，使得堆芯流量按一定的比例進入IHX，從而實現流量分配，使得相應環路中流入與流出流量平衡。當出現返流以后，可以將故障環路泵的出口邊界與入口邊界互換，從而模擬故障環路返流過程。具體處理過程如圖3所示。

（a）惰轉過程流量分配處理方法

（b）返流過程流量分配處理方法

2　計算設置及關鍵邊界條件

在本計算中，將上述不同區域網格進行一體化組裝，各區域交界處位置采用interface連接，以實現不同區域網格的耦合計算；采用精度合理的標準模型，綜合考慮計算精度及計算經濟性需要合理選擇時間步長；對于堆芯、IHX、堆內屏蔽等密集管束結構，采用多孔介質模型[11]，調節阻力系數實現堆芯及不同環路的流量比例分配，使之與實際運行狀態保持一致。計算過程中關鍵邊界條件主要包括：

（1）源項：100%穩態運行下，堆芯總功率為65 MW，IHX設置正常運行的冷卻功率16.25 MW/臺，DHX設置備用冷卻功率為0.052 5 MW/臺；瞬態計算中，以100%穩態功率運行狀態作為瞬態初始條件，堆芯熱功率及停運環路IHX冷卻功率迅速降低，經過一段時間后堆芯保持一定的衰變熱功率，停運環路IHX冷卻功率迅速降為0，正常環路IHX冷卻功率在1 000 s左右降為0，DHX冷卻功率由計算過程中監測的DHX入口溫度確定。熱功率及冷卻功率瞬態輸入如圖4所示。

（2）系統流量：100%穩態運行下，堆芯流量為301 kg/s，主容器冷卻系統流量為40 kg/s，同樣將其作為瞬態初始流量條件，兩個環路對稱設置，流量分配比例為1:1；停運工況下，鈉泵首先進行惰轉，停運環路流量在50 s左右降為0，正常環路由于存在廠內電源的支持，所以其流量在下降到一定程度（約初始流量的36%）以后會保持穩定。之后來自正常環路泵的流量一部分進入堆芯進行冷卻，一部分通過故障泵返流至停運環路（返流流量用負值表示），流量瞬態變化輸入曲線如圖5所示。入口（mass-flow-inlet）流量通過UDF進行設置，模型入口溫度通過計算過程中監測的相應位置出口（pressure-outlet）溫度確定，使得系統入口溫度與出口溫度相等，保證計算連續性。

圖4　堆內相對功率變化

圖5　不同環路相對流量變化

（3）冷卻劑物性參數：考慮冷卻劑的物理參數隨溫度的變化，主要物性參數包括鈉的密度、熱導、黏度等，對各關鍵物性參數進行隨溫度變化的分段擬合設置。

（4）阻力設置：在堆芯、柵板聯箱、堆芯圍筒、IHX等區域設置多孔介質，根據設計壓降及惰轉流量分配要求進行數值試驗確定阻力系數，瞬態下基于流量分配需要對兩個環路IHX阻力系數進行調節。

（5）計算所設涉及的各處壁面條件包括：內部壁面設置為coupled wall（導熱耦合壁面）；外部壁面采取絕熱設置；流域內的聯通處設置interior（聯通面）；連接不同區域網格交界面位置設置interface（網格連接）。

基于以上工況分析及關鍵邊界條件設置，開展CEFR一回路一臺一回路泵停運疊加失廠外電事故工況下的三維瞬態熱工計算。

3　CFD計算分析

3.1　對稱工況穩態計算

綜合考慮多種條件因素影響，為驗證系統模型的正確性和輸入數據的準確性，為后續瞬態計算過程提供基礎。根據設計參數開展獨立的額定功率穩態工況（環路對稱運行）計算。在100%額定功率下，堆本體中心截面（=0）處溫度分布、IHX中心截面溫度分布、流場分布如圖6、圖7所示。

圖6　池內溫度分布示意圖

圖7　池內流動分布示意圖

從圖6、圖7可以明顯看出，在穩態工況下，堆內溫度分布呈現明顯的熱分層現象，且兩環路保持對稱狀態。冷卻劑進入泵后，經過壓力管傳送至堆芯，經過堆芯加熱，后平均分配至兩個環路，IHX將其冷卻，其中一回路整體溫度范圍為350～540 ℃，堆內最高溫在堆芯出口處附近，最低溫位于IHX出口附近[12]。

將本研究的模型計算的熱工水力參數與中國實驗快堆穩態的設計參數進行對比，經比較得出結果吻合度較高。表1給出CEFR滿功率穩態運行計算值與設計值比較結果。

表1穩態計算結果與設計值比較

Table.1　The compare of calculated results and design values in the steady state

經過比較結果顯示，該計算方案模型合理可靠，能夠用來準確預測一回路冷卻劑流動分布。將該穩態結果作為一回路熱工水力瞬態計算輸入的初始條件。

3.2　非對稱事故瞬態計算

由穩態運行進入預計瞬態工況后，考慮主泵首先進行約50 s的惰轉，兩個環路的流量逐漸下降，停運環路流量在50 s左右降為0，當故障泵完全停止運行以后，一部分正常環路的流量通過故障泵返流至停運環路。在計算過程中需考慮兩個環路之間的流量分配及不對稱排余熱功率變化。

（1）瞬態流場分析：事故前50 s，流體通過兩個環路泵入口經壓力管進入柵板聯箱，向上經過柵板聯箱、堆芯，由于堆芯及柵板聯箱內部攪混，使得兩個環路的流量充分混合，50 s后流量逐漸穩定，停運環路出現穩定的返流現象。現分別選取50 s之前與50 s之后時刻進行分析，其中10 s及100 s時刻堆內速度場分布分別如圖8所示。

圖8　瞬態進程下池內流動分布示意圖

從圖8可知，正常環路與停運環路流速逐漸下降，而在50 s之后，停運環路由于返流的影響流速開始逐漸增大，但流動方向與正常環路方向相反。10 s時刻正常環路IHX入口流量約為153.2 kg/s，停運環路IHX入口流量約為42.6 kg/s；在流動穩定以后，正常環路IHX入口流量約為58.1 kg/s，停運環路IHX入口流量約為-20.4 kg/s（負值表示返流）。

（2）瞬態溫度場分析：在瞬態條件下，IHX冷卻功率迅速下降，特別是停運環路IHX功率迅速降為0，正常環路的IHX冷卻功率也將在1 000 s左右降為0。由圖9、圖10事故進程下的IHX相關溫度變化可知，在事故發生前期由于堆芯功率迅速下降使得堆芯出口溫度也有明顯下降，導致IHX入口溫度明顯下降，其中一環路（正常環路）的IHX入口溫度緩慢下降，而二環路（停運環路）IHX由于返流影響入口溫度先下降后上升，并在約600 s時刻達到最低約367.1 ℃；由于熱流體直接通過IHX進入冷鈉池，使得事故前期冷鈉池溫度升高，其中，一環路IHX出口溫度在600 s左右時出現極大值，約為491.9 ℃。同時由于返流影響，一環路流體先被泵吸入后一部分進入堆芯，一部分進入二環路，隨著正常環路一側的冷鈉池溫度升高，停運環路一側的冷鈉池溫度逐漸升高，1 500 s時停運環路IHX出口溫度已接近383.9 ℃。該溫度變化趨勢與系統程序計算結果相似[6]。事故后期正常環路側的冷鈉池在DHX影響下溫度開始下降，自600 s至1 500 s，正常環路IHX出口溫度下降約57.1 ℃，總體下降趨勢較緩。

由上述溫度變化分布分析可知，由于熱鈉池熱量向冷鈉池傳遞，所以冷鈉池內溫度存在上升趨勢，但幅度較緩。選取冷鈉池內中板高度處截面溫度分析，如圖 11 所示（0°-90°-180°為正常一環路側，180°-270°-360°為停運二環路側），0～1 500 s范圍內，中板上的溫度分布整體有上升趨勢，并隨著事故進程，逐漸出現冷池內的三維不對稱現象。為更好分析一回路系統的溫度分布狀況，取中板平面上沿徑向的2條圓形監測線，分析線上的溫度分布。在事故初期，冷池中板上的溫度依然分布較為均勻，800 s時已出現一定的非對稱影響，正常環路溫度逐漸高于停運環路側溫度，其中一環路側與二環路側溫差20～30 ℃。

圖9　瞬態進程下IHX截面溫度分布示意圖

圖10　瞬態進程下IHX進出口溫度變化示意圖

圖11　瞬態進程下冷池中板溫度變化示意圖

4　結論

本文以池式鈉冷快堆為研究對象，采用FLUENT，開展失廠外電下一回路一臺泵停運工況的一體化三維數值計算，以較為準確的穩態結果提供瞬態進程的初始條件，分析了事故下前1 500 s期間熱鈉池、冷鈉池內流場和溫度場非對稱瞬態變化規律，是對現有系統程序計算結果的有力補充。主要計算結論如下：

（1）通過對停運泵惰轉及返流現象的計算，較好地模擬停運事故發生后不對稱下雙環路熱工水力現象的瞬態進程，從三維角度較為準確描述堆內溫度場變化趨勢，獲得了瞬態運行下IHX進出口溫度、堆芯出口溫度等關鍵參數變化特征。

（2）在事故進程中，熱鈉池溫度有一定下降趨勢，1 500 s時堆芯出口溫度已下降至425.4 ℃；而由于熱鈉直接通過IHX進入冷鈉池，且部分冷卻劑返流至故障環路，使得冷鈉池溫度則存在一定上升趨勢，600 s左右兩個環路的IHX出口溫度出現最大溫差約125.7 ℃。

（3）通過對堆內流場分布及IHX截面、冷鈉池中板等關鍵位置溫度分析，堆內各處溫度變化較為緩慢，堆內極值溫度未超過安全限值，可以滿足該事故工況下的反應堆安全要求，此計算結果核電廠發生該嚴重事故下的安全運行提供關鍵數據參考。

［1］ 林超，馮預恒，周志偉．CEFR非對稱運行工況的研究［J］．原子能科學技術，2016，50（06）：1021-1026.

［2］ 王利霞，楊紅義，楊曉燕，等．鈉冷快堆喪失蒸汽發生器模塊后的最佳工況研究［J］．原子能科學技術，2017，51（02）：292-298.

［3］ CAHALAN J E.Modeling developments for the SAS4Aand SASSYS computer codes［R］．USA：Argonne National Laboratory，2012.

［4］ CAHALAN J E.The SAS4A/SASSYS-1safety analysis code system［R］．USA：Argonne National Laboratory，2012.

［5］ 楊紅義，徐銤．OASIS程序的開發與應用［J］．核科學與工程，2001（04）：322-325+340.

［6］ 任麗霞，王晉，胡文軍．中國實驗快堆失去廠外電后單臺主泵停運的一回路瞬態特性分析［J］．核科學與工程，2016，36（01）：35-41.

［7］ 張熙司，胡文軍，李政昕，等．中國實驗快堆1臺一回路泵切除試驗計算模擬與分析［J］．原子能科學技術，2015，49（S1）：283-287.

［8］ 楊曉燕，齊少璞，楊紅義，等．不對稱工況對池式鈉冷快堆堆芯入口溫度的影響［J］．原子能科學技術，2018，52（11）：1977-1983.

［9］ 朱桓君，許義軍，錢曉明．基于附加源項法的鈉冷快堆冷熱池三維分析［J］．核動力工程，2015，36（01）：141-143.

［10］崔滿滿．鈉冷快堆一回路主冷卻系統仿真研究［D］．哈爾濱工程大學，2012.

［11］ Zihan Xia，Daogang Lu，Jiaxuan Tang，et al. Three- dimensional thermal hydraulic transient calculation of coupled cold and hot sodium pools under a loss of feedwater accident in the China experimental fast reactor［J］．Annals of Nuclear Energy，2020，139：107217.

［12］張鈺浩，夏子涵，梁江濤，等．基于一體化三維數值模擬的中國實驗快堆冷鈉池及其堆內構件熱工特性分析［J］．核科學與工程，2020，40（03）：499-507.

Numerical Simulation of Three Dimensional Transient Thermal and Hydraulic Characteristics of CEFR under Asymmetric Conditions of One Primary Pump Trip Superimposing Loss of Off-site Power

MA Xiangfeng1,2，LU Daogang1,2，LIANG Jiangtao1,2，FU Jingpin1,2，LIU Yizhe3，YANG Jun3，GUO Zhongxiao3，ZHANG Yuhao1,2,*

（1. School of Nuclear Science and Engineering，North China Electric Power University，Beijing 102206，China；2. Beijing Key Laboratory of Passive Safety Technology for Nuclear Energy，Beijing 102206，China；3. China Institute of Atomic Energy，Beijing 102413，China）

It is very important for the design and safe operation of pool-type fast reactor to master the asymmetric three dimensional thermal hydraulic characteristics under the accident of one primary pump trip superimposing loss of off-site power. Due to the complex structure and multiple flow paths, it is difficult for the existing system codes to accurately obtain the three-dimensional thermal and hydraulic characteristics of the pool-type sodium-cooled fast reactor under the typical asymmetric accident. This study was based on the three-dimensional numerical calculation method, and a full-scale integrative CEFR cold hot sodium pool model was established. The three-dimensional transient thermal hydraulic behavior under asymmetric conditions of one primary pump trip superimposing loss of off-site power was simulated. In particular, the transient process of running down and reflux was simulated. Transient variation characteristics of three-dimensional asymmetric flow field, temperature distribution and the key thermal parameters, such as the IHX inlet and outlet temperature, were revealed.The calculation results showed that in the 1500 s, the coolant flowed upward back to the hot sodium pool through the IHX outlet of the trip loop and entered the cold sodium pool through the normal loop IHX outlet, and the average outlet temperature of the normal loop IHX appeared a maximum value of about 491.9℃ around 600 s, and then gradually decreased. However, the average outlet temperature of the trip loop IHX continued to rise and gradually tended to be consistent with the normal loop. The calculation results provided the key numerical reference for reactor safety evaluation and structural stress analysis under this working condition.

The integrative model of CEFR; One primary pump trip; Loss of off-site power; Asymmetric; Three dimensional numerical simulation

TL333

A

0258-0918（2021）06-1234-10

2021-07-11

馬翔鳳（1996—），女，山東威海人，碩士研究生，主要從事反應堆熱工水力學研究

張鈺浩，E-mail：zhangyuhao@ncepu.edu.cn