📄 压水堆热工水力设计dlg.cpp
字号:
a1=495.8,//二氧化铀导热系数
a2=4.202,
a3=8.775e-14,
V0,//计算芯块中心温度的中间式
ts=42.184*pow(P,0.223)-17.95;//饱和温度
double hm[21],//单相对流放热系数
B[21],//综合物性系数
Kzro2[21];//锆表面氧化层导热系数
for(int k=1;k<20;k++)//计算各节点的中子通量分布
shuzu0[k]=(shuzu[k-1]+shuzu[k])/2;
//对不同的工况选择，为shuzu0[]选择不同的首尾数值。
switch(nCurSel4){
case 0:
shuzu0[0]=0.06;
shuzu0[20]=0.06;
break;
case 1:
shuzu0[0]=0.75;
shuzu0[14]=0.0125;
shuzu0[15]=0;
shuzu0[16]=0;
shuzu0[17]=0;
shuzu0[18]=0;
shuzu0[19]=0;
shuzu0[20]=0;
break;
case 2:
shuzu0[0]=0.0385;
shuzu0[14]=0.680;
shuzu0[15]=0;
shuzu0[16]=0;
shuzu0[17]=0;
shuzu0[18]=0;
shuzu0[19]=0;
shuzu0[20]=0;
break;
case 3:
shuzu0[0]=0.600;
shuzu0[14]=0.0715;
shuzu0[15]=0;
shuzu0[16]=0;
shuzu0[17]=0;
shuzu0[18]=0;
shuzu0[19]=0;
shuzu0[20]=0;
break;
}
for(int i=0;i<21;i++)//计算各节点的各种温度
{ B[i]=5.7497e-3+(9.43058e-5)*tfm[i]-(1.73945e-7)*tfm[i]*tfm[i]+(8.76384e-11)*tfm[i]*tfm[i]*tfm[i]+(2.5515e-11)*exp(0.0541072*tfm[i]);
hm[i]=B[i]*pow(Gm,0.8)*pow(de,-0.2);
tcsm1[i]=tfm[i]+q*shuzu0[i]/hm[i];
tcsm2[i]=ts+0.71246*pow((q*shuzu0[i]),0.5)*exp(-P/8687);
tcsm[i]=(tcsm1[i]<=tcsm2[i])?tcsm1[i]:tcsm2[i];//最终包壳外表面温度为两者的较小的一个
Kzro2[i]=1.22e-3+(1.2e-6)*tcsm[i];
tzrm[i]=tcsm[i]+q*shuzu0[i]*8e-6/Kzro2[i];
tcim[i]=(sqrt(Ka*Ka+Kb*(2*Ka*tzrm[i]+Kb*tzrm[i]*tzrm[i]+q*shuzu0[i]*dcs*log(dcs/dci)))-Ka)/Kb;
tum[i]=tcim[i]+q*shuzu0[i]*dcs/(0.5*(dci+du)*hg);
V0=log((tum[i]+a1)/a1)+a3*pow(tum[i]+273.15,4)/(4*a2)+q*shuzu0[i]*dcs/(4*a2);
double tom1=1500,tom2;//迭代求芯块中心温度
do{
tom2=tom1;
tom1=a1*exp(V0-a3*pow(tom1+273.15,4)/(4*a2))-a1;
if(tom1>4000) break;
}while(fabs(tom1-tom2)>=5);
tom[i]=tom2;}
m_edit2.SetSel(0,-1);
m_edit2.ReplaceSel("程序正在运行，请耐心等待结果的输出····");
//平均管冷却剂空泡份额与含汽量场的计算
double g1=(1.4e-2)+(9.86923e-7)*P,//试验常数，求td[21]用
Hd,//气泡脱离点的焓
Hfs=1.018314e+3+(3.951159e-2)*P+(2.038595e-12)*P*P-(9.067458e-18)*P*P*P+(2.646485e-21)*P*P*P*P,//饱和水的焓
Hgs=3.067382e+3-(3.01551e-2)*P+(1.944121e-12)*P*P-(1.437321e-16)*P*P*P-(4.898732e-21)*P*P*P*P,//饱和汽的焓
Hfg=Hgs-Hfs,//汽化潜热
A,//计算中间式
rogs=-63.65636+(1.069103e-2)*P-(2.446517e-13)*P*P-(6.947528e-18)*P*P*P-(4.702406e-22)*P*P*P*P;//饱和汽的密度
double arfa[21],//空泡份额
x1[21],//过冷沸腾1区的含汽量
td[21],//气泡脱离点的温度
rol[21],//液态水的密度
x2[21],//过冷沸腾2区的含汽量
x[21],//最终含汽量
miu[21],//水的粘度
Cp[21],//水的比热
xe[21],//饱和沸腾放热区的含汽量
td1;//标记汽泡分离点
for (int i2=0;i2<21;i2++){//求miu[],及Cp[],rol[],td[]
miu[i2]=25.3/(pow(tfm[i2]+273.15,2)+91*(tfm[i2]+273.15)-(8.58e+4));
Cp[i2]=(1.1012+(3.136e-5)*pow(tfm[i2]-220,2)+0.00243*(tfm[i2]-220)*pow(2.2,0.1*tfm[i2]-34)+(7.95e-5)*pow(18-(0.980665e-3)*P,1.152)*pow(2.4,0.05*tfm[i2]-11))*4.1868;
if(tfm[i2]<=310)//分区求rol[]
rol[i2]=1.277808e+3-1.839807*tfm[i2]+(2.02097e-7)*tfm[i2]*tfm[i2]+(4.723564e-10)*pow(tfm[i2],3)+(8.253361e-13)*pow(tfm[i2],4);
else
rol[i2]=(-1.1755984e+6)+(8.1437361e+3)*(tfm[i2]+273.15)-21.136559*pow(tfm[i2]+273.15,2)+(2.4381598e-2)*pow(tfm[i2]+273.15,3)-(1.0549747e-5)*pow(tfm[i2]+273.15,4)+(P*1000-(1.5e+7))*exp(-14.64389+(1.1283357e-3)*exp(0.012670366*(tfm[i2]+273.15)));
td[i2]=ts-g1*q*shuzu0[i2]*rol[i2]/Gm;
}
for(int s4=0;s4<21;s4++){//寻找气泡分离点
if(tfm[s4]>=td[s4])
{td1=td[s4];
break;}
else
td1=ts;
}
if(td1<=310)//求Hd
Hd=(-2.182662e+2)+5.190077*td1-(2.384186e-7)*pow(td1,2)-(5.601945e-10)*pow(td1,3)-(9.684937e-13)*pow(td1,4);
else
Hd=exp(5.830537+(4.539119e-3)*td1-(4.547474e-11)*pow(td1,2)-(1.069343e-13)*pow(td1,3)-(2.137254e-16)*pow(td1,4))+2.1735;
for(int i1=0;i1<21;i1++){
xe[i1]=(Hfm[i1]-Hfs)/Hfg;//求xe[]
if(tcsm1[i1]<=tcsm2[i1])//求单相换热区的含汽量
{
arfa[i1]=0;
x[i1]=0;
x1[i1]=0;
x2[i1]=0;
}
else {//其他区的含气量
A=(q*shuzu0[i1]-hm[i1]*(ts-tfm[i1]))*miu[i1]*Cp[i1]/(1.07*hm[i1]*hm[i1]*(ts-tfm[i1]));
arfa[i1]=4*A/de;
x1[i1]=1/(1+((1-arfa[i1])/arfa[i1])*(rol[i1]/rogs));
x2[i1]=0.435*(Hfm[i1]-Hd)/Hfg;
x[i1]=(x1[i1]<=x2[i1])?x2[i1]:x1[i1];//取两者的较大值
}
if(x[i1]<=xe[i1])//取两者的较大值
x[i1]=xe[i1];
}
//冷却剂密度场的计算
double rom[21],//各节点平均管的密度
miuw[21],//按加热壁面计算的冷却剂的密度
qm[21],//平均管各节点的热流密度
xglm[21];//平均管各节点的线功率
for(int i4=0;i4<21;i4++){//计算以上定义的各数组
rom[i4]=1/((1-x[i4])/rol[i4]+x[i4]/rogs);
miuw[i4]=25.3/(pow(tcsm[i4]+273.15,2)+91*(tcsm[i4]+273.15)-(8.58e+4));
qm[i4]=q*shuzu0[i4];
xglm[i4]=qm[i4]*Pi*dcs;
}
//压降的计算
double Pf[20],//各步长的摩擦压降
f[20],//摩擦系数
fiso[20],//初始摩擦系数
faifo[20],//两相摩擦倍数
Re[20],//雷诺数
faiq[20],//热流修正因子
Pz[20],//提升压降
Pgd[20],//隔架压降
Kgd[20],//隔架阻力系数
pm[21];//通道内压力
double fane=1.5e-6;//壁面粗糙度
for(int i5=0;i5<20;i5++){//计算以上定义的各数组
Re[i5]=Gm*de/miu[i5+1];
faiq[i5]=1+0.004355*pow(q*shuzu0[i5+1]/Gm,0.7);
faifo[i5]=1+x[i5+1]*(rol[i5+1]-rogs)/rogs;
double fiso1=0.01,fiso2;//迭代求初始摩擦系数
do{
fiso2=fiso1;
fiso1=1/pow(-2*log10(fane*de/3.7+2.51/(Re[i5]*sqrt(fiso1))),2);
}while(fabs(fiso1-fiso2)>=0.00001);
fiso[i5]=fiso1;
f[i5]=pow(miuw[i5+1]/miu[i5+1],0.6)*fiso[i5];
Pf[i5]=f[i5]*detz*Gm*Gm*faifo[i5]*faiq[i5]*0.001/(de*2*rol[i5+1]);
Pz[i5]=9.81*detz*rom[i5+1]*0.001;
Kgd[i5]=3/pow(Re[i5],0.1);
Pgd[i5]=gejia[i5]*Kgd[i5]*Gm*Gm*0.001/(2*rom[i5+1]);
}
double Pa,//延程加速压降
Pin,//入口压降
Pex,//出口压降
Pf0=0,//总摩擦压降
Pz0=0,//总提升压降
Pgd0=0,//总隔架压降
Pm;//总压降
Pa=Gm*Gm*(1/rom[20]-1/rom[0])*0.001;
Pin=3.15*Gm*Gm*0.001/(2*rom[0]);
Pex=3.46*Gm*Gm*0.001/(2*rom[20]);
for(int i8=0;i8<20;i8++){//求各总压降
Pf0+=Pf[i8];
Pz0+=Pz[i8];
Pgd0+=Pgd[i8];
}
Pm=Pf0+Pz0+Pgd0+Pa+Pin+Pex;
pm[0]=P-Pin;//求通道内压降
for(int r1=1;r1<21;r1++){
pm[r1]=pm[r1-1]-Pf[r1-1]-Pz[r1-1]-Pgd[r1-1]-Pa*detz*r1/L;
}
m_edit2.SetSel(0,-1);
m_edit2.ReplaceSel("程序正在运行，请耐心等待结果的输出······");
//热管驱动压头的计算
double Pf1[20],Pgd1[20],ph[21];//对热管的各种总压降
double f1[20],fiso0[20],Re1[20],faiq1[20],faifo1[20];//对热管的各种系数
double Phe,//热管总压降
Kfh,//修正系数
Kah,
Pa1,//对热管的各种压降
Pin1,
Pex1,
Pf10,
Pgd10,
jdz=1;//前后两次总压降的绝对值
//计算以上定义的各种量
Kfh=pow(1-0.05,1.8);
Kah=pow(1-0.05,2);
Phe=Kfh*Pf0+Kah*(Pa+Pin+Pex+Pgd0)+Pz0;
ph[0]=P-Pin*Kah;
for(int r2=1;r2<21;r2++){
ph[r2]=ph[r2-1]-Pf[r2-1]*Kfh-Pgd[r2-1]*Kah-Pz[r2-1]-Pa*Kah*detz*r2/L;
}
double Ph,Gh,Gh1;
Gh1=1500;
do{ //迭代求热管的质量流量
Pf10=0;
Pgd10=0;
for(int i9=0;i9<20;i9++){
Re1[i9]=Gh1*de/miu[i9+1];
faiq1[i9]=1+0.004355*pow(q*shuzu0[i9+1]/Gh1,0.7);
faifo1[i9]=1+x[i9+1]*(rol[i9+1]-rogs)/rogs;
double fiso3=0.01,fiso4;
do{
fiso4=fiso3;//迭代求摩擦系数
fiso3=1/pow(-2*log10(fane*de/3.7+2.51/(Re1[i9]*sqrt(fiso3))),2);
}while(fabs(fiso3-fiso4)>=0.00001);
fiso0[i9]=fiso3;
f1[i9]=pow(miuw[i9+1]/miu[i9+1],0.6)*fiso0[i9];
Pf1[i9]=f1[i9]*detz*Gh1*Gh1*faifo1[i9]*faiq1[i9]*0.001/(de*2*rol[i9+1]);
Pgd1[i9]=gejia[i9]*Kgd[i9]*Gh1*Gh1*0.001/(2*rom[i9+1]);
}
Pa1=Gh1*Gh1*(1/rom[20]-1/rom[0])*0.001;
Pin1=3.15*Gh1*Gh1*0.001/(2*rol[0]);
Pex1=3.46*Gh1*Gh1*0.001/(2*rom[20]);
for(int j1=0;j1<20;j1++){
Pf10=Pf1[j1]+Pf10;
Pgd10=Pgd1[j1]+Pgd10;
}
Ph=Pf10+Pz0+Pgd10+Pa1+Pin1+Pex1;
jdz=fabs(Ph/Phe-1);
Gh1=Gh1+0.5;
}while(jdz>=0.01);
Gh=Gh1;
m_edit2.SetSel(0,-1);
m_edit2.ReplaceSel("程序正在运行，请耐心等待结果的输出········");
//热管温度场的计算
double Hfh[21],tfh[21];
tfh[0]=t;
Hfh[0]=Hfm[0];
double Fnr=1.476,Feq=1.03,Fnu=1.05,Fea=1.05,Feb=0.89,qh;
qh=q*Fnr*Feq*Fnu;
for(int k1=0;k1<20;k1++){
Hfh[k1+1]=Hfh[k1]+detz*Pi*dcs*qh*shuzu[k1]*Fea*Feb/(Gh*Fu*Ab);
if(Hfh[k1+1]<=1623.60)//超过饱和焓温度不再上升
tfh[k1+1]=at1+at2*Hfh[k1+1]+at3*Hfh[k1+1]*Hfh[k1+1]+at4*Hfh[k1+1]*Hfh[k1+1]*Hfh[k1+1];
else
tfh[k1+1]=ts;
}
//热管燃料元件场的计算
double tcsh1[21],tcsh2[21],tcsh[21],tzrh[21],tcih[21],tuh[21],toh[21];
double hh[21],B0[21],Kzro20[21];
for(int k2=0;k2<21;k2++){
B0[k2]=5.7497e-3+(9.43058e-5)*tfh[k2]-(1.73945e-7)*tfh[k2]*tfh[k2]+(8.76384e-11)*tfh[k2]*tfh[k2]*tfh[k2]+(2.5515e-11)*exp(0.0541072*tfh[k2]);
hh[k2]=B0[k2]*pow(Gh,0.8)*pow(de,-0.2);
tcsh1[k2]=tfh[k2]+qh*shuzu0[k2]/hh[k2];
tcsh2[k2]=ts+0.71246*pow((qh*shuzu0[k2]),0.5)*exp(-P/8687);
tcsh[k2]=(tcsh1[k2]<=tcsh2[k2])?tcsh1[k2]:tcsh2[k2];
Kzro20[k2]=1.22e-3+(1.2e-6)*tcsh[k2];
tzrh[k2]=tcsh[k2]+qh*shuzu0[k2]*(8e-6)/Kzro20[k2];
tcih[k2]=(sqrt(Ka*Ka+Kb*(2*Ka*tzrh[k2]+Kb*tzrh[k2]*tzrh[k2]+qh*shuzu0[k2]*dcs*log(dcs/dci)))-Ka)/Kb;
tuh[k2]=tcih[k2]+qh*shuzu0[k2]*dcs/(0.5*(dci+du)*hg);
V0=log((tuh[k2]+a1)/a1)+a3*pow(tuh[k2]+273.15,4)/(4*a2)+qh*shuzu0[k2]*dcs/(4*a2);
double toh1=300,toh2;
do{
toh1=toh1+1;
toh2=a1*exp(V0-a3*pow(toh1+273.15,4)/(4*a2))-a1;
if(toh1>5000)
break;
}while(fabs(toh1-toh2)>=5);
toh[k2]=toh1;
}
double tohmax=0;
for(int v1=0;v1<jiedian;v1++){
if(tohmax<=toh[v1])
tohmax=toh[v1];
}
m_edit2.SetSel(0,-1);
m_edit2.ReplaceSel("程序正在运行，请耐心等待结果的输出··········");
//空泡份额与含汽量的计算
double arfa0[21],x10[21],td0[21],rol0[21],x20[21],x0[21],miu0[21],Cp0[21],xe0[21],Hd0,td01;
int lonb;
for(int k3=0;k3<21;k3++){
miu0[k3]=25.3/(pow(tfh[k3]+273.15,2)+91*(tfh[k3]+273.15)-(8.58e+4));
Cp0[k3]=(1.1012+(3.136e-5)*pow(tfh[k3]-220,2)+0.00243*(tfh[k3]-220)*pow(2.2,0.1*tfh[k3]-34)+(7.95e-5)*pow(18-(0.980665e-3)*P,1.152)*pow(2.4,0.05*tfh[k3]-11))*4.1868;
if(tfh[k3]<=310)
rol0[k3]=(1.277808e+3)-1.839807*tfh[k3]+(2.02097e-7)*tfh[k3]*tfh[k3]+(4.723564e-10)*pow(tfh[k3],3)+(8.253361e-13)*pow(tfh[k3],4);
else
rol0[k3]=(-1.1755984e+6)+(8.1437361e+3)*(tfh[k3]+273.15)-21.136559*pow(tfh[k3]+273.15,2)+(2.4381598e-2)*pow(tfh[k3]+273.15,3)-(1.0549747e-5)*pow(tfh[k3]+273.15,4)+(P*1000-(1.5e+7))*exp(-14.64389+(1.1283357e-3)*exp(0.012670366*(tfh[k3]+273.15)));
td0[k3]=ts-g1*qh*shuzu0[k3]*rol0[k3]/Gh;
}
for(int k4=0;k4<21;k4++){
if(tfh[k4]>=td0[k4])
{td01=td0[k4];
break;}
else
td01=344.787;
}
if(td01<=310)
Hd0=(-2.182662e+2)+5.190077*td01-(2.384186e-7)*pow(td01,2)-(5.601945e-10)*pow(td01,3)-(9.684937e-13)*pow(td01,4);
else
Hd0=exp(5.830537+(4.539119e-3)*td01-(4.547474e-11)*pow(td01,2)-(1.069343e-13)*pow(td01,3)-(2.137254e-16)*pow(td01,4))+2.1735;
for(int k5=0;k5<21;k5++){
xe0[k5]=(Hfh[k5]-Hfs)/Hfg;
if(tcsh1[k5]<=tcsh2[k5])
{
arfa0[k5]=0;
lonb=k5;
x0[k5]=0;
x10[k5]=0;
x20[k5]=0;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -