Image-based semi-multiscale finite element analysis using elastic subdomain homogenization

In this paper we present a semi-multiscale methodology, where a micrograph is split into multiple independent numerical model subdomains. The purpose of this approach is to enable a controlled reduction in model fidelity at the microscale, while providing more detailed material data for component level- or more advanced finite element models. The effective anisotropic elastic properties of each subdomain are computed using periodic boundary conditions, and are subsequently mapped back to a reduced mesh of the original micrograph. Alternatively, effective isotropic properties are generated using a semi-analytical method, based on averaged Hashin–Shtrikman bounds with fractions determined via pixel summation. The chosen discretization strategy (pixelwise or partially smoothed) is shown to introduce an uncertainty in effective properties lower than 2% for the edge-case of a finite plate containing a circular hole. The methodology is applied to a aluminium alloy micrograph. It is shown that the number of elements in the aluminium model can be reduced by 99.89%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$99.89\%$$\end{document} while not deviating from the reference model effective material properties by more than 0.65%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.65\%$$\end{document}, while also retaining some of the characteristics of the stress-field. The computational time of the semi-analytical method is shown to be several orders of magnitude lower than the numerical one.


Introduction
As the computational power of modern hardware continues to increase, multi-scale modeling of materials inevitably becomes more appealing even in engineering problems of industrial scale. The academic literature on such methods is massive and in some sense also multi-disciplinary, combining efforts from fields like model reduction, high performance computing and machine learning with material modelling, just to name a few [1][2][3]. Many books are available on related topics such as multiscale modelling, homogenization methods, and micromechanics related to heterogeneous materials [4][5][6][7][8]. Of particular interest is the area of computational homogenization, where knowledge about the microscopic constituents are used in order to predict material properties at the scale of engineering products [9]. In this paper, we focus on the specific task of estimating effective mechanical properties of heterogeneous engineering materials via computational homogenization, using an uncoupled multiscale approach which we denote semimultiscale. A real, or artificial micrograph is used as input in this methodology, which is split into multiple subdomains that are evaluated independently from each other. The effective material properties are then mapped back to a reduced mesh, decreasing the level of detail in the model. This approach generates a reduced model for further analys, but also allows the evaluation of models which would otherwise be impossible to run on one machine to be evaluated in separate parts sequentially, or on multiple threads/machines in parallel.
A large amount of time and effort is spent both in academia and industry studying the microstructure of manufactured materials. In this process, micrographs are typically generated which are used in everything from microstructural characterization, process-structure-property modelling, defect identification and failure analysis to building numerical models [10,11]. Such numerical models inherently contain multiple layers of approximations. Firstly, the micrograph itself is an approximation of the real physical microstructure. Secondly, image processing operations such as segmentation only approximates the actual material distribution of individual phases. Thirdly, discretization of the segmented material description introduces a more or less jagged approximation of phase boundaries, depending on which discretization strategy that is used [10,12]. Lastly, such models might also need to be reduced depending on the available computational resources. In this paper we deal with model reduction by splitting the micrograph into multiple subdomains which are evaluated independently, meaning that regions of the same phase which were originally connected might become disconnected, and separated regions which are actually interacting might become unaware of each other. By generating artificial micrographs for shapes which have known analytical solutions for their constitutive behavior and stress under given loading conditions, the discrepancies from the first two approximations can be eliminated, and the errors introduced by discretization and model subdomain reduction can be isolated and quantified.
The paper is structured as follows; a short theoretical background on analytical expressions for effective mechanical properties, and analytical stress concentrations are given in Sect. 2. Section 3 covers implementation details, as well as the methodology used to investigate the influence of both discretization strategy and subdomain splitting. Section 4 presents results and analysis from applying the proposed subdomain methodology to the edge case of a finite plate containing a circular hole, as well as a more realistic case of a complex aluminium electron micrograph idealized as containing three interacting phases.

