Tummescheit H., Eborn J.Chemical Reaction Modeling with ThermoFluid/MF and MultiFlashChemical Reaction Modeling withThermoFluid/MF and MultiFlashHubertus Tummescheit† and Jonas Eborn‡† Department‡ Unitedof Automatic ControlLund University, [email protected] Research CenterEast Hartford, Connecticut, [email protected] free Modelica library T HERMO F LUID (see [2]and [11]) was developed for simulation of thermohydraulic applications, both for single-species applications like the water-steam cycle in a thermal powerplant and for multi-species applications with gas mixtures. It has demonstrated its flexibility for modeling thermodynamic and process applications in a variety of industrial and academic projects, see [10], [3]and [7]. This article describes how support for chemical reactions and membrane diffusion has been addedto T HERMO F LUID , thus expanding the area of possible applications to include reacting flows, chemical batch reactors, catalytic converters, etc. Anothercrucial part of the modeling work has to be spent ongetting physical property data of sufficient accuracyand with acceptable computational complexity for engineering purposes into the model. This has beenadressed in the development of a commercial interface to the industry-standard physical property package MultiFlash. The new Modelica library T HER MO F LUID /MF provides the modeler with two toolboxes. Firstly, a low-level Modelica function interfaceto MultiFlash. MultiFlash consists of a core of physical property calculation routines and a basic databaseof the most comman chemical components and a number of add-on property databases. The interface givesaccess to multi-component, multi-phase property calculations including gas, several liquid and condensedphases, wax formations and hydrates. Secondly, ahigh-level Modelica model library which is fully integrated with the T HERMO F LUID library and implements robust and efficient dynamical models for themost common process engineering equipment. In addition, reliable crossing functions for detecting phaseboundaries in multi-phase, multi-component mixturesThe Modelica Association31have been implemented for the first time in a highlevel modeling language. The crossing functions makeit possible to simulate processes correctly even at offdesign operating points and under start-up conditions.A flash volume may in such cases be filled with onlyliquid or only gas. Crossing functions for phase transitions ensure high performance simulation even in thesecases.1Flexible handling of chemical reactionsIn standard chemical textbooks, reactions are treatedas source terms in component concentration balances:dcidtcini couti(1)ri where ri are the component reaction rates, given bya kinetic expression. In a more general way, we caninclude the reaction terms in the component mass balance and total energy balancedMidtdUdtṁiniq̇in ṁoutiq̇outrZi MWi nc rZii 1f Hi(2)(3) where rZ are reaction rates in moles/s, q̇ is convectiveheat flow, ṁ mass flow and H f is component enthalpyof formation.1.1ThermoFluid balance equationsIn T HERMO F LUID, the general balance equations are implemented in the package BaseClasses.Balances. The basic balance equationsshould not be modified by the average user and thusModelica 2002, March 18 19, 2002

