All Tutorials

What’s the time? (AKA: Hello world!)

What is the time? (AKA Hello world!) is a very simple model. The model consists of a single variable (called ‘time’) and a single differential equation:

d(time)/dt = 1

This way, the value of the variable ‘time’ is equal to the elapsed time in the simulation at any moment.

This tutorial presents the basic structure of daeModel and daeSimulation classes. A typical DAETools simulation requires the following 8 tasks:

  1. Importing DAE Tools pyDAE module(s)
  2. Importing or declaration of units and variable types (unit and daeVariableType classes)
  3. Developing a model by deriving a class from the base daeModel class and:
    • Declaring domains, parameters and variables in the daeModel.__init__ function
    • Declaring equations and their residual expressions in the daeModel.DeclareEquations function
  4. Setting up a simulation by deriving a class from the base daeSimulation class and:
    • Specifying a model to be used in the simulation in the daeSimulation.__init__ function
    • Setting the values of parameters in the daeSimulation.SetUpParametersAndDomains function
    • Setting initial conditions in the daeSimulation.SetUpVariables function
  5. Declaring auxiliary objects required for the simulation
    • DAE solver object
    • Data reporter object
    • Log object
  6. Setting the run-time options for the simulation:
    • ReportingInterval
    • TimeHorizon
  7. Connecting a data reporter
  8. Initializing, running and finalizing the simulation

The ‘time’ variable plot:

_images/whats_the_time-results.png

Files

Model report whats_the_time.xml
Runtime model report whats_the_time-rt.xml
Source code whats_the_time.py

Tutorial 1

This tutorial introduces several new concepts:

  • Distribution domains
  • Distributed parameters, variables and equations
  • Setting boundary conditions (Neumann and Dirichlet type)
  • Setting initial conditions

In this example we model a simple heat conduction problem: a conduction through a very thin, rectangular copper plate.

For this problem, we need a two-dimensional Cartesian grid (x,y) (here, for simplicity, divided into 10 x 10 segments):

 y axis
    ^
    |
Ly -| L T T T T T T T T T R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
    | L + + + + + + + + + R
 0 -| L B B B B B B B B B R
    --|-------------------|-------> x axis
      0                   Lx

Points ‘B’ at the bottom edge of the plate (for y = 0), and the points ‘T’ at the top edge of the plate (for y = Ly) represent the points where the heat is applied.

The plate is considered insulated at the left (x = 0) and the right edges (x = Lx) of the plate (points ‘L’ and ‘R’). To model this type of problem, we have to write a heat balance equation for all interior points except the left, right, top and bottom edges, where we need to specify boundary conditions.

In this problem we have to define two distribution domains:

  • x (x axis, length Lx = 0.1 m)
  • y (y axis, length Ly = 0.1 m)

the following parameters:

  • rho: copper density, 8960 kg/m3
  • cp: copper specific heat capacity, 385 J/(kgK)
  • k: copper heat conductivity, 401 W/(mK)
  • Qb: heat flux at the bottom edge, 1E6 W/m2 (or 100 W/cm2)
  • Tt: temperature at the top edge, 300 K

and a single variable:

  • T: the temperature of the plate (distributed on x and y domains)

The model consists of 5 equations (1 distributed equation + 4 boundary conditions):

  1. Heat balance:

    rho * cp * dT(x,y)/dt = k * [d2T(x,y)/dx2 + d2T(x,y)/dy2];  x in (0, Lx), y in (0, Ly)
    
  2. Neumann boundary conditions at the bottom edge:

    -k * dT(x,y)/dy = Qb;  x in (0, Lx), y = 0
    
  3. Dirichlet boundary conditions at the top edge:

    T(x,y) = Tt;  x in (0, Lx), y = Ly
    
  4. Neumann boundary conditions at the left edge (insulated):

    dT(x,y)/dx = 0;  y in [0, Ly], x = 0
    
  5. Neumann boundary conditions at the right edge (insulated):

    dT(x,y)/dx = 0;  y in [0, Ly], x = Lx
    

The temperature plot (at t=100s, x=0.5, y=*):

_images/tutorial1-results.png

Files

Model report tutorial1.xml
Runtime model report tutorial1-rt.xml
Source code tutorial1.py

Tutorial 2

This tutorial introduces the following concepts:

  • Arrays (discrete distribution domains)
  • Distributed parameters
  • Making equations more readable
  • Degrees of freedom
  • Setting an initial guess for variables (used by a DAE solver during an initial phase)

The model in this example is very similar to the model used in the tutorial 1. The differences are:

  • The heat capacity is temperature dependent
  • Different boundary conditions are applied

The temperature plot (at t=100s, x=0.5, y=*):

_images/tutorial2-results.png

Files

Model report tutorial2.xml
Runtime model report tutorial2-rt.xml
Source code tutorial2.py

Tutorial 3

This tutorial introduces the following concepts:

  • Arrays of variable values
  • Functions that operate on arrays of values
  • Functions that create constants and arrays of constant values (Constant and Array)
  • Non-uniform domain grids

The model in this example is identical to the model used in the tutorial 1. Some additional equations that calculate the total flux at the bottom edge are added to illustrate the array functions.

The temperature plot (at t=100, x=0.5, y=*):

_images/tutorial3-results.png

The average temperature plot (considering the whole x-y domain):

_images/tutorial3-results2.png

Files

Model report tutorial3.xml
Runtime model report tutorial3-rt.xml
Source code tutorial3.py

Tutorial 4

This tutorial introduces the following concepts:

  • Discontinuous equations (symmetrical state transition networks: daeIF statements)

In this example we model a very simple heat transfer problem where a small piece of copper is at one side exposed to the source of heat and at the other to the surroundings.

The lumped heat balance is given by the following equation:

rho * cp * dT/dt - Qin = h * A * (T - Tsurr)

where Qin is the power of the heater, h is the heat transfer coefficient, A is the surface area and Tsurr is the temperature of the surrounding air.

The process starts at the temperature of the metal of 283K. The metal is allowed to warm up for 200 seconds, when the heat source is removed and the metal cools down slowly to the ambient temperature.

This can be modelled using the following symmetrical state transition network:

IF t < 200
  Qin = 1500 W
ELSE
  Qin = 0 W

The temperature plot:

_images/tutorial4-results.png

Files

Model report tutorial4.xml
Runtime model report tutorial4-rt.xml
Source code tutorial4.py

Tutorial 5

This tutorial introduces the following concepts:

  • Discontinuous equations (non-symmetrical state transition networks: daeSTN statements)

In this example we use the same heat transfer problem as in the tutorial 4. Again we have a piece of copper which is at one side exposed to the source of heat and at the other to the surroundings.

The process starts at the temperature of 283K. The metal is allowed to warm up, and then its temperature is kept in the interval (320K - 340K) for 350 seconds. This is performed by switching the heater on when the temperature drops to 320K and by switching the heater off when the temperature reaches 340K. After 350s the heat source is permanently switched off and the metal is allowed to slowly cool down to the ambient temperature.

This can be modelled using the following non-symmetrical state transition network:

STN Regulator
  case Heating:
    Qin = 1500 W
    on condition T > 340K switch to Regulator.Cooling
    on condition t > 350s switch to Regulator.HeaterOff

  case Cooling:
    Qin = 0 W
    on condition T < 320K switch to Regulator.Heating
    on condition t > 350s switch to Regulator.HeaterOff

  case HeaterOff:
    Qin = 0 W

The temperature plot:

_images/tutorial5-results.png

Files

Model report tutorial5.xml
Runtime model report tutorial5-rt.xml
Source code tutorial5.py

Tutorial 6

This tutorial introduces the following concepts:

  • Ports
  • Port connections
  • Units (instances of other models)

