Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2024.49.0077 Acta Polytechnica CTU Proceedings 49:77–84, 2024 © 2024 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague IMAGE-BASED RANDOM FIELDS IN NUMERICAL MODELING OF HETEROGENEOUS MATERIALS David Šilhánek, Jan Sýkora∗ Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 160 00 Prague, Czech Republic ∗ corresponding author: jan.sykora.1@fsv.cvut.cz Abstract. Modeling of spatial variability using random fields is nowadays standard in computational modeling of heterogeneous materials. The only difficulty remains in determining the input parameters of the random field in order to represent the modeled material morphology as accurately as possible. The most frequent method of constructing a random field relies on various covariance kernels that primarily use as input the values of the correlation lengths, which are typically estimated ad-hoc. However, another possibility is to extract input parameters from the modeled material morphology. Moreover, this approach enables to calculate the covariance kernel itself. The performance of the image-based procedure is compared with standard methods of random field construction on the response of the two-phase elastic material. Keywords: Karhunen-Loève expansion, random field, Monte Carlo, two-phase elastic material. 1. Introduction The stochastic finite element method is an extension of the deterministic finite element method that in- troduces randomness into the computational scheme, see [1, 2]. Mostly, the input parameters are cast in a probabilistic setting, and random fields (see [3–5]) play a crucial role in this because they allow us to incorporate randomness and spatial variability. So, our goal is to investigate the efficiency of constructing a random field using several different strategies on the material response, which is here represented by the elastic problem of the two-phase medium. 2. Estimation of the RVE cell The first step in our work is estimating the dimen- sions of the representative volume element (RVE), which means an area with statistically homogeneous properties or at least with properties insignificantly affected by a particular choice. Let us consider a two-dimensional heterogeneous composite mate- rial with two phases represented by an integer ma- trix containing only 0 (black color) and 1 (white color), see Figure 1 depicting several cut-outs of the original medium. Each matrix member corre- sponds to one pixel in the image and simultaneously to one rectangular element in the finite element mesh. The domain under the study contains different-sized circular white-labeled inclusions with the following prescribed grain-size distribution curve: diameter 2 px (10 %), 4 px (10 %), 8 px (15 %), 16 px (20 %), 24 px (20 %), and 32 px (25 %). Its overall dimensions are 10 000 × 10 000 px, and the volume fraction of inclusions is c1 = 0.5. Here, the determination of the RVE size stems from an engineer’s approach a priori estimating the dimen- (a). 150 × 150 px (b). 100 × 100 px (c). 50 × 50 px Figure 1. Samples of the windows with predefined di- mensions taken randomly from the original domain. sions based on easily observable parameters, namely the volume fraction of the white phase and the opti- mal correlation lengths. Authors are aware of a more correct approach to identifying RVE size from the statistical analysis of model responses. However, this does not have to be always possible and brings ob- 77 https://doi.org/10.14311/APP.2024.49.0077 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en David Šilhánek, Jan Sýkora Acta Polytechnica CTU Proceedings Figure 2. The convergence graph computed for the volume fraction of the inclusions as a function of the window size. The solid line denotes the mean and the light-hatched area represents the ± standard deviation. stacles, especially for the time-demanding numerical simulations. The first study addressing the optimal RVE size is built upon the parameter of volume fraction. A win- dow of predefined dimensions is randomly placed 1 000 times in the original medium. Consequently, the dif- ference in the studied properties and their statistical moments are then evaluated on the obtained dataset. The entire procedure is repeated for various dimen- sions of the placed window. The plotted graphs serve for the visual determination of optimal RVE dimen- sions. The results obtained from the first study are shown in Figure 2. The mean of the volume fraction is almost constant for all predefined window sizes, contrary to the standard deviation, which steeply converges with increasing dimensions, and the values around the window size of 100 × 100 px reach the fi- nal plateau. So, optimal RVE sizes in the view of the volume fraction are the dimensions above 100×100 px. The second study focuses on estimating the RVE size based on the optimal correlation lengths (Lx, Ly). The algorithmic procedure is similar to the previous case, except that an additional step of identifying the correlation lengths for every placed window with pre- defined dimensions is added to the calculation. This step does not significantly increase the computational costs, since the relatively small window sizes enter the optimization process of correlation lengths. The built-in Matlab function is used as an optimizer for minimizing the value of error εLopt between the corre- lation functions CG and CT P P : εLopt = 1 nx · ny nx∑ i=1 ny∑ j=1 |CT P P,ij −CG,ij(Lx, Ly)|, (1) where nx denotes the window dimension in x-direction, ny denotes the window dimension in y-direction, CT P P is the image-based correlation kernel stored as a matrix and enumerated according to Equation (4), CG is the Gaussian correlation kernel based on Equa- tion (3). The optimal correlation lengths are plotted in Fig- ure 3. As we can notice, the values of optimal correla- tion length become constant for window dimensions greater than 70 × 70 px. Moreover, the mean values of correlation length for the same threshold also cor- respond to the optimal values computed for the entire original domain (10 000×10 000 px), as one can ensure from Figure 4, where the absolute values of the error are plotted. The optimal values of correlation length are primarily affected by the spatial distribution of inclusions, the sizes of the inclusions, and the value of the volume fraction. So, the observed evolution of error curves coincides with the study findings examin- ing the effect of volume fraction on the window size. Consequently, the influence of window size on optimal correlation lengths is negligible for dimensions above 100 × 100 px and almost disappears for the window sizes above 150 × 150 px. Usually, the S2 function is computed using the discrete fast-Fourier transform algorithms (DFT) [5] to save an enormous amount of computational time. However, this procedure requires the periodic medium as an input, which is not our case of windows randomly placed into the original medium, thus the effect of the non-periodic domain on the final S2 DFT-based values needs to be carefully investigated. For this reason, the algorithm computing S2 values for the non-periodic medium has been implemented as a comparison tool with DFT-based values. The computational testing protocol is the same as in previous studies and the ob- served values of error are obtained from the following equation as: εS2,p,S2,np = 1 nx · ny nx∑ i=1 ny∑ j=1 |Speriodic 2,ij − Snonperiodic 2,ij |, (2) where nx denotes the window dimension in x-direction, 78 vol. 49/2024 Image-based RF in numerical modeling of heterogeneous materials Figure 3. The convergence of optimal correlation lengths as a function of the window size. The solid line denotes the mean and the light-hatched area represents the ± standard deviation. Figure 4. The difference between the optimal correlation length obtained from the entire original domain and mean values of the correlation length computed for the various dimensions of the placed window. ny denotes the window dimension in y-direction, Speriodic 2,ij is the i, j-th element of two-point probability matrix computed with the assumption of periodicity according to Equation (5), and i, j-th element of two-point probability matrix obtained for the non- periodic medium with the help of Equation (7). As expected, the curves in Figure 5 show good con- vergence of values of error with higher window sizes. Therefore, the effect of the periodic assumption on the optimized values of correlation lengths is negligible for the window sizes of our interest, and the fast eval- uation of the two-point probability based on the DFT algorithm can be performed in our computational scheme. Considering the outcomes mentioned above, we decided to proceed further with an RVE size of 100 × 100 px. It is a reasonable compromise between the obtained results – the variance of the volume fraction and the optimized correlation length – and computational cost. Another enlarging of the win- dow area would not bring a huge benefit since the values of mean computed for the correlation length and volume fraction are further invariable, and their standard deviations reach almost the plateau. 3. Numerical Example This numerical example is devoted to the assessment of several proposed techniques of random field con- structions in stochastic material modeling. Since the efficiency of the random fields is demonstrated mostly on the model inputs, the emphasis is here mainly on checking the quality of the prediction of mate- rial behavior. To begin with, the ideal elastic model simulating the mechanical behavior is employed for both material phases. Their material properties – E Young’s modulus and ν Poisson’s ratio – differing by an order of magnitude, are summarized in Table 1. The geometrical domain with applied boundary conditions is depicted in Figure 6. The Dirichlet boundary conditions are imposed on the left side of the domain as ux = 0, on the right side as incre- mentally increasing up to ux = 1.5 · 10−6 m, and 79 David Šilhánek, Jan Sýkora Acta Polytechnica CTU Proceedings Figure 5. The comparison of the periodic and non-periodic assumptions on the values of S2 function. Phase E [MPa] ν [-] Black phase, 0 1 325 0.43 White phase, 1 25 000 0.20 Table 1. Elastic properties of black and white phase representing the material constants of polypropylene composite materials [6]. Figure 6. A scheme of the tensile st problem. finally uy = 0 at bottom left corner. The entire do- main with dimensions of 100 × 100 px is discretized into 10 000 finite rectangular elements. The observed responses are displacements ux, uy, and average stresses σx. 3.1. Random field The spatial variability of material properties is mod- eled here by the random fields, which are essentially dependent on the auto-covariance function. Since the correct determination of the auto-covariance function plays a crucial role, several methods of constructing this key function are tested within this study. The classical approach, published in literature [4, 5, 7], is based on the approximation of the covariance function by the Gaussian correlation kernel, which is a function of a material variance σ2, correlation lengths Lx, Ly, and the vector of relative position r = (x−x0, y −y0), where (x0, y0) are coordinates of an arbitrarily chosen point (pixel in our case). The Gaussian kernel is then expressed as: CG(r) = σ2 · exp ( − (x − x0)2 2 · L2 x − (y − y0)2 2 · L2 y ) . (3) The input correlation lengths are chosen based on the expert’s guess or their values can be estimated from the optimization process minimizing the absolute difference between the Gaussian covariance kernel (see Equation (3)) and image-based covariance kernel, see [5], derived from a two-point probability function. The latter-mentioned approach represents relatively new concepts of extracting the spatial covariance from images according to the following formula: CT P P (r) = (κ0 − κ1)2 · (S1 2(r) − (c1)2), (4) where κf is the selected material property (e.g. Poisson’s ratio) of given phase f , c1 denotes the volume fraction of the white phase, S1 2 is the two-point probability function computed by the Fast Fourier Transformation, see [7], as: Sf 2 (r) = IDFT { DFT{χf (r)}DFT{χf (r)} } nx · ny , (5) where DFT is the Discrete Fourier Transformation, IDFT is the inverse Discrete Fourier Transformation, · · · stands for the complex conjugate, nx · ny is the image area, f is the given phase, χf is the characteristic function defined as: χf (x) = { 1 if the point x lies in the phase f, 0 otherwise. (6) 80 vol. 49/2024 Image-based RF in numerical modeling of heterogeneous materials If the non-periodic medium is considered, the formula used for the evaluation of a two-point probability function arises directly from its definition: Sf 2 (x, x′) = ∑(2·nx−dx−1) x=0 ∑(2·ny−dy−1) y=0 χf (x)χf (x′) (2nx − dx) · (2ny − dy) , (7) where x = (x, y) are the coordinates of the starting point, x′ = (x+dx, y+dy) are the coordinates of the ending point. Finally, the auto-covariance matrix entering the random field construction is obtained as a resulting covariance computed from Equations (3) or (4), which is further reassembled into the row vector. Then, this reshaped product determines the first row in a sym- metric square Toeplitz matrix representing the auto- covariance matrix such that Ci,j = Ci+1,j+1. Another possibility of computing the auto-covariance matrix is directly from the definition itself. Based on the as- sumption that our original domain is sufficiently large enough, many independent samples are collected via a window of predefined dimensions randomly placed in the original structure. Each of these cut-outs can be then rearranged into a row vector, and subsequently, the auto-covariance is calculated as follows: CS(X, X ′) = E[(X − E(X))(X ′ − E(X ′))]. (8) 3.1.1. Karhunen-Loève Expansion The Karhunen-Loève expansion is an extremely useful tool representing the stochastic process as an infinite linear combination of orthogonal functions. Based on the spectral decomposition of the discretized form of the auto-covariance function C(X, X ′), and the or- thogonality of eigenfunctions ϕi, the Gaussian random field of parameter κ(ω) can be written as: κ(ω) ≈ κµ + s nϕ∑ i=1 √ φiϕiξi(ω), (9) where ξi(ω) is the standard independent and identically distributed random variable, κµ is the mean of material parameters, s is the standard deviation of material parameters, φi represents a corresponding eigenvalue, nϕ is the number of eigenmodes. Smaller values of nϕ lead to more significant dimen- sionality reduction, and oppositely, higher values give better descriptions of the resulting random field. The log-normal formulation of the random field κln(ω) stemming from the Gaussian expression is given by: κln(ω) ≈ exp ( κµ,g + sg ∑nϕ i=1 √ φiϕiξi(ω)√∑nϕ i=1(√φiϕi)2 ) , (10) Tag Label Cov. Correlation Trimmed kernel length #1 Gauss (3) Lopt No #2 Gauss T (3) Lopt Yes #3 TPP (4) − No #4 TPP T (4) − Yes #5 S (8) − No #6 S T (8) − Yes #7 Gauss + (3) 2 · Lopt No #8 Gauss + T (3) 2 · Lopt Yes #9 Gauss ++ (3) 10 · Lopt No #10 Gauss ++ T (3) 10 · Lopt Yes #11 Gauss − (3) 1 2 · Lopt No #12 Gauss − T (3) 1 2 · Lopt Yes #13 Gauss −− (3) 1 10 · Lopt No #14 Gauss −− T (3) 1 10 · Lopt Yes Table 2. Description of all tested random field con- structions. The first six versions #1–#6 character- ized the main approaches of constructing random fields. The rest of them #7–#14 served for simulating the wrong initial guess of experts about correlation lengths. where κµ,g and sg are the Gaussian mean and stan- dard deviation, respectively, computed from their log- normal counterparts (κµ,ln, sln) according to the fol- lowing equations, see [2]: κµ,g = ln ( κµ,ln ) − 1 2s2 ln, (11) sg = √ ln (( sln κµ,g )2 + 1 ) . (12) 3.1.2. Description of tested random field All tested versions of random field constructions are summarized in Table 2. Overall, 14 different ver- sions have been examined with the following specifi- cations: The labels – Gauss, TPP, S – refer to the equations employed for assembling the covariance ker- nels and auto-covariance. Since the Karhunen-Loève expansion delivers the real-valued random fields, the additional label – T – is introduced as another op- tion referring to value filtering. It means that real- valued random fields are trimmed back to the bi- nary values preserving the volume fraction. Finally, the correlation lengths indicate the values entering the covariance kernels. The versions with the labels – plus and minus – serve for simulating the wrong initial guess of experts, i.e. the underestimations and overestimations of correlation lengths, respec- tively. 3.2. Homogenization bounds Since the initial guess of experts about the inputs entering the covariance kernels is missing here, the theoretical homogenization bounds are exploited for the comparison of our proposed techniques with easily 81 David Šilhánek, Jan Sýkora Acta Polytechnica CTU Proceedings Figure 7. The evolution of average stresses as a function of loading steps, nϕ = 50. Figure 8. The values of error computed for the average stress as a function of the number of eigenmodes. computed rules. The first homogenization bound is the rule of mixtures, which is given by: κhom = κ1 · c1 + κ1 · (1 − c1), (13) and the inverse rule of mixtures is defined as: κhom = ( c1 κ1 + 1 − c1 κ0 )−1 . (14) κhom is the homogenized material property, and the results stemming from these two equations are labeled as ROM and IROM, respectively. 4. Results For the comparison, it is necessary to introduce the metrics of error. They represent the absolute values calculated as a difference between the reference set of material responses, and the response set having the random fields as an input. The equations evaluating the values of error for the average stress σx and its variance var(σx) are expressed as: εσx = ∣∣∣∣σx,real − σx,RF σx,real ∣∣∣∣, (15) and: εvar(σx) = ∣∣∣∣var(σx,real) − var(σx,RF) var(σx,real) ∣∣∣∣, (16) where the subscript real represents the reference set assembled from 1 000 Monte Carlo simulations with a domain size of 100×100 px, which is randomly taken from the original structure. Then, the subscript RF refers to set collected for the 1 000 Monte Carlo simu- lations utilizing the random fields as an input. The first graph (Figure 7) illustrates the evolution of average stresses as a function of loading steps for different constructions of random fields. The results are plotted for the fixed number of eigenmodes, which is equal to 50. The true values of the average stresses computed from the reference set are shown by the solid blue line and the homogenization bounds are depicted as dashed blue lines. The first six versions #1–#6 of random field construction are only compared for the sake of clarity. From the presented graph we can see that the trimmed versions show good correspondence with the reference values, in contrast to values carried out for real-valued random fields, which shift to the upper homogenization bound. A similar conclusion is observed from Figure 8 depicting the values of error computed for the average stress depending on the number of eigenmodes. As we can see from this graph, including more eigenmodes in the Karhunen-Loève expansion does not tend to have smaller values of error; however, it has a significant impact on values of error 82 vol. 49/2024 Image-based RF in numerical modeling of heterogeneous materials Figure 9. The variance error in dependence on the number of eigenvectors included in KLE. Figure 10. The ratio of the variances computed for the average stresses of the reference set and set assembled from the random field-based simulations #1–#6, nϕ = 250. Figure 11. The average stress for RF with arbitrarily chosen correlation lengths nϕ = 50. computed for the stress variance illustrated in Figure 9. All trimmed versions achieve the minimum values for nϕ = 200, while the real-valued versions converge with the higher number of involved eigenmodes. The last comparison is finally depicted in Figure 10 showing the ratio between the variances of average stresses determined from the reference set and set assembled from the random field-based simulations, respectively. Other interesting findings arise from the study of the wrong-determined input correlation lengths caused by the expert guess. It can be seen from Figures 11 and 12 that the random fields with the incorrect corre- lation lengths as an input can have fatal consequences on the material response and can bring no benefit to the probabilistic analysis of the studied problem. From the predicted lines, the splitting into the real- valued and trimmed groups is more observed than the clustering based on the different covariance kernels. 83 David Šilhánek, Jan Sýkora Acta Polytechnica CTU Proceedings Figure 12. The ratio of the variances computed for the average stresses of the reference set and set assembled from the random field-based simulations #7-#14, nϕ = 250. Taken together, the values of error for the real-valued random fields are large in comparison to other meth- ods, and therefore they do not appear to be adequate approaches for simulating such type of material. 5. Conclusion Within this contribution, the efficiency of the image- based procedure of random field construction was studied with classical methods on the elastic prob- lem of the two-phase material. Good results have been obtained for trimmed versions of random fields. Here, our proposed techniques for identifying proper- ties from images, e.g. correlation lengths or image- based covariance, are beneficial and can be applied to such problems. Unfortunately, these random field con- structions can no longer exploit the potential of the Karhunen-Loève expansion, and their contribution to the stochastic finite method in terms of acceleration is not significant. Acknowledgements The authors are thankful for financial support from the Student Grant Competition of CTU, project No. SGS23/152/OHK1/3T/11 and the Czech Science Foundation, project No. 22-35755K. References [1] R. G. Ghanem, P. D. Spanos. Stochastic finite elements: A spectral approach. Springer New York, USA, 1991. https://doi.org/10.1007/978-1-4612-3094-6 [2] B. Rosić, H. G. Matthies. Computational approaches to inelastic media with uncertain parameters. Journal of the Serbian Society for Computational Mechanics 2(1):28–43, 2008. [3] R. J. Adler, J. E. Taylor. Random fields and geometry. Springer-Verlag New York, USA, 2007. https://doi.org/10.1007/978-0-387-48116-6 [4] J. Havelka, A. Kučerová, J. Sýkora. Dimensionality reduction in thermal tomography. Computers & Mathematics with Applications 78(9):3077–3089, 2019. https://doi.org/10.1016/j.camwa.2019.04.019 [5] M. Lombardo, J. Zeman, M. Sejnoha, G. Falsone. Stochastic modeling of chaotic masonry via mesostructural characterization. International Journal for Multiscale Computational Engineering 7(2):171–185, 2009. https: //doi.org/10.1615/IntJMultCompEng.v7.i2.70 [6] J.-H. Yun, Y.-J. Jeon, M.-S. Kang. Analysis of elastic properties of polypropylene composite materials with ultra-high molecular weight polyethylene spherical reinforcement. Materials 15(16):5602, 2022. https://doi.org/10.3390/ma15165602 [7] J. Havelka, A. Kučerová, J. Sýkora. Compression and reconstruction of random microstructures using accelerated lineal path function. Computational Materials Science 122:102–117, 2016. https://doi.org/10.1016/j.commatsci.2016.04.044 84 https://doi.org/10.1007/978-1-4612-3094-6 https://doi.org/10.1007/978-0-387-48116-6 https://doi.org/10.1016/j.camwa.2019.04.019 https://doi.org/10.1615/IntJMultCompEng.v7.i2.70 https://doi.org/10.1615/IntJMultCompEng.v7.i2.70 https://doi.org/10.3390/ma15165602 https://doi.org/10.1016/j.commatsci.2016.04.044 Acta Polytechnica CTU Proceedings 49:77–84, 2024 1 Introduction 2 Estimation of the RVE cell 3 Numerical Example 3.1 Random field 3.1.1 Karhunen-Loève Expansion 3.1.2 Description of tested random field 3.2 Homogenization bounds 4 Results 5 Conclusion Acknowledgements References