/[ascend]/trunk/models/johnpye/fprops/solve_ph.c
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/solve_ph.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2257 - (show annotations) (download) (as text)
Wed Aug 4 09:45:52 2010 UTC (9 years, 2 months ago) by jpye
File MIME type: text/x-csrc
File size: 4758 byte(s)
Converging in supercritical region now.
Fixed some nan errors in helmholtz, needs cleaning up.
1 /*
2 ASCEND modelling environment
3 Copyright (C) 2004-2010 John Pye
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License
7 as published by the Free Software Foundation; either version 2
8 of the License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18 */
19
20 #include "solve_ph.h"
21 #include "sat2.h"
22 #include "derivs.h"
23
24 #include <stdio.h>
25 #include <math.h>
26 #include <assert.h>
27
28 #define SQ(X) ((X)*(X))
29
30 int fprops_region_ph(double p, double h, const HelmholtzData *D){
31
32 double Tsat, rhof, rhog;
33
34 if(p >= D->p_c)return FPROPS_NON;
35
36 int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);
37
38 double hf = helmholtz_h(Tsat, rhof, D);
39 if(h <= hf)return FPROPS_NON;
40
41 double hg = helmholtz_h(Tsat, rhog, D);
42 if(h >= hg)return FPROPS_NON;
43
44 return FPROPS_SAT;
45 }
46
47 int fprops_solve_ph(double p, double h, double *T, double *rho, int use_guess, const HelmholtzData *D){
48 double Tsat, rhof, rhog, hf, hg;
49 double T1, rho1;
50 int subcrit = 0;
51 if(p < D->p_c){
52 int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);
53 if(res){
54 fprintf(stderr,"Unable to solve saturation state\n");
55 return res;
56 }
57 hf = helmholtz_h(Tsat, rhof, D);
58 hg = helmholtz_h(Tsat, rhog, D);
59 fprintf(stderr,"hf = %f kJ/kg, hg = %f\n",hf/1e3,hg/1e3);
60
61 if(h > hf && h < hg){
62 fprintf(stderr,"SATURATION REGION\n");
63 /* saturation region... easy */
64 double x = (h - hf)/(hg - hf);
65 *rho = 1./(x/rhog + (1.-x)/rhof);
66 *T = Tsat;
67 return 0;
68 }
69
70 subcrit = 1;
71 if(!use_guess){
72 *T = Tsat;
73 if(h <= hf)*rho = 1.1 * rhof;
74 else *rho = rhog * 0.5;
75 }
76 }else{
77 if(!use_guess){
78 *T = D->T_c;
79 *rho = D->rho_c;
80 }
81 }
82
83 fprintf(stderr,"STARTING NON-SAT ITERATION\n");
84 //*rho = 0.21584907523805483;
85 //*T = 1004.0675923338069;
86
87 T1 = *T;
88 rho1 = *rho;
89 assert(!__isnan(T1));
90 assert(!__isnan(rho1));
91
92 /* try our own home-baked newton iteration */
93 int i = 0;
94 while(i++ < 60){
95 double p1 = helmholtz_p_raw(T1,rho1,D);
96 assert(!__isnan(p1));
97 double h1 = helmholtz_h_raw(T1,rho1,D);
98 assert(!__isnan(h1));
99
100 fprintf(stderr," T = %f, rho = %f\tp = %f bar, h = %f kJ/kg\n", T1, rho1, p1/1e5, h1/1e3);
101 fprintf(stderr," p error = %f bar\n",(p1 - p)/1e5);
102 fprintf(stderr," h error = %f kJ/kg\n",(h1 - h)/1e3);
103 if(fabs(p1 - p) < 1e-6 && fabs(h1 - h) < 1e-8){
104 fprintf(stderr,"Converged to T = %f, rho = %f, in homebaked Newton solver", T1, rho1);
105 *T = T1;
106 *rho = rho1;
107 return 0;
108 }
109 /* calculate step */
110 double f = log(p1) - log(p);
111 double g = h1 - h;
112 assert(!__isnan(f));
113 assert(!__isnan(g));
114 assert(rho1 != 0);
115 assert(p1 != 0);
116 assert(!__isnan(fprops_non_dZdT_v('p', T1,rho1, D)));
117 double f_T = 1./p1 * fprops_non_dZdT_v('p', T1,rho1, D);
118 double f_rho = -1./p1/SQ(rho1) * fprops_non_dZdv_T('p', T1, rho1, D);
119 double g_T = fprops_non_dZdT_v('h', T1,rho1, D);
120 double g_rho = -1./SQ(rho1) * fprops_non_dZdv_T('h', T1, rho1, D);
121 assert(!__isnan(f_T));
122 assert(!__isnan(f_rho));
123 assert(!__isnan(g_T));
124 assert(!__isnan(g_rho));
125
126 double det = g_rho * f_T - f_rho * g_T;
127 assert(det!=0);
128 assert(!__isnan(det));
129 fprintf(stderr," ���f/���T = %e\t\t���f/���rho = %e\n",f_T, f_rho);
130 fprintf(stderr," ���g/���T = %e\t\t���g/���rho = %e\n",g_T, g_rho);
131
132 double delta_T = -1./det * (g_rho * f - f_rho * g);
133 double delta_rho = -1./det * (f_T * g - g_T * f);
134 assert(!__isnan(delta_T));
135 assert(!__isnan(delta_rho));
136 fprintf(stderr," ��T = %f\n", delta_T);
137 fprintf(stderr," ��rho = %f\n", delta_rho);
138
139 if(subcrit){
140 if(h > hg){
141 /* vapour */
142 int i = 0;
143 while(rho1 + delta_rho > rhog && i++ < 20){
144 delta_rho *= 0.5;
145 }
146 }else{
147 /* liquid */
148 while(rho1 + delta_rho < rhof && i++ < 20){
149 delta_rho *= 0.5;
150 }
151 if(T1 + delta_T < D->T_t) delta_T = 0.5 * (T1 + D->T_t);
152 if(T1 < D->T_t) delta_T = +1;
153 }
154 }
155 /* don't go too dense */
156 if(rho1 + delta_rho > 2000) delta_rho = 2000;
157
158 /* avoid huge step */
159 while(fabs(delta_T / T1) > 0.6){
160 delta_T *= 0.5;
161 }
162 while(fabs(delta_rho / rho1) > 0.2){
163 delta_rho *= 0.5;
164 }
165
166 T1 = T1 + delta_T;
167 rho1 = rho1 + delta_rho;
168 }
169
170 *T = T1;
171 *rho = rho1;
172 return 999;
173
174 int res = fprops_nonsolver('p','h',p,h,T,rho,D);
175 fprintf(stderr,"%s: Iteration failed in nonsolver\n",__func__);
176 return res;
177 }
178

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