Generalized effective dynamic constitutive relation for heterogeneous media:
Beyond the quasi-infinite and periodic limits

Heterogeneous media are considered a significant class of engineering materials, where their microstructures can be controlled to adapt to new operating environments and exhibit superior properties compared to raw materials. Metamaterials, a type of heterogeneous media consisting of building blocks large enough to be fabricated yet small enough to affect macroscopic material properties, have opened up new opportunities in the field of advanced materials. They exhibit unprecedented characteristics, particularly in manipulating dynamic material responses such as negative-index, subwavelength focusing, and invisibility cloaking for electromagnetic waves [1, 2, 3, 4, 5], collimators, superlenses, and inaudible cloaks for elastic (acoustic) waves [6, 7, 8, 9, 10, 11, 12, 13]. Extensive experimental and theoretical work has been performed to explore the rational designs of tailored metamaterials with desired macroscopic properties [14, 15, 16, 17, 18, 19].

Dynamic homogenization theories are powerful tools to describe and understand the behavior of heterogeneous media, by seeking to define effective constitutive relations and material characteristics that govern macroscopic wave propagation of a hypothetical homogenized medium. Willis first presented a linear elastodynamic homogenization theory based on the ensemble average technique and Green’s function of heterogeneous media [20, 21, 22, 23, 24, 25, 26, 27, 28, 29], from which it is discovered that the homogenized medium can exhibit non-classical features: 1) non-local material properties both in space and time; 2) tensorial density, i.e., not scalar; and 3) non-classical constitutive relation. In this non-classical constitutive relation, a coupling effect exists between linear momentum $\bm{\mu}$ and stress $\bm{\sigma}$ with velocity $\bm{v}$ and strain $\bm{\varepsilon}$, called Willis coupling, and it has been extended in the case of piezoelectric heterogeneous media to a coupling effect between $\bm{\mu}$, $\bm{\sigma}$, and electric displacement $\bm{D}$ with $\bm{v}$, $\bm{\varepsilon}$, and electric field $\bm{E}$, called electro-momentum coupling [30, 31, 32, 33]. Remarkably, these coupling effects occur even in the homogenization of deterministic media [34, 35]. In addition, this idea of Willis was particularized for periodic heterogeneous media based on Bloch-Floquet waves [36, 37, 38, 39, 40].

According to Willis’ book chapter [26], determining Green’s function of heterogeneous media is the biggest challenge when it comes to the dynamic homogenization theory. As Willis noted, “If the Green’s function could be determined, then all problems would be solved! Of course, this cannot be done.” To address this challenge, researchers have used the Hashin-Shtrikman principle [22, 24, 26, 28, 29, 41], which involves using Green’s function of a reference comparison medium that is homogeneous. However, previous works in the literature have constraints when determining Green’s function of the reference medium. For instance, Green’s function must satisfy certain boundary conditions, such as being zero along the boundaries (i.e., homogeneous boundary condition) or subject to the given types of boundary conditions of the target heterogeneous medium to be homogenized. Therefore, these approaches can only be suitable for infinite or quasi-infinite heterogeneous media large enough for boundaries to be not important. This limitation along with the Bloch-Floquet wave approach that only works for periodic cases poses challenges when it comes to predicting the effective constitutive relations of finite-size non-periodic heterogeneous media for real-world applications.

This paper presents two novel contributions to the linear elastodynamic homogenization of general heterogeneous media, including finite-size and non-periodic ones, for their effective constitutive relations. First, we introduce a new homogenization formulation based on the Hashin-Shtrikman principle, which takes into account boundary terms that have not been adequately considered previously. This formulation allows the use of the infinite-body Green’s function of the reference homogeneous medium. The effective constitutive relations obtained from this formulation display the non-classical features of the homogenized medium. Second, we demonstrate the universal applicability of the proposed formulation for both random and deterministic heterogeneous media. With these two contributions, this work provides a better understanding of the homogenization theory and enables the design of next-generation composites and metamaterials in which securing their effective characteristics is critical for dynamic wave manipulation.

