When I first started to develop the fixed bed reactor model, I quickly learned that corrections were easier to make if the difference formulas were defined as functions rather than "hard coded" directly into the program. One reason for this is because each formula is used twice in the program due to the hopscotch method. Also, the use of functions makes the program more compact and easier to read. This post shows an example of one of the difference functions.
The Differential Equation
The energy balance for the fluid will be used in the example. The model includes balances for both the fluid and the catalyst, so the heat of reaction does not appear in the fluid energy balance.
All of the variables and parameters are dimensionless. The thetas are temperatures of fluid (f) and particle (p). I will define the alpha parameters in a later post. The radial variable is x and the axial variable is y.
The difference formula
The difference form of the above equation, when solved for the fluid temperature at location x,y is show below:
The above formula is for an explicit step in the axial direction from a previously computed set of values at row y-1.
The value computed by this formula would be on the left side of the arrow. This result will be defined as the function Exyft:
The argument list
In order from left to right: weight fraction vectors for fluid and particle, temperatures of fluid and particle, pressure, x and y positions for the point being calculated, component number, radial and axial step sizes. You may wonder why n (component number) is in the argument list, but it is not used in the expression on the right of the arrow. It is needed because all of the functions to be returned from the subroutine, fb_diff, must have the same number and type of arguments. This is a Mathcad requirement.
The alpha parameters are not included in the argument list. The fb_diff routine obtains their values by calling the fb_parms routine prior to defining the difference functions.
The function naming convention
The example shown is for a general point. By that, I mean it is for a point not on the centerline or at the wall. Different formulas, and thus different functions, are needed for each type of point. In addition, another set of functions is needed for the implicit steps.
To easily identify the functions, I used the following convention for their names. Each character position in the name has a meaning (except the third). Starting with the first character and moving to the right, the meanings are as follows:
As another example of a function name, IRypn signifies the implicit value of the particle mass fraction at the wall. The component number desired is supplied in the argument list. The "n" in the function name only means the function is for a mass fraction. In hindsight, I probably should have used "m" for mass in the naming convention instead of "n".
Using the functions
The picture below shows how a function is used in the integration routine.
With the difference functions defined, the main programming task for the integration routine, fb_ss, is to keep track of the alternating explicit, implicit steps as the calculation proceeds down the reactor. In the picture above, the explicit steps are at the even radial positions. In the next row of points down the reactor, the explicit steps will occur at the odd radial positions.
Next: the alpha parameters