[發明專利]基于粒子法的核反應堆燃料早期變形分析方法

有效

• 申請號: 202110486476.0
• 申請日: 2021-04-30
• 公開（公告）號: CN113191065B
• 公開（公告）日: 2022-07-26
• 發明（設計）人: 陳榮華;蔡慶航;肖鑫坤;王金順;李勇霖;田文喜;蘇光輝;秋穗正
• 申請（專利權）人: 西安交通大學
• 主分類號: G06F30/25
• 分類號: G06F30/25;G16C10/00;G16C20/10;G06F119/08;G06F119/14
• 代理公司: 西安智大知識產權代理事務所 61215
• 代理人: 何會俠

1.一種基于粒子法的核反應堆燃料早期變形分析方法，其特征在于：步驟如下：

步驟1：根據實際的核反應堆燃料結構建立粒子幾何模型，定義材料數M區分不同的材料，M＝0,1,2分別表示核反應堆燃料元件的芯塊、包殼和結構材料，粒子根據對應的材料攜帶對應的物性和力學參數；定義ID數區分固體邊界，采用1號粒子模擬固體邊界，2號粒子模擬固體內連續介質區域；不同固體之間采用參數Object區分，對于2號粒子，只有當Object數相同時才發生相互作用；初始化所有粒子的速度、位置、溫度、應力張量和應變張量，分別為u_i、r_i、T_i、σ_i、ε_i；通過以上參數設定，建立燃料芯塊、包殼與結構材料三者之間的固體邊界，還原真實的燃料結構，并確定了燃料的初始位置、速度、溫度、應力應變狀態和相關物性參數；

步驟2：根據建立的初始粒子幾何模型，為每個粒子建立固體內粒子鄰居域，即存儲每個粒子以自身為圓心，作用范圍r_e為半徑的區域內所有Object數與該粒子相同的粒子ID，固體內鄰居域不隨計算時間的推移而發生改變，即認為固體內連續介質發生相互作用的粒子始終保持不變；除了固體內粒子鄰居域，還需要為所有1號粒子建立固體間粒子鄰居域，即存儲每個1號粒子以自身為圓心，作用范圍r_e為半徑的區域內所有Object數與該粒子不同的1號粒子ID，固體間鄰居域會隨著計算時間的推移而發生改變，即兩個固體的相對位置會不斷發生變化，當兩者表面距離小于作用范圍r_e時，兩個固體的邊界會發生相互作用，否則不會；定義粒子數密度為n，其值為粒子鄰居域內所有粒子對中心粒子的核函數值之和，如公式(1)所示；

式中

n_i——粒子i的粒子數密度；

w(r_ij,r_e)——核函數，采用的核函數形式為

r_ij——粒子i和粒子j的距離，m；

r_e——粒子作用范圍，m；

下標i——任意粒子i；

下標j——粒子i的粒子鄰居域內的任意粒子j；

由此，進一步定義固體內粒子數密度n_inner為粒子i的固體內鄰居域的粒子數密度、粒子間粒子數密度n_outer為粒子i的固體間鄰居域的粒子數密度、n⁰為參考粒子數密度，是假設粒子均勻分布情況下計算得到的粒子數密度；

步驟3：基于邊界條件和固體碰撞模型計算固體邊界粒子外力載荷，根據具體實際物理問題，為1號粒子設置外力的作用方向和大小，由公式(2)計算邊界條件外力載荷的影響，邊界條件外力載荷始終視為作用于粒子中心，因此不考慮該力所引起的粒子轉動，公式采用張量表示法和下標表示法；

式中

v_α——α方向上的速度分量，m/s；

——α方向上的邊界條件外力載荷引起的加速度，m/s²；

t——時間，s；

對于固體碰撞過程，采用固體碰撞模型計算1號粒子的碰撞載荷，如公式(3)所示，

式中

P_i^Contact——粒子i的接觸應力，N/m²；

λ——拉梅常數的第一個參數，

(ε_γγ)_i——粒子i的應變張量對角項之和，即ε_γγ＝ε₁₁+ε₂₂+ε₃₃；

