Russell inside out for Dist column 3

[dinu24]

Hi all,

I was wondering if anyone can direct me to a good reference on how to implement the Russells inside out algorithm for distillation column model calculations. What I need is to learn the theory and how it works. Also possibly get into the numeric method...

Alternatively, is there any other algorithm that you'd recommend for steady state distillation column model?

Thanks!
Dinu

 

[DickRussell]

The original article is:

Russell,R.A., "A Flexible and Reliable Method Solves Single-Tower and Crude-Distillation Column Problems," Chem. Eng., October 17, 1983, p.53.

 

[dinu24]

Thanks DickRussell. I appreciate your feedback.

 

[dinu24]

Can someone suggest a good example of a steady state distillation simulation, with all the necessary specs, to test out my implementation?

I tried simulating the deethanizer in the "Fractionation Train" but I dont have all the necessary operating specs. I am trying to get Ethane and Methane in the distillate, but the column is not converging for what I have specified. Here are the details of what I am trying;

The distillation column has 14 trays, with the inlet stream at (1.5, 24.6, 170.3, 31.0, 76.7, 76.5 moles/hr) of C1, C2, C3, iC4, nC4, C5 respectively. The column pressure is at about 2725kPa, boilup at 90`C and reflux at -4`C. I am trying a reflux ratio of 2.5. But I cant get the column to converge.

Any suggestions or thoughts? Or is there any other example you 'd suggest I try to test my implementation?

Thank you!

 

[UmeshMathur]

It is a highly non-trivial task to implement Russell's complete distillation column algorithm. Failure to converge can arise from: (a) poor convergence of thermodynamics routines, (b) poor initial guesses for column internal flows, (c) inadequate attention to the computational and convergence algorithms for the inside-out column model itself, and a number of other issues too intricate to discuss here.

Instead, you'd be well advised to consider buying Russell's excellent PD Plus, a full-blown keyword-driven process simulator that is available for a very modest fee. This inside-out implementation for columns is the easily the fastest in the industry, usually beating the major competitors by one or even two orders of magnitude. See deerhaventech.com for details.

This package includes all the major thermo options, unit operations, multiple recycle convergence, software controllers, in-line fortran, and a host of other convenience features. I have used it very happily for over ten years in petroleum refining, petrochemicals, complex chemicals, etc.

 

[dinu24]

Thank you UmeshMathur.

I agree with you on the Russell's algorithm implementation being a non trivial task. I design and implement unit operations and thermodynamic VLE calculations (based on both EOS and Activity Coef methods). And in fact, the implementation I am testing is actually a port from another language, and I have tested it with the separation of Ethanol and Water. What I want is to test it with multiple components. I want a case that should converge with PR EOS, so I can see how this column behaves. I have not implemented any short cut methods to design columns yet, so I am looking for help...

I'll check out PD Plus. It seems like a great tool from what you say. Thanks for the tip

 

[DickRussell]

In order to comment on convergence of your deethanizer test column, we need to know a few more things.

1. What are the temperature and pressure of the feed?

2. Where does the feed enter the column?

3. Is the distillate taken as a vapor product?

4. Besides reflux ratio, what other specification is there?

5. What initial estimates of product rate, internal rates, and temperatures have been given?

The usual implementation of the method is such that when a column model fails to converge the last profiles available usually indicate what the problem may be. A good implementation will not require really good initial estimates of flows and temperatures for most columns, although some systems do need something sensible, information usually available by inspection of the physical problem being solved. More often, convergence failure is due to asking for something physically infeasible.

At first glance, that reflux ratio spec could be a problem if the feed is all or partially vapor as it enters the column. The overhead vapor rate would have to be roughly (2.5 +1)*(1.5 +24.6) = 91 mols, which is only 24% of the feed. If a lot of the feed is vapor, it will run up the column to the condenser, and become either product or reflux. Assuming you want a decent split between the C2 and C3, you thus know fairly well what the product rate will be. The excess must be reflux, and the ratio of that to the product rate won't be your 2.5 spec value.

If your feed does not provide a lot of vapor upon entry, so that the reboiler must generate boilup to satisfy a reflux spec, then the problem could be something else.

Just eyeballing the expected top and bottom temperatures at the expected top and bottom compositions, they don't seem especially out of line with that pressure. You're getting up there toward the critical region, though. If you are running into physical limitations of the column and the system it is handling, you could always do a run at a high reflux ratio and lower column pressure to get a good run, use the results to initialize subsequent runs, and gradually move toward the reflux and higher pressure you think you want. At some point, you should be able to see limitations of the system become obvious in the resulting profiles.

