Design of Particulate-Reinforced Composite Materials

A microstructure-based model is developed to study the effective anisotropic properties (magnetic, dielectric or thermal) of two-phase particle-filled composites. The Green’s function technique and the effective field method are used to theoretically derive the homogenized (averaged) properties for a representative volume element containing isolated inclusion and infinite, chain-structured particles. Those results are compared with the finite element approximations conducted for the assumed representative volume element. In addition, the Maxwell–Garnett model is retrieved as a special case when particle interactions are not considered. We also give some information on the optimal design of the effective anisotropic properties taking into account the shape of magnetic particles.


Introduction
Temperature as well as magnetic and electric fields in composite materials are of interest in many engineering applications. To use them effectively in modern constructions, it is necessary to predict the effective homogenized properties. Theoretically, for particle-filled, two-phase composites, their homogenized, effective physical properties may be derived in manner similar to the linear case. The methods for searching for effective magnetic permeabilities, dielectric permittivities or thermal conductivities are analogous and they may be calculated using the same approach. Although the similarity of those three fields was noticed a long time ago, the models for the prediction of the effective composite properties are built separately for each class of problems and in addition with the use of various simplified hypothesis (assumptions).
It is worth pointing out that the analysis of an incompressible viscous fluid flow through a porous medium can be described using analogous equations such as the ones mentioned above for the two-phase composites. On the macroscopic scale, flow through the porous material is governed by Darcy's law having the permeability tensor K H , so that the methods, which are analogous to the ones discussed herein may also be successfully applied to the fluid flow problems.
There are different homogenization approaches and they can be divided into three classes, i.e., direct, indirect and variational methods. Direct methods are based on the volume average of field quantities and they can be performed by a numerical procedure, usually the finite element method (FEM) or the boundary element method (BEM). Indirect homogenization follows the idea of the equivalent inclusion method based on Eshelby's eigenstrain solution. In this area, different variants of solutions are developed and they may be divided into: the self-consistent schemes, the Mori-Tanaka method and the differential method. The variational approach can give upper and lower bounds of the effective properties. Monographs describing different homogenization methods in details have been written by Mura [1], Nemat-Nasser and Hori [2], and Qin and Yang [3]. The review presented by Wang and Pan [4] first examines the issues, difficulties and challenges in the prediction of material To demonstrate the theoretical limitations, this study firstly aims to present the possibility of theoretical predictions of effective properties for smart materials in the 3D approach. The preliminary results are discussed in the papers by Barski and Muc [27] and Barski [28]. The presented theoretical approach is based on the use of the effective field method. Having in mind the possible applications of two-phase composites as smart materials, two separate problems are discussed in detail, i.e., the case of isolated inclusions and chain-like structures. Then, the numerical method of homogenization is proposed and the results are compared with the theoretical ones. The suggested scheme of numerical homogenization is applied to optimize the effective properties varying the shape of inclusions. As it is reported in the literature, the behavior of aggregated particles has a great influence on the appropriate modeling of the smart material deformation as it is strained and especially in view of its yield stress value. For such a class of composites, the present analysis is an introduction to the global FEM modeling of rheological deformations.

The Effective Field Method
First, let us note that in engineering practice particles embedded in a matrix are usually coated by an additional material (an interface) between constituents ( Figure 1) to enhance various properties of composites, i.e., thermal, magnetic or dielectric. Surfactants are added to alleviate the settling problem. Thus, a particulate composite is made of three different phases. Such a class of problems is analysed for instance by Kamiński [29]. Since it is very difficult to estimate the physical properties of an interface, our analysis is limited to the considerations of the two-phase material demonstrated in Figure 1b. preliminary results are discussed in the papers by Barski and Muc [27] and Barski [28]. The presented theoretical approach is based on the use of the effective field method. Having in mind the possible applications of two-phase composites as smart materials, two separate problems are discussed in detail, i.e., the case of isolated inclusions and chain-like structures. Then, the numerical method of homogenization is proposed and the results are compared with the theoretical ones. The suggested scheme of numerical homogenization is applied to optimize the effective properties varying the shape of inclusions. As it is reported in the literature, the behavior of aggregated particles has a great influence on the appropriate modeling of the smart material deformation as it is strained and especially in view of its yield stress value. For such a class of composites, the present analysis is an introduction to the global FEM modeling of rheological deformations.

The Effective Field Method
First, let us note that in engineering practice particles embedded in a matrix are usually coated by an additional material (an interface) between constituents ( Figure 1) to enhance various properties of composites, i.e., thermal, magnetic or dielectric. Surfactants are added to alleviate the settling problem. Thus, a particulate composite is made of three different phases. Such a class of problems is analysed for instance by Kamiński [29]. Since it is very difficult to estimate the physical properties of an interface, our analysis is limited to the considerations of the two-phase material demonstrated in Figure 1b. Variety of physical problems dealing with two-phase composites, understood in the sense of particles embedded in matrix, are described in the similar manner, i.e., with the use of similar physical equations. The examples of such relations are given in Table 1. For two-phased composites the general mathematical formulation of the problem as well as the methodology of the derivation of the effective material properties are presented in the Appendix. The general symbols, which are used in the Appendix, are explained in the first column in Table 1. However, the further numerical investigations concern the examples of the evaluation of effective magnetic permeability.
Below, two particular cases are solved to demonstrate the effectiveness and usefulness of the approach. Variety of physical problems dealing with two-phase composites, understood in the sense of particles embedded in matrix, are described in the similar manner, i.e., with the use of similar physical equations. The examples of such relations are given in Table 1. For two-phased composites the general mathematical formulation of the problem as well as the methodology of the derivation of the effective material properties are presented in the Appendix A. The general symbols, which are used in the Appendix A, are explained in the first column in Table 1. However, the further numerical investigations concern the examples of the evaluation of effective magnetic permeability.
Below, two particular cases are solved to demonstrate the effectiveness and usefulness of the approach.