Consider a boundary-value problem of a heterogeneous medium $\Omega$ which has random or deterministic microstructures but its boundary conditions and body force are taken to be sure, in arbitrary space dimension. Here, we proceed with piezoelectric media, which possess both Willis and electro-momentum coupling, as the homogenization formulation of standard (non-piezoelectric) media only with Willis coupling can be easily reduced by taking zero piezoelectricity as shown later in this paper. Thus, our goal is the development of a general formulation to construct effective dynamic constitutive relations of the heterogeneous medium, i.e.,

with the transpose operator $(\cdot)^{T}$, in which the quantities with angle brackets mean the homogenized macroscopic ones where $\bm{\varepsilon}=\frac{1}{2}\left[\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{T}\right]$ with spatial gradient operator $\bm{\nabla}$ and displacement $\bm{u}$, and $\bm{v}=\dot{\bm{u}}$ with over-bar being time derivative, and $\bm{\nabla}\phi=-\bm{E}$ for electric potential $\phi$. The over-tilde is used to denote the effective material properties for stiffness $\bm{C}$, dielectricity $\bm{A}$, and piezoelectricity $\bm{B}$, and mass density $\rho$. We note here the effective mass density has a tensorial form. $\widetilde{\bm{S}}$ and $\widetilde{\bm{W}}$ are the Willis and electro-momentum coupling terms, respectively. As shown later in this paper, the effective material properties are non-local operators in general such that, for instance, $\left(\widetilde{\bm{C}}:\langle\bm{\varepsilon}\rangle\right)(\bm{x})=\int{\widetilde{\bm{C}}(\bm{x},\check{\bm{x}}):\langle\bm{\varepsilon}(\check{\bm{x}})\rangle}~{}d\Omega(\check{\bm{x}})$ with respect to the spatial coordinates $\bm{x}$ and $\check{\bm{x}}$.

The local mechanical and electric fields of the heterogeneous medium must satisfy the governing equations on $\Omega$

where $\bm{f}$ is body force, and $q$ is free charge density. The constitutive relations between the local kinetic and kinematic fields are written as

in which the material properties depend on the microstructure materials of the heterogeneous medium at the local position. Here, we assume the standard tensor symmetries on the properties $\bm{C}$, $\bm{A}$, and $\bm{B}$. Boundary conditions along the boundary $\partial\Omega$ are given as

where $\bm{t}$ and $\bm{u}_{bc}$ are prescribed traction and displacement, respectively, with the outward unit normal vector $\bm{n}$, and $c$ and $\phi_{bc}$ are prescribed surface charge density and electric potential, respectively. Note that $\partial\Omega_{t}\cap\partial\Omega_{u}=\partial\Omega_{c}\cap\partial\Omega_{\phi}=\varnothing$.

Now following the Hashin-Shtrikman principle [42] schematically shown in Fig. 1, let’s introduce a homogeneous reference medium with constant properties $\bm{C}_{0}$, $\bm{A}_{0}$, $\bm{B}_{0}$, and $\rho_{0}$, which may or may not be one of the microstructure materials of the heterogeneous medium. Then, field polarizations are defined as

where $\bm{\tau}$, $\bm{\chi}$, $\bm{\pi}$ are the polarization of stress, electric displacement, and linear momentum, respectively. Hereafter, we will write matrix equations for the shake of simplicity via symbolic expression based on the notations in Fig. 1 — for instance of (3) and (5), $h=Lb$ and $\zeta=\Delta Lb$ with $\Delta L=L-L_{0}$, respectively. Then, the boundary-value problem of the heterogeneous medium can be equivalently represented, as shown in Fig. 1(b), via a boundary-value problem of the homogeneous reference medium with the constant properties $L_{0}$ as

