Simulation, the game we all can play!

Dec. 20, 2004
Washington University's (St. Louis, MO) Robert Heider provides a game plan for using simulation in process control.
 By Robert Heider, PE
There are two classes of simulators, ones that require the user to program in FORTRAN, C++, Visual Basic etc. and ones that advertise that these skills are not required. Most of these programs never work the way they should with your process, but usually work well with the sample problems, which are usually very simple or elementary.

Why simulate at all? From a control systems engineer's viewpoint, it is helpful to be able to predict the dynamics thereby gaining the knowledge of how loops may behave as well as allowing the simulation to be used as a training tool. The chemical engineers want to have some data to support the plant design and a simulator will help justify equipment design and sizes.

I have been asked to sit in on many discussions and vendor presentations about simulators. The subject comes up about what type of physical property database to use or how their product interfaces to the users database or how to enter physical and chemical properties of compounds not in their database. Clearly, this is the single most important factor from the chemical engineering simulation viewpoint. They usually divert their discussion to the DIPPR database while my mind keeps wondering about the dynamics that they say they have, usually as an option. Most of these products require a large cash outlet as well as ongoing customer support from some hot line or email service.

Control Deficiencies
When you use one of these simulators for control, there are some watchouts; does you simulator actually simulate the following:

Can the dynamics simulate true deadtime? This is important because if it weren't for deadtime, most process control engineers would be out of a job. The control blocks should emulate real time as much as possible and therefore need deadtime. If your simulator does not do deadtimes, do not give up, it must do first order filtering. Just add multiple first orders in series, which looks like a first order with deadtime. This may not be sufficient, as I will explain later.

Noise simulation: This is necessary where analytical control is being implemented, noise needs to be filtered which adds lag to the process. The simulator should be able to add Gaussian noise with an adjustable standard deviation. This is very important for analytical data. Remember if the process model is only 85% accurate, it is good enough for process control. Frequently other engineers are very possessive and not willing to let others see the model because if is not accurate enough. For control purposes, accuracy is not required, just some overall dynamic behavior. The difference of 10% one way or another usually doesn't effect control simulations. Remember the objective of the simulator is to show the approximate behavior, and test various control strategies.

Simulated Physical Properties
This is the major watch out one should have with commercial packages. In order for the simulator to produce good results, you need a good estimate of the physical properties. How do they get these properties in their database? They perform regressions on experimental data or they use some equations based on molecular structure. If the data is experimental, what data did they use? Over what range was the data regressed? Frequently the data entered is only based on a few points and well outside the range where you may want to run your simulation. A good example of this is formaldehyde and water solutions.

An Actual Example
In order to minimize the size of a surge tank, a commercial package, with dynamic properties, was used to simulate a process chilled cooling system using a 50-50 mixture of ethylene glycol and water as the cooling medium. Forty percent of the system capacity is a batch process so there was some concern that a large surge tank is in place in order to allow the chiller to maintain the utility at the desired temperature.

Specifications:
Outlet Temperature: -4 degF
Cooling Media: 50-50  EG-H2O solution
Flow Rate: 600 GPM
Cooling capacity: 146 Tons
Piping: 1000 feet 6" SCH40

A commercial simulator was used to test the effect of various size surge tanks on dynamic response of the chiller's outlet temperature controller.

The first problem was the apparent inability of the simulator to calculate the viscosity of the glycol and water mixture. The simulator gave a viscosity several times the viscosity given by the vendor of the ethylene glycol and water solution. The density given was also in error. I question how these packages simulate liquid mixtures.

Once the simulation was configured, it showed that there was no control problem and that the return temperature changed instantaneously and as a result no surge capacity was required! This, of course, cannot be. The simulation failed to show the effect of the installed deadtime due to the piping. This turned out to be a very important criterion, because the liquid in the piping acted as surge capacity.

Piping should be considered both as a delay line and as surge. Consider the following:

