A numerical method for the generation of hierarchical Poisson Voronoi microstructures applied in micromechanical finite element simulations—part I: method

Poisson Voronoi (PV) tessellations as artificial microstructures are widely used in investigations of material deformation behaviors. However, a PV structure usually describes a relative homogeneous field. This work presents a simple numerical method for generating 2D/3D artificial microstructures based on hierarchical PV tessellations. If grains/particles of a phase cover a large size span, the concept of “artificial phases” can be used to create a more realistic size distribution. From case to case, detailed microstructural features cannot be directly achieved by commercial or free softwares, but they are necessary for a deep or thorough study of the material deformation behavior. PV tessellations created in our process can fulfill individual requirements from material designs. Another reason to use PV tessellations is due to the limited experimental data. Concerning the application of PV microstructures, four examples are given. The FE models and results will be presented in consecutive works, i.e. “part II: applications”.


Introduction
To investigate the material behavior, micromechanical material properties are essential, since such features, besides other factors, determine the overall material characteristics. Besides experiments, finite element (FE) simulations also play an important role in studying material behaviors. Due to limitedly measured data or other reasons, FE simulations with artificial microstructures can be as valuable as simulations with real microstructures. When a real microstructure B Y. Schneider yanling.schneider@imwf.uni-stuttgart.de 1 is applied, its phase composition may have a large discrepancy to the one of the sample. This discrepancy is caused by the limited size of a microstructure cut-out, e.g. as shown in Schneider et al. [1]. Artificial microstructures generally do not have the above mentioned composition misfitting between the applied structure and the real one. At least, it can be easily avoided. The application of the artificial microstructure is very useful for a computer-aided material design and for the structure optimization. It is more time efficient and less expensive than the application of the real microstructure. To evaluate the microstructural influences on material properties, it is important to statistically estimate the mean size, the size distribution and geometrical characteristics of each phase. A random microstructure can be approximately presented by the Voronoi tessellation [2], which is widely used in modelling and analysing cell-type structures. Voronoi tessellations are also known as Voronoi diagrams, Voronoi decompositions or Dirichlet tessellations. Tessellations are also called cells, grains, aggregates or mosaics. To be distinguished with grains in a crystalline material (e.g. Ag grains in Fig. 6), the term "grain" is not used for Voronoi tessellations in the current work. Voronoi tessellations have different classes/types, e.g. the Poisson Voronoi (PV), the Hardcore Voronoi and the Laguerre Voronoi [3][4][5][6][7]. Within the class of Voronoi tessellations, the most intensively examined model is the PV tessellations. For this type of tessellations, the most extensive set of existing results concerns the distribution of the cell quality [2,5].
Voronoi structures are used by different authors to simulate material mechanical behaviors, e.g. [8][9][10][11][12][13][14][15]. Fritzen et al. [8] reported a method of generating the 3D meshing of periodic Voronoi structures for polycrystalline aggregates. Their process was time efficient, but it required professional programming skills and no hierarchical Voronoi tessellations are included. Schneider et al. [9] used PV cells in their texture simulations for α-Fe-Cu composites, where a regular meshing was used. Falco et al. [10] applied Voronoi tessellations for the discretization of arbitrarily shaped threedimensional polycrystalline models. They also introduced an original approach to define any possible (concave or convex) shape of the final domain, independent from the initial configuration of the aggregate. To simulate the pervasive fracture of materials and structures, Bishop [11] used the randomly closed-packed Voronoi tessellations, which provide a regularized random network of facets for representing cracks. Based on the Voronoi tessellation method applied on the microscopy, an approach for the elastic property calibration of the macroscopic property is established by using the discrete-element-method [12]. This calibration approach connects the macroscopic elastic properties with microscopic properties. By applying combined real and Voronoi microstructures, a coupled approach of a crystal plasticity FE model with a phase-field model is used to simulate the grain recrystallization [13]. In this model [13], the 3D real microstructure used in an FE simulation is taken as the starting local morphology under a hot rolling deformation process. The deformed microstructure (from a real one) provided the geometric information for the generation of Voronoi tessellations used in a phase-field simulation to simulate the recrystallization during a cold rolling process. Some applications of PV structures can be found in [16]. Based on closed-pack non-overlapping circles, Jafari and Kazeminezhad [17] proposed a new Voronoi diagram in a Laguerre geometry. In their microstructures, grain sizes and fractions were created by controlling the size and distribution. Based on Laguerre Voronoi tessellations, the generation of 3D numerical models for polycrystalline microstructures was given in [14]. Falco et al. [14] used the 2D sections (by cutting 3D tessellations) to compare with real 2D surfaces, in order to achieve the most representative set of input values. Based on the minimization of the energy and distance function, a process was presented to dynamically generate the spherical Laguerre Voronoi diagram based on assumptions in the real world [18]. Such Voronoi diagram has the corresponding properties of polyhedra as described in [19]. i.e. the spherical Laguerre Voronoi diagram is valuable to analyze or model some real-world phenomena. The relative density and irregularity of Voronoi closed-cell foam structures can be used to investigate the elastic characters of the foam material [20]. With a combination of the Laguerre Voronoi tessellation method and the FE method for solving the heat transfer problem, the thermal properties were computed [21]. Luo et al. studied the optimization of the centroidal Voronoi tessellations by using a Laplacian operator, where a density function was used to control the size and distribution of Voronoi tessellations [22]. Their study results can be used in dynamic simulations of forming processes. Besides in materials science, Voronoi tessellations has also applicability in other fields, e.g. in aviation, in biology, in ecology, in geography and in astronomy as well as in telecommunications. Using the cohesive element-based numerical manifold method with Voronoi cells (grains), a fully coupled hydro-mechanical formulation is presented for the investigation of the hydraulic fracturing of rock on the micro-scale [23]. To study the rock strength, a methodology for generating Voronoi tessellations presenting a polycrystalline rock microstructure is proposed to simulate the packing processes and the crystal growth of mineral grains, where the distribution of the poly-crystals (tessellations) varies in shapes and sizes representing different mineral grains [24]. In presenting structures for individual cases, especially for detailed structural features, Voronoi tessellations have more than one drawback, if no extra constraints for cells are set. Using the Voronoi estimator for the intensity estimation, the major drawback is that it tends to paradoxically under-smooth (the over accentuating of local intensity) the data in regions with a high point density and to over-smooth in regions with low point density. To remedy this behavior, an additional smoothing operation to the Voronoi estimator is applied based on resampling the point pattern by independent random thinning [25]. Their proposed intensity estimation schemes were used to two datasets: locations of pine saplings (planar point pattern) and motor vehicle traffic accidents (linear network point pattern) [25]. Mentioned in Koufos and Dettmann [26], little is known about the properties of Voronoi cells located close to the boundaries of a compact domain. In a domain with boundaries, Voronoi cells would be naturally clipped by the boundary. They [26] used a low-complexity numerical method to compute the mean cell area of homogeneous Poisson Voronoi tessellations for seeds located along and/or close to the boundary of a quadrant. It is found that, by using the location-dependent parameters calculated by using the method of moments (second moment of the cell area), it can provide a reasonably good approximation of the wished structure, i.e. avoid the above mentioned cell area problem near to boundaries. Voronoi tessellations are ubiquitous in nature, since they possess the ability to model growth processes and equilibrium states [27]. To be simple and time efficient, users can firstly try to use softwares with open accesses, e.g. "NEPER" [28]. "NEPER" cannot fulfill some individual requirements for microstructure topologies in our cases, e.g. 3-twins in Fig. 7d-f and mixed crystals with layers in Fig. 8g. It is well known that PV tessellations result in relatively homogeneous cell sizes, which can be a defect to present structures with inhomogenous sizes. The current work presents a simple process to numerically generate 2D and 3D hierarchical PV structures, i.e. to remedy the aforementioned defect. A "hierarchical structure" is achieved by controlling the number of total nuclei in a given volume or subvolume. The size distribution of a given phase can be described by dividing the whole size spectrum into N subzones. Each subzone has its own mean size, i.e. the concept of an "artificial phase" (each subzone corresponds to an artificial phase). A further contribution of this work is to illustrate processes in order to fulfill individual microstructural requirements, e.g. twins. The method is to set constrains during the nuclei generation. In the 3D case, the tomography concept is used to present the 3D PV tessellations. Since it is only necessary to generate points (Voronoi nuclei) and write out numbers, it needs only some basic programming knowledge of any software. No deep understanding of mathematical theories is required. PV microstructures with ASCII format are directly ready for meshing. Here, PV tessellations are applied as artificial microstructures in FE simulations to investigate solid material behaviors.

