Results in Physics

We present a computational framework for two-scale asymptotic homogenization to determine the intrinsic magnetic permeability of composites. To this end, considering linear magnetostatics, both vector and scalar potential formulations are used. Our homogenization algorithm for solving the cell problem is based on the displacement method presented in Lukkassen et al. 1995, Composites Engineering, 5(5), 519-531. We propose the use of the meridional eccentricity of the permeability tensor ellipsoid as an anisotropy index quantifying the degree of directionality in the linear magnetic response. As application problems, 2D regular and random microstructures with overlapping and nonoverlapping monodisperse disks, all of which are periodic, are considered. We show that, for the vanishing corrector function, the derived effective magnetic permeability tensor gives the (lower) Reuss and (upper) Voigt bounds with the vector and scalar potential formulations, respectively. Our results with periodic boundary conditions show an excellent agreement with analytical solutions for regular composites, whereas, for random heterogeneous materials, their convergence with volume element size is fast. Predictions for material systems with monodisperse overlapping disks for a given inclusion volume fraction provide the highest magnetic permeability with the most increased inclusion interaction. In contrast, the disk arrangements in regular square lattices result in the lowest magnetic permeability and inadequate inclusion interaction. Such differences are beyond the reach of the isotropic effective medium theories, which use only the phase volume fraction and shape as mere statistical microstructural descriptors.


Introduction
Modern manufacturing processes are becoming increasingly flexible and reactive to changing process conditions [1].In the case of metal manufacturing processes such as machining, forming, and additive manufacturing, the material microstructure is a crucial factor for process performance, as it determines the local mechanical response of the material.At the same time, the microstructure evolves during production under thermo-mechanical loading.To enable the transition to Industry 4.0, there is a strong need for the development of in-line microstructural observers [2].
A promising technique in this regard is eddy current measurement.This method is used for many purposes, from inferring geometrical properties of the composites, e.g., coating thickness [3], film thickness [4] to size and location of defects [5].Several studies have shown that the electromagnetic sensor response can be related to mechanical properties such as formability, hardness, or anisotropy through extensive experimental calibration [6,7].However, a better understanding * Chair of Nonlinear Solid Mechanics, Faculty of Engineering Technology, University of Twente, 7522 NB Enschede, The Netherlands.E-mail address: c.soyarslan@utwente.nl (C.Soyarslan). of the relations between the electromagnetic and mechanical response of the materials is required for a more effective application of these techniques.Several recent studies have been devoted to the modeling of electromagnetic measurements systems to infer certain statistical moments of the material microstructure using in-line measurement data [8][9][10][11][12][13].
In the low-frequency regimes, what affects the inductive sensor output is mainly the effective magnetic permeability of the material.In contrast, effective electrical conductivity becomes important at frequencies where eddy current losses become significant [13].For most alloy steels, the constituents are distinguished by their magnetic permeabilities, whereas the property contrast in their electrical conductivities (or resistivities) are relatively small [11].An example is the two-phase ferritic-austenitic steel in which (ferromagnetic) ferrite and (paramagnetic) austenite has a magnetic permeability contrast of around  Fe ∕ Aus ≃ 1000.In contrast, their electrical resistivities are comparable within a measurement uncertainty of 7% [14].In the https://doi.org/10.1016/j.rinp.2022.106188Received 28 November 2022; Accepted 15 December 2022 determination of the effective magnetic permeability, a micro-to-macro transformation of the constituent phase properties is required.
Analytical methods developed to this end constitute various effective-medium theories such as the Maxwell, self-consistent, and differential effective-medium approximations.The (arithmetic) Voigt and the (harmonic) Reuss averages form analytical bounds in this respect.These approximations rely on the simplest microstructural descriptors, such as phase volume fraction and shape.While most effective-medium approximations mixture rules yield high accuracy for low-contrast material properties, high-contrast physical properties, and cases with intricate phase interactions, e.g., due to phase percolation, these mixture rules lose their accuracy.
A full-field numerical computation that accounts for fine microstructural details is required to remedy this.Such computational methods using field averages on material volumes are modeled either as embeddings into free space with uniform external magnetic fields or as isolated closed-magnetic circuits [15] subjected to various boundary conditions.In embedded models, due to demagnetizing fields, the resultant predictions underestimate the intrinsic magnetic permeability [12].The use of isolated microstructures with appropriate boundary conditions gives accurate results in determining the intrinsic magnetic properties of the material.In this context, besides the boundary conditions proposed in [12], those satisfying Hill-Mandel conditions are discussed in, e.g., [16][17][18][19].
In this work, we present a computational method in first-order twoscale asymptotic homogenization to determine the tensorial magnetic permeability of composites.Limiting ourselves to linear magnetostatics, our derivation of the cell problem and its solution follows the displacement method of Lukassen et al. [20], which was initially introduced for an elastomechanical problem.In our formulation of the uncoupled nondimensional magnetostatics equations, we use scalar and vector potentials where the corrector function is identical-type to the applied magnetic potential.We use both uniform Dirichlet and periodic boundary conditions.Inspired by Lamé's stress ellipsoid and considering that both the stress and the permeability tensors are of second-order, we propose the use of meridional eccentricity in quantifying the anisotropy in the effective magnetic permeability of the composites.
The developed computational framework is applied to determine the effective properties of regular and random material systems in two dimensions.The motivation behind this choice is two-fold: The first motivation is computational tractability.Secondly, with appropriate changes of variables and parameters, this permits a unified treatment of both scalar and vector potentials through a generic partial differential equation, which becomes handy in applying Dirichlet-based periodic and uniform boundary conditions.Our results are compared with the available analytical solutions and predictions of the selected effectivemedium theories and rigorous bounds.In light of our numerical computations, we proposed a simple three-parameter approximation for the effective permeability of overlapping random disk systems.

A word on notation
Let us denote by (, ) a scalar field and by (, ) and (, ) two vector fields, distributed over the domain represented by material points  at time .Considering three-dimensional space and using Einstein's summation convention, we write  =     and  =     where   for  = 1, 2, 3 represent Cartesian basis vectors and   and   the associated vector components.The gradient of a scalar field, divergence, and curl of vector fields are respectively defined as follows Here, ⋅ and × denote single-contraction (scalar) and cross products with  ⋅  =     and  ×  =         .  is the Levi-Civita symbol, such that   is 1∕ − 1 for even/odd permutations of (, , ) and 0 for repeated indices.  is the Kronecker delta with   = 1 if  =  and   = 0 otherwise.

