核反应堆热工水力课程设计解析

一、设计要求
在设计反应堆冷却系统时，为了保证反应堆运行安全可靠，针对不同的堆型，预先规定了热工设计必须遵守的要求，这些要求通常就称为堆的热工设计准则。

目前压水动力堆设计中所规定的稳态热工设计准则，一般有以下几点：
1.燃料元件芯块内最高应低于其他相应燃耗下的熔化温度；
2.燃料元件外表面不允许发生沸腾临界；
3.必须保证正常运行工况下燃料元件和堆内构件得到充分冷却；在事故工况
下能提供足够的冷却剂以排除堆芯余热；
4.在稳态额定工况和可预计的瞬态运行工况中，不发生流动不稳定性。

5.在热工设计中，通常是通过平均通道（平均管）可以估算堆芯的总功率，
而热通道（热管）则是堆芯中轴向功率最高的通道，通过它确定堆芯功率的上限，热点是堆芯中温度最高的点，代表堆芯热量密度最大的点，通过这个点来确定DNBR。

二、设计任务
某压水反应堆的冷却剂和慢化剂都是水，用二氧化铀作燃料，Zr-4作燃料包壳材料。

燃料组件无盒壁，燃料元件为棒状，正方形排列，已知下列参数：系统压力P 15.8M P a 堆芯输出热功率N t1820M W 冷却剂总流量W32500t/h 反应堆进口温度t f i n287℃堆芯高度L 3.60m 燃料组件数m121燃料组件形式n0×n017×17每个组件燃料棒数n265燃料包壳外径d c s9.5m m 燃料包壳内径d c i8.6m m 燃料包壳厚度δc0.57m m 燃料芯块直径d u8.19m m 燃料棒间距（栅距）s12.6m m 两个组件间的水隙δ0.8m m UO2芯块密度ρUO2 95%理论密度旁流系数ζ5%燃料元件发热占总发热份额F a97.4%径向核热管因子F R N 1.33轴向核热管因子F Z N 1.520
热流量核热点因子F q N F R N F Z N 2.022热流量工程热点因子F q E 1.03焓升工程热点因子FΔH E未计入交混因子) 1.142
交混因子FΔH·m
E0.95焓升核热管因子FΔH N F R N 1.085堆芯进口局部阻力系数K i n0.75堆芯出口局部阻力系数K o u t 1.0堆芯定位格架阻力系数K g r 1.05
若将堆芯自下而上分为3个控制体，其轴向归一化功率分布见下表：
通过计算，得出：
1. 堆芯流体出口温度；
2. 燃料棒表面平均热流密度以及最大热流密度，平均线功率，最大线功率；
3. 热管内的流体温度（或焓）、包壳表面温度、芯块中心温度随轴向的分布；
4. 包壳表面最高温度，芯块中心最高温度；
5. DNBR 在轴向上的变化；
6. 计算堆芯压降
三、设计正文（详细的计算过程、计算结果及分析）
1.计算过程
1.1堆芯流体出口温度（平均管）
t f,out=t f,in+
F a∙N t
W∙(1−ζ)∙C p̅̅̅
C p̅̅̅按流体平均温度t f̅=12(t f,in+t f,out)以及压力由表中查得。

假设t f,out=330℃，查表得C p̅̅̅=5.610kJ/(kg∙℃)
经过输入所查C p̅̅̅程序不断迭代得t f,out=323.9℃误差小于0.5℃。

如需更精确的值，可以继续进行迭代计算。

1.2燃料表面平均热流密度q̅
q̅=F a∙N t/F
总
式中F
总
为堆芯燃料棒的总传热面积
F
总
=m∙n∙π∙d cs∙L
代入数据得
F
总
=121×265×π×9.5×10−3×3.60=3443.40 m2
q̅=97.4%×1820×106
3443.40
=5.29×105 W/m2
燃料棒表面最大热流密度q max
q max=q̅∙F q N∙F q E
代入数据得
q max=5.29×105×2.022×1.03=1.10×106 W/m2
燃料棒平均线功率q l̅
q l̅=q̅π∙d cs∙L
L
=q̅∙π∙d cs
代入数据得
q l̅=5.29×105×π×9.5×10−3=1.57×104 W/m
燃料棒最大线功率q l,max
q l,max=q l̅∙F q N∙F q E
代入数据得
q l,max=q l̅∙F q N∙F q E=1.57×104×2.022×1.03=3.26×104 W/m
1.3平均管的情况
平均管的流速V
V=W(1−ζ)
t f
式中A t堆芯内总流通面积
A t=m∙(n0×n0)[s2−π
d cs2]+m[4(n0∙s×
δ
)]
n0为燃料组件内正方形排列时的每一排（列）的燃料元件数ρf̅̅̅由压力以及流体的平均温度t f̅查表得到:
ρf̅̅̅=1 v f
由1.1知t f̅=323.9+287
2
=305.5℃，查表得v f= 0.001397680614m3/kg
ρf̅̅̅=
1
0.001397680614
=715.471 kg/m3
A t=121×(17×17)×[(12.6×10−3)2−π
4
×(9.5×10−3)2]
+121×(4×17×12.6×10−3×0.8×10−3
2
)=3.11 m2
V =
32500×(1−5%)
3.11×715.471×3.6
=3.85 m/s
1.4为简化计算起见，假定热管内的流体流速V h 和平均管的V 相同。