where $\bm{\sigma}_{r}=\bm{\sigma}-\bm{\tau}$, $\bm{\mu}_{r}=\bm{\mu}-\bm{\pi}$, and $\bm{D}_{r}=\bm{D}-\bm{\chi}$ such that $h_{r}=L_{0}b$ via the symbolic expression, with boundary conditions

As we consider a linear elastodynamic system, this equivalent boundary-value problem can be decomposed via the superposition principle as shown in Fig. 1(c) and (d), in which Fig. 1(c) is the reference homogeneous medium subjected to the given Neumann boundary conditions of the heterogeneous medium to be homogenized. This leads to a resultant boundary-value problem building on the system of Fig. 1(d) that is independent of both the source terms $\bm{f}$ and $q$ and the Neumann boundary conditions $\bm{t}$ and $c$ as

such that $h^{\prime}=L_{0}b^{\prime}$ as the symbolic expression, with boundary condition

where $\bm{u}_{bc_{0}}$ and $\phi_{bc_{0}}$ are arbitrary Dirichlet boundary conditions on the system of Fig. 1(c). Then, the problem can be solved via the boundary-element method using the infinite-body Green’s function of the reference homogeneous medium. By employing the Laplace transform to replace the time derivative by the Laplace variable $s=-i\omega$ with angular frequency $\omega$, the solution is obtained as

as the symbolic expression with

in which the time-reduced Green’s function matrix $G(\bm{x},\check{\bm{x}})$ is defined to solve [43, 44]

where $\delta_{x\check{x}}$ is the Dirac-delta function $\delta(\bm{x}-\check{\bm{x}})$, and $\bm{\delta}_{x\check{x}}=\delta(\bm{x}-\check{\bm{x}})\bm{I}$ with the identity tensor $\bm{I}$. We note that the over-check of $\check{B}$ means the gradients are associated with $\check{\bm{x}}$ which comes from the non-local nature of Green’s function.

The solution (10) requires only the infinite-body Green’s function of a homogeneous medium, and its integration volume extends over the region occupied by the target heterogeneous medium to be homogenized. While different homogenization theories were derived based on the reference homogeneous medium in the previous literature [22, 24, 26, 28, 29, 41], the boundary terms in the Hashin-Shtrikman principle-based approach were not considered before. This results in constraints on Green’s function of the reference medium in their homogenization. Therefore, the previous formulations only work for heterogeneous media of infinite (at least quasi-infinite) size or simple geometries to which Green’s function can be determined.

The solution (10) yields the kinematic fields as

with $\Gamma(\bm{x},\check{\bm{x}})=B\Big{(}\check{B}G(\bm{x},\check{\bm{x}})\Big{)}^{T}$ where $B$ is replacing the gradients in $\check{B}$ to be associated with $\bm{x}$. We note that $\Gamma(\bm{x},\check{\bm{x}})$ is dominated by its singularity near $\bm{x}=\check{\bm{x}}$ such that it can be written as $\Gamma(\bm{x},\check{\bm{x}})=\mathcal{D}\delta_{x\check{x}}+\Gamma^{H}(\bm{x},\check{\bm{x}})$. $\mathcal{D}$ is constant which represents the contribution of an infinitesimal region around the singularity position while $\Gamma^{H}(\bm{x},\check{\bm{x}})$ is of the contribution outside of the singular region. This integral equation can be written via a compact integral operator form as $b=\langle b\rangle+\Gamma\langle\zeta\rangle-\Gamma\zeta$. From this, one can find a relation between $\langle\zeta\rangle$ and $\langle b\rangle$ that can be expressed in a general form of

in which $\mathcal{I}$ is the identity operator, and $(\cdot)^{+}$ denotes the inverse integral operator, e.g., $\Gamma^{+}\Gamma=\mathcal{I}$. Therefore, by the definition of field polarizations (i.e., $h=L_{0}b+\zeta$), the effective material property matrix $L_{eff}$, for the effective constitutive relations $\langle h\rangle=L_{eff}\langle b\rangle$ as the compact integral operator form of (1), can be finally obtained as