Nondimensional magnetostatics
Let ,  and M  denote the dimensional magnetic field, magnetic induction field and the position of the particle, respectively, and, ,  and M  = M ∕ M  ref their nondimensionalized counterparts, with M  ref representing macroscopic a reference scale length. 1 Ampère-Maxwell law for steady state conditions links the nondimensional free current density  to  with curlM   =  . ( Gauss's law for magnetism which asserts the absence of magnetic monopoles implying a solenoidal magnetic induction field  by Assumption of linear magnetostatics leads to the following constitutive relation between  and   =  ⋅  and  =  ⋅  , where,  =  −1 is the nondimensional symmetric second-order magnetic permeability tensor, where   =   and   =   for ,  = 1, 2, 3.
In the formulation of magnetostatics, there exist two possibilities.The first one is the more general vector potential formulation which links the magnetic induction field to a magnetic vector potential  with Alternatively, for vanishing free currents, i.e.,  → , using tensor calculus in view of Eq. ( 1) allows one postulate a magnetic scalar potential  from which the magnetic field  is derived With this, a scalar potential formulation is developed.In the following derivations, for the sake of completeness, we pursue both vector and scalar potential formulations.

Asymptotic homogenization
The asymptotic homogenization method used here is adapted from the displacement method presented in [20] in a mechanics context.Similar developments can be found in, e.g., [21][22][23][24].
Let  denote a periodic unit cell of size  ref ×  ref ×  ref whose tiling fills the composite material domain.Although the conclusion is not affected if we have more constituent phases, we consider two linearly magnetic phases which occupy unit cell subdomains  1 ,  2 where  1 ∪  2 = .To characterize the scale of the unit cell, another nondimensional coordinate  can be generated with Here, 0 <  ≪ 1 controls the fineness (scale separation) of the unit cell, see Fig. 1.If  is small, the result is a fine microstructure with rapidly oscillating property functions.Let   ( M ) =  ( M , ) denote a generic tensor-valued function on the composite domain   =  × , see, e.g., Ref. [25], which is −periodic in  over the unit cell  and which features slow (with M ) and fast (with ) variations.It follows that where sub-indices M  and  denote derivatives with respect to the macroscopic and microscopic coordinates M  and , respectively.
In the following, we derive expressions with vector and scalar potential formulations of asymptotic two-scale homogenization in magnetostatics in parallel.To this end, we start by assuming that the magnetic vector and scalar potentials are multiscale functions of the form   ( M ) = ( M , ) and   ( M ) = ( M , ), respectively.Then, we apply two-scale asymptotic expansion to the nondimensional potentials as the primary solution variables on   ,   ( M ) =  (0) ( M , ) +  (1) ( M , ) +  2  (2) ( M ) =  (0) ( M , ) +  (1) ( M , ) +  2  (2) ( M , ) + ( 3 ) , (10) where both  () and  () are −periodic functions in .The leading orders  (0) and  (0) respectively correspond to the homogenized magnetic vector and scalar potential fields, and as shown below, they only depend on the macroscopic coordinates M , i.e., they are constant over the unit cell domain.

Computation of components through Dirichlet boundary conditions
The cell problems given in Eqs. ( 33) and ( 34), respectively given for vector and scalar potential formulations, can be solved for   () and   () by applying Neumann boundary conditions corresponding to the right-hand-side source term at the phase boundaries, see, e.g., [20].Algorithmically, this construction, adapted from the displacement method of [20], uses only Dirichlet boundary condition associations at the cell boundary and does not require the detection of phase boundary surfaces for the imposition of boundary conditions.In what follows, we give details of our implementation.For simplicity of notation, from now on, we denote the nondimensional cell-scale fields  (0) and  (0)  with  and , respectively.

C. Soyarslan et al.
For each of the magnetic vector and scalar potential formulations, we respectively seek three solutions in the corresponding total magnetic potential denoted by   () and   () indexed with .For  = 1, 2, 3,   () and   () have the following additively decomposed forms, respectively Here,   () and   () are the unknown periodic fluctuation fields whereas, as indicated above,   () and   () are known. 2Under the above definitions and using along with Eqs. ( 35) and ( 36) for vector and scalar potential formulations, we respectively obtain 3 If we apply Eqs. ( 37) and (38), respectively, the conjugate magnetic fields   and   are obtained Averaging over the unit cell, corresponding homogenized magnetic fields M   and M   are obtained Comparison of above and Eqs. ( 43) and ( 44) yields the following expressions in matrix form for the inverse permeability and permeability tensors computed using vector and scalar potential formulations, respectively This is to say that, for the vector potential formulation, the component of the macroscopic magnetic permeability tensor  ⋆ with the indices  corresponds to the homogenized magnetic induction tensor component M   due to total magnetic vector potential influence function   computed for the imposed macroscopic magnetic field at unit cell vertices with   .The same reasoning applies to the case of the scalar potential formulation; however, this time, the components of the macroscopic magnetic permeability tensor  ⋆ are inferred.Although, for the most general case, the computation of the six linearly independent macroscopic constitutive constants making up the permeability tensor requires three load cases given in Eqs. ( 47) and ( 48) for  = 1, 2, 3, once the symmetry elements of the point group of the underlying 2 The known part is imposed at the control nodes, e.g., the periodic unit cell corners.In this case,   () and   () vanish at the periodic unit cell corners.Alternatively, the control nodes can be selected over a reference unit cube.In this case, the uniqueness of the fluctuation field can be provided by fixing it at an arbitrary node in the unit cell domain. 3Application of averaging given in Eqs.(41) to Eqs. ( 49) and ( 50) with the use of Eqs.(32) give the following unit vector representations of the macroscopic magnetic fields M   =   and M   =   for vector and scalar potential formulations, respectively , and composite geometry, which should be included in the physical property symmetry elements with Neumann's principle [27,28], allow, less number of tests can suffice.For instance, considering two-dimensional tensor formalism, the magnetic response of microstructures possessing cubic symmetry is isotropic; in consequence, their magnetic permeability tensor is spherically symmetric.Therefore, only one test can determine the corresponding nonzero main-diagonal element.The algorithmic treatment of the homogenization procedure for vector and scalar potential formulations are summarized in Algorithms 1 and 2, respectively.
Algorithm 1: FE-based computational periodic homogenization for vector potential formulation.Macroscopic loading is applied with PBC to result in M . .

Magnetic anisotropy index
In this part, we develop an index for quantification of the degree of anisotropy in the magnetic permeability of materials.To this end, we use a permeability ellipsoid (or ellipse in 2D) for the graphical representation of the second-order symmetric permeability tensor.Among other possible ellipsoids making the indicatrice surfaces geometrically representing second-order tensors, ours is the ellipsoid of the values of symmetric second-order tensor [29, p. 22].For a second-order stress tensor, this corresponds to Lamé's stress ellipsoid [30, p. 215].The corresponding meridional eccentricity of the ellipsoid, which is the eccentricity of the ellipse at the section of the longest and the shortest radii, respectively denoted by  max =  max and  min =  min , which gives max is then used as the anisotropy index with 0 ≤  ≤ 1.Here,  max and  min are the largest and the smallest eigenvalues of the permeability tensor.For  = 0, the permeability ellipsoid transforms onto a sphere which indicates isotropy in magnetic permeability.As  increases, the directionality in the material's permeability also does so.

Motivation: 2D magnetostatics
We are concerned only about two-dimensional problems for which the fields  and  have merely in-plane components with Considering two-phase composites, indexed with  = 1, 2, we let   =    and   =    where   = 1∕  with the assumption of isotropy at the microscale.Here  is the second-order identity tensor.This allows representing the governing differential equations to be solved over the unit cell, see Algorithms 1 and 2, for both the vector and the scalar potential formulations, that is curl   = 0 and div   = 0 respectively, in the following generic form in Cartesian components 4     1 Here, the vector potential formulation is recovered through replacements  ←  and  ←  3 , whereas for the scalar potential formulation, we use the substitutions  ←  and  ← .

Considered microstructures and associated model generation
In the current section, we consider regular and random arrangements of nonoverlapping and overlapping monodisperse disks in a two-dimensional matrix.All generations are periodic.In three dimensions, these correspond to unidirectional circular cylindrical inclusions in a matrix material, which extend indefinitely along the main axis of cylindrical inclusions.We take  i = 250 0 and  m =  0 for the inclusion and the matrix, respectively.We also consider phaseinterchange, which corresponds to interchanged phases contemplating the same microstructure [31].

Regular nonoverlapping square disk arrangements
In periodic material systems, the choice for the primitive unit cell is not unique.Thus, we consider three possible primitive unit cells and their  ×  spatial tilings for  = {2, 4, 8, 16, 32}, henceforth referred to as the volume elements (VEs), see Fig. 2. As a first task, we focus on the effect of boundary conditions and the choice of the potential on which the formulation is based, in conjunction with the type and 4 Besides magnetostatics in scalar potential, the elliptic differential equation div   = 0, governs dielectrics, electrical conduction, and thermal conduction.Considering two-dimensional reduction, vector potential formulation can also be treated in the same formalism.This allows one to use any available physics engine in the software to solve for magnetostatics by carefully mapping parameters.In the current study, we exploit this condition and use Abaqus coupled thermo-electrical analysis framework by omitting thermal fields.As opposed to Abaqus edge-based elements, which are generally used for modeling magnetostatics with vector potential, current use corresponds to a node-based approach.An important point to be considered is during post-processing.Although in scalar potential formulation, derivation of the magnetic field  = −   which uses the exact form in electrical conductivity as  = −   with  and  denote the electrical field and scalar electric potential, respectively, which results in  = ∕ 1  1 + ∕ 2  2 and  = ∕ 1  1 + ∕ 2  2 .In magnetic vector potential, however,  = curl  , which for the current condition where size of the used volume element on the effective properties.The type of the volume element controls its boundary property fluctuations, whereas its size controls the relative size of the emerging boundary layer.The models are subjected to periodic and uniform Dirichlet boundary conditions with both scalar and vector potential formulations of magnetostatics.Here we fix the inclusion volume fraction at  i = 0.30.A geometry-based approach, which requires mesh seeds at the phase interface, is adopted to model phase boundaries accurately.The primitive unit cell domain is discretized using 12 640 node-based linear quadrilateral elements forming a structured mesh.To eliminate any mesh-related differences, the unit cells are meshed identically.
For the studied periodic square disk arrangements, a series expansion solution in terms of the disk volume fraction is given for the effective conductivity in [32].As a second task, we compare our numerical results with this analytical solution and consider the primitive unit cell with centered disk arrangement given in the first row of Fig. 2. We created models with inclusion volume fractions of  i = {0.05,0.10, … , 0.55}.Each cell domain is finely discretized using over 10 000 node-based linear quadrilateral elements forming a structured mesh.

Random disk arrangements
In this part, we consider nonoverlapping and overlapping periodic disk arrangements which possess statistical ergodicity and uniformity.The microstructures are generated using a Monte Carlo method depending on a sequential adsorption process.The periodicity is enforced as follows: we consider a square region with size  × .We create a 3 × 3 continuous tiling of this region.Letting  denote uniform probability distribution, we generate a disk with a fixed radius , hence monodispersity, and random central coordinates (  ,   ) where   ∼  (0, ) and   ∼  (0, ).Taking the origin of each square region as a reference for each disk coordinate, we embed nine copies of this disk in our 3 × 3 square tiling.This process is repeated until the desired inclusion volume fraction is reached.Finally, we clip out the central square region, which possesses periodicity.
Generating a mesh with periodic nodal locations can be difficult for random microstructures.Although the analytical computation of inclusion volume fraction is trivial for nonoverlapping disks, for the overlapping ones, it is not.Thus, to remedy both challenges, a voxelbased meshing strategy is followed in this part.However, this introduces a relatively higher error in the representation of volumes and surfaces compared to geometry-based meshing approaches.
According to Hill [33], a representative volume element (RVE) should yield the same effective properties for the applied displacement or traction boundary conditions as far as they are macroscopically uniform.For random composites whose constituents have high-contrast physical properties, this is impractical as it results in substantial RVE sizes [34,35].In this part, we start with searching for the RVE size but adapting the definition, which defines it as the smallest material volume size that accurately represents the composite's average physical property [36].For regular periodic composites, the primitive unit cell makes up the RVE.For random disk-matrix systems, this is not the case.Such systems require a large set of unique microstructures of the same size and shape to be considered.What makes each microstructure unique is its underlying stochastic generation.These are also referred to as statistical volume elements (SVEs).The computation of effective properties requires an averaging over an ensemble of microstructures; an operation referred to as the ensemble averaging, see, e.g., [37, p. 85] and [38].The ergodic hypothesis asserts that averaging over all ensemble elements tends to volume averaging on any ensemble elements with increasing SVE size.In contrast, an equivalence is possible when the SVE size tends to infinity.For random composites, the minimal RVE size is selected by monitoring the precision of the mean value of conducted realizations [39].This size is generally not unique, and it depends on the selected physical property, and the contrast of the constituent phase properties [39].Based on these concepts, we adopt  a pragmatic methodology and consider sparsely populated ensembles with  = 15 random microstructures, VEs (or SVEs); for each selected VE size, we compute mean μ, standard deviation σ and standard error of the mean SEM = σ∕ √  in the resultant effective properties while keeping the inclusion volume fraction constant with  i = 0.30.Although not the case for overlapping disks, for nonoverlapping monodisperse disks, the volume fraction increase with each added disk, and the associated relative error is known apriori as a function of the number of disks and the VE-size.This restricts the disk volume fractions that can be obtained, especially for small cell size-to-disk radius ratios.In our generation, we use a disk diameter of  = 34 and a minimum VE size of 128, which results in a 128 × 128 pixel-based discretization.This kept the associated absolute error in volume generation smaller than 0.01.In addition, with systematic doubling, we considered square VEs with 128 × 128, 256 × 256, 512 × 512, 1024 × 1024, and 2048 × 2048 discretizations.Exemplary periodic generations for nonoverlapping and overlapping disk arrangements are demonstrated in Fig. 3.
As a second task, we keep the VE size constant and vary the inclusion volume fraction to investigate the potential difference between overlapping and nonoverlapping disk system responses and between them and regular square disk arrangements.For nonoverlapping monodisperse disks, the random sequential adsorption process has the saturation limit at  s ≃ 0.55 in two-dimensions [40] whereas, for threedimensions with spheres, it is  s ≃ 0.38 [41].This is the limit of jamming beyond which no other random disks can be added to the system without overlapping.In Fig. 4, we demonstrate five randomly generated microstructures with nonoverlapping disk arrangements with volume fractions of  i = {0.15,0.25, … , 0.55}.For each volume fraction 15 generations were realized.
As opposed to the case for nonoverlapping disks case, for overlapping disks, the attainable inclusion volume fraction is not limited.In Fig. 5, we demonstrate five randomly generated microstructures with nonoverlapping disk arrangements with volume fractions  i = {0.15,0.25, … , 0.95}.Like before, 15 generations were realized for each volume fraction.Verification of the Monte Carlo generations of the overlapping and nonoverlapping microstructures with twopoint probability functions is given in Appendix ''Two-point Probability Function''.

Regular nonoverlapping square disk arrangements
On account of the involved phase fluctuations at the cell boundaries, the  ×  tilings of the unit cell given in row 3 of Fig. 2 are selected for computation of the effective permeability plots demonstrated in Fig.
for  = {2, 4, 8, 16, 32}.The vector and scalar potential formulations, once used in conjunction with periodic boundary conditions, give results in agreement with each other with the effective permeability predictions reading 1.850 0 and 135.1 0 , considering the cases without and with phase-interchange, respectively.Moreover, these results show size and unit cell type invariance.That is, the predicted effective permeability change neither with the tiling parameter  nor the change  of the unit cell type.As can be calculated, the phase-interchange relation [42] for a two-phase composite formulated by the equation  ⋆ ( 1 ,  2 )  ⋆ ( 2 ,  1 ) =  1  2 is satisfied within a numerical tolerance.Keller derives this relation for a medium consisting of a square lattice of monodisperse disks [42].
The size invariance, however, is not the case for the results with uniform Dirichlet boundary conditions.A monotonic convergence towards the prediction of periodic boundary conditions is observed.The direction of the convergence agrees with our remark on the effective permeability prediction of vector and scalar potential formulations for vanishing corrector functions.Without and with phase-interchange, the apparent permeability predictions for the VE with 1 × 1 tiling are 45.69 0 and 147.7 0 , respectively, using scalar potential formulation.The VE with 64 × 64 tilings improves to give 3.220 0 and 135.5 0 , respectively.Using vector potential formulation, without and with phase-interchange, the apparent permeability predictions for the VE with 1 × 1 tiling are 1.692 0 and 5.471 0 , respectively.The VE with 64 × 64 tilings improves to give 1.844 0 and 77.63 0 , respectively.For the phase contrast  i ∕ m = 250∕1, Dirichlet boundary conditions give better results for vector potential formulation.In comparison, for  i ∕ m = 1∕250, scalar potential formulation predictions are the ones closer to the effective properties.This demonstrates the difficulty in locating the RVE size using uniform Dirichlet boundary conditions.In light of these computations, it can be asserted that the effective permeability predictions with uniform Dirichlet boundary conditions constitute an upper bound to periodic boundary condition predictions when used with the scalar potential formulation.Otherwise, they comprise a lower bound.The computational results are further bounded by the upper (Voigt) and lower (Reuss) bounds which, considering the current material and geometrical properties, are  ⋆ V ≃ 75.70 0 and  ⋆ R ≃ 1.426 0 , respectively.With phase-interchange giving  i ∕ m = 250∕1 these read  ⋆ V ≃ 175.3 0 and  ⋆ R ≃ 3.303 0 .It is worth noting that uniform Dirichlet boundary conditions suppress fluctuation fields only at the unit cell boundaries but not at the interior domain.Fig. 7 shows the influence of boundary conditions on magnetic flux distribution for periodic volume elements at which the high permeability phase occupies the boundary causing high property fluctuations at this region.When subjected to periodic boundary conditions, despite phase property fluctuations at the edges, the solution remains periodic within the primitive unit cell.For uniform Dirichlet boundary conditions, however, this results in the formation of a boundary layer with unusually high magnetic flux concentration at the high permeability phase.This boundary layer is diminished fast towards the interior regions of the volume element such that the field distributions for the periodic boundary conditions are recovered.This implies that the accuracy can be improved by reducing the relative size of the boundary layer, hence by increasing the VE size, as shown above.
Having shown the invariance properties of the solutions with periodic boundary conditions in conjunction with scalar and vector potential formulations, we continue the rest of this section by investigating the inclusion volume fraction effect on the effective magnetic permeability of the composite.Here, we use the primitive unit cell given in row one of Fig. 2. We vary the inclusion volume fraction within the interval of  i ∈ [0, 0.55].
A matrix with a square and periodic arrangement of nonoverlapping disks possesses magnetostatic isotropy.For the analytical solution quantifying the influence of the included content on composite's effective isotropic permeability in periodic square disk arrangements, we adopt the truncated reiteration of the series expansion solution of Godin [32] given in [43] with where the coefficients as functions of As demonstrated in Fig. 8 our computational results with scalar and vector potential formulations and periodic boundary conditions show an excellent agreement with the analytical solution given in Eq. ( 58).As before, Keller's phase-interchange relation encapsulated in the expression  ⋆ ( 1 ,  2 )  ⋆ ( 2 ,  1 ) =  1  2 is satisfied within a numerical tolerance [42].
In understanding the influence of phase contrast and disk volume fraction on the effective magnetic properties of the composite, it is helpful to analyze the magnetic induction field intensities and lines.Considering disk volume fractions of 0.05, 0.30, and 0.55, Fig. 9 demonstrates the results for the contour plots for the norm of the magnetic field  as well as the associated flux lines for the horizontally applied macroscopic magnetic field with  = [1, 0] ⊤ .Due to the microstructural symmetry, identical results are carried out for  = [0, 1] ⊤ , however, with a ∕2 in-plane rotation.As shown, the field lines are attracted towards the higher permeability material.For (a)-(c), the disconnected inclusions constitute the high permeability phase, whereas for (d)-(f), it is the continuous matrix.This explains the observed intensity difference in the flux fields.For disk volume fraction of 0.30, the magnetic flux field in the neighborhood of circular disks is comparable to Stratton's electric field solution for dilute high permittivity inclusion-low permittivity matrix systems [44,45].

Random disk arrangements
Fig. 10 demonstrates effective permeability  ⋆ predictions as a function of the VE size considering both scalar and vector potential formulations as well as nonoverlapping and overlapping random disk arrangements with 0.30 inclusion volume fraction.Here the sizes of the pixelated VEs are represented in edge sizes in pixels for {128, 256, … , 2048}, which are given as powers of two.Phase contrasts of  i ∕ e = 250∕1 and  i ∕ e = 1∕250 with phase-interchange are considered.For the periodic boundary conditions, the observed size invariance concerning the selected VE size is not the case in contrast to the regular disk arrangements.This is due to the element of randomness.Still, the relative difference between the periodic boundary condition-based predictions with vector and scalar potential formulations is below 1.3% for cases without and with phase interchange.A convergence pattern towards the corresponding effective permeability is observed in the average predictions (ensemble average) for all vector and scalar potential formulations and boundary conditions with increasing VE size.While the mean values are converging, corresponding magnitudes of the standard deviations hence the standard errors of the mean, decrease.The  convergence is strictly monotonic for Dirichlet and Neumann boundary conditions.The corresponding trends follow those indicated for size dependence investigations for the periodic regular square disk arrangements.This time, random boundary property fluctuations are responsible for the observed trends.The fastest convergence is observed in the case of periodic boundary conditions.
As a result of these computations, the effective permeability of the nonoverlapping random disk arrangements for 0.30 disk volume fraction is computed as [2.028±0.014]0 and [124.8±0.8]0 for  i ∕ e = 250∕1 and  i ∕ e = 250∕1, respectively.These are respectively [2.230±0.043]0 and [113.4 ± 2.1] 0 for the overlapping disks case.Although there is relatively low disk interaction at low volume fractions, the effective permeability prediction for the overlapping disk systems is higher for  i ∕ e = 250∕1.As before, the computational results are bounded by the upper (Voigt) and lower (Reuss) bounds which, considering the current material and geometrical properties, are  ⋆ V ≃ 75.70 0 and  ⋆ R ≃ 1.426 0 , respectively.With phase-interchange giving  i ∕ m = 250∕1 these read  ⋆ V ≃ 175.3 0 and  ⋆ R ≃ 3.303 0 .For the VE size of 2048 pixels (2 11 ), the (symmetric) off-diagonal terms of the magnetic permeability tensor are generally two to three orders of magnitude smaller than the normal components, whereas, for smaller VE sizes, their relative magnitude may be more significant.However, in-plane isotropy in the physical properties is anticipated for random disk systems.This should give a spherical two-dimensional tensor with only nonzero components being the main diagonals equal to each other.This is quantified using the indicated anisotropy index, which materializes the ellipse eccentricity as demonstrated in Fig. 11.
To this end, the eigenvalues of the computed 2 × 2 permeability matrix are determined and ordered.It is shown that with increasing VE size, the mean anisotropy index gets closer to zero.This implies a reduction in the in-plane directional dependence of the effective magnetic response.
At the VE size of 2 11 = 2048 pixels, the relative standard deviations in the computed effective permeabilities reduce below 2%.Considering the 15 random realizations used in computations, this corresponds to a relative standard error of the sum below 0.5%.The relative change of the mean effective permeability compared to the calculation for the VE size of 1024 pixels is below 1%.Moreover, this size provides considerably decreased directional dependence in the material's average in-plane magnetic permeability, with the corresponding anisotropy index falling below 0.2.With these observations, even if for Dirichlet and Neumann boundary conditions, sufficient proximity to the effective solution is not met even at this size, RVE size is determined as 2048 pixels, and the remaining analysis is conducted using this size and periodic boundary conditions.An RVE size, according to the definition of Hill, i.e., a size devoid of the Neumann and Dirichlet boundary condition effects as far as they are uniform, is to be much larger.A Comparison to Effective Medium Approximations: Effectivemedium approximations, which use simple microstructural information such as phase volume fraction and shape, prove helpful in revealing the qualitative trends in effective physical properties of the composites as a function of phase concentration.In the current study, we will consider three effective-medium approximations, namely, the Maxwell, the self-consistent, and differential effective-medium approximations  In (a) and (c), nonoverlapping disk system results are given, whereas in (b) and (d), overlapping disks are considered.The results for vector and scalar potential formulations with uniform Dirichlet and periodic boundary conditions are depicted.The upper (Voigt) and lower (Reuss) bounds correspond to  ⋆ V ≃ 75.70 0 and  ⋆ R ≃ 1.426 0 , for the phase contrast  i ∕ m = 250∕1, respectively, and  ⋆ V ≃ 175.3 0 and  ⋆ R ≃ 3.303 0 , for the phase-interchanged case.Data are presented as mean value μ (markers) and standard deviation σ from the analysis of 15 random realizations.For each volume fraction corresponding standard error of the mean SEM can be computed with SEM = σ∕ √ 15.Lines merge discrete computational results.
for which the −dimensional generalizations are known [26].For convenience, we give these −dimensional generalizations in Appendix ''Effective-Medium Theories: -dimensional Generalizations'' whereas, in the following, we present their two-dimensional specializations, i.e., for  = 2. Starting with the two-dimensional form of the Maxwell approximation given in Eq. (84) for the effective permeability  ⋆ M reads Noting that  1 = 1 −  2 , the two-component self-consistent approximation [46,47] for the effective permeability  ⋆ SC given in Eq. (85) can be reiterated for two-dimensions as follows Finally, the two-dimensional differential effective-medium approximation for the effective permeability  ⋆ DEM given in Eq. ( 86) is given with the following implicit relation The −dimensional generalizations of the effective-medium approximations given in Eqs. ( 84), ( 85) and (86), satisfy the indicated Voigt and Reuss property bounds as well as the limiting processes: For the considered nonoverlapping and overlapping random disk configurations and the selected RVE size of 2048 pixels, Figure Fig. 12 depicts the effective permeability predictions as a function of the inclusion volume fraction without and with phase-interchange in a lin-log plot.As expected, increasing high permeability phase content monotonically increases the effective magnetic permeability.The trend of the effective permeability predictions for the overlapping disk systems in the current lin-log plot resembles an unsymmetric sigmoid function whose inflection point at which the corresponding highest rate of change occurs is located in the proximity of the aforementioned percolation threshold  P ≃ 0.68.A similar trend is obtained by Helsing in the computation of effective electric properties of biphasic composites with randomly overlapping disks embedded in a homogeneous infinite medium [48].Considering the inclusion volume fractions up to 0.55, which is the range applicable to the nonoverlapping random disk microstructures, the effective permeability in overlapping disk systems is calculated to be up to 33% higher. 5The difference is even more drastic and reaches up to 56% once the results for overlapping random disk systems and regular square disk systems are compared.This is explained via the inclusion interaction, which is minimal in regular systems since the inter-disk distances in the considered monodisperse system are maximized for a selected disk volume fraction.For the case of phase interchange, identical results are valid, whereas this time, the smallest effective permeability belongs to the overlapping disk systems.For comparison purposes, the curves for the upper (Voigt) and the lower (Reuss) bounds as well as the two-dimensional Maxwell, self-consistent, and differential effective-medium approximations, respectively, given by Eqs.(64), ( 65) and (66) are also provided in Fig. 12.As anticipated, the Voigt and Reuss mixture rules bound the numerical solutions and the effective-medium approximations; nevertheless, they are far from being stringent.As given in Ref. [26], the Maxwell approximation for  2 ≥  1 and  1 ≥  2 coincide with the optimal lower and upper Hashin Shtrikman bounds, respectively.Our results are in complete agreement with this statement where for the phase contrast of  i ∕ m = 250∕1, the Maxwell approximation constitutes the optimal lower bound.For the inverted system, it acts as the optimal upper bound to all other predictions.The Maxwell Garnet approximation agrees well with Godin's analytical solution for the periodic regular square disk arrangements with a maximum difference of 5% for both phase contrasts.This is to show that the Maxwell effective-medium approximation preserves validity beyond the dilute regime with  i > 0.05.Nevertheless, this validity is limited to well-separated disk arrangements with the least interaction potential.Unlike other effective-medium approaches, the trend in self-consistent approximation shows sigmoid function features similar to the curve of the nonoverlapping disks.Nevertheless, it is symmetric with an inflection point at  i = 0.5.Thus, its validity, for the considered random monodisperse disk systems for which the phase-interchange symmetry is not applicable, is only up to the dilute limit and for  i > 0.90.The inflection point  i = 0.5 corresponds to the built-in percolation threshold prediction for the inclusions in the self-consistent approximation.The Maxwell and the differential effective-medium approximations consider that the phases remain connected until the other phase fills the space [26].Within the considered volume fraction range, the differential effective-medium approximation gives results with up to a 13% difference from our computations for the nonoverlapping random disk systems.For overlapping random disk systems with  1 ≥  2 , a correction [1 − sin(  1 ) exp(− 1 )] to the straight line [ 1 ∕ 2 ]  1 in lin-log plot results in a simple three-parameter expression given in Eq. (67) which gives a remarkably better overall agreement with the computations than the considered effective medium approximations, especially for inclusion volume fractions up to  1 = 0.6.Beyond  1 = 0.6, the observed maximum relative error is around 30%, around which occurs  1 = 0.8.Moreover, Eq. (67) satisfies the indicated Voigt and Reuss property bounds as well as the limiting processes: Finally, Keller's phase-interchange relation  ⋆ ( 1 ,  2 )  ⋆ ( 2 ,  1 ) =  1  2 is satisfied for the considered random disk arrangements within a numerical tolerance as before in agreement with Mendelson's generalization of Keller's relation [49] to two-dimensional two-phase composites and considering the macroscopic isotropy in permeability of the two-dimensional random disk arrangements.
Figs. 13 and 14 provide complementary information regarding the magnetic flux field and flux lines for several selected computations.In Fig. 13, the distribution of magnetic flux densities is given for a horizontally applied unitary macroscopic field of M  = [1, 0] ⊤ for nonoverlapping disk system.In the figure, disk volume fractions of 0.30 and 0.55 are considered with phase inversion.As indicated previously, 0.55 is the saturation volume fraction of the microstructure generation of nonoverlapping disk systems using a sequential adsorption process.It is observed that, like in the case of a regular disk system, the magnetic flux lines are attracted towards the high permeability phase.In the nonoverlapping disk system, the disk interaction increases with increasing high permeability disks.For the phase-interchanged system, at 0.30 volume fraction, the matrix phase percolation effectively hosts the magnetic flux lines.
As opposed to nonoverlapping disk arrangements, overlapping disks allow the formation of long-range connections among inclusions.This is referred to as phase percolation; see, e.g., Fig. 5.In two-dimensions, the percolation threshold is  p ≃ 0.68 [50] whereas for three-dimensions and with spheres it is  p ≃ 0.29 [51].For further information, the reader is referred to [26, p.67, p.123] and the references therein.In Fig. 14, the distribution of magnetic flux densities is given for a horizontally applied unitary macroscopic field of M  = [1, 0] ⊤ for an overlapping disk system.In the figure, disk volume fractions of 0.30 and 0.70 are considered with phase contrasting, where 0.70 is selected as a volume fraction slightly above the percolation threshold in two dimensions.It is shown that at 0.30 of phase volume fraction, which is below the percolation threshold, the behavior is similar to the nonoverlapping disk system, such that long-range high-intensity flux lines are not developed.Nevertheless, the overlapping provides the formation of higher hot spot regions.With increasing volume fraction to 0.70, the continuous pathways starting from the left face and reaching the right face of the volume element are developed.These are the percolating phase paths at which long-range connectivity is developed.The magnetic flux lines tend towards these percolating high permeability phase paths.

Conclusion
Inspired by the displacement method of Lukassen et al. [20], we developed an asymptotic homogenization method with a finite element method in the computation of the initial effective magnetic permeability of composites.We considered vector and scalar potential formulations with periodic and uniform Dirichlet boundary conditions assuming linear magnetostatics.We showed that the vector/scalar potential formulation yields the upper (Voigt)/lower (Reuss) bound in the computed effective magnetic permeability for completely suppressed perturbation fields.Using the concept of the eccentricity of an ellipse, we proposed an anisotropy index to quantify the degree of directionality in the material's permeability.We applied the developed framework to various monodisperse disk arrangements in two dimensions.We computationally demonstrated that while periodic boundary conditions with both vector and scalar potential formulations provide unit cell size-invariance for regular periodic material systems, the predictions with the uniform Dirichlet boundary conditions are size-dependent.With increasing size, the results converge to the effective behavior where the vector potential formulation results merge from lower magnitudes, and scalar potential formulation results converge from upper.For random material systems, convergence to effective properties is much faster for periodic boundary conditions.
Our results show that even for the selected limited class of material systems constituting monodisperse disk configurations, neither the analytical bounds nor the selected popular effective-medium approximations are universally applicable.Although the effective-medium approximations reveal qualitative trends in the physical properties of random systems, they fail to provide accurate quantitative predictions.This is because they only use volume fractions and inclusion shapes as descriptors.Thus, their applicability is generally limited to conditions with low contrast of properties and low inclusion volume fractions.This constitutes an obstacle in front of their generalization to random composites.On the other hand, the proposed computational asymptotic homogenization with finite element method remedies this limitation by providing an accurate account of the underlying phase geometry, which relates to the shape and size of constituent phases, and topology, which describes constituent connectivity.Using our numerical computations, we proposed a simple three-parameter approximation for the effective permeability of overlapping random disk systems.This shows a remarkable predictive capability for inclusion volume fractions up to 0.60.band are studied from 0 to 1 with a volume fraction step size of 0.05, which required 21 model generations (see Fig. E.16).The magnetic permeabilities of phases 1 and 2 are assumed to be  1 = 250  0 and  2 =  0 , leading to a phase contrast of  1 ∕ 2 = 250∕1.
As depicted in Fig. E.17(a), our computations for scalar and vector potential formulations with the application of periodic boundary conditions are in excellent agreement with each other, as well as the analytical solutions.More specifically, the magnetic permeability computed along the direction that matches that of the band agrees with the (upper) Voigt bound.In contrast, the transverse direction magnetic permeability corresponds to the (lower) Reuss bound.This directionality in the magnetic properties can be quantified using the proposed anisotropy index, whose plots computed for the stacked-layer microstructures are depicted in Fig. E.17(b) for various phase contrast magnitudes as a function of phase volume fraction.As anticipated, the index plots are symmetric concerning the point of phase balance, that is,  1 = 0.5.As the phase property contrast increases, rapid growth in the anisotropy index occurs even for small phase volume fraction differences from zero.
For the considered two-dimensional banded system as well, the phase-interchange relation encapsulated in the expression  ⋆ ( 1 ,  2 )  ⋆ ( 2 ,  1 ) =  1  2 is satisfied in agreement with Mendelson's generalization of Keller's relation [49] to two-dimensional two-phase composites and considering the principal axes of the permeability tensor for which the macroscopic tensor components are computed.

Fig. 1 .
Fig. 1.Schematic description of macroscopic and microscopic domains, respectively denoted by  and , and associated coordinates M  and  in asymptotic homogenization for a periodic composite material. corresponds to the periodic unit cell composed of two constituent domains  1 and  2 .The scale separation parameter 0 <  ≪ 1 controls the fineness of the microstructure and links domain coordinates M  and  with  = M ∕.

Fig. 3 .
Fig. 3. Example images of the voxel-based random microstructure realizations for the inclusion volume fraction of  i = 0.30 for nonoverlapping (top) overlapping (bottom) disks.The microstructures are generated using the indicated random sequential inhibition process under periodic boundary conditions.The VE sizes are 128, 256, 512, 1024, and pixels from left to right, respectively.For convenience, the image size is kept constant while demonstrating the VEs.

Fig. 4 .
Fig. 4. Example images of the voxel-based random microstructure realizations considering inclusion volume fractions of  i = {0.15,0.25, … , 0.55}, for nonoverlapping disks.The VE size is 2048 pixels.The microstructures are generated using the indicated random sequential inhibition process under periodic boundary conditions.

Fig. 5 .
Fig. 5. Example images of the voxel-based random microstructure realizations considering inclusion volume fractions of  i = {0.05,0.10, … , 0.95}, for overlapping disks.The VE size is 2048 pixels.The microstructures are generated using the indicated random sequential inhibition process under periodic boundary conditions.

Fig. 6 .
Fig. 6.Effective permeability predictions for the  ×  tilings of the unit cell with the inclusion volume fraction of  i = 0.30 given in row 3 of Fig. 2 for  = {2, 4, 8, 16, 32} for (a)  i ∕ e = 250∕1 and (b)  i ∕ e = 1∕250.The results for vector and scalar potential formulations with uniform Dirichlet and periodic boundary conditions are depicted.The upper (Voigt) and lower (Reuss) bounds correspond to  ⋆V ≃ 75.70 0 and  ⋆ R ≃ 1.426 0 , for the phase contrast  i ∕ m = 250∕1, respectively, and  ⋆ V ≃ 175.3 0 and  ⋆ R ≃ 3.303 0 , for the phase-interchanged case.The effective permeability predictions with vector and scalar potential formulations with periodic boundary conditions give 1.850 0 and 135.1 0 , considering the cases without and with phase-interchange, respectively.

Fig. 7 .
Fig. 7. Magnetic flux vector norm contour plots for 1 × 1, 2 × 2, 4 × 4, and 8 × 8 primitive unit cell tilings with the inclusion volume fraction of  i = 0.30 for (top) periodic and (bottom) uniform Dirichlet boundary conditions with scalar potential formulation for the horizontally applied macroscopic magnetic field with  = [1, 0] ⊤ .The results correspond to  i ∕ e = 250∕1.For demonstration purposes, the intervals [min, max] of the contour plots are taken as [0.026, 3.094], which corresponds to the interval of the solution with periodic boundary conditions.The uniform Dirichlet boundary conditions result in a magnetic flux distribution with a maximum of 643.3 corresponding to the dark gray region at the boundary inclusions.These results are normalized with respect to  0 .An effective permeability of 1.850 0 is predicted with vector and scalar potential formulations with periodic boundary conditions.

Fig. 8 .
Fig. 8. Effective permeability predictions for the primitive unit cell given in row 1 of Fig. 2 as a function of inclusion volume fraction for (a)  i ∕ e = 250∕1 and (b)  i ∕ e = 1∕250.The numerical results correspond to scalar potential formulation with periodic boundary conditions.Since Keller's phase-interchange relation  ⋆ ( 1 ,  2 )  ⋆ ( 2 ,  1 ) =  1  2 given in [42] is satisfied multiplying individual components of the plotted resultant effective relative permeability vectors for (a) and its phase-interchanged version (b) always yields 1 × 250.

C
.Soyarslan et al.

Fig. 10 .
Fig. 10.Effective permeability  ⋆ predictions for random disk arrangements with the inclusion volume fraction of  i = 0.30 for various VE sizes represented in terms of the powers  of 2. Thus,  = 7, 8, … , 11 corresponds to a VE size of 128, 256, … , 2048 pixels.The given predictions are for (a) and (b)  i ∕ e = 250∕1 and (c) and (d)  i ∕ e = 1∕250.In (a) and (c), nonoverlapping disk system results are given, whereas in (b) and (d), overlapping disks are considered.The results for vector and scalar potential formulations with uniform Dirichlet and periodic boundary conditions are depicted.The upper(Voigt)  and lower (Reuss) bounds correspond to  ⋆ V ≃ 75.70 0 and  ⋆ R ≃ 1.426 0 , for the phase contrast  i ∕ m = 250∕1, respectively, and  ⋆ V ≃ 175.3 0 and  ⋆ R ≃ 3.303 0 , for the phase-interchanged case.Data are presented as mean value μ (markers) and standard deviation σ from the analysis of 15 random realizations.For each volume fraction corresponding standard error of the mean SEM can be computed with SEM = σ∕ √ 15.Lines merge discrete computational results.

Fig. 11 .√ 1 − 15 .
Fig. 11.The effective anisotropy index 0 ≤  ⋆  < 1 predictions for the permeability tensor for random disk arrangements with the inclusion volume fraction of  i = 0.30 for various VE sizes represented in terms of the powers  of 2. Thus,  = 7, 8, … , 11 corresponds to a VE size of 128, 256, … , 2048 pixels.The given predictions are for (a) and (b)  i ∕ e = 250∕1 and (c) and (d)  i ∕ e = 1∕250.In (a) and (c), nonoverlapping disk system results are given, whereas in (b) and (d), overlapping disks are considered.The anisotropy index gives the eccentricity  of the 2 × 2 matrix representation of the symmetric effective permeability tensor such that  ⋆  = √ 1 −  ⋆ 2 ∕ ⋆ 1 where  ⋆  for  = 1, 2 are the ordered eigenvalues of the effective permeability tensor.Data are presented as mean value μ (markers) and standard deviation σ from the analysis of 15 random realizations.For each volume fraction corresponding standard error of the mean SEM can be computed with SEM = σ∕ √ 15.Lines merge discrete computational results.

Fig. 12 .
Fig. 12. Effective permeability predictions for the considered nonoverlapping and overlapping random disk configurations and the selected RVE size of 2048 pixels as a function of the inclusion volume fraction for (a)  i ∕ e = 250∕1 and (b)  i ∕ e = 1∕250.The numerical results correspond to scalar potential formulation with periodic boundary conditions.The simple three-parameter expression  ⋆ • =  e [ i ∕ e ]  i [1 − sin(  i ) exp(− i )] given in Eq. (67) for  i ≥  e , shows a good agreement especially for volume fractions up to  i = 0.6.In agreement with Mendelson's generalization of Keller's relation, our numerical predictions satisfy the phase-interchange relation, which gives  ⋆ ( 1 ,  2 )  ⋆ ( 2 ,  1 ) =  1  2 given in[42].The abbreviations SC and DEM stand for self-consistent and differential effective-medium approximations.
and (12)are subjected to prescribed Dirichlet and Neumann boundary conditions.For vector potential formulation, these respectively read   = Ã on    and   ×   = t on    , where   is the surface unit normal.In the same manner, we obtain   = ρ on    and   ⋅  = t on    , respectively, for scalar potential formulation.For the domain boundaries, the following conditions apply for vector potential formulation  ←   .The notations   + =   ( + ) and   − =   ( − ) apply where [ + −  − ] are position vectors of periodically located nodes on the surface .