同样，热管四根燃料元件组成的单元通道内的流量
W b =W(1−ζ)A t
A b
A b =s 2
−π4
d cs 2
代入数据得
A b =(12.6×10−3)2−
π
4
(9.5×10−3)2=0.88×10−4 m 2 W b =32500×（1−5%）3.11
×0.88×10−4=0.87 t/ℎ
1.5热管中的计算（按一个单元通道计算） （1）热管中的流体温度
t f (z)=t f,in +q ̅∙F R N
∙F ΔH E ∙F ΔH∙m E πd cs b p ∫φ(z)dz z
t f (z)=287+5.29×105×1.35×1.142×0.95π×9.5×10−30.87×10003600
∙C p
∫φ(z )dz
z
0=287+9.56×104C p ∫φ(z )dz z
其中C p 取平均温度对应的参数值，需要进行迭代计算， 下面给出第一控制体出口处温度的算法
假设t f (L
6)=300℃,查表得C p =5.3348 kJ/(kg ∙℃),带入上式
t f (L 6)=287+9.56×1045.3348×1000×0.8×3.60
6
=291.6 ℃ 与假设误差较大，进行迭代，查表知C p =5.222 kJ/(kg ∙℃)
t f (L 6)=287+9.56×1045.222×1000×0.8×3.663
=291.69 ℃ 误差|291.60−291.69|<0.5℃，可以不再进行迭代，就取t f (L
6)=291.69 ℃ 同理由程序迭代可求得
第二控制体出口处流体温度t f (2L
6)=301.38 ℃
第三控制体出口处流体温度t f (3L
6)=314.86 ℃ 第四控制体出口处流体温度t f (4L
6)=327.70 ℃ 第五控制体出口处流体温度t f (5L 6)=334.89℃ 第六控制体出口处流体温度t f (L )=338.22 ℃ 2)第一个控制体出口处的包壳外壁温度
t cs (z)=t f (z)+Δθf1(z)=t f (z)+q ̅∙F R N ∙φ(z )∙F q E
()
式中：h z)为单相水强迫对流换热系数[W/(m 2∙℃)]，可以利用以下公式来求
N u =h (z )D ε
λ
=0.023Re 0.8∙Pr 0.4
所以
h (z )=0.023Re 0.8Pr 0.4∙
λD ε
式中
Re =G ∙D ε=W b b ∙D ε
D ε=4A b U =4(s 2
−π4d cs 2)πd cs
流体的λ、μ和Pr 数根据流体的压力和温度由表查得。

