INDUSTRIAL AND ENGINEERING CHEMISTRY
November 1951
(16) where g is the acceleration due to gravity. The punched card table of thermodynamic properties may be utilized in the calculation of specific impulse in the following atepa : 1. I n the table of h(po, T),perform an inverse interpolation to determine the adiabatic flame temperature, To= T(p., HI). 2. I n the table of s(pa, T),perform a direct interpolation to determine the adiabatic, S s(pe, !Pa). 3. In the table of B(p0, T),perform an inverse interpolation to determine the exhaust temperature, TO= T ( p ,S): 4. In the table of h(p0 T),perform a &re& lnte olation to h(p0, determine the specific enthalpy, H., of the exhaust,%
-
-
To).
At the conclusion of these s t e p , the application of Equation 16 to determine the specific impulse is immediate. This method of evaluation is independent of the assumption of the validity of the polytropic expansion law.
2475
triggered by the finder card. The computer automatically collates the finder card with the table deck to determine the point of entry in the table, forms the proper differences, and carries out the required direct or inverse interpolation by means of the fourth difference Stirling formula. The-choice between direct and inverse interpolation ia effected by means of the position of a toggle - switch on the control unit. The time required to carry out the interpolations for flame temperature and specific impulse is completely determined by the time required for the reader to read the table, and it is negligible in comparison to the time required for the construction of the table of thermodynamic properties.
ACKNOWLEDGMENT This research is part of the work being done a t the U. S. Bureau of Mines on Project NA onr 25-47, supported by the Office of Naval Research and the Army Air Forces. The calculations on the electronic numerical integrator and calculator of the equilibrium composition and thermodynamic properties of combustion gases were carried out under arrangements made with the Office of the Chief of Ordnance, Department of the Army.
LITERATURE CITED
-
(1) (2)
Beattie, J. A., C h m . Reus., 44, 141-92 (1949). Brinkley, 9.R., Jr., J . C h m . Phya., 14, 663-4, 686 (1946).
(4)
Brinkley, S. R., Jr., and Lewis, B., Chem. Eng. News. 27, 2540-1
(a) IW., 15, 107-10
Anmar urd8
card
Figure 4.
ORDERS Interpolate dim! Interpolate inverse
(1947).
(1949).
Brinkley, 8. R., Jr., and Lewis, B., U. S. Bur. Mines, Rept. Inueut. 4806, in press. (6)Brinkley, S. R., Jr., and Smith, R. W., Jr., “Proc. Scientific Computation Forum 1948,” pp. 77-82, New York, International Business Machines Corp., 1950. (7) Brinkley, 8. R., Jr., and Smith, R. W., Jr., “Proc. Seminar on Scientific Computation 1949,” pp. 58-63, New York. International Business Machines Corp., 1950. (8) Brinkley, 8. R., Jr., Wagner, G. L., and Smith, R. W., Jr., “Proc. Technical Computation Forum 1950,” New York, International Business Machines Corp., 1951. (9) Hottel, H. C., and Eberhardt, J. E., C h m . Revs., 21, 439-60
(5)
CONTROLS internal, triggered by finder card
Block Diagram of Interpolation Calculator
More complicated cycles involve the same kind of operations in the punched card table of the thermodynamic properties of the combustion gases, supplemented in most cases by similar considerations involving a supplemental punched card air table. The interpolations are efficiently carried out by the card programed calculator, utilizing special control panel circuits that make the equipment an interpolation calculator. This method of operation is illustrated by the block diagram of Figure 4. In this case, the table deck, preceded by a “finder” card listing the argument for direct interpolation of the function for inverse interpolation, is read by the card reader of the program control unit. The table deck contains no control punohes of any kind. All control is internal, accomplished by control panel wiring and
(19371,
Hottel, H. C., Williams, G ,C., and Satterfield, C. N., “Thermodynamic Charta for Combustion Procemes,” New York, John Wiley & Sons, 1949. (11) Kandiner, H. J., snd Brinkley, S. R., Jr., IND. ENG.CHEM.,42, (10)
850-5 (1950). (12)
Scarborough, J. B., “Numerical Mathematical Analysis,” pp. 178, 187, Baltimore, Johns Kopkins Press, 1930.
RBOlIVlD
Apfl 12, 1961.
Monte-Carlo Method for Solving Diffusion Problems GILBERT W. KING A r t h r D. Little, he., and Massachnsetts Institate of Technology, Cambridge, Mass.
T
HE Monte-Carlo method is a new means of solving problems in physics and engineering made available by large scale computing machines. An automatic computer should not be considered as a glorified slide rule that will do engineering calculatione as they have been done in the past, but as an organism that will attack the problems in a new and powerful way. Clae&a1 mathematics is only a tool for engineera and physicists and is
not inherent in the realities with which they attempt to -,a]. It has been customary to idealize and simplify the mechanisms of the physical world in the form of differential and other types of equations of classical mathematics, because solutions or metho& of attack have been discovered during the last few hundred years with means generally available-namely, pencil, paper, and logarithm tables.
2416
INDUSTRIAL A N D ENGINEERING CHEMISTRY
Vol. 43, No. 11
Engineering problems tothe computin machines, and asolution has%een reached by day can be formulated in this T h e Monte-Carlo method is a new means of solving an entirely different computmanner but are so complex problems in physics and engineering made available by ing scheme from an that that analytical solutions are large scale computing machines. would be used by hand? The not known. With the advent main feature in the use of The problems in physics have difficulties which are automatic computing of a u t o m a t i c computing much deeper than those of engineering, whereas the machines is the simple nature machines the tendency has engineering problems are much more complex. However, of the basic routine, which is been for the differential equalarge acale computing machines have considerable capacity repeated a large number of tion to be put into the form times. and can handle a great deal of complexity. At present of difference equations and they are designed to do arithmetic only and are not parapproximate solutions found A thousand trials a r e ticularly adapted to the hightr mathematics of physics. by numerical methods. Howusually adequate. In some A very realistic model of situations involving ciaasical ever, there is no fundamental cases, however, only the tails physics or chemistry can be inserted into a computer. reason to pass through the of the final distribution are of The actual calculations used in the MontcCarlo method abstraction of the differential interest, and these consist of are extremely elementary, consisting of many additions equation. Any model of an only a fraction of the lo00 of unity in certain counters chosen by random numbers. engineering or physical procinitial particles. As a rule, a ess involves certain assump convenient procedure of biastions and idealizations which inrr the random walks can be are more or less openly implied in setting up the mathematical found so that most of the lo00 arrive in the tail, with weight factors which are easy to calculate. equation. By making other simplifications, sometimes less stringent, the situation to be studied can be put directly in the comGENERALIZATION OF THE PHYSICAL PROCESS puting machines, and a more realistic model is obtained than is permissible in the medium of differential or integral equations. In diffusion the particles do not move the same distmce, &, each time, but they move a distance following a certain distribuDIFFUSION IN ONE DIMENSION tion of free paths. Modern computing machines are constructed in such a way that it is easy to modify the calculations by the These ideas are illustrated by consideration of a simple diffusion machine to make the particle move not a constant distance each problem. Assume that we wish to find the distribution of a dye time but a distance that is chosen from a random sampling of the which is being supplied to the center of the capillary tube and distribution of the mean-free paths, which applies in this case, diffuses in time in both directions along the capillary tube. It is very often Gaussian. Thus, u more realistic picture is put directly custoniary to describe this phenomenon as a differential equation on the machines. It is true that in this simple case the more for the dye concentration, U. realistic picture can also be introduced into the diffusion equation and explicit solutions obtained with a little more difficulty. - = Dbau bX
=
GENERALIZATION OF THREE DIMENSIONS
for which, in this simple example, explicit solutions in closed form are known. The process can be described directly to the machine in the following way: The dye molecules are subject to Brownian motion and move a mean distance, Ax, in a time, At. The direction in the capillary is purely random so that Ax may be either positive or negative. Let us take an arbitrary particle and throw a coin. If i t is heads, z is increased by Ax, and if it is tails, x is decreased by Ax. This is done by assQning a counter to x. The quantity, AX (usually unity in an appropriate decimal position)l is entered into this counter every machine cycle. The counter 1s impulsed to add or subtract by a device equivalent to a relay. The latter is thrown one way or the other by a random process. A convenient source is to supply random digits. These are the numbers 0 to 9 supplied in a sequence whlch has been teated for various types of randomness. The machine examines these for odd- or evenness, and the relay is picked up or not, according to the result of this test. For example, if the random numbers supplied to the machines were 5, 7, 7, 0, 5, 1, 3, 0, 9, 4, then the particle moves backwards or forwards according to the signs -, -, -, f, -, -, -, Its successive ositions are -1, -2, -3, -2, -3, -4, -, -5, -4, -5, -4, a n t i t s final position after 10 steps (at time 10At) is x = -4. After throwin the coin for a hundred times, the particle has moved a given &stance, say xt, in the time t = 1OOAhl. If now the coin is thrown a hundred times a t random for another particle, it will follow an entirely different path (usually referred to a~a “random walk”), unrelated to the first one, and will arrive at a different point, x2. If this is done a thousand times for a thousand different particles and the number of particles are counted a t x = Ax, 2Ax, 3Ax, -&, etc., we have a 100At. This is the distribudistribution of articles a t time t tion in which t i e particles would actually have fallen had they been acted on by those particular selections of Brownian motions. The mathematical solution of the diffusion equation is an approximation of the distribution. The mathematical splution of the diffusion equation applies to the ideal situation of infinitely small steps. I n setting up the diffusion equation, more assumptions were used than were put directly, in an elementary fashion, into
+,
+.
-
Another advantage of the MonteCarlo method is that the model can be set up for three-dimensional problemn just as easily as for one. Three coordinates Can be counted simultaneously in the counters of the machine, and modern equipment usually has a large counter capacity. On the other hand, three-dimensional problems in analysis are not a simple extension of one-dimensional results but go beyond the power of present-day mathematics. There are no solutions in the form of elementary functions unless there is some symmetry (which really reduces the problems to fewer dimensions). The author has shown elsewhere ( 9 )that an efficient scheme for carrying out diffusion in three dimensions is to divide the space by a tetrahedral lattice, so that individual steps are along one of the vectors, (1, 1, l),(1, -1, -l), ( -1, 1,-l), (-1, -1, l ) , or the negatives of these. Also (8, I),it is wasteful for a step in the random walk to move the particle to a point just visited. This means that the particle need choose only one of three vectors radiating from the point a t which it is. This can be accomplished as follows: Assume the particle is at the lattice point, x, y, z, having moved there by vector vi. A single random number (modulo 3) is supplied. If i t is 0, a new vector is created by changing the sign of the x component of vi; if it is 1, g is changed; and if it is 2, z is changed. For example, let a sequence of random numbers be 5,7,7, 5,1,3,9,4,6,8 (zeros are excluded). A 1, 2, or 3 changes x; 4,5, or 6 changes y; 7,8, or 9 changes z, so the changes are to be y, z, z, y, x, x, z, y, y1z. Then starting a t the origin (0,0,O)with vector (1, 1, l), the successive vectorsare:(1,1, ~ ) , ( l , - ~ , l ~ , ( ~ , - l , - l ) , ( l , - l , l ) , (-1, l,l,l),( l , l ) , ( l , l , V , (-1, l,-l), (-1,-1,-1), (-1,1,-1),(-1, Ll). Thesuccessivelatticepointsare: (O,O,O), (1,1, l),(2,0,2), (3,-1, 1),(4,-2, 2),(5,-11,3),(4,0,4), (5,1,5),(4,2,4),(3,1,3),(2,a,2), and (1,3,3), the last value being the position of the particle after
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
November 1951
0 SAMPLE OF 1000 0 BIASSED SAMPLE OF 1000
0 EXCLUDED
VOLUME
8
tb 0% 0
ch 0 Ql 0
troduced in the Monte-Carlo method where they would be difficult in analysis. If the tube is considered to be a cylinder with absorbing walls, every time the particle strikes a wall, it is subjected to a stochastic process of being removed with the probability proportional to the absorption coefficient. Roughness of the walls, back diffusion, etc., can be introduced provided the fundamental laws of the phenomena are known. Statistical phenomena such as these, but of almost any degree of complexity, can be introduced in a relatively simple way by biasing the selection of random numbers. A great practical advantage of the Monte-Carlo method is that a priori information concerning the solution can usually be inserted easily. For example, if the form of the solution is known in a semiquantitative way, the original distribution can be modified to fit any such ideas, and the rate of convergence is markedly improved. In general, initial conditions can easily be imposed on the starting distribution.
CHANGES OF COMPOSITION
B
0 DO
0
Figure 1. Distribution of Radius Reached by a Random Process of 64 Steps of Length 4 3 in a Tetrahedral Lattice
dE;
Radius pis normalized by dividing by the density of partlcles in the a n n u l u ~between pa and p9 dpa is iven by w(p*) whioh i. divided by p to give a linear relation with pfon a logarithmic plot for a Gaussian curve
-
2477
10 steps. Figure 1 shows the distribution of the radius reached after making 64 steps. The plot is logarithmic, and the fit to the straight line shows how soon the normal distribution is set up.
BOUNDARY AND INITIAL CONDITIONS Other features of the physical situation, such as special boundary conditions, often cause considerable trouble in solving differential equations, whereas they can be given to the machine in a simple way. For inetance, if one end of the capillary is blockedfor example, at x = -a--so that the dye molecules are reflected from the surface, the diffusion is considerably modified. This boundary condition can be introduced into the programing of the successive steps of the calculation by requiring that when a particle moves to x = -a, it can no longer go to more negative z but is to be “reflected.” In this step, the sign indicated by the toss of the “coin” is disregarded, and z is always increased. It is inconvenient to use coins in this type of machine. Stochastic (trial and error) processes, such as the simple random walk described, are carried out by random numbers supplied to or generated by the computing machines. I n the one-dimensional case, the parity of the random number can be used to decide whether xis to be added or subtracted. But many other probable features can be introduced by random numbers. For example, let us consider a boundary condition in region x -a to x = -b, in which dye molecules are absorbed in proportion to their penetration-Le., -(a-b). A molecule arrivea at x = a c, a>c>b. A random number is chosen. If it is greater than c, the particle is allowed to continue to M u s e ; if it is less than c, i t is r e m o d Le., absorbed. Many other features of the physical problem can readily be in-
+
It is also evident that chemical reactions can take place during diffusion of the particles. This allows a method to be used in the study of flame propagation, for example. This kind of problem has been solved in this laboratory in another connection. Schroedinger’s equation, which is the basis for quantum mechanics, can be written as a diffusion equation
where Vis the potential, a function of the coordinates. The interpretation of Schroedinger’s equation in this form ( 1 ) is as follows: The particle diffuses in the normal fashion, but when it reaches the point (2, y, z...), instead of remaining 8 s one particle, it now propagates itself and forms e- vAt particles. This, of course, may be greater or less than unity. Thus a diffusion problem results in which the particles increase (or decrease) in concentration depending on their position in the system. This is somewhat similar to having chemical reactions proceed during the diffusion process. Another feature of the Schroedinger equation is that the wave function, u, which must be obtained in some cases, has nodes-Le., it becomes zero at certain points in space. To impose this on the machine by the Monte-Carlo method any particle is prevented from passing the place where the node is and is removed completely from the system. When this is done, the distribution of particles after a large number of steps settles down to a steady state which is the wave function of the system. This situation is exactly the same as the boundary condition of an absorbing wall. In the Schroedinger equation, solutions to the simplest problems, such as the harmonic oscillator, have been found with accuracy to a few per cent, using 1000 walks of 400 steps each. The Monte-Carlo method can be used in this way to study the diffusion of abstract particles in phase or configuration space as in quantum mechanics or as material particles in ordinary material diffusion or as in units of heat in heat flow problems.
FURTHER ELABORATION Another diffusion problem which has been studied extensively is in connection with statistical mechanics of polymers. This has been used in attempting t o find certain properties of random configurations of long chain molecules. The problem is equivalent to solving the diffusion equations of three dimensions with a peculiar condition-namely, that a particle can never go to a point in space to which i t has already been. This is the ordinary diffusion equation with a boundary condition of sinks diffusing in correlation with the particles themselves. This kind of problem has never been considered in the theory of analysis and obviously is one of great difficulty.
Vol. 43, No. 11
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
2416
Figure 1,second curve, shows the distribution of the same particles a t first curve under the condition that lattice points have been “poisoned,” preventing re-entry. The problems in physics have difficulties that are much deeper than those of engineering, whereas the engineering problems are more complex. However, large scale computing machines have considerable capacity and can handle a great deal of complexity. At present they are designed to do arithmetic only and are not particularly adapted to the higher mathematics of physics. A very realistic model of situations involving classical physics or chemistry can be inserted into the machine. The actual calculations invoived in the Monte-Carlo method are extremely elementary, involving only many additions of unity in certain counters chosen by random numbers.
If these problems of wave mechanics and statistical mechanics can be solved by the Monte-Carlo method, certainly the classical diffusion problems in chemical engineering are subject to attack by the same method. LITERATURE CITED (1) King, Gilbert W., IBM Seminar on Scientific Computation,
“Stochastic Methods in Quantum Mechanics” (November
1949). (2) King, Gilbert W., Office of Naval Research, Contract N6-onr228, T.O. 11, Tech. Rept. No. 1. (3) Ibid., Tech. Rept. No. 5. (4) King, Gilbert TV., U. S. Dept. Commerce, Natl. Bureau of Standards, Applied Mathematics Series 12.
RECEIVED February 9, 1961.
Analogical Computing Devices in the Petroleum Industry WILLIAM L. MORRIS Phillips Petroleum Co,, Bartbsoille, Okra.
T
HERE are three types Three of the numerous analogical computing devices -that is, the one whose mole now at work in the petroleum industry are the vaporfractioninthestreamiszl. of problems to which analogical computing d e liquid equilibrium computer, the Phillips “66” spectrois the absorption coefficient of the second hydrocarbon for computer, and the reservoir behavior analyzer. vices are particularly suited: 1. Problems which require These devices give results usually of greater precision this frequency of infrared than the data supplied to them; they are simple to radiation; its mole fraction in a high speed in attainingsolutions to problems if the sohoperate and maintain and are inexpensive. the stream is x2, and soon; Computing devices do not have to be of universal applitions are to be useful. b1 is the total extinction. The 2. Recurrent problems for cation, nor of tremendous size and speed to give important second equation is obtained by measuring the absorption which direct solution is not service in industry. by the stream of another frefeasible and, in consequence, require tedious cut-and-try quency of infrared radiation. GI is the absorption coefficient by the first hydrocarbon of the processes. second frequency of infrared radiation; a22 is the absorption 3. Boundary value problems in which the boundaries and/or boundary values are either extremely complex or, as is normally coefficient by the second hydrocarbon of the second frequency, the situation in design problems, the boundaries and boundary and so on. values are the quantities to be obtained from the computation. The problem of solving Equation 1 for 2 1 , z2,xs,etc., is tedious The petroleum industry has problems in all three of these cateif there are four or more variables. Also, since normally the stream is being analyzed so as to maintain a certain composition gories and various devices have been used to solve them. This paper discusses three devices in m e in the petroleum industry of the stream, it is obvious that if the solution is to have any utility it must be available a short time after the observation. and the problem which these devices have been constructed to It must, in fact, be available before there is any appreciable solve. These are typical of each of the three categories. change in the composition. The spectrocomputer, which wm DEVICES NOW IN USE IN THE PETROLEUM INDUSTRY devised and built by Morgan and Crawford (4) of Phillips PetroSolve Systems O f simultaneous equations, falls into ]cum CO. Phillips “66” Spectrocomputer and Infrared Analysis of Hydrocarbons. The analysis of hydrocarbons by the use of in+h,si O f comPukrs. The basic Principle of the Mo~gancomputer is the Gaussfrared radiation is well known and has been found to be a highly Seidel iteration method, which consists in solving the first e q u e reliable method. To the computer, however, it poses the problem of making rapid solutions of systems of simultaneous linear tion for 21by making the assumption that all the others are mro. Let this first approximation to 51be zi(l),and from Equation 1 algebraic equations as outlined in category 1. The number of equations is the same as the number of compounds in the b Z2(O) = z p = -, zl(l) , . , = Zn(0) = 0 stream and the equations will be of the form all i
%(I),
. anlzl
.
+ anzzz + . .. + anfln=
the first approximation to zs, is obtained by putting zz(l), . . .zn@)in the second equation and solving for za(1) a(‘)E bi - ~ 1 1 ~ 1 ( 1 )
za(O), z((O),
(1) Here the first equation is the one obtained by measuring the absorption by the stream of a single frequency of infrared radiation. all is the absorption coefficient of the “first” hydrocarbon bn
a22
..
Next, TI(^), z~(l), %, 26(O), (‘I), .zn(0) are inserted in the third equation which is then solved for z*(’), and so on, This process is continued until the nth equation is reached and it is solved for