Provide more details, and we can provide better advice.

 

[DickRussell]

To add to what I wrote earlier today and see if there is anything else wrong with the system, I made a quick run of your test column in PD-PLUS. I had to make some assumptions, of course. I assumed the feed to be bubblepoint liquid, halfway up the column. I also had to assume the column to be isobaric, since only one pressure was given; that is not a real issue, and a realistic profile could be entered easily enough. I fixed the vapor distillate at the sum of the C1 and C2 in the feed, plus the reflux ratio (2.5); naturally there would be some C2 lost to the bottoms and some C3 taken overhead, according to how sloppy the split is. I also used your top and bottom temperature estimates.

This first run converged readily, showing that the split was indeed sloppy; a large fraction of the ethane went to the bottoms. This warranted a step back to basics. For a reasonable split between the C2 and C3 (no more than a mole of either lost to the other product for starters), and using the relative volatility from the stage reports of the first run, a quick Fenske equation calc showed minimum stages to be 15-18, and minimum reflux ratio somewhere above 10. This led to a run with 30 stages and reflux ratio of 20. The results overshot the assumed split, with only half a mol of C2 going down and the same of C3 overhead. But now the column configuration is in the ballpark. Further refinements can be added, such as separation specs in place of product rate and reflux.

Below are the input data for the column block for the last run, plus the convergence history and final profiles. I hope line wrap doesn't foul up the appearance when this is posted.

*COLUMN DEETH 'Deethanizer column', STAGES= 30,
FEED FD STAGE 15,
REBOILER, CONDENSER,
VAPOR DISTILLATE FLOW= 26.1 KGMOL/HR, STREAM C1C2,
BOTTOMS STREAM C3PLUS,
REFLUX RATIO= 20.0,
PRESS= 2725 KPA, STAGE 30,
EST TEMP= 0 C, STAGE 30,
EST TEMP= 50 C, STAGE 25
EST TEMP= 90 C, STAGE 1,
PRINT STAGE 1 28 30,
DAMPING FACTOR 0.7

--CALCULATION HISTORY--
COLUMN BLOCK DEETH
AVG. HEAT/SPEC. ERROR= 0.02307
Calculating new matrix
Inverting the matrix
AVG. HEAT/SPEC. ERROR= 0.00123 DISTANCE FACTOR= 1.0000
ITER. 1 ERRORS: HEAT/SPECS=0.00123, VOLATILITY=0.06934, ENTHALPY=0.00000
AVG. HEAT/SPEC. ERROR= 0.00156
AVG. HEAT/SPEC. ERROR= 0.00144 DISTANCE FACTOR= -0.0625
AVG. HEAT/SPEC. ERROR= 0.00112 DISTANCE FACTOR= 0.2500
AVG. HEAT/SPEC. ERROR= 0.00088 DISTANCE FACTOR= 0.2500
ITER. 2 ERRORS: HEAT/SPECS=0.00088, VOLATILITY=0.02230, ENTHALPY=0.00000
AVG. HEAT/SPEC. ERROR= 0.00146
AVG. HEAT/SPEC. ERROR= 0.00119 DISTANCE FACTOR= 0.2500
AVG. HEAT/SPEC. ERROR= 0.00090 DISTANCE FACTOR= 1.0000
AVG. HEAT/SPEC. ERROR= 0.00008 DISTANCE FACTOR= 1.0000
AVG. HEAT/SPEC. ERROR= 0.00001 DISTANCE FACTOR= 1.0000
ITER. 3 ERRORS: HEAT/SPECS=0.00001, VOLATILITY=0.01293, ENTHALPY=0.00000
AVG. HEAT/SPEC. ERROR= 0.00082
AVG. HEAT/SPEC. ERROR= 0.00005 DISTANCE FACTOR= 1.0000
AVG. HEAT/SPEC. ERROR= 0.00001 DISTANCE FACTOR= 1.0000
ITER. 4 ERRORS: HEAT/SPECS=0.00001, VOLATILITY=0.00351, ENTHALPY=0.00000
AVG. HEAT/SPEC. ERROR= 0.00023
AVG. HEAT/SPEC. ERROR= 0.00001 DISTANCE FACTOR= 1.0000
ITER. 5 ERRORS: HEAT/SPECS=0.00001, VOLATILITY=0.00093, ENTHALPY=0.00000
AVG. HEAT/SPEC. ERROR= 0.00007
AVG. HEAT/SPEC. ERROR= 0.00000 DISTANCE FACTOR= 1.0000
ITER. 6 ERRORS: HEAT/SPECS=0.00000, VOLATILITY=0.00023, ENTHALPY=0.00000
AVG. HEAT/SPEC. ERROR= 0.00002
ITER. 7 ERRORS: HEAT/SPECS=0.00002, VOLATILITY=0.00007, ENTHALPY=0.00000
SOLUTION TIME= 0 MIN., 0.03 SEC.
--OUTPUT REPORT--