Analytical material properties
Hashin and Shtrikman [13] developed material property bounds for isotropic materials with arbitrary internal geometry, in terms of the bulk (K) and shear (G) modulus where Here, f 1 ; f 2 represent the phase volume fractions, and m 1 , m 2 Poisson's ratio for each phase. The upper bounds K U and G U can be used together with the lower bounds K L and G L as well as the well known relationships between K, G, E and m to form the averages which can be used to semi-analytically compute the effective stiffness E, and Poisson's ratio m, without any explicit information regarding the microstructural geometry. Kachanov et al. developed expressions for the effective stiffness and Poisson's ratio for a plate containing k elliptical holes, via the elastic potential of a solid with cavities of various shapes [14]. The porosity p and the hole density tensor b are given as where A is the area of the plate, a, b are the elliptical hole half-axes and n, m are the unit normals of the global coordinate system. For a circular hole the effective properties E 11 and m 12 are functions of the plate material properties E 0 , m 0 as well as p and b according to where e.g. b 11 is the first component of the hole density tensor.

Analytical stress
Kirsch derived analytical expressions for the stress in a plate containing a circular hole [15]. For general 2D loading of an infinite plate, the hoop stress r i hh;c at the hole periphery is given by r i hh;c ¼ ðr xx þ r yy Þ À 2ðr xx þ r yy Þcosð2hÞ À 4r xy sinð2hÞ; where r xx , r yy , r xy are the remotely applied stresses, and the angle h represents the polar coordinate along the circle periphery. Inglis [16] developed expressions for the maximum stress r i hh;e at an elliptical hole in an infinite plate with the ellipse half-axes a and b. The case of a rectangular notch can be modelled using an equivalent ellipse of length a with edge curvature c ¼ b 2 =a. By using this definition together with Eq. 16, the maximum stress r i hh;n can be expressed as For a finite plate of width w, containing a circular hole of radius r, the maximum hoop stress r f hh;c becomes where r 1 is the remotely applied uniaxial stress and K t is an empirical stress concentration factor [15].

Implementation
A Julia v1.5.2 [17] code denoted numerical_subdomain.jl was developed, which splits a provided micrograph into multiple individual subdomains. See Fig. 1 for an illustration of this operation. Each pixel is given one rectangular or two triangular elements, with material assignment based on the local intensity of the micrograph. Abaqus models with periodic boundary conditions (PBC) are built from each such subdomain, whose effective elastic anisotropic properties are evaluated from uniaxial-and pure shear loading using EasyPBC v1.4 [18]. Once all subdomain evaluations have been performed, the effective properties are gathered and mapped onto a reduced mesh of the same physical dimensions as the initial micrograph, by replacing each subdomain with a new element. See Fig. 2 for an illustration of this process. The input micrograph dimensions should, if possible, be chosen such that the number of pixels along each dimension is a highly composite number. Alternatively, if the number of required subdomains is known in advance, the number of pixels should have the amount of subdomains as a divisor. The overall data-flow of the numerical code is illustrated in Fig. 3. Another Julia code denoted analytical_subdomain.jl was developed, which functions in the same way, but which instead determines the effective elastic isotropic properties of each subdomain using a semianalytical method. The method is based on the average Hashin-Shtrikman bounds given by Eqs. 9 and 10, where the fractions are determined by summation over the pixels in each subdomain. It should be noted that the numerical approach operates on a finite element model of each subdomain, providing anisotropic properties, while the semi-analytical one only provides isotropic material data based on fractions.
However, the effective properties of the reduced model may still become anisotropic using both methods. Given that the application we consider later in the paper consists of three idealized phases, and the applied analytical theory only considers two phases, the following conditional logic was implemented for the three phases i, j and k: This means that only the two largest (by fraction) phases are used to determine the effective isotropic properties, neglecting the third. In other words, the phase with the largest fraction is treated as phase 1, i.e. the ''matrix'', while the phase with the second largest fraction is treated as phase 2, i.e. the ''reinforcing phase'', while the phase with the lowest fraction is neglected altogether. In practice, this becomes an increasingly good approximation for larger numbers of subdomains, since the number of pixels, and therefore also the probability of finding multiple phases in each subdomain, decreases. Using this approach, both the numerical and the semi-analytical approaches described above will converge to the same model in the limit, where each pixel would be represented by one individual subdomain.