如果流体已经达到过冷沸腾，用Jens-Lottes 公式：
Δθf2(z)=t s +25(q ̅F R N F q E φ(z
)6)0.25
∙
e −p
6.2
−t f (z)
其中t s 为气体的饱和温度，p 的单位为MPa ，p =15.8MPa 时，t s =346.38℃ 当Δθf2≥Δθf1时，用前面的式子 当Δθf2<Δθf1时，用Δθf2替换掉Δθf1 代入数据得
D ε=
4×((12.6×10−3)2−π
4×(9.5×10−3)2)
π×9.5×
10−3
=11.78×10−3 m
Re =0.86
3.60.88×10−4×11.78×10−3μ=32.35μ
h (z )=0.023×(32.35)0.8Pr 0.4λ−3=31.51λPr 0.4
0.8
Δθf1(z )=5.29×105×1.33∙φ(z )∙1.03ℎ(z )=2.25×104μ0.8φ(z
)λPr 0.4
Δθf2(z )=346.38+25(5.29×105×1.33×1.10φ(z )
6
)
0.25
∙
e −15.86.2
−t f (z )
=346.38+1.83∙φ(z )0.25−t f (z )
第一控制体出口处
t f (L
6)=291.69 ℃，查表可得λ=0.57799W/(m ∙℃) φ(L
6)=0.48 μ=9.1858×10−5 kg/(m ∙s) Pr =0.8371
Δθf1(L 3)=2.25×104
×(9.1858×10−5)0.8×0.80.57799×0..83710.4=11.82℃
Δθf2(L
3
)=346.38+1.82×0.480.25−303.91=56.21 ℃
故
t cs (L )=t f (L )+Δθf1(L
)=291.69+11.82=303.51 ℃
其余同理由程序计算得出结果如下 第二控制体出口处
t f (2L
6)=301.38 ℃，查表可得λ=0.5601W/(m ∙℃) φ(2L
6)=1.02 μ=8.8×10−5 kg/(m ∙s) Pr =0.8621
Δθf1(2L
6)=24.79 ℃
Δθf2(2L
6
)=46.84 ℃
故
t cs (2L )=t f (2L )+Δθf1(2L
)=301.38+24.79=326.15 ℃
第三控制体出口处
t f (3L 6)=314.86℃，查表可得λ=0.5333 W/(m ∙℃) φ(L )=1.50 μ=8.275×10−5 kg/(m ∙s) Pr =0.9157
Δθf1(L )=35.55 ℃ Δθf2(L )=33.55 ℃
故
t cs(3L
6
)=t f(L)+Δθf2(L)=348.41 ℃
第四控制体出口处
t f(4L
6
)=327.70℃，查表可得λ=0.5044 W/(m∙℃) φ(L)=1.56μ=7.735×10−5 kg/(m∙s)Pr=1.0051
Δθf1(L)=35.67 ℃
Δθf2(L)=20.73 ℃
故
t cs(4L
6
)=t f(L)+Δθf2(L)=348.43 ℃
第五控制体出口处
t f(4L
6
)=334.89℃，查表可得λ=0.4868W/(m∙℃) φ(L)=0.96μ=7.401×10−5 kg/(m∙s)Pr=1.0873
Δθf1(L)=21.28 ℃
Δθf2(L)=13.30℃
故
t cs(5L
6
)=t f(L)+Δθf2(L)=348.19 ℃
第六控制体出口处
t f(L)=338.22℃，查表可得λ=0.4782 W/(m∙℃) φ(L)=0.48μ=7.235×10−5 kg/(m∙s)Pr=0.9157
Δθf1(L)=11.39 ℃
Δθf2(L)=9.68℃
故
t cs(L)=t f(L)+Δθf2(L)=347.90 ℃
（3）包壳内壁温度
t ci(z)=t cs(z)+q l̅F R N F q Eφ(z)
2πk c(z)
ln
d cs
d ci
式中Zr-4的k c=0.00547×(1.8×t c̅+32)+13.8 W/m∙℃
t c̅(z)=1
2
[t cs(z)+t ci(z)]
代入数据得：
t ci(z)=t cs(z)+1.57×104×1.33×1.03φ(z)
c
()
ln
9.5
由于k c与平均温度有关，由程序迭代计算结果如下
第一控制体出口处 t ci (1L
6)=313.13 ℃ 第二控制体出口处 t ci (2L
6)=346.26 ℃ 第三控制体出口处 t ci (3L 6)=377.54 ℃ 第四控制体出口处 t ci (4L 6)=378.72℃ 第五控制体出口处 t ci (5L 6)=366.89 ℃ 第六控制体出口处 t ci (L )=357.27 ℃ 4)燃料芯块外表面温度
t u (z)=t ci (z)+q l ̅F R N F q E φ(z)
πd ci +d u 2
∙ℎg
式中ℎg 是包壳与芯块间的气隙等效传热系数，这里取ℎg =5678 W/(m 2∙℃) 代入数据得
t u (z )=t ci (z )+1.57×104×1.33×1.03φ(z)
π8.6+8.192×5678×10
−3=t ci (z )+143.70φ(z)
第一个控制体出口处
t u (L
6
)=313.13+143.70×0.48=384.03 ℃
第二个控制体出口处
t u (2L
6
)=346.26+143.70×1.02=492.83℃
第三个控制体出口处
t u (3L
)=377.54+143.70×1.50=593.09℃
第四个控制体出口处
t u (4L
6
)=378.72+143.70×1.56=602.89 ℃
第五个控制体出口处
t u (5L
6
)=366.89+143.70×0.96=504.84 ℃
t u (L )=357.27+143.70×0.48=426.25℃
5)燃料芯块中心温度用积分热导求解的方法，即
∫k u (t )dt t 0(z )0=∫k u (t )dt t u (z )0
+q l ̅F R N F q E φ(z
)4π
其中 k u =38.24
t+402.55+4.788×(t +273.15)3 代入数据得
∫
k u (t )dt t 0(z )
=∫
k u (t )dt t u (z )
0+
1.57×104×1.33×1.03φ(z )
4π
=∫
k u (t )dt t u (z )
+17.12φ(z) W/cm
∫k u (t )dt t
=∫
38.24
t +402.55
+4.788×(t +273.15)3dt
t
=38.24ln t +402.55402.55+4.788×10−13
4
×(t +273.15)4−273.154
令 f (x )=∫k u (t )dt t 0
(z )
0−∫k u (t )dt t u
(z )
0+17.12φ(z)
由于函数递增，可以通过二分法求解f x)的根得出二氧化铀中心温度
通过编程可求得结果如下 第一控制体出口处
t 0(L
6)=568.79 ℃
第二控制体出口处
t 0(2L 6)=974.83 ℃
第三控制体出口处
t 0(3L 6)=1383.57 ℃
第四控制体出口处
t 0(4L
6)=1430.38℃
第五控制体出口处
t 0(5L 6)=958.77℃
t 0(L )=619.90 ℃
1.6热管中的q DNB (z)
用w-3公式计算，同样对3个控制体都算
q DNB =3.154×106{(2.022−6.238×10−8p )+(0.1722
−1.43×10−8p)
×exp [(18.177−5.987×10−7p )x e ]}[(0.1484−1.596x e +0.1729x e |x e |)×0.2049G
106
+1.037](1.157
−0.869x e )[0.2664+0.8357exp (−124D ε)]×[0.8258+0.341×10−6(ℎfs −ℎf,in )]
式中：p 为冷却剂工作压力（Pa ），G 为冷却剂质量流密度[kg/(m 2∙ℎ)],D ε为冷却剂通道的当量直径（m ）,ℎfs 为冷却剂的饱和比焓（J/kg ）,ℎf,in 为控制体进口处冷却剂的比焓（J/kg ）,x e 为计算点z 处的平衡含气量，|x e |为其绝对值。

