(*********************************************************************\
                        simple_pot.a4c
                        by Boyd T. Safrit
                        Part of the Ascend Library

This file is part of the Ascend modeling library.

Copyright (C) 1998 Carnegie Mellon University

The Ascend modeling library is free software; you can redistribute
it and/or modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

The Ascend Language Interpreter is distributed in hope that it will be
useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License along with
the program; if not, write to the Free Software Foundation, Inc., 675
Mass Ave, Cambridge, MA 02139 USA.  Check the file named COPYING.

Use of this module is demonstrated by the associated script file
simple_pot.s.
\*********************************************************************)

(*********************************************************************\
  $Date: 97/07/03 14:51:58 $
  $Revision: 1.8 $
  $Author: mthomas $
  $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/examples/simple_pot
.asc,v $
\*********************************************************************)

(*==========================================================================*

        S I M P L E _ P O T .  A  S  C
        ------------------------------
        05/93
 *==========================================================================*)
(*
Boyd T. Safrit
May 5, 1993

This file contains the models necessary to model the dynamics of a
multicomponent still pot using the Antoine vapor pressures to calculate
the volativites of each of the components in the mixture.

SCRIPT FILE:  A script file, called simple_pot.a4s, exists that enables
        the user to quickly initialize the models and solve them properly.

MODELS DEFINED:  These models can be broken down mainly into 2 parts:
        still pot  definition, and integration definition.
        MODEL                   PURPOSE
        ------------------      ----------------------------------
        molar_hold_up           Define the hold_up, amount of liquid left
                                in the pot.
        simple_pot              Define the pot, its product streams, and
                                thermodynamics.
        pot_dynamics            Highlight the LSODE integration
                                and the ASCEND models connections.
        boiling_pot             Sets up the parameters necessary for the
                                integration and gives the values of variables
                                that are specific to the particular problem
                                the user is trying to solve.

Revised to ASCEND IV/BLSODE simplicity 1/98, Ben Allan.
This set of models could really be collapsed to *ONE* reusable
MODEL and one MODEL that specifies the final set of components
and defines the values method. As it is, this revision is a
rush job.
*)


REQUIRE "atoms.a4l";

(* **********************************+************************************** *)
(* ***********************       mixture        **************************** *)
(* **********************************+************************************** *)
MODEL mixture;

        components                      IS_A set OF symbol_constant;
        choice_component                IS_A symbol_constant;
        y[components]                   IS_A fraction;
        nc                              IS_A integer_constant;

        SUM[y[i] | i IN components] = 1.0;

  METHODS
    METHOD clear_mixture;
        y[components].fixed := FALSE;
    END clear_mixture;

    METHOD specify_mixture;
        y[components].fixed := TRUE;
        y[components].fixed := FALSE;
    END specify_mixture;

    METHOD reset_mixture;
        RUN clear_mixture;
        RUN specify_mixture;
    END reset_mixture;

END mixture;

MODEL test_mixture REFINES mixture;

        components :== ['a','b','c'];
        choice_component :== 'a';
        nc :== 3;
END test_mixture;

(* **********************************+************************************** *)
(* ***********************    molar_stream      **************************** *)
(* **********************************+************************************** *)
MODEL molar_stream;

        components                              IS_A set OF symbol_constant;
        choice_component                        IS_A symbol_constant;
        state                                   IS_A mixture;
        Ftot, f[components]                     IS_A molar_rate;
        nc                                      IS_A integer_constant;

        choice_component, state.choice_component        ARE_THE_SAME;
        nc, state.nc                            ARE_THE_SAME;

        components, state.components           ARE_THE_SAME;

        FOR i IN components CREATE
            f_def[i]: f[i] = Ftot*state.y[i];
        END FOR;

  METHODS

    METHOD clear_molar_stream;
        RUN state.clear_mixture;
        Ftot.fixed                              := FALSE;
        f[components].fixed                     := FALSE;
    END clear_molar_stream;

    METHOD seqmod_molar_stream;
        RUN state.specify_mixture;
        state.y[components].fixed               := FALSE;
    END seqmod_molar_stream;

    METHOD specify_molar_stream;
        RUN seqmod_molar_stream;
        f[components].fixed                     := TRUE;
    END specify_molar_stream;

    METHOD reset_molar_stream;
        RUN clear_molar_stream;
        RUN specify_molar_stream;
    END reset_molar_stream;

    METHOD scale_molar_stream;
        FOR i IN components DO
          f[i].nominal := f[i] + 0.1{mol/s};
        END FOR;
        Ftot.nominal := Ftot + 0.1{mol/s};
    END scale_molar_stream;