n_outer——粒子i的固體間鄰居域的粒子數密度；

n_inner——粒子i的固體內鄰居域的粒子數密度；

E——楊氏模量，N/m²；

υ——泊松比；

由碰撞載荷，計算對應的速度梯度，如公式(4)所示，

式中

v_i——粒子i的速度矢量，m/s；

t——時間，s；

ρ_i——粒子i的密度，kg/m³；

d——空間維度；

n⁰——參考粒子數密度；

r_ij——粒子i和粒子j的相對位置矢量，m；

P_i^Contact——粒子i的接觸應力，N/m²；

——粒子j的接觸應力，N/m²；

——粒子i和粒子j的相對初始位置矢量，m；

——核函數，等同于

通過步驟3的計算，模擬核反應堆燃料早期變形過程中燃料元件受到的外力載荷和碰撞力，全面考慮了核反應堆燃料早期變形過程中可能發生的外力作用、燃料元件之間的碰撞、燃料元件與結構材料之間的碰撞、燃料包殼與芯塊之間的碰撞、芯塊與芯塊之間的碰撞過程；計算得到每個粒子由外力載荷和碰撞引起的加速度，即燃料元件由外力載荷和碰撞引起的加速度；

步驟4：能量計算如公式(5)所示：

式中

h——焓值，J/kg；

t——時間，s；

ρ——密度，kg/m³；

k——熱導率，W/(m·K)；

T——溫度，K；

Q_radiation——輻射熱源，W/m³；

Q_out——外熱源，W/m³；

輻射換熱采用斯忒藩-玻爾茲曼定律，如公式(6)所示

式中

Q_radiation——輻射熱源W/m³；

ε——發射率；

σ——斯忒藩-玻耳茲曼常數；

T——溫度，K；

T_env——環境溫度K；

l₀——粒子直徑m；

溫度的拉普拉斯項采用拉普拉斯模型離散得到：

式中

d——空間維度；

n⁰——參考粒子數密度；

ρ_i——粒子i的密度，kg/m³；

——當前時刻的粒子j溫度K；

T_i^k——當前時刻的粒子i溫度K；

r_ij——粒子i和粒子j的相對位置矢量，m；

w(|r_ij|)——核函數，等同于w(|r_ij|,r_e)；

由焓值計算粒子溫度：

式中

T——溫度K；

T_s——熔點K；

h——焓值J/kg；

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

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

c_p——比熱容J/(kg·K)；

通過步驟4的計算，模擬核反應堆燃料早期變形過程中燃料元件內的導熱、燃料元件間的輻射換熱、芯塊和包殼間的換熱過程；計算得到每個粒子在不同時刻下的焓值和溫度，即得到核反應堆燃料的焓值和溫度隨時間的變化過程；

步驟5：熱膨脹應變計算如公式(8)所示，

式中

[ds^T]_i——粒子i的熱膨脹應力增量張量，N/m²；

κ_i——粒子i的熱膨脹系數；

T_i——粒子i的溫度，K；

T_i^ref——粒子i的參考溫度，K；

通過步驟5的計算，模擬核反應堆燃料早期變形過程中燃料元件的熱膨脹變形過程；計算得到每個粒子在不同時刻下的熱膨脹應變張量，即得到了燃料元件由于溫度變化而發生碰撞變形的行為規律；

步驟6：計算總應變，固體內的相對變形是由當前相對位置減去旋轉后的相對位置得到的，如公式(9)所示，

式中

u_ij——粒子i和粒子j的相對總變形矢量，m；

r_ij——粒子i和粒子j的相對位置矢量，m；

R——粒子i的旋轉矩陣，R＝R_xR_yR_z；

——粒子i和粒子j的相對初始位置矢量，m；

R_x——粒子i繞x軸旋轉的旋轉矩陣，

R_y——粒子i繞y軸旋轉的旋轉矩陣，

R_z——粒子i繞z軸旋轉的旋轉矩陣，

