燃烧室热力计算大作业

火箭发动机原理大作业

2015年5月19日

专题 燃烧室热力计算

作者 陈志远 学

号

12151175 院（系）名称

宇航学院 宇航推进系

目录

理论方法 (3)

公式推导 (4)

程序设计 (5)

原始数据 (14)

结果与分析 (15)

参考文献 (16)

理论方法

选择用平衡常数法计算平衡组分。

1.编写在给定压强、温度下计算平衡组分的函数Balance()

2.用这个函数计算出两个Tf下的两组平衡组分，然后计算出两个总焓，进行线性插值，确定Tf

3.由上一步求出的真正的Tf计算燃烧产物的平衡组分、熵值以及冻结比热比kf

4.由冻结比热比和Tf计算理论特征速度

5.由2、3步得到的结果，（用Balance函数）求出喷管出口界面在Te1、Te2条件下的平衡组分和熵值，因为是等熵流动过程，所以用熵进行“内插法”计算出Te。

6.将Te带入F函数中，求出出口的平衡组分，进而求出燃烧产物的总焓和发动机的理论比冲。

公式推导

计算平衡组分使用的公式取自《火箭发动机原理》”在给定压强和温度条件下计算平衡组分的控制方程组“一节的 (4-43) 到 (4-58) 以及 (4-43)’ 到(4-56)’.

总焓的计算使用公式

I

?m =∑I j n j N

j=1

I j 的数值由查表而来，并非来自拟合公式。由于表格数据是离散的，因此在计算初期进行了试算来确定绝热燃烧温度在2400K 和2500K 之间.

将2400K 和2500K 时的总焓求出，并用以下公式进行插值，计算出真正的绝热燃烧温度.

T f =T f1+I p ??I m1

?I m2??I m1

?(T f2?T f1)

求熵的时候使用以下公式

S

?=∑S j n j +∑S j n j N

j=L+1

L

j=1

其中S j 使用课本公式(4-162)进行计算.

理论特征速度 c 理?

=√

m

喷管中，利用Balance() 函数算出T e1和T e2时喷管出口的平衡组分，进而得到熵，然后用以下内插公式得到真正的出口温度T e 和总焓I

?e . T e =T e1+S

?0c ?S ?e1S

?e2?S ?e1(T e2?T e1)

I ?e =I ?e1+T e ?T e1

T e2?T e1

(I

?e2?I ?e1) 发动机的理论比冲为

u e =√2(I

?c ?I ?e )

程序设计

#include

#include

#include

struct mole

{

double AL2O3;

double CO2;

double H2O;

double CO;

double H2;

double N2;

double HCL;

double SO2;

double O2;

double OH;

double NO;

double H;

double CL;

}; //定义各组分摩尔数保存在结构体内

struct mole N[10]={0.0}; //N结构体数组存储不同燃烧温度下的平衡组分

int Pc=70,Pe=1;

//计算给定压强、温度下的平衡组分

void balance(struct mole *n_ph,int T,int P)

