夏 凡,劉書勇,李桃生,梅華平,汪 振,趙吉運
(1. 中國科學院合肥物質科學研究院核能安全技術研究所,合肥 230031;2. 中國科學技術大學,合肥 230026;3. 香港城市大學機械工程系,香港)
第四代鉛冷快堆具備安全性好及易小型化的優勢,在分布式供電、海島平臺等領域具有廣闊的應用前景。在安全性方面,無論是回路式還是池式設計,堵流都是重要的事故工況,必須給予考慮[1-3]。堵流的起因主要是由于氧控不良引起氧化鉛在組件內的沉積[4],或是由于腐蝕脫落的結構材料碎片隨冷卻劑帶入,造成組件內的堵流。燃料組件內堵流發生后會進而引發過熱和局部流量下降,影響堆的安全運行[5-8]。
對于液態鉛鉍冷卻的堵流實驗研究,Pacio J 等開展了塊狀堵塊實驗,通過設置低導熱系數的固體堵塊,研究了不同堵流工況下的流動傳熱[9]。該堵流實驗指出,實際堆運行中堵塊多為多孔介質,由于多孔堵塊的加工制造存在困難,因此該實驗采用的堵塊并沒有設置孔隙率。2019 年,意大利ENEA 基于鉛鉍冷卻回路(NACIE-UP)開展了堵流實驗研究,結果表明,堵流后方形成局部溫度峰值,堵流子通道活性區末端形成全局溫度峰值[10]。對于堵流模擬研究,2014 年,Di Piazza 等建立CFD 模型對無繞絲工況下燃料組件堵流現象進行了數值模擬,研究結果表明:對于多子通道堵流,在堵塊下游的回流區域出現溫度峰值;對于少子通道堵流,溫度峰值出現在活性區域末端[11]。2019 年,上海交通大學柴翔等基于KIT 實驗裝置,采用多子通道堵流方式進行了模擬,結果表明中心通道堵流包殼溫度高于角通道溫度[12]。2020 年,呂科鋒等對帶繞絲19 棒束鉛鉍多孔介質堵流進行了模擬,得出包殼存在圓周溫度變化[13]。上述研究均限于鉛鉍冷卻單相流,未涉及兩相流的流場計算。
對于鉛基堆而言,通常采用氬氣作為保護氣體或在加強自然循環中提供驅動壓力[14-16],因此存在鉛/鉛鉍-氬氣兩相流工況[17-19]。本文對于液態金屬-氬氣兩相流的數值模擬主要關注如下兩個方面的研究。一方面,對于加強自然循環工況[20-23],關注的重點在于鉛鉍空泡份額與氬氣注入速率之間的關系[21,23]。另一方面,在泵驅動的強迫循環工況下,回路運行前需先用氬氣將空氣排空,再注入氬氣將熔化的液態鉛鉍打壓進入回路。組件復雜結構內的部分區域難以被鉛鉍充滿,殘留的氬氣會滯留在回路內并形成局部“死區”,可能因導熱不良引發局部過熱,如有堵流發生,可能會加劇局部溫升。
因此,本文通過數值模擬給出鉛鉍-氬氣兩相下的氣泡行為,進而分析堵區局部過熱和壓降特征。
通常,燃料組件內液態鉛鉍的單相流動可以視為不可壓縮流體,滿足連續性、動量守恒和能量守恒方程。文獻[33]針對采用氣泡提升泵加強自然循環工況,對比了ε型及ω型湍流模型的相對誤差,得出standardk-ε計算誤差較小且節省計算資源,因此本文采用standardk-ε湍流模型,詳細公式可以參見ANSYS FLUENT theory guide[24]。對于多孔介質堵流,堵塞區域要增加考慮動量源項,見式(1) ~式(3)。其中,Sporous是多孔區域動量源項。
式中,下標hetr.可以是i,j,k三個方向的分量,表示速度在空間x,y,z上呈各向異性變化。多孔區域內的動量損失由黏性損失[即式(1)第一項]和慣性損失[即式(1)第二項]構成[25]。式(1)中變量β為不透水性,定義見式(2)。
式中,dpar是顆粒直徑,mm。 是孔隙率。該方法首次由Ergun 于1952 年提出[26],并且可以應用于鉛鉍多孔介質模擬仿真[17]。式(1)中慣性損失常數定義為摩擦系數和水利直徑的比值,如公式(3)所示。
針對液態鉛鉍低普朗特數特性,即液態鉛鉍與水相比具有更大的熱邊界層[27]。在ANSYS FLUENT 軟件中,默認的對流換熱和湍流普朗特數是針對液態水進行求解的,因此在求解鉛鉍介質時會帶來誤差。Lyu Kefeng等針對P/D為1.14、帶繞絲的鉛鉍冷卻燃料組件堵流工況,采用如公式(4)[28]所示的湍流普朗特數(Prt)關系式,獲得了較為合理的分析結果[13]。
兩相流的處理方法分為兩類:一類是均相流模型,適用于氣相含量比較低,并且兩相相對速度不大的情況;另一類是分相流模型,即對氣液兩相分別進行處理,對每一相計算其平均物理參量。ANSYS FLUENT 提供的多相流處理方法有VOF(Volume of Fraction)模型、混合物(Mixture)模型和歐拉(Eulerian)模型,前兩種模型屬于均相流模型,而Eulerian 模型屬于分相流模型。鑒于泵驅動的強迫循環和加強自然循環工況下,氬氣含量很少,鉛鉍-氬氣相對速度不大,因此可以采用VOF 模型,采用CSF(Continuum Surface Force)方法模擬表面張力[29]。
本文以德國KIT 鉛鉍冷卻燃料組件實驗流動傳熱數據為基準[9,32],進行數值計算結果的比對驗證。實驗裝置側視圖及堵塊設置如圖1 所示,燃料組件主要參數見表1。堵流實驗采用電加熱的均勻面熱源,熱功率為394±4 kW。鉛鉍入口溫度為200±0.2℃,質量流量為18.7±0.2 kg·s-1[32]。