θ_α,ij——粒子i和粒子j繞α軸旋轉的相對角度，°，α可以取x，y，z；

由總變形可以計算得到總應變：

ε_ij——粒子i和粒子j的相對總應變矢量；

通過步驟6的計算，計算得到了每個粒子在不同時刻下的總應變張量，即得到了核反應堆燃料元件的總變形狀態；

步驟7：基于彈性模型計算彈性應力，彈性模型如公式(11)所示，

式中

——彈性應力張量的分量，N/m²；

λ——拉梅常數的第一個參數，

μ——拉梅常數的第二個參數，

——彈性應變張量對角項之和，即

δ_αβ——克羅內克爾函數，

——彈性應變張量的分量；

其中記各向同性力：

式中

P_i^E——粒子i的彈性各向同性力，N/m²；

λ——拉梅常數的第一個參數，

——彈性應變張量對角項之和，即

d——空間維度；

n⁰——參考粒子數密度；

u_ij——粒子i和粒子j的相對總變形矢量，m；

r_ij——粒子i和粒子j的相對位置矢量，m；

——粒子i和粒子j的相對初始位置矢量，m；

——核函數，等同于

此外，在碰撞計算中，如果碰撞載荷P_i^Contact小于彈性各向同性力P_i^E，則令P_i^Contact＝P_i^E；

記各向異性力：

式中

——粒子i和粒子j間的彈性各向異性力，N/m²；

μ——拉梅常數的第二個參數，

——粒子i和粒子j的相對彈性應變矢量；

由彈性各向同性力和彈性各向異性力分別求得對應的速度梯度，如公式(14)和(15)所示，

式中

v_i——粒子i的速度矢量，m/s；

t——時間，s；

ρ_i——粒子i的密度，kg/m³；

d——空間維度；

n⁰——參考粒子數密度；

r_ij——粒子i和粒子j的相對位置矢量，m；

P_i^E——粒子i的彈性各向同性力，N/m²；

——粒子j的彈性各向同性力，N/m²；

——粒子i和粒子j間的彈性各向異性力，N/m²；

——粒子i和粒子j的相對初始位置矢量，m；

——核函數，等同于

通過步驟7，模擬了核反應堆燃料早期變形過程中的彈性變形；計算得到每個粒子在不同時刻下的彈性應變張量、彈性應力張量、彈性力所引起的加速度，即得到燃料元件的彈性應力應變隨時間的變化過程；

步驟8：基于角動量定理計算粒子旋轉角速度，粒子j對粒子i的切應力為：

式中

F_ij——粒子j對粒子i的切應力，N；

d——空間維度；

n⁰——參考粒子數密度；

m_i——粒子i的質量，kg；

ρ_i——粒子i的密度，kg/m³；

——粒子i和粒子j間的彈性各向異性力在與r_ij垂直方向上的投影，N/m²；

r_ij——粒子i和粒子j的相對位置矢量，m；

——粒子i和粒子j的相對初始位置矢量，m；

——核函數，等同于

粒子j對粒子i的切應力產生的扭矩為：

T_ij＝-r_ij×F_ij 公式(17)

式中

T_ij——粒子j對粒子i的切應力產生的扭矩，N·m；

r_ij——粒子i和粒子j的相對位置矢量，m；

F_ij——粒子j對粒子i的切應力，N；

粒子i的角速度由公式(18)得：

式中

ω_i——粒子i的角速度，rad/s；

t——時間，s；

I_i——粒子i的轉動慣量，kg·m²；

T_ij——粒子j對粒子i的切應力產生的扭矩，N·m；

通過步驟8的計算，計算得到每個粒子的角速度，由角速度計算得到旋轉角度，從而計算步驟6中的旋轉矩陣R；該步驟保證了粒子的運動滿足角動量定律，即得到了燃料元件變形過程中的固體旋轉或扭轉隨時間的變化過程；