{

struct mole n,ni; //n构体存储各次迭代结果

double Nh_1,No_1,A,B,ng;

int i=0;

double K[6]={0.0};

switch (T) //根据温度选择平衡常数

{

case 1000: K[0]=0.7254;K[1]=0.8728E-10;K[2]=0.3604E-11;K[3]=0.5148E-17;K[4]=4.7852E-18;K[5]=0.7302E-8; break;

case 1100: K[0]=1.0560;K[1]=0.1314E-8;K[2]=0.8519E-10;K[3]=0.6676E-15;K[4]=5.9776E-16;K[5]=0.5277E-7; break;

case 2000: K[0]=4.7410;K[1]=0.2892E-3;K[2]=0.1378E-3;K[3]=0.2631E-5;K[4]=1.9697E-6;K[5]=0.3926E-3; break;

case 2400: K[0]=6.0390;K[1]=0.3634E-2;K[2]=0.2573E-2;K[3]=0.2516E-3;K[4]=1.8055E-4;K[5]=0.2410E-2; break;

case 2500: K[0]=6.3110;K[1]=0.6037E-2;K[2]=0.4625E-2;K[3]=0.6284E-3;K[4]=4.4720E-4;K[5]=0.3391E-2; break;

case 2800: K[0]=7.0080;K[1]=0.2233E-1;K[2]=0.2091E-1;K[3]=0.6649E-2;K[4]=4.5983E-3;K[5]=0.8786E-2; break;

case 2900: K[0]=7.2020;K[1]=0.3256E-1;K[2]=0.3228E-1;K[3]=0.1312E-1;K[4]=9.0698E-3;K[5]=0.1148E-1; break;

case 3000: K[0]=7.3820;K[1]=0.4628E-1;K[2]=0.4841E-1;K[3]=0.2475E-1;K[4]=1.7070E-2;K[5]=0.1472E-1; break;

default: printf("没有所需温度的数据from function Balance()\n");

}

n.AL2O3 = 0.9266;

n.SO2 = 2.4177;

ni.AL2O3 = 0.9266;

ni.SO2 = 2.4177;

//AL2O3和SO2精确解

//初次迭代开始

n.O2 = 0;

n.OH = 0;

n.H = 0;

n.CL = 0;

n.NO = 0;

n.N2 = 0.5*(5.6773-n.NO);

n.HCL = 5.7692-n.CL;

Nh_1 = 40.8672-n.HCL-n.OH-n.H;

No_1 = 25.8632-2*n.SO2-3*n.AL2O3-2*n.O2-n.OH-n.NO;

A = No_1-11.2738;

B = 0.5*Nh_1-A;

n.CO2 = 1/(2*(K[0]-1)) * (-(K[0]*B+No_1)+sqrt(pow(K[0]*B+No_1,2)+4*(K[0]-1)*A*11.2738));

n.CO = 11.2738-n.CO2;

n.H2O = A-n.CO2;

n.H2 = B+n.CO2;

ng = n.CL+n.CO+n.CO2+n.H+n.H2+n.H2O+n.HCL+n.N2+n.NO+n.O2+n.OH+n.SO2;

//初次迭代所有平衡组分已全部解出

for (i=1;i<10;i++)

{

ni.O2 = pow((n.H2O/n.H2),2)*K[1]*K[1]*pow(P/ng,-1);

ni.OH = n.H2O*pow(n.H2,-0.5)*K[2]*pow(P/ng,-0.5);

ni.H = pow(n.H2*K[3],0.5)*pow(P/ng,-0.5);

ni.CL = n.HCL/ni.H*K[4]*pow(P/ng,-1);

ni.NO = pow(n.N2*ni.O2*K[5],0.5);

ni.N2 = 0.5*(5.6773-ni.NO);

ni.HCL = 5.7692-ni.CL;

Nh_1 = 40.8672-ni.HCL-ni.OH-ni.H;

No_1 = 25.8632-2*ni.SO2-3*ni.AL2O3-2*ni.O2-ni.OH-ni.NO;

A = No_1-11.2738;

B = 0.5*Nh_1-A;

ni.CO2 = 1/(2*(K[0]-1)) * (-(K[0]*B+No_1)+sqrt(pow(K[0]*B+No_1,2)+4*(K[0]-1)*A*11.2738));

ni.CO = 11.2738-ni.CO2;

ni.H2O = A-ni.CO2;

ni.H2 = B+ni.CO2;

ng = ni.CL+ni.CO+ni.CO2+ni.H+ni.H2+ni.H2O+ni.HCL+ni.N2+ni.NO+ni.O2+ni.OH+ni .SO2;

//至此，本次迭代所有平衡组分全部解出

if (i>1)

{

if( (fabs(ni.AL2O3-n.AL2O3)/n.AL2O3 <= 0.0001) && (fabs(ni.CL-n.CL)/n.CL <= 0.0001) && (fabs(ni.CO-n.CO)/n.CO <= 0.0001) && (fabs(ni.CO2-n.CO2)/n.CO2 <= 0.0001) && (fabs(ni.H-n.H)/n.H <= 0.0001) && (fabs(ni.H2-n.H2)/n.H2 <= 0.0001) && (fabs(ni.H2O-n.H2O)/n.H2O <= 0.0001) && (fabs(ni.HCL-n.HCL)/n.HCL <= 0.0001) && (fabs(ni.N2-n.N2)/n.N2 <= 0.0001) && (fabs(ni.NO-n.NO)/n.NO <= 0.0001) && (fabs(ni.O2-n.O2)/n.O2 <= 0.0001) && (fabs(ni.OH-n.OH)/n.OH <= 0.0001) && (fabs(ni.SO2-n.SO2)/n.SO2 <= 0.0001) )

break;

}

n.AL2O3 = ni.AL2O3; n.CL=ni.CL; n.CO=ni.CO; n.CO2=ni.CO2; n.H=ni.H; n.H2=ni.H2; n.H2O=ni.H2O; n.HCL=ni.HCL; n.N2=ni.N2; n.NO=ni.NO; n.O2=ni.O2; n.OH=ni.OH; n.SO2=ni.SO2;

}

n_ph->AL2O3 = ni.AL2O3;

n_ph->CL = ni.CL;

n_ph->CO = ni.CO;

n_ph->CO2 = ni.CO2;

n_ph->H = ni.H;

n_ph->H2 = ni.H2;

n_ph->H2O = ni.H2O;

n_ph->HCL = ni.HCL;

n_ph->N2 = ni.N2;

n_ph->NO = ni.NO;

n_ph->O2 = ni.O2;

n_ph->OH = ni.OH;

n_ph->SO2 = ni.SO2; //输出平衡组分

}

