Apart from large scale dynamic simulation models for Operator Training Simulators, we also regularly do smaller projects and interesting case studies. In this blog post I will discuss the use of a dynamic simulation model during the design phase of a steam boiler project of one of our customers. They were interested in how their new boilers would behave dynamically and how their existing control structure would be affected by these new boilers. With the aid of our dynamic simulation tool, they managed to successfully make some necessary changes to both their equipment and their control system, before the boilers were in service. Process Description and Project Goal (Note: All equipment tag names and process values have been edited on purpose, for confidentiality reasons) A Chemical Company invested in two steam boilers, which should provide steam to a 64 barg steam header and a connected steam network (including 20, 15, 12 barg headers and several lower pressure headers). The network is also kept on pressure by external suppliers. Before the company invested in the new boilers, the steam network was pressurized by the steam of a large 75 barg network of an external supplier (via a pipeline of several kilometers). The new setup would decrease the capacity (smaller holdup because of shorter distances) and therefore the company was interested in the dynamic behavior of the various steam headers in case of failure of one of the steam suppliers or other non-nominal behavior (e.g. trips in parts of the connected plants). Based on findings, the planned rules could then be adjusted. Process Modelling Using the supplied documents (PFD’s, P&ID’s, equipment details, etc.) and consultation with company process engineers, Mobatec set up a dynamic model, tested and tuned for the steam network. After the model was delivered, the company could perform a lot of different scenarios with the model. The model consisted of the following components:

• Simplified models of all boilers, based on the boiler volume and the heat input to generate steam.
• The steam network with the main steam headers and the important control valves and controls.
• The main customers of steam headers were modeled as “tunable customers”, such that they show realistic steam consumption under different conditions.
• All controllers, control loops and logics of interest were included in the model.

For the simulation only the water / steam system was important. For the process conditions on the part of the steam consumers, reasonable estimates were made, as far as was relevant, in consultation with company’s process engineers. Results The company engineers thought of a lot of different possible scenarios and it was easy for them to implement and run these and plot and discuss relevant results. As expected, the dynamic model confirmed that the steam headers reacted a lot quicker to offsets and disturbances with the new boiler setup as compared to previous setup (with the external supplier of 75 barg steam). The new boiler setup was designed to operate at roughly 50% of its capacity, such that the boilers could ramp-up (or down) in case of a disturbance (like the failure of an external boiler). It was found that for most disturbances the new design and the proposed control system were just ok and the pressures in the headers would not drop below critical values, such that no (parts of) plants should be shut down. An important change had to be implemented in the field, though, to make sure that certain plant parts would not trip after a disturbance in the steam network. The model predicted that the current size of a main “let down valve” between the 20 barg and 12 barg header was too small to cope with dynamic changes, causing one header to blow off steam, while the other dropped in pressure. Therefore, the actual valve was replaced before the new boilers where in place. All in all, it was a very nice project, with very positive results for the involved companies. Our customer gained better insight in the dynamic behavior of their steam network and we used the feedback of them to further improve our tools. Let us know what you think about using dynamic simulation for industrial cases in the comment box below, such that we can all learn something from it. To your success! Mathieu. ———————————————–