COLUMN BLOCK DEETH Deethanizer column

COLUMN SUMMARY. FLOWS ARE IN KGMOL/HR

STAGE TEMP PRESS FLOW FROM STAGE FEED PRODUCT HEAT ADDED
C KPA VAPOR LIQUID MM KCAL/HR
30 4.4 2725.0 536.1 26.1 VAP -1.606
29 10.1 2725.0 562.2 522.1
28 17.4 2725.0 548.2 497.2
27 28.0 2725.0 523.3 472.3
26 39.9 2725.0 498.4 453.2
25 50.0 2725.0 479.3 445.3
24 57.0 2725.0 471.4 444.6
23 61.2 2725.0 470.7 445.7
22 63.8 2725.0 471.8 446.6
21 65.4 2725.0 472.7 446.7
20 66.5 2725.0 472.8 445.8
19 67.6 2725.0 471.9 443.3
18 69.2 2725.0 469.4 437.9
17 71.8 2725.0 464.0 426.9
16 76.9 2725.0 453.0 406.4
15 86.4 2725.0 432.5 788.2 380.6
14 87.5 2725.0 433.7 792.3
13 88.3 2725.0 437.8 795.4
12 89.0 2725.0 440.9 797.9
11 89.6 2725.0 443.4 800.1
10 90.1 2725.0 445.6 802.0
9 90.6 2725.0 447.5 803.6
8 91.0 2725.0 449.1 804.9
7 91.3 2725.0 450.4 805.9
6 91.7 2725.0 451.4 806.6
5 92.2 2725.0 452.1 806.7
4 93.0 2725.0 452.2 805.5
3 94.5 2725.0 451.0 801.2
2 97.9 2725.0 446.7 787.6
1 105.8 2725.0 433.1 354.5 LIQ 1.779


FEED STREAM FD TO STAGE 15 = 380.6 KGMOL/HR
VAPOR STREAM C1C2 FROM STAGE 30 = 26.1 KGMOL/HR
BOTTOMS STREAM C3PLUS = 354.5 KGMOL/HR
EXTERNAL REFLUX RATIO= 20.541
INTERNAL REFLUX RATIO= 20.003

 

[dinu24]

Thank you very much, DickRussell. Your input is greatly appreciated.

What you have shown here, gives me a start, so I'll play around this a bit and see where I get to. Obviously, something like this, that seems to take you a few minutes (a couple of hours maybe) will take me days! So, I'll get to this again on thursday and see what I can make out of it.

Just the fact that you have taken the trouble to get back to me with such detail, suggests PDPlus must be a really cool product. I am ofcourse, assuming that you are a part of its development, in which case such support alone is worth much much more than the cost of the product. I have downloaded the demo copy.

Thanks again for your input, I'll keep you informed on where I get to with this!

 

[UmeshMathur]

dinu24:

Did you note that the solution time for Dick's example problem was 0.03 CPU seconds? Most other simulators will take a lot more time. Even if you change specifications to product purities, the execution time for PD Plus goes up only a few percentage points.

I have run complex many refinery fractionators with PD Plus; most of them converge in under 1 CPU second on a decent 2 GHz or faster machine. So, if you know what you're doing, this is a great way to go.

 

[dinu24]

UmeshMathur;
Yes, 0.03 CPU seconds is very impressive. Seems almost impossible without very close start values, doesn't it? I am wondering if these results were from calculations based on very good start values, either from previous iteration or user provided. Either ways, to be honest with you, the results are much better than what I have got so far and convergence is much faster too. And you are not surprised, are you?

The column I am working with is a basic tower, without condenser and reboilers attached to it. So I have to use heaters/ coolers and flash tanks to recycle reflux and boilup. This causes a time step loss, making convergence harder in my case. The results seem ok in my case, suggesting the algorithm is implemented correctly but I need to improve the numerical methods, which I have had to tweak several times, to get convergence.

Anyways, thanks for your response. I am impressed with PDPlus so far. Just wondering how I can benefit from it though. I am a chemical engineer with a graduate degree in software systems, and currently, I specialize in the development of process systems. Having said that, I do interact with people from the industry, often on a commercial level, so it is good to know of great products like this one.

 

[DickRussell]