平衡含气量x e 的计算式为
x e =ℎf (z)−ℎfs
ℎfg
其中ℎfg 为汽化潜热（J/kg ）。

ℎfs =1641.4kJ/kg ℎfg =945kJ/kg p =15.8MPa
D ε=11.78×10−3m
G =W b A b =0.87×10000.88×10
−4=9.88×106kg/(m 2∙ℎ) 通过程序计算得出结果如下 第一控制体出口处
q DNB (L
6
)=5.244×106W/m 2
第二控制体出口处
q DNB (2L
6
)=4.731×106W/m 2
第三控制体出口处
q DNB (3L
6
)=3.897×106W/m 2
第四控制体出口处
q DNB (4L
6
)=3.109×106W/m 2
第五控制体出口处
q DNB (5L
6
)=2.624×106W/m 2
第六控制体出口处
q DNB (L )=2.388×106W/m 2
1.7DNBR 的计算
DNBR(z)=q DNB (z)
R N q E
第一控制体出口处
DNBR (L 6)= 5.244×106
5.29×105×1.33×1.03×0.48
=15.077
同理可求其余段DNBR ，结果如下： 第二控制体出口处
DNBR (2L
)=6.40
第三控制体出口处
DNBR (3L
6
)=3.585
第四控制体出口处
DNBR (4L
6
)=2.750
第五控制体出口处
DNBR (5L
6
)=3.772
第六控制体出口处
DNBR (L )=6.864
1.8计算热管中的压降 单相流体的摩擦压降
∆p f =f ∙
L D ε
∙
ρ̅V 2
2
式中：
f =f iso (μw f )n
=0.1840.2∙(μw f
)0.6
使用公式编程分别计算六段控制体的摩擦压降
单相流体加速压降：
∆p a=G2(v out−v in)
同样使用公式编程分别计算六段控制体的加速压降
单相流体提升压降
Δp el=ρ̅∙g∙L
同上
局部压降，出口：
∆p out=k out ρout V out2
2
=k out
G2v out
2
代入数据得
∆p out=1.0×（9.88×106
3600）
2
×
0.00166767
2
=6.03×103 Pa
局部压降，进口：
∆p in=k in ρin V in2
2
=k in
G2v in
2
代入数据得
∆p in=0.75×（9.88×106
3600）
2
×
0.00132935
2
=5.0×103 Pa
局部压降，定位格架出口压降
∆p gr=k gr ρgr V̅2
=k gr
G2
1
2(v out+v in)
代入数据得
∆p gr=1.05×（9.88×106
3600）
2
×
0.00166767+0.00132935
4
=5.79×103 Pa
以上所使用的比热容和动力粘度都通过软件查询后输入代码中进行计算。