END molar_stream;

MODEL test_molar_stream REFINES mixture;

        components :== ['a','b','c'];
        choice_component :== 'a';
        nc :== 3;

END test_molar_stream;

(* **********************************+************************************** *)
(* ***********************    molar_hold_up     **************************** *)
(* **********************************+************************************** *)
(* This model defines the hold_up, the amount of liquid in the pot at any given
 time. *)
MODEL molar_hold_up;

        components                      IS_A set OF symbol_constant;
        choice_component                IS_A symbol_constant;
        state                           IS_A mixture;
        Mtot, m[components]             IS_A mole;
        nc                              IS_A integer_constant;

        components, state.components    ARE_THE_SAME;

        choice_component, state.choice_component        ARE_THE_SAME;

        FOR i IN components CREATE
           m[i] = Mtot * state.y[i];
        END FOR;

   METHODS
     METHOD clear_molar_hold_up;
        RUN state.clear_mixture;
        Mtot.fixed              := FALSE;
        m[components].fixed     := FALSE;
     END clear_molar_hold_up;
     METHOD specify_molar_hold_up;
        RUN state.specify_mixture;
        Mtot.fixed := TRUE;
     END specify_molar_hold_up;
     METHOD reset_molar_hold_up;
        RUN clear_molar_hold_up;
        RUN specify_molar_hold_up;
     END reset_molar_hold_up;
END molar_hold_up;

MODEL test_molar_hold_up REFINES molar_hold_up;
        components :== ['a','b','c'];
        choice_component :== 'a';
        nc :== 3;

END test_molar_hold_up;
        
(* **********************************+************************************** *)
(* ******************    simple_pot      *********************************** *)
(* **********************************+************************************** *)
(* This model defines the pot, its vapor product stream, the dynamic equations,
 * and the thermodynamics of the pot.
 *)
MODEL simple_pot;
        F_B                                     IS_A factor;
(*
 * F_B is the direction of accumulation in the ODE.
 *)
        components                              IS_A set OF symbol_constant;
        choice_component                        IS_A symbol_constant;
        vap_prod                                IS_A molar_stream;
        hold_up                                 IS_A molar_hold_up;
        alpha[components], ave_alpha            IS_A factor;
        dM_dt[components]                       IS_A molar_rate;
        T                                       IS_A temperature;
        b[components], c[components]            IS_A real;
        a[components]                           IS_A factor;
        P                                       IS_A pressure;
        Psat[components]                        IS_A pressure;
        nc                                      IS_A integer_constant;

        
        components, hold_up.components,
        vap_prod.components                     ARE_THE_SAME;

        vap_prod.choice_component, hold_up.choice_component,
                choice_component                ARE_THE_SAME;


(* Calculate the vapor pressures of all of the components using the
 * Antoine equations.   Use the vapor presssures to calculate the
 * volativities for each component. Next write the dynamic mass balance
 * for each component.  Finally calculate the equilibrium vapor mole
 * fractions.
 *)
        FOR i IN components CREATE
        VP[i]:     Psat[i] = 10^(a[i] - b[i] / (T + c[i])) * 1.0 {mmHg};
        Volatility[i]:     alpha[i] = Psat[i] / P;
        Accumulation[i]:   dM_dt[i] = F_B * vap_prod.f[i];
        Equilibrium[i]:    vap_prod.state.y[i] =
                                alpha[i]/ave_alpha * hold_up.state.y[i];
        END FOR;