//计算总焓

double I(struct mole *n,int T)

{

double I,i0[13]={0}; //i0存储的是表“某些燃烧产物的总焓”

switch (T)

{

case 1000: i0[0]=-1592666E-3; i0[1]=137347E-3; i0[2]=-88757E-3; i0[3]=-360079E-3; i0[4]=232705E-3; i0[5]=20830E-3; i0[6]=-215594E-3; i0[7]=-71140E-3; i0[8]=21620E-3; i0[9]=112785E-3; i0[10]=22861E-3; i0[11]=63170E-3; i0[12]=-262609E-3; break;

case 1100: i0[0]=-1580163E-3; i0[1]=139564E-3; i0[2]=-85410E-3; i0[3]=-354578E-3; i0[4]=234784E-3; i0[5]=23871E-3; i0[6]=-211407E-3; i0[7]=-67950E-3; i0[8]=24919E-3; i0[9]=116211E-3; i0[10]=26372E-3; i0[11]=66262E-3; i0[12]=-257122E-3; break;

case 2000: i0[0]=-1460414E-3; i0[1]=159052E-3; i0[2]=-53686E-3; i0[3]=-301847E-3; i0[4]=253500E-3; i0[5]=53120E-3; i0[6]=-169276E-3; i0[7]=-37255E-3; i0[8]=56319E-3; i0[9]=148462E-3; i0[10]=59367E-3; i0[11]=96048E-3; i0[12]=-205766E-3; break;

case 2400: i0[0]=-1378116E-3; i0[1]=167561E-3; i0[2]=-39072E-3; i0[3]=-277336E-3; i0[4]=261818E-3; i0[5]=67103E-3; i0[6]=-148768E-3; i0[7]=-22869E-3; i0[8]=70839E-3; i0[9]=163239E-3; i0[10]=74668E-3; i0[11]=110180E-3; i0[12]=-182337E-3; break;

case 2500: i0[0]=-1363026E-3; i0[1]=169680E-3; i0[2]=-35391E-3; i0[3]=-271145E-3; i0[4]=263898E-3; i0[5]=70674E-3; i0[6]=-143532E-3; i0[7]=-19225E-3; i0[8]=74479E-3; i0[9]=166959E-3; i0[10]=78550E-3; i0[11]=113774E-3; i0[12]=-176445E-3; break;

case 2800: i0[0]=-1317254E-3; i0[1]=176022E-3; i0[2]=-24294E-3; i0[3]=-252446E-3; i0[4]=270136E-3; i0[5]=81546E-3; i0[6]=-127617E-3; i0[7]=-8201E-3; i0[8]=85542E-3; i0[9]=178163E-3; i0[10]=90328E-3; i0[11]=124682E-3; i0[12]=-158694E-3; break;

case 2900: i0[0]=-1301829E-3; i0[1]=178131E-3; i0[2]=-20580E-3; i0[3]=-246174E-3; i0[4]=272216E-3; i0[5]=85219E-3; i0[6]=-122251E-3; i0[7]=-4499E-3; i0[8]=89240E-3; i0[9]=181913E-3; i0[10]=94296E-3; i0[11]=128355E-3; i0[12]=-152756E-3; break;

case 3000: i0[0]=-1286320E-3; i0[1]=180239E-3; i0[2]=-16859E-3; i0[3]=-239886E-3; i0[4]=274296E-3; i0[5]=88916E-3; i0[6]=-116858E-3; i0[7]= -786E-3; i0[8]=92946E-3; i0[9]=185669E-3; i0[10]=98284E-3; i0[11]=132046E-3; i0[12]=-146807E-3; break;

default : printf("没有所需温度的数据from function I()\n"); break;

}

I = n->AL2O3*i0[0] + n->CL*i0[1] + n->CO*i0[2] + n->CO2*i0[3] +

n->H*i0[4] + n->H2*i0[5] + n->H2O*i0[6] + n->HCL*i0[7] + n->N2*i0[8] + n->NO*i0[9] + n->O2*i0[10] + n->OH*i0[11] + n->SO2*i0[12];

return I;

}