which well exhibits the non-local feature of the effective constitutive relations of the homogenized medium. This equation can be solved by various methods such as the iteration methods, e.g., the Neumann series and the resolvent method [45], and the discretization by quadrature rules [46] that can provide a direct connection to the finite element method.

This derived formulation is general but when $\bm{C}^{T}=\bm{C}$ (i.e. $C_{ijkl}=C_{lkji}$), $L_{eff}$ is self-adjoint based on the self-adjointness of Green’s function, which delivers $\widetilde{\bm{S}}_{2}=\widetilde{\bm{S}}_{1}^{\dagger}$ and $\widetilde{\bm{W}}_{2}=\widetilde{\bm{W}}_{1}^{\dagger}$ with the formal adjoint operator $(\cdot)^{\dagger}$ with respect to the spatial variable, e.g., $\widetilde{\bm{S}}_{1}^{\dagger}(\bm{x},\check{\bm{x}})=\widetilde{\bm{S}}_{1}(\check{\bm{x}},\bm{x})$ and $\widetilde{\bm{W}}_{1}^{\dagger}(\bm{x},\check{\bm{x}})=\widetilde{\bm{W}}_{1}(\check{\bm{x}},\bm{x})$. As (15) is derived from (10), it is unique regardless of both source terms and boundary conditions on the heterogeneous medium to be homogenized. This uniqueness is a reminiscence of the source-driven homogenization [47], and this implies that the non-uniqueness issue in the absence of both source terms and boundary conditions can be mitigated via applying eigenstrains [27]. Furthermore, since we haven’t made any assumptions on the homogenization strategy to map microscopic quantities to homogenized macroscopic ones, the proposed homogenization formulation is universal. We note that this formulation approach can also be employed to obtain effective constitutive relations in different physics, e.g., electromagnetism.

To validate the proposed homogenization formulation and demonstrate its universal applicability, we compare the effective constitution obtained from our formulation to 1-D analytical solutions of problems for random and deterministic heterogeneous media.

First, consider a homogenization problem based on the ensemble average technique for a finite-size random heterogeneous medium with a periodically repeating unit cell, as shown in Fig. 2(a). This unit cell consists of PMMA–BaTiO₃–PZT4–PVDF layers with 0.5 mm, 0.5 mm, 1 mm, and 1.5 mm of each layer thickness, a total 3.5 mm of the unit cell thickness $l_{uc}$. Relying on the periodic microstructure geometry, the position of period can be considered as a uniformly distributed random variable that has $l_{uc}^{-1}$ of uniform probability density over the unit cell domain $\Omega_{uc}$. Therefore, this feature reduces the ensemble average of any $l_{uc}$-periodic function $\Upsilon$ in the homogenization to a spatial average as $\langle\Upsilon(x)\rangle=l_{uc}^{-1}\int_{\Omega_{uc}}{\Upsilon_{y_{0}}(x-y)}~{}dy$ in which the $l_{uc}$-periodic function in realization $y\in\Omega_{uc}$ is represented using $\Upsilon_{y_{0}}(x-y)$ with the subscript $y_{0}$ denoting the realization at $y=0$. Assuming the absence of free charge and external electromagnetic fields, the electric governing equation $\nabla D=0$ in (2) gives zero electric displacement across the medium domain, which leads to $\nabla\phi=(BA^{-1})\varepsilon$. Here, the bold is omitted in a 1-D manner. In this case, the boundary-value problem of the piezoelectric medium becomes a pure mechanical problem in which the local and effective constitutive relations are modified as