Tot_Pressure:   P = SUM[Psat[i] | i IN components];

METHODS

METHOD defaults;
        dM_dt[components].lower_bound   := -1e100 {kg_mole/s};
        b[components] :=         0{K};
        c[components] :=         0{K};
        ave_alpha.lower_bound := 0.0;
END defaults;

METHOD clear_simple_pot;
        RUN hold_up.clear_molar_hold_up;
        RUN vap_prod.clear_molar_stream;
        alpha[components].fixed                         := FALSE;
        ave_alpha.fixed                                 := FALSE;
        dM_dt[components].fixed                         := FALSE;
        T.fixed                                         := FALSE;
        a[components].fixed                             := FALSE;
        P.fixed                                         := FALSE;
        Psat[components].fixed                          := FALSE;
        F_B.fixed                                       := FALSE;
END clear_simple_pot;

METHOD specify_simple_pot;
        RUN hold_up.specify_molar_hold_up;
        RUN vap_prod.specify_molar_stream;
        hold_up.Mtot.fixed                              := FALSE;
        hold_up.state.y[components].fixed               := FALSE;
        hold_up.m[components].fixed                     := TRUE;
        vap_prod.f[components].fixed                    := FALSE;
        vap_prod.Ftot.fixed                             := TRUE;
        a[components].fixed                             := TRUE;
        P.fixed                                         := TRUE;
        Psat[components].fixed                          := FALSE;
        F_B.fixed                                       := TRUE;
END specify_simple_pot;

METHOD reset_simple_pot;
        RUN clear_simple_pot;
        RUN specify_simple_pot;
END reset_simple_pot;

END simple_pot;


MODEL test_simple_pot REFINES simple_pot;

        components :== ['a','b','c'];
        choice_component :== 'a';
        nc :== 3;

END test_simple_pot;
        
(* **********************************+************************************** *)
(* ****************    pot_dynamics      *********************************** *)
(* **********************************+************************************** *)

MODEL pot_dynamics;

   pot                                  IS_A simple_pot;

   choice_components[1..nc],
        choice_component                IS_A symbol_constant;

   components ALIASES  pot.components;
   nc                                   IS_A integer_constant;


(* Make the connections between the integration concept
 * and ASCEND models.
 * This is unnecessary, but it helps some people
 * understand how there must be a mapping between the
 * standard form of an ODE and the real problem.
 *
 * Note that the 'pot' MODEL is an autonomous DAE, that is
 * time does not appear in the equations. LSODE
 * still wants us to define an independent variable, however.
 * We choose to make it a time, but it could be any kind
 * (distance, etc).
 *
 * The integrator doesn't look for things named 'x' or 'dydx',
 * it instead looks at the ode_type and ode_id flags on variables.
 *)

   x                                    IS_A time;
   y[1..nc]                             IS_A solver_var;
   dydx[1..nc]                          IS_A solver_var;
   (* The total number of dynamic equations is the number of components. *)

   FOR i IN [1..nc] CREATE
        y[i], pot.hold_up.m[choice_components[i]]       ARE_THE_SAME;
        dydx[i], pot.dM_dt[choice_components[i]]        ARE_THE_SAME;
   END FOR;


METHODS

METHOD defaults;
        RUN pot.defaults;
END defaults;

METHOD clear;
        x.fixed := FALSE;
        y[1..nc].fixed := FALSE;
        dydx[1..c].fixed := FALSE;
END clear;