double Tf_get(int Tf1,int Tf2,int Pc) //内插法求Tf

{

double I1=0,I2=0,temp;

struct mole nn[2]={0.0};

balance(&nn[0],Tf1,Pc);

balance(&nn[1],Tf2,Pc);

I1=I(&nn[0],Tf1);

I2=I(&nn[1],Tf2);

temp = Tf1 + (-2400-I1)/(I2-I1)*(Tf2-Tf1);

return temp;

}

double S_get(double Tf,double p,struct mole *N)

{

double temp,R0=8.3144;

double para[13][6] =

{ {1.7612190E1,0,0,0,0,-9.5677512E1},

{2.9657881E0,-4.3004990E-4,1.6643053E-7,-2.9717571E-

11,2.0047255E-15,3.0001737E0},

{2.9511519E0,1.5525567E-3,-6.1911411E-7,1.1350336E-10,-

7.7882732E-15,6.5314450E0},

{4.4129266E0,3.1922896E-3,-1.2978230E-6,2.4147446E-10,-

1.6742986E-14,-7.2875769E-1},

{2.5000000E0,0,0,0,0,-4.6001096E-1},

{3.0436897E0,6.1187110E-4,-7.3993551E-9,-2.0331907E-

11,2.4593791E-15,-1.6481339E0},

{2.6707532E0,3.0317115E-3,-8.5351570E-7,1.1790853E-10,-

6.1973568E-15,6.8838391E0},

{2.7356497E0,1.4979644E-3,-5.0699071E-7,8.2790246E-11,-

5.1747522E-15,

6.6280129E0},

{2.8545761E0,1.5976316E-3,-6.2566254E-7,1.1315849E-10,-

7.6897070E-15,6.3902879E0},

{3.1529360E0,1.4059955E-3,-5.7078462E-7,1.0628209E-10,-

7.3720783E-15,6.9446465E0},

{0.3621950E1,0.73618260E-3,-0.19652230E-6,0.36201560E-10,-

0.28945630E-14,0.36150960E1},

{2.8885544E0,9.9835061E-4,-2.1878804E-7,1.9802785E-11,-

3.8452940E-16,5.5597016E0},

{0.52451360E1,0.19704200E-2,-0.80375770E-6,0.15149970E-9,-

0.10558000E-13,-0.10873520E1},

};

double Sj[13]={0};double P[13]={0};

double ng;

ng = N->CL+N->CO+N->CO2+N->H+N->H2+N->H2O+N->HCL+N->N2+N->N O+N->O2+N->OH+N->SO2;

P[1] = p * N->CL/ng;

P[2] = p * N->CO/ng;

P[3] = p * N->CO2/ng;

P[4] = p * N->H/ng;

P[5] = p * N->H2/ng;

P[6] = p * N->H2O/ng;

P[7] = p * N->HCL/ng;

P[8] = p * N->N2/ng;

P[9] = p * N->NO/ng;

P[10] = p * N->O2/ng;

P[11] = p * N->OH/ng;

P[12] = p * N->SO2/ng;

for (int i=1;i<13;i++)

Sj[i] = R0*(para[i][0]*log(Tf)+para[i][1]*Tf+0.5*para[i][2]*Tf*Tf+para[i][3]*pow(Tf,3)/3.0+ para[i][4]*pow(Tf,4)/4.0+para[i][5]) - R0*log(P[i]);

Sj[0]=R0*(para[0][0]*log(Tf)+para[0][1]*Tf+0.5*para[0][2]*Tf*Tf+para[0][3 ]*pow(Tf,3)/3.0+para[0][4]*pow(Tf,4)/4.0+para[0][5]); //AL2O3的Sj

temp = Sj[0]*N->AL2O3 + Sj[1]*N->CL + Sj[2]*N->CO + Sj[3]*N->CO2 + Sj[4]*N->H + Sj[5]*N->H2 + Sj[6]*N->H2O + Sj[7]*N->HCL + Sj[8]*N->N2 + Sj[9]*N->NO + Sj[10]*N->O2 + Sj[11]*N->OH + Sj[12]*N->SO2;

return temp;

}

int main()