Isolated Inclusions
Let the matrix material be isotropic (C 0 αβ = c 0 δ αβ ) and the isotropic inclusions be ellipsoids with the same sizes randomly oriented in space. Neglecting the pair interaction between inclusions, the analysis is reduced to the consideration of a single inclusion embedded in an infinite medium. Thus, the effective properties of the two-phase composite can be easily computed from the relation in Equation (A21) (see the Appendix A) with the use of the definition in Equation (A11). For instance, in the case of spherical inclusions (so-called spherical symmetry), the operator A 0 λµ takes the following form: the homogenized isotropic properties in Equation (A21) coincide with the well-known Maxwell-Garnett formula [5,6].

Chain of Inclusions
Using the above methodology, the effective properties of two-phase composites may also be found for materials with regular lattices of identical inclusions. In this case, the function Ψ(x, x ) in Equation (A16) depends on the difference x-x . From that definition one can find that: For the i-th inclusion being in the regular lattice, the integral πa 1 a 2 a 3 J(a), J(a) = (1 − 0.5|x/a| 2 )(1 + |x/a|/4), |x/a| ≤ 2 0, |x/a| > 2 (3) is the volume of the intersection of two identical ellipsoids with the center separated by the vector x. Using the above definition, the numerator in Equation (2) can be written as follows: where m is the vector of the lattice composed by the centers of inclusions. Thus, the function Ψ(x) can be expressed as follows: where the prime over the sum symbol denotes the exclusion of the term with m = 0. Inserting the above results into the relation in Equation (A17), one can find that: and the symbol p.v. means that the integral is understood in the sense of the Cauchy principal value. Finally, with using Equations (A21) and (6), one can evaluate the effective properties of the composite with a regular lattice. It is worth emphasizing that the above definition of the effective mechanical properties does not consider the physical interaction between particles but the geometrical form of the assumed elementary cells only. The analytical results of integration in Equation (6) can be obtained for specific forms of regular lattices only. For instance, such a formula can be derived for an infinite chain of spherical particles where the representative volume cell has the form of a cuboid having an infinite length in two directions ( Figure 2). In this case, the effective properties can be expressed in the following form: assuming the isotropic properties of the matrix and inclusions, and the operator A 0 is described by Equation (1).

Numerical Homogenization Strategy
All known analytical methods are valid under certain limitations and particular geometries or classes of structures. For metamaterials comprising conducting and possibly resonant elements, and for which the periodicity is not necessarily negligible relative to the free-space wavelength, analytical homogenization techniques are unreliable or not applicable. Our intent here is to justify a numerically based homogenization scheme based on Table 1, in which the local fields computed for one unit cell of a periodic structure are averaged to yield a set of macroscopic fields. Having computed the macroscopic fields, we can then determine the constitutive relationships between the

Numerical Homogenization Strategy
All known analytical methods are valid under certain limitations and particular geometries or classes of structures. For metamaterials comprising conducting and possibly resonant elements, and for which the periodicity is not necessarily negligible relative to the free-space wavelength, analytical homogenization techniques are unreliable or not applicable. Our intent here is to justify a numerically based homogenization scheme based on Table 1, in which the local fields computed for one unit cell of a periodic structure are averaged to yield a set of macroscopic fields. Having computed the macroscopic fields, we can then determine the constitutive relationships between the macroscopic fields, arriving at the effective electromagnetic parameters. We will see that there will be virtually no restrictions on the contents of the unit cell, nor will the unit cell necessarily need to be small in comparison with the wavelength. They have proven useful in many situations, e.g., low volume fraction of homogeneous spherical or ellipsoidal inclusions in a homogeneous host material, but they fail if the volume fraction is too high or if the inclusions are not spheres or ellipsoids. Then, a more accurate homogenization procedure has to be used, which includes all contributions of the interaction between the reinforcing particles.
For two-phase composites, a typical homogenization situation is depicted in Figure 3. It shows a 2D model of the composite material, which includes a matrix with inclusions. In the homogenization method, the structure of the two-phase composite materials is assumed to be periodic, and the unit cell, which is the minimum volume to represent the overall statistics, is defined. Here, it is assumed that one particle inclusion is located in the center of the cell. The unit cell is regarded as a homogeneous substance with the effective properties. The effective property is defined based on energy balance in the unit cell ( Figure 3). It is assumed that the original cell and the homogenized cell include equivalent energy when both unit cells are immersed in equivalent external field.

