Analysis of distances between inclusions in finite binary stochastic materials

Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Analysis of distances between inclusions in finite binary stochastic materials David P. Griesheimer a,n, David L. Millman b, Clarence R. Willis c a b c

Bettis Atomic Power Laboratory, Bechtel Marine Propulsion Corporation, West Mifflin, PA 15122, USA Department of Computer Science, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA Mechanical Engineering Department, University of Texas at Austin, Austin, TX 78712, USA

a r t i c l e in fo

Keywords: Stochastic mixture Monte Carlo Chord length sampling Stochastic transport Radiative transport

abstract A generalized probability density function (PDF) describing the distribution of interinclusion distances in finite, isotropic, binary stochastic materials with fixed diameter inclusions has been developed and tested. The new probability density function explicitly accounts for edge effects present in finite two- and three-dimensional stochastic materials. The generalized PDF is shown to include factors that are dependent on both the geometry of the material region as well as the statistical properties of the material. A discussion of the properties and application of this newly developed PDF is provided along with supporting numerical results for case studies in one- and two-dimensions. These numerical results demonstrate the ability of the newly derived PDF to correctly account for edge effects in finite stochastic materials, while still reproducing the expected distribution within the bulk material region. Published by Elsevier Ltd.

1. Introduction Over the past 20 years, a significant amount of research has been devoted to developing methods for solving the linear transport equation in random heterogeneous mixtures of two or more distinct materials [1–13]. The resulting methods have been successfully applied to a wide variety of problems, including the modeling of radiation transport through heterogeneous materials such as concrete, water droplets in the atmosphere, gas/dust in interstellar space, two-phase gas–liquid mixtures, biological tissue, TRISO fuel particles (commonly found in High Temperature Gas Reactor (HTGR) fuel), and certain thermal neutron poison materials used for criticality safety applications (marketed under the trade name BORALs) [3–13]. In each of these cases, the transport medium can be modeled as a stochastic mixture of homogeneous inclusions, which are randomly distributed within a uniform

n

Corresponding author. Tel.: + 1 412 4766501; fax: + 1 412 476 6305. E-mail address: [email protected] (D.P. Griesheimer).

0022-4073/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.jqsrt.2010.06.013

background material. Solving the linear transport equation for this type of mixture, referred to as a stochastic material, presents a unique set of challenges. The random distribution of inclusions within the mixture means that the properties (e.g. cross-sections) for a stochastic material have a strong spatial dependence. Furthermore, this spatial dependence is unique for every particular realization of the material, implying that the transport solution itself must be dependent on the detailed arrangement of inclusions within the material. Unfortunately, for most cases of interest, it is impractical to consider explicitly modeling every inclusion in the material. Even if the inclusions were stationary with known locations, which they typically are not, the actual transport calculation would be prohibitively expensive to solve for large numbers of inclusions. Instead, it is more practical to search for an ‘‘average’’ transport solution, which is reasonably accurate over a wide range of material realizations. Refs. [5,6] provide a thorough review of the previous methods development work in this area. In many of these previous studies, Monte Carlo methods have been employed for producing reference solutions for stochastic media transport problems.

578

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

Unfortunately, brute force approaches for calculating the ensemble average solutions with Monte Carlo methods are severely limited by the availability of computer resources. As noted earlier, the explicit representation of individual inclusions within a stochastic material is impractical for all but the smallest problems. To overcome these problems, Zimmerman and Adams [3] proposed a new technique for Monte Carlo transport calculations, which does not require inclusions to be explicitly represented. In their chord-length sampling (CLS) method, the distances between and within inclusions in a stochastic material are randomly sampled from appropriate statistical distributions during the particle transport process. By sampling the distances between material interfaces (background/inclusion or inclusion/background) on the fly, the CLS method is able to track particles through stochastic materials very efficiently. Furthermore, by only considering the distances between material interfaces that occur along the current particle flight path, the CLS method effectively reduces multi-dimensional particle tracking in a stochastic material to a series of simpler one-dimensional ray tracing calculations, which are easier to compute and analyze theoretically. In numerical experiments, variants of the CLS method have repeatedly shown excellent performance and good agreement with benchmark results [3–8,10,11]. However, while the CLS method has shown promising results for many benchmark calculations, there are still a number of open questions regarding its use in more realistic problems. While chord length distributions within slab (1-D), disk (2-D), and spherical/ellipsoidal (3-D) shaped inclusions are well known [5,11–13], the distribution of distances between adjacent inclusions along a fixed ray has not been as well characterized, especially for finite higher-dimensional problems. In infinitely thick, one-dimensional problems (with fixed-width inclusions), it is known that the distribution of distances between adjacent inclusions should be exponentially distributed with a constant rate parameter [11–14], which we will denote l. Similar results have also been derived for infinite problems in two- and three-dimensions [12,14–19]. However, for finite two- and three-dimensional problems, the actual distribution of distances between inclusions is not known, with Su and Pomraning [13] commenting that ‘‘analytic results for the [distribution] seem not to be possible’’. As a result, recent research into CLS methods has typically assumed an exponential distribution of the distances between adjacent inclusions, even for finite problems involving two- and three-dimensional inclusions. This assumption has been shown by Olsen et al. to be accurate to within 2% in the interior of dilute stochastic mixtures (containing up to 10% inclusion volume fraction) [12]. Many of the studies on the CLS method have relied on empirical fitting calculations to determine a fixed exponential rate parameter l for a particular problem of interest [4–6,11,12]. These empirically produced values have typically yielded good agreement with numerical experiments. However, requiring an empirically fit rate parameter to be calculated for each unique problem can

be viewed as a limitation of the original CLS models. These methods would be much more convenient and applicable if the rate parameter could be calculated, or at least estimated, on the fly, thus eliminating the need for auxiliary fitting calculations. To this end, recent work by Ji and Martin [10,11] has considered the distribution of distances between fuel microspheres in pebble and compact-type HTGR fuel elements. For three-dimensional problems, their result matches the theoretical distribution of distances between uniformly distributed spheres in an infinite medium [11,14]. However, subsequent analysis has shown that this limiting distribution may not be an accurate approximation near the edges of stochastic materials with finite extent [12]. This boundary effect is due to the fact that the volume fraction of inclusions is lower near the edge of a region than in the center, due to limitations on allowable particle locations near the boundary [11,12]. Previous work by Murata et al. [7,8] and Ji and Martin [11] has shown that it is possible to apply correction factors to the infinite-medium results in order to account for these boundary effects. In these studies, the boundary correction was applied through the use of a modified inclusion volume fraction, which takes into account the size and geometry of the stochastic material region. In addition, Ji and Martin [11] have developed a methodology for calculating a problem specific correction factor through the use of simple fitting calculations. In this paper we present a detailed analysis of the distribution of distances (chord lengths) between adjacent inclusions in a binary stochastic material of finite extent. This analysis is intended to account for both the random distribution of inclusions within the stochastic material, as well as the geometric constraints imposed by the outer boundaries of the finite material.

2. Derivation of the generalized probability density function 2.1. Assumptions For this analysis we will consider an arbitrary, bounded domain of finite extent in three dimensions. This volume is filled with an isotropic stochastic material, containing a number of spherical inclusions of fixed diameter Dx distributed randomly throughout. A twodimensional analog of this system is shown in Fig. 1. The background material of the volume is identified as material A, and the material in the inclusions is identified as material B. Both the number and spatial location of the inclusions is assumed to be random, subject to the conditions that no inclusion may extend beyond the boundaries of the volume and no two inclusions may overlap with one another. Before continuing, it is useful to consider the practicality of the assumptions used in this work. The assumption of identically shaped fixed-diameter inclusions is of particular relevance to engineered binary stochastic materials, such as HTGR fuel, BORALs and concrete, where the inclusion shape and size can be controlled.

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

579

2.2. Derivation of the general form for the PDF

W Δx

L Fig. 1. Illustration of a two-dimensional binary stochastic material with length L, width W, and inclusion diameter Dx.

However, this assumption is often not valid in biological, atmospheric, or other naturally occurring stochastic mixtures, where the inclusion shape and size are random variables. For one-dimensional stochastic materials, the results presented in this paper may be readily generalized to account for variable width inclusions by simply using the expected inclusion diameter, /DxS, in place of the constant diameter Dx. Unfortunately, this approach is not guaranteed to work for finite stochastic materials in twoand three-dimensions. In these cases the distribution of inclusion widths must be accounted for explicitly in the derivation of the inter-inclusion chord length distribution (specifically, the inclusion widths must be accounted for in the volume function V(x), which is described in Section 2.3). A comprehensive treatment of stochastic materials with variable-diameter inclusions is beyond the scope of this paper. Also, in many natural systems it is not reasonable to assume that the spatial locations of inclusions in the material are uncorrelated. Instead, many naturally occurring stochastic systems, such as water droplets in clouds and steam bubbles in two-phase flow, exhibit patterns of clumping among the inclusions [20]. This spatial correlation between inclusions affects the probability that an inclusion will be found in a particular location (given by function q(V), as described in Section 2.4). As noted by Olsen et al., assuming an ideally random distribution in cases where positive (or negative) pair correlations exist between inclusions is not appropriate [12]. Therefore, the results presented here will be limited to ideal (isotropic and uniformly distributed) random distributions of inclusions. With this being said, in this section we seek to derive a generalized probability density function that describes the distribution of distances that separate two adjacent inclusions along a fixed ray intersecting the volume, subject to the assumptions above. Throughout this section we will consider the detailed derivation in three dimensions, and then subsequently indicate trivial reductions necessary to reproduce the results in two or one dimensions.

To begin, let us consider a fixed point in the volume, (x0,y0,z0), located in material A. Further consider a reference ray, pointing parallel to the x-axis, which passes through the point (x0,y0,z0). The reader may consider this ray to be the direction vector of a particle or photon passing through the stochastic material. Accordingly, we will measure all distances traveled in the stochastic material along this ray, and with respect to the reference point (x0,y0,z0). This convenient choice of coordinate frame effectively reduces our problem to that of finding a probability density function in one dimension, along the x-axis. Let us now assume that the next inclusion encountered along the ray occurs at an arbitrary distance x from the reference point. From geometric considerations, we can immediately show that the center point of this next inclusion must be located on the surface of a hemispherical shell centered about the point (x0 + x,y0,z0), as illustrated in Fig. 2a. If we allow the center point of this hemispherical shell to sweep through a range of x values, from x= x0 to x =x0 +X, a three-dimensional volume is created, as shown in Fig. 2b. The resulting volume is shaped like a cylinder with convex hemispherical head, and concave hemispherical tail, both with radius Dx/2. Therefore, the volume of this region can easily be calculated to be  2 Dx : VðXÞ ¼ X p 2

ð1Þ

Again, from geometric reasoning, we can immediately see that if the center point of any inclusion falls within the volume V(X), then the distance from the reference point to the next inclusion along the x-axis must be less than or equal to X. From this observation it follows immediately that if no inclusion center point falls within the volume V(X) then the distance from the reference point to the next inclusion along the x-axis must be greater than X. We can immediately write this latter form in terms of an equality involving probabilities pðx 4 XÞ ¼ qðVðXÞÞ,