Chemical Reaction Modeling with ThermoFluid/MF and MultiFlashTummescheit H., Eborn J.they need to be attached to an object inside the volume model. This is the HeatAndMassObject,seeFigure 1, which acts as a gateway between the balance equations and possible heat- and mass transferobjects. External connectors can also be connected tothe HeatAndMassObject.InterfacesThe HeatAndMassObject interact with other objects through a number of different connectors. Thecurrently implemented connectors include the HeatFlow connector for pure heat interaction (conduction/radiation), the ChemFlow connector for chemical reactions, both kinetic and equilibrium, and a connector for membrane diffusion.Figure 1: Schematic of the HeatAndMassObject withheat interaction and reaction objects (no diffusion connector present).connector HeatFlowparameter IntegerTemperature[n]flow Power[n]end HeatFlow;n;T;q;need to be general enough to handle all cases from anconnector ChemFlowisolated gas volume to a reactor with added mass- andparameter Integerheat transfer laws. The model structure has to provideparameter Stringthe option to add any kind of heat- and mass transferTemperature[n]interaction with the control volume later, as an add-onPressure[n]Concentration[n,nc]component. The basic balance equation of a controlflow MolarFlowRate[n,nc]volume with two connectors is implemented asflow Power[n]end ChemFlowdM x a.mdot x b.mdot x rM;dU a.q conv b.q conv Q s;n, nc;MediumType;T;p;conc;rZ;q;The diffusion connector is similar to the ChemFlowconnector but has mass flow rate instead of molar flowrate since this is standard for diffusion.The flow semantics of Modelica for the molar flow rZand the heat flow q make sure that all contributionsto the mass- and energy balances are correctly takeninto account, no matter whether there are zero, oneor many connections to the HeatAndMassObject, seeFigure 1. Inside the HeatAndMassObject the contributions from the different connectors are summed1.2 The HeatAndMassObject, a gateway toup and transferred to the balance equations in the volthe balancesume.The contradiction of leaving the option open to specify production terms but not having to add a default 1.3 Objects for encapsulating reactionsvalue of zero can be handled with open flow connectors. In Modelica, all quantities which are flows are To be able to drag and drop reaction models into a volmarked with the flow-prefix. Flow variables obey ume model (a reactor), they are encapsulated in reacthe zero-sum rule (Kirchhoffs’ current law) and have tion objects. As shown in the code example below, thein unconnected connectors a zero default value. Since Basic reaction inherits interfaces and basic paramethese connectors should be internal to the volume, ters from the reaction BaseObject.This means that rM and Q s, corresponding to thesource terms in (3) should be unspecified in the general base class, and then specified at a later stage whenthe balance class is reused in the model of a specificcomponent. When there are no reactions or heat interaction with the volume there is no need for any sourceterms. In this case the model should provide a defaultvalue of zero production.Modelica 2002, March 18 19, 200232The Modelica Association

Tummescheit H., Eborn J.Chemical Reaction Modeling with ThermoFluid/MF and MultiFlashpackage Reactionspartial model BaseObjectparameter Integern, nc, nr;MolarReactionRate[n,nr] reacRate;Enthalpy[nc]compHf;parameter StoichiometricNumber[nr,nc]stoich zeros(nr,nc);equationfor i in 1:n loopr.rZ[i,:] transpose(stoich)*reacRate[i,:];end for;end BaseObject;model Basic "Simple Arrhenius reaction"extends BaseObject;parameter Rate[nr]A0;parameter Real[nr]b;parameter MolarInternalEnergy[nr] Ea;Concentration[n,nc]conc;equationfor i in 1:n loopreacRate[i,:] .;r.q[i] -compHf*r.rZ[i,:];end for;end Basic;Figure 2: Schematic of example system with H2 –O2reaction.built and reactions can graphically be added to standard reactor models.1.4end Reactions;Example, combustion of hydrogenAs an example, we consider the combustion of hydrogen and oxygen into water. In a simple setting, seeFigure 2, the system consists of a reservoir supplyingthe reactants, a reactor volume and a sink for the product flow. A heat source is added to provide the heatnecessary to ignite the mixture.The complete set of sub-reactions for this process involves a large number ( 40) of very fast reactions,see [12]. Here we only consider the 8 main reactions,involving the components O2 , H2 , H2 O, O , H , OH,Ar . Argon is included as an inert gas. The includedreactions are listed in Table 1. The corresponding stoichiometry matrix and reaction rate parameters havebeen coded into a Basic reaction object inside theGasCV reaction vessel.The result plots show clearly that the reactions are extremely fast once they started. They saturate whenall H2 is burned up and the flow through the volumereaches steady state. The mass flows in Figure 3 showa violent explosion when the mixture ignites. After the The reaction rates are calculated from standard Arrhenius expressions using the concentrations and the parameters. To use the reaction component, like in thekinetic reaction in Figure 1, the user simply needs tospecify the parameters. The stoichiometry matrix isconstructed as shown in Table 1 and the heat of formation parameters are added from the medium model.To construct models of other types of reactions the reaction BaseObject can be reused. The customizedreaction model needs to give expressions for the reaction rates, either by adding equations or by calling arate function. In this way packages of reactions can be 2.5outflowinflow2O2 H2 H2 O O H OH ArHOHOH2 OH2H2 O2H2O O2OH2HOHOArArOH OH O2OH HH2 OHH2 O H2 OHH2 ArO2 Ar 11000001 00111010 00011100 11100102 11111020 11111200Massflow [kg/s]1.5 0000000010.50 0.5 1 1.5 2 2.5012Time [s]345 3x 10Table 1: Reactions included in the H2 O2 reaction sys- Figure 3: Mass flows into and out of the control voltem and the corresponding stoichiometric matrix.ume during the ignition phaseThe Modelica Association33Modelica 2002, March 18 19, 2002

