1. Introduction
Many research methods were proposed to deal with the effective properties of composites, such as self-consistent model [1-4], variational approach [5, 6], statistical method and exact solutions [7, 8]. In these studies of composites, there are currently two types of boundary conditions: periodic boundary conditions and open boundary conditions [5, 8-13]. Periodic boundary conditions proposed by Rayleigh [14] for the periodic model of composites were developed by Nemat-Nasser et al [15] for accommodating the interactions between inclusion particles without any limitations on inclusion shapes. However, the periodic model and periodic boundary conditions still differ from real composites. Actually, the real composite medium is that inclusion particles are confined to a limited region around an observation point when the effective response of the composite medium is examined by applying an external field. At a distance from the observation point, the disturbance field induced by inclusion particles tends to zero. Therefore, as long as the field quantity is equal to the external field at the boundary or the disturbance field is taken as zero at infinity [16], it is referred to as open boundary conditions in physics. Therefore, open boundary conditions are more suitable for real composites.
For the open boundary composites of regular geometric inclusions embedded in an infinite host, some analytical solutions have been derived for the electric, elastic, piezoelectric and thermal composites, such as cylindrical, ellipsoidal and spherical composites [17-24]. The open boundary problem is more primitive and very difficult to analytically solve for the composites having complex geometric inclusions due to the complex interactions between inclusion shapes and material properties. To investigate the open boundary problem of composites, it is possible to develop Eshelby's transformation field method (TFM) to estimate the effective responses of composites having arbitrary inclusion shapes [13, 16, 25-27].
Eshelby's TFM had been proposed to calculate the elastic tensors of whole ellipsoidal composites [16, 25], where the elastic tensor can be linearly expressed by the uniform transformation strain (i.e. stress-free strain) introduced by an ellipsoid inclusion through the Eshelby's uniform tensor, which contains the geometric parameters and physical properties of ellipsoidal inclusions. Generally, Eshelby also conjectured that above property (i.e. Eshelby's property) does not hold for non-ellipsoidal inclusions. However, considering the non-uniform transformation strain tensor, Nemat-Nasser et al [15, 26] had developed TFM to estimate the overall stress and the effective elastic moduli of periodic composites having arbitrary inclusion geometry. Subsequently, TFM was extended by Gu et al [28-32] to predict the effective electric conductivity, viscosity, piezoelectric properties and dispersion relation of photonic band structure of periodic composites. For the graded periodic composites, TFM was used to calculate the effective dielectric constants of periodic composites with arbitrary gradient and geometric inclusions [33, 34]. Recently, TFM was successfully developed to estimate the effective responses of unidirectional periodic composites having open boundary condition only along one axis direction [35]. However, the open boundary problem is not completely solved by TFM for the composites having arbitrary structure inclusions embedded in an infinite host. In actual physics or engineering, the composites of one-dimensional open boundary condition are not common. If TFM can be generally used to deal with open boundary conditions in all directions, it will provide a firm theoretical basis for the study of the physical properties of composites. Therefore, along this line, we will investigate TFM to solve the fully open boundary problem of composites having arbitrary geometric inclusions.
In this paper, to break the restriction of Eshelby's uniform transformation field applied only to an ellipsoid composite, we develop the non-uniform transformation field to solve the open boundary problem of composites having complex geometric inclusions. The developed TFM copes with not only the composites of satisfying Eshelby's property (including ellipsoidal conform maps) but also non-ellipsoid composites, where the number of inclusions is more than one and the shape of inclusions is arbitrary. As an example, the two-dimensional dielectric composites having open boundary conditions in all directions are investigated by TFM. The transformation electrical fields in composite regions are introduced to set up a unified constitutive relation. To match open boundary conditions, the perturbation potential and the transformation electrical fields are expressed in terms of weighted orthogonal Hermite polynomials because the polynomial converges to a constant at infinite point and can close the equations of perturbation potentials. Furthermore, the transformation electric fields are obtained by solving a set of algebraic equations with respect to the unknown coefficients of transformation fields. The formulas of effective dielectric responses are derived by the transformation fields. Finally, the validity is verified by comparing the effective responses estimated by TFM with the exact solutions of dilute limit.
2. Formulation
2.1. Transformation field and its governing equation
Consider two-phase isotropic dielectric composites, where the inclusions of arbitrary shapes are embedded in an infinite matrix. For two-dimensional composites, we assume that the isotropic constitutive relations of both host and inclusion are linear in Cartesian coordinates $(x,y)$, where the origin of coordinates is at the center of an inclusion.
$\begin{eqnarray}{D}_{k}^{\alpha }={\varepsilon }^{\alpha }{E}_{k}^{\alpha }\left(x,y\right),\end{eqnarray}$
where $k=x,y$. $\alpha =h,i$ denote the quantities of the host and inclusion regions, respectively. ${\boldsymbol{D}}$, ${\boldsymbol{E}}$ and $\varepsilon $ are the electrical displacement, electrical field and dielectric constant, respectively. The electrostatic governing equations are ${\rm{\nabla }}\cdot {{\boldsymbol{D}}}^{\alpha }=0$ and ${\rm{\nabla }}\times {{\boldsymbol{E}}}^{\alpha }=0$. If an external uniform electric field ${E}_{k}^{0}$ is applied to the composites along the $k$-direction, the perturbation electric field ${E}_{k}^{p}$, induced by the inclusions, will occur in the whole composite region. In the inclusion region ${{\rm{\Omega }}}_{i}$ and the host region ${{\rm{\Omega }}}_{h}$, the constitutive relations are given by ${D}_{k}^{i}={\varepsilon }^{i}\left({E}_{k}^{0}+{E}_{k}^{p}\right)$ and ${D}_{k}^{h}={\varepsilon }^{h}\left({E}_{k}^{0}+{E}_{k}^{p}\right)$, respectively. Meanwhile, the governing equation requires that ${{\rm{\nabla }}}_{k}\left\{{\varepsilon }^{h}\left[{E}_{k}^{0}+{E}_{k}^{p}\left(x,y\right)\right]\right\}=0$ in host region ${{\rm{\Omega }}}_{h}$ and ${{\rm{\nabla }}}_{k}\left\{{\varepsilon }^{i}\left[{E}_{k}^{0}+{E}_{k}^{p}\left(x,y\right)\right]\right\}=0$ in the inclusion region ${{\rm{\Omega }}}_{i}$. The continuity of normal electrical displacements and electric field on the interface $\partial {{\rm{\Omega }}}_{i}$ of inclusion region ${{\rm{\Omega }}}_{i}$ is satisfied if $\left[{\varepsilon }^{h}{\left({{\boldsymbol{E}}}^{0}+{{\boldsymbol{E}}}^{p}\right)}^{h}-{\varepsilon }^{i}{\left({{\boldsymbol{E}}}^{0}+{{\boldsymbol{E}}}^{p}\right)}^{i}\right]\cdot {\boldsymbol{n}}=0$ and $\left[{\left({{\boldsymbol{E}}}^{0}+{{\boldsymbol{E}}}^{p}\right)}^{h}-{\left({{\boldsymbol{E}}}^{0}+{{\boldsymbol{E}}}^{p}\right)}^{i}\right]\times {\boldsymbol{n}}=0$ on $\partial {{\rm{\Omega }}}_{i}$, where ${\boldsymbol{n}}$ is the normal unit vector on boundary $\partial {{\rm{\Omega }}}_{i}$ of inclusion region ${{\rm{\Omega }}}_{i}$, and the superscripts $h$ and $i$ denote the quantities on the boundary of host and inclusion, respectively.Instead of solving above governing equations directly, a transformation electrical field ${E}_{k}^{* }$ is introduced to the inclusion according to the concept of Eshelby's transformation field, where the introduced transformation field can satisfy the continuous boundary conditions of electric displacement and electric field on the interface of inclusion and host [16, 26]. In this case, the transformation electrical fields obey the following equations in the host and inclusion regions, respectively [28, 33].2 ) and (3 ), a unified constitutive relation is obtained in the whole composite region ${\rm{\Omega }}={{\rm{\Omega }}}_{i}+{{\rm{\Omega }}}_{h}$.5 ) gives the relationship between the perturbation electrical field and the transformation electrical field. From the physics perspective, introducing a transformation field into the composite converts the problem of fitting complex boundary conditions on the interfaces between inclusion and host into a Poisson equation, and the source term of the Poisson equation is related to the transformation field. Therefore, to solve the perturbation electric field ${E}_{k}^{p}$ of inclusion region, the source term ${{\rm{\nabla }}}_{k}\left[{\varepsilon }^{h}{E}_{k}^{\ast }\left(x,y\right)\right]$ of equation (5 ) should be considered. Next, the perturbation potential will be derived from equation (5 ) and expressed by the transformation electrical field.
$\begin{eqnarray}{E}_{k}^{\ast }\left(x,y\right)=0\,{\rm{in}}\,{{\rm{\Omega }}}_{h},\end{eqnarray}$
$\begin{eqnarray}{\varepsilon }^{i}\left[{E}_{k}^{0}+{E}_{k}^{p}\left(x,y\right)\right]\,=\,{\varepsilon }^{h}\left[{E}_{k}^{0}+{E}_{k}^{p}\left(x,y\right)-{E}_{k}^{\ast }\left(x,y\right)\right]\,{\rm{in}}\,{{\rm{\Omega }}}_{i},\end{eqnarray}$
where $k=x,y$. Based on equations ( $\begin{eqnarray}{D}_{k}={\varepsilon }^{h}\left[{E}_{k}^{0}+{E}_{k}^{p}\left(x,y\right)-{E}_{k}^{\ast }\left(x,y\right)\right]\,{\rm{in}}\,{\rm{\Omega }}.\end{eqnarray}$
Then, the governing equation of the perturbation electrical field ${E}_{k}^{p}$ and the transformation electrical field ${E}_{k}^{* }\left(x,y\right)$ is rewritten as follows, $\begin{eqnarray}{{\rm{\nabla }}}_{k}\left\{{\varepsilon }^{h}\left[{E}_{k}^{0}+{E}_{k}^{p}\left(x,y\right)-{E}_{k}^{\ast }\left(x,y\right)\right]\right\}=0.\end{eqnarray}$
Thus, equation (2.2. Solutions of perturbation potential
For the open boundary composite, in order to match the open boundary conditions of the perturbation potential ${\rm{\Phi }}\left(x,y\right)$, i.e. the perturbation electrical field vector ${{\boldsymbol{E}}}^{p}\left(x,y\right)=-{\rm{\nabla }}{\rm{\Phi }}$ and ${{\boldsymbol{E}}}^{p}{\left.\left(x,y\right)\right|}_{\left(x,y\right)\to \infty }=0$, the perturbation potential ${\rm{\Phi }}\left(x,y\right)$ and the transformation field ${E}_{k}^{\ast }\left(x,y\right)$, (i.e. ${\left.{E}_{k}^{\ast }\left(x,y\right)\right|}_{\left(x,y\right)\to \infty }=0$ on open boundary) are expanded in terms of weighted Hermite polynomials,9 ) by substituting equations (6 )-(8 ) into equation (5 ),9 ), the following derivations of both ${u}_{m}\left(x\right)$ and ${v}_{n}\left(y\right)$ with respect to variables $x$ and $y$ are used, respectively.9 ) is an infinite series equation of both orthogonal functions ${u}_{m}\left(x\right)$ and ${v}_{n}\left(y\right)$. Considering the orthogonality of functions ${u}_{m}\left(x\right)$ and ${v}_{n}\left(y\right)$, we obtain a set of equations of the coefficients ${{\rm{\Phi }}}_{m,n}$ by truncating $m\,=\,0,1,2,\mathrm{.}.,\,M$ and $n\,=\,0,1,2,\mathrm{.}.,\,N$ in equation (9 ),
$\begin{eqnarray}{\rm{\Phi }}\left(x,y\right)=\displaystyle \sum _{m=0}^{\infty }\displaystyle \sum _{n=0}^{\infty }{{\rm{\Phi }}}_{m,n}{u}_{m}\left(x\right){v}_{n}\left(y\right),\end{eqnarray}$
$\begin{eqnarray}{E}_{x}^{\ast }\left(x,y\right)=\displaystyle \sum _{m=0}^{\infty }\displaystyle \sum _{n=0}^{\infty }{\alpha }_{m,n}{u}_{m}\left(x\right){v}_{n}\left(y\right),\end{eqnarray}$
$\begin{eqnarray}{E}_{y}^{\ast }\left(x,y\right)=\displaystyle \sum _{m=0}^{\infty }\displaystyle \sum _{n=0}^{\infty }{\beta }_{m,n}{u}_{m}\left(x\right){v}_{n}\left(y\right),\end{eqnarray}$
where ${{\rm{\Phi }}}_{m,n}$, ${\alpha }_{m,n}$ and ${\beta }_{m,n}$ are unknown coefficients. The orthogonal function ${u}_{m}\left(x\right)$ is ${u}_{m}\left(x\right)={N}_{m}{H}_{m}\left(x\right)\exp \left(-{x}^{2}/2\right)$, where ${H}_{m}\left(x\right)$ is Hermite polynomial and the coefficient ${N}_{m}={\left({2}^{m}m!\sqrt{\pi }\right)}^{-1/2}$. ${v}_{n}\left(y\right)={N}_{n}{H}_{n}\left(y\right)\exp \left(-{y}^{2}/2\right)$. Since the orthogonal function ${u}_{m}\left(x\right)$(or ${v}_{n}\left(y\right)$) tends to a constant at an infinite point (i.e. satisfies the open boundary conditions of the perturbation potential and the transformation field ), we get equation ( $\begin{eqnarray}\begin{array}{c}\displaystyle \sum _{m,n=0}^{\infty }{{\rm{\Phi }}}_{m,n}\left[{\tilde{M}}_{m}^{+}{u}_{m+2}\left(x\right)+{\tilde{M}}_{m}^{0}{u}_{m}\left(x\right)+{\tilde{M}}_{m}^{-}{u}_{m-2}\left(x\right)\right]{v}_{n}\left(y\right)\\ \,+\displaystyle \sum _{m,n=0}^{\infty }{{\rm{\Phi }}}_{m,n}\left[{\tilde{M}}_{n}{v}_{n+2}\left(y\right)+{\tilde{M}}_{n}^{0}{v}_{n}\left(y\right)+{\tilde{M}}_{n}^{-}{v}_{n-2}\left(y\right)\right]{u}_{m}\left(x\right),\\ \,=-\displaystyle \sum _{m,n=0}^{\infty }{\alpha }_{m,n}\left[{\tilde{N}}_{m}^{+}{u}_{m+1}\left(x\right)+{\tilde{N}}_{m}^{-}{u}_{m-1}\left(x\right)\right]{v}_{n}\left(y\right)\\ \,-\displaystyle \sum _{m,n=0}^{\infty }{\beta }_{m,n}\left[{\tilde{N}}_{n}^{+}{v}_{n+1}\left(x\right)+{\tilde{N}}_{n}^{-}{v}_{n-1}\left(x\right)\right]{u}_{m}\left(x\right),\end{array}\end{eqnarray}$
where ${\tilde{N}}_{m}^{+}\,=\,-\displaystyle \frac{1}{2}{N}_{m}/{N}_{m+1}$, ${\tilde{N}}_{m}^{-}\,=\,m{N}_{m}/{N}_{m+1}$, ${\tilde{M}}_{m}^{+}\,=\,{\tilde{N}}_{m}^{+}{\tilde{N}}_{m+1}^{+}$, ${\tilde{M}}_{m}^{0}\,=\,{\tilde{N}}_{m}^{+}{\tilde{N}}_{m+1}^{-}+{\tilde{N}}_{m}^{-}{\tilde{N}}_{m-1}^{+}$, ${\tilde{M}}_{m}^{-}\,=\,{\tilde{N}}_{m}^{-}{\tilde{N}}_{m-1}^{-}$. To derive equation ( $\begin{eqnarray}\displaystyle \frac{{\rm{d}}{u}_{m}\left(x\right)}{{\rm{d}}x}\,=\,{\tilde{N}}_{m}^{+}{u}_{m+1}\left(x\right)+{\tilde{N}}_{m}^{-}{u}_{m-1}\left(x\right),\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{{\rm{d}}{v}_{n}\left(y\right)}{{\rm{d}}y}\,=\,{\tilde{N}}_{n}^{+}{v}_{n+1}\left(y\right)+{\tilde{N}}_{n}^{-}{v}_{n-1}\left(y\right),\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}}^{2}{u}_{m}\left(x\right)}{{\rm{d}}{x}^{2}}\,=\,{\tilde{M}}_{m}^{+}{u}_{m+2}\left(x\right)+{\tilde{M}}_{m}^{0}{u}_{m}\left(x\right)+{\tilde{M}}_{m}^{-}{u}_{m-2}\left(x\right),\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}}^{2}{v}_{n}\left(y\right)}{{\rm{d}}{y}^{2}}\,=\,{\tilde{M}}_{n}^{+}{v}_{n+2}\left(y\right)+{\tilde{M}}_{n}^{0}{v}_{n}\left(y\right)+{\tilde{M}}_{n}^{-}{v}_{n-2}\left(y\right).\end{eqnarray}$
Clearly, equation ( $\begin{array}{c} 2 \Phi_{0,0} \tilde{M}_{0}^{0}+\Phi_{2,0} \tilde{M}_{2}^{-}+\Phi_{0,2} \tilde{M}_{2}^{-}=-\left(\alpha_{1,0} \tilde{N}_{1}^{-}+\beta_{0,1} \tilde{N}_{1}^{-}\right), \\ \Phi_{1,0}\left(\tilde{M}_{1}^{0}+\tilde{M}_{0}^{0}\right)+\Phi_{3,0} \tilde{M}_{3}^{-}+\Phi_{1,2} \tilde{M}_{2}^{-} \\ =-\left(\alpha_{0,0} \tilde{N}_{0}^{+}+\alpha_{2,0} \tilde{N}_{2}^{-}+\beta_{1,1} \tilde{N}_{1}^{-}\right), \\ \Phi_{0,1}\left(\tilde{M}_{0}^{0}+\tilde{M}_{1}^{0}\right)+\Phi_{2,1} \tilde{M}_{2}^{-}+\Phi_{0,3} \tilde{M}_{3}^{-} \\ =-\left(\alpha_{1,1} \tilde{N}_{1}^{-}+\beta_{0,0} \tilde{N}_{0}^{+}+\beta_{0,2} \tilde{N}_{2}^{-}\right), \\ 2 \Phi_{1,1} \tilde{M}_{1}^{0}+\Phi_{3,1} \tilde{M}_{3}^{-}+\Phi_{1,3} \tilde{M}_{3}^{-} \\ =-\left(\alpha_{0,1} \tilde{N}_{0}^{+}+\alpha_{2,1} \tilde{N}_{2}^{-}+\beta_{1,0} \tilde{N}_{0}^{+}+\beta_{1,2} \tilde{N}_{2}^{-}\right), \\ \Phi_{0,0} \tilde{M}_{0}^{+}+\Phi_{2,0}\left(\tilde{M}_{2}^{0}+\tilde{M}_{0}^{0}\right)+\Phi_{4,0} \tilde{M}_{4}^{-}+\Phi_{2,2} \tilde{M}_{2}^{-} \\ =-\left(\alpha_{1,0} \tilde{N}_{1}^{+}+\alpha_{3,0} \tilde{N}_{3}^{-}+\beta_{2,1} \tilde{N}_{1}^{-}\right), \\ \Phi_{0,1} \tilde{M}_{0}^{+}+\Phi_{2,1}\left(\tilde{M}_{2}^{0}+\tilde{M}_{1}^{0}\right)+\Phi_{4,1} \tilde{M}_{4}^{-}+\Phi_{2,3} \tilde{M}_{3}^{-} \\ =-\left(\alpha_{1,1} \tilde{N}_{1}^{+}+\alpha_{3,1} \tilde{N}_{3}^{-}+\beta_{2,0} \tilde{N}_{0}^{+}+\beta_{2,2} \tilde{N}_{2}^{-}\right), \\ \Phi_{0,0} \tilde{M}_{0}^{+}+\Phi_{0,2}\left(\tilde{M}_{2}^{0}+\tilde{M}_{0}^{0}\right)+\Phi_{2,2} \tilde{M}_{2}^{-}+\Phi_{0,4} \tilde{M}_{4}^{-} \\ =-\left(\alpha_{1,2} \tilde{N}_{1}^{-}+\beta_{0,1} \tilde{N}_{1}^{+}+\beta_{0,3} \tilde{N}_{3}^{-}\right), \\ \cdots \\ \Phi_{m-2, n} \tilde{M}_{m-2}^{+}+\Phi_{m, n}\left(\tilde{M}_{m}^{0}+\tilde{M}_{n}^{0}\right)+\Phi_{m+2, n} \tilde{M}_{m+2}^{-} \\ +\Phi_{m, n-2} \tilde{M}_{n-2}^{+}+\Phi_{m, n+2} \tilde{M}_{n+2}^{-} \\ =-\left(\alpha_{m-1, n} \tilde{N}_{m-1}^{+}+\alpha_{m+1, n} \tilde{N}_{m+1}^{-}+\beta_{m, n-1} \tilde{N}_{n-1}^{+}+\beta_{m, n+1} \tilde{N}_{n+1}^{-}\right) . \end{array}$
From equation (14 ), the unknown coefficients ${{\rm{\Phi }}}_{m,n}$ of the perturbation potential ${\rm{\Phi }}\left(x,y\right)$ can be linearly expressed by means of the coefficients ${\alpha }_{m,n}$ and ${\beta }_{m,n}$ of the transformation electrical field ${E}_{k}^{\ast }\left(x,y\right)$. If the coefficients ${{\rm{\Phi }}}_{m,n}$ are truncated by the $N$-th order (i.e. $m,n=0,1,2,\mathrm{.}.,N$), there are ${\left(1+N\right)}^{2}$ closed equations in equation (14 ) for solving the coefficients ${{\rm{\Phi }}}_{m,n}$. The general solutions of coefficients ${{\rm{\Phi }}}_{m,n}$ can be written as follows,
$\begin{eqnarray}{{\rm{\Phi }}}_{m,n}=\displaystyle \sum _{i,j=0}^{N+1}\left({\gamma }_{m,n,i,j}^{\alpha }{\alpha }_{i,j}+{\gamma }_{m,n,i,j}^{\beta }{\beta }_{i,j}\right),\end{eqnarray}$
where ${\gamma }_{m,n,i,j}^{\alpha }$ and ${\gamma }_{m,n,i,j}^{\beta }$ are the functions of the coefficients ${N}_{m}$.For example, the first-order approximation solutions of the perturbation potential coefficients ${{\rm{\Phi }}}_{m,n}$ (i.e. $m,n=0,1$) can be derived from the first four equations of equation (14 ) by omitting the h.c. ${{\rm{\Phi }}}_{m,n}$ (i.e. $m,n\gt 1$). Therefore, in equation (15 ), the non-zero coefficients ${\gamma }_{m,n,i,j}^{\alpha }$ and ${\gamma }_{m,n,i,j}^{\beta }$ of the first-order approximation solutions ${{\rm{\Phi }}}_{m,n}$ for $N=1$ are given as follows,14 ) by the Mathematica Software. Next, the coefficients ${\alpha }_{m,n}$ and ${\beta }_{m,n}$ of the transformation electrical field will be solved on the basis of equation (15 ).
$\begin{eqnarray*}\begin{array}{c}{\gamma }_{0,0,1,0}^{\alpha }={\gamma }_{0,0,0,1}^{\beta }=-{\tilde{N}}_{1}^{-}/\left(2{\tilde{M}}_{0}^{0}\right),\\ \,{\gamma }_{1,0,0,0}^{\alpha }={\gamma }_{0,1,0,0}^{\beta }=-{\tilde{N}}_{0}^{+}/\left({\tilde{M}}_{0}^{0}+{\tilde{M}}_{1}^{0}\right),\end{array}\end{eqnarray*}$
$\begin{eqnarray*}\begin{array}{c}{\gamma }_{1,0,2,0}^{\alpha }=-{\tilde{N}}_{0}^{-}/\left({\tilde{M}}_{0}^{0}+{\tilde{M}}_{1}^{0}\right),\,\\ {\gamma }_{0,1,1,1}^{\alpha }\,=\,{\gamma }_{1,0,1,1}^{\beta }=-{\tilde{N}}_{1}^{-}/\left({\tilde{M}}_{0}^{0}+{\tilde{M}}_{1}^{0}\right),\end{array}\end{eqnarray*}$
$\begin{eqnarray*}\begin{array}{c}{\gamma }_{0,1,0,2}^{\beta }=-{\tilde{N}}_{2}^{-}/\left({\tilde{M}}_{0}^{0}+{\tilde{M}}_{1}^{0}\right),\,\\ {\gamma }_{1,1,0,1}^{\alpha }={\gamma }_{1,1,1,0}^{\beta }\,=\,-{\tilde{N}}_{0}^{+}/\left(2{\tilde{M}}_{1}^{0}\right),\,{\gamma }_{1,1,2,1}^{\alpha }={\gamma }_{1,1,1,2}^{\beta }=-{\tilde{N}}_{2}^{-}/\left(2{\tilde{M}}_{1}^{0}\right).\end{array}\end{eqnarray*}$
Here, it is noted that the higher order approximation solutions of perturbation potential ${\rm{\Phi }}\left(x,y\right)$ can be easily derived from equation (2.3. Solutions of transformation electric fields
To solve the transformation electrical field, substituting equations (15 ), (7 ) and (8 ) into equation (3 ), we obtain the transformation field equations in the inclusion region,16 ) and (17 ) are equations of transformation electrical fields.
$\begin{eqnarray}\begin{array}{c}{\varepsilon }^{h}{E}_{x}^{* }\left(x,y\right)=\left({\varepsilon }^{h}-{\varepsilon }^{i}\right){E}_{x}^{0}\\ -\left({\varepsilon }^{h}-{\varepsilon }^{i}\right)\displaystyle \sum _{m,n=0}^{\infty }\displaystyle \sum _{i,j}^{N+1}\left({\gamma }_{m,n,i,j}^{\alpha }{\alpha }_{i,j}+{\gamma }_{m,n,i,j}^{\beta }{\beta }_{i,j}\right)\\ \left[{\tilde{N}}_{m}^{+}{u}_{m+1}\left(x\right)+{\tilde{N}}_{m}^{-}{u}_{m-1}\left(x\right)\right]{v}_{n}\left(y\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{c}{\varepsilon }^{h}{E}_{y}^{* }\left(x,y\right)=\left({\varepsilon }^{h}-{\varepsilon }^{i}\right){E}_{y}^{0}\\ -\left({\varepsilon }^{h}-{\varepsilon }^{i}\right)\displaystyle \sum _{m,n=0}^{\infty }\displaystyle \sum _{i,j}^{N+1}\left({\gamma }_{m,n,i,j}^{\alpha }{\alpha }_{i,j}+{\gamma }_{m,n,i,j}^{\beta }{\beta }_{i,j}\right)\\ \left[{\tilde{N}}_{n}^{+}{v}_{n+1}\left(y\right)+{\tilde{N}}_{n}^{-}{v}_{n-1}\left(y\right)\right]{u}_{m}\left(x\right),\end{array}\end{eqnarray}$
where ${\alpha }_{i,j}=\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}{u}_{i}\left(x\right){v}_{j}\left(y\right){E}_{x}^{* }\left(x,y\right){\rm{d}}x{\rm{d}}y$, ${\beta }_{i,j}=\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}{u}_{i}\left(x\right){v}_{j}\left(y\right){E}_{y}^{* }\left(x,y\right){\rm{d}}x{\rm{d}}y$. It is clear that equations (To obtain the solutions of the transformation fields from above equations, the transformation field ${E}_{k}^{* }\left(x,y\right)$ can be given as a series, such as a power series [15, 26, 33],18 ) into equations (16 ) and (17 ), and multiplying both sides of the resulting equations by the term ${\left(x/L\right)}^{k}{\left(y/L\right)}^{l}$, and then integrating these equations over the inclusion region ${{\rm{\Omega }}}_{i}$, we get a set of linear algebra equations about the unknown coefficient ${C}_{k}^{ij}$,19 ) and (20 ),19 ) and (20 ) if the number of equations (19 ) and (20 ), multiplied by the term ${\left(x/L\right)}^{k}{\left(y/L\right)}^{l}$ ($k,l=0,1,2,\mathrm{.}.,S$, and $S$ notes the order of the power-law for the closed equations), equals to the number of the unknown coefficient ${C}_{k}^{ij}$. From above derivations, it is clear that there is not any limitation for the inclusion shapes. The effects of the complex inclusion geometry on the transformation field are expressed by the integrations of Hermite polynomial over the inclusion regions. Thus, our formulas of equations (19 ) and (20 ) can be used to investigate the composites having the inclusions of arbitrary shapes.
$\begin{eqnarray}{E}_{k}^{* }\left(x,y\right)=\displaystyle \sum _{i,j=0}^{\infty }{C}_{k}^{ij}{\left(x/L\right)}^{i}{\left(y/L\right)}^{j},\end{eqnarray}$
where $k=x,y$. ${C}_{k}^{ij}$ is an unknown coefficient. $L$ is the characteristic length of two-dimensional inclusion. Substituting equation ( $\begin{eqnarray}\begin{array}{c}{\varepsilon }^{h}\displaystyle \sum _{i,j=0}^{\infty }{C}_{x}^{ij}{A}^{i+k,j+l}+\left({\varepsilon }^{h}-{\varepsilon }^{i}\right)\displaystyle \sum _{i,j=0}^{\infty }{C}_{x}^{ij}{\tilde{H}}_{\alpha }^{i,j,k,l}\\ \,+\,\left({\varepsilon }^{h}-{\varepsilon }^{i}\right)\displaystyle \sum _{i,j=0}^{\infty }{C}_{{\rm{y}}}^{ij}{\tilde{H}}_{\beta }^{i,j,k,l}=\left({\varepsilon }^{h}-{\varepsilon }^{i}\right){E}_{x}^{0}{A}^{k,l},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{c}{\varepsilon }^{h}\displaystyle \sum _{i,j=0}^{\infty }{C}_{y}^{ij}{A}^{i+k,j+l}+\left({\varepsilon }^{h}-{\varepsilon }^{i}\right)\displaystyle \sum _{i,j=0}^{\infty }{C}_{x}^{ij}{\tilde{F}}_{\alpha }^{i,j,k,l}\\ +\left({\varepsilon }^{h}-{\varepsilon }^{i}\right)\displaystyle \sum _{i,j=0}^{\infty }{C}_{{\rm{y}}}^{ij}{\tilde{F}}_{\beta }^{i,j,k,l}=\left({\varepsilon }^{h}-{\varepsilon }^{i}\right){E}_{y}^{0}{A}^{k,l},\end{array}\end{eqnarray}$
where ${\tilde{H}}_{\alpha }^{i,j,k,l}=\displaystyle \sum _{s,t=0}^{N+1}{H}_{\alpha }^{k,l}\left(s,t\right){G}^{i,j}\left(s,t\right)$, ${\tilde{H}}_{\beta }^{i,j,k,l}=\displaystyle \sum _{s,t=0}^{N+1}{H}_{\beta }^{k,l}\left(s,t\right){G}^{i,j}\left(s,t),\right)$ $\begin{eqnarray*}\begin{array}{c}{H}_{\alpha }^{k,l}\left(s,t\right)=\displaystyle \sum _{m,n=0}^{\infty }{\gamma }_{m,n,s,t}^{\alpha }{D}_{x,m,n}^{k,l},\,\\ {H}_{\beta }^{k,l}\left(s,t\right)\,=\,\displaystyle \sum _{m,n=0}^{\infty }{\gamma }_{m,n,s,t}^{\beta }{D}_{x,m,n}^{k,l},\end{array}\end{eqnarray*}$
$\begin{eqnarray*}\begin{array}{c}{\tilde{F}}_{\alpha }^{i,j,k,l}=\displaystyle \sum _{s,t=0}^{N+1}{G}^{i,j}\left(s,t\right){F}_{\alpha }^{k,l}\left(s,t\right),\,\\ {\tilde{F}}_{\beta }^{i,j,k,l}\,=\,\displaystyle \sum _{s,t=0}^{N+1}{G}^{i,j}\left(s,t\right){F}_{\beta }^{k,l}\left(s,t\right),\end{array}\end{eqnarray*}$
$\begin{eqnarray*}\begin{array}{c}{F}_{\alpha }^{k,l}\left(s,t\right)=\displaystyle \sum _{m,n=0}^{\infty }{\gamma }_{m,n,s,t}^{\alpha }{D}_{y,m,n}^{k,l},\,\\ {F}_{\beta }^{k,l}\left(s,t\right)\,=\,\displaystyle \sum _{m,n=0}^{\infty }{\gamma }_{m,n,s,t}^{\beta }{D}_{y,m,n}^{k,l},\end{array}\end{eqnarray*}$
$\begin{eqnarray*}\begin{array}{c}{D}_{x,m,n}^{k,l}={\tilde{N}}_{m}^{+}{G}^{k,l}\left(m+1,n\right)+{\tilde{N}}_{m}^{-}{G}^{k,l}\left(m-1,n\right),\,\\ \,D{y,m,n}^{k,l}={\tilde{N}}_{n}^{+}{G}^{k,l}\left(m,n+1\right)+{\tilde{N}}_{n}^{-}{G}^{k,l}\left(m,n-1\right),\end{array}\end{eqnarray*}$
${G}^{i,j}\left(s,t\right)=\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}{\left(x/L\right)}^{i}{\left(y/L\right)}^{j}{u}_{s}\left(x\right){v}_{t}\left(y\right){\rm{d}}x{\rm{d}}y$, ${A}^{i,j}=\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}{\left(x/L\right)}^{i}{\left(y/L\right)}^{j}{\rm{d}}x{\rm{d}}y$. Here, note that the following equations are used to derive equations( $\begin{eqnarray}{\alpha }_{m,n}=\displaystyle \sum _{i,j=0}^{\infty }{C}_{x}^{ij}{G}^{i,j}\left(m,n\right),\end{eqnarray}$
$\begin{eqnarray}{\beta }_{m,n}=\displaystyle \sum _{i,j=0}^{\infty }{C}_{y}^{ij}{G}^{i,j}\left(m,n\right).\end{eqnarray}$
Then, the unknown coefficients ${C}_{k}^{ij}$ can be solved with the closed equations (3. Effective dielectric properties
On the basis of the local solutions of perturbation potential and transformation field for a composite having an unbounded host under an external electric field, Maxwell had proposed that the effective dielectric property can be estimated by the electric field averaged over composite volumes which are large compared with the scale of the inhomogeneities [36, 37]. The effective dielectric constant ${\varepsilon }_{ij}^{e}$ of a two-dimensional linear composite can be defined as the follows [38],23 ) and (1 ) into equation (24 ), we have,25 ) shows that the effective responses can be calculated by the volume averages of the electrical field over the inclusion regions.
$\begin{eqnarray}{\bar{D}}_{k}={\varepsilon }_{kx}^{e}{\bar{E}}_{x}+{\varepsilon }_{ky}^{e}{\bar{E}}_{y},\end{eqnarray}$
where $k=x,y$. $\bar{A}=\tfrac{1}{V}\displaystyle {\int }_{\Omega }A{\rm{d}}V$ notes the spatial average of quantity $A$ over the whole composite region ${\rm{\Omega }}$ with the volume $V$. ${E}_{k}$ is the electrical field distribution in the composite system. According to Landau's effective property formula of an isotropic linear composite [38], we have the following formula, $\begin{eqnarray}\begin{array}{c}\displaystyle \frac{1}{V}\displaystyle {\int }_{\Omega }\left({D}_{k}-{\varepsilon }^{h}{E}_{k}\right){\rm{d}}V=\displaystyle \frac{1}{V}\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}\left({D}_{k}^{i}-{\varepsilon }^{h}{E}_{k}\right){\rm{d}}V\\ \,=\,{\bar{D}}_{k}-{\varepsilon }^{h}{\bar{E}}_{k}.\end{array}\end{eqnarray}$
Substituting equations ( $\begin{eqnarray}\begin{array}{c}{\varepsilon }_{kx}^{e}{\bar{E}}_{x}+{\varepsilon }_{ky}^{e}{\bar{E}}_{y}-{\varepsilon }^{h}{\bar{E}}_{k}=\displaystyle \frac{1}{V}\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}\left({\varepsilon }^{i}-{\varepsilon }^{h}\right){E}_{k}{\rm{d}}V\\ \,=\,f\displaystyle \frac{1}{{V}_{{{\rm{\Omega }}}_{i}}}\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}\left({\varepsilon }^{i}-{\varepsilon }^{h}\right){E}_{k}{\rm{d}}V,\end{array}\end{eqnarray}$
where $f={V}_{{{\rm{\Omega }}}_{i}}/V$ is the volume fraction of inclusions, and ${V}_{{{\rm{\Omega }}}_{i}}$ is the volume of inclusions. Equation (Moreover, we regard the external applied field ${E}_{k}^{0}$ as the average field ${\bar{E}}_{k}$. From equation (25 ), for $k=x$, ${\bar{E}}_{x}={E}_{x}^{0}$ and ${\bar{E}}_{y}=0$ (i.e. the external applied electric field ${E}_{x}^{0}$ along $x$-direction), the following formula of effective dielectric constant is given,25 ), we get,26 )-(28 ). Furthermore, based on the equation (3 ), the electrical fields ${E}_{k}^{i}(x,y)$ in the inclusion regions are expressed by the transformation electrical fields,30 ) into equations(26 )-(28 ), the effective dielectric constants can be calculated. In next section, these formulas will be applied to calculate the effective dielectric constants of isotropic composites.
$\begin{eqnarray}{\varepsilon }_{xx}^{e}={\varepsilon }^{h}+f\displaystyle \frac{1}{{E}_{x}^{0}{V}_{{{\rm{\Omega }}}_{i}}}\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}\left({\varepsilon }^{i}-{\varepsilon }^{h}\right){E}_{x}{\rm{d}}V.\end{eqnarray}$
In addition, for $k=x$, ${\bar{E}}_{y}={E}_{y}^{0}$ and ${\bar{E}}_{x}=0$, we have, $\begin{eqnarray}{\varepsilon }_{xy}^{e}=f\displaystyle \frac{1}{{E}_{y}^{0}{V}_{{{\rm{\Omega }}}_{i}}}\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}\left({\varepsilon }^{i}-{\varepsilon }^{h}\right){E}_{x}{\rm{d}}V.\end{eqnarray}$
Similarly, for $k=y$ in equation ( $\begin{eqnarray}{\varepsilon }_{yy}^{e}={\varepsilon }^{h}+f\displaystyle \frac{1}{{E}_{y}^{0}{V}_{{{\rm{\Omega }}}_{i}}}\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}\left({\varepsilon }^{i}-{\varepsilon }^{h}\right){E}_{y}{\rm{d}}V,\end{eqnarray}$
$\begin{eqnarray}{\varepsilon }_{yx}^{e}=f\displaystyle \frac{1}{{E}_{x}^{0}{V}_{{{\rm{\Omega }}}_{i}}}\displaystyle {\int }_{{{\rm{\Omega }}}_{i}}\left({\varepsilon }^{i}-{\varepsilon }^{h}\right){E}_{y}{\rm{d}}V.\end{eqnarray}$
Here, it is noted that both effective dielectric constants ${\varepsilon }_{yx}^{e}$ and ${\varepsilon }_{xy}^{e}$ are symmetric for usual composites. For the composites containing both isotropic host and inclusions, the effective dielectric constants are not always isotropic, and depend on the structures and distributions of inclusions. Clearly, the effects of the inclusion geometry and distribution on the effective responses are shown in equations ( $\begin{eqnarray}{E}_{k}^{i}\left(x,y\right)={E}_{k}^{0}+{E}_{k}^{p}\left(x,y\right)=\displaystyle \frac{{\varepsilon }^{h}}{{\varepsilon }^{h}-{\varepsilon }^{i}}{E}_{k}^{* }\left(x,y\right),\end{eqnarray}$
where $k=x,y$. Thus, substituting equation (4. Numerical discussion
To verify the validity of our method, as an example of two-dimensional problem, the effective dielectric responses of the isotropic cylindrical composite are estimated, where a long cylindrical inclusion with radius $a$ (i.e. radius $a$ within the $x-y$ plane) is embedded in an infinite host. In our numerical calculation, the isotropic dielectric constants of host and inclusion are ${\varepsilon }^{h}=25$ and ${\varepsilon }^{i}=15+n$, respectively, where $n$ is the parameter of the inclusion dielectric constants. The dimensionless radius $a$ of the cylinder inclusion is $a=2$, and the volume fraction of inclusions is $f=0.05$. The origin of Cartesian coordinates $\left(x,y\right)$ is located at the center of cross section of the cylindrical inclusion. Under external applied electrical field ${E}_{x}^{0}$ along the $x$-direction (or ${E}_{y}^{0}$ along the $y$-direction), our method is applied to calculate the effective responses. The results of the method are compared with the exact dilute solutions of cylindrical composites [39].
In table 1, the effective dielectric constants ${\varepsilon }_{xx}^{e}/{\varepsilon }^{h}$ are calculated by TFM on the basis of the zero-order and the first-order approximations of the perturbation potential ${\rm{\Phi }}\left(x,y\right)$ given in equation (15 ), respectively. Good agreements are obtained by comparing the results of TFM with the exact dilute solutions of the cylindrical composites [39]. It is shown that the first-order approximation is more accurate towards exact dilute solutions than the zero-order approximation. In addition, the high order approximation is needed for exactly calculating the effective properties of composites. Here, it should be noted that, in our calculation, the effective dielectric constant ${\varepsilon }_{yy}^{e}$ equals to ${\varepsilon }_{xx}^{e}$ due to the symmetry of the cross section of the cylinder in the $x-y$ plane, and both effective dielectric constant ${\varepsilon }_{xy}^{e}$ and ${\varepsilon }_{yx}^{e}$ are zero. Hence, the numerical results demonstrate that the extended TFM is valid to cope with the effective response problem of composites having arbitrary inclusion geometry since TFM does not have any limitation for the inclusion shapes.
Table 1. Effective dielectric constant contrast ${\varepsilon }_{xx}^{e}/{\varepsilon }^{h}$ estimated by the zero-order and the first-order perturbation potentials for different parameter $n$ of inclusion dielectric constant, where the dielectric constants of the host and cylindrical inclusion are ${\varepsilon }^{h}=25$ and ${\varepsilon }^{i}=15+n$ with the inclusion dielectric parameter $n$, respectively. Here the inclusion volume fraction is $f=0.05$, and 'Exact' notes the exact solution of cylindrical dielectric composite of 39. |
| $n$ | 0-order | First-order | Exact |
|---|---|---|---|
| 1 | 0.98200 | 0.98148 | 0.97805 |
| 2 | 0.98400 | 0.98360 | 0.98095 |
| 3 | 0.98600 | 0.98569 | 0.98372 |
| 4 | 0.98800 | 0.98778 | 0.98636 |
| 5 | 0.99000 | 0.98985 | 0.98889 |
| 6 | 0.99200 | 0.99190 | 0.99130 |
| 7 | 0.99400 | 0.99395 | 0.99362 |
| 8 | 0.99600 | 0.99598 | 0.99583 |
| 9 | 0.99800 | 0.99799 | 0.99796 |
| 10 | 1.00000 | 1.00000 | 1.00000 |
| 11 | 1.00200 | 1.00199 | 1.00196 |
| 12 | 1.00400 | 1.00398 | 1.00385 |
| 13 | 1.00600 | 1.00595 | 1.00566 |
| 14 | 1.00800 | 1.00791 | 1.00741 |
| 15 | 1.01000 | 1.00986 | 1.00909 |
| 16 | 1.01200 | 1.01180 | 1.01071 |
| 17 | 1.01400 | 1.01373 | 1.01228 |
| 18 | 1.01600 | 1.01565 | 1.01379 |
| 19 | 1.01800 | 1.01756 | 1.01525 |
| 20 | 1.02000 | 1.01946 | 1.01667 |
5. Summary
Eshelby's TFM has been developed to solve the fully open boundary problem of composites having arbitrary geometric inclusions. The transformation fields are introduced in the composite system to deal with the complex interface conditions between the inclusion and the host. The transformation electric fields and the perturbation potential are expressed by Hermite polynomials for coping with open boundary conditions. Moreover, a set of linear algebraic equations is given for solving the transformation field on the basis of the transformation field equations. In practical applications, TFM is applied to estimate the effective dielectric responses of the two-dimensional isotropic composites having fully open boundary conditions, and the formulas of effective dielectric constants are expressed by the transformation fields. Furthermore, the validity of TFM is verified by comparing the effective responses calculated by the zero-order and the first-order approximation solutions of perturbation potentials with the exact solutions in the dilute limit. The numerical results show that the method is valid to deal with the open boundary problem of composites. As brief conclusions, it is noted that the higher order solutions of the perturbation potentials given by Hermite polynomials will obtain more accurate results. In addition, TFM can be extended to study the open boundary problem of the anisotropic two-dimensional or three-dimensional composites, such as dielectric, permeability, elastic, thermal, piezoelectric composites, etc [40]. This method will give a new direction of exploring the problems of complex composites.