[發明專利]一種核反應堆板狀燃料熔化流固耦合無網格分析方法

有效

• 申請號: 202110484617.5
• 申請日: 2021-04-30
• 公開（公告）號: CN113192567B
• 公開（公告）日: 2022-12-09
• 發明（設計）人: 陳榮華;肖鑫坤;蔡慶航;張鵬輝;郭凱倫;田文喜;蘇光輝;秋穗正
• 申請（專利權）人: 西安交通大學
• 主分類號: G16C10/00
• 分類號: G16C10/00;G16C20/10;G16C20/90;G06F30/25;G06F30/28;G06F113/08;G06F119/14
• 代理公司: 西安智大知識產權代理事務所 61215
• 代理人: 何會俠

1.一種核反應堆板狀燃料熔化流固耦合無網格分析方法，其特征在于：步驟如下：

步驟1：對核反應堆板狀燃料熔化流固耦合問題進行粒子建模；用不同類型的粒子代表核反應堆板狀燃料的各個組成部分，0號粒子代表鋯合金材質的燃料包殼和燃料基體，1號粒子代表彌散在燃料基體中的二氧化鈾燃料微粒；對每一個粒子進行編號，所有的粒子均具有對應的物性參數，物性參數包括質量、密度、比熱、熔沸點、溫度、焓和初速度；粒子在笛卡爾坐標系內，依據核反應堆板狀燃料的具體參數均勻布置，能夠模擬核反應堆板狀燃料的初始狀態；

步驟2：使用鏈表法檢索鄰居粒子，具體方法是：將整個計算域均勻劃分成正方形表格，正方形表格的邊長與粒子最大作用半徑相等；對于任意一個粒子而言，檢索鄰居粒子時只需檢索包含粒子在內的表格以及包圍此表格的共計9個表格，三維計算時，需要檢索的就是由正方形表格圍成的正方體，需要檢索27個三維正方體；

步驟3：對時間步長進行穩定性分析，采用Courant-Friedrichs-Lewy條件簡稱CFL條件判斷；CFL條件判斷是差分方法保持計算穩定的必要條件，對于粒子法也適用，在移動粒子半隱式方法中采用下式進行判斷

式中：

Δt——能夠保持計算穩定的最小時間步長；

C——庫朗常數；

l₀——粒子初始布置時相鄰兩粒子間的距離；

u_max——各時間步長中粒子速度最大值；

計算輸入的初始條件中包含預設的時間步長大小，將輸入的時間步長和使用CFL條件計算得到的能夠保持計算穩定的最小時間步長進行比較，使用較大的時間步長作為當前時間步長；

步驟4：使用核函數表征粒子之間的相互作用程度，具體為在對動量守恒方程中的各項進行顯式計算時，先使用核函數進行離散然后積分得到動量守恒方程中各項的具體值，使用的核函數如下所示：

式中

w(r_ij)——核函數值；

r_ij——粒子i與粒子j之間的距離；

r_e——粒子作用域半徑；

步驟5：對核反應堆板狀燃料熔化后的熔融物研究對象采用如下控制方程：

式中

ρ——粒子對應物質的密度kg/m³；

t——時間s；

P——壓力Pa；

μ——動力粘度系數N·s/m²；

——速度矢量m/s；

——表面張力矢量N/kg；

——重力加速度m/s²；

公式(3)為連續性方程，針對的對象為不可壓縮流體，公式(4)為動量守恒方程，從中看出分析對象只受到重力這一外力作用，對于任意一個粒子，通過核函數先將動量守恒方程即公式(4)中的重力項、粘性項和表面張力項離散，然后對計算域內全部粒子進行積分，顯式得到粒子的估算速度和位置，其中粘性項需要使用采用如下多相粘度模型：

式中

μ——動力粘度系數N·s/m²；

μ_ij——粒子間動力粘度系數m²/s；

——速度m/s；

d——直徑m；

r_e——粒子作用半徑m；

n_G——使用高斯核函數計算得到的粒子數密度；

——粒子j速度m/s；

——粒子i速度m/s；

G——高斯核函數；

r_j——粒子j的位置矢量；

r_i——粒子i的位置矢量；

步驟6：使用如下的焓值相變模型計算粒子的相態：

式中

T——需要計算的粒子溫度K；

T_s——需要計算的粒子對應的熔點K；

h——當前粒子的焓值J/kg；

h_s0——粒子開始熔化的焓值J/kg；

h_s1——粒子結束熔化的焓值J/kg；

c_p——粒子等壓比熱容J/(kg·K)；

通過步驟6的計算，能夠模擬核反應堆板狀燃料在最高溫度超過燃料基體鋯合金的熔點時的相變過程；計算得到每個粒子在不同時刻下的種類、焓值和溫度，即得到鋯合金基體和鈾燃料微粒的相態、焓值和溫度隨時間的變化過程；