Numerical Homogenization Strategy
All known analytical methods are valid under certain limitations and particular geometries or classes of structures. For metamaterials comprising conducting and possibly resonant elements, and for which the periodicity is not necessarily negligible relative to the free-space wavelength, analytical homogenization techniques are unreliable or not applicable. Our intent here is to justify a numerically based homogenization scheme based on Table 1, in which the local fields computed for one unit cell of a periodic structure are averaged to yield a set of macroscopic fields. Having computed the macroscopic fields, we can then determine the constitutive relationships between the macroscopic fields, arriving at the effective electromagnetic parameters. We will see that there will be virtually no restrictions on the contents of the unit cell, nor will the unit cell necessarily need to be small in comparison with the wavelength. They have proven useful in many situations, e.g., low volume fraction of homogeneous spherical or ellipsoidal inclusions in a homogeneous host material, but they fail if the volume fraction is too high or if the inclusions are not spheres or ellipsoids. Then, a more accurate homogenization procedure has to be used, which includes all contributions of the interaction between the reinforcing particles.
For two-phase composites, a typical homogenization situation is depicted in Figure 3. It shows a 2D model of the composite material, which includes a matrix with inclusions. In the homogenization method, the structure of the two-phase composite materials is assumed to be periodic, and the unit cell, which is the minimum volume to represent the overall statistics, is defined. Here, it is assumed that one particle inclusion is located in the center of the cell. The unit cell is regarded as a homogeneous substance with the effective properties. The effective property is defined based on energy balance in the unit cell ( Figure 3). It is assumed that the original cell and the homogenized cell include equivalent energy when both unit cells are immersed in equivalent external field.  In an actual estimation, the solution of the Laplace equation obtained with the use of basic relation shown in Table 1 is computed by FEM. In this analysis, the potential φ is unknown, and, assuming the applied field is unidirectional at least in the cell, the boundary conditions are set as: The identical boundary conditions are formulated at each parallel boundaries of the cell shown in Figure 4 for the 2D case. Let us note that the above type of boundary conditions satisfies the periodicity of boundary conditions for the arbitrary type of regular lattices as well as for the chain of inclusions and for a single inclusion.
The identical boundary conditions are formulated at each parallel boundaries of the cell shown in Figure 4 for the 2D case. Let us note that the above type of boundary conditions satisfies the periodicity of boundary conditions for the arbitrary type of regular lattices as well as for the chain of inclusions and for a single inclusion.  Similar to the previous case for the effective field method, the FEM analysis is based on the averaging method that is carried out for the representative volume element (RVE) having the volume denoted by the symbol RVE. Thus, it is obvious that the results, understood in the sense of the average property tensor, are directly dependent on the dimensions and form of the RVE. It is worth pointing out that, even for the 2D two-phase periodic composites, the RVE may be of an arbitrary form, not necessarily rectangular as it shown in Figure 3. Since in our numerical analysis we intend to give information on the variations of the property tensor components with respect to the volume fraction p, we define it in the following way (see Figure 4): where it is assumed that the RVE has a square cross-section in the y direction. Thus, for the prescribed volume fraction p, the RVE is completely defined by the set of two parameters (geometrical ratios), i.e., gf/rp and b/rp for both 2D and 3D cases. Let us note that, for the constant volume fraction p and the constant interparticle distance gf, the geometrical dimensions of the representative cell (i.e., b and h) are uniquely determined For the selected RVE ( Figure 4) and the selected boundary conditions in the form of Equation (8) (the unidirectional external field), the average intensity and flux of the field are defined by: where TN denotes the total number of elements in the FEM mesh. Using the above relations, it is possible to compute four components of the average intensity and flux of the field for two types of boundary conditions demonstrated in Figure 4 (the 2D problem) or nine for the 3D analysis. We do not know in advance how many nonzero components in the property matrix C 1  occur. Therefore, for the linear problem, all components of the property matrix C 1  (nine for the 3D problem) can be derived directly from basic relationships presented in Table 1 where the appropriate components of the vectors E and D are replaced by their average values evaluated in the local cell (RVE). Similar to the previous case for the effective field method, the FEM analysis is based on the averaging method that is carried out for the representative volume element (RVE) having the volume denoted by the symbol Ω RVE . Thus, it is obvious that the results, understood in the sense of the average property tensor, are directly dependent on the dimensions and form of the RVE. It is worth pointing out that, even for the 2D two-phase periodic composites, the RVE may be of an arbitrary form, not necessarily rectangular as it shown in Figure 3. Since in our numerical analysis we intend to give information on the variations of the property tensor components with respect to the volume fraction p, we define it in the following way (see Figure 4): p = πr 2 p 2br p 1 + g f /r p for 2D and p = 2πr 3 p 3b 2 r p 1 + g f /r p for 3D (9) where it is assumed that the RVE has a square cross-section in the y direction. Thus, for the prescribed volume fraction p, the RVE is completely defined by the set of two parameters (geometrical ratios), i.e., g f /r p and b/r p for both 2D and 3D cases. Let us note that, for the constant volume fraction p and the constant interparticle distance g f , the geometrical dimensions of the representative cell (i.e., b and h) are uniquely determined. For the selected RVE ( Figure 4) and the selected boundary conditions in the form of Equation (8) (the unidirectional external field), the average intensity and flux of the field are defined by: where TN denotes the total number of elements in the FEM mesh. Using the above relations, it is possible to compute four components of the average intensity and flux of the field for two types of boundary conditions demonstrated in Figure 4 (the 2D problem) or nine for the 3D analysis. We do not know in advance how many nonzero components in the property matrix C 1 αβ occur. Therefore, for the linear problem, all components of the property matrix C 1 αβ (nine for the 3D problem) can be derived directly from basic relationships presented in Table 1 where the appropriate components of the vectors E and D are replaced by their average values evaluated in the local cell (RVE).
For the non-linear problem (understood in the sense of the non-linear relation, Table 1), the components of the property matrix are computed by the comparison of the magnetic energy of the homogenized and of the original (Figure 4) unit cell combined with the Newton-Raphson method. The energy is represented as follows: However, it is necessary to know in advance the D-E characteristic curve for the inclusions. Commonly, the inclusion is modeled as isotropic material but for the unit cell the local property matrix can possess anisotropic properties as it is shown for example in Section 2. Figure 5 represents the example the FEM analysis conducted for two-phase ferromagnetic composites made of phases having magnetic properties-the magnetostatic problem. The plots demonstrate the distribution of the magnetic flux density, the magnetic field and the potential φ as the external magnetic field is applied at the x direction-the boundary conditions in Equation (8) For two-phase composites, the effective property tensor has analogous properties to, e.g., the tension modulus in elasticity: if the behavior is identical in three perpendicular directions, then it is isotropic. This conclusion points out a limitation of the use of constant second order tensors for the description of behavior. Indeed, e.g., for magnetism, many experimental observations reveal that cubic single crystals are not magnetically isotropic (see for instance Webster [30] for iron and nickel or Wang et al. [31] for Terfenol-D). In fact, in the experiments, the chains of particles do not have to be aligned in the direction of the external field, namely E y . We have some more compact aggregates. The microscopic analysis on the structure of the two-phase composites revealed that there were aggregates forming rather than chains of spheres that can be approximated by ellipsoids, stripes or cylinders. Therefore, it is interesting to verify the correctness of the introduced FE model in cases when the external field is rotated with respect to coordinates defining RVE in order to consider the relationship between the orientation of inclusions and to explore symmetries in the constitutive relations in Equation (1) and their relevance to the homogenized composite medium. In the case of a generalized anisotropic structure for which principal axis of external fields and the unit cell do not coincide, the property tensor should satisfy the classical transformation rules of the second rank tensors. To interpret and investigate those effects, let us analyze the form of the property (permeability) matrix for the 3D ring structure shown in Figure 6. The geometry of unit cells is defined in the cylindrical coordinate system but the external field is directed along the line joining the center of the central spherical inclusion and the center of the ring curvature.
components of the property matrix are computed by the comparison of the magnetic energy of the homogenized and of the original (Figure 4) unit cell combined with the Newton-Raphson method. The energy is represented as follows: However, it is necessary to know in advance the D-E characteristic curve for the inclusions. Commonly, the inclusion is modeled as isotropic material but for the unit cell the local property matrix can possess anisotropic properties as it is shown for example in Section 2. Figure 5 represents the example the FEM analysis conducted for two-phase ferromagnetic composites made of phases having magnetic properties-the magnetostatic problem. The plots demonstrate the distribution of the magnetic flux density, the magnetic field and the potential ϕ as the external magnetic field is applied at the x direction-the boundary conditions in Equation (8). As can be seen, the variations of the potential inside the local cell (Figure 5a)  For two-phase composites, the effective property tensor has analogous properties to, e.g., the tension modulus in elasticity: if the behavior is identical in three perpendicular directions, then it is isotropic. This conclusion points out a limitation of the use of constant second order tensors for the description of behavior. Indeed, e.g., for magnetism, many experimental observations reveal that cubic single crystals are not magnetically isotropic (see for instance Webster [30] for iron and nickel or Wang et al. [31] for Terfenol-D). In fact, in the experiments, the chains of particles do not have to be aligned in the direction of the external field, namely Ey. We have some more compact aggregates. The microscopic analysis on the structure of the two-phase composites revealed that there were aggregates forming rather than chains of spheres that can be approximated by ellipsoids, stripes or cylinders. Therefore, it is interesting to verify the correctness of the introduced FE model in cases when the external field is rotated with respect to coordinates defining RVE in order to consider the relationship between the orientation of inclusions and to explore symmetries in the constitutive relations in Equation (1) and their relevance to the homogenized composite medium. In the case of a generalized anisotropic structure for which principal axis of external fields and the unit cell do not coincide, the property tensor should satisfy the classical transformation rules of the second rank tensors. To interpret and investigate those effects, let us analyze the form of the property (permeability) matrix for the 3D ring structure shown in Figure 6. The geometry of unit cells is defined in the cylindrical coordinate system but the external field is directed along the line joining the center of the central spherical inclusion and the center of the ring curvature.  The assumed form of the system of the unit cells reflects the situation as the external field is not always parallel to the cell edges. It may occur for instance for clusters of inclusion. Figure 7 shows the assumed boundary conditions and the values of the permeability matrix terms (the magnetostatic problem). The terms cαβ (α ≠ β) are not equal to zero, which means that it may be for instance the origin of clusters and of the inclusion aggregation at the beginning of magnetization. On the other hand, it may be easily verified that the terms of the property matrix satisfy the classical transformation rule:   The assumed form of the system of the unit cells reflects the situation as the external field is not always parallel to the cell edges. It may occur for instance for clusters of inclusion. Figure 7 shows the assumed boundary conditions and the values of the permeability matrix terms (the magnetostatic problem). The terms cαβ (α ≠ β) are not equal to zero, which means that it may be for instance the origin of clusters and of the inclusion aggregation at the beginning of magnetization. On the other hand, it may be easily verified that the terms of the property matrix satisfy the classical transformation rule:  The assumed form of the system of the unit cells reflects the situation as the external field is not always parallel to the cell edges. It may occur for instance for clusters of inclusion. Figure 7 shows the assumed boundary conditions and the values of the permeability matrix terms (the magnetostatic problem). The terms c αβ (α = β) are not equal to zero, which means that it may be for instance the origin of clusters and of the inclusion aggregation at the beginning of magnetization. On the other hand, it may be easily verified that the terms of the property matrix satisfy the classical transformation rule: and θ denotes the angle of rotation which is equal to 14.11 0 .

