控制工程师论坛

        为了更好的验证，我们决定把验证的代码公布，欢迎大家讨论和指正，发现问题可以直接在论坛发帖。

代码内容:

#include "iostream.h"
#include "stdio.h"
#include "math.h"
///////试验方式：C1——C，进口侧管道，在风管中测压////////
///////流量测试：90度弧形进口喷嘴/////////////////////////
//////////////换算条件/////////////
//大气压力：101325               //
//大气温度：20                   //
//大气湿度：50                   //
//换算转速：740                  //
///////////////////////////////////

////////测试数据：大气条件和工况点////////
/*   大气条件              |   流量测量   |  进口条件   |  出口条件  |    整机
 压力 温度 湿度         | 静压 压差    | 静压   温度 | 静压 温度 |转速  电机功率 噪声 电机效率 
 Pa   ℃    %               |   Pa      Pa   |  Pa     ℃     |  Pa   ℃      |r/min    W        dB(A)     %    
 99522 25.50 89.00| -377  377    | -205           |                     | 740   44306                0.93  

*/

void main()
{
 double pi = acos(-1.0);//
 double d  = 1.756;     //C式喉部直径
 double _k = 1.4;
 
 double Rw, pa, hu, ta; 

 pa = 99522;        //测试大气条件
 ta = 25.5;         //测试大气条件
 hu = 0.89;         //测试大气条件

 double pp3= -205;  //进口静压
 
 double mu= (17.1+0.048*ta)*0.000001;//气体粘性系数

 double psat, pv;

 //////////////////////////////计算湿空气气体常数////////////////////////
 if(ta<=30.0)
  psat = exp(17.438*ta/(239.78+ta)+6.4147);
 else if(ta>30.0 && ta<100.0)
  psat = 610.8 + 44.442*ta + 1.4133*ta*ta + 0.02768*pow(ta,3) +
         2.55667e-4*pow(ta,4) + 2.89166e-6*pow(ta,5);

 pv = hu*psat;

    Rw = 287.0/(1.0-0.378*pv/pa);
 //////////////////////////////计算湿空气气体常数////////////////////////
 

 double dp = 377;  //节流压差

 double pu = pa-dp;//当地静压值

 ///////////////////////////////////////////////////////////////////////

 double afa, e, rd, theta_u, ro_u, qm, Red,qm_2, error,aera;

 //////////////测试状态下的质量流量的计算///////////////////////////////
 aera= acos(-1.0)*d*d/4.0;   //喉部面积 
 afa = 0.98735088935932647;  //指定起始迭代的afa值,afa流量系数
 e   = 1-(0.55*dp)/pu;       //e系数
 rd  = 1-dp/pu;              //rd系数
 theta_u = ta+273.15;        //进口绝对温度  
 ro_u    = pu/(Rw*theta_u);  //密度计算

loop: qm      = (afa*e*pi*d*d*sqrt(2*ro_u*dp))/4;
  Red     = (4*qm)/(pi*mu*d);
  afa     = 1-0.004*sqrt(1000000/Red);
  qm_2    = (afa*e*pi*d*d*sqrt(2*ro_u*dp))/4;
  error   = qm-qm_2;
 if(fabs(error/qm_2)>0.001) goto loop;
 printf("%f\n",qm_2);
 qm = qm_2;
 //////////////测试状态下的质量流量的计算///////////////////////////////

 //double qv = qm/ro_u;

 double Cp;//定压热容
 Cp = Rw*_k*(_k-1.0);
 double aeredensity = pa/(Rw*(ta+273.15)); //大气密度计算

 double tsg1 = ta+(pu-pa)/(aeredensity*Cp);//进口滞止温度
 double p3 = pp3 + pa;                     //风机进口压力

 double d1p = 2.4;                         //风机流道直径
 double A3 = acos(-1.0)*d1p*d1p/4.0;
 
 double M = 0.5*qm*qm*(_k-1.0)*Rw*(tsg1+273.15)/(A3*A3*_k*p3*p3);
 double t1 = 2.0/(1.0+sqrt(1.0+4.0*M))*(tsg1+273.15)-273.15;//进口温度
 
 

 double Aqm = qm*qm/((acos(-1.0)*d1p*d1p/4.0)*(acos(-1.0)*d1p*d1p/4.0));
 double M2  = Aqm*(_k-1.0)*Rw*(tsg1+273.15)/(2.0*_k*p3*p3);
 double T   = (1.0+sqrt(1.+4.0*M2))/2.0;
 double Ma3 = sqrt(2.0*(T-1.0)/(_k-1.0));  //中间变量
 
 //Ma3 = 0.039072531167756552;

 double y3  = p3/(Rw*(t1+273.15));
 double Ve  = qm/y3/A3;
 double pd3 = 0.5*y3*Ve*Ve;               //中间变量           

 //pd3 = 106.13649099693293;

 double a    = Ma3*Ma3;
 double Fm3  = 1.0+a/4.0+a*a/40.0+a*a*a/1600.0;
 double Re3  = 4.0*qm/(acos(-1.0)*mu*d1p);
 double _LossCoefficient13 = -0.015-1.26*pow(Re3,-0.3);

 double psg1 = p3+pd3*Fm3*(1.0+_LossCoefficient13);//进口总压

 double ysg1 = psg1/(Rw*(tsg1+273.15));            //进口密度

 double qvsg1 = qm/ysg1;                           //进口体积流量

 ///////////////////////出口参数/////////////////////////////
 double p2   = pa;
 double f    = 41204.115;                          //轴功率
 //double f    = 41204.115000000005;
 //double Cp   = Rw*_k*(_k-1.0); 

 double tsg2 = tsg1 + f/qm/Cp;                     //出口总温
 double A2   = 4.08;//
 double M3   = 0.5*qm*qm*(_k-1.0)*Rw*(tsg2+273.15)/(A2*A2*_k*p2*p2);
 double t2   = 2.0/(1.0+sqrt(1.0+4.0*M3))*(tsg2+273.15)-273.15;      //出口温度

 ////////////////Ma2//////////////////
 double Aqm2 = qm*qm/A2/A2;         
 double M22  = Aqm2*(_k-1.0)*Rw*(tsg2+273.15)/(2.0*_k*p2*p2);
 double T2   = (1.0+sqrt(1.+4.0*M22))/2.0;
 double Ma2  = sqrt(2.0*(T2-1.0)/(_k-1.0));                          //出口马赫数

 ////////////////y2////////////////////
 double y2   = p2/(Rw*(t2+273.15));                                  //出口密度

 ////////////////pd2///////////////////
 double Ve2  = qm/y2/A2;
 double pd2  = 0.5*y2*Ve2*Ve2;                                       //出口动压

 /////////////////MachFactor////////////
 double a2    = (Ma2*Ma2);
 double Fm2   = 1.0+a2/4.0+a2*a2/40.0+a*a2*a2/1600.0;                //马赫数系数

 ////////////////psg2//////////////////
 double psg2   =  p2+(1.0+(Ma2*Ma2)/4.0+(Ma2*Ma2)*(Ma2*Ma2)/40.0+a2*a2*a2/1600.0)*pd2;//出口滞止压力
 double ysg2 = psg2/(Rw*(tsg2+273.15));                              //出口滞止密度
 double pf   = psg2 - psg1;                                          //通风机全压
 //////////////_fancompressibilitycorrectivecoefficient///////////
 double r    = 1.0+pf/psg1;
 double Pr   = qvsg1*pf;
 double ZK   = (_k-1.0)/_k*Pr/qvsg1/pf;
 double kp   = ZK*log10(r)/log10(1.0+ZK*(r-1.0));

 //////////////标态下的气体Rw/////////////////////////////////////
 double pao  = 101325.0;
 double tao  = 20.0;
 double huo  = 0.50;//湿度
 double psato;
 
 
 if(tao<=30.0)
  psato = exp(17.438*tao/(239.78+tao)+6.4147);
 else if(tao>30.0 && tao<100.0)
  psato = 610.8 + 44.442*tao + 1.4133*tao*tao + 0.02768*pow(tao,3) +
         2.55667e-4*pow(tao,4) + 2.89166e-6*pow(tao,5);

 double pvo = huo*psato;

    double Rwo=287.0/(1.0-0.378*pvo/pao);      //计算湿空气气体常数

 //////////////====kc====//////////////////////////////////////////
 double kcc     = Rw*(tsg1+273.15)/Rwo/(tao+273.15);
 double kc      = (1.0-kcc*(1-kp))/kp;

 double qvsg1o  = qvsg1*pow(kc,0.5);       //标准状态下的体积流量,转速相同
 //double qvsg1o2 = 61.308297042976967;

 double qvsg1operh  = qvsg1o*3600;
 //double qvsg1o2perh = qvsg1o2*3600;

 printf("%f\n",qvsg1operh);
 //printf("%f",qv);

 ///////////////////////////////标准全压计算////////////////////////
 double ysg1o = pao/(Rwo*(tao+273.15));
 double no    = 740;
 double ne    = 740;
 
 double pfo   = no*no/ne/ne*ysg1o/ysg1/kc*pf;

 ///////////////////////////////静压计算/////////////////////////////
 double pst   = p2 - psg1;

 ////计算静压系数
 r    = 1.0+pst/psg1;
 Pr   = qvsg1*pst;
 ZK   = (_k-1.0)/_k*Pr/qvsg1/pst;
 double kps   = ZK*log10(r)/log10(1.0+ZK*(r-1.0));

 ///////////////////////////////标准静压计算/////////////////////////
 double kss     = Rw*(tsg1+273.15)/Rwo/(tao+273.15);
 double ks      = (1.0-kss*(1-kps))/kps;

 double psto    = no*no/ne/ne*ysg1o/ysg1/ks*pst;

 ///////////////////////////////通风机的动压计算//////////////////////
 double VVe2    = qm/y2/A2;                   //////////////气流速度
 double pd      = 0.5*y2*VVe2*VVe2; 
 
 ////////////////////////////////动压系数计算/////////////////////////
 //
 double A1      = 3.141593;                       ////////风机进口面积
 double ysg     = psg1/(Rw*(tsg1+273.15));
 double M1      = qm*qm/(A1*A1*_k*psg1*ysg1);
 double Ma1     = sqrt(M1+1.217*M1*M1+1.369*M1*M1*M1+10.0*M1*M1*M1*M1);

 //通风机进口密度
 double y1      = ysg1*pow(1.0+0.5*(_k-1.0)*Ma1*Ma1, -1.0/(_k-1.0));
 //平均密度
 double ym      = (y1+y2)/2;
 //_DensityRatioInletToMean
 double kpy     = y1/ym;

 //计算动压系数kd
 double kds     = Rw*(tsg1+273.15)/Rwo/(tao+273.15);
 double kd      = (1.0-kds*(1-kpy))/kpy;
 //标准动压计算
 double pdo     = no*no/ne/ne*ysg1o/ysg1/kd*pd;

 ////////单位功计算//////////////////////////////////////////////////////
 double VVe1    = qm/y1/A1;                        //////////////气流速度
 double pd1     = 0.5*y1*VVe1*VVe1;               //////////////入口动压
 //马赫数系数计算
 double Ma1_a   = Ma1*Ma1;
 double Ma1_factor = 1.0+Ma1_a/4.0+Ma1_a*Ma1_a/40.0+Ma1_a*Ma1_a*Ma1_a/1600.0;

 ////单位功计算////////
 double p1      = psg1 - pd1*Ma1_factor;
 double w       = (p2-p1)/ym+pd2/y2-pd1/y1;

 ///////测试状态单位静功/////////////////////////////////////////////////////////
 double ws = w-pd2/y2;

 ////////空气功/////////
 double wu = qm*w;

 ////////空气静功///////
 double wus = qm*ws;
 
 ////////叶轮功/////////
 double wr  = f;
 double fo  = 0.0;
 double wa = f+fo;

 //全压内效率
 double er  = wu/wr;

 //静压内效率
 double esr = wus/wr;

 //轴效率
 double ea  = wu/wa;

 //静轴效率
 double eas = wus/wa;

 //////////////标准状态单位功//////////////////////////////////////////////////////
 double wo   =  pow(no/ne, 3)*ysg1o/ysg1*pow(kc, 0.5)*w;
 double wso  =  pow(no/ne, 3)*ysg1o/ysg1*pow(kc, 0.5)*ws;
 double wuo  =  pow(no/ne, 3)*ysg1o/ysg1*pow(kc, 0.5)*wu;
 double wuso =  pow(no/ne, 3)*ysg1o/ysg1*pow(kc, 0.5)*wus;
 double wro  =  pow(no/ne, 3)*ysg1o/ysg1*pow(kc, 0.5)*wr;
 double wao  =  pow(no/ne, 3)*ysg1o/ysg1*pow(kc, 0.5)*wa;

 double kpo  =  kc*kp;
 double kpso =  ks*kps;
 double kpyo =  kd*kpy;
}