Chemical Reaction Modeling with ThermoFluid/MF and MultiFlash1Molefraction [1]mole fraction of O20.5mole fraction of H20mole fraction of H200123Time [s]45 3x 10Figure 4: Molar fractions of the principal reactants andproductsinitial ignition, a steady inflow of premixed gases leadsto a steady combustion with plenty of surplus oxygen.The speed of the reactions makes the system very stiff.The whole simulation shown in Figures 4-3 spans onlya few milliseconds.2MultiFlash is the generic name for a physical property software from Infochem Ltd. It is a comprehensive system that calculates the thermophysical properties of pure substances and mixtures and carries outphase and chemical equilibrium calculations for fluidand solid phases. MultiFlash consists of several software modules: databases with the raw property data,access software to the databases and different modules for pure component property calculation, mixture models for thermodynamic properties and transport properties, handling of binary interaction parameters and phase and chemical equilibrium calculations.A number of process simulators, e. g., gPROMS fromPSE, uses an interface to MultiFlash for the calculationof physical properties. For use in a dynamic simulationprogram typically only a small fraction of the MultiFlash functions are needed. The current interface iskept as simple as possible, with all necessary interaction with the property database encapsulated into onemedium property object. The interface has been testedwith both Dymola by Dynasim AB [4] and MathModelica by MathCore AB [6].2.1The low-level interfaceThe low-level interface between Modelica models andthe MultiFlash modules, which are accessible via aModelica 2002, March 18 19, 2002Win32 Dynamic Link Library (dll) under Windows,consists of the standard Modelica foreign function interface for the C-language. This means that all callsto MultiFlash routines are provided exactly as documented in the MultiFlash programmers guide, including identical variable names. There are two minor,but necessary exceptions. Modelica does not allow tooverwrite inputs with outputs in the calling of functions. This is common practice in Fortran numericalroutines and in MultiFlash this is exclusively used toprovide estimates of the solutions in the input variables. In the Modelica interface, the estimated solutions are provided as additional input arguments to thefunction and the original MultiFlash variables are keptas outputs. The second exception is the handling of error message strings. The handling of errors and warnings is done in the C wrapper functions. Diagnosticmessages are written to the simulation log. A flag withthe number of errors is returned in the Modelica function call for error trapping purposes.2.2The MultiFlash interface34Tummescheit H., Eborn J.Computational efficiencySimulation time is an important issue and the interface library uses all available methods to make function calls computationally efficient. A simple rule isto get as many physical properties as possible fromone call to MultiFlash. All essential medium properties needed for the default dynamic model are available in one property record which is calculated withone single function call to MultiFlash. The dynamicstate model and other ThermoFluid models need everything in this standard property record, so it is computationally not efficient to slim down this functioncall. Boolean flags to the MultiFlash routines are usedto ensure that only the medium properties of interestare calculated. High level functions that return onlya single property have not been implemented in orderto close the door on unnecessarily slow models. However, all low-level MultiFlash functions are availableand thus single function calls to obtain properties canbe used if it is desired.Providing good estimates of the solution makes a bigdifference in the solution time for any nonlinear system of equations, especially for phase equilibrium calculations. It is obvious that for continuous, dynamicsimulation the result from the last time step usuallyprovides such an estimate. Internal caching of the lastsolution in the same control volume is therefore implemented in the T HERMO F LUID /MF MultiFlash interface.The Modelica Association