ð2Þ

where q(V(X)) is defined to be the probability that no inclusion center points are found within the volume V(X), and p(x 4X) is the probability that the distance to the next inclusion, x, is greater than the fixed value X. By inspection, we see that the left-hand side of Eq. (2) is related to the cumulative distribution function for distances to the next inclusion, P(X), which, in turn, is related to the probability density function for distances to the next inclusion, p(x). Thus, Eq. (2) can be rewritten to yield a general relationship between the desired probability density function (PDF) and the probability of finding zero inclusion center points within the volume V(X) 1qðVðXÞÞ ¼

Z

X

pðxÞ dx: 0

ð3Þ

580

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

V(X)

(x0,y0,z0)

Δx/2

(x0,y0,z0)

Δx/2 x

x (x0+X,y0,z0)

(x0+X,y0,z0)

Fig. 2. Two-dimensional illustration of the allowable inclusion center point locations that result in (a) a distance x =X from the reference point, or (b) a range of distances 0r xr X from the reference point.

To solve for the PDF, p(x), we simply take the derivative of both sides of Eq. (3) to give pðxÞ ¼ 

d dVðxÞ qðVðxÞÞ : dVðxÞ dx

ð4Þ

Note that, in Eq. (4), we have replaced the variable X by the variable x for notational convenience. The resulting PDF, given in Eq. (4), has a surprisingly simple form, especially considering the complexities of the problem being addressed. Unfortunately, this simple form is misleading, as the complexities of the problem have not been resolved, but simply hidden from view in the functions V(x) and q(V(x)). The benefit of the form given in Eq. (4) lies in its generality and flexibility. Although originally derived for three-dimensional problems, the form of the PDF remains unchanged for one- and two-dimensional problems. The only difference in the lower dimensions is that the volume function, V(x), has units of length and area, rather than volume, for one and two dimensions, respectively. Furthermore, the abstract interpretation of the function V(x) as the set of all inclusion center points that give a distance of less than x from a reference point to the inclusion along a fixed ray allows an analysis of edge effects, as shown in the next section. 2.3. Properties of the function V(x) As discussed previously, the function V(x) defines a volume (3-D), area (2-D), or length (1-D) defining the set of all inclusion center points where the distance to the surface of the associated inclusions is less than distance x from the reference point along the x-axis. From geometric reasoning, it is easy to show that V(x) is equal to 8 x, 1-D, > > > < xDx, 2-D,  2 VðxÞ ¼ ð5Þ > D x > > , 3-D, : xp 2 for constant-diameter inclusions and in the absence of any material boundary effects. It follows immediately that the derivative of Eq. (5) with respect to x is given by 8 1, > > > dVðxÞ < Dx,  2 ¼ > dx Dx > >p , : 2

1-D, 2-D, 3-D:

ð6Þ

In the vicinity of a material boundary, however, the function V(x) will have a more complicated form. To show this, consider a reference point, (xs,ys,zs), located on the bounding surface of the stochastic material and a reference ray passing through the point directed into the material. For convenience, we will again assume that our coordinate system is defined so that the reference ray is parallel to the x-axis. Further consider an inclusion whose (nearest) surface intersects the ray at a distance x from the reference point. As before, we can immediately show that the center point of this inclusion must lie on the surface of a hemispherical shell of radius Dx/2, centered at the point (xs + x, ys, zs). As we take the limit as x-0, effectively bringing the surface intersection of the inclusion along the reference ray into contact with the material boundary, we see that some of the possible inclusion center points (as defined by the hemispherical shell) will lie closer than Dx/2 from the material boundary. Inclusion center points that are positioned less than Dx/2 from the material boundary, in any direction, will result in inclusions that partially extend beyond the boundary of the material. By definition, these protruding inclusions are not allowed. Therefore, any potential inclusion center points that lie within a distance Dx/2 from the material surface cannot be included in the volume V(x). Thus, V(x), in its most general form, can be defined as the intersection of the set of inclusion center points that lie further than Dx/2 from any material boundary with the original set of inclusion center points that result in a distance less than or equal to x from the reference point to the next inclusion. Fig. 3 illustrates the shape of V(x) in the vicinity of material boundaries for several distinct cases. At this point we wish to emphasize that the function V(x) is based only on the geometry of the region containing the stochastic material, and not on the random nature of the material itself. If the detailed geometry of the region is known, it follows that V(x) can be calculated for any given reference ray passing through the stochastic material. Unfortunately, analytical forms of V(x) for finite materials cannot be easily derived, except when the bounding surface has a particularly simple form and orientation relative to the reference ray (e.g. a plane perpendicular to the reference ray). The primary difficulty in calculating V(x) for the general case lies in computing the intersection between the hemispherical head of the set of possible inclusion center points and the set of points located a distance Dx/2 from the nearest material

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

581

Fig. 3. Examples of edge effects on V(x) (gray shaded region) for several boundary surface shapes and orientations. In all cases, the reference point, 0, is located on the outer surface of the stochastic material boundary (heavy black line). For each example the extent of the surface boundary layer and edgeextent layer is shown.

boundary. For geometries with complicated bounding surfaces the analytical computation of this intersection can quickly become intractable. In these situations, numerical methods (possibly including stochastic methods) must be called upon to determine V(x). The complex region of intersection described above is also responsible for the presence of edge effects that have been observed in previous analyses of the distances between inclusions in finite stochastic materials [7,8,10–12]. This edge effect can be understood as follows. For any given reference point in space, there is a set of inclusion center point locations that will result in a given distance, x, from the reference point to the surface of an inclusion, along a given direction. As previously described, this set of inclusion center points defines a hemispherical shell of radius Dx/2 centered on the point of intersection

between the inclusion surface and the reference ray. As long as there are an equal number of allowed inclusion center points for any given distance x, it follows immediately that the derivative of the function V(x) will be a constant value, as shown in Eq. (6). However, in the presence of a material boundary, some members of the hemispherical set of inclusion center points for a given distance, x, may fall within the surface boundary layer, which is defined as the set of all points that lie within distance Dx/2 of any bounding surface of the stochastic material. Inclusion center points that fall within the surface boundary layer are not permitted, because they would result in inclusions that extend beyond the material boundary. Therefore, we see that some distances x (especially small distances measured from a point on the material surface) may have a truncated set of allowable

582

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

inclusion center points that will produce a particular inter-inclusion distance, x. More concisely, this statement corresponds to the observation that, near material boundaries, there are fewer ways to arrange an inclusion in space to produce a given inter-inclusion distance, without allowing the inclusion to overlap the material boundary. In situations where the set of inclusion center points has been truncated, it follows immediately that the derivative of the function V(x) will be smaller than the constant values given in Eq. (6). This observation provides a convenient metric for determining what range(s) of inter-inclusion distances will be affected by edge effects, as originally mentioned by Olsen et al. [12]. In general, distances between 0 and Dx/2, as measured from a reference point on the material surface, will suffer from edge effects. However, the actual extent of the edge-effect layer depends on the geometry of the bounding surfaces of the material and their orientation relative to the reference ray. The extent of the edge-effect layer is illustrated in Fig. 3 for several example cases. In addition to the fact that, within the edge-effect layer, the function V0 (x) will not be equal to the reference values for V0 (x) given in Eq. (6), we also note that within this region it is likely that the derivative will be a function of x rather than a constant. However, once outside of the edge-effect region, the derivative of V(x) will return to the reference value given in Eq. (6). Until this point we have considered the bounding surfaces of the material volume as the only source of edge effects in stochastic materials. All development has focused on an analysis of distances between a point located on the outer surface of the material and the surface of the nearest inclusion along a given inwardpointing direction vector. However, consideration of the general problem of finite stochastic materials leads us to conclude that edge effects will also occur near the surfaces of inclusions embedded in the material, even when these inclusions are far away from any material boundary. Since inclusions in the stochastic material are not allowed to overlap, by definition, it follows immediately that no two inclusion center points may be located closer than a distance Dx apart in any direction. Thus, any given inclusion will have a boundary layer, extending a

distance Dx/2 from all points on its surface, where center points of adjacent inclusions are prohibited. The presence of this boundary layer implies that the PDF describing distances between adjacent inclusions within the material will exhibit the same type of edge effects as the PDF that describes the distances between the material boundary and the first encountered inclusion. For the case involving adjacent inclusions it can be shown (Fig. 4a and b) that the edge-effect layer will range from a thickness of pffiffiffi ð 31ÞDx=2, for rays passing through the center point of the first inclusion, to a thickness of Dx/2, for rays that intersect tangentially with the surface of the first inclusion. This observation leads to the conclusion that the distance traveled through an inclusion along a given ray has an effect on the distance to the next inclusion along that same ray. This dependence implies that ray tracing through a stochastic material is not a Markovian process, due to the presence of edge effects between inclusions. From the preceding discussion it is apparent that the edge-effect regions clearly have an impact on the PDF that describes the distribution of inter-inclusion distances. This effect stems from the fact that, within the edge-effect region, certain inclusion center point locations are disallowed because they would produce an unphysical overlap between adjacent inclusions or the material boundary. As a result, the volume of allowable inclusion center points is smaller within the edge-effect region than it would be in the absence of the interfering boundary regions. This truncation of the set of possible inclusion center points within the edge-effect region has the effect of changing the relative probabilities of observing certain intersection chord lengths within the next inclusion. To illustrate this effect, consider a fixed ray that intersects nearly tangentially with the surface of an inclusion, as shown in Fig. 4b). Simple trigonometry shows that it is impossible for the same ray to pass through the center point of the next inclusion pffiffiffi if the spacing between the inclusions is less than ð 31ÞDx=2, unless the inclusions partially overlap. This implies that for distances pffiffiffi x o ð 31ÞDx=2 (which fall within the edge-effect layer defined by x o Dx=2) it is not possible to realize an intersection chord of length Dx through the next inclusion. Thus, we find that sampled distances between

Extent of Edge-Effect Layer

Inclusion Inclusion

Δx/2 0

V(x)

x

Δx/2 0

x V(x)

Extent of Edge-Effect Layer

Fig. 4. Illustration of two-dimensional edge effect layer surrounding a disk-shaped inclusion, and the corresponding effect on V(x) for (a) a ray that passes through the center point of the inclusion, and (b) a ray intersecting the inclusion near to the inclusion surface.

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

inclusions that fall within the edge-effect region will influence the distance traveled through the next inclusion. This result further illustrates the non-Markovian nature of the ray tracing process through stochastic materials. However, for large materials with low inclusion volume fractions, it can be shown that the significance of inclusion edge effects is relatively minor and the transition between materials A and B in the material may be assumed to be Markovian without any substantial loss in accuracy. Finally, we note that the preceding discussion on the properties of V(x) in the edge-effect region may be easily generalized from 3-D space to 2-D by simply replacing the hemispherical shell of possible inclusion center points with a semi-circular set of points in the 2-D plane. However, further reduction to 1-D provides an interesting result. The complicated nature of the edge-effect layer is due to the intersections of the curved inclusion surfaces in two- and three-dimensional problems. In these higher dimensions we must concern ourselves with the overlap of inclusions away (above, below, etc.) from the ray of interest. In one-dimensional problems, however, there are no curved surfaces to contend with and the entire problem space is confined to the ray of interest. As a result, there can be no edge interference effects in onedimensional problems. This means that the original forms for V(x) and V0 (x) given in Eqs. (5) and (6) will hold for all one-dimensional problems. The lack of edge effects in one-dimensional problems also eliminates the dependencies between observed distances between and within inclusions, which were seen in higher dimensions. Thus we find that we can treat transitions between materials in a 1-D stochastic material as a Markovian process without any loss in generality.