The modelling methodology behind Mobatec Modeller assumes that any process can be broken down into systems and connections. Systems represent a capacity, able to store mass and energy. Connections represent the transfer of mass and energy between the defined systems.
So far, however, we have dealt with “lumped” systems as a basic modelling element. Lumped systems are assumed to have uniform properties over their entire volume. A typical example of such a system is a simple (one phase) CSTR, as the whole volume is assumed to be ideally mixed. In that case, mathematically speaking we are calculating the assumed “uniform point” properties & behaviour and are using the volume of the system to quantify the system, i.e. calculate the defined system performance. From a mathematics point of view this is DAE (Differential Algebraic Equation) problem, containing Ordinary Differential Equations (ODE) and Algebraic Equations (AE), where the time is only independent variable. ODE’s are representing the balances (mass & heat) equations of the systems while AE’s are called constitutive equations and are used to define each system properties.
What is a modelling strategy towards non-uniform systems?
A cascade of lumped systems can be used to simulate (already existing) non-uniform systems successfully enough, but if the model goal is the design, optimization and model based innovation of such systems, then the model made would need to take spatial coordinates of the modelled system into account as independent variables as well. From mathematics point of view this problem then becomes PDAE – Partial Differential and Algebraic Equations problem. In this months blog we are going to talk about the Distributed Parameter System modelling in Mobatec Modeller, a solution strategy, and numerical methods used.
Distributed parameter systems are systems that cannot be considered to have uniform properties over their entire volume, hence they are called “distributed parameter” as their properties are time and space dependent.
Typical representatives of such a systems are plug flow & packed bed reactors, fixed or moving bed columns, catalytic pellets, etc.
Dominant phenomena in these systems are convective and dispersive/diffusive mass & heat transfer. In order to solve these type of problems space discretization numerical methods are required.
Space discretization numerical methods in Mobatec Modeller
Numerical discretization of a time-dependent partial differential equation (PDE) includes the following approaches: Fully-discretized in both time and space and Method of Lines. In Mobatec Modeller we employ the Method of Lines approach (MOL). (MOL) approach: discretize in space first, transform PDEs to ordinary differential equations (ODEs); and then using the time integration method of ODEs to advance in time.
There are several spatial discretization methods: Finite Difference Method (FDM), Finite Volume Method (FVM), Finite Element Method (FEM), Spectral and Pseudo-spectral Method. Finite Difference Methods (FDM) are adopted as dominant approach to numerical solution of PDE’s and are used in Mobatec Modeller.
Finite Difference Methods make derivatives approximation of univariate functions by replacing the derivatives in equations with finite difference approximations. This results in a number of ordinary differential & algebraic equations that can be solved simultaneously. The basic idea of finite difference schemes is to replace derivatives by finite differences. The most common method of obtaining the finite-difference expressions is via Taylor series expansion: The Taylor series expansions for the nearby points of Fx are used to obtain the First and later the Second derivative approximations at Fx. Depending whether we employ expression for the point before F(x-h) or the point after F(x+h) the point Fx, or combine both expressions together we will get: Forward, Backward or Central Finite Difference approximation of the derivatives. (FFDM, BFDM & CFDM)
Forward and Backward finite difference methods can be used alone, depending on the boundary conditions provided for the system, i.e. do we know the inlet or the outlet boundary condition. If we want to employ a central finite difference method then all three methods are needed to form a numerical solution strategy. This will be explained in detail later on, but let us first “get to know” each method separately: FFDM is using the point after the focus point. The arrow tip is showing the integration direction, which is from Z=Z towards Z=0. So we use FFDM if the outlet BC is known. BFDM is using the point before the focus point. The arrow tip is showing the integration direction, which is from Z=0 towards Z=Z. So we use BFDM if the inlet BC is known. “O” is representing the terms in the Taylor series expansion that are being neglected and is representing the error of the finite approximation. First order error is of the order of spatial step taken O(h). Second order error derivative approximation is obtained by including one more point in the approximation expression. For FFDM that is F(x+2h) point, and for BFDM is F(x-2h) point. The error is of the order of spatial step squared O(h2) which is more accurate solution taking that h->0. CFDM is using the point before and the point after the focus point. The arrow tip is showing the integration direction, which is from Z=0 towards Z=Z. In order to employ CFDM we use BFDM at the inlet BC point and FFDM at the outlet BC point. In the middle – throughout the distributed domain CFDM is used. This subject will be tackled in detail in a next blog once we get to know the methods properties.
Below, the graphical interpretation of the finite difference derivatives approximation by used method is given: The smaller the spatial step (h->0) the better the finite approximations of the derivatives are. The same goes with higher error order. In real case the desired solution would be an optimal number of spatial elements (spatial step size) over the distributed domain length. Too much spatial steps -> very small step size h, would take too much computational time and effort.
Let us take a look at the second order derivatives finite approximations. If we take the derived first derivative finite approximations and reuse them in the Taylor series expansion, now deriving the second derivative finite approximations we will get: The error order is of the spatial step size so they are not used that much. Main use is for purely convective terms -> first derivative finite approximation.
CFDM can be successfully used for mixed convective and (purely) diffusive/dispersive terms. Knowledge on numerical methods is essential for their successful application in real life engineering problems. In the next blog on this topic, we will explain the application of FD methods in Mobatec Modeller and why it is important that this part of the PDE solution strategy needs to remain a “white box” – available for user to interact with in order to be able to make the models to converge and not have their hands “tied together”.
What is your experience with Distributed Parameter System modelling and numerical methods used? Did you ever run into convergence problems without knowing what the cause was?
I invite to post your experiences, insights and/or suggestions in the comment box below, such that we can all learn something from it.
To your success! Vladimir.

