Ind. Eng. Chem. Res. 1996,34, 3437-3444
3437
Stochastic Simulation of Transport Phenomena Lewis E. Wedgewood* and Kevin R. Geurts Department of Chemical Engineering, University of Washington, Box 351 750, Seattle, Washington 98195-1750
In this paper, four examples are given to demonstrate how stochastic simulations can be used as a method to obtain numerical solutions to transport problems. The problems considered are two-dimensional heat conduction, mass diffusion with reaction, the start-up of Poiseuille flow, and Couette flow of a suspension of Hookean dumbbells. The first three examples are standard problems with well-known analytic solutions which can be used to verify the results of the stochastic simulation. The fourth example combines a Brownian dynamics simulation for Hookean dumbbells, a crude model of a dilute polymer suspension, and a stochastic simulation for the suspending, Newtonian fluid. These examples illustrate appropriate methods for handling sourcdsink terms and initial and boundary conditions. The stochastic simulation results compare well with the analytic solutions and other numerical solutions. The goal of this paper is to demonstrate the wide applicability of stochastic simulation as a numerical method for transport problems. 1. Introduction
The transport equations, as presented in Bird, Stewart, and Lightfoot (1960),treat the conserved quantities as continuum variables, e.g., concentration, temperature, pressure, and velocity. These equations, which govern mass, energy, and momentum transport processes, can be recast as stochastic equations which in turn can be solved numerically using stochastic simulation techniques. For example, Laso (1994)solves several one-dimensional energy transport problems using stochastic simulations. The stochastic simulation method requires that the conserved properties be discretized into small “packets” so that convection and diffusion processes can be simulated by the motion of these packets. Although each of these packets represents a small quantity of mass, energy, or momentum, the simulations do not model the system at the molecular level. Only transport properties such as diffusivity, thermal conductivity, or viscosity need to be specified in the simulations, and all the limitations of the continuum approximation apply. Transport properties and variables can be related to a molecular model by kinetic theory. For example in polymer kinetic theory, the polymer molecule can be modeled using the so-called “dumbbell”models or beadspring and bead-rod chains suspended in Newtonian solvents. Brownian dynamics is a stochastic technique which has been successfully used to simulate kinetic theory models in homogeneous flows (e.g., Zylka and Ottinger, 1989;van den Brule, 1993;Wedgewood, 1993). Recently, there have been efforts to extend Brownian dynamics to nonhomogeneous flows. This requires satisfying the momentum balance for the solvent using a numerical technique, such as numerical integration (e:g., Biller and Petruccione, 1987; Petruccione and Biller, 1988a) or the finite element method (Laso and Ottinger, 19931,or by stochastic simulation (Geurts and Wedgewood, 1994). In this paper we first consider three simple examples covering each of the transport properties. Two of the presented solutions are steady state while one is transient, and several boundary conditions are considered.
* To whom correspondence should
be addressed. E-mail:
[email protected]. 0888-5885/95/2634-3437$09.~0l0
The examples illustrate some of the features of using stochastic simulation techniques to solve transport problems, and several issues arise, such as the appropriate algorithms for boundary conditions and for source/sink terms. Finally, in a fourth example we demonstrate how Brownian dynamics can be coupled with a stochastic simulation of a Newtonian suspending fluid in order to simulate a dilute polymer suspension in plane Couette flow. 2. Stochastic Representation of the Transport
Equations The transport equations are generally in the form of convectioddiffision equations. The Fokker-Planck equation (Gardiner, 1990) is a linear convectioddiffusion equation of the form
aflr,t>- -VA(r,t) -at
+
1 flr,t) sVVB(r,t) flr,t) (2.1)
where flr,t) represents the transport variable. The quantity A is a vector, and B is a second rank tensor and is often referred t o as the diffusion tensor. Many transport problems are governed by an equation of the form of eq 2.1. The Fokker-Planck equation is, however, linear in the unknown flr,t) which is not generally true of the governing equations for all transport problems. We accept this limitation in the present work and consider only linear problems. There is a well-established equivalency between eq 2.1 and either stochastic differential equations (SDE) or the master equation (ME). For example, the stochastic process described by eq 2.1 is equivalent to the SDE d d t ) = A(r,t) dt
+ D(r,t).dW(t)
(2.2)
where the tensor D is defined by the decomposition
B =D-D~
(2.3)
and dW describes a Wiener process and has the “white noise” characteristics (dw) = 0 (dWdW) = 6dt
0 1995 American Chemical Society
(2.4)
3438 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 I .O
x=L
x=o
0.8
- Analytical Solutions 0.6 @
kvdi
Figure 1. Geometry for two-dimensional heat conduction in a thin square slab.
where 6 is the unit tensor and the angular brackets denote an ensemble average. Equation 2.2 can be integrated over finite time steps, yielding a discreteform solution for the SDE (Kloeden and Platen, 1992). Averages over a large number of consecutive realizations of this solution in a small region about r give a time evolution for the property flr,t). In the following sections, the stochastic simulation method is applied to transport problems involving each of the three conserved quantities. 3. Two-DimensionalHeat Conduction-A Dirichlet Problem
The solid, rectangular slab shown schematically in Figure 1 occupies the region 0 Ix IL , 0 Iy IH,and 0 Iz I2. The thickness 2 is taken to be very small so that the temperature is uniform in the z-direction and the problem is effectivelytwo-dimensional. The density e, thermal diffusivity a,and specific heat capacity C, of the solid are constant. The edge of the slab for which x = L is held at a temperature of Tw. The other three edges are held a t temperature TO. We wish to find the steady-state temperature distribution in the slab. The transient energy balance equation for this problem simplifies to
(3.1) Thus energy accumulation is balanced by conduction in two directions. In eq 3.1 the energy density has been rescaled so that
This form emphasizes that the quantity being conserved is the thermal energy per unit volume of solid. The definition of 0 in eq 3.2 allows us to treat three of the four slab edges as absorbing boundaries. The steady-state problem has a well-known solution in the form of an infinite series, eigenvalue expansion given by Boyce and DiPrima (1986).
o=
c-
n=1,3,5,...
nnL nn sinh(
7)
(3.3)
XIL
Figure 2. Temperature profiles for two-dimensional heat conduction in a square slab. Results of a stochastic differential equation simulation are compared with the steady-state analytical solution at various y-locations.
The results of a stochastic simulation of this example will be compared below with eq 3.3 in Figure 2. Equation 3.1 is a Fokker-Planck equation (FPE) for the distribution of thermal energy in the system. In order t o simulate this heat transfer problem, the FPE is written as the corresponding SDE. Since the FPE is two-dimensional, the SDE is a vector equation, here written in component form.
&(t) = 6 a dW,(t) dy(t) = 6dWy(t)
(3.4)
Together, the above equations define the trajectories of the thermal energy packets in the simulation ensemble. The SDE simulation proceeds as follows. First, the domain is divided into N, equally sized bins in the x-direction and Ny equally sized bins in the y-direction. The bins, as explained below, are only necessary for generating a histogram of packet locations. Then a large number of thermal energy packets are positioned in the domain with an arbitrary initial distribution of locations. Now the trajectories of all the thermal energy packets in the simulation ensemble are advanced forward one step in time by numerical integration of eq 3.4. Next, the boundary conditions must be satisfied. Any packet whose trajectory reaches a boundary where 0 = 0 is removed from the ensemble. Managing the boundary where 0 = 1 is more involved. An energy reservoir region is created outside of the problem domain and adjacent t o the hot edge of the slab as depicted in Figure 1. The properties of the reservoir are crucial to the behavior of the simulation. It is not sufficient merely to keep the total number of packets in the reservoir constant; rather it is necessary that the packets be uniformly distributed in the reservoir region as well. Thus, at the end of each time step of the simulation, after each packet is allowed to diffise, all existing packets in the wall reservoir are removed from the simulation ensemble. Then the reservoir is reconstructed by adding
(3.5) energy packets to the ensemble with uniformly distrib-
Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 3439 uted, random locations inside the reservoir. This procedure ensures that the packet density in the reservoir is truly uniform. The quantity E represents the number of packets in a single unit of thermal energy, and is a simulation parameter which is chosen large enough t o keep statistical noise low, but small enough to give acceptable computer run times. The width of the wall reservoir is given by Axwall. Many cycles of moving all packets and enforcing the boundary conditions brings the system to a steady distribution of packet locations. The thermal energy density Oij in any bin in the slab may be calculated from a histogram of packet locations, as it is proportional t o the number of packets present in the bin, nij.
of the droplet is ignored. The two boundary conditions are x2%=0
ax
at x = O
and
where the Biot number is defined by Bi =
[email protected] first condition is drawn from the symmetry of the geometry and the second condition defines the mass transfer coefficient, k , (see Bird et al. (1960), eq 17.12). The steady-state solution to this problem can be obtained by the method of separation of variables (Carslaw and Jaeger, 1959)
(3.6) The index i sums over x-partitions while the indexj sums overy-partitions of the slab. Now, eq 3.2 provides the average temperature in each bin. Results from a stochastic simulation of this heat transfer problem for a square slab are plotted in Figure 2, along with the analytical solution, eq 3.3. The plot shows temperature profiles along several lines of constant y. A total of 10 000 packets was maintained in the hot thermal reservoir, and the domain was divided into 50 equally sized bins in each coordinate direction. The reservoir region is the same length as an interior was bin. A dimensionless time step size of 1 x used, and the simulation was carried out for 13 000 time steps, which is sufficient to achieve steady state if the first term of eq 3.3 is used as the initial condition. (This initial condition was chosen to minimize computation time, although other initial conditions would lead to the same steady-state solution.) The final packet distribution was averaged over the final 1000 time steps to generate the temperature profiles in Figure 2. The statistical noise in the simulation results could be reduced by increasing the number of packets in the thermal reservoir, and thereby increasing the size of the simulation ensemble. 4. Reactionhliffusion in a Droplet
A liquid droplet of radius R and consisting initially of pure component A is surrounded by a soluble gas B at a concentration CB, The diffisivity of B into A is given by @AB. We assume first-order kinetics for B absorbed in the liquid phase with a reaction rate constant k . By symmetry, the concentration of B in the liquid phase CB depends only on the radial position r. Assuming negligible flow, the dimensionless diffusion equation which governs this problem is
.
where the dimensionless parameter and variables are defined as
where Da is the Damkohler number. To obtain a pseudo-steady solution, the effect of the accumulation of B and the reaction product on the size
and is shown in Figure 3. The diffusion equation (4.1)can alternately be expressed as the SDE,
dx = 2dt X
+ 4 d@
(4.5)
where the Wiener process has been made dimensionless according to dW = d W G I R . This SDE represents the diffusion of B in the droplet but does not account for the disappearance of B caused by chemical reaction. To simulate the reaction process, each packet has a probability of reacting given by
for a first-order reaction. Thus each time step At every packet in the system has probability h of being eliminated from the ensemble. Petruccione and Biller (1988b,c)and Wedgewood and Geurts (1995)have used a similar procedure for the destruction of network strands in network theories of polymer melts. To complete the simulation scheme, it is necessary to maintain the flux of packets constant at the droplet interface according to eq 4.3 and t o maintain a zeroflux condition at the droplet center. It should also be mentioned that because of the change of coordinates this simulation produces a contracted distribution y(x) which is related to the concentration of B by
The function y(x)& is proportional to the mass of B found in a spherical shell of thickness dx and located at the radial position x . The results of this simulation are shown in Figure 3 as open circles. The simulation results compare well with the analytic curve except a t small values of x . The volume of spherical shells decreases with decreasing radial position, and simulation results a t positions near the center of the droplet are thus necessarily based on relatively few observations. This explains the large error in the results nearest the center of the droplet. This error can be reduced, however, by either running
3440 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
a longer simulation and time averaging or by using a larger ensemble of packets.
,
1.0
i 11
0.9
5. Start-up of Poiseuille Flow To demonstrate a transient solution, we consider a Newtonian fluid flowing in a slot driven by a constant pressure drop of (Po - PL)/Lin the +x-direction. In this unidirectional flow, continuity is trivially satisfied, and the Navier-Stokes equations simplify to
I
0.8
0.7
0
n
Stochastic Simulation
- Analytical Solulion
2 0.5
0.4
0.3
for the x-component of momentum. In this example, the accumulation of fluid momentum is governed by momentum diffusion and by the pressure gradient which drives the flow. The pressure gradient has units of momentum per volume per time, and thus may be thought of as a source term for fluid momentum. It is analogous to a zero-order reaction rate term in a mass transport problem. If there is no slip a t the plates and then a t steady state the wellthe slot height is W , known solution to eq 5.1 is
0.2 0.0
0.2
0.6
0.4
I .o
0.8
rlR
Figure 3. Reactioddiffision in a droplet. Results of a stochastic differential equation simulation are compared with the steadystate analytical solution. The increasing noise at the center of the droplet is caused by the decreasing size of a volume element with decreasing radius. - Analytical Solutions o Stochastic Simulation qi/plp"=0.06 Stochastic Simulation rMplp"4.l A Stochastic Simulation Ilr/plp"=0.2 + Stochastic Simulation rNplp"d.4
-I
0.4
so that the momentum profile is parabolic across the slot. The fluid momentum is zero a t each plate and is a t its maximum a t the center of the slot. For the transient start-up problem, the fluid is initially at rest until at t = 0 when a pressure drop of (PO- PL)/Lis instantaneously applied. This start-up problem has the solution
0.2
0.0
-0.2 .0.4 .0.6 -0.8
.
.1.0
1 c (rzn 1 (-1Y
1-E/-4
,,,({[2n
0.0
lIn}{;}) X
n=O
1 J
which is shown in Figure 4 as solid curves. The trajectories of packets of fluid momentum can be represented by the stochastic differential equation d y = E d W
(5.4)
which is simply a Wiener process in which the amplitude is determined by the kinematic viscosity. The creation of fluid momentum by the pressure drop term in eq 5.1 must also be accounted for. At the beginning of each time step a set number of new momentum packets, np,is added to the simulation ensemble, (5.5)
0.1
0.2
0.3
0.4 0.S 0.6 24PV/Plp"( P,-PJ
0.7
0.8
0.9
1.0
Figure 4. Start-up of Poiseuille flow. Results of a stochastic differential equation simulation are compared with the analytical solution at four different times soon after the pressure drop is imposed.
where 2 is the slot width. These packets are placed at random y-locations in the domain such that the distribution of newly created packets is uniform across the slot. Creating and locating packets in this manner is consistent with the zero-order, spatially uniform nature of the pressure gradient in this problem. The quantity E is here the ratio of the number of packets per unit of momentum. The momentum in a bin is proportional to the number of packets in the bin. Transient simulation results for early times are plotted in Figure 4 in comparison to the analytic was used solution, eq 5.3. A time step size of 1 x for these simulations, and np was set at 250. We see that the stochastic simulation method satisfactorily approximates the unsteady solution. 6. Plane Couette Flow of a Suspension of
Hookean Dumbbells A dilute polymer solution, modeled as Hookean dumbbells suspended in an incompressible Newtonian solvent, is placed between two infinite flat plates. Both plates are moved with speed V, but in opposite directions. Also, the dumbbells are allowed to interact with
Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 3441
v+
1.0
,
I
I
,
- E- : + ;
-L% y A
+V
0
Figure 6. Schematic of plane Couette flow of a suspension of Hookean dumbbells. Geometric parameters are defined, and a sketch of the expected nonlinear velocity profile is shown.
the walls. The system is shown schematically in Figure 5, which also illustrates the general shape of the expected velocity profile. This one-dimensional flow problem has previously been studied by several researchers. For example, Goh et al. (1985) assumed a linear velocity profile which was not modified by the dumbbells and applied both Brownian dynamics and continuum theory. Biller and Petruccione (1987)relied exclusively on Brownian dynamics and accounted for the effect of the dumbbells on the solvent flow field. Bhave et al. (1991) took a phase-space kinetic theory approach and introduced assumptions about the orientation of dumbbells near the walls which were not required by the Brownian dynamics approach of Biller and Petruccione, whose method is adopted for the present work. 6.1. Problem Description. The continuity equation is trivially satisfied for this flow, and the equation of motion is
where 7;ls is the viscosity of the solvent and Q ,, is the polymer's contribution to the shear stress in the fluid. Thus the accumulation of fluid momentum is governed by momentum diffusion and by the gradient of the shear stress due t o the polymer. If interactions between the walls and the dumbbells are ignored, then the Hookean dumbbellNewtonian solvent system is equivalent to using the convected Jefieys model (Bird et al., 1987) to determine the stress tensor. Under these circumstances the steady-state velocity profile is linear as it is in the absence of polymer, although the force required to move the walls is increased by the polymer contribution to the shear stress. In the following, the walY polymer interactions are treated as hard-core repulsions. Repulsive interaction between the walls and the beads of the dumbbells leads to a nonhomogeneous flow system. The number density of polymer molecules is lower near the walls than it is in the center of the gap, as is illustrated in Figure 6. Also, molecules near the walls are more likely to be oriented parallel to the walls than t o be oriented perpendicular to the walls. As a result, the polymer contribution to the shear stress in this test problem varies with 9-location in the gap, rather than being constant everywhere as it would be if the Jefieys model were used. This leads to the nonlinear velocity profile sketched in Figure 5. At steady state, the momentum profile must be nonlinear to compensate for 9-variations in Qpyx. Biller and Petruccione (1987) calculate polymer contributions to the stress tensor with a Brownian dynamics (BD) simulation. They assume initially that the velocity profile is linear and run a BD simulation to calculate tpyx.The momentum balance, eq 6.1, is
Simulation BNnn dr Orisan (1985)
-0.2
dnJn,
Figure 6. Dumbbell bead density profile. Results from the Brownian dynamics simulation are compared with an analytical result from the literature. The dashed line shows the average density across the slot.
integrated numerically to get a more consistent velocity profile. Then, another BD simulation is run using the new velocity profile. This iterative process continues until t,, and ir, no longer change significantly. Rather than integrating the equation of motion, we couple a master equation (ME) simulation of the momentum balance with a BD simulation similar to the one described by Biller and Petruccione. The solution of this problem demonstrates that a ME treatment of the momentum balance can work not only for Newtonian fluids (Breuer and Petruccione, 1992,1993),but also for those non-Newtonian fluids that can be represented by a mechanical model from kinetic theory. 6.2. Brownian Dynamics Simulation. To start the solution algorithm, a linear velocity profile (Newtonian case) across the slot is assumed, and a BD simulation is run to provide a first estimate of %p,z/@. The SDE for a Hookean dumbbell is
[
d$ = -
+ [RBI}d2 + 2
,&?
dW (6.2)
where H is the dumbbell spring constant, 5 is the friction coefficient between the dumbbell beads and the suspending solvent, k g is the Boltzmann constant, T i s the absolute temperature, is the dumbbell connector vector, R is the transpose of the velocity gradient tensor, and d\qTis a vector of increments of independent Wiener processes. Equation 6.2 is integrated to produce dumbbell trajectories from which a shear stress profile is extracted by time and spatial averaging. The polymer shear stress profile from the BD simulation is then utilized in a ME simulation of the momentum balance, eq 6.1. This simulation is described next. 6.3. Master Equation Solvent Simulation. Dimensionless variables are introduced,
where no is the nominal polymer number density. This leads to a dimensionless momentum balance (6.4)
3442 Ind. Eng. Chem. Res., Vol. 34,No. 10,1995 M
containing one dimensionless parameter GZ
wT=cwi
n&JB G, VSV
(6.5)
which compares the magnitude of polymer stresses to stresses in the Newtonian solvent. The domain -B Iy I +B is spanned by M node points each spaced Ay apart. In a ME simulation, momentum packets exist only at the nodes, so that it is only necessary to keep a count, Ni, of the number of packets at each node. The diffusion of momentum consists of packets jumping between neighboring nodes. When a packet leaves a node, it carries its momentum (and velocity) with it. In this formulation, the displacement of packets is fured, but there is a distribution of time intervals between the jumps. If the packet count at a node is negative, then the velocity at that node is negative. Each particle carries some quantity 6u of positive or negative velocity which must be fixed for each simulation. A master equation for the many-particle probability distribution P(Ni,t) can be written for the system ( c f , Breuer and Petruccione, 1992)
The notation E indicates a step operator as defined by Van Kampen (1981). The quantity f l is equal to Ni if Ni is positive, and it is zero otherwise. Conversely, is equal to Ni if Ni is negative, and it is zero otherwise. The term -G2(atpyx/?y)in eq 6.4is treated analogously to the pressure gradients in Breuer and Petruccione (1992); that is, it is treated as a momentum creation term. In this case, though, the creation rate varies with y instead of having a constant value throughout the region between the walls. Above the center line of the domain, -azp,yx/?y < 0 so that the creation rate is negative. Thus in that region negative x-momentum is created. Conversely, in the lower half of the problem domain, the creation rate is positive and positive x momentum is created. A guess is made as to the number of velocity particles a t each node, and then the simulation is started. First, a transition rate wi is calculated for each node according to the expression
which follows from eq 6.6. It represents the probability per unit time of a transition occurring a t the ith node of the system. These individual rates are then summed to get the total transition rate W T for the entire system, which indicates the probability of a transition occurring at any node.
(6.8)
i=l
Next, the waiting time t* is determined. It is the random elapsed time between system transitions, and it depends on the total transition rate and on a real random variable 5 which is uniformly distributed on the interval [0,11. (6.9) Now it must be determined a t which node the transition will occur. The relative probability of the transition occurring at any single node is given by the ratio of the node transition rate to the total transition rate, or WJWT. A single node is then selected such that probability of that node being selected meets the calculated relative probability. One way to do this is with a rejection technique (Press et al., 1989). After selecting a node for the transition, the next step consists of choosing and implementing a specific transition involving that node. Some transitions are related to momentum diffusion and some are related to polymer shear stresses. The possible specific transitions and their relative probabilities are listed in Table 1. One of the transitions in Table 1 is randomly selected and implemented; that is, the specified changes are made in the Ni. The column marked "condition" provides the circumstances under which each transition is possible. An impossible transition has zero probability. The boundary conditions are next implemented by restoring N1 and NMt o the following specified integer values.
N l = -1Idv
(6.10)
NM = 1I6v
(6.11)
This completes one waiting time interval of the master equation solvent simulation. Then the process of selecting and implementing transitions is repeated for another interval until the total time simulated is sufficient t o obtain a next estimate of the velocity profile. At the end of the simulation, the velocity at each node is determined by the conversion vi = NiSu
(6.12)
To decrease statistical noise in the final results, often multiple simulation realizations are performed, each differing by the string of pseudo-random numbers used t o select transitions and calculate waiting times. Then final velocities are calculated as ensemble averages of all such realizations. This also allows for the calculation of standard deviations so that the statistical variance may be quantified. The velocity profile so calculated is a first improvement on the linear profile assumed a t the start of the solution algorithm. 8.4. Coupling of Simulations. The new velocity profile is utilized in a new BD simulation to provide a next estimate of &pJz/2y. Information is thus fed back and forth between the BD simulation of the polymer contribution to the stress and the ME simulation of the momentum balance. The polymer simulation passes stresses to the solvent simulation, which in turn passes velocities back to the polymer simulation. When both simulations provide equivalent (to within a certain
Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 3443 Table 1 type
transition
probability
diffusion
Ni (Ay12wi
diffusion
Ni
diffusion
diffusion
1.0
condition
Ni
2
0
Ni
2
0
Ni 5 0
Ni 5 0
polymer stress
N,-Nifl
polymer stress
Ni-Ni - 1
A
I
0 Maslcr Equation Simulation
- Biller & Petruccione (1987) __ Newtonian _.
2.
Fluid
1
4.5p, I 0.0
,
-1.0 -1.0
-0.5
0.0
,
,
0.5
,
1.O
V
Figure 7. Velocity profile for plane Couette flow of a suspension of Hookean dumbbells. Results of the master equation simulation are compared with a profile determined by numerical integration as by Biller and Petruccione (1987). The error bars show the variation between the 20 realizations of the flow. The dashed line shows the linear velocity profile of a Newtonian fluid.
tolerance) results on successive passes through them, the steady-state solution to the non-Newtonian flow problem has been determined. 6.5. Stochastic SimulationResults. We solve this flow for the same parameters used by Biller and Petruccione (1987) to generate the dotted line in their Figure 4. The gap is made very narrow, such that the root-mean-square equilibrium length of a dumbbell is equal t o the plate spacing, 2B, to emphasize the flow inhomogeneity. The previous researcher's parameters for the plate spacing, plate speed, and solvent viscosity equate to setting G2 = 1. Other required parameters for the ME solvent simulation include the number of nodes (M= 20) and the amount of velocity represented by a single packet (6u = 0.005). Twenty separate realizations are used to reduce statistical variation. A velocity profile for this non-Newtonian flow problem is shown in Figure 7. For the sake of comparison, included is a profile obtained from the numerical
Figure 8. Polymer shear stress and shear stress derivative for plane Couette flow of a suspension of Hookean dumbbells. The circles show the polymer contribution to the shear stress as determined by Brownian dynamics simulation. The squares show the spatial dependence of the momentum creation term in eq 6.4. Positive momentum is generated below the center line, and negative momentum is generated above.
integration scheme of Biller and Petruccione. Figure 8 shows the final polymer contribution to the shear stress, and its derivative, which is required in the momentum balance, eq 6.4. This problem demonstrates that a master equation simulation is possible not only for Newtonian flows but also for non-Newtonian flows. Each BD simulation takes about 5 min to run on an HP9000/720 workstation, while each ME simulation takes about 90 min. It takes four or five iterations to obtain consistent stresses and velocities. 7. Conclusions
In this paper we have illustrated the use of stochastic simulation as a method for numerically solving transport problems through four examples of increasing complexity. We believe that stochastic simulation is an increasingly useful numerical method because it is ideally structured for running on parallel computers. Laso (1994)reports, for example, that a one-dimensional
3444 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
stochastic differential equation simulation of heat conduction runs 40 times faster on a vector machine than on a scalar machine. In the extreme, a separate processor might track each packet in a stochastic differential equation transport simulation, integrate the trajectory of each dumbbell in the ensemble of a Brownian dynamics simulation, and manage each realization of a master equation transport simulation. Although all transport properties were assumed constant in the examples presented, it is a simple matter to account for dependence on temperature, pressure, concentration, etc. in a stochastic simulation. Also, the error or uncertainty of a simulation result can be decreased by increasing the size of the simulation ensemble. This means that spatial resolution can be increased without introducing such errors as might be caused by roundoff in matrix inversion operations common in other numerical methods. In section 6 we solved a non-Newtonian flow problem without reference to a constitutive equation relating stress to kinematic variables. Instead the polymeric contribution to the stress tensor was determined from a Brownian dynamics simulation of a kinetic theory model. The advantage of such an approach is that constitutive relations can be formulated at the molecular level and then tested in complex flows without additional simplifications or approximations. The master equation simulation presented in section 6 is better suited for momentum diffusion problems than is a stochastic differential equation approach. The vector nature of momentum requires that a “sign” be attached to momentum packets so that flow need not be unidirectional. In a stochastic simulation, this requires the presence of velocity/momentumpackets and “antipackets,” and it is more efficient to allow a counter to take on positive and negative values than it is t o integrate trajectories of an ensemble of packets and another ensemble of antipackets. This technique is presently being applied to complex two-dimensional flows of dumbbell suspensions (Geurts and Wedgewood, 1994).
Acknowledgment The authors are grateful to the Shell Oil Company Development Fund and the Graduate School Research Fund of the University of Washington for financial also gratefully support. One of the authors (K.R.G.) acknowledgesthe financial assistance of the Henry Gray Fellowship.
Boyce, W. E.; DiPrima, R. C. Partial Differential Equations and Fourier Series. In Elementary Differential Equations and Boundary Value Problems, 4th ed.; John Wiley & Sons: New York, 1986; pp 573-575. Breuer, H. P.; Petruccione, F. A Stochastic Approach to Computational Fluid Dynamics. Continuum Mech. Thermodyn. 1992, 4, 247-267. Breuer, H. P.; Petruccione, F. On the Numerical Integration of Burgers’ Equation by Stochastic Simulation Methods. Comput. Phys. Commun. 1993, 77,207-218. Brunn, P. 0.; Grisafi, S. Kinetic Theory of a Dilute Polymer Solution in a Small Channel: Equilibrium Results. Chem. Eng. Commun. 1985,36,367-383. Carslaw, H. C.; Jaeger, J. S. Conduction ofHeat in Solids; Oxford University Press: London, 1959. Gardiner, C. W. The Fokker-Plank Equation. In Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 2nd ed. (Corrected Printing); Springer-Verlag: Berlin, 1990; p 144. Geurts, K. R.; Wedgewood, L. E. Stochastic Simulation and Fluid Mechanics of Polymeric Liquids: Presented at the AIChE 1994 Annual Meeting, San Francisco, CA, November 1994; paper 127c. Goh, C. J.; Atkinson, J. D.; Phan-Thien, N. The Flow of Dilute Polymer Solution in a Narrow Channel. I. The Slip Effect in Simple Shear Flow. J . Chem. Phys. 1985,82 (21, 988-995. Kloeden, P. E.; Platen, E. Numerical Solution of Stochastic Differential Equations; Springer-Verlag: Berlin, 1992. Laso, M.; Ottinger, H. C. Calculation of Viscoelastic Flow Using Molecular Models: The CONNFFESSIT Approach. J . NonNewtonian Fluid Mech. 1993, 47, 1-20. Laso, M. Stochastic Dynamic Approach to Transport Phenomena. AIChE J. 1994,40, 1297-1311. Petruccione, F.; Biller, P. A Consistent Numerical Analysis of the Tube Flow of Dilute Polymer Solutions. J . Rheol. 1988a,2 (l), 1-21. Petruccione, F.; Biller, P. Rheological Properties of Network Models with Configuration-Dependent Creation and Loss Rates. Rheol. Acta. 1988b, 27, 557-560. Petruccione, F.; Biller, P. A Numerical Stochastic Approach to Network Theories of Polymeric Fluids. J . Chem. Phys. 1988c, 89, 577-582. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Random Numbers. In Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: Cambridge, 1989; pp 203-204. van den Brule, B. H. A. A. Brownian Dynamics Simulation of Finitely Extensible Bead-Spring Chains. J . Non-Newtonian Fluid Mech. 1993,47, 357-378. van Kampen, N. G. One-step Processes. In Stochastic Processes in Physics and Chemistry; North-Holland Publishing: Amsterdam, 1981; p 145. Wedgewood, L. E. Internal Viscosity in Polymer Kinetic Theory: Shear Flows. Rheol. Acta. 1993, 32 (41, 405-417. Wedgewood, L. E.; Geurts, K. R. A Non-Affine Network Model for Polymer Melts. Rheol. Acta. 1995, 34, 196-208. Zylka, W.; Ottinger, H. C. A Comparison Between Simulations and Various Approximations for Hookean Dumbbells with Hydrodynamic Interaction. J . Chem. Phys. 1989, 90, 474-480.
Literature Cited Bhave, A. V.; Armstrong, R. C.; Brown, R. A. Kinetic Theory and Rheology of Dilute, Nonhomogeneous Polymer Solutions. J . Chem. Phys. 1991,95 (4), 2988-3000. Biller, P.; Petruccione, F. The Flow of Dilute Polymer Solutions in Confined Geometries: A Consistent Numerical Approach. J. Non-Newtonian Fluid Mech. 1987,25, 347-364. Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; John Wiley & Sons: New York, 1960. Bird, R. B.; Curtiss, C. F.; Armstrong, R. C.; Hassager, 0. Elastic Dumbbell Models. In Dynamics of Polymeric Liquids Volume 2: Kinetic Theory, 2nd ed.; John Wiley & Sons: New York, 1987; pp 71-73.
Received for review January 5, 1995 Revised manuscript received July 19,1995 Accepted July 27, 1995@ IE950028X
Abstract published in Advance ACS Abstracts, September 1, 1995. @