通过程序计算结果如下
分段压降变化
总压降为
∑p=∆p out+∆p gr+∆p in+Δp el+∆p f+∆p a=4.6981×104Pa
2.计算总结
温度汇总表
临界热流密度和烧毁比汇总
所以理论上能够保证安全性。

四、课程设计感想
通过这次反应堆热工分析的课程设计，我加深了对反应堆内部传热的了解，同时我发发现了自己的很多不足之处。

拿到课程设计题目后，我首先根据题目中的问题去复习热工分析的知识，去了解每个问题该怎么算。

然后将计算思路和公式写出来，为下一步的编程做准备。

由于大部分的求解都是非线性的方程，所以采用了迭代和二分法求解方程。

当把整个问题的求解思路理清和确定计算过程中使用算法后，整个问题就在水和水蒸气热物性如何导入上了。

在最开始，我使用通过一定的数据去拟合用到的热物性在定压下与温度的关系。

但在拟合之后，通过得出的结果与接用水和水蒸气热物性查询软件得出的结果有一点差距，所以就放弃了。

然后，我想将网上通过水和蒸汽热力性质公式IAPWS_IF97计算的源码做成接口，在计算过程中直接用它计算热物性。

后来，发现太难，源码有很多地方看不懂，不知道如何接入。

于是最终也放弃了。

后来就采用了最原始的方法，就是在计算的过程中，通过水和水蒸气热物性查询软件查出结果，一步一步输进去。

最后使用程序计算时，在输入参数时耗时太久，我深深的感受到书到用处方恨少，如果原来，把编程多学一点，将整个数据导入直接通过动态链接库完成，就省事太多。

这次课设让我把原来学的编程与实际问题结合起来，给我很大感触，原来敲着书上的例子，感觉漫无目的。

通过这次自己去思考如何解决一个问题，让我理解了拿到一个问题后，该如何处理。

同时通过这次课设，我意识到反应堆热工分析是个复杂而连续的过程，每个参数都受到大量的常数参数的影响也具有很多的修正因子。

我们不应该根据自己的常识来判断数据的变化情况，相反地我们应该随时坚持以数据计算为引导，以实验作为验证。

仔细客观认真地分析堆内数据的变化，并且对堆内数据进行全程监控，防止堆内数据随时变化，对反应堆的危害性。

附录（设计流程图、程序）
1、程序说明
共7个程序由c语言编写，分别为计算流体出口温度，控制体出口流体温度，燃料包壳外壁温度，燃料包壳内壁温度，二氧化铀中心温度，qDNB,压降。

使用了迭代和二分法求解部分非线性方程。

（1）堆芯出口温度计算：
此段根据任务书给出的基本参数和热量与流量之间关系，运用迭代的算法，求出堆芯的出口温度。

（2）第一至第六控制体的各量计算：
因为六个控制体的计算过程类似，这里只说明第一个控制体的计算过程。

在现有的参数下，根据热流量与流量的关系和迭代算法，求出该控制体的出口温度。

通过流通截面积与湿周的关系求出栅元的当量直径。

再根据上面的温度，查出对应的热物性参数由雷诺数与努尔数的关系，解出控制体出口处的对流换热系数。

因为不知该处的流体状态，分别用单相强迫对流放热公式和詹斯-洛特斯传热方程算出各自的膜温压，取较小的值加上出口处的流体温度即是包壳的外表面温度。

由包壳的外表面的温度再根据圆管的传热方程运用迭代算法解出包壳内表面的温度。

芯块与包壳内表面之间的导热问题，根据间隙导热模型，即可解出芯块表面的温度，根据内热源的导热模型，依据积分热导率与温度的对应关系列出方程用二分法解出芯块中心的温度。

接下来依据冷却剂的温度，得出的控制体出口处的含汽量。

进而依据W-3公式求出该出的临界热流量qDNB，最后得出该出的烧毁比DNBR。