For the first run, I used for temperature estimates only the top and bottom estimatets you provided. Aside from the top product rate specified, no internal flow estimates were given; the program provided that from the reflux ratio spec. I usually use results of an initial run to provide better guesses for subsequent runs. Most of the time the method is quite good at finding the solution from minimal or even poor estimates. Since the engineer knows what the physical system is supposed to do, providing some decent estimates of product rate and temperatures for the first run is fairly easy to do. I advocate making use of prior calculations to help the software along, but people often don't do this, particularly if they notice that usually the software comes up with the answer anyway.

I am curious as to why your column model has no reboiler and condenser. If you have gone to the trouble of programming the I/O method, adding reboiler and condenser is easy. For a fixed duty, the duty is just another quantity in the heat balance term for that stage. In the more normal case, with duty to be calculated, the duty is just whatever heat is needed (in or out) to heat balance the stage, and some other spec takes the place of that heat balance term, such as a purity or recovery. Turning a whole column model into column plus separate flash drums and heat exchangers introduces recycle streams. This of course takes longer to converge, and it makes the use of product quality specs much more difficult, likely requiring more loops outside of the recycle loops. Having the column model solved all at once is far more efficient.

 

[dinu24]

DickRussell;
Honestly, the basic tower design just came about without a good design effort in the start. Since I was working on modeling the dynamic behaviour of the column, the basic structure of the module was taken from that. I guess, it was assumed that the design would be simpler if we limited the complexity of the column. And keeping the condenser and reboiler external was not a very smart thing to do, now that I am working on your algorithm. I will soon implement the "distillation column" and absorber and I am really excited about it!!

DickRussell/UmeshMathur;

As an aside, have you taken a look at possibilities for modeling and simulations in the "micro" process engineering? I am really intrigued by the new developments in these areas, and about the possibilities for the future. Since both of you have a lot of experience in process engineering, I thought I'll ask you

 

[dinu24]

Hi DickRussell;

I am reading your article "A flexible and reliable method solves single-tower and crude-distillation-column problems" Chem. Eng Oct17, 1983, and have a couple of questions for you;

In your article, you say;

1. "Often, distillation calculations are performed with recycle loops of a flowsheet simulation. In such cases, the jacobian inverse need be obtained for only the first time through the tower" - Do you mean after the first time the column is solved, we need not obtain the jacobian inverse, but can simply use the one obtained in the previously converged calculation?

2. "This noncorrespondence between variables and errors means that the initial Jacobian must be obtained by first perturbing the variables one at a time, then inverting. Thereafter, the Broyden formula updates the inverse after each iteration." - Are you suggesting I use the "Sherman-Morrison" formula to calculate the inverse or do you think I should approximate the Jacobian with the Broyden method and then invert it with LU or QR?

I hope you dont mind my asking these Q's. I look forward to your reply...

Thanks!
Dinu
dinu.biz

http://www.dinu.biz

 

[DickRussell]

The Jacobian inverse can be stored and reused to start off a re-calculation of the column. If the new feeds to the column are not much different and the specifications have not been changed much, then the saved matrix often is a good estimate of what it would be if a new Jacobian were to be calculated and inverted. Of course you need to check the progress toward convergence. If it seems that perhaps the old matrix inverse is not leading to steady progress toward convergence then it would be appropriate to get a new Jacobian and invert it. When the method was first developed, the computers were nowhere near as fast as today's, and the effort to calculate a new Jacobian by differencing for a column with many stages was substantial. The time to invert it (I just used a canned matrix inverter that was in the public domain) went up with the cube of the dimension of the matrix. Updating the matrix inverse from the changes made to the variables and the reductions in the set of errors was far faster. Today's computers are so fast that you can afford to be less fussy about doing every possible thing to obtain calculation efficiency and still come up with a reasonably fast program.

2. I don't know what the "Sherman-Morrison" formula is. The full Jacobian is obtained by perturbing the variables one at a time and collecting the sets of error terms corresponding to each perturbation. Of course this is not an analytical derivative, but a good approximation based on differencing (delta error/delta variable). You need to get the size of delta big enough so that the changes in errors are not down in the noise region, but not too big so as to reduce the accuracy of the derivatives. Any tool for inverting the matrix can be used, perhaps one from some math library you have. The Broyden inverse update formula is used after each successful use of the inverse to reduce the set of errors. It is not exact, of course, but it does give a very good approximation of the inverse corresponding to the new set of variables, allowing the next prediction of changes to the variables.

 

[dinu24]

Thanks DickRussell! ...Time to get my hands dirty again!!!

http://www.dinu.biz