2.4. Properties of the function q(v) Following our lengthy discussion on the properties of V(x), let us return our attention to the other unknown function present in the generalized PDF, q(v). Recall that q(v) was originally defined as the probability that no inclusion center points will be found in a given volume v. In contrast to the function V(x), which depends only on the geometry of the problem, q(v) is determined only by the statistical properties of the material. As a result, general conclusions regarding the properties of q(v) cannot be derived without knowledge of the particular stochastic material under consideration. In particular, we anticipate that q(v) will be heavily dependent on the inclusion volume fraction of the material, VB, which is defined as the ratio of the volume of material B (the inclusions) to the total volume of the material. For large inclusion volume fractions, the function q(v) will also depend on how the stochastic material was created. Stochastic materials may be formed either by an equilibrium process, which allows inclusions to interact with one another and move freely through the background material in order to find a suitable arrangement, or by random sequential addition (RSA), where inclusions are not permitted to move after their initial placement

583

within the model [21–23]. Clearly, stochastic materials that are formed by mixing a pre-determined volume of inclusions into the background material (such as concrete, HTGR fuel, and BORALs) will follow an equilibrium distribution. In addition, stochastic materials such as liquid/vapor mixtures (two-phase flow, clouds, etc.) with inclusions that are continuously in motion will follow an equilibrium distribution. However, as described earlier, these types of natural stochastic materials may exhibit correlation (clumping) between inclusion positions that may affect the equilibrium distribution [20]. By contrast, random sequential addition of inclusions can be used to model physical processes such as the adsorption of proteins onto the surface of an object [22] or the arrangement of randomly parked cars along a length of road [21]. For radiation transport applications we anticipate that most stochastic materials of interest will follow an equilibrium distribution of inclusions. However, for modeling purposes it is much easier to construct realizations of stochastic materials through RSA using rejection sampling techniques. For small inclusion volume fractions the equilibrium and RSA distributions of inclusions are nearly identical. However, as the inclusion volume fraction increases, RSA distributed inclusions will reach a jamming limit that is below the theoretical packing fraction for inclusions that are in equilibrium [21,23]. Beyond the general observations provided above, a complete discussion of the properties of q(v) for all possible types of stochastic materials and inclusion distributions is well beyond the scope of this paper. However, the statistical distribution of inclusions within stochastic mixtures has been a widely studied problem [14–19], with many analyses published on a wide variety of problems from diverse fields of study. Because of the general nature of the PDF for the distribution of distances between inclusions derived in Eq. (4), any knowledge regarding the behavior of q(V(x)) for a specific material can be directly incorporated into the formula for p(x). By necessity, such detailed analyses must be considered on a case by case basis. However, in many cases, reasonable assumptions regarding the general form of q(v) may be used to produce especially convenient solution forms for p(x). For the studies presented in this paper we will assume that q(v) has the general form qðvÞ ¼ elv ,

ð7Þ

where l is referred to as the rate parameter. Eq. (7) can be derived by assuming that the placement of inclusions in the material follows a Poisson point process. In a point process the probability of finding an inclusion center point within any small volume of space, Dv, is assumed to be constant, with the chance of success given by lDv. It follows that any arbitrarily large volume, v, can be broken up into v/Dv separate regions, each with probability lDv of containing an inclusion center point. Taking the limit as Dv-0, it can be shown that the probability of finding exactly k inclusions in the volume v will follow a Poisson distribution. Setting k= 0 in the resulting Poisson distribution gives the probability of finding exactly zero inclusion center points in volume v, and is equal to q(x), as given

584

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

in Eq. (7). Unfortunately, approximation by Poisson point process assumes that inclusions are singular points, and cannot account for the interference effects among inclusions with non-zero diameter. However, for low inclusion volume fractions the probability of overlap between Poisson distributed inclusions is small, and can often be neglected without affecting the shape of the distribution q(x) [12]. For convenience, we will use the Poisson point process assumption for the remaining derivations and examples in this paper. It is important to note that assuming an exponential shape for the function q(v), as shown in Eq. (7), is not the same as assuming that the distribution of chord lengths between inclusions will be exponentially distributed, which was made by Donovan and Dannon [5], Olsen et al. [12], and others. As described below, the presence of the derivatives in Eq. (4) will allow the probability distribution p(x) to take a non-exponential form within edge-effect regions of the material. Substituting our assumed form for q(x), Eq. (7), into the generalized PDF, Eq. (4), and evaluating yields pðxÞ ¼ lelVðxÞ

dVðxÞ : dx

ð8Þ

By inspection, we immediately recognize that, in the absence of edge effects, Eq. (8) reduces to a simple exponential distribution. This observation agrees with previously reported results, which have demonstrated that the distances between inclusions in an infinite material will be exponentially distributed [12–14]. In addition, many studies involving stochastic materials of a finite extent have started from an initial assumption that the inter-inclusion distances are exponentially distributed [4,5,10,11]. This convenient assumption provides excellent agreement with numerical experiments in regions away from the edge-effect layer of the material [7,12], but breaks down within one inclusion radius of the material boundary for two- and three-dimensional problems. This behavior is consistent with the derived form of p(x) shown in Eq. (8). Within the edge-effect region, the behavior of p(x) deviates from a pure exponential form due to the nonconstant factor dV(x)/dx, which accounts for the geometric restrictions imposed by the material boundary. Several studies have attempted to account for the presence of these edge effects by simply adjusting the rate parameter, l, in the exponential distribution [5,8,11]. The use of adjusted or empirically fitted rate parameters will result in improved agreement within the bulk material by correcting for the inability of the assumed Poisson point process distribution model to account for interference between inclusions in the material, as discussed in Section 2.4. However, adjustments to the rate parameter alone cannot account for the non-exponential behavior of p(x) within edge-effect regions. Therefore, the correct shape of the distribution p(x) can only be obtained through the inclusion of the geometry-dependent dV(x)/dx factor. Following this broad discussion on the nature of the generalized PDF p(x), and the related functions V(x) and q(v), which describe the distribution of inter-inclusion distances within a binary stochastic material, we will now present two in-depth case studies to verify the accuracy of

the derived form for the PDF and illustrate ideas from the preceding discussion. The first case study uses a one-dimensional binary stochastic material of finite length to demonstrate the applicability of p(x) to problems without edge-effect layers. Due to the simple nature of the one-dimensional problem we are able to provide additional statistical analysis of the material and even derive an analytical approximation for the rate parameter l, based on the inclusion volume fraction and the inclusion diameter. For the second case study, we consider a finite two-dimensional binary stochastic mixture with disk-shaped inclusions of fixed diameter. This study is intended to examine the accuracy of p(x) within edge-effect layers of the problem. In both cases, we are especially interested in assessing the accuracy and robustness of p(x) as a function of material properties, such as the inclusion volume fraction. To this end, the applicability of the assumed form of q(v) as the inclusion volume fraction increases is of particular interest.

3. Analysis of a one-dimensional binary stochastic material In this section we seek to develop a better understanding of the behavior of p(x) by considering the distribution of distances between adjacent inclusions in a 1-D binary stochastic media problem with finite length. For this analysis we will consider a one-dimensional slab of length L, which contains a number of inclusions of fixed length Dx distributed randomly throughout. The background material of the slab is identified as material A, and the material in the inclusions is identified as material B. Both the number and spatial location of the inclusions is assumed to be random, subject to the conditions that no inclusion may extend beyond the boundaries of the slab and no two inclusions may overlap with one another. An illustration of an example slab of this type is shown in Fig. 5. From Eqs. (4) and (6), we can immediately show that the probability distribution function for the distances between inclusions, x, in a 1-D slab is given by pðxÞ ¼ 

dqðxÞ : dx

ð9Þ

The simple form of Eq. (9) (particularly the lack of a dV(x)/dx factor) illustrates that finite 1-D stochastic materials do not suffer from the complex geometric edge effects present in 2-D and 3-D problems, as explained in Section 2.3. If we assume that the locations of the inclusion center points along the length of the slab follow a Poisson distribution,

Fig. 5. Illustration of a one-dimensional binary stochastic material with total length L, and inclusion length Dx.

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

then it is possible to substitute the exponential form for q(x), Eq. (7), into Eq. (9) to give an approximation for the probability distribution function pðxÞ ¼ lelx ,

ð10Þ

where l is a constant rate parameter, which gives the expected number of inclusion center points encountered per distance traveled. The resulting PDF for one-dimensional stochastic materials is known to be exact for ideal (infinite medium) cases [12,14], and has also been shown to produce accurate results in certain non-ideal (typically thick media) problems [5,10,11]. In 1-D infinite-medium cases the rate parameter, l, can be well approximated by the equation [12,14]

l1 ¼

VB , Dxð1VB Þ

ð11Þ

3.1. Derivation of the expected number of inclusions in a finite 1-D slab To quantify the accuracy of the infinite-medium rate parameter, Eq. (11), when applied to a finite system, we begin with the probability density function for the number of fixed-length inclusions in a slab of total length L, where the distances between inclusions are distributed according to Eq. (10) pðNÞ ¼ 8 1gð1, lDxðm1ÞÞ, > > > > < gðN, lDxðmNÞÞ gðN þ1, lDxððm1ÞNÞÞ  , ðN1Þ! N! > > > > : 0;

