Procs. Fifth SIAM Conference on
Parallel Processing for Scientific Computing
J. Dongarra, K. Kennedy, P. Messina, D.C. Sorensen and R.G. Voigt, Editors
pp. 382-387
SIAM, Philadelphia PA (1992)

Cellular Automata Modeling Isotropic Growth of Clusters of Arbitrary Morphology and their Application to the Study of Heterogeneous Reacting Systems

K. Zygourakis ^1,* and P. Markenscoff ²

¹ Department of Chemical Engineering
Rice University
Houston, TX 77251-1892

² Department of Electrical Engineering
University of Houston
Houston,TX 77204-4793

* Corresponding Author

ABSTRACT

This study evaluates a class of novel cellular automata models that can accurately describe structural transformations in physical systems characterized by isotropic growth of objects and complicated morphology. The algorithms designed for these models are highly suitable for parallel processing and their implementation on an INTEL hypercube is discussed. For a given number of processors, the speedup and efficiency increase with increasing size of the cellular array. This implies that in the practically interesting case of large cellular arrays, we can use efficiently a very large number of processors.

INTRODUCTION

Reacting systems involving solids whose structure evolves with time pose some very challenging theoretical problems. An example of great practical importance involves gas-solid systems where a porous solid reactant A undergoes an irreversible non-catalytic reaction with a gas B to produce gaseous product(s) C.

The porous solids often have very complicated structures. Some coal-derived chars, for example, have pores ranging from cavities a few hundred microns in diameter [1] to small submicropores whose size approaches molecular dimensions (diameter less than 0.8 nm). Figure 1 presents two binary images depicting the macropore structure of two randomly sectioned char particles we observed on a video microscope.

Since the chemical reaction occurs at the gas-solid interface, the rate R_s at which the solid reactant is consumed is proportional to the total pore surface area S_g. As the porous solid reacts, however, its internal pore structure changes and, consequently, the total pore surface area S_g and its reactivity will change drastically with conversion. Thus, measurements of intrinsic reaction rates and fundamental studies of reaction mechanisms must consider the temporal evolution of the solid structure.

CLUSTER GROWTH IN LATTICES.

The complicated internal structures of porous solids may be described using a lattice, that is an array of computational cells. Each cell corresponds to a small volume element of solid or void space. Clusters of cells corresponding to pores or solid reactant can then be distributed on the lattice using information obtained, for example, from binary images like those presented in Figure 1. As the chemical reaction proceeds, solid reactant cells are consumed or "gasified" and the pores enlarge. Thus, the temporal evolution of the internal pore structure can be modeled as a "cluster growth" problem.

Figure 1: Binary images of two coal particle cross-sections digitized on a video microscopy and image analysis system. The white areas correspond to the solid reactant (carbon) and the black internal areas are pore profiles.

We will consider gas-solid reactions occurring without diffusional limitations in a system with one homogeneous solid phase. In this case, the clusters (gas-filled pores) will exhibit isotropic unconstrained growth. That is, the clusters will grow at the same rate in all directions. Let dV be the volume of solid reacted over a time period dt. Then

(1)

where is the solid density, m_s(t) is the mass of solid at time t, S_g(t) is the specific pore surface area exposed to the gas and R_s is the intrinsic reaction rate. If R_s is constant over the entire surface area S_g, a solid layer with uniform thickness dz will be consumed in the time interval [t, t+dt]. Since

dV = ms(t) Sg(t) dz

equation (1) becomes

(2)

Equation (2) is the basis for relating simulation results to experimental data [2]. At every simulation step, we react a layer with uniform thickness from the solid surface area exposed to the gas. Then, every step of the simulation will correspond to a time interval

The requirement of isotropic growth precludes the use of diffusion-limited aggregation (DLA) type models [3, 4] for our problems. Time normalization and computational efficiency considerations make the Eden model inapplicable to our cluster growth problem [5, 6].

CELLULAR AUTOMATA (CA) MODELS

In order to overcome these problems, we have developed [7] a new simulation method that employs cellular automata and parallel discrete iterations. We will briefly describe here the 2-state/8-neighbor automata that model the temporal evolution of the porous structure of a solid undergoing an irreversible non-catalytic reaction. The simulation are carried out on cellular arrays with N=N_x N_y interacting cells. The state x_i of a cell can take two values.