where the subscript $a$ denotes the apparent quantities as $C_{a}=C+B^{2}A^{-1}$ for the microscopic property, and $\widetilde{C}_{a}=\widetilde{C}+\widetilde{B}^{2}\widetilde{A}^{-1}$, $\widetilde{S}_{1_{a}}=\widetilde{S}_{1}+\widetilde{B}\widetilde{W}_{1}\widetilde{A}^{-1}$, $\widetilde{S}_{2_{a}}=\widetilde{S}_{2}+\widetilde{B}\widetilde{W}_{2}\widetilde{A}^{-1}$, and $\widetilde{\rho}_{a}=\widetilde{\rho}+\widetilde{W}_{1}\widetilde{W}_{2}\widetilde{A}^{-1}$ for the effective properties. By the symbolic expression, we write (16) as $h_{a}=L_{a}b_{a}$ and $\langle h_{a}\rangle=L_{eff_{a}}\langle b_{a}\rangle$ respectively for the local and effective relation, with $h_{a}=(\sigma~{}\mu)^{T}$ and $b_{a}=(\varepsilon~{}v)^{T}$. Then, (13) is reduced, via the compact integral operator form, as

where $\zeta_{a}=(\tau~{}\pi)^{T}$, and $\Gamma_{a}(x,\check{x})=B_{a}\Big{(}\check{B}_{a}^{T}G_{a}(x,\check{x})\Big{)}$ with $\check{B}_{a}=(\check{\nabla}~{}s)^{T}$ and $B_{a}$ replacing the gradient to be associated with $x$, which can be decomposed as $\Gamma_{a}(x,\check{x})=\mathcal{D}_{a}\delta_{x\check{x}}+\Gamma_{a}^{H}(x,\check{x})$, where $\mathcal{D}_{a}$ is constant, based on its singularity near $x=\check{x}$. The apparent infinite-body time-reduced Green’s function $G_{a}(x,\check{x})$ is defined to solve

where $C_{0_{a}}$ is the apparent stiffness of the reference homogeneous medium, i.e., [48]

with $k_{0_{a}}=\omega\sqrt{\rho_{0}/C_{0_{a}}}$. Then, this leads to

where $\Delta L_{a}=L_{a}-L_{0_{a}}$ with $L_{0_{a}}$ being the apparent material property matrix of the reference homogeneous medium. We note that taking the piezoelectricity $B\rightarrow 0$ exactly recovers the effective constitutive relations of standard heterogeneous media with only Willis coupling. One can easily expand it to the multidimensional case. Fig. 2(b-e) show macroscopic stress and linear momentum fields across the medium domain, under prescribing displacements at its left and right ends respectively by zero and 0.2% of the medium size, for the different number of the unit cell at time frequency $f_{q}(=\frac{\omega}{2\pi})=0.24$ MHz. The fields obtained using (20) are compared with the ensemble average of analytical solutions of each ensemble. These results show that our proposed homogenization formulation provides an excellent estimate of the effective material properties for finite-size random heterogeneous media and displays well the necessity of the non-classical coupling terms in the effective constitutive relations of the homogenized medium, i.e., Willis and electro-momentum coupling. Fig. 2(f) and (g) visualize the apparent coupling coefficients and confirm their non-local feature in space and the adjoint relation between $\widetilde{S}_{1_{a}}$ and $\widetilde{S}_{2_{a}}$.

Second, we extract the non-classical coupling terms of the homogenized medium for a finite-size non-periodic deterministic heterogeneous medium and compare them to benchmark data of [35]. Consider a piezoelectric heterogeneous element made of 2.62 mm of PZT4, 0.19 mm of BaTiO₃, and 0.19 mm of PVDF. We here also assume the absence of free charge and external electric fields. The microstructures of this element possess wavelength much larger than those characteristic lengths, i.e., $kl\ll 1$ with wavenumber $k$ and length $l$ of the layer. When a heterogeneous medium is at the subwavelength scale, or when a single isolated element is used rather than a bulk of mutually interacting elements, the effective constitutive relations become local in space, termed Milton-Briane-Willis equation [49, 50]. It has been studied that the local coupling terms originate from broken inversion asymmetry of the microstructures [51, 47]. Then, the apparent kinematic field vector $b_{a}$ in the modified constitutive relations (16) can be related to the state vector $\varsigma=(u~{}\sigma)^{T}$ via the material properties as $b_{a}=M\varsigma$ and $\langle b_{a}\rangle=M_{eff}\langle\varsigma\rangle$ via the symbolic expression, with