（3）热管的压降计算：
热管的压降包括摩擦压降、提升压降、进出口局部压降、定位搁架出口压降。

摩擦压降可由计算单相流的达西（Darcy）公式算得。

提升压降可由根据位置的变化算得，其中参数都取平均值。

其余的压降根据形阻压降的基本公式再乘以相应的系数求得。

最后各项相加得出热管的总压降。

2、流程图
2、程序
// 3.1.cpp : 定义控制台应用程序的入口点。

//流体出口温度迭代计算
#include"stdafx.h"
#include<math.h>
#include<iostream>
using namespace std;
double tfout(double cp, double atfout){
double tfin = 287, Fa = 0.974,Nt = 1820,W=32500,plxs=0.05,tfout;
tfout = tfin + 3600*Fa*Nt / (W*cp*(1 - plxs));
return tfout;
}
int main(){
double atfout,tfin=287,tav,cp,tfoutc=0;
int panding = 1;
printf("请输入假设出口温度\n");
cin >> atfout;
while(panding>=1){
tav = (atfout + tfin) / 2;
printf("请输入%f下的Cp\n", tav);
cin >> cp;
tfoutc = tfout(cp, atfout);
if (fabs(tfoutc - atfout) <= 0.5)
panding = 0;
atfout = tfoutc;
}
printf("出口温度为%f\n",tfoutc);
while (true)
{
}
}
//3.3.1热管内流体温度计算
#include"stdafx.h"
#include<math.h>
#include<iostream>
using namespace std;
double ffout(double tfin,double atfout){
double qav, fnr=1.33, feh=1.03, fehm=0.95, dcs=9.5, gyh, cp, wh,l=3.60,tfout;
bool panding = true;
cout << "请输入平均热流密度（W/m*m） " << endl;
cin >> qav;
cout << "请输入归一化参数 " << endl;
cin >> gyh;
cout << "请输入热管单元通道流量Wh（t/h）" << endl;
cin >> wh;
while (panding){
double tav = (atfout + tfin) / 2;
cout << "请输入" << tav << "下的Cp(J/(kg.℃))" << endl;
cin >> cp;
tfout = tfin + 3.6*(qav*fnr*feh*fehm*3.14*dcs*gyh*l / 6) / (wh*cp*1000);
if (fabs(tfout - atfout) <= 0.5)
panding = false;
atfout = tfout;
}
return tfout;
}
int main(){
double tfin,atfout,tfout;
cout << "请输入流体进口温度" << endl;
cin >> tfin;
cout << "请输入流体假设流体出口温度" << endl;
cin >> atfout;
tfout = ffout(tfin, atfout);
cout << tfout;
while (true)
{
}
}
// 3.3.02.cpp : 定义控制台应用程序的入口点。

//计算包壳外壁温度
#include"stdafx.h"
#include"cmath"
#include"iostream"
using namespace std;
int main(){
double tfout[6] = { 291.69,301.38, 314.86, 327.70, 334.89, 338.22 }, gyhcs[6] = {0.48,1.02,1.50,1.56,0.96,0.48};
double pr[6],drxs[6],dlnd[6],f1[6],f2[6],f[6];
int i;
for (i = 0; i < 6; i++)
{
cout << "请输入" << tfout[i] << "下的普朗特数" << endl;
cin >> pr[i];
cout << "请输入" << tfout[i] << "下的导热系数（W/(m∙℃)）" << endl;
cin >> drxs[i];
cout << "请输入" << tfout[i] << "下的动力粘度kg/(m∙s)" << endl;
cin >> dlnd[i];
f1[i] = 2.25*pow(10, 4)*gyhcs[i] * pow(dlnd[i], 0.8)/(drxs[i]*pow(pr[i],0.4));
f2[i] = 346.38 + 1.83*pow(gyhcs[i], 0.25) - tfout[i];
if (f1[i]>f2[i])
f[i] = f2[i];
else
f[i] = f1[i];
cout << "θf1=" << f1[i] << "θf2=" << f2[i] <<"f="<< f[i] + tfout[i] << endl;
}
while (true)
{
}
}
// 3.3.2.cpp : 定义控制台应用程序的入口点。

