/[ascend]/trunk/models/steam/dsgsat3.a4c
ViewVC logotype

Diff of /trunk/models/steam/dsgsat3.a4c

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

revision 1304 by johnpye, Sat Feb 10 10:47:23 2007 UTC revision 1305 by johnpye, Sat Mar 3 04:57:04 2007 UTC
# Line 11  REQUIRE "steam/satsteamstream.a4c"; Line 11  REQUIRE "steam/satsteamstream.a4c";
11    
12  MODEL dsgsat3;  MODEL dsgsat3;
13      n IS_A integer_constant;      n IS_A integer_constant;
14      n :== 7;      n :== 10;(* with L = 10m: 5,6,7,8,9,10,11 *)
15                 (* with L =  5m: 2,3,4,5,7,9,11,12,13,1415,16 *)
16    
17        dz IS_A real_constant;
18        L IS_A real_constant;
19        L :== 5 {m};
20        dz :== L / (n-1);
21    
22        nodes,butfirst1,upwind4,central IS_A set OF integer_constant;
23        nodes :== [1..n];
24        butfirst1 :== nodes - [1];
25        upwind4 :== nodes - [1,2,n];
26        central :== nodes - [1,n];
27    
28      (* temporal derivatives *)      (* temporal derivatives *)
29      drho_dt[2..n] IS_A density_rate;      drho_dt[butfirst1] IS_A density_rate;
30      dmdot_dt[2..n] IS_A mass_rate_rate;      dmdot_dt[butfirst1] IS_A mass_rate_rate;
31      du_dt[2..n] IS_A specific_energy_rate;      du_dt[butfirst1] IS_A specific_energy_rate;
32      dTw_dt[2..n] IS_A temperature_rate;      dTw_dt[butfirst1] IS_A temperature_rate;
33    
34      (* wall properties *)      (* wall properties *)
35      rho_w IS_A mass_density;      rho_w IS_A mass_density;
# Line 28  MODEL dsgsat3; Line 40  MODEL dsgsat3;
40      h_ext IS_A heat_transfer_coefficient; (* external *)      h_ext IS_A heat_transfer_coefficient; (* external *)
41      z_A: A = 1{PI}*D^2/4;      z_A: A = 1{PI}*D^2/4;
42      z_Aw: A_w = 1{PI}*(D_2^2 - D^2)/4;      z_Aw: A_w = 1{PI}*(D_2^2 - D^2)/4;
     dz IS_A distance;  
     L IS_A distance;  
     z_dz: dz = L / (n - 1);  
