/[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 2260 - (show annotations) (download) (as text)
Thu Aug 5 06:28:01 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 5336 byte(s)
Error messages to detect temperatures below critical point.
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 = 1.1 * Tsat;
73 if(h <= hf){
74 *rho = rhof;
75 fprintf(stderr,"LIQUID GUESS: T = %f, rho = %f\n",*T, *rho);
76 }else{
77 *rho = rhog * 0.5;
78 fprintf(stderr,"GAS GUESS: T = %f, rho = %f\n",*T, *rho);
79 }
80 }
81 }else{
82 if(!use_guess){
83 *T = D->T_c * 1.03;
84 *rho = D->rho_c;
85 }
86 }
87
88 fprintf(stderr,"STARTING NON-SAT ITERATION\n");
89 //*rho = 976.82687191126922;
90 //*T = 344.80371310850518;
91
92 T1 = *T;
93 rho1 = *rho;
94 assert(!__isnan(T1));
95 assert(!__isnan(rho1));
96
97 /* try our own home-baked newton iteration */
98 int i = 0;
99 double delta_T = 0;
100 double delta_rho = 0;
101 fprintf(stderr,"STARTING ITERATION\n");
102 while(i++ < 60){
103 double p1 = helmholtz_p_raw(T1,rho1,D);
104 assert(!__isnan(p1));
105 double h1 = helmholtz_h_raw(T1,rho1,D);
106 assert(!__isnan(h1));
107
108 fprintf(stderr," T = %f, rho = %f\tp = %f bar, h = %f kJ/kg\n", T1, rho1, p1/1e5, h1/1e3);
109 fprintf(stderr," p error = %f bar\n",(p1 - p)/1e5);
110 fprintf(stderr," h error = %f kJ/kg\n",(h1 - h)/1e3);
111
112 if(p1 < 0){
113 T1 -= (delta_T *= 0.5);
114 rho1 -= (delta_rho *= 0.5);
115 continue;
116 }
117
118 if(fabs(p1 - p) < 1e-6 && fabs(h1 - h) < 1e-8){
119 fprintf(stderr,"Converged to T = %f, rho = %f, in homebaked Newton solver", T1, rho1);
120 *T = T1;
121 *rho = rho1;
122 return 0;
123 }
124 /* calculate step */
125 double f = log(p1) - log(p);
126 double g = h1 - h;
127 assert(!__isnan(f));
128 assert(!__isnan(g));
129 assert(rho1 != 0);
130 assert(p1 != 0);
131 assert(!__isnan(fprops_non_dZdT_v('p', T1,rho1, D)));
132 double f_T = 1./p1 * fprops_non_dZdT_v('p', T1,rho1, D);
133 double f_rho = -1./p1/SQ(rho1) * fprops_non_dZdv_T('p', T1, rho1, D);
134 double g_T = fprops_non_dZdT_v('h', T1,rho1, D);
135 double g_rho = -1./SQ(rho1) * fprops_non_dZdv_T('h', T1, rho1, D);
136 assert(!__isnan(f_T));
137 assert(!__isnan(f_rho));
138 assert(!__isnan(g_T));
139 assert(!__isnan(g_rho));
140
141 double det = g_rho * f_T - f_rho * g_T;
142 assert(det!=0);
143 assert(!__isnan(det));
144 fprintf(stderr," ���f/���T = %e\t\t���f/���rho = %e\n",f_T, f_rho);
145 fprintf(stderr," ���g/���T = %e\t\t���g/���rho = %e\n",g_T, g_rho);
146
147 delta_T = -1./det * (g_rho * f - f_rho * g);
148 delta_rho = -1./det * (f_T * g - g_T * f);
149 assert(!__isnan(delta_T));
150 assert(!__isnan(delta_rho));
151 fprintf(stderr," ��T = %f\n", delta_T);
152 fprintf(stderr," ��rho = %f\n", delta_rho);
153
154 if(subcrit){
155 if(h > hg){
156 /* vapour */
157 int i = 0;
158 while(rho1 + delta_rho > rhog && i++ < 20){
159 delta_rho *= 0.5;
160 }
161 }else{
162 /* liquid */
163 while(rho1 + delta_rho < rhof && i++ < 20){
164 delta_rho *= 0.5;
165 }
166 if(T1 + delta_T < D->T_t) delta_T = 0.5 * (T1 + D->T_t);
167 if(T1 < D->T_t) delta_T = +1;
168 }
169 }else{
170 /* supercritical... stay above critical temperature of density > rho_crit */
171 if(rho1 + delta_rho > D->rho_c && T1 + delta_T < D->T_c){
172 delta_T = 0.5 *(T1 + D->T_c);
173 }
174 }
175 /* don't go too dense */
176 if(rho1 + delta_rho > 2000) delta_rho = 2000;
177
178 /* don't go too hot */
179 if(T1 + delta_T > 5000) delta_T = 5000 - T1;
180
181 /* avoid huge step */
182 while(fabs(delta_T / T1) > 0.6){
183 delta_T *= 0.5;
184 }
185 while(fabs(delta_rho / rho1) > 0.2){
186 delta_rho *= 0.5;
187 }
188
189 T1 = T1 + delta_T;
190 rho1 = rho1 + delta_rho;
191 }
192
193 *T = T1;
194 *rho = rho1;
195 return 999;
196
197 int res = fprops_nonsolver('p','h',p,h,T,rho,D);
198 fprintf(stderr,"%s: Iteration failed in nonsolver\n",__func__);
199 return res;
200 }
201

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