1 #include "pfr.hpp"
2 #include "RungeKutta.cpp"
3 using namespace std;
4
pfr(stream * s1,stream * s2,double ** t,int nb_r,reaction ** rr,double u,double ta)5 pfr::pfr ( stream * s1 , stream * s2 , double ** t , int nb_r , reaction ** rr , double u , double ta ) {
6
7 F = s2;
8 F->m=0;
9 P=s1->P;
10 for ( i = 0 ; i < s1->nb ; i++ ) {
11 F->chem[i]->m = s1->chem[i]->m;
12 F->m+=F->chem[i]->m;
13 }
14 F->set(s1->P, s1->T);
15 m_in=F->m;
16 a = t;
17 rx = rr;
18 n=nb_r;
19 m= F->nb;
20 U=u;
21 Ta=ta;
22 T = F->T;
23 C = new double[m];
24 y = new double[m+1];
25 r=new double[n];
26 OK=true;
27 explode=true;
28 solver = new RungeKutta<pfr>(m+1);
29 }
30
~pfr()31 pfr::~pfr() {
32 delete [] r;
33 delete [] C;
34 delete [] y;
35 delete solver;
36 }
37
run()38 bool pfr::run() {
39
40 for ( i = 0 ; i < m ; i++ )
41 y[i]=F->chem[i]->n();
42
43 y[m]=T;
44
45 solver->set ( this , y , 0.0 , L );
46
47 dL=solver->dx();
48
49 OK=solver->run();
50
51 sum = F->m;
52 F->m = 0;
53
54 for ( i = 0 ; i < m ; i++ )
55 F->m+=F->chem[i]->m;
56 for ( i = 0 ; i < m ; i++ )
57 F->chem[i]->m *= sum/F->m;
58 // if (OK)
59 //mass balance!
60 if ( fabs(m_in-F->m) > EPS || !explode )
61 OK=false;
62
63 return OK;
64 }
65
f(int eq,double l,double * y)66 double pfr::f ( int eq , double l , double * y ) {
67
68
69 sum=F->m;
70 F->m=0;
71 for(i=0; i<m;i++)
72 {
73 if(y[i]<0) y[i]=0;
74 F->chem[i]->m = y[i]*F->chem[i]->M/1000.0;
75 F->m+=F->chem[i]->m;
76 }
77
78
79 for(i=0; i<m;i++)
80 F->chem[i]->m *= sum/F->m;
81
82
83
84 F->m=sum;
85 T=y[m];
86 if(T>MAX_TEMP)
87 {
88 cout << "ERROR 11\n\n";
89 exit(0);
90 }
91
92
93 for(i=0; i<m;i++)
94 C[i]=F->chem[i]->n()/F->v;
95
96 for(j=0;j<n;j++)
97 r[j] = rx[j]->rate(T,C);
98
99 if(0<=eq && eq<m) //return dFi/dL
100 {
101 tmp=0.0;
102 for(j=0;j<n;j++) tmp+=a[eq][j]*r[j];
103 tmp *= (pi*D*D/4.0);
104 }
105
106
107
108 if(eq==m) //return dT/dL
109 {
110
111
112 F->set(F->P,T);
113
114 tmp=0.0;
115 for(j=0;j<n;j++)
116 tmp -= r[j]*rx[j]->dHr(T);
117
118
119 tmp *= (pi*D*D/4.0);
120
121
122 tmp += (pi*D)*U*(Ta-T);
123
124
125 tmp1=0.0;
126 for(i=0;i<m;i++) {
127 tmp1+= y[i]*F->chem[i]->Cp()*0.001;
128 }
129 tmp /= tmp1;
130
131 if(fabs(tmp*dL)>500.0)
132 {
133 cout << "ERROR 13\n\n";
134 exit(0);
135 }
136 }
137
138 return tmp;
139 }
140
141
142
get_cost(void)143 double pfr::get_cost ( void ) {
144 dL=L*pi*pow(D,2)/4.0;
145 if(dL<0.3) dL=0.3; if(dL>520) dL=520;
146 sum = 3.4974+0.4485*log10(dL)+0.1074*pow(log10(dL),2);
147 sum = pow(10, sum);
148 P= (P-1)*101.325/100;
149 dL=(P+1)*D/(317.46*(850-0.6*(P+1)))+0.0315;
150 sum *=(2.25+ 1.82*dL*4.2);
151 sum = sum*MS_YEAR/MS_2001;
152 return sum;
153 }
154
get_water()155 double pfr::get_water() {
156 sum = (U>EPS && T>Ta) ? U*L*pi*pow(D,2)/4*(T-Ta)/4.185/25.0 : 0.0;
157 return sum;
158 }
159
cost()160 void pfr::cost() {
161 cout << "WRITE FILE " << RUNTIME << name << ".cost" << " :\n\tBEGIN\n";
162 cout << "\t>>" << get_cost();
163 cout << "\n\tEND\n\n";
164 }
165
water()166 void pfr::water() {
167 cout << "WRITE FILE " << RUNTIME << name << ".water" << " :\n\tBEGIN\n";
168 if (U>EPS && T>Ta) sum = (U*L*pi*pow(D,2)/4*(T-Ta)/4.185/25.0);
169 else sum = 0.0;
170 cout << "\t>>" << sum;
171 cout << "\n\tEND\n\n";
172 }
173