At time t=0, the "solid" cells neighboring on "pore" cells will be exposed to the gaseous reactant contained in the pores. The exposed solid cells may then be consumed and the state of the reacted cells will change from 0 (solid reactant) to 1 (gas product). The rules of the automaton will determine which of the cells that are candidates for reaction will actually be consumed. The cellular automata evolve by simultaneously updating the state of all the cells at discrete time levels tr (where r=0,1,2,..., t0 = 0 and t^r+1 = t^r + dt) according to

(3)

where j ranges over the eight neighbors of cell i, a_ij are the interaction (or synaptic) weights and b_i is a threshold associated with the cell i. When all cells have the same synaptic weights, the iterations can be completely defined by the threshold vector b and the array W of synaptic weights

In an earlier paper [7], we have solved the problem of selecting appropriate neighborhoods, synaptic weights and cell thresholds for our discrete iterations. These parameters were selected so that the cellular automata exhibit isotropic unconstrained growth. In order to achieve this goal, the parallel iterations were implemented by using alternating weight kernels and thresholds. Let W₁ and W₂ be the following arrays of synaptic weights:

Then the parallel iterations proceed by performing:

1. m iterations with synaptic array W₁ and b_i = 0, for all i;
2. n iterations with synaptic array W₂ and b_i = 0, for all i; and
3. k iterations with synaptic array W₁ and b_i = 2, for all i.

Note that m and n are small integer numbers and k is either 0 or 1. Iteration (3) above is performed last in every period and will change the state of a cell only if that "solid" cell has at least three "gas" neighbors. This type of iteration does not propagate vertical or horizontal straight boundaries and leads to a "rounding" of the growing regions.

The automaton with m=1, n=1, k=0 grows single cells into octagonal domains. We will refer to these automaton as P110 to indicate the number of iterations of type1, 2 or 3 performed in each period. The automata P100 and P010 grow single cells into diamonds and square regions respectively. Finally, the automata P111 and P331 grow single cells into dodecagonal and decahexagonal domains respectively, thus providing increasingly better approximations to isotropic growth behavior.

COMPARISON OF CELLULAR AUTOMATA (CA) AND ANALYTICAL MODELS

The cellular automata models are the only ones capable of describing the temporal evolution of complicated porous solid structures like the ones shown on Figure 1. In order to validate these models, however, we carried extensive comparisons of their predictions to those obtained by the few available analytical models [2, 8] that describe the temporal evolution of much simpler pore structures obtained by overlapping regular geometric entities (e.g. cylindrical or spherical pores randomply distributed in a solid). Figure 2A shows such a comparison by plotting the normalized reaction rate R_s0 vs the solid conversion defined by

where V_s0 is the volume of unreacted solid at the beginning of the reaction and V_s(t) is the volume of unreacted solid at time t after the start of the reaction.

Figure 2: Solid reaction rate vs. conversion patterns for (A) a pore structure consisting of overlapping cylindrical pores and (B) the char particle whose pore structure is shown in the binary image of Figure 1B.

 

The run of Figure 2A is for a solid whose pore structure is modeled by parallel overlapping cylinders. This run was carried out on a 1024x1024 cellular array containing 1,049 overlapping circular pore profiles and a pore volume fraction equal to 0.001. Although the simulation results from the P100 and P010 automata (diamond and square growth patterns respectively) deviate significantly from the analytical predictions, Figure 2A shows that the P110 and P111 automata can approximate very accurately the temporal evolution of this pore structure. Indeed, the approximation error of the P110 automaton is so small, that we have chosen to use this automaton for most of our simulations. After establishing the accuracy of the cellular automata models, we have used them to predict the temporal evolution of reactivity for the solid reactants whose complicated pore structures were analyzed on our image processor system. Figure 2B shows the normalized rate vs. conversion pattern for the cross-section shown in Figure 1B. This simulation was carried out on a 4096x4096 cellular array.

PARALLEL IMPLEMENTATION OF CA ALGORITHMS

The cellular automata algorithms were implemented on an INTEL hypercube. The array was partitioned into subdomains and each subdomain was assigned to a different processor. The following two domain decomposition schemes for a cellular array with N = N_x N_y cells can be easily implemented on a hypercube of dimension d:

1. Linear chain: the cellular array is partitioned into K=2^d subdomains (see Figure 3A) ; and
   