Tummescheit H., Eborn J.2.3Chemical Reaction Modeling with ThermoFluid/MF and MultiFlashUsing MultiFlash with ThermoFluid/MF temperature T and mole amounts N as dynamic states,Setting up a model that uses MultiFlash properties iscurrently done through the MultiFlash windows userinterface. The user selects the wanted components byquerying the available MultiFlash databases. For complex systems, the MultiFlash stream facility is used todefine different component streams to be used in different parts of the system. The global stream whichhas to be defined first contains the union of the components in all streams. If necessary, the thermodynamic models can be changed from the default suggested by MultiFlash through the graphical user interface. All standard equations-of-state (EOS) models,e. g., Redlich-Kwong-Soave or Peng-Robinson, can beused with a selection of mixing rules. Binary interaction parameters can be entered if needed. Transport property models are selected independently of theEOS models. The problem setup is saved in an mflfile which is then read and parsed during initialization.The file name is the main parameter to the mediummodels in the T HERMO F LUID /MF library. Problemsetup files are read from the current directory or froma repository, where all problem setup definitions canbe managed in a centralized way.2.4while p can be regarded as an algebraic variable thatcontains state information. For the standard case of aconstant volume control volume, the pressure is solvedfor iteratively to ensure that the total volume is keptconstant, V f ixed V p . This is the only non-linearequation system in the model for single phase calculations. The dynamic state model is derived from thestandard text-book form of an energy and mass balances. In block matrix notation, the inner energy andmole amount balance can be recast into temperatureand moles as dynamic states as follows (boldface forvectors and matrices, sizes follow from dimension ofN, the number of components in the mixture.): NtUtVt I dUdN T V350dUdT N V0 The most efficient method of combining the dynamicstates and the physical property calculation is tochoose the dynamic states of the model such that theyare inputs to the physical property calculation routines.That avoids the solution of non-linear equation systems during simulation. Otherwise, inputs to the property functions have to be computed from outputs ofthat functions through a non-linear equation system.This happens when the outputs are dynamic states ortime-invariant parameters, like the volume in a closedvessel. If the property functions are computationallyexpensive relative to the rest of the model, the savingin computation time by using a model which is explicit in the states is significant. When this can not beachieved, as is the case with the MultiFlash routines, itis still preferable to get non-linear equation systems ofthe lowest possible dimension. Due to the MF calling structure with pressure p, temperature T and N(mole amounts) as inputs and total volume V as output a special state model has been defined. It is incorporated in the free T HERMO F LUID library and can beused interchangeably with the simple ideal gas modelsin the free T HERMO F LUID library and the commercial MultiFlash property models. The state model uses0 NtTtVt dUdV N T01 (4)The subscript t is used for the time derivative, N standsfor mole amounts, U for total inner energy. The inverse of the jacobian is used to make this model explicit in the mole vector, the temperature and the volume as dynamic states:Dynamic state equationsThe Modelica Association I 11 dUdT N VdUdN T V0 0dUdT N V0 dUdV N T 01(5)The structure of the jacobian inverse reveals that onlythe equation for the inner energy is transformed intoone for the temperature. The mole balance equationsremain unchanged from (4).The partial derivatives occuring in the transformed dynamic and initial equations can be calculated fromderivatives that are returned by MultiFlash by settingthe appropriate flags. However, the derivatives must betransformed using thermodynamic determinants sinceMultiFlash returns derivatives at constant pressure andthe T HERMO F LUID balances are derived at constantvolume. As an example we pick the derivative of total enthalpy w. r. t. temperature (all derivatives are atconstant composition): H T V H T p H V T H T H p T V p T p (6) All derivatives on the right hand side of the equationare returned by standard MultiFlash property calls.Modelica 2002, March 18 19, 2002