Discretization strategy
The influence of different discretization strategies was investigated. A series of artificial micrographs were generated describing a homogeneous material with a circular hole. The aspect ratio of these micrographs was chosen such that a finite element model accurately reproduces the analytical stress concentration for a circular central hole using Eqs. 15 and 18, resulting in 400 Â 200 pixels for tensile loads, and 200 Â 200 pixels for shear loads. The radius of the circular hole was chosen such that 0-90% of the total crosssectional area was reduced in the center of each micrograph. Two imperfect discretization strategies were investigated relative to an explicit description of the circular hole. The first strategy yields a structured mesh with two first order plane stress right triangles (CPS3) per pixel of the micrograph, denoted ''1px'' because each catheti length is the same as one pixel. While rational, this results in a jagged and nonphysical description of any circular phase boundaries. The second strategy produces an unstructured mesh with two first order plane stress triangles (CPS3) per three pixels, denoted ''3px''. This approach yields some smoothening, but still produces jagged edges every other pixel. A Julia code was used to generate the 1px mesh, while the open source code OOF2 v2.1.19 [10] and the commercial pre-processor ANSA was used to generate and reshape the 3px mesh, ensuring that the aspect ratio and element density was approximately the same using both strategies. The main difference between these approaches is illustrated in a simplified way in Fig. 4.

Stress error
To quantify the error in stress caused by an imperfect discretization, the analytical stress at the circular hole was evaluated using finite plate theory as described in Sect. 2.2. For uni-axial loading, the micrographs were modelled as thin elastic rectangular plates of width w ¼ 2, height h ¼ 1 mm, thickness t ¼ 0:01 mm, stiffness E ¼ 1 GPa and Poisson's ratio m ¼ 0:3. A traction load was applied to the boundary of each plate such that the remote stress r 1 ¼ 1 MPa. A circular hole with radius r mm was introduced in the center of each plate, whereby the finite plate hoop stress r f hh was calculated using Eq. 18. For pure shear loading, a square elastic plate of width w ¼ 1 mm and height h ¼ 1 mm was used. Displacement controlled periodic boundary conditions (PBC) was applied to the plate using EasyPBC, with otherwise identical parameters to the uniaxial case. The maximum hoop stress due to a state of pure shear was calculated using Eq. 18. For both uniaxial and shear loading, an explicitly modelled circle was verified to reproduce the analytical solutions. The maximum principal stresses in all models were then compared to this numerical model. To circumvent the unavoidable pointwise stress singularities caused by using a jagged description of the circular hole, the principal stresses were averaged around the element with the highest principal stress. The averaged maximum principal stress r 1 from each model was defined as where e.g. r n 1 are the first principal stresses for the n neighbouring elements of the same material, around the element with the highest principal stress. The relative stress error for each modelr 1px 1 ,r 3px 1 was then defined based on these averages aŝ where e.g. r e 1 is the averaged maximum principal stress of the explicitly modelled hole, i.e. the Fig. 4 Simplified example of the different discretization strategies ''1px'', ''3px'' and ''explicit'' for a coarse artificial micrograph containing a circular hole. The circle diameter / plate width ratio (d/w-ratio) in this example is approximately 0.9 numerical solution is taken as reference. The same relative error definition is used throughout the paper with different reference models, i.e. for the quantity

Effective property error
The deviation in effective material properties was studied using several plates of width w ¼ 1, height h ¼ 1 mm, thickness t ¼ 0:01 mm, stiffness E ¼ 1 GPa and Poisson's ratio m ¼ 0:3. A circular hole with various choices of radius r mm was introduced in the center of each plate. The numerical solution for stiffness and Poisson's ratio was first validated using the analytical expressions in Eqs. 13 and 14. These properties are accurate when the ratio between the circle diameter d, and the plate width w, is d=w\0:6.
For larger d/w-ratios the numerical model starts to deviate from the analytical solution. While keeping the just mentioned accuracy criteria in mind, the error was computed using Eq. 22, with the effective material properties from the model containing an explicitly discretized hole as the reference solution.