A simple port type ‘portSimple’ is defined which contains only one variable ‘t’. Two models ‘modPortIn’ and ‘modPortOut’ are defined, each having one port of type ‘portSimple’. The wrapper model ‘modTutorial’ instantiate these two models as its units and connects them by connecting their ports.

Files

Model report tutorial6.xml
Runtime model report tutorial6-rt.xml
Source code tutorial6.py

Tutorial 7

This tutorial introduces the following concepts:

  • User-defined schedules (operating procedures)
  • Resetting of degrees of freedom
  • Resetting of initial conditions

In this example we use the same heat transfer problem as in the tutorial 4. The input power of the heater is defined as a variable. Since there is no equation defined to calculate the value of the input power, the system contains N variables but only N-1 equations. To create a well-posed DAE system one of the variable needs to be “fixed”. However the choice of variables is not arbitrary and in this example the only variable that can be fixed is Qin. Thus, the Qin variable represents a degree of freedom (DOF). Its value will be fixed at the beginning of the simulation and later manipulated in the user-defined schedule in the overloaded function daeSimulation.Run().

The plot of the inlet power:

_images/tutorial7-results.png

The temperature plot:

_images/tutorial7-results2.png

Files

Model report tutorial7.xml
Runtime model report tutorial7-rt.xml
Source code tutorial7.py

Tutorial 8

This tutorial introduces the following concepts:

  • Data reporters and exporting results into the following file formats:
    • Matlab MAT file (requires python-scipy package)
    • MS Excel .xls file (requires python-xlwt package)
    • JSON file (no third party dependencies)
    • VTK file (requires pyevtk and vtk packages)
    • XML file (requires python-lxml package)
    • HDF5 file (requires python-h5py package)
    • Pandas dataset (requires python-pandas package)
  • Implementation of user-defined data reporters
  • daeDelegateDataReporter

Some time it is not enough to send the results to the DAE Plotter but it is desirable to export them into a specified file format (i.e. for use in other programs). For that purpose, daetools provide a range of data reporters that save the simulation results in various formats. In adddition, daetools allow implementation of custom, user-defined data reporters. As an example, a user-defined data reporter is developed to save the results into a plain text file (after the simulation is finished). Obviously, the data can be processed in any other fashion. Moreover, in certain situation it is required to process the results in more than one way. The daeDelegateDataReporter can be used in those cases. It has the same interface and the functionality like all data reporters. However, it does not do any data processing itself but calls the corresponding functions of data reporters which are added to it using the function AddDataReporter. This way it is possible, at the same time, to send the results to the DAE Plotter and save them into a file (or process the data in some other ways). In this example the results are processed in 10 different ways at the same time.

The model used in this example is very similar to the model in the tutorials 4 and 5.

Files

Model report tutorial8.xml
Runtime model report tutorial8-rt.xml
Source code tutorial8.py

Tutorial 9

This tutorial introduces the following concepts:

  • Third party direct linear equations solvers

Currently there are the following linear equations solvers available:

  • SuperLU: sequential sparse direct solver defined in pySuperLU module (BSD licence)
  • SuperLU_MT: multi-threaded sparse direct solver defined in pySuperLU_MT module (BSD licence)
  • Trilinos Amesos: sequential sparse direct solver defined in pyTrilinos module (GNU Lesser GPL)
  • IntelPardiso: multi-threaded sparse direct solver defined in pyIntelPardiso module (proprietary)
  • Pardiso: multi-threaded sparse direct solver defined in pyPardiso module (proprietary)

In this example we use the same conduction problem as in the tutorial 1.

The temperature plot (at t=100s, x=0.5, y=*):

_images/tutorial9-results.png

Files

Model report tutorial9.xml
Runtime model report tutorial9-rt.xml
Source code tutorial9.py

Tutorial 10

This tutorial introduces the following concepts:

  • Initialization files
  • Domains which bounds depend on parameter values
  • Evaluation of integrals

In this example we use the same conduction problem as in the tutorial 1.

Files

Model report tutorial10.xml
Runtime model report tutorial10-rt.xml
Source code tutorial10.py

Tutorial 11

This tutorial describes the use of iterative linear solvers (AztecOO from the Trilinos project) with different preconditioners (built-in AztecOO, Ifpack or ML) and corresponding solver options. Also, the range of Trilins Amesos solver options are shown.

The model is very similar to the model in tutorial 1, except for the different boundary conditions and that the equations are written in a different way to maximise the number of items around the diagonal (creating the problem with the diagonally dominant matrix). These type of systems can be solved using very simple preconditioners such as Jacobi. To do so, the interoperability with the NumPy package has been exploited and the package itertools used to iterate through the distribution domains in x and y directions.

The equations are distributed in such a way that the following incidence matrix is obtained:

|XXX                                 |
| X     X     X                      |
|  X     X     X                     |
|   X     X     X                    |
|    X     X     X                   |
|   XXX                              |
|      XXX                           |
| X    XXX    X                      |
|  X    XXX    X                     |
|   X    XXX    X                    |
|    X    XXX    X                   |
|         XXX                        |
|            XXX                     |
|       X    XXX    X                |
|        X    XXX    X               |
|         X    XXX    X              |
|          X    XXX    X             |
|               XXX                  |
|                  XXX               |
|             X    XXX    X          |
|              X    XXX    X         |
|               X    XXX    X        |
|                X    XXX    X       |
|                     XXX            |
|                        XXX         |
|                   X    XXX    X    |
|                    X    XXX    X   |
|                     X    XXX    X  |
|                      X    XXX    X |
|                           XXX      |
|                              XXX   |
|                   X     X     X    |
|                    X     X     X   |
|                     X     X     X  |
|                      X     X     X |
|                                 XXX|

The temperature plot (at t=100s, x=0.5, y=*):

_images/tutorial11-results.png

Files

Model report tutorial11.xml
Runtime model report tutorial11-rt.xml
Source code tutorial11.py

Tutorial 12

This tutorial describes the use and available options for superLU direct linear solvers:

  • Sequential: superLU
  • Multithreaded (OpenMP/posix threads): superLU_MT

The model is the same as the model in tutorial 1, except for the different boundary conditions.

The temperature plot (at t=100s, x=0.5, y=*):

_images/tutorial12-results.png

Files

Model report tutorial12.xml
Runtime model report tutorial12-rt.xml
Source code tutorial12.py

Tutorial 13

This tutorial introduces the following concepts:

  • The event ports
  • ON_CONDITION() function illustrating the types of actions that can be executed during state transitions
  • ON_EVENT() function illustrating the types of actions that can be executed when an event is triggered
  • User defined actions

In this example we use the very similar model as in the tutorial 5.

The simulation output should show the following messages at t=100s and t=350s:

...
********************************************************
simpleUserAction2 message:
This message should be fired when the time is 100s.
********************************************************
...
********************************************************
simpleUserAction executed; input data = 427.464093129832
********************************************************
...

The plot of the ‘event’ variable:

_images/tutorial13-results.png

The temperature plot:

_images/tutorial13-results2.png

Files

Model report tutorial13.xml
Runtime model report tutorial13-rt.xml
Source code tutorial13.py

Tutorial 14

In this tutorial we introduce the external functions concept that can handle and execute functions in external libraries. The daeScalarExternalFunction-derived external function object is used to calculate the heat transferred and to interpolate a set of values using the scipy.interpolate.interp1d object.

In this example we use the same model as in the tutorial 5 with few additional equations.

The simulation output should show the following messages at the end of simulation:

...
scipy.interp1d statistics:
  interp1d called 1703 times (cache value used 770 times)

The plot of the ‘Heat_ext’ variable:

_images/tutorial14-results.png

The plot of the ‘Value_interp’ variable:

_images/tutorial14-results2.png

Files

Model report tutorial14.xml
Runtime model report tutorial14-rt.xml
Source code tutorial14.py

Tutorial 15

This tutorial introduces the following concepts:

  • Nested state transitions