43    
44      (* fluid properties *)      (* fluid properties *)
45      node[1..n] IS_A satsteamstream;      node[nodes] IS_A satsteamstream;
46            
47      (* flow properties *)      (* flow properties *)
48      vel[1..n] IS_A speed;      vel[nodes] IS_A speed;
49      T_w[2..n] IS_A temperature;      T_w[butfirst1] IS_A temperature;
50      T[1..n] IS_A temperature;      T[nodes] IS_A temperature;
51    
52      (* constants, for the moment: *)      (* constant, for the moment: *)
53      f IS_A positive_factor;      f IS_A positive_factor;
54      (* mu_f IS_A viscosity; *)      (* mu_f IS_A viscosity; *)
55      T_amb IS_A temperature;      T_amb IS_A temperature;
56    
57      (* system dynamics *)        (* system dynamics *)  
58      qdot_t[2..n], qdot_l[2..n] IS_A power_per_length;      qdot_t[butfirst1], qdot_l[butfirst1] IS_A power_per_length;
59      qdot_s IS_A power_per_length;      qdot_s IS_A power_per_length;
60    
61      FOR i IN [1..n] CREATE      FOR i IN nodes CREATE
62          z_vel[i]: vel[i] = node[i].v*node[i].mdot/A;          z_vel[i]: vel[i] = v[i]*mdot[i]/A;
63      END FOR;      END FOR;
64    
65      (* some aliases just for easier review of the state of the model *)      (* some aliases just for easier review of the state of the model *)
66      x[1..n] IS_A fraction;      x[nodes] IS_A fraction;
67      mdot[1..n] IS_A mass_rate;      mdot[nodes] IS_A mass_rate;
68      p[1..n] IS_A pressure;      p[nodes] IS_A pressure;
69      FOR i IN [1..n] CREATE      rho[nodes] IS_A mass_density;
70          x[i], node[i].x ARE_THE_SAME;      u[nodes] IS_A specific_energy;
71        v[nodes] IS_A specific_volume;
72        FOR i IN nodes CREATE
73            x[i],    node[i].x ARE_THE_SAME;
74          mdot[i], node[i].mdot ARE_THE_SAME;          mdot[i], node[i].mdot ARE_THE_SAME;
75          p[i], node[i].p ARE_THE_SAME;          p[i],    node[i].p ARE_THE_SAME;
76          T[i], node[i].T ARE_THE_SAME;          T[i],    node[i].T ARE_THE_SAME;
77            rho[i],  node[i].rho ARE_THE_SAME;
78            u[i],    node[i].u ARE_THE_SAME;
79            v[i],    node[i].v ARE_THE_SAME;
80      END FOR;      END FOR;
81    
82      (* differential equations *)      (* mass conservation *)
83            FOR i IN upwind4 CREATE (* 4-pt upwind biased *)
84      (* use central difference for several derivatives *)          z_massbal1[i]: A * drho_dt[i] * dz =
85      FOR i IN [2..n-1] CREATE                  - (mdot[i+1] + 6.*mdot[i] - 3.*mdot[i-1] - 2.*mdot[i-2]) / 6.;
86          z_massbal[i]: A * drho_dt[i] * dz = - (node[i+1].mdot - node[i-1].mdot) / 2.0;      END FOR;
87        FOR i IN [] CREATE
88          z_enbal[i]: dz * (qdot_t[i] - node[i].rho * A * du_dt[i]) =          z_massbal2[i]: A * drho_dt[i] * dz =
89               + node[i].mdot * (node[i].u - node[i-1].u)                  - (mdot[i+1] - mdot[i-1]) / 2.;
90               + (node[i+1].p*node[i+1].v*node[i+1].mdot - node[i-1].p*node[i-1].v*node[i-1].mdot) / 2;      END FOR;
91        FOR i IN [2,n] CREATE
92          z_mombal[i]:  - dz/A*dmdot_dt[i] =          z_massbal[i]: A * drho_dt[i] * dz = - (mdot[i] - mdot[i-1]);
                         (node[i+1].p-node[i-1].p) / 2.0  
                         + dz * f/D/2 * node[i].rho * vel[i]^2  
                         + (node[i+1].rho*vel[i+1]^2 - node[i-1].rho*vel[i-1]^2) / 2.0;  