圖1 燃料組件實驗裝置側視圖(黃色矩形區域為堵塊位置) [11]Fig.1 Side view of the testing platform (blockages are highlighted in yellow) [11]

表1 燃料組件主要參數[32]Table 1 Main parameters of the fuel assembly[32]
本課題組已完成了正常運行工況和全堵工況的流動傳熱及網格無關性驗證,結果顯示,溫度的模擬結果與實驗相比,最大相對誤差在±4.9%,摩擦系數與實驗推薦經驗關系式的相對誤差在±9.2%[33]。相對誤差的計算方法為:(數值模型計算值-實驗參考值)/實驗參考值。
基于表1 的燃料組件結構,本文采用FLUENT 軟件中自帶的FLUENT Meshing 模塊,針對表1 的中心6 子通道堵流(簡稱C6工況)和邊1 子通道堵流(簡稱E1工況)兩種工況,分別進行幾何結構建模及網格劃分,并在堵塊上游設置一個圓形氣泡的起始位置,如圖2所示。

圖2 堵塊及氣泡起始位置局部網格劃分Fig.2 Mesh display for the blockage regions and gas bubbles
加熱棒的壁面邊界條件采用與圖1 基準實驗一致的均勻熱流密度,對于小堵塊E1工況,為0.93 MW·m-2;對于大堵塊C6工況,考慮降功率保護,設置為E1工況的1/4,即0.23 MW·m-2。鉛鉍入口條件采用不同雷諾數下的質量流量進口,溫度為473.15 K;出口處設置為壓力出口,壓強取0 Pa。氬氣氣泡的起始位置設置在堵塊區域前方10 mm 處,直徑為1.0~2.0 mm(見表3 第2 列),并對氣泡可能流經的區域進行網格加密。本文采用FLUENT 軟件對流場進行初始化,在圓形氣泡區域通過patch 方式使其充滿氬氣。在瞬態計算時,設置連續性方程收斂條件取10-4量級,動量及能量方程收斂條件取10-5量級。考慮到流速和最小網格尺寸,時間步長取1×10-3s。
為驗證兩相流模型,本文同步通過速度擬合法[34]和氣泡穩定形態圖示法[16]進行驗證。其中,前者屬于經驗關系式,后者屬于經驗圖示,均來源于對大量實驗數據和結果的總結歸納。而本文數值模型采用VOF 方法屬于兩相求解方法,其與連續性方程、動量方程和能量方程一并用于兩相速度、能量的計算求解。采用VOF 方法得到的數值模擬結果若與經驗圖示及經驗關系式的求解結果相一致,則表明VOF 方法用于鉛鉍-氬氣工況計算是適合的。
對于速度擬合法,研究指出,通過計算無量綱數Fr,Eo,Mo和Nf數(Nf數為Eo、Mo的組合)可以得到兩相間作用力的影響情況。當Eo>70 時,表面張力可以不考慮;當Nf>550時,黏性作用影響很小;當Fr<0.05 時,慣性作用可以不考慮[34]。根據ANSYS 理論手冊[24],對于Re>>1 的工況,We數(慣性力和張力之比)不可忽略,當We>>1 時,張力影響可以忽略。
這些無量綱數的表達式(見表2 第2 行)與本文工況點1~6(見圖4)的計算結果各參數列于表2 第3 行。由計算結果可知,慣性力和張力的影響是主要影響因素,慣性力占主導。當慣性力的影響占主導時,根據文獻[34]可知,垂直管內流體速度應當滿足u=Const.×Ul+Ug,c,其中Ul為液相流速,m/s;Ug,c為最終穩定時氣泡頭部中心線上的速度,m/s。對于湍流,Const.的經驗值取1.2。由于VOF 模型只能獲得流體速度u,因此,將流體速度u與鉛鉍入口速度Ul進行擬合(見圖3),得到Const.約為1.05,該值稍小于經驗值[34],相關系數R2為0.9997,這與其他研究者的結論一致,因此可以認為模擬的結果是可以接受的。