2. Mesh: The cellular array is divided into K=K_x K_y subdomains (Figure 3B), where K_x = 2^dx and K_y = 2^dy are the number of subdomains in the x- and y-directions and d_x, d_y are integers such that d=d_x+d_y.

Figure 3: Linear chain and mesh partitioning schemes. Shaded areas indicate amount of data communications among processors.

Let us assume an NxN cellular grid and K processors with K = k². For the linear chain scheme, the amount of data communication among adjacent subdomains is proportional to N and independent of the number of processors, yielding communication and execution costs for each processor of O(N) and O ( N² / k² ) respectively. For the mesh partitioning scheme on the other hand, the amount of data communication among adjacent subdomains decreases with increasing k and the communication and executions costs for each processor are O(N/k) and O ( N² / k² ) respectively. If we define R as the ratio of the communication cost over the execution cost of a processor, we have

R = O ( k/N ) for the mesh partitioning scheme and

 

R = O ( 1/N ) for the linear chain scheme.

For fixed N, the mesh partitioning scheme exhibits better performance than the linear chain scheme as the number of processors K increases. Finally, for fixed K << N the efficiency of both schemes will increase with increasing size N of the cellular array.

These results were verified by simulation on an INTEL IPSC-1 hypercube. Table 1 shows the speedup we obtained for various sizes of the cellular array as a function of the number of processors p (p=1,2,4,8). The efficiency calculated for these runs is shown on Table 2.

TABLE 1

Speedup obtained for an INTEL IPSC-1 hypercube

# of Processors	32x32 array	64x64 array	128x128 array	256x256 array
1	1.0	1.0	1.0	1.0
2	1.75	1.78	1.89	1.97
4	2.33	2.73	3.35	3.76
8	3.50	4.10	5.61	6.88

TABLE 2

Efficiency obtained for an INTEL IPSC-1 hypercube

# of Processors	32x32 array	64x64 array	128x128 array	256x256 array
1	1.00	1.00	1.00	1.00
2	0.88	0.89	0.95	0.99
4	0.58	0.68	0.84	0.94
8	0.44	0.51	0.70	0.86

CONCLUSIONS

We developed and evaluated a class of novel discrete models that can accurately treat structural transformations in physical systems (e.g. reacting solids) characterized by isotropic growth of objects (e.g. pores) and complicated morphology. The cellular automata simulations provide structural information for these systems that cannot be obtained via standard analytical techniques. The algorithms designed for these models are highly suitable for parallel processing and were implemented on an INTEL hypercube. For a given number of processors, the speedup and efficiency increase with increasing size of the cellular array. This implies that in the practically important case of large cellular arrays, we can use efficiently a very large number of processors.

REFERENCES

1. G. BALLAL and K. ZYGOURAKIS, Evolution of Pore Surface Area During Noncatalytic Gas-Solid Reactions: 2. Experimental Results and Model Validation., Ind. Eng. Chem. Research. 26 (1987) pp. 1787-1796.
   
2. K. ZYGOURAKIS and C.W. SANDMANN, Pore Structure Evolution During Noncatalytic Gas Solid Reactions: 3. Development and Applications of Discrete Structural Models, AIChE J. 34 (1988) pp. 2030-2040.
   
3. T.A. WITTEN and L. M. SANDER, Diffusion-limited aggregation., Phys. Rev. B. 27 ( 1983) pp. 5686-5697.
   
4. P. MEAKIN, Diffusion-controlled cluster formation in 2-6-dimensional space, Phys. Rev. A. 27 (1983) pp. 1495-1507.
   
5. M. EDEN, Fourth Berkeley Symposium on Mathematical Statistics and Probability. IV (1961) p. 223.
   
6. H.P. PETERS, D. STAUFFER, H.P. HOLTERS and K. LOEWENICH, Radius, Perimeter, and Density Profile for Percolation Clusters and Lattice Animals, Z. Physik B, 34 (1979) pp. 399-408.
   
7. P. MARKENSCOFF and K. ZYGOURAKIS, Parallel Iterations and Cellular Automata Models for Simulating Pore Structural Transformations, Procs. 1991 ACM Computer Conference, (1991) pp. 409-419.
   
8. G. BALLAL and K. ZYGOURAKIS, Evolution of Pore Surface Area During Noncatalytic Gas-Solid Reactions: 1. Model Development., Ind. Eng. Chem. Res., 26 (1987) pp. 911-921.

Top

[Return to List of Publications]