{

double Tf,Te,Ic,Ie,m,gama,c_star,Sc,ypxl,ng_e,m_e,u_e,Is; int Te1,Te2;

double R0=8314.4;

//平衡组分计算

balance(&N[0],2400,Pc);

Ic=I(&N[0],2400);

printf("2400℃时的平衡组分：\n");

printf("AL2O3 : %lf\nCL : %lf\nCO : %lf\nCO2 : %lf\nH : %lf\nH2 : %lf\nH2O : %lf\nHCL : %lf\nN2 : %lf\nNO : %lf\nO2 : %lf\nOH : %lf\nSO2 : %lf\n",N[0].AL2O3,N[0 ].CL,N[0].CO,N[0].CO2,N[0].H,N[0].H2,N[0].H2O,N[0].HCL,N[0].N2,N[0].NO,N[0].O 2,N[0].OH,N[0].SO2);

printf("总焓I = %lf\n",Ic);

Tf=Tf_get(2400,2500,Pc);

printf("绝热燃烧温度Tf = %lf\n",Tf);

m = 1000.0/(N[0].AL2O3+N[0].CL+N[0].CO+N[0].CO2+N[0].H+N[0].H2+N[0].H2O+N [0].HCL+N[0].N2+N[0].NO+N[0].O2+N[0].OH+N[0].SO2);

//使用2400℃时的定压比热来计算冻结比热比

double Cp_f = (N[0].AL2O3*35.59 + N[0].CL*5.0654 + N[0].CO*8.784 + N[0].CO2*14.763 + N[0].H*4.9680 + N[0].H2*8.4977 + N[0].H2O*12.463 + N[0].HCL*8.685 + N[0].N2*8.733 + N[0].NO*8.874 + N[0].O2*9.248 + N[0].OH*8.561 + N[0].SO2*14.063)*4.1858518; //单位：J/mol.K

double Cv_f = Cp_f - R0/m;

double kf = Cp_f/Cv_f;

//特征速度计算

gama = sqrt(kf)*pow((2/(kf+1)),(kf+1)/(2*kf-2));

c_star = sqrt(R0*Tf/m)/gama;

//喷管计算

Sc = S_get(Tf,Pc,&N[0]);

Te = Tf * pow(Pe*1.0/Pc,(kf-1)/kf);

Te1=1000; Te2=1100;

balance(&N[1],Te1,Pe);

balance(&N[2],Te2,Pe);

double S_Te1 = S_get(Te1,Pe,&N[1]);

double S_Te2 = S_get(Te2,Pe,&N[2]);

Te = Te1 + (Sc-S_Te1)/(S_Te2-S_Te1)*(Te2-Te1); //内插法求Te

double I_Te1 = I(&N[1],Te1);

double I_Te2 = I(&N[2],Te2);

Ie = I_Te1 + (Te-Te1)/(Te2-Te1)*(I_Te2-I_Te1); //内插法求Ie

u_e = sqrt(2*(Ic-Ie)*1000);

Is = u_e; //设计状态

//输出命令组

printf("喷管出口的平衡组分：\n");

printf("AL2O3 : %lf\nCL : %lf\nCO : %lf\nCO2 : %lf\nH : %lf\nH2 : %lf\nH2O : %lf\nHCL : %lf\nN2 : %lf\nNO : %lf\nO2 : %lf\nOH : %lf\nSO2 : %lf\n",N[2].AL2O3,N[2 ].CL,N[2].CO,N[2].CO2,N[2].H,N[2].H2,N[2].H2O,N[2].HCL,N[2].N2,N[2].NO,N[2].O 2,N[2].OH,N[2].SO2);

printf("理论特征速度c* = %lf\n",c_star);

printf("燃烧室的熵Sc = %lf\n",Sc);

printf("冻结比热比kf = %lf\n",kf);

printf("喷管出口截面燃烧产物的温度Te = %lf K\n",Te);

printf("喷管出口截面温度Te = %lf\n",Te);

printf("燃烧产物的总焓Ie = %lf\n",Ie);

printf("发动机的理论比冲Is= %lf\n",Is);

system("pause");

return 0;

}

原始数据

假定化学式C11.2738H40.8672O25.8632N5.6773Cl5.7692S2.4177Al1.8532

推进剂总焓

p 2400kJ/kg

I=-

燃烧室压强p c=70.924×105Pa 喷管出口压强p e=1atm

结果与分析

程序运行结果：

程序给出的是2400℃时的平衡组分，原因找不到2412℃（即得到的Tf）时的平衡常数的数据，因此给出了较为接近的2400℃的平衡组分，温度误差不到1%，精确度应该能得到保证。

通过运行结果可以看出，燃烧产物中容易反应的组分基本已经消耗殆尽，燃烧产物的总焓下降非常大，用来对外做功，推进火箭。比冲为2150.83 m/s，并不是十分优秀。

参考文献

[1]李宜敏,张中钦,张远君. 固体火箭发动机原理[M],2010.

[2]王新月. 气体动力学基础[M],2006.

[3]NASA. SP-273, 1972.