Numerical Results
The homogenization method is applied to various 2D and 3D test problems to evaluate the distributions of the terms of the property matrix and to compare theoretical predictions with numerical ones that take into account the finite dimensions of the unit cell. In the test problems, a sample two-phase composite material is composed of an isotropic matrix and the inclusions have the following material properties: cp = 2000 and cf = 1. The analysis is conducted for 2D and 3D unit cells to test the capability and limitations of the proposed model. Similar to the previous case, the numerical model corresponds to the analysis of MR fluids.

2D Problems
Let us consider a single circular particle surrounded by a nonmagnetic carrier fluid-the planar problem in Figure 4. This is a typical homogenization problem analyzed for the MR fluids (see, e.g., Simon et al. [32] for MR fluids). However, on the contrary to the cited work, we compute the four (the planar problem) permeability matrix coefficients. The off-axis terms are equal to zero, and two others are plotted in Figure 8 for various volume fractions.

Numerical Results
The homogenization method is applied to various 2D and 3D test problems to evaluate the distributions of the terms of the property matrix and to compare theoretical predictions with numerical ones that take into account the finite dimensions of the unit cell. In the test problems, a sample two-phase composite material is composed of an isotropic matrix and the inclusions have the following material properties: c p = 2000 and c f = 1. The analysis is conducted for 2D and 3D unit cells to test the capability and limitations of the proposed model. Similar to the previous case, the numerical model corresponds to the analysis of MR fluids.