Chemical Reaction Modeling with ThermoFluid/MF and MultiFlash2.5InitializationIn the Modelica 2.0 specification and Dymola version4.2 the possibility to separate the equations for steadystate initialization from the dynamic states was introduced. Due to that separation, a control volume cannow easily be initialized at steady state pressure, evenwhen the pressure is not a dynamic state. In complexflow sheets, the calculation of an initial steady state isusually numerically much more challenging than thesubsequent dynamic simulation. A helpful way aroundthat problem is to use a suitable “pseudo steady state”instead which avoids harsh initial transients. One possibility to do so is to use steady state initialization onlyfor the states with relatively fast eigenvalues. Settingthe pressure gradient to zero, but supplying initial estimates for temperature and composition is one suchsuitable choice of a “pseudo steady state”. The fastmodes of the system (pressure and, if a dynamic momentum balance is used, mass flows) are initialized insteady state, while the much slower modes of temperature and composition are set by a non-steady state initial guess. The pressure gradient becomes:dpdtdpdT dTN dtn i 1 dpdNi TdNidt0(7)This equation together with given initial compositionand temperature is much easier to solve than a fullsteady state, especially for large networks, but the initial transients due to errors in the initial guesses areorders of magnitude smaller than the ones obtainedfrom non steady state pressures. The new initializationmethod has been implemented for all state models inthe T HERMO F LUID and T HERMO F LUID /MF librariesand has improved the handling of model initialization considerably. Before implementation of the improved initialization, computation time for small problems was dominated by the time to simulate past theinitial transients. With that obstacle removed, typical simulation times for small systems are an order ofmagnitude faster than before.This initialization is the default setup when the T HER MO F LUID /MF high-level models are used. The initialstate is defined by given temperature and mass fractions and an initial pressure estimate. The initialization then solves for the mole amount states.2.6Tummescheit H., Eborn J.identifier is passed to the wrapper functions calling theMultiFlash property routines. Using the unique control volume identifier, it is possible to connect errorand warning messages from the MultiFlash routines tothe location in the flow sheet where the error occured,e. g., if a temperature rises above the range of validityof the property function. All error and warning messages from MultiFlash are written to the Dymola simulation log. Information about the version, the configuration, the number and composition of streams etc. isalso included in the log.2.7ThermoFluid/MF high level modelsModeling of process engineering problems can notbe cast into fixed, unchangeable model library components as for example multibody systems. Insteadflexibility is needed to have basic building blocks taking care of the standard parts of any dynamic model.These basic models need to be easy to adapt to a specific problem. A large part of the physical propertycalculations is identical for all modeling problems.The T HERMO F LUID /MF library provides such basicmodels and building blocks for control volume modelsbased on MultiFlash properties. Extensions are simple to add by using elements of the T HERMO F LUIDor T HERMO F LUID /MF libraries. Some examples oflumped and distributed models demonstrate how tobuild components and larger systems from the building blocks in the library. A mixture which is typical forfuel cell reformer systems is used to demonstrate howthe minmimal physical property model is used and alsohow to add transport properties. Transport propertiesare not included in the T HERMO F LUID library exceptfor water, but MultiFlash includes several models forviscosity, thermal conductivity and surface tension forpure components and mixtures.DebuggingIn order to improve feedback and error messages fordebugging, an identifier for each control volume isallocated during the initialization of the model. TheModelica 2002, March 18 19, 200236Figure 5: Example models from the library.The Modelica Association