步驟7：得到粒子的相態后，通過非連續介質模型的分析方法離散單元法DEM對固體之間的碰撞作用進行如下的計算：

非連續介質模型的分析方法離散單元法DEM的計算對象為剛性固體，當計算對象中存在大固體時，需要先解決將離散顆粒聚集成剛性物體統一運動的問題；采用如下公式初始化大固體運動速度；

式中

——標記為ii的剛體的速度矢量m/s；

——粒子i當前的速度矢量m/s；

——標記為ii的剛體的角速度矢量rad/s；

m_i——粒子i的質量kg；

——粒子i的位置矢量；

——標記為ii的剛體的位置矢量；

I_ii——標記為ii的剛體的轉動慣量；

對于任意固體單元m_ii，采用如下公式對單元受力作用后產生的運動進行分析；

式中

——標記為ii的剛體的加速度矢量；

——標記為ii的剛體的角加速度矢量；

——標記為ii的剛體的速度矢量；

——標記為ii的剛體的角速度矢量；

——標記為ii的剛體所受合力矢量；

——標記為ii的剛體所受合力矩矢量；

公式(8)和公式(9)中固體單元所受的合力與和力矩通過下式計算

式中

F_X,jj→ii——標記為jj的剛體對標記為ii的剛體施加的法向合力；

e_n——法向彈性力；

d_n——法向阻尼力；

k_n——法向彈性系數；

c_n——法向阻尼系數；

F_Y,jj→ii——標記為jj的剛體對標記為ii的剛體施加的切向合力；

e_s——切向彈性力；

d_s——切向阻尼力；

k_s——切向彈性系數；

c_s——切向阻尼系數；

這樣，固體單元m_ii就移動到了一個新的位置，并產生新的接觸力和接觸力矩，再次計算固體單元m_ii所受的合力和合力矩重復計算流程，即得到每個固體單元及整個散粒體的運動形態，此時固體粒子的速度和位置均為估算值；

步驟7：通過核反應堆板狀燃料熔化后的熔融物的不可壓縮特性和動量守恒方程中的壓力項聯立得到如下壓力泊松方程，隱式求解如下的壓力泊松方程

式中

ρ——粒子對應物質密度kg/m³；

P——壓力Pa；

n⁰——初始粒子數密度常數；

n^*——臨時粒子數密度；

n^k——當前時層的粒子數密度；

n^k-1——上一時層的粒子數密度；

n^k-1——下一時層的粒子數密度；

β——壓力源項修正系數；

γ——壓力源項修正系數；

通過步驟8的計算得到動量守恒方程中的壓力項，使用該壓力項所得值對步驟5和步驟7中得到的粒子估算速度和位置進行修正，得到當前時間步長全部粒子的真實速度和位置，這樣就完成了一個時間步長的核反應堆板狀燃料熔化流固耦合的模擬；

步驟9：為了準確分析模擬對象中固體和液體之間的耦合作用，還需要對步驟5中的粘性項計算和步驟8中的壓力項計算進行處理，這是因為，固體和液體之間的作用通過下式表現

式中

——流體對固體的作用合力；

——流體對固體的壓力作用；

——流體對固體的粘性作用；

在步驟5的粘性項計算中，通過多相粘度模型中使用的調和粘性系數來表現固體的物性，在計算時，這個粘性項只發生在固體和液體接觸時，對固體內部不進行粘性項的計算；在步驟8的壓力項計算中，流體粒子只和接觸的固體粒子施加壓力作用，后續壓力作用將在非連續介質模型的分析方法離散單元法DEM的計算中擴展至全部固體粒子；

綜上，通過步驟1對核反應堆板狀燃料熔化流固耦合問題中，鋯合金材質的燃料包殼和燃料基體和二氧化鈾燃料微粒的位置、速度和初始物性參數進行設定；通過步驟6模擬板燃料的鋯合金基體在溫度超過其熔點后發生相變形成熔融物，計算得到每個粒子在不同時刻下的種類、焓值和溫度，得到熔融物的相態、焓值和溫度隨時間的變化過程；通過步驟5至步驟9模擬核反應堆板狀燃料熔化流固耦合過程，模擬熔融物和未熔化的燃料微粒之間的相互作用，模擬燃料微粒之間的相互作用過程，計算得到熔融物和燃料微粒粒子的速度、位置和壓力，即得到熔融物遷徙過程中的運動和壓力變化過程；綜合以上步驟，分析了核反應堆板狀燃料熔化流固耦合問題，得到熔化流固耦合過程中熔融物的位置、速度、壓力、相態、溫度、焓值隨時間的變化，通過以上步驟對熔融物的遷徙行為展開機理性分析。