2D Problems
Let us consider a single circular particle surrounded by a nonmagnetic carrier fluid-the planar problem in Figure 4. This is a typical homogenization problem analyzed for the MR fluids (see, e.g., Simon et al. [32] for MR fluids). However, on the contrary to the cited work, we compute the four (the planar problem) permeability matrix coefficients. The off-axis terms are equal to zero, and two others are plotted in Figure 8 for various volume fractions.
As can be seen, the effective permeability c zz decreases as the interparticle distance increases and it is the highest for the highest volume fraction p = 0.5. The decrease of the effective permeability c xx is associated with the increase of the effective permeability c zz . Let us note that, for the constant volume fraction p, the variations of the ratio r p /h (or g f /r p ) results in the change of the ratio b/r p (see Equation (9)). Using the single unit cell (Figure 4), one can observe that, for the external magnetic field having the non-zero component H z only, the chains of ferromagnetic particles are completely isolated since there is no interaction at the x direction. In fact, the experiments evidently demonstrate evidently that they form clusters of different shapes (see Bossis et al. [33]). In Figure 9, a relationship between the volume fraction and the effective permeability for arbitrary assumed value of ratio r p /h is depicted. As can be observed in the case of c zz , these relations can be considered as almost linear. In the case of c xx for small values of volume fraction, they are also linear. However, for the values of volume fraction greater then p > 0.35, the effective permeability values increases exponentially. For both parameters, c zz and c xx , the values are achieved for volume fraction p = 0.5. numerical ones that take into account the finite dimensions of the unit cell. In the test problems, a sample two-phase composite material is composed of an isotropic matrix and the inclusions have the following material properties: cp = 2000 and cf = 1. The analysis is conducted for 2D and 3D unit cells to test the capability and limitations of the proposed model. Similar to the previous case, the numerical model corresponds to the analysis of MR fluids.

2D Problems
Let us consider a single circular particle surrounded by a nonmagnetic carrier fluid-the planar problem in Figure 4. This is a typical homogenization problem analyzed for the MR fluids (see, e.g., Simon et al. [32] for MR fluids). However, on the contrary to the cited work, we compute the four (the planar problem) permeability matrix coefficients. The off-axis terms are equal to zero, and two others are plotted in Figure 8 for various volume fractions.  As can be seen, the effective permeability czz decreases as the interparticle distance increases and it is the highest for the highest volume fraction p = 0.5. The decrease of the effective permeability cxx is associated with the increase of the effective permeability czz. Let us note that, for the constant volume fraction p, the variations of the ratio rp/h (or gf/rp) results in the change of the ratio b/rp (see Equation (9)). Using the single unit cell (Figure 4), one can observe that, for the external magnetic field having the non-zero component Hz only, the chains of ferromagnetic particles are completely isolated since there is no interaction at the x direction. In fact, the experiments evidently demonstrate evidently that they form clusters of different shapes (see Bossis et al. [33]). In Figure 9, a relationship between the volume fraction and the effective permeability for arbitrary assumed value of ratio rp/h is depicted. As can be observed in the case of czz, these relations can be considered as almost linear. In the case of cxx for small values of volume fraction, they are also linear. However, for the values of volume fraction greater then p > 0.35, the effective permeability values increases exponentially. For both parameters, czz and cxx, the values are achieved for volume fraction p = 0.5.

3D Problems
Now, the 3D Laplace relationships shown in Table 1 have been solved for a 3D unit cell using the FEM package. The results for the variations of the property tensor components are demonstrated in Figures 10 and 11 considering the assumption in Equation (9), which leads directly to the transversely isotropic properties of composites (i.e., cxx = cyy). The plots are drawn both for spherical (a = ax = ay = az) and spheroid (a = ax = ay ≠ az) inclusions. For spherical inclusions, the distributions of the homogenized properties for the 3D case are similar to those for the 2D case (see Figures 8 and 9). However, assuming the identical geometric ratios of the RVE, for the 3D case, the averaged values in the x direction are higher than those evaluated for the planar case, whereas the values in the z direction are almost identical. The shape of inclusions has a significant influence on the averaged values. For the increasing az/h parameter, the effective property czz decreases, and cxx increases since the x axis corresponds to the longer axis of spheroids. Figures 10 and 11 show that the proposed model provides a transversely isotropic effective permeability, whereas the Maxwell-Garnett model gives an isotropic one. The Maxwell-Garnett model always yields the same estimates for any rp/h because it is insensitive to the microstructure. Thus, the Maxwell-Garnett model cannot be used for these composites because, even when the overall volume fraction is very small, the distance between particles of the same chain is so small that particle interactions cannot be disregarded. For the low volume fractions, the Maxwell-Garnett gives very good estimations in the x directions only (Figure 11).