While I was reading back through an earlier blog post on automatic component distribution, I noticed the following sentence: “To enhance the definition of the species topology, permeability and directionality are introduced as properties of mass connections. They constrain the mass exchange between systems by making the species transfer respectively selective or uni-directional.”
I realised it may not be entirely clear what is meant with the terms “permeability” and “directionality”, so I decided to spend some more lines on that:
The modelling methodology implemented in Mobatec Modeller makes use of several topologies, the Physical Topology, Species Topology and Equation Topology.
The idea of the construction of the species topology is to facilitate the definition and identification of species and reactions occurring in various parts of the processing plant. The species topology defines which species are present in which part of the process and thus defines the basics of the chemistry/biology of a processing plant model. This enables the generation of the relevant balances of extensive quantities for the systems.
The species topology is superimposed on the physical topology. Its construction efficiently aids the model designer in the definition of the chemistry/biology of the processing plant because it generates knowledge about which species are present in every system in a physical topology and about which reactions occur in every reactive system.
It should be noted that at this stage of the model definition process the quantities of the species in the systems are of no relevance. The association of a certain species with a certain system only implies that this system has a capacity for this particular species. If, for example, a system Σ is associated with only one species, S, then this implies that species S is the only species that can be present in system Σ. It does not imply that S has to be present in Σ. At certain moments or maybe all the time during an evolution/simulation of the process, species S may even be absent. At these particular moments the system has no mass, but of course still has a capacity to store the species S.
The existence of certain species and reactions in one part of the system does not imply that they are the same in all parts of the system. In principle, however, every chemical or biological species could be present in any part of the system, so one could decide not to construct the species topology and only generate a global set of species. This would mean that one defines only one set of species, which comprises all species that appear in the whole plant. Every system of the physical topology would then consist of all these species. For each species in each system the species mass balances could then be generated, which would of course result in a flood of equations. In most cases, many of these equations would be obsolete, because species are usually constrained to a group of systems and therefore do not appear in every system. The species of the fluids in the two parts of a shell and tube heat exchanger, for example, are usually not the same.
To aid in the definition of the species topology, species and reaction databases are used and three concepts are introduced, “permeability”, “directionality” and “reactions”, which support the idea of locality of chemical or biological species. These three concepts will facilitate the definition of the species topology and will be illustrated in the next sections. The permeability and directionality concepts are associated with connections and therefore with the common boundary element which the two connected systems have in common.
Permeability of Mass Connections
The permeability is defined as a property of a mass connection, which constrains the mass exchange between systems by making the species transfer selective. The mass connection is then defined as a semi-permeable thermodynamic wall. An example of a process where the permeability concept could come in handy is a simple diffusion process. The application of the selective constraint on the transportation of species can also assist in the representation of separation processes, such as ultra-filtration or reverse osmosis.
———————————————————————————-
Example: Diffusion of a Component from one Phase to Another Consider a tank which contains two immiscible phases A and B (figure a). Phase A consists of species P and Q. Phase B consists of species P and R. The simplest way to describe the diffusion of species P from phase A to phase B (or the other way around) is shown in figure b, which represents a physical topology of the diffusion process.
The two phases are each modelled as a simple lumped system and between the systems exists a mass connection which describes the transfer of all species between those systems. Note that the direction of the connection is only a convention, defining a reference direction for the actual stream. The actual direction of the flow depends on the driving force, which in this case would be the difference in chemical potentials of P in each phase.
To prevent the species Q and R from being transferred through the connection, the permeability concept is introduced: The connection between the two systems is defined as a semi-permeable thermodynamic wall, which is only permeable for species P. This means that only the species P will be transferred through the connection and that, consequently, the transfer terms of the other species (Q and R) will not appear in the (component) mass balances of the systems because they are not transferred.
———————————————————————————-
Note that in principle the permeability concept does not have to be limited to mass c.q. species. The concept can be readily extended to any extensive quantity. In practice, though, this is not necessary as one normally distinguishes between mass, heat, momentum and work transfer, which are all handled individually.
Directionality of Mass Connections
The concept of (uni- or bi-) directionality is also a property that is related with mass connections. The species topology, as it was described so far, allows for bi- directional flow in any mass connection. Since the flow is (normally) dependent on the state of the two connected systems, nothing can be said a priori about the direction of the actual flow.
In fact, in many situations the flow would change in the course of the time interval over which the process is being observed. Species that are not inhibited by a defined permeability may thus appear on both sides of the connection, when computing the species topology. In many applications, though, one-directional flows are assumed and thus species propagate in a constraint flow pattern across the network. This desired behaviour calls for defining two classes of mass connections, namely uni-directional and bi-directional, where the latter should be the normal case as it leaves the system unconstraint. Inversion of the flow for a uni-directional flow would have to be handled as an exception as it contradicts the validity range of the model.
Reactions
The reaction concept more or less speaks for itself. A chemical or biological reaction is a mechanism that transforms a set of species, the reactants, into a set of different species, the products. At this stage of defining the species topology it is only necessary to know which species are transformed and which species are generated in every (reactive) thermodynamic system. The properties of reactions (kinetic laws, heat of reaction, etc.) and the conditions under which the reaction may take place are introduced at a later stage.
I invite to post your experiences, insights and/or suggestions in the comment box below, such that we can all learn something from it.
To your success!
Mathieu.
———————————————–