METHOD specify;

        RUN pot.specify_simple_pot;

 (* Always TRUE for calling LSODE:
  * The state variables and
  * independent variable must be set
  * TRUE, while their derivatives are FALSE.
  *)
        x.fixed := TRUE;
        y[1..nc].fixed := TRUE;
        dydx[1..nc].fixed := FALSE;

END specify;

METHOD set_ode_id;
        x.ode_type := -1;               (* independent *)
        y[1..nc].ode_type := 1;         (* state *)
        dydx[1..nc].ode_type := 2;      (* derivative *)
        FOR i IN [1..nc] DO
                y[i].ode_id := i;       (* pair up y,yprime for all i *)
                dydx[i].ode_id := i;
        END FOR;
END set_ode_id;

METHOD set_obs;
        (* you can watch any variable you like while integrating,
         * and then change what you watch and integrate again.
         *)
        pot.hold_up.m[components].obs_id := 1;
        pot.hold_up.Mtot.obs_id := 1;
END set_obs;

METHOD reset;
        RUN clear;
        RUN specify;
END reset;

END pot_dynamics;


MODEL test_pot_dynamics REFINES pot_dynamics;

        components :== ['a','b','c'];
        choice_component :== 'a';
        nc :== 3;
        choice_components[1] :== 'a';
        choice_components[2] :== 'b';
        choice_components[3] :== 'c';
        
END test_pot_dynamics;
                
(* **********************************+************************************** *)
(* ******************    bOiling_pot     *********************************** *)
(* **********************************+************************************** *)
(* This model is the model that actually gets integrated.
 * It defines the necessary values of the integration
 * parameters, as well as giving values to the variables
 * that defines the actual problem the user is trying to
 * solve.
 *)

MODEL boiling_pot REFINES pot_dynamics;

(* Define what the actual components are to be. *)
   components :== ['n_pentane', 'n_hexane', 'n_heptane'];
   nc :== CARD[components];
   choice_components[1] :== 'n_pentane';
   choice_components[2] :== 'n_hexane';
   choice_components[3] :== 'n_heptane';


METHODS

METHOD values;
        FOR i IN [1..nc] DO
(* Define the absolute and arelative error tolerance user wants LSODE to have.
 * Again, this is to help see the connection to LSODE by the
 * flags (ode_rtol,ode_atol) on the state variables y[i].
 *)
                y[i].ode_rtol := 1.0e-6;
                y[i].ode_atol := 1.0e-6;
        END FOR;

        pot.alpha['n_pentane']          := 1.44880;
        pot.alpha['n_hexane']           := 0.97472;
        pot.alpha['n_heptane']          := 0.70000;
        pot.ave_alpha                   := 1.0;
        pot.a['n_pentane']              := 6.85221;
        pot.a['n_hexane']               := 6.87776;
        pot.a['n_heptane']              := 6.90240;
        pot.b['n_pentane']              := 791.48 {K};
        pot.b['n_hexane']               := 898.38 {K};
        pot.b['n_heptane']              := 994.965 {K};
        pot.c['n_pentane']              :=-41.15 {K};
        pot.c['n_hexane']               :=-48.784 {K};
        pot.c['n_heptane']              :=-56.25 {K};
        pot.vap_prod.Ftot               := 0.8 {kg_mole/s};
        pot.hold_up.m['n_pentane']      := 50.0 {kg_mole};
        pot.hold_up.m['n_hexane']       := 50.0 {kg_mole};
        pot.hold_up.m['n_heptane']      := 50.0 {kg_mole};
        pot.dM_dt['n_pentane']          := 0.0 {kg_mole/s};
        pot.dM_dt['n_hexane']           := 0.0 {kg_mole/s};
        pot.dM_dt['n_heptane']          := 0.0 {kg_mole/s};
        pot.P                           := 1.0 {atm};
        pot.T                           := 300 {K};
        x                               := 0 {s};
        pot.F_B                         := -1.0; (* forward *)
END values;

METHOD clear;
        RUN pot.clear_simple_pot;
END clear;


END boiling_pot;