93      END FOR;      END FOR;
94    
95      (* the equations for the last node use backwards difference *)      (* energy conservation *)
96        FOR i IN [] CREATE
97            z_enbal2[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
98                 + mdot[i] * (node[i+1].u + 6.*u[i] - 3.*u[i-1] - 2.*u[i-1]) / 6.
99                 + (p[i+1]*node[i+1].v*mdot[i+1] - p[i-1]*v[i-1]*mdot[i-1]) / 2.;
100        END FOR;
101        FOR i IN central CREATE
102            z_enbal1[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
103                 + mdot[i] * (u[i] - u[i-1]) (* NOTE: not central *)
104                 + (p[i+1]*v[i+1]*mdot[i+1] - p[i-1]*v[i-1]*mdot[i-1]) / 2.;
105        END FOR;
106      FOR i IN [n] CREATE      FOR i IN [n] CREATE
107          z_massbal1[i]: A * drho_dt[i] * dz = - (node[i].mdot - node[i-1].mdot);          z_enbal[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
108                 + mdot[i] * (u[i] - u[i-1])
109          z_enbal1[i]: dz * (qdot_t[i] - node[i].rho * A * du_dt[i]) =               + (p[i]*v[i]*mdot[i] - p[i-1]*v[i-1]*mdot[i-1]);
110               + node[i].mdot * (node[i].u - node[i-1].u)      END FOR;
              + (node[i].p*node[i].v*node[i].mdot - node[i-1].p*node[i-1].v*node[i-1].mdot);  
111    
112          z_mombal1[i]:  - dz/A*dmdot_dt[i] =      (* momentum conservation *)
113                          (node[i].p-node[i-1].p)      FOR i IN upwind4 CREATE
114                          + dz * f/D/2 * node[i].rho * vel[i]^2          z_mombal2[i]:  - dz/A * dmdot_dt[i]
115                          + (node[i].rho*vel[i]^2 - node[i-1].rho*vel[i-1]^2);               =  (p[i]-p[i-1])  (* backdiff for pressure *)
116                    + dz * f/D/2 * rho[i] * vel[i]^2
117                    + (rho[i+1]*vel[i+1]^2 + 6.*rho[i]*vel[i]^2 - 3.*rho[i-1]*vel[i-1]^2 - 2.*rho[i-2]*vel[i-2]^2) / 6.;
118        END FOR;
119        FOR i IN [] CREATE
120            z_mombal1[i]:  - dz/A * dmdot_dt[i]
121                 =  (p[i+1]-p[i-1]) / 2.
122                    + dz * f/D/2 * rho[i] * vel[i]^2
123                    + (rho[i+1]*vel[i+1]^2 - rho[i-1]*vel[i-1]^2) / 2.;
124        END FOR;
125        FOR i IN [2,n] CREATE
126            z_mombal[i]:  - dz/A * dmdot_dt[i]
127                 =  (p[i]-p[i-1])
128                    + dz * f/D/2 * rho[i] * vel[i]^2
129                    + (rho[i]*vel[i]^2 - rho[i-1]*vel[i-1]^2);
130      END FOR;      END FOR;
131    
132      (* these equation have no spatial derivatives *)      (* internal/external convection, and thermal mass of wall -- no spatial derivs here *)
133      FOR i IN [2..n] CREATE      FOR i IN butfirst1 CREATE
134          z_wall[i]:  rho_w*A_w*c_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];          z_wall[i]:  rho_w*A_w*c_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
135          z_loss[i]:  qdot_l[i] = h_ext*(1{PI}*D_2)*(T_w[i] - T_amb);          z_loss[i]:  qdot_l[i] = h_ext*(1{PI}*D_2)*(T_w[i] - T_amb);
136          z_trans[i]: qdot_t[i] = h_int*(1{PI}*D)  *(T_w[i] - node[i].T);          z_trans[i]: qdot_t[i] = h_int*(1{PI}*D)  *(T_w[i] - T[i]);
137      END FOR;      END FOR;
138    
139      t IS_A time;      t IS_A time;
140  METHODS  METHODS
141      METHOD bound_self;      METHOD bound_self;
142          vel[1..n].upper_bound := 100 {m/s};          vel[nodes].upper_bound := 100 {m/s};
143          qdot_l[2..n].lower_bound := 0 {W/m};          qdot_l[butfirst1].lower_bound := 0 {W/m};
144          FOR i IN [1..n] DO          FOR i IN nodes DO
145              RUN node[i].bound_self;              RUN node[i].bound_self;
146          END FOR;          END FOR;
147      END bound_self;      END bound_self;
148      METHOD default_self;      METHOD default;
149          (* these are DEFAULT values and are overwritten in METHOD values below *)          (* these are initial guesses only; fixed parameters are overwritten by 'values' below *)
150          D := 0.06 {m};          D := 0.06 {m};
151          D_2 := 0.07 {m};          D_2 := 0.07 {m};
152          A_w := 0.25{PI}*D_2^2 -0.25{PI}*D^2;          A_w := 0.25{PI}*D_2^2 -0.25{PI}*D^2;
153          A := 1 {m^2};          A := 1 {m^2};
         L := 50 {m};  
154          T_amb := 298 {K};          T_amb := 298 {K};
155          dz := 3.08 {m};          f := 0.0;
         f := 0.01;  
156          h_ext := 10 {W/m^2/K};          h_ext := 10 {W/m^2/K};
157          h_int := 100 {W/m^2/K};          h_int := 100 {W/m^2/K};
158          qdot_s := 1000 {W/m};          qdot_s := 1000 {W/m};
159          rho_w := 1000 {kg/m^3};          rho_w := 1000 {kg/m^3};
160          t := 0 {s};          t := 0 {s};
161          FOR i IN [1..n] DO          FOR i IN nodes DO
162              T[i] := 298 {K};              T[i] := 298 {K};
163              vel[i] := 1 {m/s};              vel[i] := 1 {m/s};
164              RUN node[i].default_self;              RUN node[i].default_self;
165          END FOR;          END FOR;
166          FOR i IN [2..n] DO          FOR i IN butfirst1 DO
167              T_w[i] := 298 {K};              T_w[i] := 298 {K};
168              drho_dt[i] := 0 {kg/m^3/s};              drho_dt[i] := 0 {kg/m^3/s};
169              dmdot_dt[i] := 0 {kg/s/s};              dmdot_dt[i] := 0 {kg/s/s};
# Line 138  METHODS Line 171  METHODS
171              dTw_dt[i] := 0 {K/s};              dTw_dt[i] := 0 {K/s};
172              qdot_t[i] := 0 {W/m};              qdot_t[i] := 0 {W/m};
173              qdot_l[i] := 0 {W/m};              qdot_l[i] := 0 {W/m};
174                x[i] := x[1];
175          END FOR;          END FOR;
176      END default_self;      END default;
177      METHOD specify;      METHOD specify;
178          (* change to a proper steady-state problem, with fluid properties FREEd *)          (* change to a proper steady-state problem, with fluid properties FREEd *)
179          FOR i IN [1..n] DO          FOR i IN nodes DO
180              RUN node[i].specify;              RUN node[i].specify;
181              FIX dTw_dt[i];   FREE T_w[i];              FIX dTw_dt[i];   FREE T_w[i];
182          END FOR;          END FOR;
183          FIX node[1].p;          FIX p[1];
184          FREE node[1].T;          FREE T[1];
185          FIX qdot_s;          FIX qdot_s;
186          FIX D, D_2, L;          FIX D, D_2;
187          FIX h_int, c_w, rho_w, h_ext;          FIX h_int, c_w, rho_w, h_ext;
188          FIX f;          FIX f;
189          (* FIX mu_f; *)          (* FIX mu_f; *)
190          FIX T_amb;          FIX T_amb;
191          (* fix derivatives to zero *)          (* fix derivatives to zero *)
192          FOR i IN [2..n] DO          FOR i IN [2..n] DO
193              (* FIX dmdot_dt[i]; FREE node[i].mdot; *)              FREE x[i]; FIX p[i];
194              FREE node[i].x; FIX node[i].p;              FIX drho_dt[i];  FREE p[i];
195              FIX drho_dt[i];  FREE node[i].p;              FIX du_dt[i]; FREE T[i];
             FIX du_dt[i]; FREE node[i].T;  
196              FREE mdot[i]; FIX dmdot_dt[i];              FREE mdot[i]; FIX dmdot_dt[i];
197          END FOR;          END FOR;
198      END specify;      END specify;
199      METHOD values;      METHOD values;
         L := 10 {m};  
200          h_int := 100 {W/m^2/K};          h_int := 100 {W/m^2/K};
201          h_ext := 20 {W/m^2/K};          h_ext := 20 {W/m^2/K};
202          f := 0.03;          f := 0.0;
203          node[1].mdot := 0.26 {kg/s};          mdot[1] := 0.26 {kg/s};
204          node[1].p := 10 {bar};          p[1]    := 10 {bar};
205          node[1].x := 0.23;          x[1]    := 0.1;
206          qdot_s := 1000 {W/m^2} * D_2 * 10;          qdot_s := 1000 {W/m^2} * D_2 * 10;
207          FOR i IN [2..n] DO          FOR i IN [2..n] DO
208              dmdot_dt[i] := 0.0 {kg/s/s};              dmdot_dt[i] := 0.0 {kg/s/s};
209              du_dt[i] := 0 {kJ/kg/s};              du_dt[i] := 0 {kJ/kg/s};
210              node[i].v := 0.2 {L/kg};              v[i] := 0.2 {L/kg};
211              node[i].rho := 6 {kg/L};              rho[i] := 6 {kg/L};
212              node[i].dp_dT := +0.5 {kPa/K};              node[i].dp_dT := +0.5 {kPa/K};
213              node[i].p := 5 {bar};              p[i] := 5 {bar};
214          END FOR;          END FOR;
215      END values;      END values;
216      METHOD on_load;      METHOD on_load;
# Line 200  METHODS Line 232  METHODS
232      (*------------------------- the dynamic problem ------------------------------*)      (*------------------------- the dynamic problem ------------------------------*)
233      METHOD configure_dynamic;      METHOD configure_dynamic;
234          FOR i IN [2..n] DO          FOR i IN [2..n] DO
235              FREE drho_dt[i];  FIX node[i].rho;              FREE drho_dt[i];  FIX rho[i];
236              FREE dmdot_dt[i]; FIX node[i].mdot;              FREE dmdot_dt[i]; FIX mdot[i];
237              FREE du_dt[i];    FIX node[i].u;              FREE du_dt[i];    FIX u[i];
238              FREE dTw_dt[i];   FIX T_w[i];              FREE dTw_dt[i];   FIX T_w[i];
239              FREE node[i].x;              FREE x[i];
240              FREE node[i].T;              FREE T[i];
241          END FOR;          END FOR;
242          t := 0 {s};          t := 0 {s};
243      END configure_dynamic;      END configure_dynamic;
244      METHOD free_states;      METHOD free_states;
245          FOR i IN [2..n] DO          FOR i IN [2..n] DO
246              FREE node[i].rho;              FREE rho[i];
247              FREE node[i].mdot;              FREE mdot[i];
248              FREE node[i].u;              FREE u[i];
249              FREE T_w[i];              FREE T_w[i];
250          END FOR;          END FOR;
251      END free_states;          END free_states;    
# Line 222  METHODS Line 254  METHODS
254          t.ode_type := -1;          t.ode_type := -1;
255    
256          FOR i IN [2..n] DO          FOR i IN [2..n] DO
257              drho_dt[i].ode_id := 4*i;     node[i].rho.ode_id := 4*i;              drho_dt[i].ode_id := 4*i;     rho[i].ode_id := 4*i;
258              drho_dt[i].ode_type := 2;     node[i].rho.ode_type := 1;              drho_dt[i].ode_type := 2;     rho[i].ode_type := 1;
259    
260              dmdot_dt[i].ode_id := 4*i+1;  node[i].mdot.ode_id := 4*i+1;              dmdot_dt[i].ode_id := 4*i+1;  mdot[i].ode_id := 4*i+1;
261              dmdot_dt[i].ode_type := 2;    node[i].mdot.ode_type := 1;              dmdot_dt[i].ode_type := 2;    mdot[i].ode_type := 1;
262                            
263              du_dt[i].ode_id := 4*i+2;     node[i].u.ode_id := 4*i+2;              du_dt[i].ode_id := 4*i+2;     u[i].ode_id := 4*i+2;
264              du_dt[i].ode_type := 2;       node[i].u.ode_type := 1;              du_dt[i].ode_type := 2;       u[i].ode_type := 1;
265    
266              dTw_dt[i].ode_id := 4*i+3;    T_w[i].ode_id := 4*i+3;              dTw_dt[i].ode_id := 4*i+3;    T_w[i].ode_id := 4*i+3;
267              dTw_dt[i].ode_type := 2;      T_w[i].ode_type := 1;              dTw_dt[i].ode_type := 2;      T_w[i].ode_type := 1;
268          END FOR;          END FOR;
269    
270          FOR i IN [1..n] DO          FOR i IN nodes DO
271              (* p[i].obs_id :=         1 + 10*i; *)              (* p[i].obs_id :=         1 + 10*i; *)
272              (* x[i].obs_id :=         2 + 10*i; *)              (* x[i].obs_id :=         2 + 10*i; *)
273              (* node[i].mdot.obs_id := 4 + 10*i; *)              (* mdot[i].obs_id := 4 + 10*i; *)
274              (* node[i].h.obs_id :=    3 + 10*i; *)              (* node[i].h.obs_id :=    3 + 10*i; *)
275          END FOR;          END FOR;
276          FOR i IN [2..n] DO          FOR i IN [2..n] DO
# Line 247  METHODS Line 279  METHODS
279              (* T[i].obs_id :=         6 + 10*i;*)              (* T[i].obs_id :=         6 + 10*i;*)
280          END FOR;          END FOR;
281    
282          node[n].mdot.obs_id := 1;          mdot[n].obs_id := 1;
283          node[n].x.obs_id := 1;          x[n].obs_id := 1;
284          p[n].obs_id := 1;          p[n].obs_id := 1;
285          vel[n].obs_id := 1;          vel[n].obs_id := 1;
286    
287      END ode_init;      END ode_init;
288      METHOD fix_outlet_quality;      METHOD fix_outlet_quality;
289          FIX x[n];          FIX x[n];
290          FREE node[1].mdot;          FREE mdot[1];
291      END fix_outlet_quality;      END fix_outlet_quality;
292    
293  END dsgsat3;  END dsgsat3;

Legend:
Removed from v.1304  
changed lines
  Added in v.1305

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