Efficient 3D reconstruction of random heterogeneous media via random process theory and stochastic reconstruction procedure

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 354 (2019) 1–15 www.elsevier.com/locate/cma

Efficient 3D reconstruction of random heterogeneous media via random process theory and stochastic reconstruction procedure Wenliang Zhanga , Lei Songa ,∗, Juanjuan Lib a

State Key Laboratory of Geomechanics and Deep Underground Engineering, School of Mechanics & Civil Engineering, China University of Mining and Technology, Xuzhou 221116, China b IoT Perception Mine Research Center, China University of Mining and Technology, Xuzhou, 221116, China Received 27 December 2018; received in revised form 26 April 2019; accepted 22 May 2019 Available online 27 May 2019

Highlights • Three-dimensional reconstruction of random heterogeneous media. • An enhanced autocorrelation function to generate more detailed features. • An improved algorithm by integrating random process theory and SAA.

Abstract This paper introduces an algorithm for reconstructing three-dimensional (3D) heterogeneous media, in which the randomly distributed individual entities, such as aggregates, pores or other material phases with different sizes and shapes can be explicitly represented. The reconstruction procedure achieves high efficiency and reliability by integrating random process theory and stochastic reconstruction procedure. Specifically, the random process theory is applied to generate an initial model based on the statistics obtained from the digital image of heterogeneous material sample and then this coarse model is evolved by the simulated annealing algorithm to match the same statistical descriptors of sample structure. Two 3D models of different sizes are reconstructed. The autocorrelation function is applied to validate the accuracy of the reconstructed models. The comparison between the reconstructed models and the real model shows good agreement. The proposed method is suitable for reconstructing large size 3D random media for its high computational efficiency. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: 3D reconstruction; Random heterogeneous media; Random process theory; Stochastic reconstruction procedure; Simulated annealing

1. Introduction Random heterogeneous media are ubiquitous in practical engineering application. Evaluating or predicting the macroscopic properties of these materials is the primary task for material design or service status assessment and prediction. Macroscopic properties, such as strength, elastic modulus, permeability, conductivity and wave impedance, depend mainly on the media structures. Take concrete, which bears the functions of sealing, supporting, ∗ Corresponding author.

E-mail address: [email protected] (L. Song). https://doi.org/10.1016/j.cma.2019.05.033 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

bearing and waterproof in modern infrastructure, as an example, the size of aggregate has a certain impact on the attenuation [1] and velocity [2] of the ultrasonic wave, which are the key parameters of ultrasonic non-destructive testing; the shape of aggregate also has a significant effect on the microstructure of concrete by influencing the diffusion rate of mortar [3]. The permeability and strength of concrete are determined mainly by the type of aggregate [4] and mixing ratio [5]. Whether predicting the macroscopic properties of materials to be designed or estimating that of structures in site, the primary task is to obtain reliable samples with detailed structural information. Although Scanning Electron Microscopy and computed micro-tomography (CT) [6] are currently available to achieve this goal with providing a direct 3D description of the internal structure, they are usually timeconsuming and economically expensive. More importantly, both methods can only provide a single small structural model, which is clearly not enough for material design and obtaining statistical macroscopic properties. It is of great significance to reconstruct large-size and large number of media models based on the limited morphological information of a few or even a single scanning sample. To achieve the purpose, digital media model reconstruction can be an efficient and economical alternative. Various reconstruction methods [7] based on limited morphological and statistical information extracted from 2D/3D imaging data have been proposed for a wide range of random media. The reconstruction models include particulate structures [8,9], porous structures [10,11] and composite materials [12]. Such reconstruction methods include the Gaussian random field method [13], random-set method [14] and stochastic reconstruction procedure [15]. In the Gaussian random field method, the goal is to generate a random field to satisfy the two-point correlation requirement of the target field. In general, a field–field correlation function is used to constructed the Gaussian field and level-cut is applied to obtain a binary microstructure. The performance and applicability of this method depend on whether the field–field correlation function is appropriate, which is often difficult to meet [13]. In the random-set method, the random morphology of inclusions is constructed by Boolean operation [14]. This method can be fast for composites with sphere or polygon inclusions. However, for the inclusions with arbitrary shapes, the calculation is complex and difficult to realize. The main idea of stochastic reconstruction procedure is to quantitatively characterize the initial and reference model with appropriate statistical functions, and then iteratively optimize the initial model with simulated annealing (SA) [16] until the reconstructed model is statistically equivalent to the reference model. Combining with the dilation–erosion method we can improve the accuracy of the reconstruction [17]. This method can handle a wide range of heterogeneous materials because it can incorporate any type and number of correlation functions [18]. However, incorporating multiple functions will multiply the computational expenses of optimization iterative. A few works have demonstrated that computational costs can be significantly reduced from the aspects of effective initial model [19] and efficient updating procedure of statistical functions [20]. In this paper, an efficient reconstruction algorithm is proposed, which can be directly used to generate 3D model of two-phase random media. Firstly, an initial model is generated by random process theory and the reconstruction in this step exhibits high efficiency as the main computation is done by fast Fourier transform. The stochastic reconstruction procedure is applied to further refine the initial structure model. Because the initial structure model has preliminarily possessed the statistical characteristics of the reference model, the number of interchanges in the iterative optimization is greatly reduced and the convergence of reconstruction process is accelerated. The remainder of the paper is structured as follows: in Section 2, several statistical functions used for random heterogeneous media reconstruction are introduced. In Section 3, the algorithm flow of random media reconstruction is described in detail. In Section 4, two 3D random media of different sizes are reconstructed by the proposed algorithm with a 3D concrete structure as the reference sample. Final concluding remarks are summarized in Section 5. 2. Statistical and morphological functions of random media With the development of hardware and software, digital image processing technique has been an innovative way to accurately measure the spatial distribution and morphology of inhomogeneities in random media [21]. A digital image consists of a set of discrete pixels and each pixel has an integer value that is called gray level. A digital image can be expressed as a discrete function g [22]. For a Cartesian coordinate system, a 2D image of m × n can

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

