Simulation, the game we all can play!

Washington University's (St. Louis, MO) Robert Heider provides a game plan for using simulation in process control.

Share Print Related RSS
Page 2 of 2 1 | 2 Next » View on one page

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)
' 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
End If
If j = 2 Then
    For i = 1 To ssv
    k2(i) = step * x_dot(i)
    xnew(i) = x(i) + k2(i) / 2
End If
If j = 3 Then
    For i = 1 To ssv
    k3(i) = step * x_dot(i)
    xnew(i) = x(i) + k3(i)
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
End If

' 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:
(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
Page 2 of 2 1 | 2 Next » View on one page
Share Print Reprints Permissions

What are your comments?

You cannot post comments until you have logged in. Login Here.


No one has commented on this page yet.

RSS feed for comments on this page | RSS feed for all comments