The pipeline, 1000 feet in length (500-foot supply and return branches) of 6" SCH40 pipe
contains 1.5 gal/foot. When flowing at 600 GPM or 6.6 feet/second, the total delay per branch would be 500/6.6 or about 75 seconds. With a sample time of 5 seconds this would be 75/5 = 15 pipe segments with 500/15 = 33 feet/segment which is 33.33 feet/segment * 1.5 gal/foot = 50 gal/ segment.

In this case, it became necessary to simulate the pipeline as a series of 50 gallon tanks, the idea that deadtime can be simulated as a series of first order lags. Once this change was made, to the simulation, the simulator showed the first order effect, but still failed to show the effect of true deadtime. The system was finally simulated in EXCELâ to control the integration scan time, which was necessary to see the effect of the piping deadtime. Once the simulation was complete, it demonstrated that there is sufficient reserve of coolant in the supply and return piping to eliminate the surge tank. The following plots show the correct results of the simulation. Notice the dead time between the load temperatures.

This simple simulation demonstrates the potential problems with commercial process simulator packages.

Integration Algorithms
So how can you integrate with a spreadsheet? I’ll bet it is too complicated.

Writing differential equations is not a very complicated exercise. The integration algorithms are found in most college math texts or on the Internet. Writing a differential equation is just a matter of writing a difference equation and solving it with an algorithm. Remember, difference equations first, then given to solve through the algorithm.

Integrators come in two varieties, relative to step size, variable or fixed. Examples of variable sized are Gear’s Stiff and Runge-Kutta-Fehlberg. These algorithms reduce the integration interval to minimize the error. There is the Simpson’s rule that states that the error in an algorithm is related to the size of the integration interval raised to the power of the number of evaluations. Variable size algorithms keep reducing the interval until the error is less than that specified by the user. These algorithms are frequently used by those who are making a space shot or in determining the reaction rate constants for some process chemistry. They are of little use for control simulations because in control, fixed intervals are required. Even if the process is static, control must continue and mixing different algorithms is seldom worth the effort. The Runge-Kutta level 4, meaning four levels or intermediate solutions per iteration, is the best for most control purposes. Even a simple Euler algorithm can give good results.

To write a differential equation, just remember you only need to write the difference, which is the input minus the output. Just remember solving a differential equation required the initial value, so you have to start somewhere. The following is a Visual Basic code for a differential equation program solving a simple pressure control loop using the ideal gas law and the Universal Gas Sizing Equations. In this case, a tank is blanketed with nitrogen. The nitrogen is supplied through a regulator upstream of an orifice plate. The tank is vented through a control valve. x is a three element array, x(1) = mass in the tank, x(2) = valve time constant, x(3) = integrated error term. x is set to the initial condition values. The differential equations solve for a new x term based on the differences. In this example, it is best to determine the mass in the tank by assuming an initial pressure.

Runge-Kutta level 4 integration example
‘ first set the existing variables to a xnew
For j = 1 To ssv
   xnew(j) = x(j)