Subdomain split
The influence of dividing a micrograph into individual subdomains was investigated. A 1px discretization strategy was used, since the edge case of one subdomain per pixel will reproduce the same model without any reduction. Also, while the explicit discretization strategy is applicable to the simple case of a circular hole where there exists a known analytical solution, it is not so in general for micrographs containing complex microstructures. For this reason, the reference solution will be different in this case, when compared to the investigation described in Sect. 3.2. The effective material behavior of each subdomain was evaluated using both the numerical and semi-analytical approaches described in Figs. 3 and 5, and the resulting effective elastic properties were mapped back to a geometrically homogeneous solid plate of the same dimensions. In the semi-analytical case the subdomains are isotropic, while for the numerical case they may become anisotropic. However, the resulting reduced model may still become anisotropic using both methods due to the locally heterogeneous material distribution. The effective properties of this plate were then re-evaluated using EasyPBC, and the relative error in these were computed using Eq. 22. In the numerical case, the reference solution was taken as the model with no subdomain reduction. In the semi-analytical case, the model with pixelwise reduction. Both of which are identical.

Finite plate
The artificial micrographs from Sect. 3.2 were considered in order to emulate a worst-case scenario, i.e. finite plates with central circular holes covering 10-90% of the plate width. A micrograph resolution of 200 Á 200 ¼ 40000 pixels was considered. The solid was modelled using E = 1 GPa, m = 0.3, while the hole was modelled using E = 0.000001 GPa and m = 0.3.

Aluminium micrograph
A complex aluminium microstructure was studied, where the real structure was idealized into a system of three numerical phases A, B and C. The initial model contains 30 million triangular plane stress elements (CPS3) based on a input micrograph of 15 million The micrograph describes a square physical area of approximately 1 mm 2 and was constructed by stitching 7 Á 8 ¼ 56 individual electron micrographs together. The stitching operation results in local contrast fluctuations, which were normalized using the normalize local contrast plugin of ImageJ [19]. A histogram-based segmentation was performed, defining three numerical phases using the following ranges of grayscale values; A 2 ½7501; 34999, B 2 ½35000; 65535 and C 2 ½0; 7500. The original micrograph as well as the segmented result are given in Fig. 6. The phase specific material data used were where the numerical phase A represents a-aluminium, C represents porosity and B represents all other metallurgical phases present in the microstructure. Once all subdomain evaluations were completed, the effective material properties were mapped back to a 1 mm 2 square plate for re-evaluation and analysis.