3

Fig. 1. (a) Grayscale image of concrete; (b) variation of gray level along the line; (c) binary image.

be expressed as: ⎡

f (1, 1) ⎢ f (2, 1) ⎢ g(i, j) = ⎢ . ⎣ .. f (n, 1)

f (1, 2) f (2, 2) .. .

... ...

f (n, 2) ...

⎤ f (1, m) f (2, m)⎥ ⎥ .. ⎥ . ⎦

(1)

f (n, m)

The statistical characteristics used for random media reconstruction in this paper are derived from images obtained by CT scanning and digital camera. Fig. 1(a) is a 2D grayscale image of a concrete cross-section, it can be clearly seen that there is a significant difference in gray level between aggregate and mortar matrix: the aggregate appears brighter and the cement matrix is darker. The gray level is extracted along the blue line in Fig. 1(a) and draw in Fig. 1(b). Obviously, the area with gray level higher than a certain threshold belongs to aggregate, while the area with gray level lower than this threshold belongs to cement matrix. According to this difference, the separation of aggregate and cement matrix can be achieved by applying threshold segmentation and a digital model of two-phase random media can then be obtained as shown in Fig. 1(c). It can be concluded that the random heterogeneity of concrete is consistent with the spatial distribution of gray level. Therefore, the statistical and morphological characteristics of individual entities embedded in matrix can be obtained by statistical description of the image. 2.1. The characterization of random field The gray level of a random media image can be regarded as a random variable and the distribution of gray level in space can be modeled as a random field. Under the ergodic hypothesis, the statistical parameters regarding the random field, such as mean g, variance σ 2 , root mean square gr ms are all invariant when shifted in space. g, σ 2 and gr ms are given as follows: g=

