-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathSeparator.c
273 lines (248 loc) · 9.03 KB
/
Separator.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
//this is a file in the VAc dynamic model software package
//this is the main file for the separator
//parent-routine: VAmodel.c
//subroutines involve: beta_separator.c enthalpy_7.c gamma_wilson_7.c
//Copyright: Rong Chen and Kedar Dave, June 2002
//VERSION 1.0
//if a user wants to build a DLL for the separator only, follow the instruction below:
//1. Add this line to the head of this file: #include "mex.h"
//2. Enable mexFunction(), which is the second routine included in this file (currently disabled)
//3. Use the following line to generate a DLL file: separator.dll
// mex separator.c beta_separator.c enthalpy_7.c gamma_wilson_7.c
//4. The name of the generated function is "separator()". This function should be called in MATLAB.
//The separator is modeled as a partial condenser
//A static flash calculation is first carried out and then the two outlet streams are fed into the gas header and the liquid buffer tank respectively
//The most tricky thing in this model is that we assume the flash temperature is always 5 degC higher than the cooling jacket temperature
#include <stdio.h>
#include <math.h>
#include "debug.h"
//#include "mex.h"
void Separator(double *dstatedt, double states[], double Feed, double Feed_comp[], double Feed_T, double F_SepLiquidOUT,
double T_SepCooler, double F_SepGasOUT, double Working_Level_Volume, double Total_Gas_Loop_Volume, double UA,
double VIJ[][7], double aij[][7], double MW[], double SpG[],double HVAP[], double HCAPLa[], double HCAPLb[],
double HCAPVa[], double HCAPVb[], double A[], double B[], double C[], double R, double Dw, double T_ref)
{
//--------------------------------------------------------------------------
//physical properties: VIJ, aij, MW, SpG, HVAP, HCAPLa, HCAPLb, HCAPVa, HCAPVb, A, B, C, R, Dw, T_ref
//--------------------------------------------------------------------------
/*output:
dstatedt: state derivatives
F_SepLiquidOUT_out: separator liqiud outlet flowrate (duplicated when no perfect control is selected)
--------------------------------------------------------------------------*/
/*input:
states: states
Feed: reactor feed flowrate
Feed_comp: reactor feed composition
Feed_T: reactor feed temperature
F_SepLiquidOUT: liqiud outlet flowrate
T_SepCooler: cooling jacket temperature
F_SepGasOUT: gas outlet flowrate
------------------------------------------------------------------------*/
/*Equipment Specification
Working_Level_Volume: liquid volume
Total_Gas_Loop_Volume: gas header volume
UA: Inverse of the Total Thermal Resistance
------------------------------------------------------------------------*/
//local variables
//states
double xL[1][7], Volume, TL[1], yV[1][7], P[1], TV[1];
//intermediate variables
double liquid_concentration, liquid_density, Y_Feed_comp[1][7];
double T[1]; //flash temperature
double ps[7]; //partial pressure
double Q, hl[1], hv[1], hl_feed[1], hv_feed[1];
double VRT, error[1], Feed_V, Feed_L, Feed_xL[7], Feed_xL_i[7], Feed_yV[1][7], beta[1]; //flash calculation variables
//derivatives
double dxdt[7], dydt[7], dMdt, dLdt, dTLdt, dPdt, dTVdt;
//others
int i, j;
double sum1, sum2;
double dummy[1], temp_HL[1], temp_HV[1];
#ifdef DEBUG
debug("Separator entered.\n");
#endif
/*get states*/
sum1=0.;
for (i=0;i<3;i++)
{
xL[0][i]=states[i];
xL[0][i+4]=states[i+3];
sum1+=xL[0][i] + xL[0][i+4];
};
xL[0][3]=1-sum1;
sum1=0.;
sum2=0.;
for (i=0;i<7;i++)
{
sum1+=xL[0][i]*MW[i]/SpG[i];
sum2+=xL[0][i]*MW[i];
};
liquid_concentration=Dw/sum1;
liquid_density=sum2*liquid_concentration;
Volume=states[6]*Working_Level_Volume; /*m3*/
TL[0]=states[7];
enthalpy_7(temp_HL, dummy, xL, xL, TL, MW, HVAP, HCAPLa, HCAPLb, HCAPVa, HCAPVb);
hl[0] = temp_HL[0];
sum1=0.;
for (i=0;i<3;i++)
{
yV[0][i]=states[i+8];
yV[0][i+4]=states[i+11];
sum1+=yV[0][i] + yV[0][i+4];
};
yV[0][3]=1-sum1;
P[0]=states[14];
TV[0]=states[15];
enthalpy_7(dummy, temp_HV, yV, yV, TV, MW, HVAP, HCAPLa, HCAPLb, HCAPVa,HCAPVb);
hv[0] = temp_HV[0];
for (j=0;j<7;j++)
Y_Feed_comp[0][j] = Feed_comp[j];
enthalpy_7(dummy, temp_HV, Y_Feed_comp, Y_Feed_comp, &Feed_T, MW, HVAP, HCAPLa, HCAPLb, HCAPVa, HCAPVb);
hl_feed[0] = temp_HV[0];
/*-----------------------------------------------------------------------------
/*Do a flash calculation under known P and T*/
T[0]=T_SepCooler+5;
error[0]=1.0;
for (i=0;i<7;i++)
{
Feed_xL[i]=xL[0][i];
ps[i]=exp(A[i]+B[i]/(T[0]+C[i]));
};
beta[0]=F_SepGasOUT/Feed;
while (fabs(error[0])>1e-10)
{
for (j=0;j<7;j++)
Feed_xL_i[j] = Feed_xL[j];
Beta_Separator(error, Feed_xL, beta, ps, P, T, Feed_xL_i, Feed_comp, VIJ, aij, R, T_ref);
beta[0]=beta[0]-0.1*error[0];
}
Feed_V=Feed*beta[0];
Feed_L=Feed-Feed_V;
for (i=0;i<7;i++)
Feed_yV[0][i]=(Feed*Feed_comp[i]-Feed_L*Feed_xL[i])/Feed_V;
enthalpy_7(dummy, temp_HV, Feed_yV,Feed_yV, T, MW, HVAP, HCAPLa, HCAPLb, HCAPVa,HCAPVb);
hv_feed[0] = temp_HV[0];
/*-----------------------------------------------------------------------------
/*liquid dynamics*/
sum1=0.;
sum2=0.;
for (i=0;i<7;i++)
{
sum1+= Feed_xL[i]*MW[i];
sum2+= xL[0][i]*MW[i];
};
dMdt=Feed_L-F_SepLiquidOUT;
dLdt=dMdt/liquid_concentration/Working_Level_Volume; /*volume percentage*/
/*composition update*/
sum1=0.;
for (i=0;i<7;i++)
{
dxdt[i]=(Feed_L*Feed_xL[i]-F_SepLiquidOUT*xL[0][i] - xL[0][i]*dMdt)/liquid_concentration/Volume;
sum1+=xL[0][i]*MW[i]*(HCAPLa[i]+HCAPLb[i]*TL[0]);
};
Q=UA*(T[0]-T_SepCooler);
dTLdt=(Feed*(hl_feed[0]-hl[0])-Feed_V*(hv_feed[0]-hl[0]) - Q)/(liquid_concentration*Volume)/sum1;
/*-----------------------------------------------------------------------------
/*vapor dynamics*/
VRT=Total_Gas_Loop_Volume/8.314/(TV[0]+273.15);
dPdt=(Feed_V-F_SepGasOUT)/VRT*14.696/101.325;
sum1=0.;
for (i=0;i<7;i++)
{
dydt[i]=(Feed_V*Feed_yV[0][i]-F_SepGasOUT*yV[0][i]-VRT*yV[0][i]*dPdt)/VRT/(P[0]/14.696*101.325);
sum1+=yV[0][i]*MW[i]*(HCAPLa[i]+HCAPLb[i]*TV[0]);
};
dTVdt=Feed_V*(hv_feed[0]-hv[0])/(VRT*P[0])/sum1;
/*-----------------------------------------------------------------------------*/
//derivatives
dstatedt[0] = dxdt[0];
dstatedt[1] = dxdt[1];
dstatedt[2] = dxdt[2];
dstatedt[3] = dxdt[4];
dstatedt[4] = dxdt[5];
dstatedt[5] = dxdt[6];
dstatedt[6] = dLdt;
dstatedt[7] = dTLdt;
dstatedt[8] = dydt[0];
dstatedt[9] = dydt[1];
dstatedt[10] = dydt[2];
dstatedt[11] = dydt[4];
dstatedt[12] = dydt[5];
dstatedt[13] = dydt[6];
dstatedt[14] = dPdt;
dstatedt[15] = dTVdt;
#ifdef DEBUG
debug("Separator left.\n");
#endif
}
/*
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *states;
double Feed;
double *Feed_comp;
double Feed_T;
double F_SepLiquidOUT;
double T_SepCooler;
double F_SepGasOUT;
double Working_Level_Volume;
double Total_Gas_Loop_Volume;
double UA;
double *VIJ;
double *aij;
double *MW;
double *SpG;
double *HVAP;
double *HCAPLa;
double *HCAPLb;
double *HCAPVa;
double *HCAPVb;
double *A;
double *B;
double *C;
double R;
double Dw;
double T_ref;
int mrows;
int ncols;
double *dstatedt_Sep;
double *F_SepLiquidOUT_out;
if(nrhs!=25)
mexErrMsgTxt("inputs not correct");
if(nlhs!=1)
mexErrMsgTxt("outputs not correct");
states = mxGetPr(prhs[0]);
Feed = mxGetScalar(prhs[1]);
Feed_comp = mxGetPr(prhs[2]);
Feed_T = mxGetScalar(prhs[3]);
F_SepLiquidOUT = mxGetScalar(prhs[4]);
T_SepCooler = mxGetScalar(prhs[5]);
F_SepGasOUT = mxGetScalar(prhs[6]);
Working_Level_Volume = mxGetScalar(prhs[7]);
Total_Gas_Loop_Volume = mxGetScalar(prhs[8]);
UA = mxGetScalar(prhs[9]);
VIJ = mxGetPr(prhs[10]);
aij = mxGetPr(prhs[11]);
MW = mxGetPr(prhs[12]);
SpG = mxGetPr(prhs[13]);
HVAP = mxGetPr(prhs[14]);
HCAPLa = mxGetPr(prhs[15]);
HCAPLb = mxGetPr(prhs[16]);
HCAPVa = mxGetPr(prhs[17]);
HCAPVb = mxGetPr(prhs[18]);
A = mxGetPr(prhs[19]);
B = mxGetPr(prhs[20]);
C = mxGetPr(prhs[21]);
R = mxGetScalar(prhs[22]);
Dw = mxGetScalar(prhs[23]);
T_ref = mxGetScalar(prhs[24]);
//-------------------------------------------
mrows = mxGetM(prhs[0]);
ncols = mxGetN(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(mrows,ncols, mxREAL);
dstatedt_Sep = mxGetPr(plhs[0]);
Separator(dstatedt_Sep,states,Feed,Feed_comp,Feed_T,
F_SepLiquidOUT,T_SepCooler,F_SepGasOUT, Working_Level_Volume,Total_Gas_Loop_Volume,UA,
VIJ, aij, MW, SpG, HVAP, HCAPLa, HCAPLb, HCAPVa, HCAPVb, A, B, C, R, Dw, T_ref);
}
*/