3D Problems
Now, the 3D Laplace relationships shown in Table 1 have been solved for a 3D unit cell using the FEM package. The results for the variations of the property tensor components are demonstrated in Figures 10 and 11 considering the assumption in Equation (9), which leads directly to the transversely isotropic properties of composites (i.e., c xx = c yy ). The plots are drawn both for spherical (a = a x = a y = a z ) and spheroid (a = a x = a y = a z ) inclusions. For spherical inclusions, the distributions of the homogenized properties for the 3D case are similar to those for the 2D case (see Figures 8 and 9). However, assuming the identical geometric ratios of the RVE, for the 3D case, the averaged values in the x direction are higher than those evaluated for the planar case, whereas the values in the z direction are almost identical. The shape of inclusions has a significant influence on the averaged values. For the increasing a z /h parameter, the effective property c zz decreases, and c xx increases since the x axis corresponds to the longer axis of spheroids. Figures 10 and 11 show that the proposed model provides a transversely isotropic effective permeability, whereas the Maxwell-Garnett model gives an isotropic one. The Maxwell-Garnett model always yields the same estimates for any r p /h because it is insensitive to the microstructure. Thus, the Maxwell-Garnett model cannot be used for these composites because, even when the overall volume fraction is very small, the distance between particles of the same chain is so small that particle interactions cannot be disregarded. For the low volume fractions, the Maxwell-Garnett gives very good estimations in the x directions only (Figure 11). For spherical inclusions, as can be observed in Figures 10 and 11, the effective field model gives much better approximations of the effective values evaluated with the use of the FE model than the Maxwell-Garnett model. However, the effectiveness of the effective field model decreases for the high particle interactions (a z /h < 0.1) and for the high volume fractions (p > 0.25). It is obvious that theoretical estimations, i.e., with the use of the Maxwell-Garnett model and the effective field model have limited applications in comparison with the FEM model since the first corresponds to the random (quasi-isotropic) structure of reinforced particles, and the second to the chain-like structure, i.e., in two directions the dimensions of the RVE tend to infinity.

Optimal Design
In the theoretical and numerical analysis, it is commonly assumed that the reinforcement particle has an ideal spherical form. However, as demonstrated in Figures 10 and 11, the shape of the particle can significantly affect the effective material properties. Therefore, the development of optimized multifunctional composite materials becomes of great interest from the technological and theoretical viewpoints to all engineering fields. In general, such an approach belongs to the class of

Optimal Design
In the theoretical and numerical analysis, it is commonly assumed that the reinforcement particle has an ideal spherical form. However, as demonstrated in Figures 10 and 11, the shape of the particle can significantly affect the effective material properties. Therefore, the development of optimized multifunctional composite materials becomes of great interest from the technological and theoretical viewpoints to all engineering fields. In general, such an approach belongs to the class of Figure 11. Variations of the permeability coefficients with volume fraction p (transversely-isotropic body).

Optimal Design
In the theoretical and numerical analysis, it is commonly assumed that the reinforcement particle has an ideal spherical form. However, as demonstrated in Figures 10 and 11, the shape of the particle can significantly affect the effective material properties. Therefore, the development of optimized multifunctional composite materials becomes of great interest from the technological and theoretical viewpoints to all engineering fields. In general, such an approach belongs to the class of so-called shape optimization problems (see Muc et al. [34]). The broader description of optimization algorithms that can be applied herein is discussed by Muc [35].
This section designs such materials computationally using the method of parametric optimization. In particular, two-dimensional periodic two-phase composite materials are optimized for the optimal effective properties. To analyze the effects of different shapes, for simplicity, it is assumed that the particle is modeled as a superellipse (Figure 4) having the following form: where a and b denote the superellipse semi-axis, and n is a parameter greater than 1. Having the constant volume fraction in the RVE and varying the n value, one can observe the change of the terms in the effective property matrix.
For higher values of the n parameter (n > 10), the shape resembles a rectangle and, in this case, both components of the effective permeabilities, i.e., c yy and c xx reach their optimal values (see Figures 12 and 13). However, the optimal values of the effective properties are strongly dependent on the values of the geometrical ratios r p /h and the volume fractions p. Let us note that the maximal value of the term c zz is much higher than those plotted previously in Section 4, and the values of c zz are much lower. Therefore, it seems to be reasonable to conclude that the optimal rectangular form of the particles can prevent their aggregation in ellipsoids or cylinders instead of linear chains and in this sense, the theoretical effective field model may be applicable in the estimations of the effective properties. It is worth mentioning that the obtained optimal designs completely resemble those obtained by Guest and Prevost [36] for fluid transportation problem (the Darcy law). They concluded that the Schwartz P minimal surface is believed to be the maximum permeability structure in the 3D case. However, the authors of the cited paper assumed in advance the isotropic properties of the permeability matrix. for the optimal effective properties. To analyze the effects of different shapes, for simplicity, it is assumed that the particle is modeled as a superellipse (Figure 4) having the following form: where a and b denote the superellipse semi-axis, and n is a parameter greater than 1. Having the constant volume fraction in the RVE and varying the n value, one can observe the change of the terms in the effective property matrix.
For higher values of the n parameter (n > 10), the shape resembles a rectangle and, in this case, both components of the effective permeabilities, i.e., cyy and cxx reach their optimal values (see Figures 12 and 13). However, the optimal values of the effective properties are strongly dependent on the values of the geometrical ratios rp/h and the volume fractions p. Let us note that the maximal value of the term czz is much higher than those plotted previously in Section 4, and the values of czz are much lower. Therefore, it seems to be reasonable to conclude that the optimal rectangular form of the particles can prevent their aggregation in ellipsoids or cylinders instead of linear chains and in this sense, the theoretical effective field model may be applicable in the estimations of the effective properties. It is worth mentioning that the obtained optimal designs completely resemble those obtained by Guest and Prevost [36] for fluid transportation problem (the Darcy law). They concluded that the Schwartz P minimal surface is believed to be the maximum permeability structure in the 3D case. However, the authors of the cited paper assumed in advance the isotropic properties of the permeability matrix. (a)