A simple numerical method for the generation of Poisson Voronoi tessellations
As mentioned, the current work is aimed at presenting a simple numerical way to generate PV tessellations for material simulations. For a detailed knowledge about Poisson field and Voronoi cells, it can be referred to [3][4][5]. Here our method will be directly introduced. Tessellations in PV structures should match characters of material microstructures.

Basic concepts
where · indicates the distance between two points in E. In our case, {x i } points are generated by the subdivision of E in a regular way (equal subvolumes S x i with ∪S x i = E), and {x i } points lie in the center of each S x i . {x i } points are uniformly distributed, i.e. a homogeneous Poisson's distribution. The above generation process is also suitable for 2D cases, i.e. within a controlled area. Figure 1 schematically presents the growth process of 2D Voronoi nuclei (generation points). As a result of the aforementioned growth, Voronoi cells/tessellations are created. In a 2D controlled area, 10 Voronoi nuclei given as black points are randomly generated (Fig. 1a). During the growth, the speed is identical in both the x-and the y-direction. In Fig. 1b, each circle (including dashed and solid lines) has a Voronoi nuclei as its center point. The two Voronoi nuclei with dashed circles will firstly meet each other, since the distance between these two neighboring nuclei is the minimum one among all. Fig. 1c denotes the final result of such a growing process of Voronoi nuclei, i.e. Voronoi tessellations are created. More simply, the aforementioned Voronoi tessellation structure can be obtained by: (I) supposing each Voronoi nuclei has M neighbors, drawing the middle line between each two neighboring Voronoi nuclei, (II) connecting intersection points of these middle lines in sequence; (III) an individual Voronoi tessellation being a polygon with M edges. Figure 2 schematically illustrates the generation process of 2D PV tessellations.   (Fig. 2b). After searching for the first layer, it goes to the next neighboring layer (second layer). The search still begins with the left subzone and so on. Figure 2d means that 14 subzones of the Poisson field (pixels with inclined lines) belong to a certain Voronoi cell after the whole searching process. The area fraction of this tessellation with 14% ( 14 10×10 ) is slightly higher than the mean value 12.5%. Each tessellation will be assigned an individual pixel color number for the graphical presentation. Figure 2 is suitable for one phase material, which can be polycrystalline. Figure 3 schematically gives the generation process of a two-phase material with different phase mean sizes. This PV tessellation generation process (Fig. 3) can be analogously applied for multi-phase (N-phase) materials with non-identical mean sizes of phases. A detailed description to generate Fig. 3 is referred to the appendix. Here, the general method to create a hierarchical PV structure with N phases is given in Table 1 for the 2D case.
For the generation of 3D PV tessellations, the process is analogous to the one for the 2D case. Some precal- } are identical. Concerning the Poisson field, one can imagine that the volume E is regularly divided into layers in the z-direction. In Fig. 4b, the volume E is divided into 4 equal layers by the dashed lines. 4 images go through center points of the above mentioned equal layers and they are perpendicular to the z-axis. These 4 images compose an artificial tomography. Each above mentioned image will be subdivided into n x × n y pixels. It is user dependent how fine the net is. The refinement degree of the net depends on the number of the total layers and the total pixels of each image. In Fig. 4b, each image in the artificial tomography is subdivided into n x × n y = 5 × 5 pixels. These pixels are presented by squares in blue. "Squares" indicates that the length of a pixel is identical in the x and in the y direction (Fig. 4b, c). A small blue rectangle point locates at the center point of each pixel. If the length of each layer in z direction, the length of each pixel in the x and the y direction are the same (Fig. 4b, c), then the controlled volume is regularly subdi- Step Description Calculation of two types of nuclei numbers for each phase, i.e. P j−total and P j−real (generator points) with j ∈ [1, N ] is only an intermediate tool to achieve a geometrical hierarchical structure, where P j−real ≤ P j−total and P j−total ≤ P total as well as { p For each phase j with j ∈ [1, N ], assignation of remaining {x i } to { p j−total }. The condition for the assignment is After the assignment for the current phase, then j = j + 1 and the data should be updated. Then, it starts the assignment for the next phase; 5 At last, a hierarchical PV structure with N phases will be created with ∪S spectrum is [0, 255]. A colored ASCII format is out of the discussion range of the current work. For a multi-phase microstructure, phases can have non-overlapping subranges of the above mentioned color spectrum. If a microstructure includes more than 256 grains and particles (more than 256 pixel groups), our method is still applicable. Strictly, an extra variable should be used in the code to trace each pixel group. The usage of this variable is to guarantee that no neighboring groups will be assigned an identical pixel color. It is very rare that a grain or a particle cluster has more than 255 neighbors in the engineering applications. Practically, the code can only assign each pixel group a color number one after another. The possibility is very low for neighboring pixel groups to have the same pixel color, since the cells are randomly generated in the control volume and there are large enough data to guarantee a representative microstructure. Besides, each phase has its own color number range. It means that pixel groups belonging to different phases will have non-identical pixel color numbers. Even though two or a few neighboring pixel groups have the identical pixel color, the percent ratio is still near to zero among a few hundreds of PV cells (total number of PV cells at least larger than 256). In FE simulations of polycrystalline materials, the maximum number of considered cells are within 1000, e.g. [1,8,9].

Generation of 2D and 3D multi-phase PV microstructures
In the generation process of the PV structure, the Voronoi nuclei are firstly randomly generated in the range of [0, 1] (ABC O in Fig. 5a). Then, coordinates of nuclei are scaled into the desired dimensions ( Fig. 5a). To keep the same nuclei growing rate in all the directions, the above mentioned desired dimensions should be a square in 2D cases. To create a rectangle-shaped PV microstructure, the involved area should be firstly selected as a square. The edge length of this square is the same as the longer edge length of an arbitrary size distribution of a phase divide into "m" subzones ; b division of the whole size spectrum into "m" subzones to realize a better size distribution in a PV microstructure, where each subzone is taken as an extra artificial phase and has its own mean size and volume fraction the rectangle. Then, the PV structure should be cut into the desired size. An example is presented in Fig. 8g, h. In 3D cases, the controlled volume should be a cube. In usual, a PV tessellation diagram gives a relative homogeneous size distribution. This can be a defect for the representation of real microstructures. A simple way, i.e. adding artificial phases, can be used to compensate this defect. If a phase (phase-A) possesses a relative large size spectrum in reality, "m" (m > 1) phases in PV structure can be used to describe phase-A (Fig. 5b). Each artificial phase would have its own mean size, but the mean size of "m" phases is the same as the one of phase-A, where the volume weighting factor should be considered according to the phase-A size distribution. It is pointed out that the concept of an "artificial phase" is a way to improve the shortage of a relative homogeneous cell size distribution. But this concept increases the total amount of constraints for the nuclei, which makes the finding of all nuclei coordinates more and more difficult. A good mapping of a very irregular real size distribution of a phase might be impossible by using this concept, since the generation of all desired nuclei might be too difficult to be successful. Figure 5b schematically describes the way to present the size distribution spectrum of a phase by using PV tessellations. Each subzone in Fig. 5b will have its own mean size. As a final result, the size distribution can be better described. For a given phase in a PV structure, a volume deviation exists between the phase volume fraction in the final structure and the prescribed volume fraction as an input parameter. The larger the number of total nuclei, the less the volume deviation. Statistically, a larger number of total nuclei depicts usually a more representative microstructure. Generally, some basic physical data, like the volume fraction, mean size and size spectrum of each phase as well as cell/grain/particle distributions, should be obtained from each individual material and applied in the PV structure generation. To obtain a representative structure, the total nuclei number should be big enough for each phase. The minimum distance between nuclei must be considered to avoid unrealistic tessellation concentrations, e.g. islands. If islands exist, some phases would be completely surrounded by the other phases. In our code, an algorithm according to Ahrens et al. [29] is used to generate pseudo-random numbers in the range of (0, 1). By using this multiplicator [29], the advantage is that the nuclei distribution quality is high, i.e. there are no islands observed in our PV structures. A minimum distance control helps to ensure that nuclei for each phase are distributed in the whole region, i.e. no complete concentration of all nuclei of a phase in a small subregion of the whole image exists. This control also avoids that nuclei of secondary phases (small volume fraction and/or less nuclei) would solely locate on boundaries of the major phase. In the current work, there are more than 100 nuclei for each phase in all the presented PV structures (not the one to show geometrical adaptive meshing). In combination with the applied nuclei distance control, there exist neither islands nor the phase concentration on the boundary of another phase. After each generation, the resulted geometry should be reviewed critically. If there is anything far away from the reality, more constraints should be set for the generation of nuclei.

2D two-phase polycrystalline microstructures
For a 2D image, a pixel length can be taken as a unit one. This unit length is identical in the x and in the y direction (coordinate axes see Fig. 2). Here, a polycrystalline two-phase Ag/SnO 2 oxide dispersion strengthend (ODS) composite is given as an example. A controlled area with 500 × 500 units is supposed, i.e. an image with 500 × 500 pixels. If absolute dimensions are necessary in FE simulations, node coordinates can be scaled into desired dimensions after meshing. Based on the experimental data in Table 2 [30], the mean size of SnO 2 particles is smaller than the one of the Ag phase, i.e.
In generating 2D PV microstructures, the volume fraction is also taken as the area fraction for 2D cases ( }. It is worthy to mention that all the data for the generation of PV tessellations come from material design or physically measured data. Based on the statistical data (Table 2), several material input parameter sets are used to generate the PV artificial microstructures. As an example, a PV structure with hierarchical tessellations is shown in Fig. 6b. It has 800 SnO 2 particles and 150 Ag grains. Figure 6a presents a real microstructure cut-out obtained from an electron backscatter diffraction (EBSD) test. In comparison with Fig. 6a, b, it is obvious that the mean size ratio of SnO 2 to Ag in Fig. 6b is smaller than the one in Fig. 6a. The reason is that the calculation of the statistical data in Table 2 is based on a much larger size than the size given in Fig. 6a. If an artificial microstructure should be more matchable to the real cut-out (Fig. 6a), one can enlarge the mean size of the particle phase by keeping the Ag mean size constant. Compared to Fig. 6b, c is taken as a better mapping for the real microstructure. A particle cluster marked by a red circle in Fig. 6a shows a character of the real microstructure, i.e. some small particles compose a cluster and are located in Ag grains. In Fig. 6a, a large particle surrounded with a purple circle locates among Ag grains. Such a particle can be taken as a neighbor of Ag grains. Particles in the yellow circle (Fig. 6a) are partly inside certain Ag grains and partly along Ag grain boundaries. The PV tessellations in Fig. 6c can catch the above mentioned three microstructure characteristics marked in corresponding colored circles as in Fig. 6a.
In Fig. 6a, some very small particles are sparsely distributed in the whole structure. If such small particles should be embodied in the PV structure, a third phase (an artificial phase) can be added. It means that the particle phase is described by phase 1 and 2 with a mean diameter of d phase−1 and d phase−2 , respectively. Since the above mentioned artificial phase has the smallest mean size, it should be the phase 1. This means that its mean size d phase−1 is used to calculate the { p Table 1). The Ag phase will be described by phase 3 with d Ag = d phase−3 . Figure 6d possesses 3% phase 1, 14% phase 2 and 83% phase 3 with 1 2 d phase−1 : 1 2 d phase−2 : 1 2 d phase−3 = 2:5:10. It is worth mentioning that the very small particles in Fig. 6d will lead to more elements in meshing. In such a microstructure with many small particles, the total number of elements may exceed the calculation capacity limit in 3D cases.
Our Ag/SnO 2 ODS is heat treated, i.e. two passes of the hot extrusion [30]. After hot extrusion, the Ag phase has about 20-25 vol.% 3-twins. Figure 7a is identical to Fig. 6a. Figure 7b shows the twin boundaries of Fig. 7a. Figure 7c schematically shows a simple way to generate nuclei pairs for twins. In Fig. 7c, point P 1 (a, b) in black is randomly generated. Points P 1 (a, b) and P 2 (x, y) (in red) compose a pair of twins. In the current work, point P 1 (a, b) and point P 2 (x, y) are called a master nuclei and a slave nuclei, respectively. The length c 1 denotes the distance between the twin boundary AB and the master nuclei. The distance between the twin boundary and the slave nuclei is c 2 . The angle between the xaxis and the twin boundary is presented as α. For the sake of simplicity, it is taken as c 1 = c 2 = c in Fig. 7d-f. For a given distance c and an angle α, the coordinates of P 2 (x, y) can be determined. It also means that the slave nuclei depends on the master nuclei. In Fig. 7d-f, master nuclei generate the black grains. Grains in grey color are from slave nuclei. The area fraction of twins is about 20% in Fig. 7d-f.  Fig. 7f. The value of α and c can be statistically obtained from experimental data. If the boundary length also should be taken into consideration, more constraints should be set to generate master and slave nuclei.

Overlapping: a simple way to create multi-phase microstructures
Following the generation process as given in Fig. 6b } with j ∈ [1, N ] need to be updated for each phase. The advantage of this process is that the volume fraction of each phase is guaranteed. For an N-phase material, a simpler way would be the overlapping: (I) 1-N phases are arranged according to their vol.%, i.e. vol.% phase−1 ≤ vol.% phase−2 · · · ≤ vol.% phase−N ; (II) N images are generated, i.e. each phase one image; (II) the first image (for phase 1) overlaps the second image. The overlapping resulting image overlaps the third one and so on. As a final result, an image with N-phases will be created. In this way, the deviation of the volume fraction before and after overlapping is kept possibly small for phases with small volume fractions. The advantage of the overlapping process is that no calculation of { p phase− j−total i } is necessary with j ∈ [1, N ]. The disadvantage is that the phase volume fraction is only guaranteed for two-phase (N = 2) materials, but not for N ≥ 3. Practically, even though N ≥ 3, the final image for given a, b, c: The PV tessellations as artificial microstructures are also widely used for ceramic materials. Zirconia ceramic ZrO 2 is a kind of a brittle technical ceramic for industrial and biomedical applications. It possesses many favorable properties, like the high strength, high toughness, very good thermal insulation and shape memory ability [31,32]. Figure 8a presents a microstructure of such a material [33]. Such a ceramic is polycrystalline and belongs to the high-performance ceramics. For some usages, subregions with different phases might appear inside one microstructure cut-out, e.g. a layered or coated material. Figure 8b presents a zirconia-alumina ceramic composite, where zirconia is shown in light grey color [34]. Such zirconia-alumina ceramic composite also can be manufactured as layered material and coated materials. Figure 8c illustrates a small microstructure cut-out from a SEM image of a polished surface from a ZrO 2 + 30 wt.%Al 2 O 3 composite [35]. Zone A 1 in Fig. 8c located near the interface between the zirconia phase and the alumina phase, which presented the grain growth (in alumina phase). Zone B 1 in Fig. 8c referred to the formation of a protective alumina layer on the zirconia-rich surface. The interface in Fig. 8c can be approximately taken as vertically oriented. Figure 8d presents a composite with the zirconia as a coating material, where the interface is horizontally oriented [36]. It is pointed out that the current work concentrates on the presentation of microstructures by PV tessellations. If the interface is non-linear in a larger scale (larger than the corresponding microstructure's scale), such a case is out of the topic of this work. The characteristics illustrated by Fig. 8ad can be summarized as: (I) polycrystalline composites; (II) non-uniform mean sizes of different phases; (III) approximately linear interface on the given microscales; (IV) the orientation of the interface varies. In order to create a PV microstructure with the aforementioned (I)-(IV) characteristics, Fig. 8e, f schematically shows the approach. In Fig. 8e, the microstructure is divided into 2 subzones by the straight line CD, which has an angle α to the x-axis (Fig. 8f). The above mentioned two subzones correspond to subzone A and B in Fig. 8f. Both subzones are supposed to be polycrystalline. Subzone A has only one phase (phase-1). Subzone B is composed of two different phases (phase-1 and phase-2), each possessing its own mean size. In subzone B, the mean size ratio between phase-1 and phase-2 is taken as 1:5 and the area ratio 1:1. It is simple to create such a PV microstruc- PV artificial microstructure applied for ceramic layered/coated materials: a a measured microstructure cut-out of a polycrystalline zirconium dioxide (ZrO 2 ) ceramics [33]; b zirconia (light grey color) with alumina (dark grey color) [34]; c zirconia with a protective alumina layer (zone A 1 : alumina grain growth, zone B 1 : alumina layer in zirconia-rich zone) [35]; d zirconia as coating material with a horizontally oriented interface (EB-PVD: electron beam physical vapor deposition) [36]; e a sketch for a layered material with the interface orientation in a general case, where the microstructure includes multiphase polycrystals in subregions and different phase mean sizes; f a condition to fulfill the requirements of the individual phase distribution given in e; g a PV microstructure for polycrystalline layered materials with different phases and phase mean sizes in subzones; h if necessary, cutting the image into desired dimensions ture as prescribed by Fig. 8e, f: (I) define two points on edges of the controlled area to get two regions A and B in Fig. 8f. The aforementioned two points give a straight line with a slope α; (II) for a randomly generated Voronoi nuclei with coordinates (P x , P y ), the y-coordinate P α y can be calculated (Fig. 8f). If P α y > P y , it is a Voronoi nuclei in region B. Otherwise, in region A. By using the overlapping method, Fig. 8g shows a PV two-phase layered hierarchical microstructure with subzones. It includes two polycrystalline phases and fulfills individual requirements (characteristics shown in Fig. 8a-d) and assumptions given in Fig. 8e, f) for subregions. Since there are only two phases (Fig. 8g), the area fraction remains for each phase, i.e. not changed by the overlapping. If a rectangle shape is necessary, as an example, Fig. 8h is the cutting result of Fig. 8g to fulfill the size requirement of the whole structure.

Generation of 3D microstructures
Like aforementioned, the concept of the tomography is used for generating 3D PV microstructures. The procedure is the same as for 2D ones. Simply, the images of the tomography is created one after another in the third direction. Due to the identical cell growing rate in all directions, the 3D controlled volume should be a cube for PV tessellations. It implies that: (I) a pixel (unit one) has a square shape; (II) the distance of two neighboring images equals the edge length of a pixel; (III) a voxel is a cube. For a PV structure, if a cuboid is required (edge length not equal), a cube with the longer edge length should be firstly created. Then, it is cut into the desired cuboid. 3D PV microstructures will be shown together with meshing in Sect. 2.3.

An example of a two-phase polycrystalline material: Ag/SnO 2 ODS
The meshing can be done by softwares available for all, e.g. the free software "OOF2D" and "OOF3D" with open accesses [37,38]. Software "NEPER" [28] with an open access also has the meshing function. In the current work, the commercial software "Simpleware Scan IP" [39] is used. This software provides much more functions than "OOF3D" and meshes only 3D structures. If a 2D meshing is necessary, it can be easily created in 3 steps: (I) a 3D structure can be created by copying the given 2D image in manifolds (Fig. 9a); (II) 3D meshing is generated (Fig. 9b); (III) 2D meshing can be obtained by extracting the nodes and elements on the z = 0 surface. Figure 9c is identical to a single image in Fig. 9a, but in color. As a result, a geometrical adaptive meshing is done for a polycrystalline 2D microstructure. The 2D geometrical adaptive meshing marked in a rectangle area ABCD (Fig. 9c) is presented in Fig. 9d. It is worthy to mention that a 2D meshing can be directly applied in axisymmetric analyses ("CAX" element type in ABAQUS). Based on the statistical data in Table 2 Figure 10a shows 3 artificial tomographic images. Practically, a volume reduction occurs for the particle phase from the tomographic images (voxel vol.%) to the final FE meshing (elements vol.%). In order to fulfill the meshing conditions, e.g. no overlapping space among cells and no pores in the controlled volume E, a closed surface needs to be created for each tessellation. During this surface creation, "shrinkage" occurs. As a result, the volumes of such tessellations with "shrunk" closed surfaces are reduced. Such a volume reduction occurs also in other meshings, e.g. the one for an Al/Al 2 O 3 real microstructure (commercial software AVIZO) [40]. Some trial meshings are performed for PV structures with an enlarged volume fraction of the SnO 2 phase. Figure 10c illustrates the geometrical adaptive 3D meshing for a microstructure with about 16.2vol.% SnO 2 particles ("EVOL" in ABAQUS). The basic characteristics of the particle geometry and distribution can be captured, e.g. the non-uniform shape and size of particles. Some small spherical SnO 2 particles locate in the larger Ag grain, like SnO 2 -A in Fig. 10c. Some particle clusters can be taken as neighbors of Ag grains, like the particle SnO 2 -B in Fig. 10c. Still some clusters locate among Ag grains as well as go through/into Ag grains, like the particle SnO 2 -C in Fig. 10c. In the example, the reduced volume of the particle phase is integrated to the matrix phase.

An example of a three-phase composite: Co/WC/diamond
To investigate the deformation and damage behavior of Co/WC/diamond metal matrix composite (MMC), 2D and 3D real microstructures are applied. Figure 11a with 54 diamonds illustrates a 2D real microstructure [41]. After the pixel handling, the image is given in Fig. 11b, which also shows a part of the meshing along the particle boundaries. After the meshing, the WC and diamond phase have an area fraction of 5.23% and 8.64%, respectively. For the 3D real microstructure cut-out, some 2D images from a tomography and the final geometrical adaptive meshing is presented in Fig. 11c. There are 5 diamond particles, two of which are located on the boundaries of the cut-out and much bigger than the other three. This means that the 3D real structure (Fig. 11c) is less representative than the 2D one (Fig. 11b).
There is no possibility to create a real 3D microstructure with about 50 diamond particles, since the experimental tomography covers a much smaller size than the 2D image (Fig. 11a). Comparing FE results of these two models, some conflicts arose for the phase stress. There is no testing data for the phase stress. To solve this conflict, an artificial 3D PV microstructure is preferred. It should have 54 diamond particles. After only two trial meshings (with different WC vol.%), a suitable structure is found and shown in Fig. 12.
The element volume fractions of the WC and diamond phase are kept the same as in the 2D real microstructure (Fig. 11b).
Here, "the same" volume fraction means within a tolerable deviation. Figure 12a presents a stack view of some selected artificial tomographic images. The whole size is 50×50×50, in which 3580 WC particles are included. The mean size ratio of the WC and the diamond particles is 1:3.6. The mean size ratio in the 2D real microstructure (Fig. 11b) is taken into consideration for the creation of Fig. 12. It is pointed out that the statistical mean size of the WC phase cannot be directly used here, since WC particles should be partly invisible to obtain a good matchable 3D PV structure with the real 2D one (Fig. 11b). In the FE simulation, the matrix has homogenized properties from the pure Co and the invisible WC particles. Figure 12b is an overview of the geometrical adaptive meshing. Figure 12c, d present the diamond (8.00 vol.%) and the WC phase (5.89 vol.%) after the meshing, respectively.
To generate the tomographic images like given in Figs. 10a and 11a, it took only about 1 min. It takes about 40 min from the creation of polycrystalline PV images Fig. 10a to the final applicable meshing Fig. 10c. From Fig. 11a to the final applicable meshing Fig. 11b, it needs only a few minutes, since the manual work of the pixel group assortment is much less. Fig. 11 a A 2D real microstructure cut-out presented by a SEM image [41]; b the corresponding image of a after the pixel selection, where the regular mesh is created by a FORTRAN file due to a large amount of small WC particles, after which the meshing is improved by the software Hypermesh for some sharp corners (WC: 5.23%; diamond: 8.64%); c measured tomographic images, selected regions and the meshing for the FE simulation

Boundary conditions
Artificial boundary effects exist for the applied PV structure, which is introduced by tessellations (materials) outside the resulting geometry. A periodic PV structure can reduce these artificial boundary effects, which requires periodic boundary conditions (BCs). For micromechanical FE simulations, BCs are one point which must be considered. BCs given as constraints for node degree of freedoms are always a disturbance for the calculated material behavior, i.e. it causes, more or less, inaccuracy of the predicted results. There are usually two types of BCs, i.e. homogeneous and periodic ones. Homogeneous BCs set too strong constraints for the degree of freedoms of boundary nodes, where nodes will move parallel with each other, including the zero displacement case. Comparatively, periodic BCs soften these strict constraints, which allow non-uniform displacements of boundary nodes and in a periodic formulation for counterparts of edges/surfaces. In many cases, it is a major task to calculate the stress-strain relations for FE simulations. If so, constraints for the other (stress or strain) should be also considered, when constraints are set to one of them. Taking strain controlled loading as an example, periodic BCs require periodic structures and antiperiodic conditions for corresponding stresses, e.g. as shown in [8,9]. Our code to generate hierarchical PV microstructures will be extended to cover the case for periodic microstructures in our consecutive works. Some of our previous works used Voronoi (not hierarchical) structures with periodicity x y block-2 block-3 block-4 block-8 block-9 A B D C Fig. 13 A sketch for the generation of periodic Voronoi nuclei: (i) randomly generate the nuclei in "block-1" (range ABCD); (ii) copy the "block-1" nuclei to "block-2"-"block-9", i.e. "block-1"-"block-9" with the identical nuclei distribution; (iii) "block-1" will be an applicable periodic Voronoi structure [9,42]. Figure 13 illustrates a sketch for the generation of 2D periodic Voronoi structure. In the 2D case, 3 2 = 9 samearea-sized blocks of PV structures (Fig. 13), which with the identical Voronoi nuclei distribution are neighbors, can result in an applicable periodic Voronoi structure for the FE simulation, i.e. the center block. For the 3D case, it is analogous, i.e. the above mentioned blocks should be 3 3 = 27. Based on experiences of the current authors, there is a way to apply non-periodic microstructures and without the strict constraints of homogeneous BCs applied on this microstructure, i.e. simultaneous micro-macro two-scale FE simulations. In our case, BCs are applied on the macrostructure, which means the microstructure is free of BCs. This two-scale simulation method is also suitable for the application of real microstructures, since, in reality and for most cases, there exists no periodicity for real microstructure cutouts. Figure 10c for the Ag/SnO 2 ODS is used to prove the correctness of a simulation process, which is found by axisymmetric FE simulations. 2D and 3D PV microstructures will also be applied to investigate the influence of 3-twins on the texture evolution.

Further works
Besides the application mentioned in Sect. 2.3.2, PV tessellations will also be used to simulate the crack path, the debonding and the damage evolution for the Co/WC/diamond MMC. Since the generation of a PV structure is time efficient, it is preferable to make parametric studies and to optimize the microstructure. FE results predicted by PV structures will be presented in consecutive reports later on.

Conclusion and outlook
In the numerical investigation of material behaviors, the application of microstructures in one-scale or multi-scale is getting more and more essential. Real microstructures may be unavailable, especially 3D polycrystalline ones. Still, it can be difficult to achieve a well representative microstructure in a limited small cut-out of the real microstructure. This work presents a simple numerical method to generate multi-phase and polycrystalline artificial microstructures by using hierarchical Poisson Voronoi tessellations. By using softwares with an open access, e.g. "Neper" [28], some individual microstructural characters cannot be described. The basic concept of our method is to create images with an ASCII format. For 3D PV microstructures, "artificial tomographic images" are used, i.e. more than one 2D images aligned in order in the third direction. In this work, four examples are presented for applications of artificial 2D/3D PV microstructures. The meaning of our method is: • to fulfill individual requirements from material designs by setting constraints for the Voronoi nuclei generation, e.g. twins and mixed phases in selected regions; • to make it easy by simply writing out numbers (pixel color) for pictures and applicable for engineers with only some basic programming knowledge of any software; • to obtain a PV image with a better match to the real one, e.g. a better description of the size distribution by using the concept of artificial phases; • to have potentials for further studies, e.g. generation of 3D twins for the Ag phase in Ag/SnO 2 ODS in our example.

Appendices
To describe the generation process of Fig. 3, the content listed in Table 1 is considered. For Fig. 3, the assumption is that a two-phase material has phase-1 and phase-2. It has a composition of 20% area fraction for phase-1 and 80% area fraction for phase-2. The mean size ratio (diameter in 2D) between "d 1 " for phase-1 and "d 2 " for phase-2 is d 1 : d 2 = 1 : 2, i.e. (d 1 ) 2 : (d 2 ) 2 = 1 : 4. It implies that PV tessellations should be hierarchical. To create an artificial microstructure for such a material composition, supposing the quadratic microstructure with an edge length D 0 , it is as follows: • Calculation of the Voronoi nuclei numbers for each phase: -by considering the mean size and the volume fraction of each phase, the cell number of each phase can be fixed, if there are no other influential factors: • in our example (Fig. 3), the cell number for phase-1 is P 1−real = 20% × D 2 0 1 4 π(d 1 ) 2 and the one for phase-2 is P 2−real = 80% × • the ratio of the cell number between phase-1 and phase-2 is P 1−real : P 2−real = 1 : 1 ; • supposing that the mean size of phase-1 tessellations is 5% as an area fraction, then this value for phase-2 will be (d 2 ) 2 (d 1 ) 2 × 5% = 4 × 5% = 20%; the number of Voronoi tessellations for phase- • Calculation of the number of total Voronoi nuclei to achieve hierarchical tessellations, i.e. the maximum value P total for { p total i }: -by taking the phase with the smallest mean size (d small ) among all phases, the value of P total can be fixed. In 3D case, P total = vol total 4 3 π( 1 2 d small ) 3 . • for the 2D case as given in Fig. 3, it is phase-1 with a mean size d 1 . Since the mean size of phase-1 tessellations has 5% area fraction, P total = 100% 5% = 20; • the above 20 Voronoi nuclei (P total ) include the 4 Voronoi cells for phase-1 in { p 1−real i } and the 4 for phase-2 in { p 2−real i } (Fig. 3a); • not all the P total = 20 Voronoi nuclei present a tessellation in the final structure. Only the 8 nuclei (real ones) for phase-1 and phase-2 will result in 8 tessellations in the final image. The other 12 nuclei are used only for the creation of different phase mean sizes; • Assignment of each subzone in Poisson field to a certain Voronoi nuclei: -arranging phases according to their mean sizes d, i.e. d 1 ≤ d 2 ≤ · · · ≤ d N ; -beginning with the phase with the smallest mean size (phase-1, n=1); • for example in Fig. 3, it begins with phase-1, and the minimum distance search goes through all the 10 × 10 = 100 Poisson's subzones and the 20 Voronoi nuclei ({ p 1−total i }, i ∈ [1,20]); • after the search, 20 pixels with inclined lines (Fig. 3b) are assigned to 4 Voronoi nuclei of phase-1 (in { p 1−real i }). It implies that the mean size of phase-1 tessellations is 20 4 = 5 (pixels); • in the final images, pixel groups present tessellations and each pixel group will be assigned a pixel color number. The grey scale for the pixel color is applied here. It means that the pixel color number is in the range of [0, 255] for the ASCII format; -going to the next phase with the second smallest mean size (phase-2, n = 2); • it is phase-2 in Fig. 3. The search goes through the remaining 100 − 20 = 80 subzones in the Poisson field (pixels without inclined lines) and the 4 Voronoi points for phase-2. It means that the mean size of phase-2 tessellations is 80 4 = 20 (pixels). It is to mention that i ∈ [1,4] is identical for { p 2−total i } and { p 2−real i }; • in the PV diagram, the mean size ratio between phase-1 and phase-2 is 5% : 20% = 1 : 4 = (d 1 ) 2 : (d 2 ) 2 ; -Going to the next phase (phase-3 upto phase-n), if necessary; • for the next phase (phase-n), the number of subzones in the Poisson field is the number of remaining pixels. Here, "remaining" pixels refers to the ones not yet assigned to any Voronoi tessellation • for a representative PV microstructure, the number of the Voronoi nuclei and of the subzones in the Poisson field needs to be large enough.