步驟9：總應變為彈性應變加塑性應變加熱膨脹應變，求解塑性應變基于應變增量理論；假設在計算塑性應變之前，先假設無塑性應變，由該假設得彈性應變計算等效應力，并將等效應力與屈服極限比較，如果等效應力大于屈服極限，則認為發生塑性變形，等效應力求解如公式(19)至公式(21)所示，

[ds]_n＝2μ[dε^e]_n 公式(19)

[s]_n＝[s]_n-1+[ds]_n 公式(20)

式中

[ds]_n——第n時間步下的應力增量張量，N/m²；

μ——拉梅常數的第二個參數，

[dε^e]_n——第n時間步下的彈性應變增量張量；

[s]_n——第n時間步下的應變張量，N/m²；

[s]_n-1——第n-1時間步下的應變張量，N/m²；

——粒子i的等效應力，N/m²；

s_ij——應力矢量，N/m²；

通過步驟9的計算，得到了每個粒子在不同時刻下是否發生塑性變形的判定，即判定燃料元件的不同位置是否發生塑性變形；

步驟10：塑性應力應變計算，塑性應變增量如公式(22)和(23)所示，

[dε]_n＝[dε^p]_n+[dε^e]_n+[dε^T]_n 公式(23)

式中

[ds^p]_n——第n時間步下的塑性應力增量張量，N/m²；

[ds]_n——第n時間步下的應力增量張量，N/m²；

[ds^T]_n——第n時間步下的熱膨脹應力增量張量，N/m²；

[ds^e]_n——第n時間步下的彈性應力增量張量，N/m²；

[s]_n-1——第n-1時間步下的應變張量，N/m²；

μ——拉梅常數的第二個參數，

dλ_n——增量參數，由材料的力學特性決定；

式中dλ_n的取值影響塑性應變和彈性應變的關系，不同特性的材料，其dλ_n的求解不同，以下針對線性彈塑性，dλ_n的求解如公式(24)所示，

式中

——假設所有應變均為彈性應變條件下的等效應力，N/m²；

H′——塑性硬化系數；

Y——屈服極限，N/m²；

——第n-1時間步下的等效塑性應變；

μ——拉梅常數的第二個參數，

由塑性模型獲得塑性應變后，由公式(23)計算得到彈性應變，再采用步驟7中的彈性模型計算彈性應力；

通過步驟10的計算，模擬了核反應堆燃料早期變形過程中的塑性變形過程；計算得到了每個粒子在不同時刻下的塑性應變和彈性應變，即得到了燃料元件的塑性變形和彈性變形隨時間的變化過程；

步驟11：斷裂判斷，當粒子間距大于或小于斷裂閾值時，則認為發生斷裂，斷裂后斷裂部分的2號內部粒子轉化為1號邊界粒子，且斷裂粒子之間不再發生固體連續介質內的相互作用，僅作為固體邊界參與邊界載荷和碰撞的計算；根據外力載荷、碰撞和彈性力更新粒子速度、位置、角速度和旋轉角度；

通過步驟11的計算，模擬了核反應堆燃料早期變形過程中的斷裂過程；計算得到了每個粒子在不同時刻下的是否發生斷裂的判定，即獲得了燃料元件發生斷裂的位置；

通過步驟1對實際的核反應堆燃料結構建立粒子幾何模型，設定粒子的位置、速度、溫度、物性參數和力學參數；通過步驟2存儲粒子鄰居域，為后續計算提供基礎；通過步驟3至步驟10計算固體的溫度和受力，受力包括外力、碰撞、彈性應力應變、塑性應變和熱膨脹應變，更新粒子的速度和位置；通過步驟11判斷粒子是否發生斷裂；由步驟3至步驟10計算得到的加速度，更新所有粒子的速度和位置；對計算是否結束進行判定，若是，輸出結果，若否，更新固體間粒子鄰居域，再返回步驟3；綜合以上步驟，模擬了核反應堆燃料早期變形過程，得到了變形過程中燃料元件的位置、速度、溫度、應力狀態、應變狀態，通過以上數據對核反應堆燃料早期變形過程中的傳熱、熱膨脹、彈性變形、塑性變形和受力情況展開機理性分析。