In this example we use the same model as in the tutorial 4 with the more complex STN:

IF t < 200
  IF 0 <= t < 100
     IF 0 <= t < 50
       Qin = 1600 W
     ELSE
       Qin = 1500 W
  ELSE
    Qin = 1400 W

ELSE IF 200 <= t < 300
  Qin = 1300 W

ELSE
  Qin = 0 W

The plot of the ‘Qin’ variable:

_images/tutorial15-results.png

The temperature plot:

_images/tutorial15-results2.png

Files

Model report tutorial15.xml
Runtime model report tutorial15-rt.xml
Source code tutorial15.py

Tutorial 16

This tutorial shows how to use DAE Tools objects with NumPy arrays to solve a simple stationary heat conduction in one dimension using the Finite Elements method with linear elements and two ways of manually assembling a stiffness matrix/load vector:

d2T(x)/dx2 = F(x);  x in (0, Lx)

Linear finite elements discretisation and simple FE matrix assembly:

                  phi                 phi
                     (k-1)               (k)
                     
                    *                   *             
                  * | *               * | *            
                *   |   *           *   |   *          
              *     |     *       *     |     *        
            *       |       *   *       |       *      
          *         |         *         |         *        
        *           |       *   *       |           *      
      *             |     *       *     |             *    
    *               |   *           *   |               *  
  *                 | *  element (k)  * |                 *