The dynamic behaviour of a process can be abstracted to the point at which it can be represented by the smooth, continuous change of a series of state variables. This allows such a process to be represented mathematically as a set of DAEs.
The majority of processes, however, cannot be considered entirely continuous, but also experience significant discrete changes or discontinuities superimposed on their predominantly continuous behaviour. Such discrete changes, often referred to as events, can either be planned changes, called time events (for example, planned operational changes, such as start-up and shut-down, feed stock and/or product changes, process maintenance, switching a valve) or non-planned changes, called state events (overflow of a tank, (dis)appearance of a phase, inversion of a flow).
The exact time of occurrence of time events is either known a priori or it can be computed from the occurrence of a previous event. The time of occurrence of state events is not known in advance because it depends on the system fulfilling certain state conditions.
In addition to the two types of events we can classify events in another way, namely in events that do not change the (dynamic) model structure and events that change the (dynamic) model structure. With events that do not change the model structure, the balance equations are not altered, but some algebraic equations are changed. This usually results in a non-smooth continuation of the primary state variables.
———————————————————————————-
Example: Tank with Fast and Equilibrium Reaction Consider a tank with an inlet (consisting of two components A and B) and an outlet. The contents of the tank consist of four components (A, B, C and D) among which two fast reactions occur, namely:
A + B → C  (1)
A ↔ D        (2)
The first reaction (1) is a very fast, one way reaction. This means that there is no measurable inverse reaction and one of the two components A or B will always be completely consumed (i.e.: n_A =0 or n_B =0). The other reaction (2) is an equilibrium reaction, which means that there will always be a certain ratio between components A and D (c_A = k ∗ c_D). Because both reactions are fast in the context of the viewed time scale, the reaction terms that appear in the balance equations can be removed be using the Simple Index Reduction Algorithm (as discussed in previous blogs). The figure shows the result of a simulation of the tank. In this figure we can see three events taking place.
At the start of the simulation there is an excess of A in the tank, which means that there will be no B, since all B that is available directly reacts with A to form C. The inlet flow contains more B than A and this results in diminishing quantities of A (and, consequently, D) in the tank until the amount of A (and D) in the tank reaches zero. At this point we encounter the first event, for now there will be an excess of B available in the tank and A (and D) will be non-existent.
At a certain point in time the input is changed such that the inlet flow will contain more A than B. A second event is the result. The amount of B will decrease until it reaches zero. This will initiate the third event, such that A (and D) will be present again.
Obviously, the first and last event were state events, while the second event was a time event.
———————————————————————————-
In our Modelling Methodology the model structure is typically never changed when an event occurs. Special IF-THEN-ELSE equations are used to describe state-dependent events. As an example, consider the overflow equation that was mentioned in our blog post on the Overflowing Tank:
F2:  IF V > Vmax THEN
         F2 = K * (V – Vmax)
       ELSE
         F2 = 0
       END
Instead of changing the balance equations when there is no overflow (if the volume in the tank is lower than the maximum volume Vmax), the value of the overflow is simply set to zero.
For those who are interested in the simultaneous solving of (large) mathematical sets of differential and algebraic equations (DAE) that also include conditional statements (IF-THEN-ELSE type equations), it should be noted that the conditions are checked before a calculation step is performed. So, with the information of the previous time step it is determined which equations are used for solving the DAE. The effect of a state-dependent conditional statement is, therefore, typically “too late” if the calculation step is too large. This effect can give rise to “interesting” modelling results and should be well understood when making dynamic models. I will most likely discuss it in another blog post ;).
What is your experience with events and discontinuities in your models? Did you ever run into unexpected model behaviour because of these events?
I invite to post your experiences, insights and/or suggestions in the comment box below, such that we can all learn something from it.
To your success!
Mathieu.
———————————————–