Concluding Remarks
For two-phase particulate composites, the compact, unique formulations and methods of solutions (theoretical and numerical) are presented in the paper. Although numerical examples deal with the analysis of magnetic phenomena, the methods and results can be easily extended to various other classes of problems.
The properties of a heterogeneous medium (two-phase composites) made of inclusions distributed in a locally periodic way in a matrix have been derived and studied. A uniform test external field is applied on the boundary of the composite, and then the averaged fields of the particles and matrix are derived theoretically by the Green's function technique and then compared with FE results based on the numerical homogenization technique. An anisotropic effective property tensor is further provided. The effective property tensor of the composite medium is symmetric, positive definite, generally anisotropic, and depends on the microstructure both for 2D and 3D cases.
The results of the theoretical and computational (FEM) approximations are strongly dependent on the construction of the computational (representative) cell using specified dimensions and distributions of the reinforced particles, i.e., their location with respect to the centers of representative volume elements. The proposed methodology can be enriched by adding the interfacial elements between particles and matrices.
The proposed method can be successfully applied to the analysis of the non-linear problems, considering the non-linearity of the characteristic curves (e.g., B-H). From these models, it is found that the averaged property tensor components are strongly dependent on the dimensionless interparticle distance and the volume fraction.
In addition, this paper proposes a shape optimization methodology for designing multifunctional two-phase composite materials optimized for tensor property components. For the 2D problems, the optimal shape resembles rectangle with rounded edges. It is verified in this study that the optimal design based on the finite element analysis is a valid method for the output improvement of constructions.
The validity of the model should be also estimated by the inclusion of the uncertainties dealing with the shape of reinforcements or representative cells, e.g., using fuzzy sets [37].