*-------------------*+++++++++++++++++++*-------------------*-
                    x                   x
                     (k-i                (k)
                   
                    \_________ _________/
                              |
                              dx

The comparison of the analytical solution and two ways of assembling the system is given in the following plot:

_images/tutorial16-results.png

Files

Model report tutorial16.xml
Runtime model report tutorial16-rt.xml
Source code tutorial16.py

Tutorial 17

This tutorial introduces the following concepts:

  • TCPIP Log and TCPIPLogServer

In this example we use the same heat transfer problem as in the tutorial 7.

The screenshot of the TCP/IP log server:

_images/tutorial17-screenshot.png

The temperature plot:

_images/tutorial17-results.png

Files

Model report tutorial17.xml
Runtime model report tutorial17-rt.xml
Source code tutorial17.py

Tutorial 18

This tutorial shows one more problem solved using the NumPy arrays that operate on DAE Tools variables. The model is taken from the Sundials ARKODE (ark_analytic_sys.cpp). The ODE system is defined by the following system of equations:

dy/dt = A*y

where:

A = V * D * Vi
V = [1 -1 1; -1 2 1; 0 -1 2];
Vi = 0.25 * [5 1 -3; 2 2 -2; 1 1 1];
D = [-0.5 0 0; 0 -0.1 0; 0 0 lam];

lam is a large negative number.

The analytical solution to this problem is:

Y(t) = V*exp(D*t)*Vi*Y0

for t in the interval [0.0, 0.05], with initial condition y(0) = [1,1,1]’.

The stiffness of the problem is directly proportional to the value of “lamda”. The value of lamda should be negative to result in a well-posed ODE; for values with magnitude larger than 100 the problem becomes quite stiff.

In this example, we choose lamda = -100.

The solution:

lamda = -100
reltol = 1e-06
abstol = 1e-10

--------------------------------------
   t        y0        y1        y2
--------------------------------------
 0.0050   0.70327   0.70627   0.41004
 0.0100   0.52267   0.52865   0.05231
 0.0150   0.41249   0.42145  -0.16456
 0.0200   0.34504   0.35696  -0.29600
 0.0250   0.30349   0.31838  -0.37563
 0.0300   0.27767   0.29551  -0.42383
 0.0350   0.26138   0.28216  -0.45296
 0.0400   0.25088   0.27459  -0.47053
 0.0450   0.24389   0.27053  -0.48109
 0.0500   0.23903   0.26858  -0.48740
--------------------------------------

The plot of the ‘y0’, ‘y1’, ‘y2’ variables:

_images/tutorial18-results.png

Files

Model report tutorial18.xml
Runtime model report tutorial18-rt.xml
Source code tutorial18.py

Tutorial 19

This tutorial introduces the thermo physical property packages.

Since there are many thermo packages with a very different API the CapeOpen standard has been adopted in daetools. This way, all thermo packages implementing the CapeOpen thermo interfaces are automatically vailable to daetools. Those which do not are wrapped by the class with the CapeOpen conforming API. At the moment, two types of thermophysical property packages are implemented:

  • Any CapeOpen v1.1 thermo package (available only in Windows)
  • CoolProp thermo package (available for all platforms) wrapped in the class with the CapeOpen interface.

The central point is the daeThermoPhysicalPropertyPackage class. It can load any COM component that implements CapeOpen 1.1 ICapeThermoPropertyPackageManager interface or the CoolProp thermo package.

The framework provides low-level functions (specified in the CapeOpen standard) in the daeThermoPhysicalPropertyPackage class and the higher-level functions in the auxiliary daeThermoPackage class defined in the daetools/pyDAE/thermo_packages.py file. The low-level functions are defined in the ICapeThermoCoumpounds and ICapeThermoPropertyRoutine CapeOpen interfaces. These functions come in two flavours:

  1. The ordinary functions return adouble/adouble_array objects and can only be used to specify equations:
    • GetCompoundConstant (from ICapeThermoCoumpounds interface)
    • GetTDependentProperty (from ICapeThermoCoumpounds interface)
    • GetPDependentProperty (from ICapeThermoCoumpounds interface)
    • CalcSinglePhaseScalarProperty (from ICapeThermoPropertyRoutine interface: scalar version)
    • CalcSinglePhaseVectorProperty (from ICapeThermoPropertyRoutine interface: vector version)
    • CalcTwoPhaseScalarProperty (from ICapeThermoPropertyRoutine interface: scalar version)
    • CalcTwoPhaseVectorProperty (from ICapeThermoPropertyRoutine interface: vector version)
  2. The functions starting with the underscores can be used for calculations (they use and return float values):
    • _GetCompoundConstant
    • _GetTDependentProperty
    • _GetPDependentProperty
    • _CalcSinglePhaseScalarProperty
    • _CalcSinglePhaseVectorProperty
    • _CalcTwoPhaseScalarProperty
    • _CalcTwoPhaseVectorProperty

The daeThermoPackage auxiliary class offers functions to calculate specified properties, for instance:

  • Transport properties:
    • cp, kappa, mu, Dab (heat capacity, thermal conductivity, dynamic viscosity, diffusion coefficient)
  • Thermodynamic properties:
    • rho
    • h, s, G, H, I (enthalpy, entropy, gibbs/helmholtz/internal energy)
    • h_E, s_E, G_E, H_E, I_E, V_E (excess enthalpy, entropy, gibbs/helmholtz/internal energy, volume)
    • f and phi (fugacity and coefficient of fugacity)
    • a and gamma (activity and the coefficient of activity)
    • z (compressibility factor)
    • K, surfaceTension (ratio of fugacity coefficients and the surface tension)
Nota bene:
Some of the above functions return scalars while the others return arrays of values. Check the thermo_packages.py file for details.

All functions return properties in the SI units (as specified in the CapeOpen 1.1 standard).

Known issues:

  • Many properties from the CapeOpen standard are not supported by all thermo packages.
  • CalcEquilibrium from the ICapeThermoEquilibriumRoutine is not supported.
  • CoolProp does not provide transport models for many compounds.
  • The function calls are NOT thread safe.
  • The code generation will NOT work for models using the thermo packages.
  • Some CapeOpen thermo packags refuse to return properties for mass basis (i.e. density).

In this tutorial, we use a very simple model: a quantity of liquid (water + ethanol mixture) is heated using the constant input power. The model uses a thermo package to calculate the commonly used transport properties such as specific heat capacity, thermal conductivity, dynamic viscosity and binary diffusion coefficients. First, the low-level functions are tested for CapeOpen and CoolProp packages in the test_single_phase, test_two_phase, test_coolprop_single_phase functions. The results depend on the options selected in the CapeOpen package (equation of state, etc.). Then, the model that uses a thermo package is simulated.

The plot of the specific heat capacity as a function of temperature:

_images/tutorial19-results.png
Nota bene:
There is a difference between results in Windows and other platforms since the CapeOpen thermo packages are available only in Windows.

Files

Model report tutorial19.xml
Runtime model report tutorial19-rt.xml
Source code tutorial19.py

Tutorial 20

This tutorial illustrates the sensitivity analysis features in DAE Tools.

This model has one state variable (T) and one degree of freedom (Qin). Qin is set as a parameter for sensitivity analysis.

The sensitivity analysis is enabled and the sensitivities can be reported to the data reporter like any ordinary variable by setting the boolean property simulation.ReportSensitivities to True.

Raw sensitivity matrices can be saved into a specified directory using the simulation.SensitivityDataDirectory property (before a call to Initialize). The sensitivity matrics are saved in .mmx coordinate format where the first dimensions is Nparameters and second Nvariables: S[Np, Nvars].

The plot of the sensitivity of T per Qin:

_images/tutorial20-results.png

Files

Model report tutorial20.xml
Runtime model report tutorial20-rt.xml
Source code tutorial20.py

Tutorial 21

This tutorial illustrates the additional sensitivity analysis features. Here, the numerical sensitivities for the Constant coefficient first order equations are compared to the available analytical solution.

The sensitivity analysis is enabled and the sensitivities are reported to the data reporter. The sensitivity data can be obtained in two ways:

  • Directly from the DAE solver in the user-defined Run function using the DAESolver.SensitivityMatrix property.
  • From the data reporter as any ordinary variable.

The comparison between the numerical and the analytical sensitivities:

_images/tutorial21-results.png

Files

Model report tutorial21.xml
Runtime model report tutorial21-rt.xml
Source code tutorial21.py

Advanced Tutorial 1

This tutorial presents a user-defined simulation which instead of simply integrating the system shows the pyQt graphical user interface (GUI) where the simulation can be manipulated (a sort of interactive operating procedure).

The model in this example is the same as in the tutorial 4.

The simulation.Run() function is modifed to show the graphical user interface (GUI) that allows to specify the input power of the heater (degree of freedom), a time period for integration, and a reporting interval. The GUI also contains the temperature plot updated in real time, as the simulation progresses.

The screenshot of the pyQt GUI:

_images/tutorial_adv_1-screenshot.png

Files

Model report tutorial_adv_1.xml
Runtime model report tutorial_adv_1-rt.xml
Source code tutorial_adv_1.py

Advanced Tutorial 2

This tutorial demonstrates a solution of a discretized population balance using high resolution upwind schemes with flux limiter.

Reference: Qamar S., Elsner M.P., Angelov I.A., Warnecke G., Seidel-Morgenstern A. (2006) A comparative study of high resolution schemes for solving population balances in crystallization. Computers and Chemical Engineering 30(6-7):1119-1131. doi:10.1016/j.compchemeng.2006.02.012

It shows a comparison between the analytical results and various discretization schemes:

  • I order upwind scheme
  • II order central scheme
  • Cell centered finite volume method solved using the high resolution upwind scheme (Van Leer k-interpolation with k = 1/3 and Koren flux limiter)

The problem is from the section 3.1 Size-independent growth.

Could be also found in: Motz S., Mitrović A., Gilles E.-D. (2002) Comparison of numerical methods for the simulation of dispersed phase systems. Chemical Engineering Science 57(20):4329-4344. doi:10.1016/S0009-2509(02)00349-4

The comparison of number density functions between the analytical solution and the high-resolution scheme:

_images/tutorial_adv_2-results.png

The comparison of number density functions between the analytical solution and the I order upwind, II order upwind and II order central difference schemes:

_images/tutorial_adv_2-results2.png

Files

Model report tutorial_adv_2.xml
Runtime model report tutorial_adv_2-rt.xml
Source code tutorial_adv_2.py

Advanced Tutorial 3

This tutorial introduces the following concepts:

  • DAE Tools code-generators
    • Modelica code-generator
    • gPROMS code-generator
    • FMI code-generator (for Co-Simulation)
  • DAE Tools model-exchange capabilities:
    • Scilab/GNU_Octave/Matlab MEX functions
    • Simulink S-functions

The model represent a simple multiplier block. It contains two inlet and two outlet ports. The outlets values are equal to inputs values multiplied by a multiplier “m”:

out1.y   = m1   x in1.y
out2.y[] = m2[] x in2.y[]

where multipliers m1 and m2[] are:

STN Multipliers
   case variableMultipliers:
      dm1/dt   = p1
      dm2[]/dt = p2
   case constantMultipliers:
      dm1/dt   = 0
      dm2[]/dt = 0

(that is the multipliers can be constant or variable).

The ports in1 and out1 are scalar (width = 1). The ports in2 and out2 are vectors (width = 1).

Achtung, Achtung!! Notate bene:

  1. Inlet ports must be DOFs (that is to have their values asssigned), for they can’t be connected when the model is simulated outside of daetools context.
  2. Only scalar output ports are supported at the moment!! (Simulink issue)

The plot of the inlet ‘y’ variable and the multiplied outlet ‘y’ variable for the constant multipliers (m1 = 2):

_images/tutorial_adv_3-results.png

The plot of the inlet ‘y’ variable and the multiplied outlet ‘y’ variable for the variable multipliers (dm1/dt = 10, m1(t=0) = 2):

_images/tutorial_adv_3-results2.png

FMI Cross-Check results:

  1. Start the DAEPlotter (or change the data reporter below).

  2. Execute:

    ./fmuCheck.linux64 -n 10 tutorial_adv_3.fmu
    
  3. Results:

    [INFO][FMUCHK] FMI compliance checker 2.0 [FMILibrary: 2.0] build date: Aug 22 2014
    [INFO][FMUCHK] Will process FMU tutorial_adv_3.fmu
    [INFO][FMILIB] XML specifies FMI standard version 2.0
    [INFO][FMUCHK] Model name: tutorial_adv_3
    [INFO][FMUCHK] Model GUID: e9654532-0998-11e6-957b-9cb70d5dfdfc
    [INFO][FMUCHK] Model version:
    [INFO][FMUCHK] FMU kind: CoSimulation
    [INFO][FMUCHK] The FMU contains:
    0 constants
    3 parameters
    0 discrete variables
    4 continuous variables
    2 inputs
    2 outputs
    0 local variables
    0 independent variables
    0 calculated parameters
    6 real variables
    0 integer variables
    0 enumeration variables
    0 boolean variables
    1 string variables
    
    [INFO][FMUCHK] Printing output file header
    time,out_1.y,out_2.y
    [INFO][FMUCHK] Model identifier for CoSimulation: tutorial_adv_3
    [INFO][FMILIB] Loading 'linux64' binary with 'default' platform types
    [INFO][FMUCHK] Version returned from CS FMU:   2.0
    ***********************************************************************
                                     Version:   1.5.0
                                     Copyright: Dragan Nikolic
                                     Homepage:  http://www.daetools.com
           @                       @
           @   @@@@@     @@@@@   @@@@@    @@@@@    @@@@@   @      @@@@@
      @@@@@@        @   @     @    @     @     @  @     @  @     @
     @     @   @@@@@@   @@@@@@     @     @     @  @     @  @      @@@@@
     @     @  @     @   @          @     @     @  @     @  @           @
      @@@@@@   @@@@@@    @@@@@      @@@   @@@@@    @@@@@    @@@@  @@@@@
    ***********************************************************************
    DAE Tools is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License version 3
    as published by the Free Software Foundation.
    DAE Tools is distributed in the 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 this program. If not, see <http://www.gnu.org/licenses>.
    ***********************************************************************
    Creating the system...
    The system created successfully in: 0.002 s
    Starting the initialization of the system... Done.
    [INFO][FMUCHK] Initialized FMU for simulation starting at time 0
    0,1.0000000010000001E+00,2.0000000039999999E+00
    10,1.0100000000000001E+02,4.0200000000000006E+02
    20,2.0099999999999994E+02,8.0199999999999977E+02
    30,3.0100000000000000E+02,1.2020000000000000E+03
    40,4.0100000000000000E+02,1.6020000000000000E+03
    50,5.0099999999999989E+02,2.0019999999999995E+03
    60,6.0100000000000000E+02,2.4020000000000000E+03
    70,7.0100000000000000E+02,2.8020000000000000E+03
    80,8.0100000000000000E+02,3.2020000000000000E+03
    90,9.0099999999999977E+02,3.6019999999999991E+03
    100,1.0009999999999998E+03,4.0019999999999991E+03
    [INFO][FMUCHK] Simulation finished successfully at time 100
    FMU check summary:
    FMU reported:
            0 warning(s) and error(s)
    Checker reported:
            0 Warning(s)
            0 Error(s)
    

Files

Model report tutorial_adv_3.xml
Runtime model report tutorial_adv_3-rt.xml
Source code tutorial_adv_3.py

Advanced Tutorial 4

This tutorial illustrates the C++ MPI code generator. The model is identical to the model in the tutorial 11.

The temperature plot (at t=100s, x=0.5128, y=*):

_images/tutorial_adv_4-results.png

Files

Model report tutorial_adv_4.xml
Runtime model report tutorial_adv_4-rt.xml
Source code tutorial_adv_4.py

Tutorial deal.II 1

An introductory example of the support for Finite Elements in daetools. The basic idea is to use an external library to perform all low-level tasks such as management of mesh elements, degrees of freedom, matrix assembly, management of boundary conditions etc. deal.II library (www.dealii.org) is employed for these tasks. The mass and stiffness matrices and the load vector assembled in deal.II library are used to generate a set of algebraic/differential equations in the following form: [Mij]{dx/dt} + [Aij]{x} = {Fi}. Specification of additional equations such as surface/volume integrals are also available. The numerical solution of the resulting ODA/DAE system is performed in daetools together with the rest of the model equations.

The unique feature of this approach is a capability to use daetools variables to specify boundary conditions, time varying coefficients and non-linear terms, and evaluate quantities such as surface/volume integrals. This way, the finite element model is fully integrated with the rest of the model and multiple FE systems can be created and coupled together. In addition, non-linear and DAE finite element systems are automatically supported.

In this tutorial the simple transient heat conduction problem is solved using the finite element method:

dT/dt - kappa/(rho*cp)*nabla^2(T) = g(T) in Omega

The mesh is rectangular with two holes, similar to the mesh in step-49 deal.II example:

Dirichlet boundary conditions are set to 300 K on the outer rectangle, 350 K on the inner ellipse and 250 K on the inner diamond.

The temperature plot at t = 500s (generated in VisIt):

_images/tutorial_dealii_1-results.png

Files

Model report tutorial_dealii_1.xml
Runtime model report tutorial_dealii_1-rt.xml
Source code tutorial_dealii_1.py

Tutorial deal.II 2

In this example a simple transient heat convection-diffusion equation is solved.

dT/dt - kappa/(rho*cp)*nabla^2(T) + nabla.(uT) = g(T) in Omega

The fluid flows from the left side to the right with constant velocity of 0.01 m/s. The inlet temperature for 0.2 <= y <= 0.3 is iven by the following expression:

T_left = T_base + T_offset*|sin(pi*t/25)| on dOmega

creating a bubble-like regions of higher temperature that flow towards the right end and slowly diffuse into the bulk flow of the fluid due to the heat conduction.

The mesh is rectangular with the refined elements close to the left/right ends:

_images/rect(1.5,0.5)-100x50.png

The temperature plot at t = 500s:

_images/tutorial_dealii_2-results.png

Files

Model report tutorial_dealii_2.xml
Runtime model report tutorial_dealii_2-rt.xml
Source code tutorial_dealii_2.py

Tutorial deal.II 3

In this example the Cahn-Hilliard equation is solved using the finite element method. This equation describes the process of phase separation, where two components of a binary mixture separate and form domains pure in each component.

dc/dt - D*nabla^2(mu) = 0, in Omega
mu = c^3 - c - gamma*nabla^2(c)

The mesh is a simple square (0-100)x(0-100):

_images/square.png

The concentration plot at t = 500s:

_images/tutorial_dealii_3-results.png

Files

Model report tutorial_dealii_3.xml
Runtime model report tutorial_dealii_3-rt.xml
Source code tutorial_dealii_3.py

Tutorial deal.II 4

In this tutorial the transient heat conduction problem is solved using the finite element method:

dT/dt - kappa * nabla^2(Τ) = g in Omega

The mesh is a pipe submerged into water receiving the heat of the sun at the 45 degrees from the top-left direction, the heat from the suroundings and having the constant temperature of the inner tube. The boundary conditions are thus:

  • [at boundary ID=0] Sun shine at 45 degrees, gradient heat flux = 2 kW/m**2 in direction n = (1,-1)
  • [at boundary ID=1] Outer surface where y <= -0.5, constant flux of 2 kW/m**2
  • [at boundary ID=2] Inner tube: constant temperature of 300 K

The pipe mesh is given below:

_images/pipe.png

The temperature plot at t = 3600s:

_images/tutorial_dealii_4-results.png

Files

Model report tutorial_dealii_4.xml
Runtime model report tutorial_dealii_4-rt.xml
Source code tutorial_dealii_4.py

Tutorial deal.II 5

In this example a simple flow through porous media is solved (deal.II step-20).

K^{-1} u + nabla(p) = 0 in Omega
-div(u) = -f in Omega
p = g on dOmega

The mesh is a simple square:

_images/square(-1,1)x(-1,1)-50x50.png

The velocity plot at t = 500s:

_images/tutorial_dealii_5-results.png

Files

Model report tutorial_dealii_5.xml
Runtime model report tutorial_dealii_5-rt.xml
Source code tutorial_dealii_5.py

Tutorial deal.II 6

A simple steady-state diffusion and first-order reaction in an irregular catalyst shape (Proc. 6th Int. Conf. on Mathematical Modelling, Math. Comput. Modelling, Vol. 11, 375-319, 1988) applying Dirichlet and Robin type of boundary conditions.

D_eA * nabla^2(C_A) - k_r * C_A = 0 in Omega
D_eA * nabla(C_A) = k_m * (C_A - C_Ab) on dOmega1
C_A = C_Ab on dOmega2

The catalyst pellet mesh:

_images/ssdr.png

The concentration plot:

_images/tutorial_dealii_6-results1.png

The concentration plot for Ca=Cab on all boundaries:

_images/tutorial_dealii_6-results2.png

Files

Model report tutorial_dealii_6.xml
Runtime model report tutorial_dealii_6-rt.xml
Source code tutorial_dealii_6.py

Chem. Eng. Example 1

Continuously Stirred Tank Reactor with energy balance and Van de Vusse reactions:

 A -> B -> C
2A -> D

Reference: G.A. Ridlehoover, R.C. Seagrave. Optimization of Van de Vusse Reaction Kinetics Using Semibatch Reactor Operation, Ind. Eng. Chem. Fundamen. 1973;12(4):444-447. doi:10.1021/i160048a700

The concentrations plot:

_images/tutorial_che_1-results.png

The temperatures plot:

_images/tutorial_che_1-results2.png

Files

Model report tutorial_che_1.xml
Runtime model report tutorial_che_1-rt.xml
Source code tutorial_che_1.py

Chem. Eng. Example 2

Binary distillation column model.

Reference: J. Hahn, T.F. Edgar. An improved method for nonlinear model reduction using balancing of empirical gramians. Computers and Chemical Engineering 2002; 26:1379-1397. doi:10.1016/S0098-1354(02)00120-5

The liquid fraction after 120 min (x(reboiler)=0.935420, x(condenser)=0.064581):

_images/tutorial_che_2-results.png

The liquid fraction in the reboiler (tray 1) and in the condenser (tray 32):

_images/tutorial_che_2-results2.png

Files

Model report tutorial_che_2.xml
Runtime model report tutorial_che_2-rt.xml
Source code tutorial_che_2.py

Chem. Eng. Example 3

Batch reactor seeded crystallisation using the method of moments.

References (model equations and input parameters):

  • Nikolic D.D., Frawley P.J. (2016) Application of the Lagrangian Meshfree Approach to Modelling of Batch Crystallisation: Part I – Modelling of Stirred Tank Hydrodynamics. Chemical Engineering Science 145:317–328. doi:10.1016/j.ces.2015.08.052
  • Mitchell N.A., O’Ciardha C.T., Frawley P.J. (2011) Estimation of the growth kinetics for the cooling crystallisation of paracetamol and ethanol solutions. Journal of Crystal Growth 328:39–49. doi:10.1016/j.jcrysgro.2011.06.016

The main assumptions:

  • Seeded crystallisation
  • Ideal mixing
  • Fixed cooling rate
  • Size independent growth

Solubility of Paracetamol in ethanol:

---------------------------------------------------------------------------------
Temperature, C   Solubility, kg Parac./kg EtOH    Solubility, mol Parac./m3 EtOH
---------------------------------------------------------------------------------
0                0.11362                          593.0387
10               0.14128                          737.4215
20               0.17568                          916.9562
30               0.21845                          1140.2008
40               0.27163                          1417.7972
50               0.33777                          1762.9779
60               0.42000                          2192.1973
---------------------------------------------------------------------------------

The supersaturation plot:

_images/tutorial_che_3-results.png

The concentration plot:

_images/tutorial_che_3-results2.png

The recovery plot:

_images/tutorial_che_3-results3.png

The yield plot:

_images/tutorial_che_3-results4.png

The total number of crystals plot:

_images/tutorial_che_3-results5.png

Files

Model report tutorial_che_3.xml
Runtime model report tutorial_che_3-rt.xml
Source code tutorial_che_3.py

Chem. Eng. Example 4

This example shows a comparison between the analytical results and the discretised population balance equations results solved using the cell centered finite volume method employing the high resolution upwind scheme (Van Leer k-interpolation with k = 1/3) and a range of flux limiters.

This tutorial can be run from the console only.

The problem is from the section 4.1.1 Size-independent growth I of the following article:

  • Nikolic D.D., Frawley P.J. Application of the Lagrangian Meshfree Approach to Modelling of Batch Crystallisation: Part II – An Efficient Solution of Integrated CFD and Population Balance Equations. Preprints 2016, 20161100128. doi:10.20944/preprints201611.0012.v1

and also from the section 3.1 Size-independent growth of the following article:

  • Qamar S., Elsner M.P., Angelov I.A., Warnecke G., Seidel-Morgenstern A. (2006) A comparative study of high resolution schemes for solving population balances in crystallization. Computers and Chemical Engineering 30(6-7):1119-1131. doi:10.1016/j.compchemeng.2006.02.012

The growth-only crystallisation process was considered with the constant growth rate of 1μm/s and the following initial number density function:

n(L,0): 1E10, if 10μm < L < 20μm
           0, otherwise

The crystal size in the range of [0, 100]μm was discretised into 100 elements. The analytical solution in this case is equal to the initial profile translated right in time by a distance Gt (the growth rate multiplied by the time elapsed in the process).

The flux limiters used in the model are:

  • HCUS
  • HQUICK
  • Koren
  • monotinized_central
  • minmod
  • Osher
  • ospre
  • smart
  • superbee
  • Sweby
  • UMIST
  • vanAlbada1
  • vanAlbada2
  • vanLeer
  • vanLeer_minmod

Comparison of L1- and L2-norms (ni_HR - ni_analytical):

--------------------------------------
         Scheme  L1         L2
--------------------------------------
       superbee  1.786e+10  7.016e+09
          Sweby  2.817e+10  8.614e+09
          Koren  3.015e+10  9.293e+09
          smart  2.961e+10  9.326e+09
             MC  3.258e+10  9.807e+09
           HCUS  3.638e+10  1.001e+10
         HQUICK  3.622e+10  1.005e+10
  vanLeerMinmod  3.581e+10  1.011e+10
        vanLeer  3.874e+10  1.059e+10
          ospre  4.139e+10  1.094e+10
          UMIST  4.363e+10  1.136e+10
          Osher  4.579e+10  1.156e+10
     vanAlbada1  4.574e+10  1.157e+10
         minmod  5.653e+10  1.325e+10
     vanAlbada2  5.456e+10  1.331e+10
 -------------------------------------

The comparison of number density functions between the analytical solution and the solution obtained using high-resolution scheme with the Superbee flux limiter at t=60s:

_images/tutorial_che_4-results.png

The comparison of number density functions between the analytical solution and the solution obtained using high-resolution scheme with the Koren flux limiter at t=60s:

_images/tutorial_che_4-results2.png

Files

Source code tutorial_che_4.py
Analytical solution fl_analytical.py
Flux limiters flux_limiters.py

Chem. Eng. Example 5

Similar to the chem. eng. example 4, this example shows a comparison between the analytical results and the discretised population balance equations results solved using the cell centered finite volume method employing the high resolution upwind scheme (Van Leer k-interpolation with k = 1/3) and a range of flux limiters.

This tutorial can be run from the console only.

The problem is from the section 4.1.2 Size-independent growth II of the following article:

  • Nikolic D.D., Frawley P.J. Application of the Lagrangian Meshfree Approach to Modelling of Batch Crystallisation: Part II – An Efficient Solution of Integrated CFD and Population Balance Equations. Preprints 2016, 20161100128. doi:10.20944/preprints201611.0012.v1

and also from the section 3.2 Size-independent growth of the following article:

  • Qamar S., Elsner M.P., Angelov I.A., Warnecke G., Seidel-Morgenstern A. (2006) A comparative study of high resolution schemes for solving population balances in crystallization. Computers and Chemical Engineering 30(6-7):1119-1131. doi:10.1016/j.compchemeng.2006.02.012

Again, the growth-only crystallisation process was considered with the constant growth rate of 0.1μm/s and with the different initial number density function:

n(L,0):                      0, if        L <= 2.0μm
                          1E10, if  2μm < L <= 10μm (region I)
                             0, if 10μm < L <= 18μm
      1E10*cos^2(pi*(L-26)/64), if 18μm < L <= 34μm (region II)
                             0, if 34μm < L <= 42μm
      1E10*sqrt(1-(L-50)^2/64), if 42μm < L <= 58μm (region III)
                             0, if 58μm < L <= 66μm
1E10*exp(-(L-70)^2/(2sigma^2)), if 66μm < L <= 74μm (region IV)
                             0, if 74μm < L

The crystal size in the range of [0, 100]μm was discretised into 200 elements. The analytical solution in this case is equal to the initial profile translated right in time by a distance Gt (the growth rate multiplied by the time elapsed in the process).

Comparison of L1- and L2-norms (ni_HR - ni_analytical):

-------------------------------------
        Scheme  L1         L2
-------------------------------------
      superbee  4.464e+10  1.015e+10
         smart  4.727e+10  1.120e+10
         Koren  4.861e+10  1.141e+10
         Sweby  5.435e+10  1.142e+10
            MC  5.129e+10  1.162e+10
        HQUICK  5.531e+10  1.194e+10
          HCUS  5.528e+10  1.194e+10
 vanLeerMinmod  5.600e+10  1.202e+10
       vanLeer  5.814e+10  1.225e+10
         ospre  6.131e+10  1.252e+10
         UMIST  6.181e+10  1.259e+10
         Osher  6.690e+10  1.275e+10
    vanAlbada1  6.600e+10  1.281e+10
        minmod  7.751e+10  1.360e+10
    vanAlbada2  7.901e+10  1.413e+10
-------------------------------------

The comparison of number density functions between the analytical solution and the solution obtained using high-resolution scheme with the Superbee flux limiter at t=100s:

_images/tutorial_che_5-results.png

The comparison of number density functions between the analytical solution and the solution obtained using high-resolution scheme with the Koren flux limiter at t=100s:

_images/tutorial_che_5-results2.png

Files

Source code tutorial_che_5.py
Analytical solution fl_analytical.py
Flux limiters flux_limiters.py

Chem. Eng. Example 6

Model of a lithium-ion battery based on porous electrode theory as developed by John Newman and coworkers. In particular, the equations here are based on a summary of the methodology by Karen E. Thomas, John Newman, and Robert M. Darling,

Thomas K., Newman J., Darling R. (2002). Mathematical Modeling of Lithium Batteries in Advances in Lithium-ion Batteries. Springer US. 345-392. doi:10.1007/0-306-47508-1_13

A few simplifications have been made rather than implementing the more complete model described there. For example, the following assumptions have (currently) been made:

  • two porous electrodes are used rather than providing the option for a “half cell” in which one electrode is lithium foil.
  • conductivity in the electron-conducting phase is infinite
  • constant exchange current density in Butler-Volmer reaction expression
  • no electrolyte convection
  • constant and uniform solvent concentration (ions vary according to concentrated solution theory)
  • monodisperse particles in electrode
  • no volume occupied by binder, filler, etc. in the electrode

The up to date version of the model is available at Raymond’s GitHub repository: https://github.com/raybsmith/daetools-example-battery.

The voltage plot:

_images/tutorial_che_6-results.png

The current plot:

_images/tutorial_che_6-results2.png

Files

Model report tutorial_che_6.xml
Runtime model report tutorial_che_6-rt.xml
Source code tutorial_che_6.py

Chem. Eng. Example 7

Steady-state Plug Flow Reactor with energy balance and first order reaction:

A -> B

The problem is example 9.4.3 from the section 9.4 Nonisothermal Plug Flow Reactor from the following book:

  • Davis M.E., Davis R.J. (2003) Fundamentals of Chemical Reaction Engineering. McGraw Hill, New York, US. ISBN 007245007X.

The dimensionless concentration plot:

_images/tutorial_che_7-results.png

The dimensionless temperature plot (adiabatic and nonisothermal cases):

_images/tutorial_che_7-results2.png

Files

Model report tutorial_che_7.xml
Runtime model report tutorial_che_7-rt.xml
Source code tutorial_che_7.py

Chem. Eng. Example 8

Model of a gas separation on a porous membrane with a metal support. The model employs the Generalised Maxwell-Stefan (GMS) equations to predict fluxes and selectivities. The membrane unit model represents a generic two-dimensonal model of a porous membrane and consists of four models:

  • Retentate compartment (isothermal axially dispersed plug flow)
  • Micro-porous membrane
  • Macro-porous support layer
  • Permeate compartment (the same transport phenomena as in the retentate compartment)

The retentate compartment, the porous membrane, the support layer and the permeate compartment are coupled via molar flux, temperature, pressure and gas composition at the interfaces. The model is described in the section 2.2 Membrane modelling of the following article:

  • Nikolic D.D., Kikkinides E.S. (2015) Modelling and optimization of PSA/Membrane separation processes. Adsorption 21(4):283-305. doi:10.1007/s10450-015-9670-z

and in the original Krishna article:

  • Krishna R. (1993) A unified approach to the modeling of intraparticle diffusion in adsorption processes. Gas Sep. Purif. 7(2):91-104. doi:10.1016/0950-4214(93)85006-H

This version is somewhat simplified for it only offers an extended Langmuir isotherm. The Ideal Adsorption Solution theory (IAS) and the Real Adsorption Solution theory (RAS) described in the articles are not implemented here.

The problem modelled is separation of hydrocarbons (CH4+C2H6) mixture on a zeolite (silicalite-1) membrane with a metal support from the section ‘Binary mixture permeation’ of the following article:

  • van de Graaf J.M., Kapteijn F., Moulijn J.A. (1999) Modeling Permeation of Binary Mixtures Through Zeolite Membranes. AIChE J. 45:497–511. doi:10.1002/aic.690450307

The CH4 and C2H6 fluxes, and CH4/C2H6 selectivity plots for two cases: GMS and GMS(Dij=∞), 1:1 mixture, and T = 303 K:

_images/tutorial_che_8-results.png

Files

Model report tutorial_che_8.xml
Runtime model report tutorial_che_8-rt.xml
Source code tutorial_che_8.py
Membrane unit membrane_unit.py
Variable types membrane_variable_types.py
Membrane model membrane.py
Support model support.py
In/out compartment compartment.py

Chem. Eng. Example 9

Chemical reaction network from the Dow Chemical Company described in the following article:

  • Caracotsios M., Stewart W.E. (1985) Sensitivity analysis of initial value problems with mixed odes and algebraic equations. Computers & Chemical Engineering 9(4):359-365. doi:10.1016/0098-1354(85)85014-6

The sensitivity analysis is enabled and the sensitivities are reported to the data reporter. The sensitivity data can be obtained in two ways:

  • Directly from the DAE solver in the user-defined Run function using the DAESolver.SensitivityMatrix property.
  • From the data reporter as any ordinary variable.

The concentrations plot (u1, u3, u4):

_images/tutorial_che_9-results1.png

The concentrations plot (u6, u8):

_images/tutorial_che_9-results2.png

The sensitivities plot (k2*du1/dk2, k2*du2/dk2, k2*du3/dk2, k2*du4/dk2, k2*du5/dk2):

_images/tutorial_che_9-results3.png

Files

Model report tutorial_che_9.xml
Runtime model report tutorial_che_9-rt.xml
Source code tutorial_che_9.py

The implementations of the COPS tests differ from the original ones in following:

  • The Direct Sequential Approach has been applied while the original tests use the Direct Simultaneous Approach
  • The analytical sensitivity Hessian matrix is not available. The limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm from IPOPT is used.

As a consequence, the results slightly differ from the published results. In addition, the solver options should be tuned to achieve faster convergence.

Chem. Eng. Optimisation Example 1

Optimisation of the CSTR model and Van de Vusse reactions given in tutorial_che_1:

Not fully implemented yet.

Reference: G.A. Ridlehoover, R.C. Seagrave. Optimization of Van de Vusse Reaction Kinetics Using Semibatch Reactor Operation, Ind. Eng. Chem. Fundamen. 1973;12(4):444-447. doi:10.1021/i160048a700

Files

Model report tutorial_che_opt_1.xml
Runtime model report tutorial_che_opt_1-rt.xml
Source code tutorial_che_opt_1.py

Chem. Eng. Optimisation Example 2

COPS test 5: Isomerization of α-pinene (parameter estimation of a dynamic system).

Very slow convergence.

Determine the reaction coefficients in the thermal isometrization of α-pinene (y1) to dipentene (y2) and allo-ocimen (y3) which in turn produces α- and β-pyronene (y4) and a dimer (y5).

Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS-273, 2004. PDF

Experimental data taken from: Rocha A.M.A.C., Martins M.C., Costa M.F.P., Fernandes, E.M.G.P. (2016) Direct sequential based firefly algorithm for the α-pinene isomerization problem. 16th International Conference on Computational Science and Its Applications, ICCSA 2016, Beijing, China. doi:10.1007/978-3-319-42085-1_30

Run options:

  • Simulation with optimal parameters: python tutorial_che_opt_2.py simulation
  • Parameter estimation console run: python tutorial_che_opt_2.py console
  • Parameter estimation GUI run: python tutorial_che_opt_2.py gui

Currently, the parameter estimation results are (solver options/scaling should be tuned):

Fobj  57.83097
p1    5.63514e-05
p2    2.89711e-05
p3    1.39979e-05
p4   18.67874e-05
p5    2.23770e-05

The concentration plots (for optimal ‘p’ from the literature):

_images/tutorial_che_opt_2-results.png

Files

Model report tutorial_che_opt_2.xml
Runtime model report tutorial_che_opt_2-rt.xml
Source code tutorial_che_opt_2.py

Chem. Eng. Optimisation Example 3

COPS test 6: Marine Population Dynamics. (Not working properly)

Given estimates of the abundance of the population of a marine species at each stage (for example, nauplius, juvenile, adult) as a function of time, determine stage specific growth and mortality rates.

Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS-273, 2004. PDF

Experimental data generated following the procedure described in the COPS test.

Run options:

  • Simulation with optimal parameters: python tutorial_che_opt_3.py simulation
  • Parameter estimation console run: python tutorial_che_opt_3.py console
  • Parameter estimation GUI run: python tutorial_che_opt_3.py gui

Currently, the parameter estimation results are (suboptimal results obtained, solver options/scaling should be tuned):

Fobj =  1.920139e+8
m(0)    3.358765e-01
m(1)    4.711709e-01
m(2)    1.120524e-01
m(3)    8.509170e-02
m(4)    9.683579e-02
m(5)    1.919142e-01
m(6)    2.418778e-01
m(7)    2.421000e-01
g(0)    1.152995e+00
g(1)    7.529383e-01
g(2)    5.024174e-01
g(3)    5.704327e-01
g(4)    4.180333e-01
g(5)    3.185407e-01
g(6)    2.250250e-01

The distribution moments 1,2,5,6 plots (for optimal results from the literature):

_images/tutorial_che_opt_3-results.png

The distribution moments 3,4,7,8 plots (for optimal results from the literature):

_images/tutorial_che_opt_3-results2.png

Files

Model report tutorial_che_opt_3.xml
Runtime model report tutorial_che_opt_3-rt.xml
Source code tutorial_che_opt_3.py

Chem. Eng. Optimisation Example 4

COPS test 12: Catalytic Cracking of Gas Oil.

Determine the reaction coefficients for the catalytic cracking of gas oil into gas and other byproducts.

Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS-273, 2004. PDF

Experimental data generated following the procedure described in the COPS test.

Run options:

  • Simulation with optimal parameters: python tutorial_che_opt_4.py simulation
  • Parameter estimation console run: python tutorial_che_opt_4.py console
  • Parameter estimation GUI run: python tutorial_che_opt_4.py gui

Currently, the parameter estimation results are (solver options/scaling should be tuned):

Fobj = 4.841995e-3
p1   = 10.95289
p2   =  7.70601
p3   =  2.89625

The concentration plots (for optimal ‘p’ from the literature):

_images/tutorial_che_opt_4-results.png

Files

Model report tutorial_che_opt_4.xml
Runtime model report tutorial_che_opt_4-rt.xml
Source code tutorial_che_opt_4.py

Chem. Eng. Optimisation Example 5

COPS test 13: Methanol to Hydrocarbons.

Determine the reaction coefficients for the conversion of methanol into various hydrocarbons.

Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS-273, 2004. PDF

Experimental data generated following the procedure described in the COPS test.

Run options:

  • Simulation with optimal parameters: python tutorial_che_opt_5.py simulation
  • Parameter estimation console run: python tutorial_che_opt_5.py console
  • Parameter estimation GUI run: python tutorial_che_opt_5.py gui

Currently, the parameter estimation results are (solver options/scaling should be tuned):

Fobj = 1.274997e-2
p1 = 2.641769
p2 = 1.466245
p3 = 1.884254
p4 = 1.023885
p5 = 0.471067

The concentration plots (for optimal ‘p’ from the literature):

_images/tutorial_che_opt_5-results.png

The concentration plots (for optimal ‘p’ from this optimisation):

_images/tutorial_che_opt_5-results2.png

Files

Model report tutorial_che_opt_5.xml
Runtime model report tutorial_che_opt_5-rt.xml
Source code tutorial_che_opt_5.py

Chem. Eng. Optimisation Example 6

COPS optimisation test 14: Catalyst Mixing.

Determine the optimal mixing policy of two catalysts along the length of a tubular plug flow reactor involving several reactions.

Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS-273, 2004. PDF

In DAE Tools numerical solution of dynamic optimisation problems is obtained using the Direct Sequential Approach where, given a set of values for the decision variables, the system of ODEs are accurately integrated over the entire time interval using specific numerical integration formulae so that the objective functional can be evaluated. Therefore, the differential equations are satisfied at each iteration of the optimisation procedure.

In the COPS test, the problem is solved using the Direct Simultaneous Approach where the equations that result from a discretisation of an ODE model using orthogonal collocation on finite elements (OCFE), are incorporated directly into the optimisation problem, and the combined problem is then solved using a large-scale optimisation strategy.

The results: fobj = -4.79479E-2 (for Nh = 100) and -4.78676E-02 (for Nh = 200).

The control variables plot (for Nh = 100):

_images/tutorial_che_opt_6-results.png

The control variables plot (for Nh = 200):

_images/tutorial_che_opt_6-results2.png

Files

Model report tutorial_che_opt_6.xml
Runtime model report tutorial_che_opt_6-rt.xml
Source code tutorial_che_opt_6.py

Optimisation tutorial 1

This tutorial introduces IPOPT NLP solver, its setup and options.

Files

Model report opt_tutorial1.xml
Runtime model report opt_tutorial1-rt.xml
Source code opt_tutorial1.py

Optimisation tutorial 2

This tutorial introduces Bonmin MINLP solver, its setup and options.

Files

Model report opt_tutorial2.xml
Runtime model report opt_tutorial2-rt.xml
Source code opt_tutorial2.py

Optimisation tutorial 3

This tutorial introduces NLOPT NLP solver, its setup and options.

Files

Model report opt_tutorial3.xml
Runtime model report opt_tutorial3-rt.xml
Source code opt_tutorial3.py

Optimisation tutorial 4

This tutorial shows the interoperability between DAE Tools and 3rd party optimization software (scipy.optimize) used to minimize the Rosenbrock function.

DAE Tools simulation is used to calculate the objective function and its gradients, while scipy.optimize.fmin function (Nelder-Mead Simplex algorithm) to find the minimum of the Rosenbrock function.

Files

Model report opt_tutorial4.xml
Runtime model report opt_tutorial4-rt.xml
Source code opt_tutorial4.py

Optimisation tutorial 5

This tutorial shows the interoperability between DAE Tools and 3rd party optimization software (scipy.optimize) used to fit the simple function with experimental data.

DAE Tools simulation object is used to calculate the objective function and its gradients, while scipy.optimize.leastsq function (a wrapper around MINPACK’s lmdif and lmder) implementing Levenberg-Marquardt algorithm is used to estimate the parameters.

Files

Model report opt_tutorial5.xml
Runtime model report opt_tutorial5-rt.xml
Source code opt_tutorial5.py

Optimisation tutorial 6

daeMinpackLeastSq module test.

Files

Model report opt_tutorial6.xml
Runtime model report opt_tutorial6-rt.xml
Source code opt_tutorial6.py

Optimisation tutorial 7

This tutorial introduces monitoring optimization progress.

Files

Model report opt_tutorial7.xml
Runtime model report opt_tutorial7-rt.xml
Source code opt_tutorial7.py