/[ascend]/branches/sid/models/johnpye/fprops/ttse.c
ViewVC logotype

Diff of /branches/sid/models/johnpye/fprops/ttse.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3023 by sid, Wed Jul 1 02:21:00 2015 UTC revision 3024 by sid, Thu Jul 23 04:16:35 2015 UTC
# Line 25  Line 25 
25    
26  #define INIT -23456789.123  #define INIT -23456789.123
27    
28  #define FO  //#define FO
29  //#define SO  #define SO
30  inline TtseMatrix alloc_matrix(int tp, int rhop) {  inline TtseMatrix alloc_matrix(int tp, int rhop) {
31    
32      TtseMatrix matrix =  ASC_NEW_ARRAY (TtseMatrix, tp*rhop);  //tp rows and rhop columns      TtseMatrix matrix =  ASC_NEW_ARRAY (TtseMatrix, tp*rhop);  //tp rows and rhop columns
# Line 51  inline void remove_matrix(TtseMatrix mat Line 51  inline void remove_matrix(TtseMatrix mat
51  void alloc_tables(Ttse * table)  void alloc_tables(Ttse * table)
52  {  {
53    
54        table->satFRho =  ASC_NEW_ARRAY ( double , NSAT);
55        table->satFdRhodt =  ASC_NEW_ARRAY (double, NSAT);
56        table->satFd2RhodT2 =  ASC_NEW_ARRAY (double, NSAT);
57        table->satGRho =  ASC_NEW_ARRAY (double, NSAT);
58        table->satGdRhodt =  ASC_NEW_ARRAY (double, NSAT);
59        table->satGd2RhodT2 =  ASC_NEW_ARRAY (double, NSAT);
60    
61    
62    
63      table->s = alloc_matrix(NTP,NRHOP);      table->s = alloc_matrix(NTP,NRHOP);
64      table->dsdt = alloc_matrix(NTP,NRHOP);      table->dsdt = alloc_matrix(NTP,NRHOP);
65      table->d2sdt2 = alloc_matrix(NTP,NRHOP);      table->d2sdt2 = alloc_matrix(NTP,NRHOP);
# Line 95  void alloc_tables(Ttse * table) Line 104  void alloc_tables(Ttse * table)
104    
105  void remove_tables(Ttse *table)  void remove_tables(Ttse *table)
106  {  {
107    
108        ASC_FREE(table->satFRho );
109        ASC_FREE(table->satFdRhodt );
110        ASC_FREE(table->satFd2RhodT2 );
111    
112        ASC_FREE(table->satGRho );
113        ASC_FREE(table->satGdRhodt );
114        ASC_FREE(table->satGd2RhodT2 );
115    
116      remove_matrix(table->dsdt,NTP);      remove_matrix(table->dsdt,NTP);
117      remove_matrix(table->d2sdt2,NTP);      remove_matrix(table->d2sdt2,NTP);
118      remove_matrix(table->dsdrho,NTP);      remove_matrix(table->dsdrho,NTP);
# Line 159  void build_tables(PureFluid *P){ Line 177  void build_tables(PureFluid *P){
177      int i,j;      int i,j;
178      FpropsError err = FPROPS_NO_ERROR;      FpropsError err = FPROPS_NO_ERROR;
179    
180        double Tt = P->data->T_t;
181        double Tc = P->data->T_c;
182        double dt = (Tc - Tt)/(NSAT-1);
183    
184        MSG("triple point and critical temperature -->  %f  %f",Tt,Tc);
185    
186        for(i=0; i<NSAT; ++i)
187        {
188    
189            double T = Tt + i*dt;
190    
191            double  rhof,rhog;
192            P->sat_fn(T,&rhof,&rhog,P->data,&err);
193    
194    
195            //The fluid saturation line rho's and 1st & 2nd  derivatives of rho with respect to T
196    
197            double dpdT_rho  = P->dpdT_rho_fn(T,rhof,P->data,&err);
198            double dpdrho_T  = P->dpdrho_T_fn(T,rhof,P->data,&err);
199    
200            double d2pdrho2_T = P->d2pdrho2_T_fn(T,rhof,P->data,&err);
201            double d2pdrhodT = P->d2pdTdrho_fn(T,rhof,P->data,&err);
202            double d2pdT2_rho = P->d2pdT2_rho_fn(T,rhof,P->data,&err);
203    
204    
205            double ddrho_drhodT_p_constT = ( dpdT_rho*d2pdrho2_T - dpdrho_T*d2pdrhodT ) / pow(dpdrho_T,2);
206            double ddT_drhodT_p_constrho = ( dpdT_rho*d2pdrhodT - dpdrho_T*d2pdT2_rho ) / pow(dpdrho_T,2);
207    
208            double drhodT_p =  (-dpdT_rho )/(dpdrho_T);
209    
210    
211            PT->satFRho[i] = rhof;
212            PT->satFdRhodt[i] = drhodT_p;
213            PT->satFd2RhodT2[i] =  ddT_drhodT_p_constrho  +  ddrho_drhodT_p_constT * drhodT_p;
214    
215    
216            //The Vapour saturation line rho's and 1st & 2nd derivatives of rho with respect to T
217    
218            dpdT_rho  = P->dpdT_rho_fn(T,rhog,P->data,&err);
219            dpdrho_T  = P->dpdrho_T_fn(T,rhog,P->data,&err);
220    
221            d2pdrho2_T = P->d2pdrho2_T_fn(T,rhog,P->data,&err);
222            d2pdrhodT = P->d2pdTdrho_fn(T,rhog,P->data,&err);
223            d2pdT2_rho = P->d2pdT2_rho_fn(T,rhog,P->data,&err);
224    
225            ddrho_drhodT_p_constT = ( dpdT_rho*d2pdrho2_T - dpdrho_T*d2pdrhodT ) / pow(dpdrho_T,2);
226            ddT_drhodT_p_constrho = ( dpdT_rho*d2pdrhodT - dpdrho_T*d2pdT2_rho ) / pow(dpdrho_T,2);
227    
228            drhodT_p =  (-dpdT_rho )/(dpdrho_T);
229    
230    
231            PT->satGRho[i] = rhog;
232            PT->satGdRhodt[i] = drhodT_p;
233            PT->satGd2RhodT2[i] =  ddT_drhodT_p_constrho  +  ddrho_drhodT_p_constT * drhodT_p;
234    
235    
236        }
237    
238    
239      double tmin,tmax,rhomin,rhomax;      double tmin,tmax,rhomin,rhomax;
240    
241  //Pseudo values for water  //Pseudo values for water
242  //Should be implemented else where per fluid  //Should be implemented elsewhere per fluid
243      PT->tmin = 200;      PT->tmin = 200;
244      PT->tmax = 4200;      PT->tmax = 4200;
245      PT->rhomin = 400;      PT->rhomin = 400;
# Line 175  void build_tables(PureFluid *P){ Line 252  void build_tables(PureFluid *P){
252      rhomin = PT->rhomin;      rhomin = PT->rhomin;
253      rhomax = PT->rhomax;      rhomax = PT->rhomax;
254    
255      double dt = (tmax-tmin)/NTP;      dt = (tmax-tmin)/NTP;
256      double drho = (rhomax-rhomin)/NRHOP;      double drho = (rhomax-rhomin)/NRHOP;
257    
258    
# Line 253  void build_tables(PureFluid *P){ Line 330  void build_tables(PureFluid *P){
330  #ifdef SO  #ifdef SO
331  #define EVALTTSEFN(VAR) \  #define EVALTTSEFN(VAR) \
332      double evaluate_ttse_##VAR(PureFluid *P , double t, double rho){\      double evaluate_ttse_##VAR(PureFluid *P , double t, double rho){\
333              int i,j;\          double rho_f, rho_g;\
334              double tmin = P->table->tmin; double tmax = P->table->tmax;\          FpropsError err;\
335              double rhomin  = P->table->rhomin; double rhomax= P->table->rhomax;\          int i,j;\
336              double dt = (tmax-tmin)/NTP;\          double tmin = P->data->T_t;\
337              double drho = (rhomax-rhomin)/NRHOP;\          double tmax = P->data->T_c;\
338              i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));\          if(t >= tmin  && t< tmax) {\
339              j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));\              double dt = (tmax-tmin)/NSAT;\
340                i = (int)round(((t - P->data->T_t)/(P->data->T_c - P->data->T_t)*(NSAT-1)));\
341              double delt = t - ( tmin + i*dt);\              double delt = t - ( tmin + i*dt);\
342              double delrho = rho - ( rhomin + j*drho);\              rho_f =  P->table->satFRho[i] + delt*P->table->satFdRhodt[i] + 0.5*delt*delt*P->table->satFd2RhodT2[i];\
343              double ttse##VAR = P->table->VAR[i][j]\              rho_g =  P->table->satGRho[i] + delt*P->table->satGdRhodt[i] + 0.5*delt*delt*P->table->satGd2RhodT2[i];\
344                   + delt*P->table->d##VAR##dt[i][j] + 0.5*delt*delt*P->table->d2##VAR##dt2[i][j]\              if(rho_g < rho && rho < rho_f){\
345                   + delrho*P->table->d##VAR##drho[i][j] + 0.5*delrho*delrho*P->table->d2##VAR##drho2[i][j]\                      double x = rho_g*(rho_f/rho - 1)/(rho_f - rho_g);\
346                   + delrho*delt*P->table->d2##VAR##dtdrho[i][j];\                      double Qf = P->VAR##_fn( t,rho_f,P->data,&err);\
347              return ttse##VAR;\                      double Qg = P->VAR##_fn( t,rho_g,P->data,&err);\
348                        return x*Qg + (1-x)*Qf;\
349                    }\
350                }\
351            tmin = P->table->tmin;\
352            tmax = P->table->tmax;\
353            double rhomin  = P->table->rhomin;\
354            double rhomax = P->table->rhomax;\
355            double dt = (tmax-tmin)/NTP;\
356            double drho = (rhomax-rhomin)/NRHOP;\
357            i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));\
358            j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));\
359            double delt = t - ( tmin + i*dt);\
360            double delrho = rho - ( rhomin + j*drho);\
361            double ttse##VAR = P->table->VAR[i][j]\
362                 + delt*P->table->d##VAR##dt[i][j] + 0.5*delt*delt*P->table->d2##VAR##dt2[i][j]\
363                 + delrho*P->table->d##VAR##drho[i][j] + 0.5*delrho*delrho*P->table->d2##VAR##drho2[i][j]\
364                 + delrho*delt*P->table->d2##VAR##dtdrho[i][j];\
365            return ttse##VAR;\
366          }          }
367  #endif  #endif
368    

Legend:
Removed from v.3023  
changed lines
  Added in v.3024

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22