Stress error
Both the jagged ''1px'' and the smoother ''3px'' meshing strategies lead to incorrect stress predictions as seen in Fig. 7, especially for the case of small circles (d=w ¼ 0:01) where the error is on the order of 76%. However a large error is to be expected here, since the circle goes towards a square when limited by the image resolution. In fact, at the data point with the highest error in Fig. 7, the circle is represented by a perfect square which has been rotated p=4 radians around its center. For larger circles which are more accurately represented, i.e. d=w [ 0:3 with the considered pixel-density, the shear stress deviation stays below 30%. The error in uniaxial stress is in general lower than that of shear, and stays below 15% for the considered geometry. It is noteworthy that the shear stress error is 5-10% lower in the case of a jagged 1px mesh, when compared to the partially smooth 3px mesh. In general, the errors due to 1px or 3px discretization strategies become unpredictable for d=w [ 0:6. It can also be noted that stresses are underestimated in general when compared to the explicit circular hole. This could potentially be explained by the analytical expression for a rectangular notch in Eq. 17, where for a more square representation of the hole c ! 1 which leads to a reduction in the stress concentration.

Effective property error
As can be seen in Fig. 8, a jagged mesh tends to underestimate all measured material properties. The smoother mesh tends to over-estimate stiffness, but under-estimates Poisson's ratio. The errors are in general low for both cases; below 0:5% for the isotropic stiffness E, and below 0:25% for Poisson's ratio m regardless of discretization strategy. The shear modulus G ¼ E 2ð1þmÞ which depends on both properties only deviates at most 2% from the reference model. It should be highlighted once again that the error measure is taken in relation to the numerical solution, where for d=w [ 0:6 the analytical model results are no longer retained.

Finite plate: numerical method
The relative error in effective material properties determined from the numerical approach is illustrated in Figs. 9 and 10, with dashed lines representing the threshold of d=w [ 0:6. While the error is one order of magnitude higher than that coming from the choice of discretization strategy for higher d/w ratios, the error in effective stiffness is lower than 6% for d=w\0:6.
Errors in Poisson's ratio are below 25% for the same interval. With the exception of the large deviation at d/ w = 0.9, when the circle is split into four quarters at the

Finite plate: semi-analytical method
The relative error in effective material properties determined from the semi-analytical approach is illustrated in Figs. 11 and 12, with dashed lines representing the threshold of d=w [ 0:6. In contrast to the numerical method, the errors start off high and decreases for an increasing number of subdomains. The large initial error could be thought of as the shapediscrepancy, since no consideration of geometry is taken here. As an increasing number of subdomains are introduced, the shape is indirectly resolved. Once again, the plates with large hole fractions, i.e. d=w [ 0:6 are more unpredictable. As expected, it can also be noted that the errors are higher in general than for the numerical approach.

Aluminium micrograph: numerical method
Reducing the aluminium micrograph using any investigated number of subdomains does not lead to a significant error in homogenized effective material properties, as can be seen in Fig. 13. The slight increase in error for larger numbers of subdomains can be explained by the fact that the micrograph  Fig. 14 should be compared to the reference stress field in figure 15, where it can be seen that a certain level of detail can still be retained while reducing the number of elements by 99.89%. For an increasing number of subdomains, the reduced stress field tends towards the reference stress field. Like for the finite plate containing a circular hole, stresses around holes are in general underestimated in the reduced model, which becomes clear by comparing the maximum values in the legends of Figs. 14 and 15. The maximum vonMises stress is one order of magnitude higher in the reference model than in the reduced one. The computational times required for several values of subdomains are given in Fig. 19.

Aluminium micrograph: semi-analytical method
Reducing the aluminium micrograph using the semianalytical approach also does not introduce any significant error, as can be seen in Fig. 16. With the exception of the initial shape-discrepancy, the reason for the slight error fluctuation can once again be explained by the micrograph dimensions, which are not divisible by all investigated numbers of subdomains. The stress-field visualized in Fig. 17 also captures the characteristics of the reference plate for the same level of reduction. The maximum vonMises   stress is still higher in the reference model than in the reduced one, however the difference is smaller using the semi-analytical method, when compared to the numerical approach. The numerical and the semianalytical stress fields are similar, with slight local deviations. This deviation shrinks when considering even higher numbers of subdomains. By comparing the computational time of the numerical and the semianalytical methods in Figs. 18 and 19, it is clear that the semi-analytical method is several orders of magnitude faster than the numerical one.

Conclusions
While some of the investigated situations are unlikely to occur by choice in any numerical analysis, they are unavoidable when dealing with large micrographs of real microstructures. Features which are circular/ elliptical in reality will be represented in a more ''blocky'' way regardless of discretization strategy, as they fail to be completely resolved by the imaging technique used. Some microstructural features will inevitably take up large parts of certain subdomains, introducing large local errors in the analysis. However, the above results indicate that the edge-case of a circular hole in a finite plate represents a worst-case situation. The deviation in effective material properties for the applied case of an aluminium micrograph is several orders of magnitude lower than that of the finite plate. The numerical approach, while more computationally intensive than the semi-analytical one, provides effective anisotropic subdomain-level properties, as well as an explicit model of the microstructure which can be used for a more sophisticated analysis. To summarize, the major findings from this work are:

Discretization strategy
• Regardless of the chosen discretization strategy (pixelwise or smoothed), local averaged shear stress deviations of À76%, and uniaxial averaged  • Numerically determined stiffness deviates by less than 25%, and Poisson's ratio by less than 75% for the case of finite plates containing circular holes of various sizes • Semi-analytically determined stiffness deviates by less than 60%, and Poisson's ratio between À60% and 175% when determining fractions by counting pixels. The error decreases with an increasing number of subdomains

Subdomain split of aluminium micrograph
• The numerically determined errors in effective material properties were found to be DE ii 2 ½0; 0:2%; Dm 12 2 ½À0:5; 0% and DG 12 2 ½À0:2; 0% for all investigated levels of subdomain reduction • The semi-analytically determined errors in effective material properties were found to be DE 2 ½À0:2; 0:65% and Dm 2 ½À0:4; 0% for all investigated levels of subdomain reduction • The numerical approach underestimates the maximum vonMises stress the most, followed by a slightly more accurate semi-analytical approach. • Both the numerical and semi-analytical approaches are able to reproduce the characteristics of the uniaxial reference stress field, while reducing the original mesh size by 99.89% • The computational time can be reduced by several orders of magnitude by using the semi-analytical approach