m n 1 ∑∑ gi j m × n i=1 j=1 m n )2 1 ∑∑( gi j − g m × n i=1 j=1   m ∑ n  1 ∑ =√ g2 m × n i=1 j=1 i j

(2)

σ2 =

(3)

gr ms

(4)

Another important statistical function for characterizing random field is autocorrelation function G(x, y), which is an important parameter to describe the characteristics of stationary random media structure [23]. The autocorrelation

4

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

Fig. 2. Second-order statistics of the random media in Fig. 1(a): (a) normalized autocorrelation function; (b) power spectrum function.

function is defined as [24]: 1 L 1 →∞ 4L 1 L 2

G (x, y) = lim

L 2 →∞

∫

L1

−L 1

∫

L2

g (x1 , y1 ) g (x1 + x, y1 + y) dx1 dy1

(5)

−L 2

Following the isotropic assumption, the autocorrelation function only depends on the distance r = (x 2 + y 2 )1/2 between two points. Other statistics, such as the distribution of crests heights, the density of crests and the mean curvature of crests are all related to the power spectrum of the random field. According to the Wiener–Khinchin theorem, the power spectrum of a stationary random field is the Fourier transform of the corresponding autocorrelation function: ∫ ∞∫ ∞ ( ) [ ] 1 Φ k x , k y = F [G (x, y)] = G(x, y)exp −i(xk x + yk y ) dxdy (6) 2 4π −∞ −∞ Similarly, the autocorrelation function is the inverse Fourier transform of the power spectrum: ∫ ∞∫ ∞ ( ) [ ] [ ( )] −1 Φ k x , k y exp i(xk x + yk y ) dk x dk y G (x, y) = F Φ kx , k y = −∞

(7)

−∞

Fig. 2(a) is the normalized autocorrelation function of the grayscale image in Fig. 1(a) and the corresponding power spectrum is shown in Fig. 2(b). 2.2. Lineal-path function A two-phase random media can be simply expressed by a unit step function: { 1 if r belongs to white phase I (r) = 0 if r belongs to black phase

(8)

where vector r is spatial position of a pixel (in 2D case) or voxel (in 3D case) in the image. For a digital image with size m × n, let p represent the volume fraction of the white phase that is interested in reconstructing, then: p=

m ∑ n ∑

I (i, j)/(m × n)

(9)

i=1 j=1

The two-phase random media shown in Fig. 1(c), where the white pixels represent the aggregate and black pixels represent the cement matrix, according to Eq. (9), the aggregate volume fraction is p = 0.6.

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

5

Fig. 3. The lineal-path function of the model in Fig. 1(c).

Lineal-path function L ( j) (r1 , r2 ) contains substantial structural information, such as the size and connectedness information of the heterogeneities embedded in matrix, as its definition is the probability of finding a line segment between r1 and r2 lies entirely in the phase j and a detailed description can be found in [25]. For an isotropic media, the lineal-path function only depends on the distance r. It should be noted that the general lineal-path function is not sufficient to characterize the anisotropic microstructure. In this case, the directional correlation functions should be applied. Interested readers are referred to Ref. [26] for detailed descriptions of directional correlation functions and their properties. Because of its simplicity and practicability, lineal-path function is widely used for the characterization of clusters in random media structures [27]. In this paper, the heterogeneity of random media is statistically quantified by this function. Hereafter, we primarily assign the subscript to the phase that we are interested in reconstructing. For all media with a volume fraction of p, clearly, L (0) = p

(10)

The lineal-path function of the sample in Fig. 1(c) is plotted in Fig. 3. The function has a value of 0.6 at 0, that is, the volume fraction of aggregate in the sample is about 0.6. As can be seen where r is greater than 60, the corresponding function values are equal to 0 which indicates that aggregate particle size range is within 60 pixels along the sampling direction. 3. Random media reconstruction 3.1. Gray-level random distribution model reconstruction In this section, the Monte Carlo method is used to simulate the randomly distributed gray level in space. Suppose the size of the space in which the gray level to be generated is L x × L y × L z , the discrete points with equal spacing are L, M and N, and the corresponding intervals are ∆x, ∆y and ∆z, respectively. Based on Kuga’s method [28], the gray level g(x, y, z) in three-dimensional space is given as follows: ⏐ L/2 M/2 N /2 ∑ ∑ ∑ ⏐ ( ) [ ( 1 g (x, y, z) = ⏐⏐ F k xl , k ym , k zn exp i k xl x + k ym y L x L y L z l=−L/2+1 m=−M/2+1 n=−N /2+1 (11) ⏐ )] ⏐ +k xl z ⏐⏐

6

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15 Table 1 Parameters of the reference model shown in Fig. 1. Statistical characteristic values

Parameters of autocorrelation function

Volume fraction

Mean

Standard deviation

a (cm)

b (cm)

σ (x, y)

Aggregate

130.4

44.4

0.3

0.3

N(0.75, 0.1)

0.6

where ( ) [ ( )]1/2 F k xl , k ym , k zn = 2π L x 2π L y 2π L z S k xl , k ym , k zn { [N (0, 1) + i N (0, 1)] , l ̸= 0, L/2 and m ̸= 0, M/2 and n ̸= 0, N /2 × N (0, 1), otherwise √ k xl = 2πl/L x k ym = 2π m/L y k zn = 2π n/L z i = −1

(12)

where k xl , k ym and k zn are the discrete set of spatial frequencies; S(k x , k y , k z ) is the power spectrum of a random process; N (0, 1) denotes a sequence of normally distributed numbers in [0, 1] with zero mean and unity standard deviation. Therefore, an appropriate power spectrum that can characterize the random process must be selected first before generating randomly distributed gray level according to Eq. (12). A variety of power spectra is available to simulate random processes, notably in the fields of composite materials, electromagnetism, computational mechanics, machinery manufacturing, geoscience, and so on. For example, the Pierson–Moskowitz (PM) spectrum for sea surface simulation [29]; the fractal spectrum which is suitable for simulating rough surface with fractal features [30]; some other power spectrum is obtained indirectly by Fourier transform of autocorrelation function based on Wiener–Khinchin theorem [31,32]. Although it applies only to stationary random process, it is used extensively due to its maturity, ubiquity and performance. Ellipsoidal autocorrelation function is a common statistic characterizing random process in space [33], the expression is given [34]: [ ] 1/r G (x, y) = exp −(x 2 /a 2 + y 2 /b2 ) (13) where a and b are the autocorrelation lengths, r is a constant, called roughness factor. The ellipsoidal autocorrelation function can be further divided into Gaussian, exponential and hybrid type, the corresponding r values are 1, 2 and other values, respectively. In this paper, an improved 3D autocorrelation function is proposed as follows: [ ( )1/σ (x,y,z) ] G (x, y, z) = exp − x 2 /a 2 + y 2 /b2 + z 2 /c2 (14) where a, b and c are the autocorrelation lengths in the x, y and z directions, respectively. σ (x, y, z) is a custom roughness function and its value is a variable changing with the variation of space position. The roughness function can be used to characterize the heterogeneity of random media more flexibly, and more detailed structures can be generated than a single roughness value. For a better understanding of the performance of the roughness function σ , taking one-dimensional random gray level as an example, let the autocorrelation length a = 1, and substitute the power spectrum obtained by Fourier transform of Eq. (14) into Eq. (12), and then use Eq. (11) to generate gray curve with different roughness from the same set of random numbers, as shown in Fig. 4. According to the normalized autocorrelation function of the sample, the appropriate autocorrelation length can be selected, whose value is approximately equal to the distance from the maximum value 1 to exp−1 [35], more detailed descriptions can be found in the literature [33]. The roughness function adopts empirical value, and after several tests, the normal distribution sequence with 0.75 mean and 0.1 standard deviation is finally selected. The statistical parameters of the model shown in Fig. 1 are given in Table 1: 3.2. Gray-level random distribution model binarization To obtain a two-phase model, an image segmentation algorithm should be employed for the gray-level random distribution model to separate the target phase from the background matrix. This process includes two steps:

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

7

Fig. 4. One-dimensional random distribution model of gray level: (a) a = 1, σ (x) = 1; (b) a = 1, σ (x) = N(1,0.1) (a normal distribution with a mean of 1 and a variance of 0); (c) a = 1, σ (x) = 0.5; (d) a = 1, σ (x) = N(0.5,0.1) (a normal distribution with a mean of 0.5 and a variance of 0.1).

transform the gray model into a binary model containing only black and white pixels; adjust the volume fraction of the target phase to that of the reference model. Considering the operation efficiency, we chose threshold segmentation algorithm for binary processing. Both global threshold algorithm [36] and local threshold algorithm [37] can realize binarization of the gray model. Even one can customize a certain gray value to divide the whole model pixels. In this step, we are more concerned about the statistical characteristics of the target phase similar to that of the reference model. Therefore, the determination of threshold value depends on whether the obtained binary model is consistent with the reference model in terms of the number, spatial distribution and size of the aggregate. After the binarization processing, a two-phase model is obtained, in which the volume fraction of the target phase cannot be guaranteed to be the same as that of the sample. Therefore, further adjustment need to be done to match that of the reference model. Increasing or decreasing voxels of the target phase randomly in the whole space would introduce noise (i.e., isolated black or white voxels), which would worsen the model. The corresponding solution is given in literature [19], that is, increasing or decreasing voxels only on the edge of an existing phase cluster. In this step, it is still important to note that the cluster size of reconstruction should not exceed that of the sample. 3.3. Two-phase model optimization The theoretical autocorrelation function used to reconstruct the gray level distribution is the approximation of the actual autocorrelation function of the reference sample. Therefore, further processing is needed to reduce the

8

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

Fig. 5. Illustrative example of the pixel selection process: (a) object of random perturbation with arbitrary morphological heterogeneities; (b) white pixel selection process; (c) black pixel selection process.

Fig. 6. (a) Reconstruction of gray level model. (b) Intermediate two-phase model. (c) Final optimized model.

approximation error. In this step, simulated annealing is chosen to further adjust the reconstruction model. The binary image of sample is taken as the reference model, and the intermediate two-phase model obtained in Section 3.2 is to be optimized. The lineal-path function is used as input variables in the optimization process. Let f 0 (ri ) represent the function value of the reference model and f (ri ) represent the function value of the reconstruction model. To

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

9

Fig. 7. Comparisons of the lineal-path function before SA.

evolve the intermediate model towards the reference model, a stochastic perturbation process should be performed and the new function value f ′ (ri ) after perturbing is calculated. The “energy” in the simulated annealing is defined as: ∑ E= (15) [ f (ri ) − f 0 (ri )]2 i ′

E =

∑[

]2 f ′ (ri ) − f 0 (ri )

(16)

i

∆E = E ′ − E The random perturbation is accepted with probability p(∆E) via the Metropolis method: { 1, ∆E ≤ 0 p(∆E) = exp(−∆E/T ), ∆E > 0

(17)

(18)

where T is temperature and determines the cooling schedule. The random perturbation in this case refers to the process of randomly selecting a voxel from the black pixel region and the white voxel region respectively and then exchanging the states of these two voxels. In the traditional stochastic reconstruction algorithm, the pair of voxels used for interchange are arbitrarily selected from the whole model [18], or to accelerate the evolutionary process, literature [38] proposed a multiple interchanging process (MIP) that can exchange multiple voxel pairs simultaneously. However, both perturbation methods are applicable to the randomly distributed points models, which are not applicable to the intermediate model in this step. To solve this problem, the surface optimization rule should be employed, where the voxel pairs are randomly selected from the edge of clusters. The selection process of black and white voxel pairs is briefly illustrated in Fig. 5. Firstly, randomly select a cluster and find its inner edge pixels, then a white pixel is picked from the inner edge pixels. The black pixel is picked from the outer edge pixels in the same way. The two selected pixels are then swapped. The random perturbation in each optimization iteration must go through such a process. Such a random perturbation process can not only avoid the introduction of noise (i.e., isolated black or white voxels), but also reduce the unnecessary interchange process and accelerate the evolution of the intermediate model. The algorithm flow of random media reconstruction is summarized as follows: Step 1. Given a grayscale image of the reference sample, calculate the statistics such as mean value and root mean square. Obtain the power spectrum by fast Fourier transform and then the autocorrelation function is achieved by Wiener–Khinchin theorem. Estimate parameters of autocorrelation function.

10

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

Fig. 8. Optimization histories of different reconstruction methods.

Fig. 9. 3D reference model of concrete: (a) gray level model; (b) two-phase model.

Step 2. Threshold segmentation is carried out on the grayscale image to obtain a binary image of the sample. Calculate the volume fraction and plot the lineal-path function curve of the sample. Step 3. Plug estimated autocorrelation function (14) into Eq. (12) and reconstruct the gray-level random distribution model according to Eq. (11). Step 4. Perform a threshold cut to the gray-level random distribution model to generate a two-phase model and adjust the volume fraction of the target phase to match that of the sample. Step 5. A stochastic optimization algorithm is employed to refine the reconstruction model. For a better understanding of the entire reconstruction process, Fig. 6 shows several representative structures from a 2D random media reconstruction with Fig. 1 as the reference model.

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

11

Fig. 10. Reconstruction model R1: (a1) reconstruction of gray level model; (a2) intermediate two-phase model; (a3) final optimized model; (b1) XY section. (b2) XZ section. (b3) YZ section.

To demonstrate the capability of the proposed method, the original simulated annealing method based on the randomly distributed points model [18] is also tested. Both methods adopt the same temperature T in the process of simulated annealing optimization. Only the lineal-path function is incorporated in the “energy” E in Eq. (16). The lineal-path function before optimization is shown in Fig. 7. It is observed that the reconstructed model with the proposed method is already close to the target structure. The change of the energy E with the number of iterations is shown in Fig. 8. Obviously, the proposed method has absolute advantages in minimizing the energy to some small tolerance value with less iteration which means that this method can reconstruct the random media with higher efficiency. 4. 3D random media reconstruction In this section, the proposed algorithm is employed to reconstruct a 3D aggregate composite structure with a volume fraction of V = 50% based on the statistical information from the 3D CT samples shown in Fig. 9. For the limitation of the instrument on the scanning sample, the reference model is smaller than that get from digital camera in Fig. 1. Its dimensions are 3.5 cm × 3.5 cm × 3.5 cm containing 64 × 64 × 64 voxels. Nonetheless, the sample has enough statistics to achieve statistically equivalent reconstruction by assuming it is ergodic and stable stochastic processes. To demonstrate the capability of the proposed method, two 3D reconstruction random media with different size are presented. One is the same size as the reference model, hereafter denoted by R1. The other is eight times the size of the reference model, which is 7.0 cm × 7.0 cm × 7.0 cm, denoted by R2. Representative models obtained in reconstruction process of R1 and R2 are shown in Figs. 10 and 11, respectively. Following the algorithm flow, the initial gray level random distribution models are established based on Eq. (11). The smaller is shown in Fig. 10(a1) and the bigger one is shown in Fig. 11(a1). Then threshold value segmentation is used for binary processing of the gray level models and the intermediate two-phase models can be achieved

12

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

Fig. 11. Reconstruction model R2: (a1) reconstruction of gray level model; (a2) intermediate two-phase model; (a3) final optimized model; (b1) XY section. (b2) XZ section. (b3) YZ section.

Fig. 12. Comparisons of lineal-path function between the reference model and two reconstruction models before and after optimization. (a) reconstruction model R1 in Fig. 10; (b) reconstruction model R2 in Fig. 11.

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

13

Fig. 13. Autocorrelation functions of reference and reconstruction models in Figs. 7, 8, 9: (a1–b1) The autocorrelation along x direction; (a2–b2) the autocorrelation along y direction; (a3–b3) the autocorrelation along z direction.

by adjusting the white pixels to the same volume ratio as that of the reference model as shown in Figs. 10(a2) and 11(a2). The stochastic reconstruction procedure is used to further refine the intermediate model and only the lineal-path function is incorporated into E in Eq. (15). The final reconstructed two-phase models are shown in Figs. 10(a3) and 11(a3). Cross-sectional views shown in Figs. 10(b1–b3) and 11(b1–b3) demonstrate the internal structure of the reconstructed model is visually similar to the reference model. The lineal-path functions of the reference model and the reconstructed model before and after optimization are illustrated in Fig. 12. As mentioned in Section 3.3, the autocorrelation function used to reconstruct the gray scale model is an approximation of the actual function of the sample. Moreover, the autocorrelation function of the reconstructed model has changed after threshold segmentation. Therefore, it is necessary to verify whether the autocorrelation function of the final reconstructed model is basically consistent with the sample or not. The autocorrelation functions are calculated directly from the reference and reconstructed models based on the Wiener–Khinchin theorem and plotted along three directions as shown in Fig. 13. It should be noted that for two-phase model, the value of the autocorrelation function is also related to the size of model. Therefore, a portion of R2 with the same size as sample is involved in comparison. Under the ergodic hypothesis, this part used to calculate the autocorrelation function can be considered to have the same statistical characteristics as the whole reconstructed model. The autocorrelation functions of the reference model and reconstruction models are in good agreement, which indicates the reliability of the proposed method. 5. Conclusions In this paper, an efficient algorithm has been developed for reconstructing two-phase heterogeneous materials with random morphological characteristics. The new method is based on the theory of random process and stochastic optimization algorithm. Random process theory mainly applies fast Fourier transform to generate a 3D initial gray

14

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

level random distribution model with high computational efficiency. Meanwhile a new autocorrelation function form is introduced to the procedure to characterize more detailed features of random media. The stochastic perturbation is improved to accelerate convergence of simulated annealing algorithm by greatly reducing unnecessary interchanges during optimization procedure. The proposed method was applied in reconstructing 3D concrete structure with two different sizes. The good match of the autocorrelation function with that of the sample model suggests the validity of this approach. It should be noted that the reference model is assumed to be ergodic. However, in practice, for random heterogeneous media, the inclusions embed in the matrix vary in size, shape and random distribution. Therefore, the volume fractions of the two phases in different samples may not be equal. Furthermore, the binary process can also change the volume fraction. But this is not a fundamental deficiency of the method and it can be very trustworthy if the volume fraction is in error range. Acknowledgment This work was supported by the Fundamental Research Funds for the Central Universities, China, No.: 2018ZDPY08. References [1] L.J. Jacobs, J.O. Owino, Effect of aggregate size on attenuation of Rayleigh surface waves in cement-based materials, J. Eng. Mech. 126 (11) (2000) 1124–1130. [2] A.B. Bouda, S. Lebaili, A. Benchaala, Grain size influence on ultrasonic velocities and attenuation, Ndt E Int. 36 (1) (2003) 1–5. [3] L. Liu, D. Shen, H. Chen, et al., Aggregate shape effect on the diffusivity of mortar: a 3D numerical investigation by random packing models of ellipsoidal particles and of convex polyhedral particles, Comput. Struct. 144 (2014) 40–51. [4] J.A. Bogas, M.G. Gomes, A. Gomes, Compressive strength evaluation of structural lightweight concrete by non-destructive ultrasonic pulse velocity method, Ultrasonics 53 (5) (2013) 962–972. [5] A. Benouis, A. Grini, Effect of Concrete Mixtures on Estimation of Porosity by Ultrasonic Velocity, in: Nondestructive Testing of Materials and Structures, Springer, Dordrecht, 2013, pp. 309–316. [6] Y. Li, Y. Li, Z. Guan, et al., Elastic modulus damage model of cement mortar under salt freezing circumstance based on X-ray CT scanning, Constr. Build. Mater. 191 (2018) 1201–1209. [7] R. Bostanabad, Y. Zhang, X. Li, et al., Computational microstructure characterization and reconstruction: Review of the state-of-the-art techniques, Prog. Mater. Sci. (2018). [8] M.D. Rintoul, S. Torquato, Reconstruction of the structure of dispersions, J. Colloid Interface Sci. 186 (2) (1997) 467–476. [9] L. Liu, Z. Li, S. Arcone, et al., Radar wave scattering loss in a densely packed discrete random medium: Numerical modeling of a box-of-boulders experiment in the Mie regime, J. Appl. Geophys. 99 (2013) 68–75. [10] H. Li, P.E. Chen, Y. Jiao, Accurate reconstruction of porous materials via stochastic fusion of limited bimodal microstructural data, Transp. Porous Media 125 (1) (2018) 5–22. [11] Y. Ju, Y. Huang, J. Zheng, et al., Multi-thread parallel algorithm for reconstructing 3D large-scale porous structures, Comput. Geosci. 101 (2017) 10–20. [12] X.F. Wang, Z.J. Yang, J.R. Yates, et al., Monte Carlo simulations of mesoscale fracture modelling of concrete with random aggregates and pores, Constr. Build. Mater. 75 (2015) 35–45. [13] J.W. Feng, S. Cen, C.F. Li, et al., Statistical reconstruction and Karhunen–Loève expansion for multiphase random media, Internat. J. Numer. Methods Engrg. 105 (1) (2016) 3–32. [14] D. Stoyan, Random sets: models and statistics, Internat. Statist. Rev. 66 (1) (1998) 1–27. [15] L.M. Pant, S.K. Mitra, M. Secanell, Stochastic reconstruction using multiple correlation functions with different-phase-neighbor-based pixel selection, Phys. Rev. E 90 (2) (2014) 023306. [16] Y. Jiao, F.H. Stillinger, S. Torquato, A superior descriptor of random textures and its predictive capacity, Proc. Natl. Acad. Sci. 106 (42) (2009) 17634–17639. [17] E.Y. Guo, N. Chawla, T. Jing, et al., Accurate modeling and reconstruction of three-dimensional percolating filamentary microstructures from two-dimensional micrographs via dilation-erosion method, Mater. Charact. 89 (2014) 33–42. [18] C.L. Yeong, S. Torquato, Reconstructing random media, Phys. Rev. E 57 (1) (1998) 495–506. [19] Z. Jiang, W. Chen, C. Burkhart, Efficient 3D porous microstructure reconstruction via Gaussian random field and hybrid optimization, J. Microsc. 252 (2) (2013) 135–148. [20] L.M. Pant, S.K. Mitra, M. Secanell, Stochastic reconstruction using multiple correlation functions with different-phase-neighbor-based pixel selection, Phys. Rev. E 90 (2) (2014) 023306. [21] A. Maiti, D. Chakravarty, K. Biswas, et al., Development of a mass model in estimating weight-wise particle size distribution using digital image processing, Int. J. Min. Sci. Technol. 27 (3) (2017) 435–443. [22] Z.Q. Yue, S. Chen, L.G. Tham, Finite element modeling of geomaterials using digital image processing, Comput. Geotech. 30 (5) (2003) 375–397. [23] L. Klimeš, Correlation functions of random media, Pure Appl. Geophys. 159 (7–8) (2002) 1811–1831.

W. Zhang, L. Song and J. Li / Computer Methods in Applied Mechanics and Engineering 354 (2019) 1–15

15

[24] P.R. Nayak, Random process model of rough surfaces, J. Lubr. Technol. 93 (3) (1971) 398–407. [25] S. Torquato, B. Lu, Chord-length distribution function for two-phase random media, Phys. Rev. E 47 (4) (1993) 2950. [26] Y. Jiao, N. Chawla, Modeling and characterizing anisotropic inclusion orientation in heterogeneous material via directional cluster functions and stochastic microstructure reconstruction, J. Appl. Phys. 115 (9) (2014) 093511. [27] D. Li, Review of structure representation and reconstruction on mesoscale and microscale, Jom 66 (3) (2014) 444–454. [28] Y. Kuga, P. Phu, Experimental studies of millimeter-wave scattering in discrete random media and from rough surfaces, Prog. Electromagn. Res. 14 (1996) 37–88. [29] S.L. Broschat, The phase perturbation approximation for rough surface scattering from a pierson-moskowitz sea surface, IEEE Trans. Geosci. Remote Sens. 31 (1) (1993) 278–283. [30] S. Kumar, G.S. Bodvarsson, Fractal study and simulation of fracture roughness, Geophys. Res. Lett. 17 (6) (1990) 701–704. [31] C.F. Li, Y.T. Feng, D.R.J. Owen, et al., A Fourier–Karhunen–Loève discretization scheme for stationary random material properties in SFEM, Internat. J. Numer. Methods Engrg. 73 (13) (2008) 1942–1965. [32] H.M. Stanley, T. Kato, An FFT-based method for rough surface contact, J. Tribol. 119 (3) (1997) 481–485. [33] L.T. Ikelle, S.K. Yung, F. Daube, 2-D random media with ellipsoidal autocorrelation functions, Geophysics 58 (9) (1993) 1359–1372. [34] X. Xi, Y. Yao, Simulations of random medium model and intermixed random medium, Earth Sci.-J. China Univ. Geosci. 27 (1) (2002) 67–71, (in Chinese). [35] G.U.O. Shi-Li, J.I. Meng-En, Z.H.U. Pei-Min, et al., Study on multiphase discrete random medium model and its GPR wave field characteristics, Chin. J. Geophys. 58 (4) (2015) 422–435. [36] Nobuyuki OHTSU, A threshold selection method from gray-level histograms, IEEE Trans. Syst. Man Cybern. 9 (1) (1979) 62–66. [37] W.B. Lindquist, A. Venkatarangan, Investigating 3D geometry of porous media from high resolution images, Phys. Chem. Earth A 24 (7) (2007) 593–599. [38] Y. Ju, J. Zheng, M. Epstein, et al., 3D numerical reconstruction of well-connected porous structure of rock using fractal algorithms, Comput. Methods Appl. Mech. Engrg. 279 (2014) 212–226.