本發(fā)明涉及反應(yīng)堆堆芯設(shè)計和安全分析
技術(shù)領(lǐng)域：
，具體涉及一種穩(wěn)定精確的反應(yīng)堆物理熱工耦合計算方法，稱之為全耦合法。
背景技術(shù)：
：反應(yīng)堆是一個多物理場耦合的復(fù)雜系統(tǒng)，由于其復(fù)雜性，針對各個物理場形成了相對獨立的學(xué)科，其中核反應(yīng)堆物理分析以及核反應(yīng)堆熱工分析是反應(yīng)堆設(shè)計和安全分析的兩大重要領(lǐng)域。所謂物理熱工耦合，是指堆內(nèi)的物理參量間互相存在反饋效應(yīng)，這些反饋效應(yīng)在反應(yīng)堆發(fā)生瞬態(tài)安全性事故時會更加顯著，因此準確地考慮這些反饋效應(yīng)則需要物理熱工的耦合分析計算。目前廣泛采用的物理熱工耦合計算方法是一種簡單耦合，國際上亦稱之為“代碼耦合”或“松耦合”，物理方程和熱工方程的求解彼此是解耦進行的，其實現(xiàn)方式是通過外部腳本或代碼合并，將一個專門求解反應(yīng)堆中子學(xué)問題的程序和一個專門求解反應(yīng)堆熱工學(xué)問題的程序聯(lián)接，物理程序?qū)⒂嬎愕玫降墓β史植紓鬟f給熱工程序，熱工程序?qū)⒂嬎愕玫降臒峁⒘糠答伣o物理程序，這一過程中僅存在數(shù)據(jù)交互，而核心求解過程彼此獨立，互為“黑匣子”，在此基礎(chǔ)上，在瞬態(tài)分析計算中又分為顯示耦合方法和隱式耦合方法，其區(qū)別在于顯示耦合在一個時間步內(nèi)只進行一次物理計算和一次熱工計算，而隱式耦合需要進行物理熱工間的迭代，直到計算結(jié)果在一個時間步內(nèi)滿足一定的收斂判據(jù)?，F(xiàn)有的物理熱工耦合方法，最大的優(yōu)勢在于易于實現(xiàn)，因為它幾乎不需要程序再開發(fā)，然而，由于這種“松耦合”的算法本質(zhì)，其收斂速度有限，若要獲得高收斂精度的結(jié)果需要進行很大的物理熱工迭代次數(shù)，因此需要在計算精度和計算效率間平衡，此外現(xiàn)有耦合方法在某些情況下亦會出現(xiàn)計算震蕩甚至發(fā)散的問題，雖然通過增加物理熱工耦合計算間的阻尼因子可以顯著緩解該現(xiàn)象，但依舊不能跳出“代碼耦合”的思想框架，物理方程和熱工方程的依舊是解耦求解的。毫無疑問，物理方程和熱工方程的解耦是導(dǎo)致現(xiàn)有耦合計算方法存在問題的重要源頭，只有將物理熱工方程聯(lián)立，一次計算中同時求解中子學(xué)方程和熱工方程才能從根本上解決上述問題。技術(shù)實現(xiàn)要素：為了解決上述現(xiàn)有技術(shù)存在的問題，本發(fā)明的目的在于提供一種穩(wěn)定精確的反應(yīng)堆物理熱工耦合計算方法，全耦合法，該方法利用牛頓迭代法，下文簡稱為牛頓法，求解聯(lián)立的物理熱工耦合的非線性方程組，實現(xiàn)了物理熱工的耦合計算分析，是一種與傳統(tǒng)耦合方法迥異的新的耦合計算方法，理論上具有良好的數(shù)值穩(wěn)定性和數(shù)值求解精度，同時牛頓法的超線性收斂速度亦能保證計算效率。牛頓法求解非線性方程組的一般過程如下：設(shè)F(x)＝0為一非線性方程組，求F(x)＝0的解。給定初值x＝x0作do循環(huán)JF(xk)δxk＝-F(xk)xk+1＝xk+δxk當滿足收斂判據(jù)時循環(huán)結(jié)束其中，的雅閣比矩陣，k--牛頓迭代步數(shù)。為了實現(xiàn)全耦合方法，本發(fā)明采用如下技術(shù)方案：一種穩(wěn)定精確的反應(yīng)堆物理熱工耦合計算方法，包括如下步驟：步驟1：確定中子學(xué)方程，即時空中子擴散方程的離散形式，對時空中子擴散方程采用非線性迭代半解析節(jié)塊法進行空間上的離散，時間項采用隱式差分，能量取兩群近似，得到帶節(jié)塊耦合修正因子項的節(jié)塊平衡方程，即帶修正因子的CMFD方程，其表達式如下：1vgφgm,n+1-φgm,nΔtn+1=Σg′=1G(χpgβpνΣfg′m,n+1+Σg′gm,n+1)φg′m,n+1+χdgΣlλlClm,n+1-Σtgm,n+1φgm,n+1-Lgm,n+1,g=1,2,...,G]]>Cln+1=e-λlΔtnCln+Fln0Σg′=1G(νΣfg′φg′)n+Fln1Σg′=1G(νΣfg′φg′)n+1]]>Lgm=Σu=x,y,z1Δum(Jgu+m-Jgu-m)]]>Jgu+m=-Dgm,FDM(fgu-m+1φgm+1-fgu+mφgm)-Dgm,NOD(fgu-m+1φgm+1+fgu+mφgm)]]>式中，G——能群總數(shù)；Δt——時間步長；——第g群n時刻節(jié)塊m內(nèi)的平均中子通量；——第g群n+1時刻節(jié)塊m內(nèi)的平均中子通量；vg——第g群中子速度；χpg——瞬發(fā)中子裂變譜；βp——瞬發(fā)中子份額；χdg——緩發(fā)中子裂變譜；νΣf——中子產(chǎn)生截面；Σg'g——g’群到g群宏觀散射截面；Σt——宏觀總截面；χpg——瞬發(fā)中子裂變譜；λl——瞬發(fā)中子份額；Cl——第l群緩發(fā)中子先驅(qū)核濃度；——由先驅(qū)核濃度線性近似推導(dǎo)得到的常數(shù)項系數(shù)；——第g群節(jié)塊m內(nèi)的中子泄漏項；——第g群節(jié)塊m內(nèi)u方向左右表面的中子靜流；Δu——u方向的節(jié)塊寬度；——節(jié)塊等效擴散系數(shù)；——非線性迭代節(jié)塊耦合修正因子；——m節(jié)塊在u方向左側(cè)面的不連續(xù)因子；確定熱工方程的離散形式，對熱工方程采用單通道模型，只考慮徑向?qū)幔臻g上采用有限體積法離散，時間項采用隱式差分，得到有限差分格式的熱工離散方程：AΔXjΔt(ρj-ρjn)+mj-mj-1=0]]>ΔXjΔtAρj(hj-hjn)-mj-1(hj-1-hj)=HPrΔX(Tw-Tb)]]>(ρfcpV)iTi-TinΔt+(Ki+Ki-1)Ti-Ki-1Ti-1-KiTi+1=qi′′′]]>式中A——通道流動面積；ρj——第j網(wǎng)格冷卻劑密度；mj——第j網(wǎng)格質(zhì)量流量；hj——第j網(wǎng)格冷卻劑比焓；H——傳熱系數(shù)；Pr——釋熱周長；Tw——壁面溫度；Tb——冷卻劑溫度；ρf——燃料密度；cp——燃料比熱；V——傳熱網(wǎng)格體積；K——等效導(dǎo)熱率；q″′——體積釋熱率；步驟2：由于中子學(xué)方程與熱工方程在數(shù)學(xué)表達式中并無顯示的聯(lián)系，因此需要對求解的節(jié)塊平衡方程與熱工離散方程建立數(shù)學(xué)關(guān)系，在節(jié)塊平衡方程中截面是慢化劑密度和燃料溫度的隱式函數(shù)Σx＝g(ρ,Tf)，熱工離散方程中體積釋熱項是關(guān)于中子通量和裂變截面的函數(shù)從而根據(jù)這兩個關(guān)系，便能構(gòu)建封閉的中子學(xué)方程和熱工方程耦合的非線性方程組，形式如下：F(x)=fφ(x)fh(x)fm(x)fTf(x)=1vgφgm,n+1-φgm,nΔtn+1-Σg′=1G(χpgβpνΣfg′m,n+1+Σg′gm,n+1)φg′m,n+1-χdgΣlλlClm,n+1+Σtgm,n+1φgm,n+1-Lgm,n+1ΔXjΔtAρj(hj-hjn)-mj-1(hj-1-hj)-Hn+1PrΔX(Tw-Tb)AΔXjΔt(ρj-ρjn)+mj-mj-1(ρfcpV)in+1Tin+1-TinΔt+(Ki+Ki-1)Ti-Ki-1Ti-1-KiTi+1-qi′′′=0]]>式中，fφ(x),fh(x),fm(x)分別為殘差形式的中子方程，傳熱方程，能量方程和質(zhì)量方程；F(x)便是聯(lián)立的物理熱工方程組，x為解向量，分別為中子通量，冷卻劑比焓，冷卻劑質(zhì)量流量，燃料溫度；步驟3：建立了物理熱工耦合的非線性方程組F(x)，在用牛頓法求解F(x)時，需在每一個牛頓迭代步構(gòu)造F(x)的雅閣比矩陣，如下所示，JF(x)=∂fφ∂φ∂fφ∂h∂fφ∂m∂fφ∂T∂fh∂φ∂fh∂h∂fh∂m∂fh∂T∂fm∂φ∂fm∂h∂fm∂m∂fm∂T∂fT∂φ∂fT∂h∂fT∂m∂fT∂T]]>在雅閣比矩陣元素的計算中，需要求解中子殘差方程對比焓的偏導(dǎo)數(shù)這實際上是要計算截面關(guān)于比焓的偏導(dǎo)數(shù)由于截面是關(guān)于慢化劑密度和燃料溫度的隱式函數(shù)Σx＝g(ρ,Tf)，因此使用求導(dǎo)的鏈式法則計算得到，∂Σx∂h=∂Σx∂ρ∂ρ∂h=∂g∂ρ∂ρ∂h]]>對于雅閣比矩陣中剩下的元素的計算，均能夠通過各個殘差方程得到解析表達式；通過以上運算便得到了雅閣比矩陣各個元素的數(shù)值，為了獲得每個牛頓迭代步的更新量δx，需要求解以下代數(shù)方程組，JFδx＝-F(x)JF是一個大型的稀疏矩陣，可以使用許多成熟的代數(shù)方程組迭代方法，如GMRES算法；以k作為牛頓迭代步記號，當||F(xk)||≤εNewton時認為牛頓迭代收斂，||F(xk)||為F(xk)的范數(shù)，εNewton為牛頓迭代收斂判據(jù)；步驟4：物理熱工耦合的非線性方程組F(x)求解完畢后，需要根據(jù)最新解得的節(jié)塊中子通量更新節(jié)塊耦合修正因子Dk,nod，因此是在牛頓迭代法外層再涵蓋一層用于更新Dk,nod的非線性迭代步；節(jié)塊耦合修正因子Dk,nod的計算與半解析節(jié)塊法相同。與現(xiàn)有技術(shù)相比，本發(fā)明具有如下幾個突出特點：1.統(tǒng)一求解聯(lián)立的物理熱工方程組，同時解得中子通量，慢化劑比焓，燃料溫度。2.將理論上具備二階收斂速度的牛頓法用于物理熱工耦合計算分析。3.顯示構(gòu)造求物理熱工非線性方程組的雅各比矩陣，避免了近似計算雅閣比矩陣對迭代收斂性的影響。4.三重迭代格式，最外層為非線性迭代，計算節(jié)塊耦合修正因子，中間層為牛頓迭代，計算物理熱工耦合非線性方程組，內(nèi)層為雅閣比矩陣求解的線性方程組迭代，非線性迭代可以保證中子方程的解收斂于高階節(jié)塊法，牛頓迭代則保證了耦合計算的精度和效率。附圖說明圖1三重迭代格式結(jié)構(gòu)示意圖。具體實施方式首先確定所需聯(lián)立的方程形式，并建立封閉的物理熱工聯(lián)立方程組。進行空間時間能量離散，緩發(fā)中子先驅(qū)核濃度進行時間上的線性近似后，節(jié)塊內(nèi)的2群時空中子擴散方程為1vgφgm,n+1-φgm,nΔtn+1=Σg′=1G(χpgβpνΣfg′m,n+1+Σg′gm,n+1)φg′m,n+1+χdgΣlλlClm,n+1-Σtgm,n+1φgm,n+1-Lgm,n+1,,g=1,2---(1)]]>Cln+1=e-λlΔtnCln+Fln0Σg′=1G(νΣfg′φg′)n+Fln1Σg′=1G(νΣfg′φg′)n+1---(2)]]>Lgm=Σu=x,y,z1Δum(Jgu+m-Jgu-m)]]>Jgu+m=-Dgm,FDM(fgu-m+1φgm+1-fgu+mφgm)-Dgm,NOD(fgu-m+1φgm+1+fgu+mφgm)]]>式中：G——能群總數(shù)；Δt——時間步長；——第g群n時刻節(jié)塊m內(nèi)的平均中子通量；——第g群n+1時刻節(jié)塊m內(nèi)的平均中子通量；vg——第g群中子速度；χpg——瞬發(fā)中子裂變譜；βp——瞬發(fā)中子份額；χdg——緩發(fā)中子裂變譜；νΣf——中子產(chǎn)生截面；Σg'g——g’群到g群宏觀散射截面；Σt——宏觀總截面；χpg——瞬發(fā)中子裂變譜；λl——瞬發(fā)中子份額；Cl——第l群緩發(fā)中子先驅(qū)核濃度；——由先驅(qū)核濃度線性近似推導(dǎo)得到的常數(shù)項系數(shù)；——第g群節(jié)塊m內(nèi)的中子泄漏項；——第g群節(jié)塊m內(nèi)u方向左右表面的中子靜流；Δu——u方向的節(jié)塊寬度；——節(jié)塊等效擴散系數(shù)；——非線性迭代節(jié)塊耦合修正因子；——m節(jié)塊在u方向左側(cè)面的不連續(xù)因子；由緩發(fā)先驅(qū)核濃度方程(2)式可以看出，新一時刻的緩發(fā)先驅(qū)核濃度可在中子通量求出后代入得到，因此全耦合方法中無需將緩發(fā)先驅(qū)核濃度方程放入最終的大方程組中聯(lián)立求解，所需求解的未知量為節(jié)塊中子平均通量此處需特別指出，對于非線性迭代節(jié)塊耦合修正因子僅需知道它是通過對每個節(jié)塊求解一個局部兩節(jié)塊問題得到的，其相關(guān)理論十分成熟，與非線性迭代半解析節(jié)塊法完全一致。基于單通道模型的離散后的熱工方程組如下AΔXjΔt(ρj-ρjn)+mj-mj-1=0---(3)]]>ΔXjΔt(mj-mjn)+mjUj′-mj-1Uj-1′=-A(Pj-Pj-1)-gAΔXjρj-12(ΔXfφ2Dhρl)j|mj|mjA---(4)]]>ΔXjΔtAρj(hj-hjn)-mj-1(hj-1-hj)=Hn+1PrΔX(Tw-Tb)---(5)]]>(ρfcpV)in+1Tin+1-TinΔt+(Ki+Ki-1)Ti-Ki-1Ti-1-KiTi+1=qi′′′---(6)]]>(3)～(6)式依次為質(zhì)量守恒方程，動量守恒方程，能量守恒方程，傳熱方程，式中j——流體場網(wǎng)格標識；i——傳熱場網(wǎng)格標識A——通道流動面積；ρj——第j網(wǎng)格冷卻劑密度；mj——第j網(wǎng)格質(zhì)量流量；hj——第j網(wǎng)格冷卻劑比焓；U'——冷卻劑流速；P——壓力；g——重力加速度；f——摩擦系數(shù)；φ2——兩相摩擦倍增因子；Dh——水力學(xué)直徑；ρl——冷卻劑液相密度；H——傳熱系數(shù)；Pr——釋熱周長；T——壁面溫度；Tb——冷卻劑溫度；ρf——燃料密度；cp——燃料比熱；V——傳熱網(wǎng)格體積；K——等效導(dǎo)熱率；q″′——體積釋熱率；熱工方程中未知量為h,m,P,T，若忽略壓力對流體物性的影響，則由動量守恒方程(4)式可以看出，當求得h,m后，便可由(4)式求得壓力P，從而全耦合方法最終聯(lián)立的大方程組只包含(3)，(5)，(6)式，所需求解未知量為h,m,T。經(jīng)過以上討論，將(2)，(3)，(5)，(6)式聯(lián)立，并將每個方程的所有項移到等式同一端得到各個方程的殘差形式，如下F(x,DNOD)=fφ(φ,h,m,T,DNOD)fm(φ,h,m,T)fh(φ,h,m,T)fT(φ,h,m,T)=0,x=[φ,h,m,T]T]]>如圖1所示，全耦合方法中，最外層非線性迭代更新DNOD，得到新的DNOD值后，內(nèi)層的牛頓迭代實際上求解了如下問題F(x)=fφ(φ,h,m,T)fm(φ,h,m,T)fh(φ,h,m,T)fT(φ,h,m,T)=0]]>其雅閣比矩陣為JF(x)=∂fφ∂φ∂fφ∂h∂fφ∂m∂fφ∂T∂fm∂φ∂fm∂h∂fm∂m∂fm∂T∂fh∂φ∂fh∂h∂fh∂m∂fh∂T∂fT∂φ∂fT∂h∂fT∂m∂fT∂T]]>根據(jù)方程性質(zhì)，雅閣比矩陣實際上是稀疏的，這是由于中子平衡方程與質(zhì)量流量m并無直接關(guān)系，因此雅閣比矩陣中元素值為零，同理，值為零的元素還有因此雅閣比矩陣中實際上需要計算的元素為JF(x)=∂fφ∂φ∂fφ∂h0∂fφ∂T0∂fm∂h∂fm∂m0∂fh∂φ∂fh∂h∂fh∂m∂fh∂T∂fT∂φ∂fT∂h0∂fT∂T]]>對于某一個牛頓迭代步，其計算邏輯為JF(xk)δxk＝-F(xk)xk+1＝xk+δxk求解上式是一個線性方程組求解問題，當問題規(guī)模不大時可采用直接法，如LU分解，但一般說來對于堆芯量級的問題，雅閣比矩陣規(guī)模會特別巨大，可采用GMRES，雙共軛梯度法等一系列數(shù)學(xué)上成熟的大型線性方程組數(shù)值方法求解。當滿足如下收斂判據(jù)時，牛頓迭代結(jié)束||F(x)||＜εNewton當滿足非線性迭代層收斂判據(jù)時，非線性迭代結(jié)束，至此全耦合法完成一次完整計算。當前第1頁1 2 3