圖3 慣性力主導下流體速度和液相速度擬合結果Fig. 3 The fitting curve for inertial force dominated velocities between two-phase flow and the liquid LBE

圖4 本模型工況點1~6 氣泡穩定上升時的形態Fig.`4 The numerical results of rising bubbles for E1 blockage conditions

表2 無量綱數的表達式及本模型中計算結果Table 2 Expressions for dimensionless numbers and calculation results
對于氣泡穩定形態圖示法[30,31],本文選取邊1 子通道堵流工況下的6 個工況點作為研究對象,其中1~2 研究氣泡直徑的影響,3~6 研究雷諾數的影響,工況點各參數取值見表3。網格加密區域的兩相流流速和氣相雷諾數由FLUENT 求解器輸出,分別見表3 第4 列及第5列。Mo數、Eo數、Re數定義式分別見式(5) ~式(7)[16],通過計算獲取這三個無量綱數列于表3后3 列之中。由FLUENT 軟件計算后,獲得氣泡隨無量綱數變化穩定時的形態如圖4 所示。考慮到計算結果的可靠性,選取Grace 經驗圖示(見圖5)作為基準形態,將表3 參數對照圖5 坐標,畫出工況點1~6 各點位置,見圖5 黃線所示。根據圖5 形態劃分,工況點1 落在球狀(Spherical)區域,其余5 個工況點均落在抖動狀(Wobbling)區域。對比圖4 中本模型中輸出的穩定形態,以及圖5 中氣泡穩定時的理論形態,可知本論文采用VOF 方法的模擬結果與經驗結果相吻合。

圖5 Grace 經驗圖示下的理論形態Fig.5 The stable shapes for gas bubbles under Grace correlation

表3 氬氣氣泡在液態鉛鉍中穩定上升時各參數取值 (E1 工況下,孔隙率0.2)Table 3 Values of parameters for steady rising bubbles for E1 blockage
本節考慮鉛鉍夾帶氬氣工況下,基于兩相流模型驗證結果,對C6工況和E1工況的局部過熱和壓降分別進行計算分析,并對入口雷諾數、孔隙率及氣泡直徑進行了參數敏感性分析。堵塊邊界條件設置分為兩種:一種是以外邊界為壁面來模擬實心堵塊情況;另一種是外邊界為內部面并設置孔隙率來模擬多孔介質堵塊情況。
3.1.1 入口雷諾數的影響
(1) 堵塊外邊界為壁面條件
本文通過設置氣相等體積分率的監測面,觀察含氣區域因導熱不良造成的溫升。選取組件軸向位置348~552 mm 的體平均溫度為參考溫度,堵塊在492~546.6 mm 區域,相對于參考節段位置為70.6%~97.4%,氣泡直徑保持2 mm不變。圖6 給出了不同雷諾數下氣相等體積分率(VOFAr)在[0.6~0.9]范圍內的過熱情況。由圖6 可知,隨流動時間的增加,VOFAr較高的等值面不一定存在連續性。其中圖6(a)至圖6(b)表明雷諾數較低時,VOFAr為0.9 的等值面大概率能夠在0.08 s 內監測到;而雷諾數較高時,且VOFAr值較大時,氣泡存在時間較短,僅約為0.03 s,如圖6(c)所示。且VOFAr值越高,引起的局部過熱越明顯。壁面邊界與內部面邊界相比,在瞬態條件下氣相體積分率可以保持在較高的水平(CAr>0.5),從而導致引起的過熱更嚴重。

圖6 不同雷諾數下氣相等值面溫度變化(C6 工況,外邊界為壁面)Fig.6 Temperature of iso-surfaces for gas phase at C6 blockage conditions(wall boundary conditions)
圖7 給出了當流動時間為0.01 s,Re為15000 時入口條件下的局部VOFAr分布及溫度云圖。由圖6 的計算結果可知,盡管VOFAr監測面溫度均值約在550 K,尚在可接受的范圍,然而圖7 局部最高溫度接近1000~1500 K,這可能對局部結構完整性造成威脅。

圖7 Re =15000 時堵區局部過熱,t=0.01s (C6 工況)Fig.7 Local overheating for diff. Re of C6 blockage conditions
圖8 給出了中低雷諾數下氣泡被封閉的現象,由圖可知雷諾數較低時,氣泡可能會卡在堵塊附近形成“死區”,并伴隨明顯的過熱。圖9 顯示了高雷諾數下氣泡逸出現象,由于在高雷諾數下高速流動的液態鉛鉍可能將氣泡帶出堵流區域,引起氣泡逃逸,從而使高含氣量的維持時間很短。即便在0.1 s 后再次監測到高等值面,逸出氣體引起的局部過熱也微乎其微。

圖8 中低雷諾數下氣泡被封閉現象(C6 工況)Fig.8 The sealed bubble for low and immediate Re of C6 blockage conditions

圖9 高雷諾數下氣泡逸出現象(C6 工況)Fig.9 The escape bubbles for high Re of C6 blockage conditions
(2)堵塊外邊界為內部面條件
在該類邊界條件下,氣泡可以貫穿多孔堵流區域,最終逃逸出該區域。為了研究不同入口Re影響,本文取孔隙率為0.8 和氣泡直徑為2 mm 且保持不變,觀察不同時刻氣相等值面溫度變化。
圖10 給出了不同雷諾數下VOFAr等值面溫度計算結果,由圖示可知,高VOFAr值面存留的時間非常短,且低VOFAr值時溫度隨流動時間近似線性升高。例如,Re取5000 的條件下,VOFAr取0.8 的等值面存留時間僅為50 ms,如圖10(a)所示。由圖10(b)可以看出,隨著雷諾數增加,VOFAr等值面存在時間縮短,且最終引起的溫升有所降低。然而,堵塊在內部面邊界條件下引起的過熱并不顯著。例如,流動時間為0.15 s 時,以堵塊所在節段的流體的體平均溫度為參考,在Re為5000 時VOFAr取0.6 的等值面過熱為5.4 K,而壁面條件的相同雷諾數下,其過熱溫度已達22.6 K。

圖10 不同雷諾數下VOFAr 等值面溫度計算結果(C6 工況,外邊界為內部面)Fig.10 Temperature of iso-surfaces for gas phase at C6 blockage conditions (internal boundary conditions)
為監測氬氣及堵流區域的靜壓瞬態變化,研究者在FLUENT 軟件中創建了一條監測線貫穿多孔堵流區域。該直線的x,y坐標與氣泡起始位置的坐標保持一致,得到不同雷諾數下沿軸向方向靜壓變化曲線。圖11 顯示隨著堵塊位置軸向坐標的增加,高VOFAr區域可能引起局部靜壓的小幅上升(微正壓現象)。特別是在氣泡進入多孔區域瞬間,引起的相對靜壓增量約為0.3%。由圖11(a) ~圖11(b)的變化可知,雷諾數越低,微正壓現象越明顯。這一現象從兩相流穩定條件出發可以解釋為:氣泡在液體中平衡存在除了需要具有一定的過熱度外,氣相和液相壓差還需滿足pg-pl=2σ/r*條件,其中pg為氣泡內的壓力,Pa;pl為液相壓力,Pa;σ為表面張力,N/m;r*是界面曲率,1/m[35]。因此,對于不透明介質,可以利用這一現象來監測氣泡的位置。

圖11 不同雷諾數下兩相區域壓降計算結果,Re 5000~15000(C6 工況)Fig.11 The pressure drops of diff. inlet Re for C6 blockage conditions
3.1.2 孔隙率的影響
為研究孔隙率對局部過熱和壓降特征的影響,本文鉛鉍入口雷諾數取臨界Re,即過渡流向湍流轉變的RebT=15700,且氣泡直徑取2 mm 保持不變,計算得到不同時刻孔隙率[0.2,0.8]的VOFAr值監測面溫度曲線,如圖12 所示。

圖12 不同孔隙率下溫度計算結果(C6 工況,Re=15700,氣泡直徑為2.0 mm)Fig.12 Temperature of iso-surfaces for gas phase at C6 blockage conditions
由圖12(a)可知,孔隙率取0.2 時,等值面溫度比該節段流體的體平均溫度過熱顯著,溫度曲線和全堵工況趨勢類似。這是由于氣泡被限制在堵塊局部形成 “死區”[見圖13(a)],與全堵工況時形成的“死區”相比,低孔隙率下氣泡并未完全被限制,而是分裂成一個較大的主氣泡和一個較小的子氣泡。主氣泡未被限制而逃逸,而子氣泡被限制在堵塊區域,如圖13(b)所示。

圖13 孔隙率0.2 時氣泡行為,t =0.1 s(C6 工況)Fig.13 The behavior of gas bubbles when porosity is 0.2 at C6 blockage conditions,t =0.1 s
對于孔隙率為0.4[見圖12(b)]的情況,在0.03 s 后僅低VOFAr值能夠被監測到,這種工況下,氣相也會引起比較明顯的過熱。氣相隨流動時間增加分裂成多個子泡,緩慢逃逸出多孔堵流區域。當孔隙率增至0.6 時,除了可以觀察到氣泡逃逸[圖14(a)]之外,氣相還會隨流動時間增加緩慢耗散在多孔堵流內部,如圖14(b)所示。當孔隙率增至0.8 時,氣泡能夠完整地隨流體從堵塊出口逸出多孔介質區域,并不會引起明顯的過熱現象。

圖14 孔隙率0.6 時氣泡逃逸和耗散現象(C6 工況)Fig.14 The behavior of gas phase when porosity is 0.6 at C6 blockage conditions
3.1.3 氣泡直徑的影響
為探究氣泡直徑的影響,研究者分別設置直徑為2.0 mm、1.5 mm 和1.0 mm 的圓形氣泡,其余參數取Re=RebT=15700,由于小孔隙率可能引起較為嚴重的過熱,故選取孔隙率為0.2 并對所有工況保持一致。
圖15 為不同氣泡直徑下VOFAr值監測面溫度計算結果。由圖15(a)可知,僅在直徑為2.0 mm 時,高VOFAr值監測面的存留時間較長,在0.02 s 內引起的過熱較為明顯。圖15(b)顯示直徑為1.5 mm 的氣泡在0.015 s 左右存在約為2 K 的短期過熱。這是由于僅直徑為2.0 mm 分裂出子泡并伴隨有“死區”形成(見圖13),其他氣泡直徑下均未分裂出子泡。因此,在鉛鉍-氬氣的夾帶工況中,應當關注氣泡尺寸(直徑≥2.0 mm)較大的情況;在氬氣加強自然循環工況中應當合理控制氣泡的注入尺寸。

圖15 不同氣泡直徑下氣相溫度計算結果(C6 工況,孔隙率0.2,Re=15700)Fig.15 Temperature of gas iso-surfaces for different bubble diameters( =0.2,Re =15700,C6 blockage conditions)
3.2.1 入口雷諾數影響
與上述探討C6工況時控制變量法的思想類似,為研究鉛鉍入口雷諾數影響,研究者控制氣泡直徑為2.0 mm,孔隙率為0.8,設置入口雷諾數以5000 為增量遞增,變化范圍為5000~40000,觀察堵流工況下氣泡行為、壓降及過熱特征,結果如下。
(1) 堵塊外邊界為壁面條件
在邊1 子通道堵流(E1工況)下,氣泡相較中心通道堵流有更大的逃逸趨勢。圖16(a)(b)分別顯示了Re為10000 和15000 時氣泡的逃逸現象。模擬結果顯示,在Re≥10000 時,氣體在多孔區域內的存留時長不超過0.01 s 便緩慢逸出。因此,除了Re為5000 時氣相的過熱較為明顯外,其余雷諾數下氣相引起的過熱皆不明顯。

圖16 不同雷諾數下氣泡逃逸現象(E1 工況,外邊界為壁面)Fig.16 Different behavior of gas bubbles for different inlet Re of E1 blockage
(2)堵塊外邊界為內部面條件
在E1工況下,堵塊位于軸向710.7~765.3 mm段,參考節段為組件加熱段軸向696~870 mm時,堵塞位置所在該段8.4%~39.8%的中上游位置。圖17 顯示了Re為20000 的條件下,堵塊區域靜壓的變化情況。由圖示可知,當流動時間由0.01 s 增至0.03 s 時,高氣相含量區域存在與中心6 子通道工況類似的微正壓現象(圖中箭頭所指位置為氣泡位置)。

圖17 Re=20000 時靜壓變化情況(E1 工況,外邊界為內部面)Fig.17 Pressure changes for E1 blockage,Re=20000
溫升方面,在0.03 s 內,VOFAr等值面溫度近似指數上升,然而其引起的過熱并不明顯,因此本文不進行展開闡述。
3.2.2 孔隙率影響
邊通道堵流下,孔隙率大小可能會對氣泡在多孔介質內的存留時長造成影響。在圖18中,研究者對不同孔隙率下氣泡是否逃逸進行了對比,由模擬結果可知,當孔隙率不太大( <0.8)時,氣泡容易在短期內逃逸,VOFAr等值面溫度上升到一個峰值溫度然后下降;當孔隙率比較大( ≥0.8)時,氣泡不易從多孔介質內逸出,因此能夠被加熱較長的時間,VOFAr等值面溫度近似指數形式上升,溫升也比小孔隙率時更為顯著。

圖18 不同孔隙率氣泡行為,Re=RebT,氣泡直徑為2.0 mmFig.18 Different behavior of gas bubbles for different porosities,diameter 2.0 mm
3.2.3 氣泡直徑影響
流動時間為0.01 s 時,氣泡直徑分別為1.0 mm、1.5 mm 和2.0 mm 的氣相參數和兩相流動的模擬結果已經在模型驗證部分進行了闡述(見圖4)。對于不同直徑,氣泡皆能夠存留在多孔區域較長時間,但引起的溫升均在1 K 以內,并不顯著。隨流動時間的增加,不同直徑的氣泡緩慢耗散在液態鉛鉍中,VOFAr監測值也逐漸降低。
本文對帶繞絲的19 棒束六邊形燃料組件中心通道和邊通道堵流工況下,鉛鉍-氬氣夾帶兩相流的局部過熱和壓降特征進行了闡述,主要結論如下:
氣泡行為可以概括為三類,即逃逸、耗散以及受限形成“死區”。入口雷諾數、孔隙率、氣泡直徑三者共同作用對氣泡行為造成影響,進而影響堵塊區域的傳熱和壓降特征。
(1)孔隙率較小,入口雷諾數較大時,較易發生逃逸現象。邊子通道堵流條件下發生的可能性大于相同參數下的中心通道堵流。
(2)在孔隙率較大時,較易發生耗散現象,可能伴隨氣泡逃逸。
(3)鉛鉍入口雷諾數較小(Re≤RebT),孔隙率較小時,較易發生受限現象。直徑較大(2 mm)的氣泡可能分裂出一個小的子泡,引起受限過熱。
在過熱特征方面,氣相更易在堵塊附近存留形成高體積分率區域,并引起過熱。與堵塊外邊界為內部面時相比,堵塊外邊界為壁面時所引起的過熱現象更為顯著。
在壓降特征方面,高VOFAr區域會引起局部微正壓現象,造成堵塊軸向監測線上靜壓的小幅上升,微正壓的出現位置與流動時間和鉛鉍入口雷諾數有關。
致謝:作者對所有為本文工作提供幫助指正的人員表示感謝。本文是在科技部國家重點研發計劃,Grant No. 2022YFB1902503 項目的支持下完成的,本文數值模擬得到了合肥先進計算中心的支持,在此特別表示感謝。