N ¼ 0, ( 0o N r m1, 0 o N r bmc,

m2Z

m= 2Z ,

Other:

ð12Þ In Eq. (12) N is the number of inclusions in the slab and m is the dimensionless slab length (defined as the ratio of the total slab length, L, to the inclusion length, Dx), and g(n, ly) is the lower incomplete gamma function, which is defined by Z ly gðn, lyÞ ¼ t n1 et dt: ð13Þ 0

A complete derivation of the probability distribution function given in Eq. (12) is provided in Ref. [24] and has been reproduced in Appendix A for the reader’s convenience. The probability distribution function given in Eq. (12) is a complicated, non-intuitive, function of both slab length, L, and inclusion width, Dx. In order to gain understanding about the function p(N) it is useful to consider the limiting behavior of the function for both large and small inclusion widths. It can be readily shown [24] that in the limit as the inclusion width becomes small, lim pðNÞ ¼

Dx-0

ðlLÞN elL ¼ ppois ðN; lLÞ, N!

where ppois(N; lL) is the Poisson distribution. Thus we see that as the inclusion width approaches zero, the number of inclusions in the slab approaches a Poisson distribution. This result should be expected due to the fact that the distribution of inter-inclusion distances, p(x), was derived by assuming that the inclusion center points follow a Poisson process. For the second limiting case, we will consider the behavior of p(N) as the inclusion width Dx becomes large (Dx-L). For values of Dx4L/2, Eq. (12) simplifies to a single-inclusion distribution function given by ( 1gð1, lDxðm1ÞÞ, N ¼ 0, pðNÞ ¼ ð15Þ gð1, lDxðm1ÞÞ, N ¼ 1: The expected value of the distribution given in Eq. (15) can be easily calculated, R  E½pðNÞ ¼ gð1, lDxðm1ÞÞ,

where Dx is the inclusion length and VB is the fraction of the slab length containing material B. However, for finitemedia problems, the rate parameter approximation given in Eq. (11) tends to under-predict the expected number of inclusions occurring in the slab [24].

ð14Þ

585

ð16Þ

where R is defined as the probability of finding an inclusion in the slab. Substituting Eq. (16) into Eq. (15) allows the single-inclusion probability distribution function to be written in terms of R ( 1R, N ¼ 0, pðNÞ ¼ ð17Þ R, N ¼ 1, which is the Bernoulli distribution with the variable R as the probability of success. Thus we see that, in cases where the inclusion width is greater than half the total slab length, the observation of an inclusion in the slab reduces to a Bernoulli experiment; for each trial there is some fixed probability of finding an inclusion in the slab. Now that we have developed an analytical formula for the probability of finding exactly N inclusions in a slab of total length L, we can calculate the expected number of inclusions in the slab /NS ¼

1 X

n pðnÞ:

ð18Þ

n¼0

Substituting Eq. (12) into Eq. (18) and evaluating yields a final expression for the expected number of inclusions in a slab of length L /NS ¼

bm c1 X n¼0

gðn þ 1, lDxððm1ÞnÞÞ n!

,

ð19Þ

as shown in Ref. [24] and Appendix B. A simple Monte Carlo simulation was used to confirm the accuracy of Eq. (19). In the numerical experiment, multiple independent realizations of a 1-D binary stochastic material were created and used to estimate the ensemble average number of inclusions as a function of m and l. For each realization, inclusions of width Dx were placed at random intervals within a slab of length L= 1. The spacing between adjacent inclusions was determined by sampling from an exponential distribution (Eq. (10)) with a given (fixed) value of l. When no more inclusions would fit within the slab, the total number of inclusions for the realization was recorded and the process began again for the next realization. After many realizations the population mean and variance was used to estimate the expected (ensemble average) number of inclusions for the

586

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

Table 1 Comparison of analytical and statistical estimates for the expected number of inclusions in binary stochastic materials of length L =1. Statistical estimates were generated with 50,000 independent realizations. Uncertainty values shown represent 95% confidence interval on the result.

l

m

hN i

^ N

1 10 100 100 100 100 100

10 10 1.01 2.25 5.5 10 100

0.8223 4.6250 0.6285 1.9998 4.9480 8.6664 49.6250

0.8209 7 0.0082 4.62937 0.0192 0.6267 7 0.0070 1.99987 0.0126 4.94827 0.0199 8.6660 7 0.0262 49.64267 0.0630

slab. Table 1 shows the results of 50,000 realization numerical experiments performed with different values of m and l. In all cases the calculated ensemble average was found to agree to within statistical uncertainty with the expected value given by Eq. (19).

3.2. Approximating the rate parameter l The results in Section 3.1 provide an analytical equation for the expected number of inclusions as a function of the geometric properties, L and m= L/Dx, and the rate parameter l. Ideally, we would like to be able to directly solve for the rate parameter as a function of the material properties. Unfortunately, Eq. (19) contains two unknown quantities, /NS and l, so it will be necessary to find a second relationship between these quantities before a solution can be attempted. The required condition can be obtained from a volume balance equation, which relates the expected number of inclusions in the slab to the material properties L, Dx, and VB, defined as the total volume fraction occupied by material B in the slab. In this case, VB is treated as a known material property, which may be obtained directly from manufacturing data (for man-made systems such as concrete, HTGR fuel elements, or BORALs), or estimated through auxiliary measurements of the bulk material density and average inclusion size (for naturally occurring stochastic mixtures). In either case, for a fixed inclusion width of Dx, the expected number of inclusions is given by [12,24] VB L ¼ VB m, N¼ Dx

ð20Þ

where again, VB is the volume fraction of material B, as determined from physical measurements or manufacturing data. Eq. (20) provides the necessary relationship required to solve for the unknown rate parameter, l. Assuming that the distribution of distances between inclusions is well represented by an exponential distribution, then the unknown rate parameter, l, must satisfy the relationship /NS ¼ N, where /NS is the expected number of inclusions (Eq. (19)) obtained by calculating the expectation value of p(N), and N is the expected number of inclusions based on physical (mass/volume balance) conservation properties of the problem (Eq. (20)).

3.2.1. Single-inclusion case For cases where Dx4L/2 (i.e., mo2) the expected value of p(N) has a particularly simple form, as shown in Eq. (15), /NS ¼ gð1, lðLDxÞÞ ¼ 1elðLDxÞ

ðfor Dx 4 L=2Þ:

ð21Þ

Because of the simple form it is possible to analytically solve Eq. (21) for l in terms of the expected number of inclusions /NS. Applying the condition /NS ¼ N, and using Eq. (20), yields a simple expression for the rate parameter,

l¼

ln½1VB mDx ðm1Þ

ðfor m r 2Þ:

ð22Þ

Again, we remind the reader that the expression given in Eq. (22) is exact if the distribution of distances between inclusions (or in this case the distance between the inclusion and the edges of the slab) follow an exponential distribution. However, the denominator of Eq. (22) provides some indication about the limitations of this original assumption. As the inclusion width, Dx, approaches the total length of the slab, L, (i.e., m-1) the denominator goes to zero, which causes the rate parameter to approach N. This result is in keeping with the physical realities of the problem, which dictate that as Dx-L, the position of the inclusion within the slab becomes deterministic rather than stochastic. In this case, the use of an exponential distribution to describe the position of the inclusion within the slab is no longer appropriate. Instead, the Bernoulli distribution should be used to describe the probability that an inclusion will (or will not) be found in a particular observation of the stochastic material. 3.2.2. Multiple-inclusion case Inspection of Eq. (19) reveals that the equation cannot be easily inverted to solve for l in terms of /NS for cases where mZ2. This eliminates the possibility of finding a simple, analytic, solution for l that will cover all cases of interest. Instead, we return our attention to the reasonability of using the original, infinite-medium, estimate for the rate parameter given in Eq. (11),

l1 ¼

VB : ð1VB ÞDx

ð11Þ

To determine this accuracy of the infinite-medium rate parameter for finite 1-D problems, we will consider   the residual difference between the true mean Nl1 calculated by Eq. (19) using the infinite-medium rate parameter (Eq. (11)), and the predicted mean N as defined in Eq. (20)

e ¼ /Nl1 SN ¼

bm c1 X n¼0

gðn þ 1,ðVB =ð1VB ÞÞððm1ÞnÞÞ n!

! VB m:

ð23Þ Numerical results for the accuracy of the infinitemedium rate parameter as a function of volume fraction (VB) and dimensionless slab length (m) are provided in the next section. Previous studies by Murata et al. [8], and Ji and Martin [11] have shown that the accuracy of the infinite-medium

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

rate parameter can be further improved by using an effective inclusion volume fraction, VnB, in place of the traditional VB in Eq. (11). In Ref. [11] Murata et al., proposes using an effective volume fraction given by VB ¼ VB

Vregion , Vregion DVedge

ð24Þ

where Vregion is the volume of the entire stochastic material region and DVedge is the volume of the edge boundary layer extending one inclusion radius into the material from the outer surfaces. In a finite one-dimensional system Eq. (24) simplifies to VB ¼ VB

m : m1

ð25Þ

Inspection of Eq. (25) reveals that the ratio of Murata’s effective volume fraction, VnB, to the true volume fraction, VB, is a constant value, which depends only on the dimensionless slab length, m. The volume correction shown in Eq. (25) can be applied to any stochastic material, so long as the resulting effective inclusion volume fraction, VnB, does not exceed 1. This implies that the effective volume fraction correction can only be used for stochastic materials with true inclusion volume fractions VB,max o(m 1)/m. An improved, empirically based formula for the  effective inclusion volume fraction, V~ B , was developed in Ref. [24], and is given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  V~ B ¼ 1m þ ðm1Þ2 þ 2VB m, ð26Þ where again, VB is the true inclusion volume fraction and m is the dimensionless slab length. Unlike Murata’s formula, which predicts that the effective volume fraction is linearly proportional to the true volume fraction, Eq. (26) provides for a non-linear relationship between the effective and true inclusion volume fractions.

587

In particular, the relationship given in Eq. (26) predicts a  decrease in the ratio V~ B =VB as the true inclusion volume fraction, VB, increases [24]. However, in the limit as the true inclusion volume fraction approaches zero (i.e., highly dilute systems) the effective volume fraction predicted by Eq. (26) approaches the value predicted by Murata’s formula, as given in Eq. (25). As before, we see that Eq. (26) cannot be applied to stochastic materials that would result in an effective inclusion volume fraction, V*B, greater than 1. This condition prevents use of the improved effective inclusion volume fraction for stochastic materials with a true inclusion volume fraction, VB,max 41 1/(2m), as shown in Ref. [24].

3.3. Numerical results The accuracy of the rate parameters calculated from Eqs. (22) and (11) for finite stochastic media problems was determined through a series of numerical experiments that were designed to cover a wide range of material properties (specifically, dimensionless slab lengths, m, and inclusion loading densities, VB). In each experiment, either Eq. (11) or (22) was used to estimate the rate parameter for the material, depending on the particular values for VB and m, and assuming L=1 in all cases. For each set of experiments, the same type of Monte Carlo simulation used to verify Eq. (19) was again used to estimate the expected (ensemble average) number of ^ l , for the calculated rate inclusions in the slab, N parameter. Each experiment used 50,000 independent ^ l . The realizations to produce the ensemble average, N ^ l , was then compared calculated number of inclusions, N against the desired number of inclusions, N, as calculated by Eq. (20). The relative accuracy of the rate parameter

Table 2 Comparison of analytical and statistical estimates for the expected number of inclusions in binary stochastic materials of length L= 1. Statistical estimates were generated with 50,000 independent realizations, where distances between inclusions were sampled from an exponential distribution using a calculated rate parameter. For materials with mr 2 the rate parameter was calculated from the exact formula for the rate parameter in a single-inclusion slab. For materials with m42, the rate parameter was calculated with the infinite-medium approximation formula using the true (VB), effective (V*B), and modified effective (VB ) inclusion volume fractions for the material. Uncertainty values shown represent 95% confidence the result. VB

m

N

^l N

mr2 0.8 0.6 0.4 VB

1.001 1.50 2.00

0.8008 0.900 0.800

0.80247 0.0080 0.90067 0.0084 0.7986 7 0.0080

m

N

^l N 1

^  N l

^ ~ N l

2.0001 2.0001 2.0001 3 4 10 10 20

0.2 0.6 1.0 1.5 2.0 8.0 9.0 19

0.10657 0.0030 0.3510 7 0.0027 0.6330 7 0.0072 1.1310 7 0.0084 1.6201 7 0.0114 7.5240 7 0.0246 8.54157 0.0262 18.53087 0.0386

0.2202 7 0.0021 0.7770 7 0.0039 1.5930 7 0.0056 1.79587 0.0060 2.23387 0.0067 8.39967 0.0130 9.0000 7 0.0134 N/A

0.2026 7 0.0040 0.6088 70.0035 0.9356 70.0086 1.517370.0110 1.9959 70.0126 7.9848 70.0252 8.9958 70.0268 18.9999 70.0390

m 42 0.1 0.3 0.5 0.5 0.5 0.8 0.9 0.95

588

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

calculated from Eq. (11) or (22) can be determined by ^ l approximates N. measuring how closely N Table 2 shows the results for eleven stochastic mixtures, which cover a wide range of inclusion volume fractions and dimensionless slab lengths. For the three materials with m r2 (i.e., single-inclusion slabs), the rate parameter l was calculated from the exact single-inclusion formula given by Eq. (22). As anticipated, the results for these materials (Table 2) demonstrate that the calculated ensemble average number of inclusions matches the predicted expected number of inclusions to within statistical uncertainty for all three cases. For the eight stochastic materials with m42, three separate rate parameters, lN, ln, and l~ * were calculated from the infinite-medium approximation given in Eq. (11). The rate parameter lN was calculated directly from Eq. (11) using the true volume fraction, VB, of the material, while the corrected rate parameters, ln, and l~ * were calculated from  Eq. (11) using the effective volume fractions V*B and V~ B , respectively, as given in Eqs. (25) and (26). The results for the slabs with m 42 (Table 2) show that the use of the uncorrected infinite-medium formula for the rate parameter, Eq. (11), produces poor results for finite-length slabs. In all cases considered, the ensemble average number of observed inclusions was lower than the expected number of inclusions by a statistically significant amount. The use of the infinite-medium rate parameter was especially poor for small dimensionless slab lengths (m= 2.0001), where the number of inclusions in the slab was under-predicted by nearly half (36–47%). As the dimensionless slab length, m, increases, the accuracy of the infinite-medium rate parameter improves, with the under-prediction in the number of inclusions dropping to less than 5% for slabs with m Z10, even for large inclusion volume fractions. When the infinite-medium formula for rate parameter is used in conjunction with the effective inclusion volume fraction proposed by Murata et al., the agreement between the observed (ensemble average) number of inclusions and the predicted number of inclusions improves for most of the cases looked at. However, even with Murata’s effective inclusion volume fraction the differences between the expected and observed number of inclusions remain outside of the statistical uncertainty band, indicating that a systematic bias remains. In contrast to the uncorrected infinite-medium rate parameter, Murata’s volume correction tends to over-predict the number of inclusions in the slab. In the case where m=2.0001 and VB =0.01 (representing a dilute stochastic material in a very short slab, which can hold at most two inclusions side by side), the use of Murata’s effective volume fraction over-predicts the number of inclusions in the slab by 10%. This is a significant improvement over the 47% under-prediction obtained with the uncorrected infinite-medium rate parameter. However, for longer slabs with higher inclusion volume fractions, the impact of the effective volume fraction is reduced relative to the uncorrected infinite-medium rate parameter, providing only a 1–7% improvement in accuracy. In addition, the effective volume fraction proposed by Murata et al., behaves unpredictably for cases near the VB =(m–1)/m volume

fraction limit for the formula. In one case (m=2.0001, VB =0.5) the use of the effective volume fraction over-predicts the number of inclusions by nearly 60%, producing a result that is less accurate than produced by the uncorrected infinite-medium rate parameter approximation. However, as the volume fraction and dimensionless slab length increase, the use of the effective volume fraction improves, yielding a nearly exact result for the case where m=10 and VB =0.9. Finally, we wish to point out that the volume fraction for the last stochastic material (m=20, VB =0.95) listed in Table 2 is larger than the maximum volume fraction that can be used in Eq. (25) to calculate VnB. Therefore, no results are shown for this case. When the infinite-medium formula for rate parameter is used in conjunction with the modified effective  inclusion volume fraction, V~ B , given in Eq. (26), a significant improvement between observed (ensemble average) and predicted number of inclusions is seen. With the modified effective volume fraction the relative agreement between the observed and predicted number of inclusions in the slab improves to better than 98.5% in most cases looked at. In addition, the difference between the predicted and observed number of inclusions for 6 of the 8 cases considered falls within the reported 95% confidence interval. However, the case with m =2.0001 and VB =0.5 shows a difference of nearly 6.5% between predicted and calculated results. For this particular case the volume fraction VB = 0.5 is considerably larger than the recommended maximum volume fraction, VB,recommended o1

3 1 þ , 2m 2m3

ð27Þ

as described in Ref. [24]. Table 2 also shows that the agreement between the calculated and predicted values appears to improve as the values of m increase, with agreement for the m =20 case reaching 99.9995%. In order to further quantify the error of the approx imation for l~ as a function of the material properties VB and m, an additional set of numerical experiments were conducted. These experiments were designed to determine the largest (uncorrected) inclusion volume fraction that would result in the absolute value of the approximation error e (Eq. (23)) being less than or equal to a constant value, as a function of the scale parameter m. The resulting iso-error plots from these experiments are shown in Fig. 6, for 6 different absolute errors. Each curve was generated by first calculating the absolute error of the approximation error (computed from Eq. (23)) at 200 equally spaced values of VB between 0 and 1, for a fixed value of m. The 100 resulting values for e were then searched in ascending order by volume fraction to find the first value that exceeded the reference error for the curve. This process was then repeated for 490 equally spaced values of m, ranging from 2 to 100, in order to produce the curves shown in Fig. 6. The results shown in Fig. 6 indicate a well defined relationship between the accuracy of the approximation and the material properties VB and m. As expected, the accuracy tends to decrease as m becomes small and VB becomes large. For values of m less than 10, the accuracy

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

1− 1

1 2m

1=0.1 1=0.01 1=1×10-3 1=1×10-4 1=1×10-5 1=1×10-6

0.8

1−

VB

0.6

589

1 3 + 2m 2m 3

0.4

0.2

20

40

60

80

100

m Fig. 6. Iso-error plots showing the largest allowable inclusion volume fraction for a given level of absolute accuracy (e), as a function of the dimensionless scale parameter m.

of the method begins to fall off sharply for inclusion volume fractions greater than roughly 50%. However, for larger values of m (m 425), the approximation shows excellent agreement with theoretical results, with errors on the order of 1  10-6 (absolute) for inclusion volume fractions up to 80%.

4. Analysis of a two-dimensional binary stochastic material The progression from one-dimensional to two-dimensional problems increases the difficulty of the analysis of stochastic materials substantially. As discussed in Section 2.3, stochastic mixtures in two and three dimensions exhibit edge-effect layers, which are not present in onedimensional problems. The presence of these edge-effect layers presents two main difficulties. Firstly, the volume (or area, in two dimensions) of allowable center points as a function of distance traveled through the material, V(x), is non-trivial to calculate and depends heavily on the specific geometry of the problem, as described in Section 2.3. Secondly, the presence of a variable dV(x)/dx factor in the PDF p(x) means that inter-inclusion distances within the material will not follow an exponential distribution. Given that it was the especially convenient properties of the exponential distribution that enabled a thorough analysis of the one-dimensional problem, it follows that a similar set of comprehensive theoretical results will not be possible for higher-dimensional problems. In particular, we find that the presence of edge-effect layers within the material prevents any straightforward derivation for a PDF describing the number of inclusions that intersect an arbitrary ray intersecting the material volume.

Recognizing that it will not be possible to produce the same level of comprehensive analysis or the number of theoretical results for a two-dimensional stochastic material, we will restrict our attention to the behavior of p(x) within edge-effect regions of the material. During the discussion of edge-effect layers in Section 2.3, it was stated that the distribution of inter-inclusion distances, p(x), may not be the same for all reference points within a stochastic material. This observation is confirmed by Murata et al., who notes that three separate types of probability density functions (referred to as nearest neighbor distributions) are required for analysis of a finite two- or three-dimensional problem [7]. These PDFs are required to accurately describe the different distributions of inter-inclusion distances relative to a reference point that is located: (1) on the outer surface of the material, (2) on the surface of an inclusion within the material, and (3) within the background material itself, away from any material or inclusion boundaries. For our case study, we will consider the behavior of the probability distribution function, p(x), only for the first situation described above, and for a particularly simple problem geometry. However, our method of analysis should be directly applicable to other types (or shapes) of edge-effect layers as well as three-dimensional problems. For our two-dimensional test case we will consider a rectangular box filled with a stochastic mixture of homogeneous, disk-shaped, inclusions of material B uniformly distributed within a uniform background of material A. All inclusions are assumed to have a fixed diameter, Dx, and the length and width of the bounding box is denoted by the variables L and W, respectively. An illustration of the problem geometry is provided in Fig. 1. We wish to consider the distribution of distances to the

590

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

nearest inclusion as measured from the left-hand surface of the material region, along a ray perpendicular to the surface. This case is of particular interest because it includes the presence of an edge effect layer of thickness Dx/2, as illustrated in Fig. (3a). To simplify the analysis, we will also assume that the point of intersection between the reference ray and the left surface of the region is located more than Dx/2 from the top and bottom of the material. This precludes the presence of edge effects due to the top and bottom surfaces along the ray of interest.

0 x

V(x)

4.1. Derivation of V(x) and p(x) We begin our analysis by deriving the expression for V(x), where x denotes the distance along the reference ray from its intersection point with the left-hand surface. As discussed previously, V(x) will have a different form in the edge effect layer (0rx r Dx/2) than in the bulk material. Thus, it will be necessary to compute V(x) independently in the two regimes and then combine the two solutions in a piecewise definition for the final V(x). We begin by considering the behavior of V(x) in the bulk material. Outside of the edge effect layer (x4 Dx/2), V(x) will have a bullet shape, with a rectangular body of length x  Dx/2 and width Dx, capped with a semi-circular head of radius Dx/2. Therefore, we can immediately show that the area of this subregion is given by   Dx pDx2 þ VðxÞ ¼ Dx x 2 8

Dx

for

2

o x oL

Dx 2

:

ð28Þ

The derivative of Eq. (28) can be readily computed to yield VuðxÞ ¼ Dx

for

Dx 2

o x o L

Dx 2

,

ð29Þ

which matches the expected result (for the bulk material region), previously reported in Eq. (6). Within the edge effect layer, xr Dx/2, the determination of V(x) becomes more difficult due to the need to calculate the area within the intersection between a semicircle and a line, as shown in Fig. 7. Through simple trigonometry, we can show that this area can be expressed as the definite integral

VðxÞ ¼ 2

Z

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 Dx x~ 2 dx~ 2 Dx=2x Dx=2

for

xr

Dx 2

,

ð30Þ

Surface Boundary Layer

Δx/2 Fig. 7. Illustration of the intersection between the set of allowable inclusion center points and the surface boundary layer, in two dimensions. The shaded region defines V(x) for a point x, located within the edge-effect region, as shown.

which can be evaluated to yield " pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pDx2 ðDx2xÞ xðDxxÞ VðxÞ ¼  8 2 " ## Dx2 ðDx2xÞ arctan pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ 4 2 xðDxxÞ

for

Evaluation of V(x) for x= 0 and x= Dx/2   Dx pDx2 ¼ , Vð0Þ ¼ 0, V 2 8

xr

Dx 2

:

ð31Þ

ð32Þ

confirms the expected behavior of V(x) , with the function ranging from 0 to pDx2/8 within the edge effect layer. In particular, the value of V(Dx/2) matches the value of Eq. (28), providing the anticipated continuity of V(x) at the edge effect boundary. Finally, with an analytical form of V(x) within the boundary layer, we can evaluate the derivative of the function pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dx VuðxÞ ¼ 2 xðDxxÞ for x r : ð33Þ 2 Having derived the V(x) functions for the edge effect and bulk material layers independently, we can now combine Eqs. (28) and (31) to produce a piecewise definition for V(x) for all values of xoL  Dx/2 (up to the edge effect boundary for the right-hand side of the box),

" " ## pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 > pDx2 ðDx2xÞ xðDxxÞ Dx2 ðDx2xÞ > > þ ArcTan pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , > < 8  2 4 2 xðDxxÞ VðxÞ ¼   2 > pDx Dx > > > , þ Dx x : 2 8

xr

Dx 2

Dx 2

,

ox o L

Dx 2

ð34Þ :

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

Similarly, we can create a piecewise definition for the derivative of V(x) over the same domain 8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dx > > 2 xðDxxÞ, x r , dVðxÞ < 2 ð35Þ ¼ Dx Dx > dx > : Dx, o x oL : 2 2 Eqs. (34) and (35) provide the geometry-dependent information needed to determine p(x), the probability density function describing the distribution of distances to the first inclusion as measured from the left-hand surface of the material. Substituting Eqs. (34) and (35) for this geometric arrangement into the generalized formula for p(x), as derived in Section 2, gives the desired probability distribution function

disposal, we will rely on numerical experiments to determine reference values for the rate parameter l. One such technique for experimentally determining l was used by Donovan and Danon in their study of the Limited Chord Length Sampling (LCLS) method [5]. In order to calculate l, they first performed a Monte Carlo experiment to determine the probability, p(N =0), that a ray of fixed length could pass through a given stochastic material without intersecting any inclusions. By assuming an exponential distribution (with constant rate parameter l) for the inter-inclusion distances and integrating this distribution over all distances greater than L an expression for l as a function of p(N =0) can be easily derived. However, the required assumption of a pure exponential

" " ##! 8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > pDx2 ðDx2xÞ xðDxxÞ Dx2 ðDx2xÞ > > þ  ArcTan pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l > 2 8 4 2 xðDxxÞ > < 2l xðDxxÞe , pðxÞ ¼     2 > p D x D x > > > þ Dx x l > > 2 8 > , : lDxe

where l is a constant rate parameter, as described in Section 2.4. A plot of the relative shape of the derived p(x) is shown in Fig. 8 for a binary stochastic material with parameters l = 1 and Dx= 1. The general shape of the probability density function shows qualitative agreement with results of numerical experiments conducted on finite-media problems by Murata et al. [7] and Ji and Martin [10]. Fig. 8 also demonstrates that the probability distribution function given in Eq. (36) differs significantly from an exponential distribution within the edge effect layer. Therefore, any assumption that the distances to the first inclusion are exponentially distributed will result in an over-prediction of the number of inclusions (or probability of finding an inclusion) within the edge effect layer. The magnitude of this over-prediction will depend on the properties (l and Dx) of the material. For a binary stochastic material with l =1 and Dx = 1, the exponential distribution will over-predict the probability of finding an inclusion within the edge effect layer (Dx/2 of the material boundary) by 21.25%.

591

xr

Dx 2

Dx 2

o x o L

ð36Þ

Dx 2

,

form for p(x) will not correctly account for the edge-effect layer in the material, which makes this technique unsuitable for our application. To develop an alternate approach for experimentally determining l, we return to our original definition for l, given in Section 2.4 qðvÞ ¼ elv ,

ð7Þ

we see that l is related to the probability that no inclusion center points will be found within a volume v of the material. This observation suggests a slight modification to the technique developed by Donovan and Danon. Instead of calculating the probability that no inclusion intersections will occur along a ray of fixed length, we can instead use a similar Monte Carlo simulation to determine the probability that no inclusion center points fall within a region of fixed volume (or area in the two-dimensional case). The resulting measured probability is, by definition, equivalent to q(v). Therefore, we can easily invert Eq. (7) to express the rate parameter as a function of the measured value of q(v)

4.2. Approximating the rate parameter l

l¼ For the one-dimensional case it was possible to derive an exact analytical relationship between the expected number of inclusions in the slab and the rate parameter, l. This relationship served as a convenient metric for quantifying the accuracy of different approximations for l. Unfortunately, in higher dimensions we do not know of any method for deriving a similar test function. Thus, we will be somewhat limited in our ability to quantitatively assess the accuracy of existing approximations for l, especially in the limit of large inclusion volume fractions (near the jamming or packing fraction for the material). Without sufficient analytical tools at our

ln½qðvÞ : v

ð37Þ

Ideally, for any given stochastic material, the rate parameter will be a constant value, regardless of the size or shape of the test volume v used in Eq. (37). As noted in Section 2.4, the properties of q(v) in stochastic materials is a well-studied problem, related to the development of Nearest Neighbor Distribution (NND) functions. In particular, Torquato and Lu have derived analytical solutions for q(v) in one-, two-, and three-dimensional problems under certain conditions [14–18]. For infinite-medium problems with uniformly distributed inclusions of fixed diameter Dx, these analy-

592

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

4.3. Numerical results

e–xΔx

1

0.8

p (x)

0.6

0.4

0.2

1

2

3

4

x

Edge-Effect Layer

Fig. 8. Probability distribution of distances x between a plane and a randomly placed disk-shaped inclusion, along a ray perpendicular to the plane. Solid line illustrates the derived generalized probability distribution function, p(x), for a stochastic material with l = 1 and Dx = 1. The dashed line shows an exponential distribution with the same values for l and Dx for comparison.

tical results can be used to show that the rate parameter l can be well approximated by [14]

l1 ¼

8 1 VB > > , > > Dx ð1VB Þ > > > < 4 VB

, pDx2 ð1VB Þ > > > > > VB > 3 > , : 2Dx3 ð1VB Þ

1-D, 2-D ,

ð38Þ

3-D:

The approximations given in Eq. (38) are equivalent to the expressions for l used by Ji and Martin [10,11], and Olson et al. [12] for studies on stochastic materials in two- and three-dimensions. As noted in the analysis of a one-dimensional slab, Murata et al. has proposed modifying the inclusion volume fraction, VB, to account for edge effects in finite regions and achieve better agreement with experimental results [8]. The effective inclusion volume fraction, V*B, may be calculated by applying an analytical correction factor to VB, giving VB ¼ VB

Vregion , Vregion DVedge

ð24Þ

where, again, Vregion is the volume of the entire stochastic material region and DVedge is the volume of the edge boundary layer extending one inclusion radius into the material from the outer surface. A subsequent study by Ji and Martin indicates that the effectiveness of this correction in two and three-dimensions is highly dependent on the geometry of the material and cannot be expected to provide more than a first-order improvement for general cases [11]. For our two-dimensional study we will find it worthwhile to consider the infinite-medium approximation for l given in Eq. (38) both with and without the volume correction factor proposed by Murata.

The accuracy of the derived probability density function p(x), given in Eq. (36), and the corresponding infinitemedium approximation for the rate parameter, lN, given in Eq. (38), were verified through a series of Monte Carlobased numerical experiments. These experiments follow the general procedure used by Donovan and Danon [5]. All simulations use a common material geometry, defined as a two-dimensional rectangular region, as illustrated in Fig. 1, with length L= 50 and width W=20. A stochastic mixture is created in this region by randomly placing a fixed number, N, of disk-shaped inclusions, subject to the condition that no inclusion may overlap with the material boundary or another inclusion within the material. Each random set of inclusion positions within the rectangular region then results in a unique realization of a twodimensional stochastic material. Placement of the inclusions is achieved by the process of random sequential addition. For each inclusion added to the material, a candidate center point location is sampled from a uniform distribution defined over the rectangular region. If the sampled center point is located within a distance Dx/2 of the material boundary or the surface of any adjacent (previously emplaced) inclusions, the candidate center point is rejected and resampled. Once an allowable center point has been sampled, an inclusion is positioned at that location and the process continues until all N inclusions have been placed. For each stochastic material realization, the distance from the edge of the material region to the surface of the nearest inclusion can be computed. All distances, x, are measured from the left-hand bounding surface of the material region (defined as x =0), along a ray oriented parallel to the length of the region and located a distance W/2= 10 from the bottom surface of the region. The computed distance, x, for each material realization is itself a random variable, which should theoretically follow the probability distribution p(x), as described in Section 4.1. Agreement between the theoretical and experimentally measured distributions for p(x) was verified for stochastic materials covering a range of inclusion volume fractions (VB). In each experiment, the distance from the boundary to the surface of the first inclusion was measured for ten million independent material realizations. A fixed inclusion diameter of Dx=1.5 was used for all simulations. The inclusion volume fraction of the stochastic material was Table 3 Comparison of analytical and statistical estimates for the rate parameter l for stochastic materials with a range of inclusion volume fractions, VB. Reference results lref are reported for a distance X = 1.5 from the material boundary. For each material type, reference results were generated from one million independent realizations of the stochastic material. N

VB

V*B

lref

l1

l

6 28 57 85 113

0.010603 0.049480 0.100727 0.150207 0.199687

0.011817 0.055146 0.112262 0.167409 0.222555

0.00678 0.03293 0.07125 0.11337 0.16162

0.00606 0.02946 0.06339 0.10002 0.14120

0.00677 0.03303 0.07156 0.11378 0.16199

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

varied between experiments by adjusting N, the number of inclusions placed in the material. A total of five experiments were run, with N=6, 28, 57, 85, 113. Given the inclusion diameter Dx and the overall dimensions (L and W) of the material region, these values of N correspond to inclusion volume fractions ranging between of VB =0.01060 and 0.19969, as shown in Table 3. Before a fair comparison between the theoretical and measured distributions of p(x), can be performed, we must determine the true rate parameter l for each of the different stochastic materials under consideration. As discussed in Section 4.2, the rate parameter can be obtained from measurements of q(V(X)), which corresponds to the probability of finding zero inclusion center points within a volume V(X). Because the volume V(X) corresponds to the set of allowable inclusion center points that are located less than a distance X from the reference point, it follows that q(V(X)) is equal to the probability that the distance to the first inclusion, x, is greater than X, as shown in Eq. (2). We can therefore estimate q(V(X)) from a set of M random material realizations as ~ qðVðXÞÞ ¼

M 1 X HðXxi Þ, Mi¼1

ð39Þ

where xi are the measured distances for each material realization and H(X xi) is the Heaviside step function. ~ Given the estimated qðVðXÞÞ determined from the numerical experiment, we can then calculate the rate parameter l using Eqs. (37) and (34) for each stochastic material. Given that Eq. (39) allows us to estimate q(V(X)) for many values of X, using only a single set of material realizations, we see that we can actually compute the rate parameter as a function of distance X. A plot of the resulting rate

593

parameters, lN(X), is shown in Fig. 9 for stochastic materials with the five different volume fractions listed above. As expected, Fig. 9 illustrates that the true rate parameter l is relatively constant over all distances X. For larger volume fractions, the rate parameter tends to increase slightly for larger distances (410 Dx). It is worth noting that this behavior has been observed in previous studies, and is discussed briefly in the paper by Olsen et al. [12]. We believe that this result indicates a limitation in the applicability of the Poisson point process approximation used in Section 2.4 to derive q(V(X)) as X and VB become large. We conjecture that this effect stems from the inability of the Poisson point process model to account for overlap between inclusions. In particular, this model cannot account for the effects of inclusions whose center point do not lie in V(X), but whose area (and associated edge effect layer) extends partially into the V(X) region. The presence of these partially overlapping inclusions will provide an additional limitation on allowable center point locations within the material, and, therefore, the true shape of V(X). Based on an analysis of the results shown in Fig. 9, a reference value for the true rate parameter was calculated for each type of stochastic material. In all cases the reference rate parameter was taken to be the calculated rate parameter at X= Dx= 1.5. This value was judged to be representative of a wide range of distances for each of the materials considered. Table 3 lists the measured (reference) rate parameter, lref, along with a comparison to the theoretical infinite-medium rate parameter, lN, computed from Eq. (38). In addition, the effective inclusion volume fraction, VnB for each material was computed using the formula provided by Murata, given in Eq. (24). The theoretical infinite-medium rate parameter was then recalculated with the effective inclusion

0.2

VB=19.97%

0.18 0.16 0.14

VB=15.02%

λ (X)

0.12 0.1

VB=10.07%

0.08 0.06

VB=4.95%

0.04 0.02 0 0.01

0.1

1 X

10

VB=1.06%

100

Fig. 9. Measured rate parameter, l, as a function of distance, X, for five stochastic materials with inclusion volume fractions as indicated.

594

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

0.25

0.2

p (x)

0.15

0.1

0.05

0 0.01

01

1 x

10

100

Fig. 10. Comparison between observed (points) and predicted (lines) distribution of distances between the material surface and the nearest inclusion boundary for five stochastic materials with inclusion volume fractions as indicated. Error bars on the data points show the 95% confidence band.

volume fraction and the resulting ln values are shown in Table 3. Inspection of Table 3 indicates good agreement between the measured rate parameter and the theoretically predicted values. For the range of material properties considered in this study it appears that the use of the effective inclusion volume fraction VnB produces a more accurate result. Now that we have methods for analytically and experimentally calculating values for l, it is possible to directly compare the measured distribution of distances, p(x), against the theoretically predicted distribution given in Eq. (36). Fig. 10 shows the comparison between measured and theoretical results for all 5 stochastic materials under consideration. In all cases the corrected infinite-medium rate parameter, ln, for the particular material was used to generate the theoretically predicted PDF (Eq. (36)). The results indicate excellent agreement between measurement and prediction for all distances 0 ox o48.25Dx. Inspection of Fig. 10 shows that the theoretically derived p(x) correctly predicts the distribution of distances within the edge effect layer, which extends from 0 ox r Dx/2. These results provide a compelling experimental verification for the accuracy of the generalized probability density function for the distribution of inter-inclusion distances provided in Eq. (4).

5. Conclusions A detailed analysis of the distribution of distances (chord lengths) between adjacent inclusions in a finite, isotropic, binary stochastic material with fixed diameter inclusions has been performed. This analysis is intended

to account for both the random distribution of inclusions within the stochastic material, as well as the geometric constraints imposed by the outer boundaries of the finite material. As a result of this analysis, a generalized probability density function (PDF) for describing the distribution of inter-inclusion distances in a binary stochastic material has been developed and tested. The new probability density function explicitly accounts for edge effects present in finite two- and three-dimensional stochastic materials. The generalized PDF is shown to include factors that are dependent on both the geometry of the material region as well as the statistical properties of the material. The first term gives a measure of the volume of allowable center points for the next inclusion along a one-dimensional ray, accounting for all edge effects due to the material boundary. This term can be solved analytically or numerically, if the shape of the material domain is known. The second term describes the probability of finding zero inclusions within a given volume. When expressed on a normalized, per volume basis, this term is shown to be related to the familiar rate parameter that characterizes a Poisson process. Several methods for approximating this rate parameter in one-, two-, and three-dimensional systems have been identified and their accuracy assessed. In one-dimensional problems, it has been shown that the generalized form of the PDF reduces to the exponential distribution, under certain simplifying assumptions. This result supports the use of an exponential distribution when sampling the distances between fixed-width inclusions in one-dimensional problems. However, in two- and three-dimensional problems, the newly developed generalized PDF shows a marked deviation from the

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

exponential distribution near material interfaces occurring within or at the edge of a stochastic material. To this end, we have demonstrated that the use of an exponential distribution to approximate the distribution of distances between adjacent inclusions will systematically overpredict the probability of distances less than one inclusion diameter. This over-prediction of short distances has been observed in the past, and is often attributed to edge effects related to the outer boundary of the material [7,8,11,12]. However, we further show that this nonexponential behavior for small distances can occur far away from material boundaries (or even in infinite problems), due to boundary-layer effects from the inclusions themselves. A discussion of the application of this newly developed generalized PDF has been provided along with supporting numerical results for case studies in one- and twodimensions. These numerical results demonstrate the ability of the newly derived PDF to correctly account for edge effects in finite stochastic materials, while still reproducing the expected distribution within the bulk material region.

Acknowledgements The authors would like to acknowledge Dr. Tom Sutton and Dr. Tim Donovan of the Knolls Atomic Power Laboratory and Prof. Ziya Akcasu of the University of Michigan for their helpful comments and suggestions regarding the detailed analysis of one-dimensional materials presented in this paper.

space. Consequently, Eq. (A.1) indicates that the distances between inclusions should be exponentially distributed. This prediction for one-dimensional stochastic materials is known to be exact for ideal (infinite medium) cases [12,14], and has also been shown to produce accurate results in certain non-ideal (typically thick media) problems [5,10,11]. In this analysis we will consider the validity of this assumed distribution for different combinations of slab length (L) and inclusion thickness (Dx). We proceed by deriving a probability density function for the number of inclusions in a slab of total length L, based on the assumed distribution of distances between adjacent inclusions, as given in Eq. (A.1).

A.1. Single-inclusion case Before deriving the distribution for a slab with multiple inclusions, let us first consider a simpler two-region problem that contains a random-length segment of material A followed by a single inclusion of material B. In this case, the length of the material A segment, x1, is a realization of the random variable X, which is governed by the probability distribution p(x; l). The total length of the simple two-region problem is equal to l1 ¼ x1 þ Dx:

ðA:2Þ

The probability of finding no inclusions in a fixed length L is equal to the probability that l1 4L, or, correspondingly, the probability that x1 4L  Dx. From the probability distribution given in Eq. (A.1) we can show pðx 4 LDxÞ ¼

Appendix A. Derivation of p(N) for a 1-D slab

595

Z1

lelx dx:

ðA:3Þ

LDx

For this analysis we will consider a one-dimensional slab of length L, which contains a number of inclusions of fixed length Dx distributed randomly throughout. The background material of the slab is identified as material A, and the material in the inclusions is identified as material B. Both the number and spatial location of the inclusions is assumed to be random, subject to the conditions that no inclusion may extend beyond the boundaries of the slab and no two inclusions may overlap with one another. An illustration of an example slab of this type is shown in Fig. 5. In this subsection we seek to derive statistical distributions that describe both the number of inclusions in a slab, and also the distribution of distances that separate two adjacent inclusions. We begin the derivations by assuming that the distance between any two adjacent inclusions can be well-described by the 1-D form of the generalized PDF derived in Section 2.4, Eq. (8) pðx; lÞ ¼ lelx ,

ðA:1Þ

where x is the distance between inclusions and l is a constant rate parameter, which gives the expected number of inclusion center points encountered per distance traveled. As discussed in Section 2.4, the derivation of Eq. (A.1) relies on the initial assumption that the inclusion center points follow a Poisson distribution in

Evaluating the integral in Eq. (A.3) yields an expression for the probability pðN ¼ 0Þ ¼ elðLDxÞ ,

ðA:4Þ

where N is the number of inclusions in the slab. We note that Eq. (A.4) agrees with the N =0 result previously derived by Donovan and Danon [5]. In cases where the total length of the slab is so narrow that more than one inclusion cannot fit (Lo2 Dx), we can immediately calculate the probability of finding an inclusion in the slab by the relationship, p(N = 1)= 1 p(N= 0). However, for cases where L42 Dx we must consider the possibility of multiple inclusions occurring within the slab.

A.2. Multiple-inclusion case For slabs containing more than one inclusion, we immediately recognize that any slab containing multiple inclusions may be broken up into a series of two-region segments of the form described in Section 3.2.1. Under this partitioning, the length of the material A region in each segment is an independent realization from the random variable X. Therefore, the total length of n segments, Ln, can be calculated by summing Eq. (A.2)

596

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

Substituting Eq. (A.9) and (A.10) into Eq. (A.8) yields

over n independent realizations from X Ln ¼

n X

ðxi þ DxÞ:

ðA:5Þ

Z1

pðN rn1Þ ¼ pðyn 4 LnDxÞ ¼

i¼1

pEr ðy; n, lÞ dy:

ðA:13Þ

LnDx

Since the length of the inclusions is a constant value, we may factor this term out of the summation and rewrite Eq. (A.5) as n X

Ln ¼ nDx þ

xi :

The right-hand side of Eq. (A.13) can be rewritten in terms of the cumulative probability distribution, PEr, by Z LnDx Z 1 pEr ðy; n, lÞ dy pEr ðy; n, lÞ dy pðN rn1Þ ¼

ðA:6Þ

i¼1

In cases where the inclusion widths are not constant, an equivalent result to Eq. (A.6) can be obtained by replacing Dx with the expected inclusion width for the material. Because the realizations of X are guaranteed to be positive, Ln must be a monotonically increasing function of n. Therefore, for any given sequence of realizations {x1, x2, x3,y} there exists a value n0 such that Ln0 + j 4L for all j Z1. Thus n0 gives the number of inclusions that occur in a slab of length L for a given sequence of realizations {x1, x2, x3,y}. Additional inclusions are not possible because they will not fit within the given slab length, L. Using this information, we wish to determine the probability that a given number of inclusions, N, will occur within a slab of length L. We begin by determining the probability that Ln 4L, implying that n inclusions will not fit in a slab of length L. From Eq. (A.6) we write ! n X xi 4L : ðA:7Þ pðLn 4LÞ ¼ p nDx þ

0

0

pðN rn1Þ ¼ 1PEr ðy; n, lðLnDxÞÞ:

ðA:14Þ

Substituting Eq. (A.11) into Eq. (A.14) gives the cumulative distribution function Pðn1Þ ¼ pðN rn1Þ ¼ 1

gðn, lðLnDxÞÞ ðn1Þ!

:

ðA:15Þ

Now that we have determined the probability that there are n 1 or fewer inclusions we can determine the probability that there are exactly n inclusions by noting that pðN ¼ nÞ ¼ pðN rnÞpðN rn1Þ:

ðA:16Þ

Substituting Eq. (A.15) for N rn and N rn  1 into Eq. (A.16) gives the final result pðNÞ ¼

gðN, lDxðmNÞÞ gðN þ 1, lDxððm1ÞNÞÞ 

ðN1Þ!

N!

,

ðA:17Þ

where m is the dimensionless slab length, which is defined as the ratio of the slab length to the inclusion width, or m  L/Dx.

i¼1

From the left-hand side of Eq. (A.7) we note that when Ln 4L, the actual number of inclusions in the slab, N, must be less than or equal to n  1 ! n X xi 4LnDx : ðA:8Þ pðN rn1Þ ¼ p i¼1

The right-hand side of Eq. (A.8) is the probability that the sum of n independent identically distributed random variables will be greater than a constant value. Writing this sum as a new random variable Yn, we have yn ¼

n X

xi :

ðA:9Þ

i¼1

In the case where the xi realizations are all taken from an identical exponential distribution with rate parameter l, it can be shown that the yn values are Erlang distributed with probability pEr ðy; n, lÞ ¼

ln yn1 ely ðn1Þ!

,

gðn, lyÞ ðn1Þ!

,

pðN ¼ 0Þ ¼ elDxðm1Þ ,

ðA:11Þ

ðA:4Þ

which can be rewritten in terms of the incomplete gamma function, Eq. (A.12), to yield pðN ¼ 0Þ ¼ 1gð1, lDxðm1ÞÞ:

ðA:18Þ

Understanding the behavior of p(N) for N 4m 1 requires a little more effort. Returning to the cumulative distribution function defined in Eq. (A.15) we can show that, for cases where the slab length is an integer multiple of the inclusion thickness (m 2 Z þ , where Z þ is the set of positive integers) pðN rm1Þ ¼ 1

where g(n, ly) is the lower incomplete gamma function, which is defined by Z ly gðn, lyÞ ¼ t n1 et dt: ðA:12Þ 0

By inspection we see that the probability distribution function given in Eq. (A.17) is not well defined for cases where N = 0 (due to N  1 factorial in first term) or N 4m  1 (due to the argument of the gamma function in the second term). The expression for p(N = 0) was derived in Section A.1 (Eq. (A.4)), and shown to be

ðA:10Þ

and cumulative probability PEr ðy; n, lÞ ¼

A.3. General solution for p(N)

gðm,0Þ m!

¼1

ðA:19Þ

by the definition of the lower incomplete gamma function, given in Eq. (A12). The result shown in Eq. (A.19) demonstrates that, for cases where the slab length, L, is an integer multiple of the inclusion length, Dx, the probability of observing N 4m 1 inclusions is zero. For cases where the slab length is not an integer multiple of the inclusion length, the corresponding analysis becomes more complicated. Because the probability distribution for

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

and g to be rewritten as

the number of inclusions (Eq. (A.15)) is a discrete valued function, we cannot calculate the probability that N is less than or equal to a non-integer value. However, we can calculate the probability that N is less than or equal to the nearest integer below m 1. For this case we can show that pðN r bmc1Þ ¼ 1 pðN r bmc1Þ ¼ 1

gðbmc, lDxðmbmcÞÞ ðbmcÞ!

gðbmc, lDxf Þ ðbmcÞ!

k X

ðA:20Þ

bmc X

bmc X

Given an analytical formula for the probability of finding exactly N inclusions in a slab of total length L, Eq. (12), the expected number of inclusions in the slab can be calculated as /NS ¼

1 X

n pðnÞ:

ðB:1Þ

n¼0

Applying the condition that p(n)= 0 for n 4 bmc allows Eq. (B.1) to be rewritten as a finite sum /NS ¼

bmc X

n pðnÞ:

ðB:2Þ

n¼0

Eq. (B.2) is in a convenient form for a summation by parts transformation, which allows a finite sum of two factors f

n pðnÞ ¼ bmc

bmc X

pðjÞ

gj :

ðB:3Þ

j¼0

bm c1 X

ðn þ1nÞ

n¼0

j¼0

n pðnÞ ¼ bmc

n X

pðjÞ:

ðB:4Þ

j¼0

bm c1 X

PðnÞ:

ðB:5Þ

n¼0

n pðnÞ ¼

n¼0

bm c1 X

1PðnÞ:

ðB:6Þ

n¼0

Substituting the definition for the cumulative distribution function, Eq. (A.15), into Eq. (B.6) gives a compact formula for the expected number of inclusions in the slab

when m 2 Z. Using the results given in Eqs. (A.17), (A.18), (A.21) and (A.22), we can now write the probability distribution for the number of inclusions in a slab of total length L as

Appendix B. Derivation of expected number of inclusions in a 1-D slab

n X

Factoring the first term on the right-hand side of Eq. (B.5) into the summation yields

ðA:22Þ

Before continuing, we wish to take a moment to remind the reader that the preceding derivations and resulting probability distribution function are exact if the distances between adjacent inclusions are exponentially distributed with constant rate parameter l.

n¼0

n¼0

for all cases, and

ðA:23Þ

ðfn þ 1 fn Þ

By noting that the first summation over j is equal to 1, by Eq. (A.21), and the second summation over j is the cumulative distribution function P(n), we can rewrite Eq. (B.4) as

ðA:21Þ

8 1gð1, lDxðm1ÞÞ, N ¼ 0, > > > < gðN, lDxðmNÞÞ gðN þ 1, lDxððm1ÞNÞÞ ( 0 o N r m1, m 2 Z,  , pðNÞ ¼ 0 oN r bmc, m= 2Z, ðN1Þ! N! > > > : 0, Other

k1 X

gj 

j¼0

n¼0

where bmc denotes the closest integer below the real value m, and f ¼ ðmbmcÞ is the fractional part of m, which is guaranteed to fall between 0 and 1. The result shown in Eq. (A.20) demonstrates that there is a non-zero probability of observing N 4 bmc1 inclusions in the slab. In addition, we can also recognize that it is physically impossible to allow N Z bmc þ 1 inclusions in the slab, as the total length of the bmc þ 1 inclusions is Dxðbmc þ 1Þ ¼ L þ Dxð1f Þ, which is greater than the length of the slab, L. Therefore, we conclude that the probability of observing N 4 bmc inclusions must be zero. From the above discussion, it is clear that

pðN 4 bmc1Þ ¼ 0

k X

Applying Eq. (B.3) with fn = n and gn = p(n) allows Eq. (B.2) to be rewritten as bmc X

pðN 4 bmcÞ ¼ 0

fn gn ¼ fk

n¼0

,

o1,

597

hN i ¼

bm c1 X

gðn þ1, lDxðmðn þ 1ÞÞÞ

n¼0

n!

:

ðB:7Þ

References

:

[1] Adams ML, Larsen EW, Pomraning GC. Benchmark results for particle transport in a binary Markov statistical medium. J Quant Spectrosc Radiat Transfer 1989;42(4):253–66. [2] Pomraning GC. Transport theory in stochastic mixtures. Trans Am Nucl Soc 1991;64:286–7. [3] Zimmerman GB, Adams ML. Algorithms for Monte Carlo particle transport in binary stochastic mixtures. Trans Am Nucl Soc 1991;64:287–8. [4] Donovan TJ, Danon Y. HTGR unit fuel pebble k-infinity results using chord length sampling. Trans Am Nucl Soc 2003;89:291–2. [5] Donovan TJ, Danon Y. Application of Monte Carlo chord-length sampling algorithms to transport through a two-dimensional binary stochastic mixture. Nucl Sci Eng 2003;143:226–39. [6] Reinert DR. Investigation of stochastic radiation transport methods in random heterogeneous mixtures. PhD dissertation, University of Texas, Austin, 2008. [7] Murata I, Mori T, Nakagawa M. Continuous energy Monte Carlo calculations of randomly distributed spherical fuels in hightemperature gas-cooled reactors based on a statistical geometry model. Nucl Sci Eng 1996;123:96–109. [8] Murata I, et al. New sampling method in continuous energy Monte Carlo calculation for pebble bed reactors. J Nucl Sci Technol 1997;34:734–44. [9] Brown F, Martin WR. Stochastic geometry capability in MCNP5 for the analysis of particle fuel. Ann Nucl Energy 2004;31:2039–47. [10] Ji W, Martin WR. Monte Carlo simulation of VHTR particle fuel with chord length sampling. In: Proceedings of the joint international topical meeting on mathematics & computation and supercomputing in nuclear applications (M&C + SNA 2007), American Nuclear Society, CD-ROM, 2007. [11] Ji W, Martin WR. Application of chord length sampling to VHTR unit cell analysis. In: Proceedings of the international conference on the physics of reactors (PHYSOR 2008), American Nuclear Society, CDROM, 2008. [12] Olson GL, Miller DS, Larsen EW, Morel JE. Chord length distributions in binary stochastic media in two and three dimensions. J Quant Spectrosc Radiat Transfer 2006;101:269–83.

598

D.P. Griesheimer et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 577–598

[13] Su B, Pomraning GC. A stochastic description of a broken cloud field. J Atmos Sci 1994;51(13):1969–77. [14] Torquato S, Lu B. Chord-length distribution function for two-phase random media. Phys Rev E 1993;47(4):2950–3. [15] Lu B, Torquato S. Lineal-path function for random heterogeneous materials. Phys Rev A 1992;45(2):922–9. [16] Torquato S, Lu B. Nearest-neighbor distribution functions in manybody systems. Phys Rev A 1990;41(4):2059–75. [17] Torquato S, Lee SB. Computer simulations of nearest-neighbor distribution functions and related quantities for hard-sphere systems. Physica A 1990;167:361–83. [18] Torquato S, Lu B, Rubinstein J. Nearest-neighbor distribution function for systems of interacting particles. J Phys A 1990;23:L103–7. [19] Gue´ron D, Mazzolo A. Properties of chord length distributions across ordered and disordered packing of hard disks. Phys Rev E 2003;68:066117.

[20] Davis AB, Marshak A. Photon propagation in heterogeneous optical media with spatial correlations: enhanced mean-free-paths and wider-than-expected free-path distributions. J Quant Spectrosc Radiat Transfer 2004;84:3–34. [21] Widom B. Random sequential addition of hard spheres to a volume. J Chem Phys 1966;44(10):3888–94. [22] Schaaf P, Talbot J. Kinetics of random sequential adsorption. Phys Rev Lett 1989;62(2):175–8. [23] Tarjus G, Schaaf P, Talbot J. Random sequential addition: a distribution function approach. J Statist Phys 1991;63(1/2): 167–202. [24] Griesheimer DP, Millman DL. Analysis of distances between inclusions in finite one-dimensional binary stochastic materials. In: Proceedings of the international conference on mathematics, computational methods & reactor physics (M&C 2009), American Nuclear Society, CD-ROM, 2007.