where the derivative of the state vector can also be calculated using the material properties and itself as $\nabla\varsigma=KM\varsigma$ and $\nabla\langle\varsigma\rangle=K_{eff}M_{eff}\langle\varsigma\rangle$, with

For the $i$-th layer, the local state vectors at its two ends $x_{L}^{(i)}$ and $x_{R}^{(i)}$ ($x_{L}^{(i)}<x_{R}^{(i)}$) are related via $\varsigma(x_{R}^{(i)})=T_{l}^{(i)}\varsigma(x_{L}^{(i)})$ with the layer’s transfer matrix $T_{l}^{(i)}$ as

with $k_{a}^{(i)}=\omega\sqrt{\rho^{(i)}/C_{a}^{(i)}}$. By the subwavelength characteristic, $\Gamma_{a}^{H}(x,\check{x})$ can be expressed by expanding it around the element ends upto the linear term such that

with the element total length $l_{t}$. Then, based on the continuity of the state vector at the layer intersection and $\langle\zeta_{a}\rangle=(L_{eff_{a}}-L_{0_{a}})\langle b_{a}\rangle$, integrating (17) over $\Omega(x)$ gives an equation to find $L_{eff_{a}}$ of the deterministic heterogeneous element at the subwavelength scale as

in terms of the arbitrary $\varsigma(x_{L}^{(1)})$ with $n_{l}$ being the number of layer, in which the fact that the state vector is sure along the boundary is used, i.e., $\langle\varsigma(x_{L}^{(1)})\rangle=\varsigma(x_{L}^{(1)})$, with

where $T^{(i)}=T_{l}^{(i)}T_{l}^{(i-1)}\ldots T_{l}^{(1)}$ with $T_{l}^{(0)}=I$, $\Gamma_{a+}^{H}=\Gamma_{a}^{H}(l_{t})$ and $\Gamma_{a-}^{H}=\Gamma_{a}^{H}(-l_{t})+\nabla\Gamma_{a}^{H}(-l_{t})l_{t}$, and $I$ is the identity matrix. This equation can be easily solved by using numerical methods, such as the multivariate Newton–Raphson method, with respect to the components in $L_{eff_{a}}$, i.e., $\widetilde{C}_{a}$, $\widetilde{S}_{1_{a}}$, $\widetilde{S}_{2_{a}}$, and $\widetilde{\rho}_{a}$. Fig. 3 shows the pure imaginary apparent coupling coefficients of the target heterogeneous element and its inversion, i.e., when the first and third layers are interchanged. The coupling coefficients obtained from our formulation have $\widetilde{S}_{1_{a}}=\widetilde{S}_{2_{a}}$ and are opposite in sign for the regular and inverted elements, i.e., a direction-dependent response, which is coincident with the observation in [35]. The obtained coupling coefficient values are also well agreed with the benchmark data. These results demonstrate the universal applicability of our proposed formulation not only for ensemble average-based homogenization on random heterogeneous media but also for deterministic cases.

In conclusion, we have derived a formulation to find the dynamic effective constitutive relations for general heterogeneous media, including finite-size and non-periodic ones, by carefully taking into account boundary terms in the Hashin-Shtrikman principle. Our formulation relies on the infinite-body Green’s function of a reference homogeneous medium, making it free from the difficulty of determining Green’s function even for the homogenization of finite-size media. We note that as our theory generalizes the homogenization theory on two-phase disordered heterogeneous media [52, 53, 54], our formulation satisfies rapid convergence at least up to two-phase disordered heterogeneous media, with the decomposition of $\Gamma(\bm{x},\check{\bm{x}})$. It is envisioned that this study will enable more accurate predictions of the effective material properties of heterogeneous media and expand to more real-world applications.