Next
' begin the four step process
For j = 1 To 4
   ' calculate the initial mass in tank  p=nRT/V
    P_tank = (xnew(1) / MW) * R_gas * t / vol
   ' calculate the change in pressure
   ' first we need the outlet flow, n_out
   dP = P_tank - p_atm
   If (dP > P1_orifice - p_atm) Then
   dP = P1_orifice - p_atm
   End If
   If (dP < 0#) Then
   dP = 0#
   End If
    Qout = Cg_valve(Index) * P_tank * _
        ((520 / (spgr * t)) ^ 0.5) * _
        Sin((59.64 / C1_valve) * ((dP / P_tank) ^ 0.5))
    Qout = Qout / minperhr
    n_out = Qout * spgr / 13.1
   ' then calculate the flow across the orifice plate, n_in
   dP = P1_orifice - P_tank
   If (dP > P1_orifice - p_atm) Then
    dP = P1_orifice - p_atm
    End If
   If (dP < 0#) Then
    dP = 0#
    End If
   Qin = Cg_orifice * P1_orifice * _
      ((520 / (spgr * t)) ^ 0.5) * _
      Sin((59.64 / C1_orifice) * ((dP / P1_orifice) ^ 0.5))
   Qin = Qin / minperhr
    n_in = Qin * spgr / 13.1
 
‘ for x(1), the mass difference is the inlet minus the outlet
‘ for x(2), the first order lag simulating valve travel
‘ for x(3), the integration of the control error

   x_dot(1) = n_in - n_out
   x_dot(2) = (1 / tau_v) * (u(Index) - xnew(2))
   x_dot(3) = (Psp(Index) - P_tank)

‘ These are the integration equations     
If j = 1 Then
    For i = 1 To ssv
    k1(i) = step * x_dot(i)
    xnew(i) = x(i) + k1(i) / 2
    Next
End If
If j = 2 Then
    For i = 1 To ssv
    k2(i) = step * x_dot(i)
    xnew(i) = x(i) + k2(i) / 2
    Next
End If
If j = 3 Then
    For i = 1 To ssv
    k3(i) = step * x_dot(i)
    xnew(i) = x(i) + k3(i)
    Next
End If
If j = 4 Then
    For i = 1 To ssv
    k4(i) = step * x_dot(i)
    x(i) = x(i) + (1 / 6) * (k1(i) + 2 * k2(i) + 2 * k3(i) + k4(i)) ' Calculated next state
    Next
End If
Next

' calculate the new outlet pressure by calculating the mass in tank  p=nRT/V
P_tank = (x(1) / MW) * R_gas * t / vol

Using the Control System as a Simulator
In many cases where the user is only interested in simulating hydraulic or thermal systems, or where chemical reactions are simplified or ignored, the control system itself can be modified to provide the simulated process. The following example illustrates this, the controller is diagramed in figure 2 and the simulation is shown in figure 3. Assume that using an internal cooling coil cools a heat-jacketed vessel. The temperature is controlled by the amount of cooling water through the coil and the heat is controlled to keep the flow it a high enough value it facilitate good heat transfer.

Modifying the control program can simply be done adding the required function blocks to calculate the heat transfer. The cooling heat value can be subtracted from the heating heat value and be totalized. The resultant total functions as an integrator, which is the heat value in BTU/minute, This can be converted to a temperature reading.

In the simulation blocks the following functions are simply calculated:
The cooling flow, F_COIL, is equal to K*(TV-104) where TV-104 in percent.
The heat transfer coefficient, U_COIL, is equal to a K0 + K1*(TV-104).
The coil outlet temperature is calculated by a heat balance across the coil is equal to the heat flowing through the coil:
Q_COIL1 = U_COIL * A *( (TI-104R) – (TOUTCOIL + Tincoil) /2 ) =
 F_COIL * (TOUTCOIL – Tincoil)
Solving the above for TOUTCOIL:
TOUTCOIL =
(U_COIL * A ((TI-104R) – (Tincoil/2) + F_COIL*Tincoil) / (F_COIL + U_COIL*A/2)      
The heat transferred through the coil, Q_COIL1, is equal to F_COIL*(TCOILOUT – Tincoil) / 60.
Tincoil is a constant of 529.7 DegR.
The electrical power percentage, P-104, is converted to heat in BTUs. This value is as lagged by a first order lag equal to the thermal time constant across the jacket, tau.
Q_ELECT = K*(P-104)* (1.0 - e-t/tau  )
The Q_COIL is a lagged value of the calculated Q_COIL1.

Once the electrical heat and cooling heat values are calculated, the subtracted element, D, subtracts the cooling heat from the electrical heat. This heat value is then totalized. The controlled variable, TI-104, is finally calculated in degrees C. The temperature in degrees R, TI-104R, is calculated also.

Robert L. Heider, PE, is Adjunct Professor in the Chemical Engineering Department of Washington University, St. Louis, MO. He can be reached by phone at 314-935-6070; by fax at 314- 935-7211; or by e-mail at [email protected].