In our last blog post we discussed about fundamental time scale assumptions and the implication these assumptions can have on the final model. This month we discuss types of assumptions that have less effect on the structure and contents of the model, but that can certainly lead to (computational) problems if they are not recognized.
Constraints and assumptions describe all kinds of algebraic relationships between process quantities which have to hold at any time. A volume constraint, for example, restricts the volume of a simple system or the sum of volumes of a set of systems. Assumptions may result in models which are cumbersome to solve. Potential problems can very often be avoided at an early stage of the model development by keeping clear of certain assumptions or by directly dealing with the cause of the potential problems.
As discussed in the previous blog post one can distinguish between several types of assumptions: Structural assumptions (i.e. the construction of the physical topology), order of magnitude assumptions (very small versus very large) and assumptions on relative time scale (very slow versus very fast). These assumptions are usually introduced with a goal, namely to simplify the description of the behaviour by neglecting what is considered insignificant in the view of the application one has in mind for the model. Whilst indeed such assumptions do simplify the description, they are also the source of numerous problems, such as “index problems”, which make the solving of the equations very difficult.
 

High-index Models

Dynamic process models, as derived with our modelling methodology, consist of differential and algebraic equations (DAEs). Unfortunately, most engineers have little knowledge of the theory of DAEs, since most of the calculations that have to be performed during education are steady-state simulations. In the rare case that dynamics are considered, a mathematical description of the (usually very simple) model is derived in the form of ordinary differential equations (ODEs).
One of the major advantages of writing a model in DAE form as opposed to ODE form, is that a modeller does not have to perform a set of often cumbersome mathematical manipulations, such as substitution and symbolic differentiation. Dynamic process models can be divided into “low index” models (index 0 and 1) and “high index” models (index two and higher and some index one models). The differential index or index is a measure of the problems related to initialisation and integration of dynamic process models. The problems related to solve a dynamic process model increase with increasing index.
We will not go into details nor give a mathematically correct definition of the index of a model at this stage, but will just state that the index of a DAE is the number of times that all or part of a DAE must be differentiated (with respect to time) in order to get an ODE.
It is not recommend, by the way, to actually perform this series of differentiations as a general solution procedure for DAEs. Rather the number of such differentiation steps that would be required in theory turns out to be an important quantity in understanding the behaviour of numerical methods.
According to the definition an ODE (either explicit or implicit) has index zero. DAEs with index zero and one are generally much simpler to understand (and much simpler to solve) than DAEs with index two or higher.
With our modelling method we strive to produce semi-explicit index one models, since these can be easily used for simulation by any DAE-solver. Higher index models must either be simulated directly by a special integrator which tackles high index DAEs, or be transformed to semi-explicit index one and integrated. Two simple tests that guarantee structurally semi-explicit index one models are:

• All algebraic variables must be present in the algebraic equations.
• It must be possible to assign each algebraic equation to an algebraic variable. Assignment of all equations and variables must be possible.

If one or several algebraic variables are absent from the algebraic equations, then the model index is two or higher.
 

Assumptions Leading to High-index Models

Many high-index problems are caused by a model purpose that is not carefully considered, because the modeller does not want to include some of the rapid dynamics in the model, by assumptions that may not be essential or because the modeller wants to include certain variables in the model.
Our modelling method forces a model designer to be more aware of the assumptions he makes. Therefore, potential high index model formulations can be detected and/or avoided at an early stage of the model development. Knowing the causes of high index formulations helps the model designer in carefully considering the modelling purpose and the assumptions he wants to make.
Assumptions that impose direct or indirect constraints on the differential variables lead to high index models. But a modeller cannot simply impose some constraints on the differential variables. A constraint is always imposed by some “driving force” (i.e. a flow or reaction, since these are the only forces that appear in the differential equations), which forces the differential variables to adhere to the constraint. This means that instead of giving a description for the rate, a flow or reaction remains “unmodelled” and a (direct or indirect) constraint on the differential variables is given.
There can be several reasons why a modeller wants to introduce assumptions:

• Only slow dynamics of the process are of interest. In this case the rapid dynamics can be neglected.
• Difficulties in finding reliable rate equations may force a modeller to make quasi steady-state assumptions.
• In order to perform model reduction, simplifying assumptions may be introduced.

In our previous blog post the term time scale was introduced. Modelling systems in a range of time scales, the capacity terms are chosen accordingly but also the transport and the production terms. For parts being outside of the time scale in which the dynamics are being modelled, a pseudo steady-state assumption is made. For example, (very) fast reactions – fast in the measure of the considered range of time scale – are assumed to reach the equilibrium (for all practical purposes) instantaneously, and very slow ones do not appreciably occur and may be simply ignored.
In a next blog post we will zoom into “steady-state assumptions” and “events and discontinuities”.
What is your experience with constraints and assumptions? Did you ever run into computational problems without knowing what the cause was?
I invite to post your experiences, insights and/or suggestions in the comment box below, such that we can all learn something from it.
To your success!
Mathieu.
———————————————–