核反應堆熱工分析課程設計

PAGE南京工程學院課程設計報告設計題目某壓水核反應堆熱工水力設計課程名稱核反應堆課程設計系部能源與動力工程學院專業(yè)熱能與動力工程（核電）班級核電集控081學號1111111姓名qq起止日期2011.6.27～2011.7.8指導教師22222日期2011年7月6日目錄設計任務課程設計目的計算過程及分析程序程序設計框圖代碼說明書課程設計總結參考資料代碼

一、設計任務某壓水反應堆的冷卻劑及慢化劑都是水，用二氧化鈾作燃料，用鋯-4作包殼材料。燃料組件無盒壁，燃料元件為棒狀，正方形排列。已知下列參數(shù)：系統(tǒng)壓力16MPa堆芯輸出功率1840MW冷卻劑總流量32600t/h反應堆進口溫度288℃堆芯高度3.8m燃料組件數(shù)121燃料組件形式17×17每個組件燃料棒數(shù)265燃料包殼直徑9.6mm燃料包殼內(nèi)徑8.5mm燃料包殼厚度0.55mm燃料芯塊直徑8.35mm燃料棒間距（柵距）12.5mm芯塊密度95%理論密度旁流系數(shù)5%燃料元件發(fā)熱占總發(fā)熱的份額97.4%徑向核熱管因子1.35軸向核熱管因子1.53局部峰核熱管因子1.12交混因子0.95熱流量工程熱點因子1.08焓升工程熱管因子1.08堆芯入口局部阻力系數(shù)0.75堆芯出口局部阻力系數(shù)1.0堆芯定位隔架局部阻力系數(shù)1.05燃料元件中心最高溫度不超過 2200℃若將堆芯自上而下劃分為6個控制體，則其軸向歸一化功率分布如下表：堆芯軸向歸一化功率分布（軸向等分6個控制體）自上而下控制體號123456歸一化功率分布0.481.021.501.560.960.48通過計算，得出堆芯出口溫度；燃料棒表面平均熱流及最大熱流密度，平均線功率，最大線功率；熱管的焓，包殼表面溫度，芯塊中心溫度隨軸向大的分布；包殼表面最高溫度，芯塊中心最高溫度；DNBR在軸向的變化；計算堆芯壓降；二、課程設計的目的1、深入理解壓水堆熱工設計準則；2、深入理解單通道模型的基本概念，基本原理。包括了解平均通道（平均管），熱通道（熱管）,熱點等在反應堆熱工水力設計中的作用；3、掌握堆芯焓場的計算并求出體現(xiàn)反應堆安全性的主要參數(shù)；燒毀比DNBR,最小燒毀比MDNBR,燃料元件中心最高溫度t0及其最高溫度t0,max,包殼表面溫度tcs4、求出體現(xiàn)反應堆先進性的主要參數(shù)；堆芯流量功率比，堆芯功率密度，燃料元件平均密度(熱通量)，最大熱流密度，冷卻劑平均流速，冷卻劑出口溫度等；5、通過本課程設計，掌握壓水堆熱工校核的具體工具；6、掌握壓降的計算；7、掌握單項及沸騰時的傳熱計算。三、計算過程及過程目前，壓水核反應堆的穩(wěn)態(tài)熱工設計準則有：（1）燃料元件芯塊內(nèi)最高溫度應低于其相應燃耗下的熔化溫度。目前，壓水堆大多采用UO2作為燃料。二氧化鈾的熔點約為2805±15℃，經(jīng)輻照后，其熔點會有所降低。燃耗每增加104兆瓦·日/噸鈾，其熔點下降32℃。在通常所達到的燃耗深度下，熔點將降至2650℃左右。在穩(wěn)態(tài)熱工設計中，一般將燃料元件中心最高溫度限制在2200～2450℃之間。（2）燃料元件外表面不允許發(fā)生沸騰臨界。通常用臨界熱流密度比DNBR來定量地表示這個限制條件。DNBR是根據(jù)堆內(nèi)某處燃料元件周圍的冷卻劑狀態(tài)使用專門的計算公式而得到的臨界熱流密度與該處燃料元件表面的實際熱流密度的比值。DNBR隨堆芯通道的長度是變化的，在整個堆芯內(nèi)，DNBR的最小值稱為最小DNBR，用MDNBR或DNBRmin表示。為了確保燃料元件不燒毀，當計算的最大熱功率下，MDNBR不應低于某一規(guī)定值。如果計算熱流密度的公式?jīng)]有誤差，則當MDNBR=1時，表示燃料元件表面要發(fā)生沸騰臨界。若該公式存在誤差，則MDNBR就要大于1。例如，W-3公式的誤差為23%，所以當使用W-3公式計算DNBR時，就要求MDNBR≥1.3。（3）必須保證正常運行工況下燃料元件和堆內(nèi)構件能得到充分冷卻；在事故工況下能提供足夠的冷卻劑以排出堆芯余熱。（4）在穩(wěn)態(tài)額定工況下，要求在計算的最大熱功率下，不發(fā)生流動不穩(wěn)定性。對于壓水堆，只要在堆芯最熱通道出口附近冷卻劑中的含氣量不大于某一數(shù)值，就不會發(fā)生流動不穩(wěn)定性。在反應堆內(nèi)，即使燃料元件的形狀、尺寸、密度和裂變物質(zhì)的濃縮度都相同，堆芯內(nèi)的中子通量分布也是不均勻的，再加上堆芯內(nèi)存在控制棒，水隙、空泡及反射層的影響，中子通量的分布更是不均勻的。從而，堆芯內(nèi)的熱功率分布也是不均勻的。而燃料元件在加工、安裝及運行中的各類工程因素也能造成實際值與設計值之間產(chǎn)生偏差。為了表示有關的熱工參數(shù)的最大值偏離平均值（或名義值）的程度，引入了熱管因子的概念。分兩類：核熱管因子和工程熱管因子。本計算根據(jù)《核反應堆熱工分析》課程設計指導書中的計算提示，采用簡單的matlab語言編程計算。將堆芯沿軸向劃分為六個等分控制體進行計算以下是計算結果：1.流體堆芯出口溫度tf,out=324.0686℃2.堆芯內(nèi)燃料棒的總傳熱面積Ft=3.6748e+03m3.燃料棒表面平均熱流密度q=4.8768e+05w/m24.燃料棒表面最大熱流密度qmax=1.0879e+06w/m5.燃料棒平均線功率ql=1.4708e+046.燃料棒最大線功率ql,max=3.281e+047.堆芯內(nèi)總流通面積Af=2.7715m8.平均管流速V=4.3436m/s；9.單元通道內(nèi)流量Wu=0.10.單元通道面積Ab=8.3868e-511.第一控制體出口流體溫度tf,h（L1）=292.3695℃12.第一控制體出口處的包殼外壁溫tcs,h（L1）=303.2159℃13.第一控制體出口處的包殼內(nèi)壁溫tci,h(L1)=307.399114.第一控制體出口處的UO2芯塊外表面溫度tu,h（L1）=15.第一控制體出口處的UO2芯塊中心溫度to,h(L1)=16.熱管中的qDNB,h(L1)=5.8949e+6w/m17.DNBR(L1)=17.271918.第二控制體出口流體溫度tf,h（L2）=300.7738℃19.第二控制體出口處的包殼外壁溫tcs,h（L2）=321.9306℃20.第二控制體出口處的包殼內(nèi)壁溫tci,h(L2)=330.138421.第二控制體出口處的UO2芯塊外表面溫度tu,h（L2）=22.第二控制體出口處的UO2芯塊中心溫度to,h(L2)=23.熱管中的qDNB,h(L2)=5.3775e+6w/m24.DNBR(L2)=7.877925.第三控制體出口流體溫度tf,h（L3）=314.0155℃26.第三控制體出口處的包殼外壁溫tcs,h（L3）=347.2633℃27.第三控制體出口處的包殼內(nèi)壁溫tci,h(L3)=360.267228.第三控制體出口處的UO2芯塊外表面溫度tu,h（L3）=29.第三控制體出口處的UO2芯塊中心溫度to,h(L3)=30.熱管中的qDNB,h(L3)=4.5508e+6w/m31.DNBR(L3)=4.102632.第四控制體出口流體溫度tf,h（L4）=326.2418℃33.第四控制體出口處的包殼外壁溫tcs,h（L4）=349.2519℃34.第四控制體出口處的包殼內(nèi)壁溫tci,h(L4)=361.731235.第四控制體出口處的UO2芯塊外表面溫度tu,h（L4）=36.第四控制體出口處的UO2芯塊中心溫度to,h(L4)=37.熱管中的qDNB,h(L4)=3.7716e+6w/m38.DNBR(L4)=3.536139.第五控制體出口流體溫度tf,h（40.第五控制體出口處的包殼外壁溫tcs,h（L5）=349.075℃41.第五控制體出口處的包殼內(nèi)壁溫tci,h(L5)=357.310642.第五控制體出口處的UO2芯塊外表面溫度tu,h（L5）=43.第五控制體出口處的UO2芯塊中心溫度to,h(L5)=44.熱管中的qDNB,h(L5)=3.2719e+6w/m45.DNBR(L5)=4.511346.第六控制體出口流體溫度tf,h（47.第六控制體出口處的包殼外壁溫tcs,h（L6）=346.1839℃48.第六控制體出口處的包殼內(nèi)壁溫tci,h(L6)=350.070249.第六控制體出口處的UO2芯塊外表面溫度tu,h（L6）=50.第六控制體出口處的UO2芯塊中心溫度to,h(L6)=51.熱管中的qDNB,h(L6)=3.058e+6w/m52.DNBR6=8.959853.單相流體的摩擦壓降△Pf=2.6511e+454.單相流體加速壓降△Pa=055.單相流體提升壓降△Pel=2.6612e+456.堆芯出口局部壓降△Pout=7.6260e+357.堆芯進口局部壓降△Pin=4.8090e+358.定位格架出口壓降△Pgr=7.37e+359.總的壓降△P=7.2928e+4Pa計算結果誤差分析：由于采用的是W-3公式，且該設計中的給出參數(shù)與該公式的適用范圍有些偏差，且在計算物性時粗糙地采用了線性插值的方法，更是帶來了較大誤差。但是其算出的結果還是能客觀反映出熱管中各量的變化趨勢的。臨界熱流與燒毀比的匯總表項目臨界熱流×10^6w/m^2DNBRL5.894917.27192L5.37757.87793L4.55084.10264L3.77163.53615L3.27194.51136L3.0588.9598臨界熱流與燒毀比隨著高度增加的變化趨勢圖從圖上可以觀察出，燒毀比在中間段時有最小值，即堆芯中間段是最危險的地方。而臨界熱流隨著高度的增加逐漸減小。各溫度的匯總表項目控制體出口溫度包殼外表面溫度包殼內(nèi)表面溫度芯塊表面溫度芯塊中心溫度L292.3695303.2159307.3991375.892560.6872L300.7738321.9306330.1384467.1242923.62423L314.0155347.2633360.2672582.7261513.94L326.2418349.2519361.7312575.61671456.65L333.2833349.075357.3106502.8581014.36L336.468346.1839350.0702418.5631613各溫度隨高度變化趨勢圖從這幅反應各溫度隨軸向高度變化情況的表中可以知道：冷卻劑的溫度隨著高度的增加不斷上升；包殼外表面、內(nèi)表面、芯塊表面、芯塊中心的溫度都有一個最值，且最大值所在的高度逐漸的往堆芯半高處移動。四、程序開始1、程序設計框圖開始讀輸入?yún)?shù)讀輸入?yún)?shù)計算有關堆參數(shù)計算有關堆參數(shù)估算控制體出口溫度估算控制體出口溫度tf計算控制體出口溫度計算控制體出口溫度θ|1-|1-θtf 重估tf 重估tf計算該處含汽量計算該處含汽量 是計算包殼外表面溫度計算包殼外表面溫度根據(jù)根據(jù)W-3算臨界熱流包殼內(nèi)表面溫度包殼內(nèi)表面溫度芯塊表面溫度計算燒毀比芯塊表面溫度計算燒毀比芯塊中心溫度芯塊中心溫度打印輸出值打印輸出值停機停機2、代碼說明書本代碼主要由五個小部分組成。堆芯出口溫度計算、堆熱流量計算、堆平均參數(shù)計算、第一至第六控制體各量計算、熱管的壓降計算。（1）堆芯出口溫度計算：此段根據(jù)任務書給出的基本參數(shù)和熱量與流量之間關系，運用迭代的算法，求出堆芯的出口溫度。（2）堆熱流量計算：先根據(jù)堆芯的輸出功率和釋熱率以及總的傳熱面積，求出燃料元件表面平均熱流量，再根據(jù)熱管因子求出最大熱流量。再求出平均線功率和最大線功率。（3）堆平均參數(shù)計算：根據(jù)基本的尺寸，求出堆體的流通截面積和一個柵元的流通截面積。然后再求出流經(jīng)柵元的流量。依據(jù)上面的溫度結果，查出熱物性參數(shù)，再求出冷卻劑的流速。（4）第一至第六控制體的各量計算：因為六個控制體的計算過程類似，這里只說明第一個控制體的計算過程。在現(xiàn)有的參數(shù)下，根據(jù)熱流量與流量的關系和迭代算法，求出該控制體的出口溫度。通過流通截面積與濕周的關系求出柵元的當量直徑。再根據(jù)上面的溫度，查出對應的熱物性參數(shù)由雷諾數(shù)與努爾數(shù)的關系，解出控制體出口處的對流換熱系數(shù)。因為不知該處的流體狀態(tài)，分別用單相強迫對流放熱公式和詹斯-洛特斯傳熱方程算出各自的膜溫壓，取較小的值加上出口處的流體溫度即是包殼的外表面溫度。由包殼的外表面的溫度再根據(jù)圓管的傳熱方程運用迭代算法解出包殼內(nèi)表面的溫度。芯塊與包殼內(nèi)表面之間的導熱問題，根據(jù)間隙導熱模型，即可解出芯塊表面的溫度，根據(jù)內(nèi)熱源的導熱模型，依據(jù)積分熱導率與溫度的對應關系和插值方法，解出芯塊中心的溫度。接下來依據(jù)冷卻劑的溫度，得出的控制體出口處的含汽量。進而依據(jù)W-3公式求出該出的臨界熱流量qDNB（5）熱管的壓降計算：熱管的壓降包括摩擦壓降、提升壓降、進出口局部壓降、定位擱架出口壓降。摩擦壓降可由計算單相流的達西（Darcy）公式算得。提升壓降可由根據(jù)位置的變化算得，其中參數(shù)都取平均值。其余的壓降根據(jù)形阻壓降的基本公式再乘以相應的系數(shù)求得。最后各項相加得出熱管的總壓降。五、課程設計小結為期兩周的課程設計結束了，這兩周里可謂受益匪淺。剛開始的時候，對于這個新鮮的設計感到一籌莫展，真不知道從哪里下手。后來在老師的指導下，先靜下心來把書看了一遍。在看書的同時，我也結合課程設計指導書，將一些沒有頭緒的問題慢慢理順了。漸漸地，在和老師、同學的交流中我把問題弄明白了。知道了什么是熱管因子，怎么理解歸一化系數(shù)，學會了如何從冷卻劑的溫度逐步算出芯塊中心溫度。在我們的這個小組中，經(jīng)過不斷的討論還提高了我們的團隊解決問題的能力。最后還要感謝老師的耐心指導，才能使我順利的完成課程設計。六、參考資料1．《核反應堆熱工分析》，于平安著，上海交通大學出版社2．《反應堆熱工流體力學》，趙兆頤,朱瑞安，清華大學出版社。3．《反應堆熱工水力學》，俞冀陽，賈寶山，清華大學出版社。4．《壓水反應堆熱工分析》，(美)湯糧孫，J.韋斯曼，袁乃駒譯，原子能出版社。5．《MATLAB語言與控制系統(tǒng)仿真》，孫亮主編，北京工業(yè)大學出版讓，2001。七、程序代碼%流體堆芯出口溫度計算tfin=288;Fa=0.974;Nt=1840e+6;Wt=9055.56;b=0.05;tfout=323;e0=0.01;whilee0>0.001t0_=0.5*(tfout+tfin);Cp_=1000*(0.0265*(t0_-300)+5.63);xi=tfin+Fa*Nt/(Wt*(1-b)*Cp_);e0=(tfout-xi)/tfout;tfout=xi%堆芯出口處溫度end%熱流密度計算m=121;n=265;dcs=9.6e-3;L=3.8;q_=Fa*Nt/(m*n*pi*dcs*L)%燃料元件表面平均熱流量FRN=1.35;FZN=1.53;FqN=FRN*FZN;FqE=1.08;FDHE=1.08;FDHmE=0.95;qmax=q_*FqN*FqE%最大熱流量ql_=q_*pi*dcs%平均線功率qlmax=ql_*FqN*FqE%最大線功率%平均管情況B=17;S=12.5e-3;dx=0.8e-3;Af=m*n*(S^2-pi/4*dcs^2)+m*4*B*S*dx;%總的流通截面積tf_=0.5*(tfout+tfin)%熱管平均溫度vf_=4.2e-6*(tf_-300)+0.001374;pf_=1/vf_;%平均密度v=Wt*(1-b)/(Af*pf_);%平均流速Ab=S^2-pi/4*dcs^2;%單元流通截面積Wu=Wt*(1-b)*Ab/Af;%單元截面流量%第一控制體溫度計算e11=0.01;tf1=300;L1=3.8/6;fai1=0.48;whilee11>0.001t11_=0.5*(tf1+tfin);Cp1_=1000*(0.0269*(t11_-280)+5.068);x1i=tfin+q_*FRN*FDHE*FDHmE*pi*dcs*L1*fai1/(Wu*Cp1_);e11=(x1i-tf1)/tf1;tf1=x1i%求出該控制體出口處的溫度endDe=4*(S^2-pi/4*dcs^2)/(pi*dcs);%單元通道當量直徑u1=944e-7;Pr1=0.85;k1=575.5e-3;%查得該溫度下的熱物性Re1=Wu*De/(Ab*u1);h1=0.023*Re1^0.8*Pr1^0.4*k1/De;%該處的對流換熱系數(shù)dtf11=q_*FRN*fai1*FqE/h1;%單相強迫對流放熱公式算得的溫壓ts=347.328;P=16;dtf12=25*(q_*FRN*fai1*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf1;%采用詹斯-洛特斯傳熱方程算得的過冷沸騰膜溫壓ifdtf11<dtf12%膜溫壓取兩個中較小值，算得包殼外表面溫度tcs1=tf1+dtf11elsetcs1=tf1+dtf12enddci=8.5e-3;tci1=349;e12=0.01;whilee12>0.001t12_=0.5*(tci1+tcs1);kc1=0.0547*(1.8*t12_+32)+13.8;yi=tcs1+ql_*FRN*fai1*FqE/(2*pi*kc1)*log(dcs/dci);e12=(yi-tci1)/yi;tci1=yi%采用迭代算法求得包殼內(nèi)表面溫度endhg=5678;du=8.35e-3;tu1=tci1+ql_*FRN*FqE*fai1*2/(pi*(dci+du)*hg)%燃料芯塊表面溫度d1_ku=ql_*FRN*FqE*fai1/(4*pi*100);tu1_ku=(26.42-21.32)/(400-300)*(tu1-300)+21.32;to1_ku=tu1_ku+d1_ku;to1=(600-500)/(34.97-30.93)*(to1_ku-30.93)+500%根據(jù)積分熱導率圖表查得芯塊中心溫度p=16e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3600;h1=1296.4746e+3;x1=(h1-hfs)/(hgs-hfs);%該點含汽量qDNB1=3.154e6*((2.022-6.238e-8*p)+...%根據(jù)W-3公式計算出臨界熱流量(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x1))*...((0.1484-1.596*x1+0.1729*x1*abs(x1))*0.2049*G/10^6+1.037)*...(1.157-0.869*x1)*...(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))DNBR1=qDNB1/(q_*FRN*FqE*fai1)%計算燒毀比%第二控制體溫度計算fai2=0.96;L2=3.8/6;e21=0.01;tf2=310;whilee21>0.001t21_=0.5*(tf1+tf2);Cp2_=1000*(0.0265*(t21_-300)+5.63);x2i=tf1+q_*FRN*FDHE*FDHmE*pi*dcs*L2*fai2/(Wu*Cp2_);e21=(x2i-tf2)/tf2;tf2=x2i%求出該控制體出口處的溫度endDe=4*(S^2-pi/4*dcs^2)/(pi*dcs);u2=919e-7;Pr2=0.91;k2=562e-3;%查得該溫度下的熱物性Re2=Wu*De/(Ab*u2);h2=0.023*Re2^0.8*Pr2^0.4*k2/De;%該處的對流換熱系數(shù)dtf21=q_*FRN*fai2*FqE/h2;%單相強迫對流放熱公式算得的溫壓ts=347.328;P=16;dtf22=25*(q_*FRN*fai2*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf2;%采用詹斯-洛特斯傳熱方程算得的過冷沸騰膜溫壓ifdtf21<dtf22%膜溫壓取兩個中較小值，算得包殼外表面溫度tcs2=tf2+dtf21elsetcs2=tf2+dtf22enddci=8.5e-3;tci2=349;e22=0.01;whilee22>0.001t22_=0.5*(tci2+tcs2);kc2=0.0547*(1.8*t22_+32)+13.8;zi=tcs2+ql_*FRN*fai2*FqE/(2*pi*kc2)*log(dcs/dci);e22=(zi-tci2)/zi;tci2=zi%采用迭代算法求得包殼內(nèi)表面溫度endhg=5678;du=8.35e-3;tu2=tci2+ql_*FRN*FqE*fai2*2/(pi*(dci+du)*hg)%燃料芯塊表面溫度d2_ku=ql_*FRN*FqE*fai2/(4*pi*100);tu2_ku=(30.93-26.42)/(500-400)*(tu2-400)+26.42;to2_ku=tu2_ku+d2_ku;to2=(1000-900)/(48.06-45.14)*(to2_ku-45.14)+900%根據(jù)積分熱導率圖表查得芯塊中心溫度p=16e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3600;h2=1341.5988e+3;x2=(h2-hfs)/(hgs-hfs);%該點含汽量qDNB2=3.154e6*((2.022-6.238e-8*p)+...%根據(jù)W-3公式計算出臨界熱流量(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x2))*...((0.1484-1.596*x2+0.1729*x2*abs(x2))*0.2049*G/10^6+1.037)*...(1.157-0.869*x2)*...(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))DNBR2=qDNB2/(q_*FRN*FqE*fai2)%計算燒毀比%第三控制體溫度計算fai3=1.56;L3=3.8/6;e31=0.01;tf3=315;whilee31>0.001t31_=0.5*(tf3+tf2);Cp3_=1000*(0.0265*(t31_-300)+5.63);x3i=tf2+q_*FRN*FDHE*FDHmE*pi*dcs*L3*fai3/(Wu*Cp3_);e31=(x3i-tf3)/tf3;tf3=x3i%求出該控制體出口處的溫度endDe=4*(S^2-pi/4*dcs^2)/(pi*dcs);u3=869e-7;Pr3=1.01;k3=533e-3;%查得該溫度下的熱物性Re3=Wu*De/(Ab*u3);h3=0.023*Re3^0.8*Pr3^0.4*k3/De;%該處的對流換熱系數(shù)dtf31=q_*FRN*fai3*FqE/h3;%單相強迫對流放熱公式算得的溫壓ts=347.328;P=16;dtf32=25*(q_*FRN*fai3*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf3;%采用詹斯-洛特斯傳熱方程算得的過冷沸騰膜溫壓ifdtf31<dtf32%膜溫壓取兩個中較小值，算得包殼外表面溫度tcs3=tf3+dtf31elsetcs3=tf3+dtf32enddci=8.5e-3;tci3=349;e32=0.01;whilee32>0.001t32_=0.5*(tci3+tcs3);kc3=0.0547*(1.8*t32_+32)+13.8;ai=tcs3+ql_*FRN*fai3*FqE/(2*pi*kc3)*log(dcs/dci);e32=(ai-tci3)/ai;tci3=ai%采用迭代算法求得包殼內(nèi)表面溫度endhg=5678;du=8.35e-3;tu3=tci3+ql_*FRN*FqE*fai3*2/(pi*(dci+du)*hg)%燃料芯塊表面溫度d3_ku=ql_*FRN*FqE*fai3/(4*pi*100);tu3_ku=(34.97-30.93)/(600-500)*(tu3-500)+30.93;to3_ku=tu3_ku+d3_ku;to3=(1560-1405)/(61.95-58.4)*(to3_ku-58.4)+1405%根據(jù)積分熱導率圖表查得芯塊中心溫度p=16e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3600;h3=1416.5e+3;x3=(h3-hfs)/(hgs-hfs);%該點含汽量qDNB3=3.154e6*((2.022-6.238e-8*p)+...%根據(jù)W-3公式計算出臨界熱流量(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x3))*...((0.1484-1.596*x3+0.1729*x3*abs(x3))*0.2049*G/10^6+1.037)*...(1.157-0.869*x3)*...(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))DNBR3=qDNB3/(q_*FRN*FqE*fai3)%計算燒毀比%第4控制體溫度計算 fai4=1.5;L4=3.8/6;e41=0.01;tf4=320;whilee41>0.001t41_=0.5*(tf3+tf4);Cp4_=1000*(0.0265*(t41_-300)+5.63);x4i=tf3+q_*FRN*FDHE*FDHmE*pi*dcs*L4*fai4/(Wu*Cp4_);e41=(x4i-tf4)/tf4;tf4=x4i%求出該控制體出口處的溫度endDe=4*(S^2-pi/4*dcs^2)/(pi*dcs);u4=821e-7;Pr4=1.08;k4=505e-3;%查得該溫度下的熱物性Re4=Wu*De/(Ab*u4);h4=0.023*Re4^0.8*Pr4^0.4*k4/De;%該處的對流換熱系數(shù)dtf41=q_*FRN*fai4*FqE/h4;%單相強迫對流放熱公式算得的溫壓ts=347.328;P=16;dtf42=25*(q_*FRN*fai4*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf4;%采用詹斯-洛特斯傳熱方程算得的過冷沸騰膜溫壓ifdtf41<dtf42%膜溫壓取兩個中較小值，算得包殼外表面溫度tcs4=tf4+dtf41elsetcs4=tf4+dtf42enddci=8.5e-3;tci4=349;e42=0.01;whilee42>0.001t42_=0.5*(tci4+tcs4);kc4=0.0547*(1.8*t42_+32)+13.8;bi=tcs4+ql_*FRN*fai4*FqE/(2*pi*kc4)*log(dcs/dci);e42=(bi-tci4)/bi;tci4=bi%采用迭代算法求得包殼內(nèi)表面溫度endhg=5678;du=8.35e-3;tu4=tci4+ql_*FRN*FqE*fai4*2/(pi*(dci+du)*hg)%燃料芯塊表面溫度d4_ku=ql_*FRN*FqE*fai4/(4*pi*100);tu4_ku=(34.97-30.93)/(600-500)*(tu4-500)+30.93;to4_ku=tu4_ku+d4_ku;to4=(1560-1405)/(61.95-58.4)*(to4_ku-58.4)+1405%根據(jù)積分熱導率圖表查得芯塊中心溫度p=16e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3600;h4=1490.58e+3;x4=(h4-hfs)/(hgs-hfs);%該點含汽量qDNB4=3.154e6*((2.022-6.238e-8*p)+...%根據(jù)W-3公式計算出臨界熱流量(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x4))*...((0.1484-1.596*x4+0.1729*x4*abs(x4))*0.2049*G/10^6+1.037)*...(1.157-0.869*x4)*...(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))DNBR4=qDNB4/(q_*FRN*FqE*fai4)%計算燒毀比%第5控制體溫度計算 fai5=1.02;L5=3.8/6;e51=0.01;tf5=330;whilee51>0.001t51_=0.5*(tf5+tf4);Cp5_=1000*(0.0877*(t51_-320)+6.16);x5i=tf4+q_*FRN*FDHE*FDHmE*pi*dcs*L5*fai5/(Wu*Cp5_);e51=(x5i-tf5)/tf5;tf5=x5i%求出該控制體出口處的溫度endDe=4*(S^2-pi/4*dcs^2)/(pi*dcs);u5=788e-7;Pr5=1.17;k5=485e-3;%查得該溫度下的熱物性Re5=Wu*De/(Ab*u5);h5=0.023*Re5^0.8*Pr5^0.4*k5/De;%該處的對流換熱系數(shù)dtf51=q_*FRN*fai5*FqE/h5;%單相強迫對流放熱公式算得的溫壓ts=347.328;P=16;dtf52=25*(q_*FRN*fai5*FqE/10^6)^0.25*exp(-P/6.2)+ts-tf5;%采用詹斯-洛特斯傳熱方程算得的過冷沸騰膜溫壓ifdtf51<dtf52%膜溫壓取兩個中較小值，算得包殼外表面溫度tcs5=tf5+dtf51elsetcs5=tf5+dtf52enddci=8.5e-3;tci5=380;e52=0.01;whilee52>0.001t52_=0.5*(tci5+tcs5);kc5=0.0547*(1.8*t52_+32)+13.8;ci=tcs5+ql_*FRN*fai5*FqE/(2*pi*kc5)*log(dcs/dci);e52=(ci-tci5)/ci;tci5=ci%采用迭代算法求得包殼內(nèi)表面溫度endhg=5678;du=8.35e-3;tu5=tci5+ql_*FRN*FqE*fai5*2/(pi*(dci+du)*hg)%燃料芯塊表面溫度d5_ku=ql_*FRN*FqE*fai5/(4*pi*100);tu5_ku=(34.97-30.93)/(600-500)*(tu5-500)+30.93;to5_ku=tu5_ku+d5_ku;to5=(1100-1000)/(50.81-48.06)*(to5_ku-48.06)+1000%根據(jù)積分熱導率圖表查得芯塊中心溫度p=16e+6;hfin=1273.59e+3;hfs=1650.54e+3;hgs=2584.84e+3;G=pf_*v*3600;h5=1540.03e+3;x5=(h5-hfs)/(hgs-hfs);%該點含汽量qDNB5=3.154e6*((2.022-6.238e-8*p)+...%根據(jù)W-3公式計算出臨界熱流量(0.1722-1.43e-8*p)*exp((18.177-5.987e-7*p)*x5))*...((0.1484-1.596*x5+0.1729*x5*abs(x5))*0.2049*G/10^6+1.037)*...(1.157-0.869*x5)*...(0.2664+0.8357*exp(-124*De))*(0.8258+0.341e-6*(hfs-hfin))DNBR5=qDNB5/(q_*FRN*FqE*fai5)%計算燒毀比%第6控制體溫度計算 fai6=0.48;L6=3.8/6;e61=0.01;tf6=330;whilee61>0.001t61_=0.5*(tf6+tf5);Cp6_=1000*(0.0877*(t61_-320)+6.16);x6