Thermal hydraulic power.rar_压水堆_压水堆堆芯平均管热工水力计算_堆芯_流动 热_热处理

#include<stdio.h> #include<WASPCN32.h> #include<math.h> static double pi=3.14159; void main() { void xierewenjian(int n,double *shu,char **name); double moody(double re,double k); void xiewenjian(int n,double *delt_p); double nt=2890000,w=32100000,p=15.5,t_init=289,t_in=289,t_out,t_av,cp=0,py[3]={0.8,1.5,0.7},ql,qs,z=3.66; double h_in=1279,h,wb,ab,t_cs,t_cs_up,t_ci,t_us,t_u,di,DNBR; double re,delta_p_in,delta_p_out,delta_p_ga,x_e,G,delta_p[4][3]; double shu[10][3]; char *name[10]; int i,range; w=w/3600; /*****************************堆芯出口温度计算*****************************************/ { double av_cp,av_cha_t=0,t_out_av=344; do { t_av=(t_init+t_out_av)/2; PT2CP97(p,t_av,av_cp,range); { double t_out1_av; t_out1_av=(t_init+0.974*nt/(av_cp*w*(1-0.05))); av_cha_t=fabs(t_out_av-t_out1_av); t_out_av=t_out1_av; } }while(av_cha_t>0.2); printf("出口温度:%f\n比容:%f\n",t_out_av,av_cp); /**********************平均管的热流密度和单通道流量计算**************************/ { double a_tr,l,qsm,qlm; a_tr=121*265*pi*0.0095*3.66; qs=0.974*nt/a_tr; qsm=qs*1.35*1.528*1.03; l=121*265*3.66; ql=nt*0.974/l; qlm=ql*1.35*1.528*1.03; printf("堆芯总的换热面积:%f\n表面平均热流密度:%f\n表面最大热流密度:%f\n",a_tr,qs,qsm); printf("元件总长:%f\n表面平均线热流密度:%f\n表面最大线热流密度:%f\n",l,ql,qlm); { double af,pf_av,V,v_init,v_out_av; af=121*17*17*(0.0126*0.0126-pi*0.0095*0.0095/4)+121*(4*17*0.0126*0.0008/2); PT2V(p,t_init,v_init,range); PT2V(p,t_out_av,v_out_av,range); pf_av=1/((v_init+v_out_av)/2); V=w*(1-0.05)/(af*pf_av); ab=0.0126*0.0126-pi*0.0095*0.0095/4; wb=w*(1-0.05)*ab/af; printf("堆芯流通面积:%f\n平均管平均密度:%f\n",af,pf_av); printf("平均管流速:%f\n单通道面积:%f\n单通道流量:%f\n\n",V,ab,wb); } } } h=h_in; //控制体进口焓初始化 /**************************控制体，循环体****************************/ for(i=1;i<4;i++) { /**************************流体流过控制体后温度****************************/ { double cha_f=2,t_out1; t_out=t_in+10; for(;cha_f>0.01;) { t_av=(t_in+t_out)/2; PT2CP97(p,t_av,cp,range); t_out1=t_in+qs*1.35*1.085*0.95*pi*0.0095*py[i-1]/(wb*cp);//0.242 cha_f=fabs(t_out-t_out1); t_out=t_out1; } h+=ql*1.35*1.03*py[i-1]*z/(3*wb); //控制体焓值每步叠加 shu[0][i-1]=t_out; name[0]="t_f"; shu[1][i-1]=h; name[1]="h"; printf("第%d步长：\n平均温度%f，出口温度:%f，定压比容:%f\n",i,t_av,t_out,cp); } /*************************包壳外壁面温度*********************************/ { double H,k,pr,t_s,delta_theta,u_cs; P2T97(p,t_s,range); //15.5MP下的饱和温度 di=4*(0.0126*0.0126-pi*0.0095*0.0095/4)/(pi*0.0095); if(i==1) { double delta_tcs_init; PT2ETA97(p,t_init,u_cs,range); re=wb*di/(u_cs*ab); PT2RAMD(p,t_init,k,range); PT2PRN(p,t_init,pr,range); H=0.023*k*pow(re,0.8)*pow(pr,0.4)/(1000*di); //包壳表面换热系数h delta_tcs_init=qs*1.35*1.03*0.01/H; //未沸腾下的温升 t_cs_up=t_init+delta_tcs_init; } PT2ETA97(p,t_av,u_cs,range); re=wb*di/(u_cs*ab); printf("动力粘度=%e，",u_cs); // printf("t_av=%f,u_cs=%e，",t_av,u_cs); PT2RAMD(p,t_av,k,range); PT2PRN(p,t_av,pr,range); H=0.023*k*pow(re,0.8)*pow(pr,0.4)/(1000*di); //包壳表面换热系数h { double delta_theta1,delta_theta2; delta_theta1=qs*1.35*1.03*py[i-1]/H; //未沸腾下的温升 delta_theta2=t_s-t_out+25*pow(qs*1.35*1.03*py[i-1]/1000,0.25)*exp(-p/6.2); //沸腾下的温升 printf("delta_theta1=%f，delta_theta2=%f\n",delta_theta1,delta_theta2); //沸腾时温升的选用 if(delta_theta2>=delta_theta1) { delta_theta=delta_theta1; } else { delta_theta=delta_theta2; } t_cs=t_out+delta_theta; shu[2][i-1]=t_cs; name[2]="t_cs"; } } /***********************包壳内表面和燃料芯块外表面温度******************************/ { double kc,hg=5678,t_avci,t_ci1,cha_tci=5; t_ci=(t_cs+50); for(;cha_tci>0.1;) { t_avci=(t_cs+t_ci)/2; kc=(7.73+0.0315*t_avci-2.87*pow(t_avci,2)*pow(10,-5)+1.552*(t_avci,3)*pow(10,-8))/1000; t_ci1=t_cs+ql*1.35*1.03*py[i-1]*log(9.5/8.6)/(2*pi*kc); cha_tci=fabs(t_ci1-t_ci); t_ci=t_ci1; } t_us=t_ci+ql*1000*1.35*1.03*2*py[i-1]/(pi*(0.0086+0.00819)*hg); printf("包壳表面温度:%f，包壳内表面温度:%f，芯块外表温度:%f\n",t_cs,t_ci,t_us); shu[3][i-1]=t_ci; name[3]="t_ci"; shu[4][i-1]=t_us; name[4]="t_us"; } /*****************************燃料芯块中心温度***********************************/ { double t_u0=2600,qv_a,qv_a1,integ_b,integ_a,cha_qva=1000; qv_a=ql*1.35*1.03*py[i-1]*1000/(4*pi); integ_b=3824*log(402.55+t_us)+4.788*pow(10,-11)*0.25*pow(t_us+273.15,4); //ku对t的积分下限 for(;cha_qva>0.001;) { t_u0-=0.0001; integ_a=3824*log(402.55+t_u0)+4.788*pow(10,-11)*0.25*pow(t_u0+273.15,4); //ku对t的积分上限 qv_a1=integ_a-integ_b; cha_qva=fabs(qv_a1-qv_a); if(t_u0<1000)break; } printf("燃料芯块中心温度:%f，",t_u0); shu[5][i-1]=t_u0; name[5]="t_u0"; } /*******************************DNBR的计算*************************************/ { double h_fg,h_gs,h_ls,q_c; P2HG(p,h_gs,range); //饱和气的焓 P2HL(p,h_ls,range); //饱和水的焓 printf("h_gs=%f，h_ls=%f",h_gs,h_ls); h_fg=h_gs-h_ls; //汽化潜热 x_e=(h-h_ls)/h_fg; G=wb/ab; {//W-3公式 double oper1,oper2,oper3; oper1=3.1546*pow(10,6)*(2.022-6.24*pow(10,-2)*p+(0.1722-1.43*pow(10,-2)*p)*exp((18.177-5.99*pow(10,-1)*p)*x_e)); oper2=((0.1484-1.596*x_e+0.1729*x_e*fabs(x_e))*7.374*pow(10,-4)*G+1.037)*(1.157-0.869*x_e); oper3=(0.2644+0.8357*exp(-124.1*di))*(0.8258+3.41*pow(10,-7)*(1000*(h_ls-h_in))); q_c=oper1*oper2*oper3; } DNBR=q_c/(qs*1000*1.35*1.03*py[i-1]); printf("出口焓:%f，x:%f，q_c:%f，DNBR:%f\n",h,x_e,q_c,DNBR); shu[6][i-1]=q_c; name[6]="q_c"; shu[7][i-1]=DNBR; name[7]="DNBR"; } /***********************************水利计算************************************/ { double delta_p_f,delta_p_h,delta_p_a; double v_init,v_out_end; double g_w=9.8,G,loss_v1,loss_v2,l_v,high; high=z/3; G=wb/ab; PT2V(p,t_init,v_init,range); // printf("\nv_in=%f，",v_init); delta_p_in=0.75*pow(G,2)*v_init/2; //流道进口阻力 delta_p[0][0]=delta_p_in; if(i==3) { PT2V(p,t_out,v_out_end,range); delta_p_out=pow(G,2)*v_out_end/2; //流道出口阻力 // printf("\nv_=%f，G=%f\n",v_out_end,G); printf("di=%f，",di); delta_p_ga=1.05*pow(G,2)*(v_init+v_out_end)/2; //局部阻力，定位格架 delta_p[0][1]=delta_p_out; delta_p[0][2]=delta_p_ga; } PT2V(p,t_out,loss_v2,range); PT2V(p,t_in,loss_v1,range); delta_p_a=pow(G,2)*(loss_v2-loss_v1); //加速压降 delta_p_h=2*g_w*high/(loss_v2+loss_v1); //位能 { double u_w,u_f,f,moody_f; PT2ETA(p,(t_in+t_out)/2,u_f,range); PT2ETA(p,(t_cs+t_cs_up)/2,u_w,range); // printf("动力粘度=%f，"u_w); moody_f=moody(re,0.0000015/di); f=moody_f*pow((u_w/u_f),0.6); // printf("re=%f，moody_f=%f，f=%f，",re,moody_f,f); l_v=(loss_v2+loss_v1)/2; // printf("\nl_v=%f，di=%f，",l_v,di); delta_p_f=f*high*pow(G,2)*l_v/(2*di); //摩擦压降 } delta_p[i][0]=delta_p_a; delta_p[i][1]=delta_p_h; delta_p[i][2]=delta_p_f; // printf("\nfffff %f\n",delta_p[i][3]); printf("加速压降:%f，位能:%f\n摩擦压降:%f\n\n",delta_p_a,delta_p_h,delta_p_f); } /*************************************************************************************************************/ t_in=t_out; //主循环一次结束，更改下一控制体进口温度 t_cs_up=t_cs; //主循环一次结束，更改下一控制体进口壁面温度 } printf("进口阻力:%f，出口阻力:%f，格架阻力:%f\n\n",delta_p