//计算内壁温度
#include"stdafx.h"
#include"iostream"
#include"math.h"
using namespace std;
double ftci(double tcs,double atci){
double frn = 1.33, feq = 1.03, gyhcs, dcs=9.5, dci=8.6, kc, tav,q,tci;
bool panding = true;
cout << "请输入归一化参数" << endl;
cin >> gyhcs;
cout << "请输入平均线功率，单位w/m" << endl;
cin >> q;
while (panding)
{
tav = (tcs + atci) / 2;
kc = 0.00547*(1.8*tav + 32) + 13.8;
tci = tcs + log(dcs/dci)*(q*frn*feq*gyhcs) / (2 * 3.14*kc);
if (tci - atci <= 0.5)
panding = false;
atci = tci;
}
return tci;
}
int main(){
double tcs, tci,atci;
cout << "请输入外壁温度℃" << endl;
cin >> tcs;
cout << "请输入内壁假设温度℃"<<endl;
cin >> atci;
tci = ftci(tcs, atci);
cout << tci;
while (true)
{
}
}
// 3.3.3.cpp : 定义控制台应用程序的入口点。

//计算二氧化优中心温度
#include"stdafx.h"
#include"iostream"
#include"math.h"
#include<cmath>
using namespace std;
double f1(double t){
double fout;
fout = 38.24*log((t+402.55)/402.55)+4.788*pow(10, -13)*(pow((t + 273.15), 4) - pow(273.15, 4));
return fout;
}
double f2( double tci, double atu)
{
double tu, fnr = 1.33, feq = 1.03, gyhcs, q, tmid, jieguo;
bool panding = true;
cout << "请输入归一化参数" << endl;
cin >> gyhcs;
cout << "请输入线热流密度/w" << endl;
cin >> q;
if ((f1(atu) - f1(tci) - (q*fnr*feq*gyhcs) / (4 * 314)) > 0)
{
tu = tci;
panding = false;
}
while (panding)
{
if ((f1(atu) - f1(tci) - (q*fnr*feq*gyhcs) / (4 * 314)) > 0)
{
panding = false;
tu = atu - 100;
}
else
atu = atu + 100;
}
while (atu-tu>0.0000001)
{
tmid = (atu + tu) / 2;
jieguo = f1(tmid) - f1(tci) - (q*fnr*feq*gyhcs) / (4 * 314);
if (jieguo == 0)
tu = tmid;
else
if (jieguo > 0)
atu = tmid;
else
tu = tmid;
}
return tu;
}
int main()
{
double tci, atu, tu;
cout << "请输入燃料芯块外壁温度/℃" << endl;
cin >> tci;
cout << "请输入假设二氧化铀中心温度/℃" << endl;
cin >> atu;
tu=f2(tci, atu);
cout << tu;
while (true)
{
}
}
// 3.5.cpp : 定义控制台应用程序的入口点。

//计算DNBR
#include"stdafx.h"
#include"iostream"
#include"cmath"
#define pi 3.14
#define P 15800000
using namespace std;
double fqDNB(double hfout, double hfin, double hfs, double hfg, double g, double dlzj){ double x, a1, a2, a3, a, b, c, d, e, qDNB;
x = (hfout - hfs) / hfg;
a1 = 2.022 - 6.238*pow(10, -8)*P;
a2 = 0.1722 - 1.43*pow(10, -8)*P;
a3 = exp((18.177 - 5.987*pow(10, -7)*P)*x);
a = 3.154 * pow(10, 6) * (a1 + a2*a3);
b = (0.1484 - 1.596*x + 0.1729*x*fabs(x))*0.2049*g / 1000000 + 1.037;
c = 1.157 - 0.869*x;
d = 0.2664 + 0.8357*exp(-124 * dlzj);
e = 0.8258 + 0.341*pow(10, -6)*(hfs - hfin);
qDNB = a*b*c*d*e;
return qDNB;
}
int main(){
double hfout[6], hfin[6], gyhcs[6],qDNB[6],DNBR[6], hfs, hfg, dcs = 0.0095, g, s = 0.0126, q, frn=1.33, fqe=1.03,dlzj;
int i;
dlzj = 4 * (s*s - pi*dcs*dcs / 4) / (pi*dcs);
for (i = 0; i < 6; i++)
{
cout << "请输入第" << i+1 << "段控制体出口冷却剂比焓（J/kg）" << endl;
cin >> hfout[i];
cout << "请输入第" << i +1<< "段进口温度下的冷却剂比焓（J/kg）" << endl;
cin >> hfin[i];
cout << "请输入第" << i+1 << "段控制体归一化参数" << endl;
cin >> gyhcs[i];
}
cout << "请输入15.8MPa下的冷却剂的饱和比焓（J/kg）" << endl;
cin >> hfs;
cout << "请输入15.8MPa下的冷却剂的汽化潜热（J/kg）" << endl;
cin >> hfg;
cout << "请输入热管中冷却剂质量流量密度[kg/(m*m*h)]" << endl;
cin >> g;
cout << "请输入线热流密度(W/m*m*m)" << endl;
cin >> q;
for (i = 0; i <6; i++)
{
qDNB[i] = fqDNB(hfout[i], hfin[i], hfs, hfg, g, dlzj);
DNBR[i] = qDNB[i]/ (q*fqe*frn*gyhcs[i]);
cout <<"第"<<i+1<<"段控制体qDNB为"<< qDNB[i] << endl;
cout << "第" << i + 1 << "段控制体DNBR为" << DNBR[i] << endl;
}
while (true)
{
}
}
// 3.6.cpp : 定义控制台应用程序的入口点。

//计算堆芯压降
#include"stdafx.h"
#include"iostream"
#include"cmath"
using namespace std;
const double g= 9.8;
const double Ko = 1.0;
const double Ki = 0.75;
const double Kgr = 1.05;
const double L = 3.60;
const double V = 3.85;
const double G = 9.88e+6/3600;
const double dlzj = 1.178e-2;
const double tfout[6] = { 291.69, 301.38, 314.86, 327.70, 334.89, 338.22 }; const double tfin[6] = { 287, 291.69, 301.38, 314.86, 327.70, 334.89 }; const double tcs[6] = { 303.51, 326.15, 348.41, 348.43, 348.19, 347.90 }; double fpel(double midu) //计算提升压降
{
return midu*L*g/6;
}
double fpf(double f, double dlzj, double midu) //计算摩擦压降
{
return f*L*midu*V*V / (2 * dlzj*6);
}
double fpa( double vi, double vo) //计算加速压降
{
return G*G*(vo - vi);
}
double fpout(double vo) //计算出口局部压降
{
return Ko*G*G*vo / 2;
}
double fpin( double vi) //计算进口局部压降
{
return Ko*G*G*vi / 2;
}
double fpgj(double vi,double vo) //计算格架出口压降
{
return Kgr*G*G*(vo + vi) / 4;
}
double fmxs(double re, double uw, double uf) //计算摩擦系数
{
return 0.184*pow(uw / uf, 0.6) / pow(re, 0.2);
}
double fre(double dlzj, double midu, double u) //计算雷诺数
{
return dlzj*midu*V / u;
}
int main(){
double pel[6], pf[6], pa[6], pout, pin, pgj, mxs[6], re[6], midu[6], p[6], pt = 0;
int i;
double birong[6] = { 0.001337049414, 0.001362247294, 0.001409055895,
0.001475639644, 0.001541851039, 0.001585457905 };
double vi[6] = { 0.001329349786, 0.001344993020, 0.001380774324, 0.001440939652, 0.001515991823, 0.001570780783 };
double vo[6] = { 0.001344993020, 0.001380774324, 0.001440939652, 0.001515991823, 0.001570780783, 0.001601107797 };
double uw[6] = { 8.724724e-005, 7.803018e-005, 2.328906e-005, 2.328904e-005,
2.328945e-005, 2.329043e-005 };
double uf[6] = { 9.277614e-005, 8.996754e-005, 8.543740e-005, 8.011150e-005,
7.572227e-005, 7.319442e-005 };
for (i = 0; i < 6; i++)
{
midu[i] = 1 / birong[i];
re[i] = midu[i]*V*dlzj/uf[i];
mxs[i] = 0.184*pow(uw[i] / uf[i], 0.6)/pow(re[i],0.2);
}
pout = fpout(vo[5]);
pin = fpin(vi[0]);
pgj = fpgj(vi[0], vo[5]);
for (i = 0; i < 6; i++)
{
pel[i] = fpel(midu[i]);
pf[i] = mxs[i] * L*midu[i] * V / (dlzj * 6 * 2);
pa[i] = fpa(vi[i],vo[i]);
p[i] = pel[i] + pf[i] + pa[i];
pt = pt + p[i];
}
pt = pt + +pout + pin + pgj;
for (i = 0; i < 6; i++)
{
cout <<i+1<<"段提升"<< pel[i] << endl;
cout << i+1<< "段摩擦"<<pf[i] << endl;
cout << i+1 << "段加速"<< pa[i] << endl;
}
cout <<"出口"<< pout << endl;
cout <<"入口"<< pin << endl;
cout <<"格架"<< pgj << endl;
cout << pt << endl;
while (true)
{
}
}。