(1)
where the first term corresponds to the interaction with an external magnetic field (h_i) and the second term to a spin-spin interaction, which favors ferromagnetic (antiferromagnetic) ordering in case that the coupling constant (J_{ij}) is negative (positive). The last, spin-independent term (H_0) is only an irrelevant additive constant. From the comparison of the two above expressions we identify (h_i = - 2mu _text{eff} x varepsilon _0 d) and (J_{ij} = 2 mu _text{eff} varepsilon _0^2 d ^2). First, we note that the external deformation is here analogous to the magnetic field in the Ising description. Second, the spin-spin interaction term (J_{ij}) is positive, hence favoring antiferromagnetic ordering. Also, this term is independent of the spin numbers i,j, which means that this interaction does not depend on the distance between the grains. In other words, the elastic interaction depends only on the averaged magnetization (N_+-N_-), which implies a mean field interaction.
The goal of the formulation is to minimize the elastic energy and to find the optimal spin or variant configuration ({s_i}). To this end, we use three different numerical approaches (see methods section), and the results are compared to the analytical solution above: First, a brute force approach iterates over all spin configurations to find the energetic minimum exactly, second we use simulated annealing as probabilistic ground state finder, and finally the quantum annealing approach. Fig.1a shows the resulting magnetization ((N_+-N_-)/N), i.e.the average variant orientation, as function of applied displacement (x/d Nvarepsilon _0), which corresponds to the magnetic field in the Ising model.
Results of the one-dimensional model comparing different numerical and analytical methods. (a) Mean variant orientation ((N_+-N_-) / N) as function of the displacement (x / d N varepsilon _0). Comparison between the results obtained by numerical minimization (solid lines) versus the analytical theory for an infinite and continuous system (dashed line). For large displacement, all spins are aligned and therefore the magnetization saturates. The inset shows a sketch of the one-dimensional arrangement of martensite variants (s_i=+1) (red) and (s_i=-1) (green). The bottom row is fixed to position (x=0), whereas the top grain has a mean position (x_0) in the stress free state. If an additional external stress or strain is applied, the top layer is moved to position x, and the entire microstructure is sheared to the dashed configuration. (b) Elapsed computation time as a function of the number of grains. Different methods and algorithms are compared. Dashed parts of the QA curve belong to the regime of chain breaks. For large system sizes, only the hybrid quantum annealing approach remains feasible, showing an almost constant computing time need for less than 1000 spin variables (inset).
As expected, the results agree with the analytical theory up to the aforementioned discretization effect, which becomes less pronounced for large grain numbers. For high displacements saturation sets in when all variants are de-twinned, which means that all spins are either in the state (+1) or (-1). We note that for the investigated number of spins all used algorithms lead to the same energy minimum, which confirms that also the probabilistic approaches indeed find the global minimum states.
Fig.1b shows the required computation time for the different methods and algorithms as a function of the number of grains N. All conventional algorithm implementations are based on single core computations without parallelization and are mainly shown for a qualitative comparison, as the focus of the investigations is on the quantum annealing approach. For the latter, we use quantum processing unit (QPU) implementations up to the highest possible number of spins (typically (Napprox 170) for the Pegasus architecture34 of a D-Wave machine). The brute force approach, where iterations over all spin configurations are run, has the highest computation time. Even at small spin systems of around (Napprox 40) the elapsed user time was too large for practical applications due to the simulation time scaling (sim {{mathcal {O}}}(2^N)). The pure quantum annealing method produces the fastest results and ends up with an almost constant elapsed QPU access time. Overall, the computations for (Napprox 150) are roughly three orders of magnitude faster than for the other classical approaches. Beyond around (Napprox 50) spins, so called chain breaks35 occur occasionally. They result from the need to encode strongly coupled spins as a single logical spin. Ideally, these spins should represent the same state as the individual spins, but in practise this identity can be violated. To avoid this issue and to simulate even larger systems, which are essential for higher dimensional modeling in the following sections, hybrid classical and quantum annealing approaches can be used, which combine pure QA with conventional minimization approaches36. The numerical results in Fig.1b show an increase of the computation time of the hybrid solver compared to the pure QA, but the relative acceleration compared to the classical algorithms becomes even more striking. For the hybrid solver, the elapsed computation time is essentially independent of the number of spin variables and increases only beyond (10^3) grains to several seconds. Altogether, the hybrid QA is clearly the fastest approach for large grain numbers and is therefore used in the following two-dimensional simulations.
For the determination of the linear elastic energy beyond one dimension, we consider coherent precipitates of different variants which form inside the matrix. In this way, the entire material can be considered to consist of small entities (in the following denoted as grains), which can be in one of the different martensitic states. The simplest possible (cartesian) discretization is to use small cuboidal grains with edge length a. All grains are assumed to be coherent (the elastic displacements and tractions are continuous at the interfaces between the grains), and we use homogeneous elasticity, i.e.we ignore differences in the elastic constants between the different phases or variants. This has the consequence that the elastic energy reduces to combinations of pairwise interactions between all grains37.
For demonstrational purposes we perform here two-dimensional simulations in a plane strain setup, but a transfer to three dimensions is straightforward. In particular, the annealer part does not depend on the dimensionality of the description. The qualitatively new aspect beyond 1D is the appearance of distance and orientation dependent spin-spin interactions, which decay only slowly with the distance between the grains, and therefore leads to fully populated matrices (J_{ij}). As it turns out that an accurate determination of the elastic interaction energy is essential for a precise prediction of the equilibrium microstructure, we use Fourier transformation approaches with periodic boundary conditions as outlined in the methods section. As boundary conditions, we use either vanishing average stress in the periodic volume V, (langle sigma _{ij} rangle = frac{1}{V} int sigma _{ij}(textbf{r}),dtextbf{r} = 0), or, similarly to the 1D description, a given average strain (langle varepsilon _{ij} rangle). We employ in the following for simplicity isotropic elasticity, which is e.g.described by the Lam coefficient (lambda) and the shear modulus (mu), i.e.the stress-strain relationship reads (sigma _{ij} = 2mu (varepsilon _{ij}-varepsilon _{ij}^{(0)}) + lambda delta _{ij} (varepsilon _{kk}-varepsilon _{kk}^{(0)})), where implicit summation over repeated indices is used. The position dependent eigenstrain (varepsilon _{ij}^{(0)}(textbf{r})) is known for a given microstructure with fixed phase dependent stress free strains (relative to the austenitic mother phase), (varepsilon _{ij}^{(0)}(textbf{r}) = theta (textbf{r}) varepsilon _{ij}^0), where the indicator function (theta) is zero in the austenite and either (+1) or (-1) in the two considered martensite variants. For a given microstructure, the elastic energy can then be computed in reciprocal space, as shown in the methods section. For the formulation as Ising model we discretize our microstructure using small non-overlapping cuboidal grains as discussed above and assign a spin (s_i) to each of them like before, such that the indicator field becomes a superposition (theta (textbf{r}) = sum _i s_i theta _i(textbf{r})), where (theta _i) equals one inside the corresponding square and is zero outside. Therefore, the elastic energy decomposes into pairwise interactions (for (ine j)) and self-energy terms (for (i=j))
$$begin{aligned} E_{i,j} = s_i s_j frac{1}{2V} int dtextbf{r} int dmathbf {r'} B(textbf{r}-mathbf {r'}) theta _i(textbf{r}) theta _j(mathbf {r'}), end{aligned}$$
(2)
where the integral kernel (B(textbf{r})) is defined through the inverse of the elastic Greens function. Hence, it is sufficient to perform the Fourier transform energy calculations for all pairs of the same martensite variant (s_i=s_j=1) on the discrete lattice sites in the volume V; for periodic boundary conditions and identical grain shapes, it is sufficient to calculate the elastic interaction energy between a reference grain and all the other grains, due to translation invariance. In case of fixed average strain boundary conditions, an additional homogeneous term appears (see methods section), contributing both to the spin-spin interaction (J_{ij}) as well as to the magnetic field term (h_i), which is absent for zero average stress boundary conditions. The resulting fully populated matrix of coupling constants with both positive and negative entries has similarities to spin glass systems with random couplings, which have been investigated in the literature with conventional approaches, see e.g.38.
For the simplest case that the eigenstrain is purely dilatational and isotropic the Bitter-Crum theorem applies and the total energy depends only on the volume fraction of the martensite variant, where no interaction between the grains is present and only a self energy term remains39.
For a nontrivial elastic interaction and the link to the previous 1D description, we consider a shear transformation strain with (varepsilon _{xy}^0=varepsilon _{yx}^0=varepsilon _0), where all other components vanish. In this case, we obtain a distance and orientation dependent interaction as depicted in Fig.2a, which is computed here for the case of vanishing average stress, (langle sigma _{ik}rangle =0). Here and in the following parts the Poisson ratio is chosen as (nu =1/4) (i.e.(lambda =mu)).
Interaction energies of two grains of equal variant type ((mathbf {s_i=s_j})). Interaction energies in the case of (a) shear eigenstrain and vanishing average stress and (b) tetragonal eigenstrain. The interaction energy per length is given in units of (lambda a^3 varepsilon _0^2), and the computations were done using a system size of (L_x/a=L_y/a=50), where a is the edge length of the grains. At distance (r/a=0) the grains touch each other. The symbols on the continuous curves indicate the information for the interaction at discrete lattice sites, which is actually used in the annealer simulations.
The interaction energy is obtained by subtracting the elastic self energies (E_text{self}) for each of the two (isolated) martensite grains inside the austenitic matrix from the total elastic energy (E_text{el}) of the two-grain arrangement, i.e.(E_text{int}=E_text{el}-2E_text{self}), to normalize the interaction energy such that it decays to zero for large grain separations. For short distances, a transition between attraction and repulsion is found for the (langle 100rangle) direction, whereas a purely repulsive interaction results for the diagonal (langle 110rangle) directions. Due to the periodic boundary conditions, the result depends on the system size (V=L_xtimes L_y), as the grains also interact with their periodic images, hence (rll L_x, L_y) is required to observe the decay of the interaction.
We note that in two dimensions the interaction energy decays asymptotically as (r^{-2}), whereas in three dimensions it scales as (r^{-3}) in large systems, which follows from the elastic Greens function40. For the quantum annealer implementation, the interaction energies are needed only for the discrete lattice points (symbols on the curves). Although the decay of the elastic interaction may suggest to cut it off beyond a certain distance in real space, it turns out that such an approach is inappropriate, as it leads in the end to invalid equilibrium microstructures, and it is therefore essential to keep all interaction terms (J_{ij}) with high precision to avoid spurious effects. We note that the formulation on the quantum annealer does not depend on the dimensionality, therefore the scaling plot in Fig.1b applies here as well.
Based on the calculation of the elastic interactions, we obtain from the Ising model implementation on the quantum annealer with hybrid solver stripe patterns in (langle 100rangle) direction as equilibrium structures. These patterns are irregular in the sense that the widths of the stripes are not uniform. This is in analogy to the 1D model, which was discussed above, where we found that the arrangement of the two variants is not determined. This coincidence, which is physically expected, is nontrivial from the model formulation, as (i) in the 1D model we had a distance independent interaction in the discretized model, where here the interaction is significantly more complex, but adds up to the same effective descriptions for the periodic arrangement; (ii) a rotation of the pattern by 90 degree is possible and sometimes obtained from the optimal configuration due to the discrete rotational symmetry; (iii) the fixing of the average stress compared to the given average strain in the 1D formulation can lead to unequal distributions of the different variants. In particular, for the presently considered absence of an external strain (implying a vanishing magnetic field in the Ising terminology), there is no constraint of the sort (langle s_irangle = 0) for the average spin alignment. All stripe configurations are energetically equivalent, which includes the possibility of a single variant configuration. These results therefore confirm simultaneously the accuracy of the elastic interaction calculation with the pairwise decomposition as well as the ability of the quantum annealer to identify the true ground state configurations.
As next example, we use a tetragonal eigenstrain with the only nonvanishing components (varepsilon _{xx}^0=-varepsilon _{yy}^0=varepsilon _{zz}^0=varepsilon _0). First, we consider again the situation with vanishing average stress, (langle sigma _{ij}rangle = 0). The corresponding interaction energy is shown in Fig.2b for (nu =1/4). In this case, the equilibrium microstructure is trivial and consists of a single variant, as in this case the elastic energy is zero for the periodic system. Therefore, the situation differs from the previous shear transformation, where also lamellar arrangements with both variants lead to stress free situations. The reason is that any interface between two variants leads to a mismatch between adjacent variants for the tetragonal transformation, and therefore such a situation is energetically unfavorable here. A change of boundary conditions to vanishing average strain, (langle varepsilon _{ij}rangle =0), alters the situation, since then arrangements with equal amounts of both variants are preferred, as this lowers the volumic part of the elastic energy. In this case, we find regular inclined stripes as equilibrium pattern, as shown in Fig.3a.
Resulting stripe patterns for tetragonal eigenstrain. (a) Equilibrium structure with three stripe pairs (counted along the horizontal axis) in a system consisting of (50times 50) cuboidal grains. A vanishing mean strain, (langle varepsilon _{ij}rangle =0), is imposed. The width of the stripes is uniform, consisting of grains with configuration (s_i=+1) (red) and (s_i=-1) (green). (b) Elastic energy of stripe patterns with different inclination angles (phi .) The solid curves correspond to smooth stripes (the grain size (a/L_x, a/L_yll 1) is negligible) and show a pronounced stationary point for inclinations for which the pattern repeats periodically without kinks at the boundaries. The squares correspond to situations with the same number of stripes, where the system is discretized by (50times 50) quadratic grains, leading to pronounced aliasing effects, and the resulting elastic energy is higher than for the smooth stripes. This shifts the energetic minimum for 6 stripe pairs at (phi approx 40^circ) to a lower angle (phi approx 33^circ) with 3 stripe pairs. The infinite system size limit for smooth stripes is depicted as black dotted curve.
Again, the solution is not unique; in particular, due to translation invariance, the annealer returns also configurations where the stripes are shifted. Also, a switch of the sign of the inclination angle (phi) (see definition in the figure) leads to energetically equivalent solutions. However, we do not find ground state configurations which lead to different (absolute) inclination angles or strip widths or even irregular variations of the latter, contrary to the shear transformation case before.
The reason for the observed ground state morphologies is a combination of continuum elasticity effects, the granular structure of the material and constraints induced by periodic boundary conditions. Figure3b shows the computed elastic energy for different numbers of regular arrangements of stripes in the periodic system as function of the inclination angle (phi). Here we see a pronounced influence of the grain size, as the elastic energy of configurations with regular stripe pairs with a discretization by (50times 50) grains (squares in the figure) is higher than for corresponding cases with very fine grains, where discretization effects do not play a role anymore (smooth curves). The oscillating nature is due to the periodic boundary conditions, as improper choices of the inclination angle lead to discontinuities of the stripe patterns at the boundaries, which is energetically unfavorable. Therefore, continuous patterns correspond to the stationary points of the curves. For specific angles, the curves for three and six stripe pairs meet at local minima, which is a consequence of the scale invariance of linear elasticity. From the smooth, continuum limit curves one would conclude that an angle of about (phi approx 40^circ) should lead to the energetically lowest configuration (absolute minimum of the smooth red curve). Moreover, in the limit of infinite systems, where periodic boundary conditions do not play a role anymore, an analytical treatment is possible, leading to the energy expression (E_text{el}^infty = V B(n)/2) for equal volume fraction of the two variants with
$$B(n) = frac{{4mu }}{{lambda + 2mu }}varepsilon _{0}^{2} {text{ }}left[ {(3lambda + 2mu ) - 2(3lambda + 2mu )n^{2} + 4(lambda + mu )n^{4} } right]$$
with (n=cos phi). Energy minimization gives the optimal angle (phi =cos ^{-1}sqrt{5/8}approx 37.8^circ), see Fig.3b (minimum of the black dotted curve).
However, these predictions disagree with the finding from the quantum annealer, which favors a configuration with three stripe pairs at a lower angle of (phi approx 33^circ). This observation can be understood by consideration of the granular structure of the patterns investigated here, as the microstructure in the annealer simulations consists of (50times 50) square grains. First, the explicit appearance of the length scale a breaks the scale invariance of the periodic pattern, and therefore the minima of the energy curves belonging to the discrete microstructures (squares in Fig.3b) do not coincide anymore at the local minima. Additionally, with increasing inclination antialiasing effects of the patterns become more relevant, and therefore the energy curves show an increasing disagreement with the continuum limit curves. As a result, the energetic minimum in the discrete microstructure indeed shifts toward a configuration with three stripe pairs at (phi approx 33^circ) (absolute minimum of the blue squares in Fig.3b), which is in agreement with the prediction of the quantum annealer. Consequently, details of the granular structure can change the energetics compared to a full continuum approximation, especially since many local minima of the elastic energy are located close to each other.
The approach presented above is not limited to mutually interacting cuboidal grains, but can also be applied to realistic microstructures. To illustrate the procedures, we have generated a microstructure consisting of (N=400) grains using a Voronoi tesselation41. Each grain is allowed to take one out of two martensite variants with the tetragonal eigenstrain tensor, and we pre-compute all mutual elastic interactions between them. We note that contrary to the case with the cuboidal grains in a periodic array here we cannot exploit translational invarince due to the different shapes of the grains, and hence these elastic interaction energy calculations scale here as ({{mathcal {O}}}(N^2)) instead of ({{mathcal {O}}}(N)) before, although we still use periodic boundary conditions. Additionally, we consider now arbitrary given external strains (langle varepsilon _{ij}rangle), which leads to the appearance of a magnetic term like in the one dimensional description. With that, we can predict the equilibrium variant distribution within the microstructure using the hybrid quantum annealer, and this step is typically executed within a few seconds of runtime.
Examples for the equilibrium microstructures are shown in Fig.4 as function of the externally applied strain (langle varepsilon _{xx}rangle), whereas the other average strain components vanish.
Resulting equilibrium variant distribution with uniform grain orientation. The microstructures consist of 400 grains and tensile strain is applied in horizontal (x) direction. Red (green) grains correspond to variant (s_i=+1) ((s_i=-1)). The tensile strain is (a) (langle varepsilon _{xx}rangle /varepsilon _0 = 0), (b) (langle varepsilon _{xx}rangle /varepsilon _0 = 0.1), (c) (langle varepsilon _{xx}rangle /varepsilon _0 = 0.5), (d) (langle varepsilon _{xx}rangle /varepsilon _0 = 0.9), (e) (langle varepsilon _{xx}rangle /varepsilon _0 = 1.1) and (f) (langle varepsilon _{xx}rangle /varepsilon _0 = 1.3).
The observed microstructures are indeed similar to what we have found before using the square discretization, although here the band widths and orientation deviate from the previous case due to microstructural details and the smaller number of grains, and these effects can be explained using an analysis similar to the one done for Fig.3b. We note that in these microstructures all grains have the same orientation, and therefore the application of a tensile strain strongly favors the selection of the grain variant (s_i=+1) (for a compressive situation we observe the opposite behavior), and we find a full alignment of all variants in the last snapshot.
Additionally, we have performed the same analysis for grains with uniformly distributed random orientation, which implies a rotation of the local transformation strain tensor, see Fig.5 for the grain orientations and for the variant selection.
Resulting equilibrium variant distribution with random grain orientation. (a) Grain orientation map corresponding to the microstructures. In the color bar the grain rotation angle is given in radian (modulo (pi) due to symmetry). The rotation axis is along the [001] direction. The microstructures consist of 400 grains and tensile strain is applied in horizontal (x) direction. The grains have a random orientation, which is the same for all cases, based on a uniform distribution. The tensile strain in horizontal direction is (b) (langle varepsilon _{xx}rangle /varepsilon _0 = 0) and (c) (langle varepsilon _{xx}rangle /varepsilon _0 = 2.1). Red (green) grains correspond to variant (s_i=+1) ((s_i=-1)).
Here, also the equilibrated spatial distribution of the variants appears to be irregular. Application of a tensile strain again favors the alignment of the variant, but this time even for high strains not all grains select the same variant, which is due to the local rotation. In fact, a grain, which is rotated by (90^circ) with respect to the straining direction has a preference to be in variant state (s_i=-1), as then the direction of expansion is aligned with the external tensile strain. This can be clearly seen e.g.in Fig.5(c) for the highest tensile strain in x direction, where the remaining patches with spin (s_i=-1) correspond to the grains with orientation close to (pi /2) (or (3pi /2)). We emphasize that for a given microstructure (shapes of all grains) the mutual elastic grain-grain interactions have to be computed only once. As mentioned before, this step has to be done with high precision, and consequently this is the step which demands the highest amount of computing time. After that, all changes of the external boundary conditions affect only the (k=0) mode contributing to the interactions (J_{ij}) and (h_i), and these terms can be calculated analytically (see methods section). As each hybrid quantum annealing calculation typically requires only a few seconds, the entire microstructural change during mechanical loading can be calculated extremely fast.
See original here:
Quantum annealing for microstructure equilibration with long-range ... - Nature.com