Concluding Remarks
For two-phase particulate composites, the compact, unique formulations and methods of solutions (theoretical and numerical) are presented in the paper. Although numerical examples deal with the analysis of magnetic phenomena, the methods and results can be easily extended to various other classes of problems.
The properties of a heterogeneous medium (two-phase composites) made of inclusions distributed in a locally periodic way in a matrix have been derived and studied. A uniform test external field is applied on the boundary of the composite, and then the averaged fields of the particles and matrix are derived theoretically by the Green's function technique and then compared with FE results based on the numerical homogenization technique. An anisotropic effective property tensor is further provided. The effective property tensor of the composite medium is symmetric, positive definite, generally anisotropic, and depends on the microstructure both for 2D and 3D cases.
The results of the theoretical and computational (FEM) approximations are strongly dependent on the construction of the computational (representative) cell using specified dimensions and distributions of the reinforced particles, i.e., their location with respect to the centers of representative volume elements. The proposed methodology can be enriched by adding the interfacial elements between particles and matrices.
The proposed method can be successfully applied to the analysis of the non-linear problems, considering the non-linearity of the characteristic curves (e.g., B-H). From these models, it is found that the averaged property tensor components are strongly dependent on the dimensionless interparticle distance and the volume fraction.
In addition, this paper proposes a shape optimization methodology for designing multifunctional two-phase composite materials optimized for tensor property components. For the 2D problems, the optimal shape resembles rectangle with rounded edges. It is verified in this study that the optimal design based on the finite element analysis is a valid method for the output improvement of constructions.
The validity of the model should be also estimated by the inclusion of the uncertainties dealing with the shape of reinforcements or representative cells, e.g., using fuzzy sets [37].
It is important to emphasize that the underlying methodology of homogenization and optimization is quite general and can be applied to the design of arbitrary composite materials in different scales, i.e., micro-, meso-and even nanoscales.
The present analysis is conducted for single fields separately. In the global formulation of homogenization problems for the composite considered, the next step can be achieved by the compact formulation for coupled fields, i.e., the thermal and magneto-electric effects will be analyzed simultaneously.
Author Contributions: Aleksander Muc prepared the analytical solution of homogenization problem, wrote the text, Marek Barski performed numerical computations and prepared the graphs.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Mathematical Formulation and Derivation of the Effective (Homogenized) Material Properties
Let us consider a homogeneous medium with properties described by a tensor C 0 αβ . The medium contains the region Ω (1) (inhomogeneity) with another property tensor C αβ . The intensity E(x) and the flux D(x) of the field in the medium are described by the following system of differential equations: where q(x) is a scalar density of the field sources. The third equation is satisfied automatically if the field E α (x) is the gradient of a scalar function ϕ(x) called the potential of the field: The physical meaning of the potential ϕ(x) depends on the problem under considerations. It means the temperature field for heat transfer problems, the potential of the electric (magnetic) field for the electrostatics (magnetostatics, respectively) fields.
The tensor C αβ represents the local properties of the inclusion (Figure 1b) and is defined as the sum: Thus, if we consider a homogeneous medium with the property tensor C 0 containing a set of inclusions with the property tensor C, the system of differential Equations (A1) and (A2) may be reduced to integral equations for the fields E(x) inside the inclusions, i.e., where V(x ) is a characteristic function for the kth inclusion that is equal to one if the variable x belongs to the region Ω 1 (k) occupied by the kth inclusion and to zero if x does not belong to the domain Ω 1 (k) (see Figure 1). Here, E 0 is the external field in the medium without the inclusion (C 1 (x) = 0) by the action of the same sources of the field. The kernel K(x) is determined in the classical manner: where G(x) is the Green function for the infinite homogeneous medium with the property tensor C 0 . The Green function satisfies the equation: The equivalence between the relations in Equations (A4) and (A1)-(A3) can be proven under two additional assumptions: -With the use of the potential ϕ(x), the solution of Equation (A1) can be decomposed as follows: ϕ(x) = ϕ 0 (x) + ϕ 1 (x), where ϕ 0 (x) is the potential in the medium without the inclusion, and ϕ 1 (x) is the perturbation of the potential due to the presence of the inclusion that tends to zero when |x| → ∞ . - The potential ϕ 0 (x) satisfies the relation analogous to Equation (A1), i.e., ∇ α C 0 αβ ∇ β ϕ 0 (x) = −q(x).
The solution of this equation for an arbitrary anisotropic medium takes the following form: The simplest version of the effective field method is based on the hypothesis that the local external field E * β (x), which acts on each inclusion, is the same for all inclusions, i.e., the field E α (x) inside each inclusion can be presented in the following way: where the tensor Λ E αβ (x) is determined from the solution of the problem for isolated inclusion in the medium with the property tensor C 0 by the action of the field E * β (x), and, for ellipsoidal inclusions, it takes the following form: where A 0 αβ = A 0 1 e 1 α e 1 β + A 0 2 e 2 α e 2 β + A 0 3 e 3 α e 3 β , A 0 n = a 1 a 2 a 3 2c 0 ∞ 0 dσ (a 2 n +σ) (a 2 1 +σ)(a 2 2 +σ)(a 2 3 +σ) , n = 1, 2, 3 (A10) a 1 , a 2 , and a 3 are the semiaxes of an ellipsoid, and e n α are the unit vectors of the ellipsoid principal axes, the orientation of which is given by the normal m. For a spheroid inclusion with the semiaxes a 1 = a 2 = a, a 3 , the tensor A 0 takes the following form: A 0 αβ = A 0 1 θ αβ + A 0 3 mαm β , θ αβ = δ αβ − mαm β , A 0 1 = 1 c0 f 0 (γ), A 0 3 = 1 c0 (1 − 2f 0 (γ)), γ = a a3 > 1, f 0 (γ) = 1−g 2(1−γ 2 ) , g = γ 2 √ γ 2 −1 arctg γ 2 − 1 (A11) The effective local exciting field E * β (x) acting on the kth inclusion is the sum of the external field E 0 β (x) applied to the medium and the field induced by others surrounding inclusions. The field E β (x) inside the isolated kth inclusion being in the background matrix and caused by the action of the field E * β (x) is defined by Equation (A4) where E 0 β (x) is replaced by E * β (x). The field induced in the region Ω 1 (k) of the kth inclusion by all surrounding inclusions can also be represented using Equation (A4) in the following form: Averaging Equation (A12), assuming that the field E 0 α (x) is fixed in the problem, i.e., and identifying the conditional mean E * α (x) |x with the effective field E * α , we finally obtain the following relation: The symbol . denotes the averaging, and the function Ψ is defined as follows: p is the volume concentration of inclusions and it is equal to V(x) . The symbol · |x denotes the averaging over the realizations of the random set of inclusions by the condition x ∈ V. The solution of Equation (A15) can be expressed in the following way: Multiplying Equation (A4) by the property tensor C 0 and using the definition in Equation (A1) of the flux tensor D(x) (for the inclusion and the medium) after a set of transformation of the result, it is possible to derive the average of the flux tensor, i.e., D α (x) = C 0 αβ + P 0 αλ δ λβ + pK Ψ λµ P 0 On the other hand, using Equations (A1) and (A3), the tensor of the effective (homogenized) properties of the composite may be defined as follows: Combining Equations (A18), (A19) and (A14), one can find that the tensor of the effective properties of two-phase composite is represented by the relation: C * αβ = C αβ (x) = C 0 αβ + pP 0 αλ δ λβ + pK Ψ λµ P 0 µβ −1 (A20) Taking into account Equations (A12) and (A15) as well as assuming that the tensor C 1 βµ is independent on the x variable, the above relation takes the following form: If the ratio γ has the same order as the inclusion aspect ratio, K Ψ λµ = −A 0 λµ and the above equation is further reduced. Thus, as may be easily observed, the two-phase composite has isotropic, transversely-isotropic or anisotropic properties and they are directly dependent on the form of the operator A 0 λµ -under the assumption that both the matrix and the inclusion have isotropic properties.