Tummescheit H., Eborn J.Chemical Reaction Modeling with ThermoFluid/MF and MultiFlashThe number of high level models in the T HER MO F LUID /MF library is fairly small, because moststandard models can be used from the T HERMO F LUIDlibrary. The only base models that are differentare control volume models. For flash models newflash control volumes are introduced. They determinewhich phase is flowing in or out at a connector fromthe position of the flow connector and the liquid levelin the volume. Tray models for destillation columnswill be added later.Temperature 3Use that property model in a suitble control volume model from the T HERMO F LUID /MF library.Crossing Functions for multicomponent multi-phase mixturesIn Modelica, crossing functions are usually automatically generated from all equations that contain statements which indicate that a function f x is discontinuous at a certain point x0 . For example, in the following equation:370365 360phase if h hliq or h hvap or p pcrit then 1 else 2;355350Liquid Phase FractionDefine a T HERMO F LUID-compatible propertymodel, following the examples in the T HER MO F LUID /MF library. 2468Time2468Time10120. crossing functions are introduced to monitor thestates of the boolean conditions. This is necessary because numerical integration routines assume continuity of their right-hand side functions. This assumptionis violated in most if-clauses. This can not be automated for external functions that are discontinuous ata point x0 . Thus crossing functions have to be provided by the user in order to make the simulator detect the discontinuity. These crossing functions have tobe consistent with the actual discontinuities, otherwisethey will not work. In the context of phase equilibrium calculations for multi-component fluid mixturesthis means that a unique function of composition, pressure and temperature N p T must be returned fromthe phase equilibrium calculations which has a signchange at the point where a new phase is formed or onephase ceases to exist. At a phase boundary thermodynamic variables have discontinuous first derivativesor are discontinuous by itself, like the heat capacityat constant volume cv . For a mixture with n compoℜ.nents the crossing function is a function ℜ n 1It calculates a measure for the distance to the phaseboundary surface which is in ℜn . Figure 6: Change of temperature and liquid phase fraction in a water-ehtanol mix during a pressure transient. The simulation result from the depressurization ofa flash volume filled with a water-ethanol mixturein thermodynamic equilibrium of the two phases isshown in Figure 6. A ramp from 1 bar to 0 5 barsis imposed on the volume which has a feed flow ofconstant composition. The jump in the liquid phasefraction at the start and end of the transient is due tothe changing in- and outflow phase fractions.Using MultiFlash properties with the T HER MO F LUID /MF library is very simple and requiresonly few steps of setup: 3.1 Deviation indexCollaboration with Infochem Ltd. [5] brought forthan implementation of such a function to increase efficiency and reliability of phase equilibrium calculations in dynamic simulations. It is available in thelatest release of MultiFlash, version 3.1. This is thefirst time that a multi-phase property package has beenDefine the components, phases and models to be equiped with this feature, which is indispensable forused in the MultiFlash user interface and save the being able to reliably simulate the formation or disresult in a model setup file.appearance of phases in a control volume with highThe Modelica Association37Modelica 2002, March 18 19, 2002

Chemical Reaction Modeling with ThermoFluid/MF and MultiFlashquality integrators with event detection and error control. Infochem calls the new function the deviationindex. The calculation of the deviation index is numerically much more efficient than other possibilitiesto determine the number of phases at a given N p Tduring dynamic simulation. Geometrically, the deviation index can be interpreted as a normalized length ofthe normal vector from the current point in the spacespanned by composition, pressure and temperature tothe n-dimensional tangent hyperplane to the phase separation surface. The tangent plane is also known asGibbs' tangent plane. It has been used for stabilityanalysis in phase equilibrium calculations before, see[8], [9] and [1]. The new feature is to assign a value tothe distance from the hyperplane which allows a solverusing interpolation to exactly locate the point in timewhen the simulation trajectory in the N p T - spacewill pass through the hyperplane. At equilibrium, theGibbs energy of the system is at a minimum. Thiscondition may be expressed as the equality of fugacities for each component in all phases or equivalently([1]) as 1 np(8)where nc is the number of components, n p is the number of phases, Ki j is the K-value for component i inphase j, Fi j is the fugacity coefficent for component iin phase j and ri is the the reference phase for component i. The K-values are defined asln Ki j ln Fi j ln FiriKi j0 i yi jyiri1 nc ; j(9)where yi j is the mole fraction of component i in phasej. In the vicinity of the phase split surface, the lefthand side of (8) gives the value of the desired crossingfunction, the deviation index.Furthermore, the function needs to reliably calculatethe properties used in equation (8) of a phase whichis unstable at the current N p T . Considering forsimplicity single component mixtures and the calculation of thermodynamic properties for a phase whichis unstable inside the 2-phase dome (i. e., superheatedliquid or subcooled vapour), it becomes clear that thenumerical computation is only possible to the limit ofthe so called spinoidal lines. For simple cubic EOSthe spinoidal lines are defined by the connection linesof the maxima and minima of the theorectical isothermes inside the two-phase dome. An implementationthus has to guard against erroneous results far fromthe phase boundary. Modelica 2002, March 18 19, 2002 384Tummescheit H., Eborn J.ConclusionsThe inclusion of reaction calculations and the interface to the physical property database MultiFlash intothe T HERMO F LUID library opens new possibilities ofmodeling process systems and combustion processeswhich

The speed of the reactions makes the system very stiff. The wholesimulationshownin Figures 4-3 spans only a few milliseconds. 2 The MultiFlash interface MultiFlash is the generic name for a physical prop-erty software from Infochem Ltd. It is a comprehen-sive system that calculates the thermophysical proper-