RSSInterested in linking to "Simulation, the game we all can play!"?
You may use the Headline, Deck, Byline and URL of this article on your Web site. To link to this article, select and copy the HTML code below and paste it on your own Web site.
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.
12/20/2004
To cast your vote, log in or become a member. This quick, one-time registration gives you access to members-only site benefits.
1 vote
Page 2 of 2
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 heider@wuche.che.wustl.edu.
Page 2 of 2
To cast your vote, log in or become a member. This quick, one-time registration gives you access to members-only site benefits.
1 vote
