- Email: [email protected]
Qualitative analysis of fiber composite microstructure: Influence of boundary conditions
Probabilistic Engineering Mechanics 21 (2006) 317–329 www.elsevier.com/locate/probengmech
Qualitative analysis of fiber composite microstructure: Influence of boundary conditions ˇ Jan Gajdoˇs´ık, Jan Zeman ∗ , Michal Sejnoha Department of Structural Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Th´akurova 7, 166 29 Prague 6, Czech Republic Received 16 August 2005; received in revised form 11 November 2005; accepted 21 November 2005 Available online 9 January 2006
Abstract The principal objective of this work is to classify the influence of different boundary conditions that can be applied when computing various statistical descriptors for generally random microstructures. Although applicable to most statistical descriptors, attention is limited to twopoint probability functions. Binary images of random fibrous composites are assumed in the present study. Here, the periodic, mirror and no (plain) boundary conditions are addressed when processing the binary images. Also, the minimum number of samples (Representative Volume Elements—RVEs) in the ensemble that still provide statistically representative results is sought. While the available results promote the use of periodic boundary conditions, particularly from the computational point of view regardless of the RVE size, the number of samples required is largely affected by their size. If available, the most efficient procedure results from the use of a large binary image of a real microstructure combined with periodic boundary conditions. c 2005 Elsevier Ltd. All rights reserved.
Keywords: Fibrous composites; Quantification of microstructure morphology; Two-point probability function; Boundary conditions; Representative volume element
1. Introduction to problem of microstructure analysis The reliable and accurate analysis of heterogeneous material systems increasingly often takes advantage of knowledge of real microstructures or even nanostructures; see, e.g., [1–3]. A specific example is given in Fig. 1, illustrating a random distribution of graphite fibers within a single bundle of large woven composite tube [4]. The periodic character of such systems, though often assumed, is generally not available. The random nature of most material systems at higher levels of resolution is therefore naturally accepted. Such a step, however, calls for a suitable method of microstructure1 quantification. A number of contributions have advocated the use of various statistical descriptors such as the one- and two-point probability functions or the lineal path function (see [5–8], to cite a few).
∗ Corresponding author. Tel.: +420 2 2435 4482; fax: +420 2 2431 0775.
E-mail addresses: [email protected] (J. Gajdoˇs´ık), ˇ [email protected] (J. Zeman), [email protected] (M. Sejnoha). 1 For simplicity, all structural levels, starting from the nanoscale up to the meso-scale, that resemble some kind of randomness, including random imperfections, will be jointly referred to as microstructure. c 2005 Elsevier Ltd. All rights reserved. 0266-8920/$ - see front matter doi:10.1016/j.probengmech.2005.11.006
In their recent works, the authors exercised the use of loworder statistical descriptors in a number of problems linked to the analysis of material systems with random or imperfect geometries. The formulation of statistically equivalent periodic unit cells, derived by matching the material statistics up to two-point probability function, has been offered not only for classical fibrous composites [9,10] but also for complex threedimensional textile composites with various sources of random imperfections [11], and even for masonry structures with an irregular arrangement of stone blocks [12]. Direct introduction of the two-point probability function into variational principles of Hashin and Shtrikman allowed a simplified yet reasonably accurate and extremely efficient analysis of complex material systems, not only in elastic [13] but also in visco- and nonlinear visco-elastic regimes [8,14]. In all applications, the periodic boundary conditions were intuitively selected when processing the binary images of various sizes taken as random cuts from images of true microstructures; Fig. 1(a). The resulting anisotropic response of such systems gave rise to various objections that question the assumed approach. Clearly, a number of errors may occur in the process of transforming the original color image into its binary
318
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Fig. 1. (a) A fiber tow, (b) a real micrograph of a selected cut.
counterpart. While relatively high resolution color images are available through the application of various image input devices, including the direct use of a digital camera or a camera connected to a microscope, the quality of the image itself depends on the quality of the original sample. Quite often, such specimens are of relatively poor quality, requiring manual adjustments of the resulting image for further analysis. The most critical region quite often appears in the vicinity of the sample boundary where the images are usually incomplete, e.g. “cut” fibers in the case of fibrous composites. This is particularly important when studying a large ensemble of small samples with imperfect boundaries created when taking random cuts. Investigating the image within this region leads to a number of problems. Note that the two-point probability function is defined as the probability of finding two specified points randomly thrown into a medium in specified phases. The finite dimensions of the image pose a clear question of how to treat cases, when one of the points is found out of the image boundary. Such a pair can be either directly disposed or the distribution of microstructure outside the image can be predicted in some way. Comparison of various techniques dealing with image boundaries supplies the major part of this contribution. In addition, the effect of the image size on the analyzed statistical descriptor is also studied. This problem is addressed particularly from the efficiency point of view, comparing the analysis on a single large sample (with
reference to ergodicity assumption) or a large ensemble of small samples. The minimum number of required samples given their size is provided. An additional constraint on the minimum size of RVE linked to the issue of statistical independence applied in [8] can also be put forward. This condition, however, was not examined in the present contribution. The present study is closely linked to the analysis of apparent linear mechanical properties of random composite materials initiated by Huet [15]. In this approach, to compensate for a limited size of an analyzed material sample, the microstructure is subjected to kinematic, static and periodic boundary conditions. With increasing size of the material sample, the response corresponding to these boundary conditions converges to a common value determining the size of the RVE; see also [16] for a detailed discussion and [17] for a recent exhaustive review. This procedure has been performed numerically for different material systems (see, e.g., [18–21] to cite a few) and all these studies confirm the superiority of periodic boundary conditions, even for non-periodic material systems. This conjecture lead to the recent introduction of a concept of periodization of random media [22]. This paper can be seen as an extension of these studies to quantification of the morphology of random microstructures. Finally, note that recent works on the definition of RVE for heterogeneous materials explicitly involve spatial statistics (up to two-point probability functions) for the definition of RVE [22,23]. For a moderate contrast between constituents’ properties, it seems that the condition of convergence in terms of a two-point probability function is more restrictive than the convergence of apparent material properties. The paper is organized as follows. The introductory part, Section 1, is followed by the specification of various statistical descriptors generally used in microstructure quantification. Several possible approaches for their evaluation, based on specific boundary conditions, are offered. Image analysis of a section of a graphite fiber tow is briefly discussed in Section 3. The results comparing various boundary conditions also from the RVE-size stand point are summarized in Section 4. 2. Microstructure quantification—statistical measures A number of statistical measures have been examined in the last few decades for addressing the problem of heterogeneous media with either a random distribution of material phases or the presence of various types of material or geometrical imperfections of more or less random nature. Starting from the characteristic function as the basic descriptor for microstructure quantification, the literature offers a number of derived statistical functions, including the radial path function [24], the lineal path function [25] or the general n-point probability function [6]. Of them all, the two-point probability function appears to be the most prominent, owing to the ease and efficiency of its evaluation, particularly when taking advantage of the Fast Fourier transform (FFT) applied to the binary images of real microstructures [26]. To prevent the reader from searching a vast body of literature devoted to this function, we now briefly summarize its most important features and the
319
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
essential statistical assumptions allowing its relatively efficient evaluation even for non-periodic media.
Further simplification arises when assuming the medium to be statistically isotropic [26]. Then Sr s (x1 , x2 ) reduces to
2.1. Fundamental function and statistical moments
Sr s (x1 − x2 ) = Sr s (k x1 − x2 k).
Consider an ensemble of a two-phase random medium. To provide a general statistical description of such a system, it is useful to characterize each member of the ensemble by a stochastic function — characteristic function χr (x, α), which is equal to one when point x lies in the phase r of the sample α and equal to zero otherwise [27,6,28], 1, if x ∈ Dr (α), χr (x, α) = (1) 0, otherwise, where Dr (α) denotes the domain occupied by the r -th phase. Except where noted, composites consisting of a clearly distinguishable continuous matrix phase are considered. Therefore, r = m, f , is further assumed to take values m for the matrix phase, while symbol f is reserved for the second phase. For such a system the characteristic functions χm (x, α) and χ f (x, α) are related by χm (x, α) + χ f (x, α) = 1.
(2)
Following [27,28,6], we write the ensemble average of the product of characteristic functions Sr1 ,...,rn (x1 , . . . , xn ) = χr1 (x1 , α) · · · χrn (xn , α),
(3)
where function Sr 1 , . . . , rn , referred to as the general npoint probability, gives the probability of finding n points x1 , . . . , xn randomly thrown into a medium located in the phases r1 , . . . , rn . 2.2. Functions of the first and second order Hereafter, we limit our attention to functions of the order of one and two, since higher-order functions are quite difficult to determine in practice.2 Therefore, the description of a random medium will be provided by the one-point probability function Sr (x): Sr (x) = χr (x, α),
(4)
which simply gives the probability of finding the phase r at x, and by the two-point probability function Sr s (x1 , x2 ) Sr s (x1 , x2 ) = χr (x1 , α)χs (x2 , α),
(5)
which denotes the probability of finding simultaneously the phase r at x1 and the phase s at x2 . In general, evaluation of these characteristics may prove to be prohibitively difficult. Fortunately, a substantial simplification applies when accepting an assumption regarding the material as being statistically homogeneous, so that Sr (x) = Sr , Sr s (x1 , x2 ) = Sr s (x1 − x2 ).
(6) (7)
2 Note, however, that relatively efficient procedures for the approximation of higher-order probability functions for ergodic and statistically isotropic media were proposed in [29].
(8)
Finally, making an ergodic assumption allows a substitution of the one-point correlation function by its volume average [26], i.e., volume concentration or volume fraction of the r -th phase cr : Sr = cr .
(9)
2.3. Limiting values In addition, the two-point probability function Sr s incorporates the one-point probability function Sr for certain values of its arguments such that for x1 = x2 : Sr s (x1 , x2 ) = δr s Sr (x1 ), for kx1 − x2 k → ∞ : lim Sr s (x1 , x2 )
(10)
kx1 −x2 k→∞
= Sr (x1 )Ss (x2 ),
(11)
where symbol δr s stands for Kronecker’s delta. Relation (10) states that the probability of finding two different phases at a single point is equal to 0 (see also Eq. (2)) or is given by the one-point probability function if phases are identical. Eq. (11) manifests that, for large distances, points x1 and x2 are statistically independent. This relation is often denoted as the no-long-range-order hypothesis (see e.g. [30,5]). Finally, according to Eq. (2), we may determine one- and two-point probability functions for all phases, provided that these functions are given for one arbitrary phase. For the onepoint probability function of a statistically homogeneous and ergodic medium, this relation assumes a trivial form: cm = 1 − c f .
(12)
Relations for the two-point probability functions of statistically uniform and ergodic medium are summarized in Table 1.3 2.4. Numerical evaluation To determine Sr1 ,...,rn , we recall that the general n-point probability gives the probability of finding n points x1 , . . . , xn randomly thrown into a medium located in the phases r1 , . . . , rn . In view of Table 1, we further consider only the fiber probability functions. To follow the above definition, the one-point fiber probability function S f gives the chance of finding a randomly placed point located in the fiber phase. To determine this quantity, a simple Monte Carlo-like simulation can be utilized — we randomly throw a point into the microstructure and count successful “hits” into the matrix phase. Then, the value of function S f can be estimated as Sf ≈
n0 , n
(13)
3 Note that, by definition (5) and the assumption of statistical homogeneity, Sr s (x) = Ssr (x).
320
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Table 1 Relations between two-point probability functions Known function
Smm (x)
Sm f (x)
S f f (x)
Smm (x) Sm f (x) S f f (x)
Smm (x) cm − Smm (x) c f − cm + Smm (x)
cm − Sm f (x) Sm f (x) c f − cm − Sm f (x)
cm − c f + S f f (x) c f − S f f (x) S f f (x)
Fig. 2. (a) Microstructure with periodic boundary conditions, (b) microstructure with mirror boundary conditions, and (c) microstructure with no (or “plain”) boundary conditions.
Fig. 3. (a) Original binary image with no modifications. (b) Modified binary image. (c) Regenerated “ideal” binary image.
where n 0 is the number of successful hits and n denotes the total number of throws. An entirely similar procedure can be employed to determine the values of S f f (x).4 Another more attractive approach is available when the real microstructure is replaced by its binary image. A sample of a binary image is shown in Fig. 3(c). Such a digitized micrograph can be imagined as a discretization of the characteristic function χr (x), usually presented in terms of a W × H bitmap. Denoting the value of χr for the pixel located in the i-th row and jth column as χr (i, j) allows the writing of the first moments of function χr for an ergodic and statistically homogeneous
4 For a statistically isotropic microstructure, Smith and Torquato [31] proposed a more efficient procedure for the determination of Smm (k x12 k). Instead of tossing a line corresponding to x into a medium, a sampling template is used for the determination of the two-point probability function. See also [9] for a comparison of this method with bitmap-based approaches presented hereafter and Section 4.2 for further application of this idea.
medium in the form5 Sr =
−1 H −1 X 1 WX χr (i, j). W H i=0 j=0
Formally, one can write the two-point probability function Sr s using the similar relation Sr s (m, n) =
−1 −1 H X 1 WX χr (i, j)χs (i + m, j + n), W H i=0 j=0
(14)
which, however, can generate indexes in the characteristic function χs exceeding the bitmap dimensions. To treat such a case, the boundary conditions depicted in Fig. 2 can be applied.
5 Throughout the text, the C-language type of array indexing is consistently used, i.e., we denote the first element of an array α as α0 and the last element of the array as α L−1 , where L is the array length.
321
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
2.4.1. Periodic boundary conditions In the case of periodic boundary conditions, Fig. 2(a), the starting image of a microstructure, e.g. the one displayed in Fig. 3(c), is assumed to be periodically translated in all directions. The two-point probability function then assumes the form Sr s (m, n) =
−1 H −1 X 1 WX χr (i, j)χs ((i + m)%W, W H i=0 j=0
( j + n)%H ),
Table 2 Bitmap dimensions in pixels Height
Width
Size
No. of bitmaps
1000 850 700 550 400 250
1000 850 700 550 400 250
1 000 000 722 500 490 000 302 500 160 000 62 500
18 × 1 18 × 10 18 × 20 18 × 30 18 × 40 18 × 50
(15)
where the symbol “%” stands for modulo. The above equation, usually termed the cyclic correlation [32], readily implies periodicity of function Sr s . Note that the correlation property of DFT holds for cyclic correlation. Referring to [32], it is given by the following relation DFT{Sr s (m, n)} = DFT{χr (m, n)}DFT{χs (m, n)},
(16)
where · now denotes the complex conjugate. The inverse DFT, denoted IDFT, then serves to derive function Sr s at the final set of discrete points as [33] 1 IDFT{DFT{χr (m, n)}DFT{χs (m, n)}}. (17) WH This method is very economical and its accuracy depends only on the selected resolution of the digitized medium. Usually, the Fast Fourier Transform, which needs only O(W H log(W H ) + W H ) operations, is called on to carry out the numerical computation.6 This approach is extremely efficient, especially when comparing with traditional “Monte Carlo”-like or “trial and error” methods, which need to be employed for the remaining boundary conditions. Sr s (m, n) =
2.4.2. Mirror boundary conditions In this particular case, the microstructure is expected to be mirrored; Fig. 2(b). When taking this route we are bound to modify Eq. (14) in the form Sr s (m, n) =
−1 H −1 X 1 WX χr (i, j)χs (i\ + m, j[ + n), W H i=0 j=0
(18)
where i +m for 0 ≤ (i + m) < W i\ +m = 2W − 1 − (i + m) for W ≤ (i + m) < 2W j +n for 0 ≤ ( j + n) < H j[ +n = 2H − 1 − ( j + n) for H ≤ ( j + n) < 2H. Observe that to O((W H )2 ) operations are needed to evaluate the function Sr s . This might be computationally demanding, particularly for a large micrograph. These boundary conditions are therefore selected only for the purpose of comparison. Nevertheless, there are still other fields of research where this approach may prove to be useful. 6 The public-domain package FFTW version 2.1.3 [34] was used for the evaluation of Eq. (17).
Fig. 4. Graphical representation of the process of bitmap selection.
2.4.3. Plain boundary conditions In the third case, here termed the plain boundary condition, we essentially disregard the cases for which χs is located outside the image; see Fig. 2(c). The particular form of the twopoint probability function follows from Sr s (m, n) =
(i M
iX N −1 M −1 jX 1 χr (i, j) − i m )( j N − jn ) i=i j= j
χs (i + m, j + n),
m
n
(19)
where i m = max(0, −m), i M = min(W, W − m) and jn = max(0, −n), j N = min(H, H − n). Although the computational complexity of this approach is comparable to the mirror boundary condition, it is the most “suitable” one, as it introduces not a prior assumption about the sample boundary7 but significantly reduces the amount of information obtained. It is evident that long vectors (m, n) can be positioned just a few times to be totally contained in the image. If the vector with the length of the rectangle’s diagonal is used, it can be thrown just once. It becomes clear that results for large spacing of the two points yield lower reliability when compared with those with small spacing. On the other hand, it has been demonstrated in, see e.g. [8,9], that, for distances exceeding approximately 10 times the fiber diameter, a random throw of the two points becomes a statistically independent event and therefore the results for large spacing can be safely disregarded. 7 More precisely, it can be shown that these boundary conditions deliver an unbiased estimator of two-point probability function; see, e.g., [35, Section 5.1.3].
322
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Fig. 5. Exemplary graph of measures M.
3. Image analysis of fibrous composites Image processing represents an important part of the micromechanical analysis of real world material systems allowing scientists, through computer based image acquisition, to visualize essential material data. Unfortunately, the original data, usually supplied in the form of color images, are often of poor quality, so that further treatment is needed. Such a tedious step typically involves manual marking of object borders, completion of distorted objects, extracting individual color spaces, filtration and final conversion into binary image. Although limiting the subsequent analysis to two-phase material systems, the two-dimensional binary data provide detailed information about microstructure configuration and may reveal various sources of imperfections for relatively complex and even three-dimensional structures such as textiles [11]. Note that evaluation of the basic statistical descriptors described in Section 2 is not limited to two-phase systems. At the present time the color images, transformed into 256 colors in gray scale, can be analyzed by applying the same tools as for the binary image data [36]. Nevertheless, owing to the complexity of the process of retrieving the images of a real composite, we settle in the present study for conventional binary images. Such a simulation of real microstructure is more than sufficient when studying the effect of boundary conditions on the evaluation of material statistics from Section 2. In view of our previous work, we continue with the analysis of graphite fiber tow impregnated by a polymer matrix [37,38,8,14]. All specimens subjected to individual steps of image processing were extracted from a single-layer woven composite tube supplied by Compo-Tech, Ltd. [4]. To adjust the original color image of Fig. 1, the image analyzer LUCIA G [39] was employed to obtain the desired binary (monochromatic) image data. The result of the first step appears in Fig. 3(a). The respective binary image was obtained through the process of gradual thresholding of the original color image by setting the range of color scales or intensity and saturation for objects, for which the output binary image should be black (or white). Clearly, such an image is not very suitable for further analysis as it displays, as a result of sample preparation, a number of partially damaged and
possibly connected fibers. Manual separation and completion of individual fibers applicable in the software LUCIA G resulted in an improved image, plotted in Fig. 3(b). The subsequent analysis of this image provided the basic information about the original microstructure, such as the number of objects and positions of their gravity centers. Although this image is relatively satisfactory and applicable for the evaluation of the desired statistical descriptors, it still shows fibers of various sizes and shapes that are far from circular. The “ideal” image in Fig. 3(c) was thus created by generating circular fibers placed in the gravity centers of original fibers. The fiber radius common for all fibers was selected such as to keep the original fiber volume fraction. Further details can be found in [36]. This final image was then exploited in Section 4 to compute the two-point probability function introduced in Section 2. The above approach was exercised to retrieve about 25 images from various sections of the fiber tow cross-section of Fig. 1(a). Each image of approximately 200 fibers is assumed to be representative of the actual material system, at least when accepting the ergodicity assumption tested, e.g. in [9]. Eventually, only 18 images out of all 25 were selected for the analysis. In particular, images with at least one fifth of the area being fiber free were eliminated. Those images clearly experienced some error, which was not captured when processing the image. Finally, all remaining images were adjusted by taking a smaller cut-out to reintroduce partial fibers along the “ideal” image boundary, Fig. 3(c), which were originally present in the image, Fig. 1(b) and Fig. 3(a), but removed when applying the modification steps. 4. Two-point probability function for various boundary conditions As stated in Section 2, probability functions of infinite order are needed to fully describe the material statistics of a random medium. This step, however, is computationally infeasible. Note that even the three-point probability function that appears in some variational formulations [27,30] requires enormous computational time for its evaluation. Consequently, only the two-point probability function will be considered hereafter, as it is still relatively easy to evaluate, especially when assuming a generally random medium. Recall that, by definition, the two-point probability function can be interpreted as throwing a needle of a certain length into the image and counting the number of successful hits, out of all of them, when both ends of the needle fall into the selected colors (the number of colors is equal to the number of assumed phases). It remains to determine what happens when one end of the needle falls out of the image. This event will be studied here, assuming three possible boundary conditions introduced in Section 2.4. 4.1. Two-point probability function via periodic boundary conditions The basic steps for the evaluation of the two-point probability function for binary and periodic media were
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
323
Fig. 6. Graphs of measures M derived for periodic boundary condition.
described in Section 2. The most straightforward approach is to perform an averaging of the two-point probability function for all 18 statistically similar samples forming the ensemble. It can be advantageous, however, to check whether evaluating the two-point probability function on a larger set of smaller samples could be more efficient. It remains to determine how many samples of smaller size, taken as random cuts from the original “large” samples, are then required to achieve the desired accuracy provided by the large samples, if at all possible. To that end, the following strategy was followed. First, the largest sample of size 1000 × 1000 pixels was examined. In the next step the size was reduced by 150 pixels in both directions to arrive at a bitmap of 850 × 850 pixels. This step appeared reasonable, particularly with reference to the average fiber size of about 70 pixels. The smaller image was selected 10 times from random positions of the largest image. In each subsequent step the size was reduced again by 150 pixels in both directions and the random selection was increased by 10. The resulting bitmaps subjected to further study are listed in Table 2. Recall that this process was carried out 18 times for each of the largest bitmap. A graphical representation of this process appears in Fig. 4. The size of the pools was selected from the point of view of computational feasibility, especially when referring to Eqs. (18) and (19) for the two-point probability function evaluation. To compare individual results, the following measure was adopted:
M=
1 b·H b W
·
b−1 b −1 H W X X i=0 j=0
|Sr s (i, j) − Pr s (i, j)| . |Sr s (i, j)|
(20)
b and H b represent the number of points used in In Eq. (20), W the width and height directions of the bitmap, respectively. If every point of the bitmap is taken into account for the measure b and H b correspond exactly to the width and evaluation, then W height of a given bitmap. Sr s (i, j) is the value of the twopoint probability function of an image evaluated at the point with coordinates i and j. Pr s (i, j) is the value of the twopoint probability function of the largest bitmap that is assumed to provide the “exact” solution. This enabled comparison of measures taken from bitmaps with variable sizes. The results of the evaluation of the two-point probability function were stored in the matrix of the same dimensions as the bitmap. Note that function Sr s (i, j) corresponds to an average value evaluated for a given number of bitmaps selected from each pool. In other words, not only the total set of bitmaps (e.g., 20 samples for the bitmap size equal to 700 × 700 pixels) for each pool was used, but also partial sets, selected sequentially to get smaller pools of 2, 3, 4, . . . N (with N being the number of bitmaps of a given size), were tested. The random selection of bitmaps from each pool for averaging purposes was omitted, as the pool itself was already generated randomly. The specific average value of Sr s (i, j) was obtained as an average of independently computed two-point probability functions at the position i j.
324
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Fig. 7. Evaluation of function S f f from Eq. (19): (a) selection of a size and vector orientation; (b) positioning of a vector into a bitmap; (c) effect of vector positioning; and (d) checking isotropic assumption.
Fig. 8. Graphs of measures M derived for mirror boundary condition.
The results appear in Fig. 6. For detailed explanation of the results presented, we refer the reader to Fig. 5. The selected
title for each graph corresponds to the results for the smallest bitmap size plotted in a given graph. Each line is constructed
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
325
Fig. 9. Graphs of measures M derived for plain boundary condition.
Fig. 10. Comparison of measure M for various boundary conditions.
from averages of the measure M computed for smaller sets of bitmaps available for a given pool. Therefore, for the first graph
of bitmap size 850 × 850 in Fig. 6, the first point represents a value of M derived from a single bitmap (more exactly an
326
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
average of 18 bitmaps, as there are 18 large bitmaps at hand to create smaller pools) while the last point stands for the value of measure M found as an average from 10 × 18 bitmaps (the largest number of bitmaps available for this bitmap size; see Table 2). Each of the remaining graphs thus stores the results provided by all bitmaps of larger sizes. It is evident that larger bitmaps enable results to be obtained with better reliability. Quite surprisingly, even for smaller bitmaps there is no common value where all measures completely stabilize. Also note that averages of the measure M for various bitmap sizes oscillate around different values. Nevertheless, all the graphs clearly show that the smaller the bitmap size, the larger the number of bitmaps needed to obtain acceptable results. For most of the assumed bitmap sizes, the number of 20 bitmaps appears sufficient, at least from the selected measure point of view. 4.2. Two-point probability function via mirror and plain boundary conditions Recall that, for the mirror and plain boundary conditions, the lack of periodicity prevents the use of FFT techniques for the evaluation of the two-point probability function, which leaves us with Eq. (19) to deal with. It has been suggested in [9,26] that the medium under consideration can be regarded as ergodic and statistically homogeneous. In such a case, the two-point probability function can be evaluated irrespective of the choice of the coordinate system origin. As such, it depends only on the absolute distance of the two points and the orientation of the line connecting the two points. In view of this, the polar coordinate system r, φ is typically adopted when plotting the function. In addition, due to the symmetry of the function, the angle φ sweeping the area from 0◦ to 180◦ is sufficient, which considerably reduces the computational time. Nevertheless, when compared to the solution of Eq. (17) that employs the FFT method, the solution of Eq. (19) for a statistically homogeneous medium is still extremely demanding computationally. Owing to the fact that our primary goal is not to analyze the current medium exactly, we settled, in view of the above comments, for the assumption of statistical isotropy, so that the solution no longer depends on the angle orientation. Such an assumption renders a relatively fast and simple algorithm for evaluation of the function. Individual steps can be inferred from Fig. 7. As shown in Fig. 7(a), the same grid of points corresponding to the resolution of the analyzed bitmap was created. The two-point probability function was then evaluated for a set of vectors connecting the coordinate system origin placed in the upper left corner and each point in the grid. The starting point of each vector was then successively placed into every point in the bitmap, as depicted in Fig. 7(b). The location for the end point was then checked, and a number of successful hits that met the required condition were recorded, e.g. for the fiber–fiber two-point probability function S f f both points are required to be found in the fiber phase (see Fig. 7(c); clearly, among all four vectors, only vector No. 4 is acceptable — the remaining three vectors must be disregarded). Dividing this value by the total number of hits (the number of positionings of the vector of a given length into the bitmap) gives the desired
Fig. 11. Distribution of function S f f : (a) derived from a single bitmap, (b) averaged over the maximum number of sample sets for a given bitmap.
value of the two-point probability function for a given vector length. See also Fig. 7(c) for a graphical explanation of the condition of statistical isotropy (vectors a, b and c are identical, and the value of S f f does not depend on their orientation; if only the condition of statistical homogeneity applies, then the two vectors a and b are identical but differ from vector c for the S f f function evaluation). Additional simplifying assumptions can be introduced to improve the efficiency of function evaluation. Two possibilities are at hand. First, the value of the two-point probability function can be computed for all available vectors of the function generated in Fig. 7(a), but not all the pixels in the bitmap are used for the function evaluation, Fig. 7(b). In other words, the needles of every available length are used, but just a limited number of throws is done. Another possibility is to evaluate the function just for a few specific vectors, but all data stored in the bitmap are used for the function evaluation. In other words, we throw the needles of just a few chosen lengths, but the needles are thrown into every point of the bitmap. In the present work the second possibility was adopted, as it allows direct comparison with results obtained by the Fast Fourier Transform method that evaluates the function using all data stored in the image. The selection of step size was based on the average size of a single fiber equal to approximately 70 pixels. This value promoted use of the step of 15 pixels, which is less than one quarter of the fiber diameter. The obtained results appear in Figs. 7(d)–9. Fig. 7(d) compares the variation of the two-point probability function
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
327
Fig. 12. Plane view of function S f f for various boundary conditions: (a) periodic boundary conditions, (b) mirror boundary conditions, (c) plain boundary conditions.
for mirror and plain boundary conditions derived from Eq. (19), assuming statistical isotropy with the results derived for the periodic boundary conditions under the assumption of statistical homogeneity and then averaged over all angles for a given distance r to obtain the necessary isotropized values. Evidently, the results are almost identical, particularly when comparing the periodic and plain boundary conditions. Figs. 8 and 9 display plots of the measure M derived for the same sample sizes as presented in Section 4.1. Similar concluding remarks to those already suggested in the previous section can be drawn from these graphs. 4.3. Comparison of individual boundary conditions Except for Fig. 7(d), the results discussed so far were presented with no reference to each other. To get better insight into the effect of the selected boundary conditions, the results from Figs. 6, 8 and 9 are re-plotted in Fig. 10 for the selected sample sizes. On each graph, only the results corresponding to the selected sample size are shown. Each line in the graph displays the measure M found for a specific boundary condition. The results are rather inconclusive, and it appears that the same number of samples of smaller bitmaps are needed to stabilize the solution. A valuable piece of information is also available from Fig. 11, showing the variation of the fiber–fiber two-point probability function derived from selected sample sets. The solid line represents the results derived from one of the largest 1000 × 1000 bitmaps selected as an average from the original ensemble of 18 bitmaps. The dashed and dotted lines were then found from sample sets of 20 (700 × 700) and 50 (400 × 400) bitmaps, respectively. Evidently, the smaller the sample size the shorter the corresponding line, as longer vectors do not fit into smaller cut-outs. As can be already deduced from the previous results, if samples of smaller sizes are used, it is strongly recommended to perform averaging over the number of samples in the ensemble to arrive at meaningful results; compare Figs. 11(a) and (b). The final comparison is presented in the form of twodimensional plots of the fiber–fiber two-point probability function. It follows from Figs. 12 and 13 that the difference between periodic and mirror or plain boundary conditions is random and does not follow any regular pattern. The assumed
Fig. 13. Difference in function S f f plotted in absolute values: (a) comparison of periodic and mirror boundary conditions, (b) comparison of periodic and plain boundary conditions.
isotropy is also evident. Clearly, the visible difference is shifted away from the origin where statistical dependency no longer applies. Close to the origin, where the distribution of the two-point probability function provides the most significant statistical information about the microstructure configuration, the difference is almost negligible. These final plots thus further support the fact that neither boundary condition is superior to the other, apart from the statement that the plain boundary conditions are capable of delivering statistically unbiased results. 5. Conclusion It can be concluded from the results presented that all boundary conditions can be considered as equal. Such a
328
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
conclusion follows not only from Fig. 10, showing a close match of measure M for various boundary conditions, but also from Fig. 13, suggesting negligible difference in the twopoint probability function, especially between periodic and plain boundary conditions. Recall also Fig. 7(d), advocating no difference between the periodic and plain boundary conditions. The similarity in the behavior of the selected measure M can also be inferred from Figs. 6 and 8–10, which further supports the suggested conclusion. Due to the fact that no boundary conditions are superior to any other, the only objective that remains is the computational complexity. While the use of periodic boundary conditions allows implementation of a very efficient algorithm based on the Fast Fourier Transform (FFT) method resulting in negligible computational time, the solution of Eq. (19) for mirror and plain boundary conditions is much more demanding. For the largest bitmap of 1000×1000 pixels, the FFT method is about 125 000 faster on a standard PC with a 1.6 GHz Athlon 3000+ processor. The final conclusion is therefore twofold: • As the computational time needed for the evaluation of Eq. (19) depends heavily on the bitmap dimensions, it is recommended that, if the plain boundary conditions are needed for some reason, it is more efficient to compute the two-point probability function from a larger set of smaller bitmaps. With reference to Fig. 11 their size, however, should be sufficiently large to receive enough statistical information about microstructure configuration. • If, on the other hand, the applicability of periodic boundary conditions is possible, it is reasonable to use a single bitmap of as large a size as possible. Although the results for a single large bitmap and sets of smaller bitmaps are at least qualitatively comparable, see Fig. 11, the results for larger bitmaps are still more accurate. A rather naive but illustrative explanation of all the presented results assumes that the two-point probability function behaves a little bit like the surface of water in a pond after being hit by a stone in the middle of the pond. If stones of different shapes, but quite similar in size, are thrown into a pond, the waves close to the place of impact are quite different but further from the place of impact the waves are small. If the wave is small, the difference is small too. If the observer is monitoring just the small area of the surface close to the place of impact, he will see that the differences between the shapes of waves made by stones of different shapes are bigger than if he monitors a larger area. The results from this research are the same. The average values of measure M are always smaller if determined from functions computed for larger bitmaps; compare Figs. 10(a) and (e). Acknowledgments ˇ No. 106/03/H150 and The support provided by the GACR CEZ MSM 6840770003 is gratefully acknowledged. References [1] Babuˇska I, Andersson B, Smith PJ, Levin K. Damage analysis of fiber composites Part I: Statistical analysis on fiber scale. Computer Methods in Applied Mechanics and Engineering 1999;172(1–4):27–77.
[2] Torquato S. Statistical description of microstructures. Annual Review of Materials Research 2002;32:77–111. [3] Raghavan P, Ghosh S. Adaptive multi-scale computational modeling of composite materials. CMES-Computer Modeling in Engineering & Sciences 2004;5(2):151–70. [4] Compo Tech Plus webpage. http://www.compotech.com. [5] Willis JR. Bounds and self-consistent estimates for the overall properties of anisotropic composites. Journal of the Mechanics and Physics of Solids 1977;25:185–202. [6] Torquato S, Stell G. Microstructure of two-phase random media. I. The n-point probability functions. Journal of Chemical Physics 1982;77(4): 2071–7. [7] Torquato S. Random heterogeneous materials: Microstructure and macroscopic properties. Springer-Verlag; 2002. ˇ [8] Sejnoha M, Zeman J. Overall viscoelastic response of random fibrous composites with statistically quasi uniform distribution of reinforcements. Computer Methods in Applied Mechanics and Engineering 2002;191(44): 5027–44. ˇ [9] Sejnoha M, Zeman J. Micromechanical analysis of random composites. CTU Reports, vol. 6. Prague: Czech Technical University; 2002. ˇ [10] Proch´azka P, Sejnoha J. A BEM formulation for homogenization of composites with randomly distributed fibers. Engineering Analysis with Boundary Elements 2003;27(2):137–44. ˇ [11] Zeman J, Sejnoha M. Homogenization of balanced plain weave composites with imperfect microstructure: Part I–theoretical formulation. International Journal for Solids and Structures 2004;41(22–23):6549–71. ˇ [12] Sejnoha M, Zeman J, Nov´ak J. Homogenization of random masonry structures—comparison of numerical methods. In: EM 2004–17th ASCE engineering mechanics division conference. Newark: University of Delaware; 2004. p. 1–8 [on CD ROM]. [13] Drugan W, Willis J. A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites. Journal of the Mechanics and Physics of Solids 1996;44(4): 497–524. ˇ [14] Sejnoha M, Valenta R, Zeman J. Nonlinear viscoelastic analysis of statistically homogeneous random composites. International Journal for Multiscale Computational Engineering 2004;2(4):645–74. [15] Huet C. Application of variational concepts to size effects in elastic heterogeneous bodies. Journal of the Mechanics and Physics of Solids 1990;38(6):813–41. [16] Sab K. On the homogenization and the simulation of random materials. European Journal of Mechanics A-Solids 1992;11(5):585–607. [17] Ostoja-Starzewski M. Material spatial randomness: From statistical to representative volume element. Probabilistic Engineering Mechanics (in press). doi:10.1016/j.probengmech.2005.07.007. [18] Pecullan S, Gibiansky LV, Torquato S. Scale effects on the elastic behavior of periodic and hierarchical two-dimensional composites. Journal of the Mechanics and Physics of Solids 1999;47(7):1509–42. [19] Terada K, Hori M, Kyoya T, Kikuchi N. Simulation of the multi-scale convergence in computational homogenization approaches. International Journal of Solids and Structures 2000;37(16):2285–311. [20] Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D. Determination of the size of the representative volume element for random composites: Statistical and numerical approach. International Journal of Solids and Structures 2003;40(13–14):3647–79. [21] Cluni F, Gusella V. Homogenization of non-periodic masonry structures. International Journal of Solids and Structures 2004;41(7):1911–23. [22] Sab K, Nedjar B. Periodization of random media and representative volume element size for linear composites. Comptes Rendus Mecanique 2005;333(2):187–95. [23] Gusella V, Cluni F. Random field and homogenization for masonry with non-periodic microstructure [submitted for publication]. [24] Pyrz R. Correlation of microstructure variability and local stress field in two-phase materials. Materials Science and Engineering 1994;A177: 253–9. [25] Lu B, Torquato S. Lineal-path function for random heterogeneous materials. Physical Review E 1992;45(2):922–9.
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329 [26] Zeman J. Analysis of composite materials with random microstructure. CTU reports, vol. 7. Prague: CTU; 2003. p. 177. [27] Beran MJ. Statistical continuum theories. Monographs in statistical physics. Interscience Publishers; 1968. [28] Stoyan D, Kendall W, Mecke W. Stochastic geometry and its applications. Berlin: Akademie-Verlag; 1987. [29] Derr R. Statistical modeling of microstructure with applications to effective property computation in materials science, Ph.D. thesis. University of North Carolina. Also available as http://citeseer.nj.nec.com/derr99statistical.html; 1999. [30] Markov KZ. On the cluster bounds for the effective properties of microcraked solids. Journal of the Mechanics and Physics of Solids 1998; 46(2):357–88. [31] Smith P, Torquato S. Computer simulation results for the two-point probability function of composite media. Journal of Computational Physics 1988;76:176–91. [32] Burrus C, Parks TW. DFT/FFT and convolution algorithms: Theory and implementation. In: Topics in digital signal processing. A WileyInterscience Publication; 1985. [33] Berryman JG. Measurement of spatial correlation functions using
[34]
[35] [36]
[37]
[38]
[39]
329
image processing techniques. Journal of Applied Physics 1984;57(7): 2374–84. Frigo M, Johnson S. FFTW: An adaptive software architecture for the FFT. Proceedings of the 1998 IEEE international conference on acoustics, speech and signal processing, vol. 3. New York: IEEE; 1998. p. 1381–4. See also http://www.fftw.org. Ohser J, M¨ucklich F. Statistics in practice. Chichester (England): John Wiley & Sons; 2000. Gajdoˇs´ık J. Qualitative analysis of fiber composite microstructure, Master’s thesis. Czech Technical University in Prague, Faculty of Civil Engineering; 2004. ˇ Matouˇs K, Lepˇs M, Zeman J, Sejnoha M. Applying genetic algorithms to selected topics commonly encountered in engineering practice. Computer Methods in Applied Mechanics and Engineering 2000;190(13–14): 1629–50. ˇ Zeman J, Sejnoha M. Numerical evaluation of effective properties of graphite fiber tow impregnated by polymer matrix. Journal of the Mechanics and Physics of Solids 2001;49(1):69–90. Laboratory Imaging L. LUCIA—Laboratory Universal Computer Image Analysis. http://www.laboratory-imaging.com/index.php.
Qualitative analysis of fiber composite microstructure: Influence of boundary conditions ˇ Jan Gajdoˇs´ık, Jan Zeman ∗ , Michal Sejnoha Department of Structural Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Th´akurova 7, 166 29 Prague 6, Czech Republic Received 16 August 2005; received in revised form 11 November 2005; accepted 21 November 2005 Available online 9 January 2006
Abstract The principal objective of this work is to classify the influence of different boundary conditions that can be applied when computing various statistical descriptors for generally random microstructures. Although applicable to most statistical descriptors, attention is limited to twopoint probability functions. Binary images of random fibrous composites are assumed in the present study. Here, the periodic, mirror and no (plain) boundary conditions are addressed when processing the binary images. Also, the minimum number of samples (Representative Volume Elements—RVEs) in the ensemble that still provide statistically representative results is sought. While the available results promote the use of periodic boundary conditions, particularly from the computational point of view regardless of the RVE size, the number of samples required is largely affected by their size. If available, the most efficient procedure results from the use of a large binary image of a real microstructure combined with periodic boundary conditions. c 2005 Elsevier Ltd. All rights reserved.
Keywords: Fibrous composites; Quantification of microstructure morphology; Two-point probability function; Boundary conditions; Representative volume element
1. Introduction to problem of microstructure analysis The reliable and accurate analysis of heterogeneous material systems increasingly often takes advantage of knowledge of real microstructures or even nanostructures; see, e.g., [1–3]. A specific example is given in Fig. 1, illustrating a random distribution of graphite fibers within a single bundle of large woven composite tube [4]. The periodic character of such systems, though often assumed, is generally not available. The random nature of most material systems at higher levels of resolution is therefore naturally accepted. Such a step, however, calls for a suitable method of microstructure1 quantification. A number of contributions have advocated the use of various statistical descriptors such as the one- and two-point probability functions or the lineal path function (see [5–8], to cite a few).
∗ Corresponding author. Tel.: +420 2 2435 4482; fax: +420 2 2431 0775.
E-mail addresses: [email protected] (J. Gajdoˇs´ık), ˇ [email protected] (J. Zeman), [email protected] (M. Sejnoha). 1 For simplicity, all structural levels, starting from the nanoscale up to the meso-scale, that resemble some kind of randomness, including random imperfections, will be jointly referred to as microstructure. c 2005 Elsevier Ltd. All rights reserved. 0266-8920/$ - see front matter doi:10.1016/j.probengmech.2005.11.006
In their recent works, the authors exercised the use of loworder statistical descriptors in a number of problems linked to the analysis of material systems with random or imperfect geometries. The formulation of statistically equivalent periodic unit cells, derived by matching the material statistics up to two-point probability function, has been offered not only for classical fibrous composites [9,10] but also for complex threedimensional textile composites with various sources of random imperfections [11], and even for masonry structures with an irregular arrangement of stone blocks [12]. Direct introduction of the two-point probability function into variational principles of Hashin and Shtrikman allowed a simplified yet reasonably accurate and extremely efficient analysis of complex material systems, not only in elastic [13] but also in visco- and nonlinear visco-elastic regimes [8,14]. In all applications, the periodic boundary conditions were intuitively selected when processing the binary images of various sizes taken as random cuts from images of true microstructures; Fig. 1(a). The resulting anisotropic response of such systems gave rise to various objections that question the assumed approach. Clearly, a number of errors may occur in the process of transforming the original color image into its binary
318
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Fig. 1. (a) A fiber tow, (b) a real micrograph of a selected cut.
counterpart. While relatively high resolution color images are available through the application of various image input devices, including the direct use of a digital camera or a camera connected to a microscope, the quality of the image itself depends on the quality of the original sample. Quite often, such specimens are of relatively poor quality, requiring manual adjustments of the resulting image for further analysis. The most critical region quite often appears in the vicinity of the sample boundary where the images are usually incomplete, e.g. “cut” fibers in the case of fibrous composites. This is particularly important when studying a large ensemble of small samples with imperfect boundaries created when taking random cuts. Investigating the image within this region leads to a number of problems. Note that the two-point probability function is defined as the probability of finding two specified points randomly thrown into a medium in specified phases. The finite dimensions of the image pose a clear question of how to treat cases, when one of the points is found out of the image boundary. Such a pair can be either directly disposed or the distribution of microstructure outside the image can be predicted in some way. Comparison of various techniques dealing with image boundaries supplies the major part of this contribution. In addition, the effect of the image size on the analyzed statistical descriptor is also studied. This problem is addressed particularly from the efficiency point of view, comparing the analysis on a single large sample (with
reference to ergodicity assumption) or a large ensemble of small samples. The minimum number of required samples given their size is provided. An additional constraint on the minimum size of RVE linked to the issue of statistical independence applied in [8] can also be put forward. This condition, however, was not examined in the present contribution. The present study is closely linked to the analysis of apparent linear mechanical properties of random composite materials initiated by Huet [15]. In this approach, to compensate for a limited size of an analyzed material sample, the microstructure is subjected to kinematic, static and periodic boundary conditions. With increasing size of the material sample, the response corresponding to these boundary conditions converges to a common value determining the size of the RVE; see also [16] for a detailed discussion and [17] for a recent exhaustive review. This procedure has been performed numerically for different material systems (see, e.g., [18–21] to cite a few) and all these studies confirm the superiority of periodic boundary conditions, even for non-periodic material systems. This conjecture lead to the recent introduction of a concept of periodization of random media [22]. This paper can be seen as an extension of these studies to quantification of the morphology of random microstructures. Finally, note that recent works on the definition of RVE for heterogeneous materials explicitly involve spatial statistics (up to two-point probability functions) for the definition of RVE [22,23]. For a moderate contrast between constituents’ properties, it seems that the condition of convergence in terms of a two-point probability function is more restrictive than the convergence of apparent material properties. The paper is organized as follows. The introductory part, Section 1, is followed by the specification of various statistical descriptors generally used in microstructure quantification. Several possible approaches for their evaluation, based on specific boundary conditions, are offered. Image analysis of a section of a graphite fiber tow is briefly discussed in Section 3. The results comparing various boundary conditions also from the RVE-size stand point are summarized in Section 4. 2. Microstructure quantification—statistical measures A number of statistical measures have been examined in the last few decades for addressing the problem of heterogeneous media with either a random distribution of material phases or the presence of various types of material or geometrical imperfections of more or less random nature. Starting from the characteristic function as the basic descriptor for microstructure quantification, the literature offers a number of derived statistical functions, including the radial path function [24], the lineal path function [25] or the general n-point probability function [6]. Of them all, the two-point probability function appears to be the most prominent, owing to the ease and efficiency of its evaluation, particularly when taking advantage of the Fast Fourier transform (FFT) applied to the binary images of real microstructures [26]. To prevent the reader from searching a vast body of literature devoted to this function, we now briefly summarize its most important features and the
319
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
essential statistical assumptions allowing its relatively efficient evaluation even for non-periodic media.
Further simplification arises when assuming the medium to be statistically isotropic [26]. Then Sr s (x1 , x2 ) reduces to
2.1. Fundamental function and statistical moments
Sr s (x1 − x2 ) = Sr s (k x1 − x2 k).
Consider an ensemble of a two-phase random medium. To provide a general statistical description of such a system, it is useful to characterize each member of the ensemble by a stochastic function — characteristic function χr (x, α), which is equal to one when point x lies in the phase r of the sample α and equal to zero otherwise [27,6,28], 1, if x ∈ Dr (α), χr (x, α) = (1) 0, otherwise, where Dr (α) denotes the domain occupied by the r -th phase. Except where noted, composites consisting of a clearly distinguishable continuous matrix phase are considered. Therefore, r = m, f , is further assumed to take values m for the matrix phase, while symbol f is reserved for the second phase. For such a system the characteristic functions χm (x, α) and χ f (x, α) are related by χm (x, α) + χ f (x, α) = 1.
(2)
Following [27,28,6], we write the ensemble average of the product of characteristic functions Sr1 ,...,rn (x1 , . . . , xn ) = χr1 (x1 , α) · · · χrn (xn , α),
(3)
where function Sr 1 , . . . , rn , referred to as the general npoint probability, gives the probability of finding n points x1 , . . . , xn randomly thrown into a medium located in the phases r1 , . . . , rn . 2.2. Functions of the first and second order Hereafter, we limit our attention to functions of the order of one and two, since higher-order functions are quite difficult to determine in practice.2 Therefore, the description of a random medium will be provided by the one-point probability function Sr (x): Sr (x) = χr (x, α),
(4)
which simply gives the probability of finding the phase r at x, and by the two-point probability function Sr s (x1 , x2 ) Sr s (x1 , x2 ) = χr (x1 , α)χs (x2 , α),
(5)
which denotes the probability of finding simultaneously the phase r at x1 and the phase s at x2 . In general, evaluation of these characteristics may prove to be prohibitively difficult. Fortunately, a substantial simplification applies when accepting an assumption regarding the material as being statistically homogeneous, so that Sr (x) = Sr , Sr s (x1 , x2 ) = Sr s (x1 − x2 ).
(6) (7)
2 Note, however, that relatively efficient procedures for the approximation of higher-order probability functions for ergodic and statistically isotropic media were proposed in [29].
(8)
Finally, making an ergodic assumption allows a substitution of the one-point correlation function by its volume average [26], i.e., volume concentration or volume fraction of the r -th phase cr : Sr = cr .
(9)
2.3. Limiting values In addition, the two-point probability function Sr s incorporates the one-point probability function Sr for certain values of its arguments such that for x1 = x2 : Sr s (x1 , x2 ) = δr s Sr (x1 ), for kx1 − x2 k → ∞ : lim Sr s (x1 , x2 )
(10)
kx1 −x2 k→∞
= Sr (x1 )Ss (x2 ),
(11)
where symbol δr s stands for Kronecker’s delta. Relation (10) states that the probability of finding two different phases at a single point is equal to 0 (see also Eq. (2)) or is given by the one-point probability function if phases are identical. Eq. (11) manifests that, for large distances, points x1 and x2 are statistically independent. This relation is often denoted as the no-long-range-order hypothesis (see e.g. [30,5]). Finally, according to Eq. (2), we may determine one- and two-point probability functions for all phases, provided that these functions are given for one arbitrary phase. For the onepoint probability function of a statistically homogeneous and ergodic medium, this relation assumes a trivial form: cm = 1 − c f .
(12)
Relations for the two-point probability functions of statistically uniform and ergodic medium are summarized in Table 1.3 2.4. Numerical evaluation To determine Sr1 ,...,rn , we recall that the general n-point probability gives the probability of finding n points x1 , . . . , xn randomly thrown into a medium located in the phases r1 , . . . , rn . In view of Table 1, we further consider only the fiber probability functions. To follow the above definition, the one-point fiber probability function S f gives the chance of finding a randomly placed point located in the fiber phase. To determine this quantity, a simple Monte Carlo-like simulation can be utilized — we randomly throw a point into the microstructure and count successful “hits” into the matrix phase. Then, the value of function S f can be estimated as Sf ≈
n0 , n
(13)
3 Note that, by definition (5) and the assumption of statistical homogeneity, Sr s (x) = Ssr (x).
320
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Table 1 Relations between two-point probability functions Known function
Smm (x)
Sm f (x)
S f f (x)
Smm (x) Sm f (x) S f f (x)
Smm (x) cm − Smm (x) c f − cm + Smm (x)
cm − Sm f (x) Sm f (x) c f − cm − Sm f (x)
cm − c f + S f f (x) c f − S f f (x) S f f (x)
Fig. 2. (a) Microstructure with periodic boundary conditions, (b) microstructure with mirror boundary conditions, and (c) microstructure with no (or “plain”) boundary conditions.
Fig. 3. (a) Original binary image with no modifications. (b) Modified binary image. (c) Regenerated “ideal” binary image.
where n 0 is the number of successful hits and n denotes the total number of throws. An entirely similar procedure can be employed to determine the values of S f f (x).4 Another more attractive approach is available when the real microstructure is replaced by its binary image. A sample of a binary image is shown in Fig. 3(c). Such a digitized micrograph can be imagined as a discretization of the characteristic function χr (x), usually presented in terms of a W × H bitmap. Denoting the value of χr for the pixel located in the i-th row and jth column as χr (i, j) allows the writing of the first moments of function χr for an ergodic and statistically homogeneous
4 For a statistically isotropic microstructure, Smith and Torquato [31] proposed a more efficient procedure for the determination of Smm (k x12 k). Instead of tossing a line corresponding to x into a medium, a sampling template is used for the determination of the two-point probability function. See also [9] for a comparison of this method with bitmap-based approaches presented hereafter and Section 4.2 for further application of this idea.
medium in the form5 Sr =
−1 H −1 X 1 WX χr (i, j). W H i=0 j=0
Formally, one can write the two-point probability function Sr s using the similar relation Sr s (m, n) =
−1 −1 H X 1 WX χr (i, j)χs (i + m, j + n), W H i=0 j=0
(14)
which, however, can generate indexes in the characteristic function χs exceeding the bitmap dimensions. To treat such a case, the boundary conditions depicted in Fig. 2 can be applied.
5 Throughout the text, the C-language type of array indexing is consistently used, i.e., we denote the first element of an array α as α0 and the last element of the array as α L−1 , where L is the array length.
321
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
2.4.1. Periodic boundary conditions In the case of periodic boundary conditions, Fig. 2(a), the starting image of a microstructure, e.g. the one displayed in Fig. 3(c), is assumed to be periodically translated in all directions. The two-point probability function then assumes the form Sr s (m, n) =
−1 H −1 X 1 WX χr (i, j)χs ((i + m)%W, W H i=0 j=0
( j + n)%H ),
Table 2 Bitmap dimensions in pixels Height
Width
Size
No. of bitmaps
1000 850 700 550 400 250
1000 850 700 550 400 250
1 000 000 722 500 490 000 302 500 160 000 62 500
18 × 1 18 × 10 18 × 20 18 × 30 18 × 40 18 × 50
(15)
where the symbol “%” stands for modulo. The above equation, usually termed the cyclic correlation [32], readily implies periodicity of function Sr s . Note that the correlation property of DFT holds for cyclic correlation. Referring to [32], it is given by the following relation DFT{Sr s (m, n)} = DFT{χr (m, n)}DFT{χs (m, n)},
(16)
where · now denotes the complex conjugate. The inverse DFT, denoted IDFT, then serves to derive function Sr s at the final set of discrete points as [33] 1 IDFT{DFT{χr (m, n)}DFT{χs (m, n)}}. (17) WH This method is very economical and its accuracy depends only on the selected resolution of the digitized medium. Usually, the Fast Fourier Transform, which needs only O(W H log(W H ) + W H ) operations, is called on to carry out the numerical computation.6 This approach is extremely efficient, especially when comparing with traditional “Monte Carlo”-like or “trial and error” methods, which need to be employed for the remaining boundary conditions. Sr s (m, n) =
2.4.2. Mirror boundary conditions In this particular case, the microstructure is expected to be mirrored; Fig. 2(b). When taking this route we are bound to modify Eq. (14) in the form Sr s (m, n) =
−1 H −1 X 1 WX χr (i, j)χs (i\ + m, j[ + n), W H i=0 j=0
(18)
where i +m for 0 ≤ (i + m) < W i\ +m = 2W − 1 − (i + m) for W ≤ (i + m) < 2W j +n for 0 ≤ ( j + n) < H j[ +n = 2H − 1 − ( j + n) for H ≤ ( j + n) < 2H. Observe that to O((W H )2 ) operations are needed to evaluate the function Sr s . This might be computationally demanding, particularly for a large micrograph. These boundary conditions are therefore selected only for the purpose of comparison. Nevertheless, there are still other fields of research where this approach may prove to be useful. 6 The public-domain package FFTW version 2.1.3 [34] was used for the evaluation of Eq. (17).
Fig. 4. Graphical representation of the process of bitmap selection.
2.4.3. Plain boundary conditions In the third case, here termed the plain boundary condition, we essentially disregard the cases for which χs is located outside the image; see Fig. 2(c). The particular form of the twopoint probability function follows from Sr s (m, n) =
(i M
iX N −1 M −1 jX 1 χr (i, j) − i m )( j N − jn ) i=i j= j
χs (i + m, j + n),
m
n
(19)
where i m = max(0, −m), i M = min(W, W − m) and jn = max(0, −n), j N = min(H, H − n). Although the computational complexity of this approach is comparable to the mirror boundary condition, it is the most “suitable” one, as it introduces not a prior assumption about the sample boundary7 but significantly reduces the amount of information obtained. It is evident that long vectors (m, n) can be positioned just a few times to be totally contained in the image. If the vector with the length of the rectangle’s diagonal is used, it can be thrown just once. It becomes clear that results for large spacing of the two points yield lower reliability when compared with those with small spacing. On the other hand, it has been demonstrated in, see e.g. [8,9], that, for distances exceeding approximately 10 times the fiber diameter, a random throw of the two points becomes a statistically independent event and therefore the results for large spacing can be safely disregarded. 7 More precisely, it can be shown that these boundary conditions deliver an unbiased estimator of two-point probability function; see, e.g., [35, Section 5.1.3].
322
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Fig. 5. Exemplary graph of measures M.
3. Image analysis of fibrous composites Image processing represents an important part of the micromechanical analysis of real world material systems allowing scientists, through computer based image acquisition, to visualize essential material data. Unfortunately, the original data, usually supplied in the form of color images, are often of poor quality, so that further treatment is needed. Such a tedious step typically involves manual marking of object borders, completion of distorted objects, extracting individual color spaces, filtration and final conversion into binary image. Although limiting the subsequent analysis to two-phase material systems, the two-dimensional binary data provide detailed information about microstructure configuration and may reveal various sources of imperfections for relatively complex and even three-dimensional structures such as textiles [11]. Note that evaluation of the basic statistical descriptors described in Section 2 is not limited to two-phase systems. At the present time the color images, transformed into 256 colors in gray scale, can be analyzed by applying the same tools as for the binary image data [36]. Nevertheless, owing to the complexity of the process of retrieving the images of a real composite, we settle in the present study for conventional binary images. Such a simulation of real microstructure is more than sufficient when studying the effect of boundary conditions on the evaluation of material statistics from Section 2. In view of our previous work, we continue with the analysis of graphite fiber tow impregnated by a polymer matrix [37,38,8,14]. All specimens subjected to individual steps of image processing were extracted from a single-layer woven composite tube supplied by Compo-Tech, Ltd. [4]. To adjust the original color image of Fig. 1, the image analyzer LUCIA G [39] was employed to obtain the desired binary (monochromatic) image data. The result of the first step appears in Fig. 3(a). The respective binary image was obtained through the process of gradual thresholding of the original color image by setting the range of color scales or intensity and saturation for objects, for which the output binary image should be black (or white). Clearly, such an image is not very suitable for further analysis as it displays, as a result of sample preparation, a number of partially damaged and
possibly connected fibers. Manual separation and completion of individual fibers applicable in the software LUCIA G resulted in an improved image, plotted in Fig. 3(b). The subsequent analysis of this image provided the basic information about the original microstructure, such as the number of objects and positions of their gravity centers. Although this image is relatively satisfactory and applicable for the evaluation of the desired statistical descriptors, it still shows fibers of various sizes and shapes that are far from circular. The “ideal” image in Fig. 3(c) was thus created by generating circular fibers placed in the gravity centers of original fibers. The fiber radius common for all fibers was selected such as to keep the original fiber volume fraction. Further details can be found in [36]. This final image was then exploited in Section 4 to compute the two-point probability function introduced in Section 2. The above approach was exercised to retrieve about 25 images from various sections of the fiber tow cross-section of Fig. 1(a). Each image of approximately 200 fibers is assumed to be representative of the actual material system, at least when accepting the ergodicity assumption tested, e.g. in [9]. Eventually, only 18 images out of all 25 were selected for the analysis. In particular, images with at least one fifth of the area being fiber free were eliminated. Those images clearly experienced some error, which was not captured when processing the image. Finally, all remaining images were adjusted by taking a smaller cut-out to reintroduce partial fibers along the “ideal” image boundary, Fig. 3(c), which were originally present in the image, Fig. 1(b) and Fig. 3(a), but removed when applying the modification steps. 4. Two-point probability function for various boundary conditions As stated in Section 2, probability functions of infinite order are needed to fully describe the material statistics of a random medium. This step, however, is computationally infeasible. Note that even the three-point probability function that appears in some variational formulations [27,30] requires enormous computational time for its evaluation. Consequently, only the two-point probability function will be considered hereafter, as it is still relatively easy to evaluate, especially when assuming a generally random medium. Recall that, by definition, the two-point probability function can be interpreted as throwing a needle of a certain length into the image and counting the number of successful hits, out of all of them, when both ends of the needle fall into the selected colors (the number of colors is equal to the number of assumed phases). It remains to determine what happens when one end of the needle falls out of the image. This event will be studied here, assuming three possible boundary conditions introduced in Section 2.4. 4.1. Two-point probability function via periodic boundary conditions The basic steps for the evaluation of the two-point probability function for binary and periodic media were
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
323
Fig. 6. Graphs of measures M derived for periodic boundary condition.
described in Section 2. The most straightforward approach is to perform an averaging of the two-point probability function for all 18 statistically similar samples forming the ensemble. It can be advantageous, however, to check whether evaluating the two-point probability function on a larger set of smaller samples could be more efficient. It remains to determine how many samples of smaller size, taken as random cuts from the original “large” samples, are then required to achieve the desired accuracy provided by the large samples, if at all possible. To that end, the following strategy was followed. First, the largest sample of size 1000 × 1000 pixels was examined. In the next step the size was reduced by 150 pixels in both directions to arrive at a bitmap of 850 × 850 pixels. This step appeared reasonable, particularly with reference to the average fiber size of about 70 pixels. The smaller image was selected 10 times from random positions of the largest image. In each subsequent step the size was reduced again by 150 pixels in both directions and the random selection was increased by 10. The resulting bitmaps subjected to further study are listed in Table 2. Recall that this process was carried out 18 times for each of the largest bitmap. A graphical representation of this process appears in Fig. 4. The size of the pools was selected from the point of view of computational feasibility, especially when referring to Eqs. (18) and (19) for the two-point probability function evaluation. To compare individual results, the following measure was adopted:
M=
1 b·H b W
·
b−1 b −1 H W X X i=0 j=0
|Sr s (i, j) − Pr s (i, j)| . |Sr s (i, j)|
(20)
b and H b represent the number of points used in In Eq. (20), W the width and height directions of the bitmap, respectively. If every point of the bitmap is taken into account for the measure b and H b correspond exactly to the width and evaluation, then W height of a given bitmap. Sr s (i, j) is the value of the twopoint probability function of an image evaluated at the point with coordinates i and j. Pr s (i, j) is the value of the twopoint probability function of the largest bitmap that is assumed to provide the “exact” solution. This enabled comparison of measures taken from bitmaps with variable sizes. The results of the evaluation of the two-point probability function were stored in the matrix of the same dimensions as the bitmap. Note that function Sr s (i, j) corresponds to an average value evaluated for a given number of bitmaps selected from each pool. In other words, not only the total set of bitmaps (e.g., 20 samples for the bitmap size equal to 700 × 700 pixels) for each pool was used, but also partial sets, selected sequentially to get smaller pools of 2, 3, 4, . . . N (with N being the number of bitmaps of a given size), were tested. The random selection of bitmaps from each pool for averaging purposes was omitted, as the pool itself was already generated randomly. The specific average value of Sr s (i, j) was obtained as an average of independently computed two-point probability functions at the position i j.
324
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
Fig. 7. Evaluation of function S f f from Eq. (19): (a) selection of a size and vector orientation; (b) positioning of a vector into a bitmap; (c) effect of vector positioning; and (d) checking isotropic assumption.
Fig. 8. Graphs of measures M derived for mirror boundary condition.
The results appear in Fig. 6. For detailed explanation of the results presented, we refer the reader to Fig. 5. The selected
title for each graph corresponds to the results for the smallest bitmap size plotted in a given graph. Each line is constructed
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
325
Fig. 9. Graphs of measures M derived for plain boundary condition.
Fig. 10. Comparison of measure M for various boundary conditions.
from averages of the measure M computed for smaller sets of bitmaps available for a given pool. Therefore, for the first graph
of bitmap size 850 × 850 in Fig. 6, the first point represents a value of M derived from a single bitmap (more exactly an
326
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
average of 18 bitmaps, as there are 18 large bitmaps at hand to create smaller pools) while the last point stands for the value of measure M found as an average from 10 × 18 bitmaps (the largest number of bitmaps available for this bitmap size; see Table 2). Each of the remaining graphs thus stores the results provided by all bitmaps of larger sizes. It is evident that larger bitmaps enable results to be obtained with better reliability. Quite surprisingly, even for smaller bitmaps there is no common value where all measures completely stabilize. Also note that averages of the measure M for various bitmap sizes oscillate around different values. Nevertheless, all the graphs clearly show that the smaller the bitmap size, the larger the number of bitmaps needed to obtain acceptable results. For most of the assumed bitmap sizes, the number of 20 bitmaps appears sufficient, at least from the selected measure point of view. 4.2. Two-point probability function via mirror and plain boundary conditions Recall that, for the mirror and plain boundary conditions, the lack of periodicity prevents the use of FFT techniques for the evaluation of the two-point probability function, which leaves us with Eq. (19) to deal with. It has been suggested in [9,26] that the medium under consideration can be regarded as ergodic and statistically homogeneous. In such a case, the two-point probability function can be evaluated irrespective of the choice of the coordinate system origin. As such, it depends only on the absolute distance of the two points and the orientation of the line connecting the two points. In view of this, the polar coordinate system r, φ is typically adopted when plotting the function. In addition, due to the symmetry of the function, the angle φ sweeping the area from 0◦ to 180◦ is sufficient, which considerably reduces the computational time. Nevertheless, when compared to the solution of Eq. (17) that employs the FFT method, the solution of Eq. (19) for a statistically homogeneous medium is still extremely demanding computationally. Owing to the fact that our primary goal is not to analyze the current medium exactly, we settled, in view of the above comments, for the assumption of statistical isotropy, so that the solution no longer depends on the angle orientation. Such an assumption renders a relatively fast and simple algorithm for evaluation of the function. Individual steps can be inferred from Fig. 7. As shown in Fig. 7(a), the same grid of points corresponding to the resolution of the analyzed bitmap was created. The two-point probability function was then evaluated for a set of vectors connecting the coordinate system origin placed in the upper left corner and each point in the grid. The starting point of each vector was then successively placed into every point in the bitmap, as depicted in Fig. 7(b). The location for the end point was then checked, and a number of successful hits that met the required condition were recorded, e.g. for the fiber–fiber two-point probability function S f f both points are required to be found in the fiber phase (see Fig. 7(c); clearly, among all four vectors, only vector No. 4 is acceptable — the remaining three vectors must be disregarded). Dividing this value by the total number of hits (the number of positionings of the vector of a given length into the bitmap) gives the desired
Fig. 11. Distribution of function S f f : (a) derived from a single bitmap, (b) averaged over the maximum number of sample sets for a given bitmap.
value of the two-point probability function for a given vector length. See also Fig. 7(c) for a graphical explanation of the condition of statistical isotropy (vectors a, b and c are identical, and the value of S f f does not depend on their orientation; if only the condition of statistical homogeneity applies, then the two vectors a and b are identical but differ from vector c for the S f f function evaluation). Additional simplifying assumptions can be introduced to improve the efficiency of function evaluation. Two possibilities are at hand. First, the value of the two-point probability function can be computed for all available vectors of the function generated in Fig. 7(a), but not all the pixels in the bitmap are used for the function evaluation, Fig. 7(b). In other words, the needles of every available length are used, but just a limited number of throws is done. Another possibility is to evaluate the function just for a few specific vectors, but all data stored in the bitmap are used for the function evaluation. In other words, we throw the needles of just a few chosen lengths, but the needles are thrown into every point of the bitmap. In the present work the second possibility was adopted, as it allows direct comparison with results obtained by the Fast Fourier Transform method that evaluates the function using all data stored in the image. The selection of step size was based on the average size of a single fiber equal to approximately 70 pixels. This value promoted use of the step of 15 pixels, which is less than one quarter of the fiber diameter. The obtained results appear in Figs. 7(d)–9. Fig. 7(d) compares the variation of the two-point probability function
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
327
Fig. 12. Plane view of function S f f for various boundary conditions: (a) periodic boundary conditions, (b) mirror boundary conditions, (c) plain boundary conditions.
for mirror and plain boundary conditions derived from Eq. (19), assuming statistical isotropy with the results derived for the periodic boundary conditions under the assumption of statistical homogeneity and then averaged over all angles for a given distance r to obtain the necessary isotropized values. Evidently, the results are almost identical, particularly when comparing the periodic and plain boundary conditions. Figs. 8 and 9 display plots of the measure M derived for the same sample sizes as presented in Section 4.1. Similar concluding remarks to those already suggested in the previous section can be drawn from these graphs. 4.3. Comparison of individual boundary conditions Except for Fig. 7(d), the results discussed so far were presented with no reference to each other. To get better insight into the effect of the selected boundary conditions, the results from Figs. 6, 8 and 9 are re-plotted in Fig. 10 for the selected sample sizes. On each graph, only the results corresponding to the selected sample size are shown. Each line in the graph displays the measure M found for a specific boundary condition. The results are rather inconclusive, and it appears that the same number of samples of smaller bitmaps are needed to stabilize the solution. A valuable piece of information is also available from Fig. 11, showing the variation of the fiber–fiber two-point probability function derived from selected sample sets. The solid line represents the results derived from one of the largest 1000 × 1000 bitmaps selected as an average from the original ensemble of 18 bitmaps. The dashed and dotted lines were then found from sample sets of 20 (700 × 700) and 50 (400 × 400) bitmaps, respectively. Evidently, the smaller the sample size the shorter the corresponding line, as longer vectors do not fit into smaller cut-outs. As can be already deduced from the previous results, if samples of smaller sizes are used, it is strongly recommended to perform averaging over the number of samples in the ensemble to arrive at meaningful results; compare Figs. 11(a) and (b). The final comparison is presented in the form of twodimensional plots of the fiber–fiber two-point probability function. It follows from Figs. 12 and 13 that the difference between periodic and mirror or plain boundary conditions is random and does not follow any regular pattern. The assumed
Fig. 13. Difference in function S f f plotted in absolute values: (a) comparison of periodic and mirror boundary conditions, (b) comparison of periodic and plain boundary conditions.
isotropy is also evident. Clearly, the visible difference is shifted away from the origin where statistical dependency no longer applies. Close to the origin, where the distribution of the two-point probability function provides the most significant statistical information about the microstructure configuration, the difference is almost negligible. These final plots thus further support the fact that neither boundary condition is superior to the other, apart from the statement that the plain boundary conditions are capable of delivering statistically unbiased results. 5. Conclusion It can be concluded from the results presented that all boundary conditions can be considered as equal. Such a
328
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329
conclusion follows not only from Fig. 10, showing a close match of measure M for various boundary conditions, but also from Fig. 13, suggesting negligible difference in the twopoint probability function, especially between periodic and plain boundary conditions. Recall also Fig. 7(d), advocating no difference between the periodic and plain boundary conditions. The similarity in the behavior of the selected measure M can also be inferred from Figs. 6 and 8–10, which further supports the suggested conclusion. Due to the fact that no boundary conditions are superior to any other, the only objective that remains is the computational complexity. While the use of periodic boundary conditions allows implementation of a very efficient algorithm based on the Fast Fourier Transform (FFT) method resulting in negligible computational time, the solution of Eq. (19) for mirror and plain boundary conditions is much more demanding. For the largest bitmap of 1000×1000 pixels, the FFT method is about 125 000 faster on a standard PC with a 1.6 GHz Athlon 3000+ processor. The final conclusion is therefore twofold: • As the computational time needed for the evaluation of Eq. (19) depends heavily on the bitmap dimensions, it is recommended that, if the plain boundary conditions are needed for some reason, it is more efficient to compute the two-point probability function from a larger set of smaller bitmaps. With reference to Fig. 11 their size, however, should be sufficiently large to receive enough statistical information about microstructure configuration. • If, on the other hand, the applicability of periodic boundary conditions is possible, it is reasonable to use a single bitmap of as large a size as possible. Although the results for a single large bitmap and sets of smaller bitmaps are at least qualitatively comparable, see Fig. 11, the results for larger bitmaps are still more accurate. A rather naive but illustrative explanation of all the presented results assumes that the two-point probability function behaves a little bit like the surface of water in a pond after being hit by a stone in the middle of the pond. If stones of different shapes, but quite similar in size, are thrown into a pond, the waves close to the place of impact are quite different but further from the place of impact the waves are small. If the wave is small, the difference is small too. If the observer is monitoring just the small area of the surface close to the place of impact, he will see that the differences between the shapes of waves made by stones of different shapes are bigger than if he monitors a larger area. The results from this research are the same. The average values of measure M are always smaller if determined from functions computed for larger bitmaps; compare Figs. 10(a) and (e). Acknowledgments ˇ No. 106/03/H150 and The support provided by the GACR CEZ MSM 6840770003 is gratefully acknowledged. References [1] Babuˇska I, Andersson B, Smith PJ, Levin K. Damage analysis of fiber composites Part I: Statistical analysis on fiber scale. Computer Methods in Applied Mechanics and Engineering 1999;172(1–4):27–77.
[2] Torquato S. Statistical description of microstructures. Annual Review of Materials Research 2002;32:77–111. [3] Raghavan P, Ghosh S. Adaptive multi-scale computational modeling of composite materials. CMES-Computer Modeling in Engineering & Sciences 2004;5(2):151–70. [4] Compo Tech Plus webpage. http://www.compotech.com. [5] Willis JR. Bounds and self-consistent estimates for the overall properties of anisotropic composites. Journal of the Mechanics and Physics of Solids 1977;25:185–202. [6] Torquato S, Stell G. Microstructure of two-phase random media. I. The n-point probability functions. Journal of Chemical Physics 1982;77(4): 2071–7. [7] Torquato S. Random heterogeneous materials: Microstructure and macroscopic properties. Springer-Verlag; 2002. ˇ [8] Sejnoha M, Zeman J. Overall viscoelastic response of random fibrous composites with statistically quasi uniform distribution of reinforcements. Computer Methods in Applied Mechanics and Engineering 2002;191(44): 5027–44. ˇ [9] Sejnoha M, Zeman J. Micromechanical analysis of random composites. CTU Reports, vol. 6. Prague: Czech Technical University; 2002. ˇ [10] Proch´azka P, Sejnoha J. A BEM formulation for homogenization of composites with randomly distributed fibers. Engineering Analysis with Boundary Elements 2003;27(2):137–44. ˇ [11] Zeman J, Sejnoha M. Homogenization of balanced plain weave composites with imperfect microstructure: Part I–theoretical formulation. International Journal for Solids and Structures 2004;41(22–23):6549–71. ˇ [12] Sejnoha M, Zeman J, Nov´ak J. Homogenization of random masonry structures—comparison of numerical methods. In: EM 2004–17th ASCE engineering mechanics division conference. Newark: University of Delaware; 2004. p. 1–8 [on CD ROM]. [13] Drugan W, Willis J. A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites. Journal of the Mechanics and Physics of Solids 1996;44(4): 497–524. ˇ [14] Sejnoha M, Valenta R, Zeman J. Nonlinear viscoelastic analysis of statistically homogeneous random composites. International Journal for Multiscale Computational Engineering 2004;2(4):645–74. [15] Huet C. Application of variational concepts to size effects in elastic heterogeneous bodies. Journal of the Mechanics and Physics of Solids 1990;38(6):813–41. [16] Sab K. On the homogenization and the simulation of random materials. European Journal of Mechanics A-Solids 1992;11(5):585–607. [17] Ostoja-Starzewski M. Material spatial randomness: From statistical to representative volume element. Probabilistic Engineering Mechanics (in press). doi:10.1016/j.probengmech.2005.07.007. [18] Pecullan S, Gibiansky LV, Torquato S. Scale effects on the elastic behavior of periodic and hierarchical two-dimensional composites. Journal of the Mechanics and Physics of Solids 1999;47(7):1509–42. [19] Terada K, Hori M, Kyoya T, Kikuchi N. Simulation of the multi-scale convergence in computational homogenization approaches. International Journal of Solids and Structures 2000;37(16):2285–311. [20] Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D. Determination of the size of the representative volume element for random composites: Statistical and numerical approach. International Journal of Solids and Structures 2003;40(13–14):3647–79. [21] Cluni F, Gusella V. Homogenization of non-periodic masonry structures. International Journal of Solids and Structures 2004;41(7):1911–23. [22] Sab K, Nedjar B. Periodization of random media and representative volume element size for linear composites. Comptes Rendus Mecanique 2005;333(2):187–95. [23] Gusella V, Cluni F. Random field and homogenization for masonry with non-periodic microstructure [submitted for publication]. [24] Pyrz R. Correlation of microstructure variability and local stress field in two-phase materials. Materials Science and Engineering 1994;A177: 253–9. [25] Lu B, Torquato S. Lineal-path function for random heterogeneous materials. Physical Review E 1992;45(2):922–9.
J. Gajdoˇs´ık et al. / Probabilistic Engineering Mechanics 21 (2006) 317–329 [26] Zeman J. Analysis of composite materials with random microstructure. CTU reports, vol. 7. Prague: CTU; 2003. p. 177. [27] Beran MJ. Statistical continuum theories. Monographs in statistical physics. Interscience Publishers; 1968. [28] Stoyan D, Kendall W, Mecke W. Stochastic geometry and its applications. Berlin: Akademie-Verlag; 1987. [29] Derr R. Statistical modeling of microstructure with applications to effective property computation in materials science, Ph.D. thesis. University of North Carolina. Also available as http://citeseer.nj.nec.com/derr99statistical.html; 1999. [30] Markov KZ. On the cluster bounds for the effective properties of microcraked solids. Journal of the Mechanics and Physics of Solids 1998; 46(2):357–88. [31] Smith P, Torquato S. Computer simulation results for the two-point probability function of composite media. Journal of Computational Physics 1988;76:176–91. [32] Burrus C, Parks TW. DFT/FFT and convolution algorithms: Theory and implementation. In: Topics in digital signal processing. A WileyInterscience Publication; 1985. [33] Berryman JG. Measurement of spatial correlation functions using
[34]
[35] [36]
[37]
[38]
[39]
329
image processing techniques. Journal of Applied Physics 1984;57(7): 2374–84. Frigo M, Johnson S. FFTW: An adaptive software architecture for the FFT. Proceedings of the 1998 IEEE international conference on acoustics, speech and signal processing, vol. 3. New York: IEEE; 1998. p. 1381–4. See also http://www.fftw.org. Ohser J, M¨ucklich F. Statistics in practice. Chichester (England): John Wiley & Sons; 2000. Gajdoˇs´ık J. Qualitative analysis of fiber composite microstructure, Master’s thesis. Czech Technical University in Prague, Faculty of Civil Engineering; 2004. ˇ Matouˇs K, Lepˇs M, Zeman J, Sejnoha M. Applying genetic algorithms to selected topics commonly encountered in engineering practice. Computer Methods in Applied Mechanics and Engineering 2000;190(13–14): 1629–50. ˇ Zeman J, Sejnoha M. Numerical evaluation of effective properties of graphite fiber tow impregnated by polymer matrix. Journal of the Mechanics and Physics of Solids 2001;49(1):69–90. Laboratory Imaging L. LUCIA—Laboratory Universal Computer Image Analysis. http://www.laboratory-imaging.com/index.php.