Shape-dependent initial pressure effects of fluid inclusions in an elastic solid under plane deformation

Abstract

We reconsider the plane deformation of an isotropic elastic solid matrix enclosing macroscale compressible fluid inclusions. The initial pressure inside the fluid inclusions proves to play a significant role in tuning the elastic response of the solid–fluid composite to external loadings. Such initial pressure effects cannot be captured via the classical treatment based on linear elasticity. We follow a modified boundary condition proposed earlier in the literature to examine the influence of the initial pressure in the inclusions on the external loading-induced local and overall elastic behavior of the composite for general shapes of the inclusions. Specifically, we derive initial pressure-dependent closed-form solutions for incremental stress distributions (caused by an in-plane far-field loading) around an isolated fluid inclusion of practically arbitrary shape and then attain initial pressure-dependent in-plane effective properties of the composite containing randomly distributed fluid inclusions (for a given volume fraction) based on certain homogenization methods. Numerical results are presented mainly for a soft elastic solid containing, respectively, approximately regular polygonal liquid inclusions, approximately rectangular liquid inclusions and elliptical liquid inclusions. We show that the initial pressure in the inclusions always enhances the effective moduli of the corresponding composites and that for given volume fraction of the inclusions such enhancement effects are reduced with an increasing number of sides of the inclusions (in terms of regular polygonal cases) while they are amplified with increasing slenderness of the inclusions (in terms of rectangular and elliptical cases).

1. Introduction

How internal inclusions (including pores and rigid inclusions) disturb the stress distribution in the elastic hosting solid is one of the fundamental problems in the micromechanics of porous materials, perforated structures or particle composite materials. In terms of three-dimensional deformations, Sadowsky and Sternberg [1] and Edwards [2] obtained exact closed-form solutions for the stress field in an isotropic matrix containing a spheroidal void or elastic inclusion under uniform remote loading conditions. Their approach employed three harmonic displacement functions formulated in orthogonal curvilinear coordinates. In contrast, Eshelby [3,4] introduced a groundbreaking conceptual framework, now known as Eshelby’s formalism, to analytically address a broader class of problems involving ellipsoidal inclusions. In terms of two-dimensional deformations, Kirsch [5] and Inglis [6] derived explicit analytical solutions for the stress concentration around a small hole in an infinite elastic plane subjected to uniform remote in-plane loadings, addressing the cases of a circular hole and an elliptical hole, respectively. Muskhelishvili [7] developed a powerful complex-variable approach, now referred to as Muskhelishvili’s formalism, and introduced conformal mapping techniques to solve a broader class of elastostatic problems involving holes or inclusions in elastic solids under plane deformation. For two-dimensional inclusion problems (especially for inclusions of general shapes), many analytical solutions are obtained through complex variable methods (e.g., elliptical inclusions [8–10] and arbitrarily shaped inclusions [11–14]).

Fluid inclusions are prevalent in nature and are commonly found within soft or stiff solids, including biological tissues [15], crystals [16], rocks [17], and gels [18]. In recent years, extensive research has been conducted on the elastic behavior of fluid inclusion-filled solids. In terms of macroscale fluid inclusions, Kachanov et al. [19], Shafiro and Kachanov [20], and Wang and Schiavone [21,22] examined the external loading-caused internal pressure in a fluid inclusion and/or the effective moduli of the fluid–solid composite relative to the aspect ratio, orientation and shape of the fluid inclusion; Chen et al. [23] demonstrated the role of the compressibility of a fluid inclusion in tuning the local stress concentration in the surrounding solid; Huang et al. [24] established geometrical nonlinear solutions for the local and overall elastic behavior of a soft, weakly nonlinear solid containing cylindrical liquid inclusions; and most recently Ardeshiri Jouneghani et al. [25] presented an analytic treatment on the case of an isolated liquid inclusion of general shape which is defined by a polynomial mapping function with an arbitrary number of terms. Regarding microscale fluid inclusions, Style et al. [18] provided both experimental and theoretical evidence that small-scale fluid inclusions with surface tension can enhance the stiffness of the hosting soft solid, thereby inducing a hardening effect. Ti et al. [26] further discovered that the surface tension effects of fluid inclusions increase their resistance to deformation while simultaneously amplifying the surrounding stress concentrations (induced by an external loading), while in contrast, Dai et al. [27] identified an almost contrary result that in the existing presence of surface elasticity of fluid inclusions, additional incorporation of surface tension of fluid inclusions would reduce the surrounding stress concentrations.

In general, fluid inclusions typically possess an initial pressure. For example, gas inclusions inherently possess an initial pressure, while liquid inclusions may also develop an initial pressure when inserted into a hosting solid during the preparation process. From a fundamental perspective, the initial pressure inevitably influences the overall mechanical behavior of corresponding fluid-solid composite materials. For example, a basketball would be stiffer when inflated somewhat (this is mainly not attributed to the fact that the internal air has a bulk modulus since the bulk modulus of the internal air is identical to the internal pressure which is, however, much smaller than the modulus of the external rubber or leather material). Dai et al. [28] proposed modified boundary conditions capable of capturing the effects of initial pressure on the incremental elastic behavior of fluid inclusion-filled solids (induced by external loadings) within the framework of linear elasticity (conventional boundary conditions incorporating initial pressure in linear elasticity fail to capture any initial-pressure effect on the incremental local or overall behavior). These modified boundary conditions or the essential ideas behind them were later used in exploring the initial pressure effects in specific cases of spherical fluid inclusions [29] and elliptical cylindrical liquid inclusions [30], respectively.

This paper aims to further investigate the initial pressure effects using the modified boundary conditions proposed in Dai et al. [28], by extending the analysis to fluid inclusions of arbitrary geometry. In fact, this work is stimulated by the consideration that realistic shapes of fluid inclusions existing in an elastic solid are usually out of the common regular groups of inclusion shapes (such as circular, elliptical, and spherical) leading to some unique initial pressure effects (which are inherently shape-dependent). Here, we particularly mention that an accurate analysis of the initial pressure effects for fluid inclusions in an elastic solid (especially a soft one) usually needs to resort to the finite deformation framework or at least small-on-large theories [31–33]. In terms of fluid inclusions of general shapes, however, the use of finite deformation or small-on-large theories precludes almost any analytical treatment of corresponding boundary value problems and any possibility of attaining concise analytic solutions. In fact, to obtain explicit referential results for fluid inclusions in a soft elastic solid, the linear elasticity theory with certain slight modifications has remained the most preferred approximate framework [21–23,26–30,34] despite the fact that the hosting soft solid may behave nonlinearly and undergo relatively large deformations.

We arrange our paper as follows. In section 2, we formulate the boundary value problem for the incremental plane deformation of a macroscale, compressible, arbitrarily shaped fluid inclusion with initial pressure in an elastic solid. In section 3, we determine explicitly the full incremental elastic field in the solid when an in-plane constant far-field loading is applied, while in section 4 we derive the effective properties of the fluid-solid composite system based on certain effective medium theories. In section 5, we discuss the local incremental stress concentration and overall effective moduli for several common non-circular inclusion shapes via a few numerical examples. Finally, we summarize the main results in section 6.

2. Problem description

We consider a macroscale linearly compressible fluid inclusion embedded in an elastic solid matrix subjected to plane deformation. The fluid inclusion has an intrinsic initial pressure after the fluid–solid composite structure is produced but before any external loading is imposed on the structure. The (initial) shape of the inclusion is assumed to be practically arbitrary as defined by a polynomial mapping function (see equation (16)). The purpose here is then to analytically determine the local incremental stress distribution around the inclusion when a constant far-field loading is applied to the matrix, as well as the effective in-plane response of the fluid-solid composite to the applied far-field loading based on some homogenization schemes. Rigorously speaking, an exact determination of the incremental stress field induced by the far-field loading requires the use of finite deformation theories which, however, would reduce significantly the possibility of an analytic study. Consequently, we still retain the overall framework of small deformation theories but will use a modified boundary condition on the fluid–solid interface to capture some of the initial pressure effects on the incremental response of the structure to the far-field loading.

We refer to the rectangular x-y coordinate system and illustrate the essential problem in Figure 1. The interface curve between the inclusion and matrix is denoted by L, while the unbounded area encompassed by the matrix is represented by S. Specifically, T and N signify the counterclockwise tangent and outward-pointing normal to the curve L, respectively, while $\alpha$ represents the angle between the T-direction and the positive x-axis. The in-plane rectangular components of the stress and displacement of the matrix are represented, respectively, by ( ${\sigma }_{\mathrm{xx}}$ , ${\sigma }_{\mathrm{yy}}$ , ${\sigma }_{\mathrm{xy}}$ ) and ( $u_x$ , $u_y$ ), while the in-plane components of the local stress and displacement of the matrix on the interface, with reference to the N-T coordinate system, are denoted by ( ${\sigma }_{\mathrm{NN}}$ , ${\sigma }_{\mathrm{TT}}$ , ${\sigma }_{\mathrm{NT}}$ ) and ( $u_N$ , $u_T$ ). The initial pressure and bulk modulus of the inclusion are denoted by $P_0$ and K. The rectangular components of the far-field loading are represented by ( ${\sigma }_{\mathrm{xx}}^{\infty }$ , ${\sigma }_{\mathrm{yy}}^{\infty }$ , ${\sigma }_{\mathrm{xy}}^{\infty }$ ).

Figure 1.

A fluid inclusion with initial pressure in an elastic solid under remote loading.

The interface between the inclusion and matrix is deformed immediately after the far-field loading is imposed on the matrix remotely, leading to a certain change in the volume of the inclusion and therefore a change in the fluid pressure. The classical boundary condition describing the fluid pressure of the inclusion acting on the internal boundary of the matrix is expressed, in terms of the normal and tangential stresses ( ${\sigma }_{\mathrm{NN}}$ and ${\sigma }_{\mathrm{NT}}$ ) of the matrix, as

$${\sigma }_{\mathrm{NN}}=-(P_0+\mathrm{\Delta }P),{\sigma }_{\mathrm{NT}}=0,\mathrm{on}L,$$

(1)

where $P_0$ is the initial pressure of the inclusion (before the far-field loading is applied), while $\mathrm{\Delta }P$ is the change in the pressure of the inclusion after deformation defined by ( $\mathrm{\Delta }V/V$ refers to the rate of change in the volume of the inclusion)

$$\mathrm{\Delta }P=-K\mathrm{\Delta }V/V.$$

(2)

We note from equation (1) that the initial pressure-induced traction on the matrix is specific and fixed independent of the far-field loading, indicating the fact that with the use of the classical boundary condition the initial pressure of the inclusion would play no part in determining both the incremental elastic field (induced by the far-field loading) in the matrix and the effective response of the entire composite structure to the far-field loading. However, such a vanishing role of initial pressure of the inclusion, related to the use of classical boundary condition (1), in determining the incremental elastic response of a perforated or porous structure is inconsistent with common sense. For example, a basketball would be stiffer when inflated somewhat (this is mainly not attributed to the fact that the internal air has a bulk modulus since the bulk modulus of the internal air is identical to the internal pressure which is, however, much smaller than the modulus of the external rubber or leather material). To address this deficiency to some extent for the current case of a fluid inclusion with initial pressure, a modified boundary condition was suggested in the work by Dai et al. [28] to account for the effects of the directional change in the internal pressure of the inclusion during deformation since, from a mathematical point of view, the change in the volume of the inclusion during deformation is of the same order of that in the normal direction of the interface curve (in which the internal pressure acts on the surrounding matrix). This modified boundary condition may be taken as [28]

$${\sigma }_{\mathrm{NN}}=-(P_0+P_0{\epsilon }_{\mathrm{TT}}+\mathrm{\Delta }P),{\sigma }_{\mathrm{NT}}=-P_0\omega ,\mathrm{on}L,$$

(3)

where $\mathrm{\Delta }P$ is defined in equation (2), while ${\epsilon }_{\mathrm{TT}}$ and $\omega$ denote the stretching strain and the rotation angle of an infinitesimal element of the curve L during deformation, respectively. Here, the ${\epsilon }_{\mathrm{TT}}$ -related term may alternatively be removed from equation (3) (in this case, the first part of equation (3) becomes identical to that of equation (1)), but for the sake of establishing a more mathematically consistent boundary value problem, we retain the term in equation (3) since for a general inclusion shape (especially irregular shapes) the resultant traction (induced by the initial pressure) imposed on the matrix following equation (3) would be possibly not self-equilibrated if the ${\epsilon }_{\mathrm{TT}}$ -related term is disregarded. Clearly, since the term $P_0$ on the right side of the first part of equation (3), which is independent of the incremental deformation induced by the far-field loading, plays no part in determining the incremental elastic field in the matrix caused by the far-field loading, it is justified to use a simplified version of equation (3) (by removing this term),

$${\sigma }_{\mathrm{NN}}=-(P_0{\epsilon }_{\mathrm{TT}}+\mathrm{\Delta }P),{\sigma }_{\mathrm{NT}}=-P_0\omega ,\mathrm{on}L,$$

(4)

for our purpose of studying the initial pressure effects on the incremental elastic behavior of the fluid-solid structure.

We introduce the Muskhelishvili’s formalism of plane elasticity to formulate the current boundary value problem. In this formalism, the rectangular components of the stress and displacement in the matrix can be represented via two complex analytic functions $\varphi (z)$ and $\psi (z)$ (z is a complex variable describing the position of a point in the matrix), as [7]

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\mathrm{xx}}+{\sigma }_{\mathrm{yy}}=2[{\varphi }^{\prime }(z)+\overline{{\varphi }^{\prime }(z)}],{\sigma }_{\mathrm{yy}}-{\sigma }_{\mathrm{xx}}+2i{\sigma }_{\mathrm{xy}}=2[\overline{z}{\varphi }^{\prime \prime }(z)+{\psi }^{\prime }(z)],}\end{array}$$

(5)

$$2G(u_x+i{u}_y)=\kappa \varphi (z)-z\overline{{\varphi }^{\prime }(z)}-\overline{\psi (z)},$$

(6)

with

$$\kappa =\{\begin{array}{cc}3-4\nu ,& planestrain,\\ \left(3-\nu \right)/\left(1+\nu \right),& planestress,\end{array}$$

(7)

where $z=x+\mathrm{iy}$ and i is the imaginary unit, while G and $\nu$ are, respectively, the shear modulus and Poisson’s ratio of the matrix. In the context of the Muskhelishvili’s formalism, we are able to represent equations (4) and (2) in an integral form as (refer to Sections 2 and 3 of Dai et al. [28] for details of the complete mathematical formulation)

$$\varphi (t)+t\overline{{\varphi }^{\prime }(t)}+\overline{\psi (t)}=-P^{*}[\kappa \varphi (t)-t\overline{{\varphi }^{\prime }(t)}-\overline{\psi (t)}]-\mathrm{\Delta }\mathrm{Pt}+C,\forall t\in L,$$

(8)

and

$$\mathrm{\Delta }P=-K^{*}\frac{\mathrm{Im}\left({\int }_L[\kappa \overline{\varphi (t)}-\overline{t}{\varphi }^{\prime }(t)-\psi (t)]\mathrm{dt}\right)}{\mathrm{Im}\left({\int }_L\overline{t}\mathrm{dt}\right)},$$

(9)

where C is the constant of integration and the dimensionless parameters $P^{*}$ and $K^{*}$ are defined by

$$P^{*}=\frac{P_0}{2G},K^{*}=\frac{K}{G}.$$

(10)

In terms of a constant far-field loading applied at infinity, the analytic functions $\varphi (z)$ and $\psi (z)$ can be decomposed into

$$\varphi (z)={\mathrm{\Gamma }}_{1}z+{\varphi }_0(z),\psi (z)={\mathrm{\Gamma }}_{2}z+{\psi }_0(z),$$

(11)

where

$${\mathrm{\Gamma }}_1=\frac{{\sigma }_{\mathrm{xx}}^{\infty }+{\sigma }_{\mathrm{yy}}^{\infty }}{4},{\mathrm{\Gamma }}_2=\frac{{\sigma }_{\mathrm{yy}}^{\infty }-{\sigma }_{\mathrm{xx}}^{\infty }}{2}+i{\sigma }_{\mathrm{xy}}^{\infty },$$

(12)

while ${\varphi }_0(z)$ and ${\psi }_0(z)$ are two new complex analytic functions still defined in the matrix domain S (specifically we may stipulate, without loss of generality, that they vanish at infinity). Substituting equation (11) into equations (8) and (9), we obtain

$$(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)t+A{\varphi }_0(t)+B[t\overline{{\varphi }_0^{\prime }(t)}+\overline{{\mathrm{\Gamma }}_2}\overline{t}+\overline{{\psi }_0(t)}]-C=0,\forall t\in L,$$

(13)

and

$$\mathrm{\Delta }P=-K^{*}\frac{\mathrm{Im}\left({\int }_L[(\kappa -1){\mathrm{\Gamma }}_1\overline{t}+\kappa \overline{{\varphi }_0(t)}-\overline{t}{{\varphi }_0}^{\prime }(t)-{\psi }_0(t)]\mathrm{dt}\right)}{\mathrm{Im}\left({\int }_L\overline{t}\mathrm{dt}\right)},$$

(14)

where the dimensionless parameters A, B, and D are defined by

$$A=1+\kappa P^{*},B=1-P^{*},D=2+(\kappa -1)P^{*}.$$

(15)

In what follows, we present how to determine ${\varphi }_0(z)$ and ${\psi }_0(z)$ analytically from equations (13) with (14) for a general shape of the fluid inclusion which is defined by a polynomial mapping function.

3. Solution procedure

We introduce a conformal mapping function to describe the geometry of the inclusion, i.e., Muskhelishvili [7]

$$z=w(\xi )=R{e}^{i\theta }(\xi +\sum _{m=1}^n{a}_m{\xi }^{-m}),\left|\xi \right|\ge 1,z\in S,$$

(16)

where the real constant R and the complex constants $a_m$ characterize the size and shape of the inclusion, respectively, while $\theta$ determines the orientation of the inclusion. Here, we are mainly concerned with some smooth shapes of the inclusion which have no strictly sharp corners. Theoretical inclusion shapes with strictly sharp corners may be approximately treated as their smooth counterparts with chamfered corners of sufficient sharpness (from an engineering point of view). Equation (16) maps not only the curve L in the z-plane onto the unit circle $\gamma$ (denoted by $\sigma$ ) in the $\xi$ -plane, but also the entire S in the z-plane to the entire exterior of the unit circle $\gamma$ in the $\xi$ -plane. Particularly, the bijection of the conformal mapping equation (16) (for a smooth inclusion shape without strictly sharp corners) requires that its derivative $w^{\prime }(\xi )$ has no roots on or outside the unit circle in the $\xi$ -plane, i.e.,

$$w^{\prime }(\xi )\ne 0,\forall \left|\xi \right|\ge 1.$$

(17)

In the context of equation (16), ${\varphi }_0(z)$ and ${\psi }_0(z)$ in S can be regarded as two specific functions of $\xi$ that are analytic outside the unit circle in the $\xi$ -plane, allowing us to denote

$$f_1(\xi )\stackrel{\mathrm{\Delta }}{=}{\varphi }_0(w(\xi ))\equiv {\varphi }_0(z),f_2(\xi )\stackrel{\mathrm{\Delta }}{=}{\psi }_0(w(\xi ))\equiv {\psi }_0(z),z\in S,\left|\xi \right|\ge 1.$$

(18)

Furthermore, we can expand them into the Laurent series in $\xi$ as follows:

$$f_1(\xi )=\sum _{m=1}^{\infty }{\alpha }_{-m}{\xi }^{-m},f_2(\xi )=\sum _{m=1}^{\infty }{\beta }_{-m}{\xi }^{-m},\left|\xi \right|\ge 1.$$

(19)

Here, as explained hereinabove (see the text following equation (12)), $f_1(\xi )$ and $f_2(\xi )$ need to vanish at infinity so that equation (19) contains no positive-exponent or constant terms. Using equations (16) and (18), equations (13) and (14) are rewritten on the unit circle $\gamma$ in the $\xi$ -plane as

$$\begin{array}{c}\multicolumn{1}{c}{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)w(\sigma )+A{f}_1(\sigma )+}\\ \multicolumn{1}{c}{B[\frac{w(\sigma )}{\overline{w^{\prime }(\sigma )}}\overline{{f_1}^{\prime }(\sigma )}+\overline{{\mathrm{\Gamma }}_2}\overline{w(\sigma )}+\overline{f_2(\sigma )}]-C=0,\forall \sigma \in \gamma ,}\end{array}$$

(20)

with

$$\mathrm{\Delta }P=-K^{*}\frac{\mathrm{Im}({\int }_{\gamma }[(\kappa -1){\mathrm{\Gamma }}_1\overline{w(\sigma )}+\kappa \overline{f_1(\sigma )}-\frac{\overline{w(\sigma )}}{w^{\prime }(\sigma )}{f_1}^{\prime }(\sigma )-f_2(\sigma )]w^{\prime }(\sigma )d\sigma )}{\mathrm{Im}({\int }_{\gamma }\overline{w(\sigma )}{w}^{\prime }(\sigma )d\sigma )},$$

(21)

where $\sigma =e^{i\beta }(0\le \beta <2\pi )$ represents $\xi$ when $\xi$ falls on the unit circle $\gamma$ .

We will apply Cauchy integral techniques to obtain analytic solutions for $f_1(\xi )$ and $f_2(\xi )$ from the boundary condition (20). First, we need to explore the series expansion of $\overline{{f_1}^{\prime }(\sigma )}w(\sigma )/\overline{w^{\prime }(\sigma )}$ on the unit circle. To this end, we simplify $w(\sigma )/\overline{w^{\prime }(\sigma )}$ into

$$\frac{w(\sigma )}{\overline{w^{\prime }(\sigma )}}=e^{2i\theta }\frac{(\sigma +a_1{\sigma }^{-1}+\cdots +a_n{\sigma }^{-n})}{(1-\overline{a_1}{\sigma }^2-\cdots -n\overline{a_n}{\sigma }^{n+1})}=e^{2i\theta }{\sigma }^{-n}\frac{({\sigma }^{n+1}+a_1{\sigma }^{n-1}+\cdots +a_n)}{(1-\overline{a_1}{\sigma }^2-\cdots -n\overline{a_n}{\sigma }^{n+1})}.$$

(22)

Here, in equation (22), $1/\overline{w^{\prime }(\sigma )}$ can be treated as the boundary value of $1/{\overline{w}\prime }^{-1}$ (on the unit circle) which is analytic in the entire disk defined by $\left|\eta \right|\le 1$ (noting from equation (17) that ${\overline{w}\prime }^{-1}$ in the entire disk) and therefore can be always expanded into a standard Taylor series in $\eta$ with no negative-exponent terms. As a result, the Fourier series expansion of $1/\overline{w\prime (\sigma )}$ on the unit circle contains certainly no negative-exponent terms even if $|1-\overline{\omega \prime (\sigma )}|=\left|\overline{a_1}{\sigma }^2-\cdots -n\overline{a_n}{\sigma }^{n+1}\right|\ge 1$ for some points on the unit circle. Consequently, the general form of the Fourier series expansion for equation (22) turns out to be

$$\frac{w(\sigma )}{\overline{w^{\prime }(\sigma )}}=e^{2i\theta }(b_{-n}{\sigma }^{-n}+b_{-n+1}{\sigma }^{-n+1}+\cdots +b_{-1}{\sigma }^{-1}+b_0+\sum _{m=1}^{\infty }{b}_m{\sigma }^m).$$

(23)

We further have

$$\frac{w(\sigma )}{\overline{w^{\prime }(\sigma )}}\overline{{f_1}^{\prime }(\sigma )}=e^{2i\theta }(K_{-1}{\sigma }^{-1}+K_{-2}{\sigma }^{-2}+\cdots +K_{-n+2}{\sigma }^{-n+2}+K_0+F_1(\sigma )),$$

(24)

where

$$\begin{array}{c}\multicolumn{1}{c}{K_0=-\overline{{\alpha }_{-1}}{b}_{-2}-2\overline{{\alpha }_{-2}}{b}_{-3}-\cdots -(n-1)\overline{{\alpha }_{-n+1}}{b}_{-n},}\\ \multicolumn{1}{c}{K_{-1}=-\overline{{\alpha }_{-1}}{b}_{-3}-2\overline{{\alpha }_{-2}}{b}_{-4}-\cdots -(n-2)\overline{{\alpha }_{-n+2}}{b}_{-n},}\\ \multicolumn{1}{c}{K_{-2}=-\overline{{\alpha }_{-1}}{b}_{-4}-2\overline{{\alpha }_{-2}}{b}_{-5}-\cdots -(n-3)\overline{{\alpha }_{-n+3}}{b}_{-n},}\\ \multicolumn{1}{c}{...}\\ \multicolumn{1}{c}{K_{-n+3}=-\overline{{\alpha }_{-1}}{b}_{-n+1}-2\overline{{\alpha }_{-2}}{b}_{-n},}\\ \multicolumn{1}{c}{K_{-n+2}=-\overline{{\alpha }_{-1}}{b}_{-n},}\end{array}$$

(25)

while $F_1(\sigma )$ is a function involving only positive-exponent terms of $\sigma$ . We then apply the operator $\frac{1}{2\pi i}\underset{\gamma }{\oint }d\sigma /(\sigma -\xi )$ ( $\left|\xi \right|>1$ ) to both sides of equation (20) and obtain

$$\begin{array}{c}\multicolumn{1}{c}{B\overline{{\mathrm{\Gamma }}_2}R{e}^{-i\theta }{\xi }^{-1}+(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)R{e}^{i\theta }\sum _{m=1}^n{a}_m{\xi }^{-m}}\\ \multicolumn{1}{c}{+A\sum _{m=1}^{\infty }{\alpha }_{-m}{\xi }^{-m}+(1-P^{*})e^{2i\theta }\sum _{m=1}^{n-2}{K}_{-m}{\xi }^{-m}=0,}\end{array}$$

(26)

indicating that

$$\begin{array}{c}\multicolumn{1}{c}{A{\alpha }_{-1}+B{e}^{2i\theta }{K}_{-1}=-B\overline{{\mathrm{\Gamma }}_2}R{e}^{-i\theta }-(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)R{e}^{i\theta }{a}_1}\\ \multicolumn{1}{c}{A{\alpha }_{-m}+B{e}^{2i\theta }{K}_{-m}=-(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)R{e}^{i\theta }{a}_m(m=2,\dots ,n-2),}\\ \multicolumn{1}{c}{A{\alpha }_{-m}=-(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)R{e}^{i\theta }{a}_m(m=n-1,n),}\\ \multicolumn{1}{c}{A{\alpha }_{-m}=0(m>n),}\end{array}$$

(27)

and $f_1(\xi )$ reduces to a finite-term series as

$$f_1(\xi )=\sum _{m=1}^n{\alpha }_{-m}{\xi }^{-m}.$$

(28)

Here, all the non-vanishing ${\alpha }_{-m}$ may be determined from solving the combination of equations (25) and (27). After determining $f_1(\xi )$ we may similarly attain $f_2(\xi )$ from the conjugate of equation (20) which is

$$\begin{array}{c}\multicolumn{1}{c}{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)\overline{w(\sigma )}+A\overline{f_1(\sigma )}+}\\ \multicolumn{1}{c}{B[\frac{\overline{w(\sigma )}}{w^{\prime }(\sigma )}{f_1}^{\prime }(\sigma )+{\mathrm{\Gamma }}_{2}w(\sigma )+f_2(\sigma )]-C=0,\forall \sigma \in \gamma ,}\end{array}$$

(29)

We see from the conjugate of equation (24) that

$$\frac{\overline{w(\sigma )}}{w^{\prime }(\sigma )}{f_1}^{\prime }(\sigma )=e^{-2i\theta }(\overline{K_{-1}}\sigma +\overline{K_{-2}}{\sigma }^2+\cdots +\overline{K_{-n+2}}{\sigma }^{n-2}+\overline{K_0}+\overline{F_1(\sigma )}),$$

(30)

where $\overline{F_1(\sigma )}$ is a function involving only negative-exponent terms of $\sigma$ . It is implied from equation (30) that $\overline{w}({\xi }^{-1}){f_1}^{\prime }(\xi )/w^{\prime }(\xi )$ ( $\left|\xi \right|\ge 1$ ), whose boundary value on the unit circle in the $\xi$ -plane is $\overline{w(\sigma )}{f_1}^{\prime }(\sigma )/w^{\prime }(\sigma )$ , is analytic in the entire exterior of the unit circle in the $\xi$ -plane (except for the infinity point where it has a (n-2)th-order pole). Consequently, applying the Cauchy integral operator to the equation (29) yields

$$\begin{array}{c}\multicolumn{1}{c}{f_2(\xi )=-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}R{e}^{-i\theta }{\xi }^{-1}-{\mathrm{\Gamma }}_{2}R{e}^{i\theta }\sum _{m=1}^n{a}_m{\xi }^{-m}}\\ \multicolumn{1}{c}{-[\frac{\overline{w(\xi )}}{w^{\prime }(\xi )}{f_1}^{\prime }(\xi )-e^{-2i\theta }\sum _{m=0}^{n-2}\overline{K_{-m}}{\xi }^m].}\end{array}$$

(31)

Now we have determined analytically $f_1(\xi )$ and $f_2(\xi )$ in the entire exterior of the unit circle in the $\xi$ -plane (or equivalently ${\varphi }_0(z)$ and ${\psi }_0(z)$ in the entire S).

In addition, substituting $f_1(\sigma )$ , $f_2(\sigma )$ , and $w(\sigma )$ into equation (21), we readily determine the internal pressure inside the inclusion as

$$\mathrm{\Delta }P=\frac{-K^{*}\left(\begin{array}{c}(\kappa -1){\mathrm{\Gamma }}_1{R}^2(1-{\sum }_{m=1}^{n}m{a}_m\overline{a_m})+\frac{\left(D{\mathrm{\Gamma }}_1{R}^2\right)}{B}-\mathrm{Re}\left(\kappa R{e}^{i\theta }{\sum }_{m=1}^{n}m{a}_m\overline{{\alpha }_{-m}}\right)\hfill \\ \multicolumn{1}{c}{+\mathrm{Re}\left({\mathrm{\Gamma }}_2{R}^2{e}^{2i\theta }{a}_1\right)+\mathrm{Re}\left(R{e}^{-i\theta }{\sum }_{m=1}^{n-2}\overline{K_{-m}}m{a}_m\right)}\end{array}\right)}{[(1-{\sum }_{m=1}^{n}m{a}_m\overline{a_m})+\frac{K^{*}}{1-P^{*}}]R^2}$$

(32)

where ${\alpha }_{-m}$ and $K_{-m}$ are already determined (in advance) from the combination of equations (25) and (27).

3.1. Hypotrochoidal and elliptical fluid inclusions

The current solution degenerates readily into some completely explicit results for hypotrochoidal fluid inclusions and elliptical fluid inclusions.

For a hypotrochoidal fluid inclusion, the corresponding mapping function (16) reduces to

$$z=w(\xi )=R{e}^{i\theta }(\xi +a_n{\xi }^{-n}),a_n\in \mathbb{R},n\ge 2,\left|\xi \right|\ge 1.$$

(33)

Substituting equation (33) into equations (22)–(28), the function $f_1(\xi )$ is simplified into

$$f_1(\xi )=-\frac{E\overline{{\mathrm{\Gamma }}_2}R{e}^{-i\theta }}{(1-E^2(n-2)a_n^2)}{\xi }^{-1}-\frac{E^2{\mathrm{\Gamma }}_2{a}_{n}R{e}^{3i\theta }}{(1-E^2(n-2)a_n^2)}{\xi }^{-n+2}-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}R{e}^{i\theta }{a}_n{\xi }^{-n},$$

(34)

where the dimensionless parameter $E$ is defined by

$$E=\frac{1-P^{*}}{1+\kappa P^{*}}.$$

(35)

Applying equations (33) and (34) to equation (32) and using ${\sigma }_{\mathrm{TT}}\approx 4\mathrm{Re}[{\mathrm{\Gamma }}_1+\frac{{f_1}^{\prime }(\xi )}{w^{\prime }(\xi )}]+\mathrm{\Delta }P$ , we identify the incremental hoop stress in the matrix around the inclusion and the far-field loading-induced pressure change in the inclusion as

$${\sigma }_{\mathrm{TT}}\approx 4\mathrm{Re}[{\mathrm{\Gamma }}_1+\frac{\begin{array}{c}E\overline{{\mathrm{\Gamma }}_2}{e}^{-2i\theta }{\sigma }^{-2}+(n-2)E^2{\mathrm{\Gamma }}_2{a}_n{e}^{2i\theta }{\sigma }^{-n+1}+\hfill \\ \multicolumn{1}{c}{\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}(1-E^2(n-2)a_n^2)n{a}_n{\sigma }^{-n-1}}\end{array}}{(1-n{a}_n{\sigma }^{-n-1})(1-E^2(n-2)a_n^2)}]+\mathrm{\Delta }P,$$

(36)

$$\mathrm{\Delta }P=\frac{-{\mathrm{\Gamma }}_1(\kappa +\frac{1}{E})(1+{\mathrm{Ena}}_n^2)}{[\frac{1}{B}+\frac{(1-{\mathrm{na}}_n^2)}{K^{*}}+\kappa {\mathrm{na}}_n^2\frac{1}{A}]},$$

(37)

where A and B are defined from equation (15). Surprisingly, we find that the modified results here for a hypotrochoidal fluid inclusion incorporating initial pressure effects have the same structure as the classical counterparts derived in Wang and Schiavone [21] (the only differences lie in that the parameters A, B, and E are related to the normalized initial pressure parameter $P^{*}$ in the modified case while they all reduce to 1, i.e., A = B = E = 1, in the classical case). For example, letting A = B = E = 1, we obtain immediately

$$\mathrm{\Delta }P=\frac{-{\mathrm{\Gamma }}_1(\kappa +1)(1+{\mathrm{na}}_n^2)}{[1+\frac{(1-{\mathrm{na}}_n^2)}{K^{*}}+\kappa {\mathrm{na}}_n^2]},$$

(38)

which is consistent with corresponding classical result for pressure change (see equation (37) in the work by Wang and Schiavone [21]).

For an elliptical fluid inclusion, the corresponding mapping function (16) reduces to

$$z=w(\xi )=R{e}^{i\theta }(\xi +a_1{\xi }^{-1}),a_1=\frac{a-b}{a+b},R=\frac{a+b}{2},\left|\xi \right|\ge 1,$$

(39)

where a and b are the semi-major and semi-minor axes of the inclusion, respectively. Following similar procedures of simplification, we attain $f_1(\xi )$ , ${\sigma }_{\mathrm{TT}}$ and $\mathrm{\Delta }P$ as

$$f_1(\xi )=-[E\overline{{\mathrm{\Gamma }}_2}R{e}^{-i\theta }+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}R{e}^{i\theta }{a}_1]{\xi }^{-1},$$

(40)

$${\sigma }_{\mathrm{TT}}\approx 4\mathrm{Re}[{\mathrm{\Gamma }}_1+\frac{[E\overline{{\mathrm{\Gamma }}_2}{e}^{-2i\theta }+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}{a}_1]{\sigma }^{-2}}{(1-a_1{\sigma }^{-2})}]+\mathrm{\Delta }P$$

(41)

$$\mathrm{\Delta }P=-\frac{[(\kappa -1){\mathrm{\Gamma }}_1(1-a_1^2)+(\kappa E+1)a_1\mathrm{Re}\left({\mathrm{\Gamma }}_2{e}^{2i\theta }\right)+(\frac{\kappa a_1^2}{A}+\frac{1}{B})D{\mathrm{\Gamma }}_1]}{(\frac{1}{B}+\frac{(1-a_1^2)}{K^{*}}+\frac{\kappa a_1^2}{A})},$$

(42)

the last part of which agrees essentially with the result derived in the work by Sun and Dai [30] for an initially pressurized elliptical inclusion (see equation (26) there).

4. Effective properties of the composite material

We consider a solid-based composite material containing a large number of fluid inclusions and derive its macroscopic effective properties based on the effective medium theories. The effective compliance matrix of the composite material is represented as [35]

$$H_{\mathrm{eff}}=H+f_{\mathrm{in}}(A_c-H{B}_c),$$

(43)

where H is the compliance matrix of the solid surrounding the fluid inclusions defined as

$$H=\left[\begin{array}{ccc}\frac{1+\kappa }{8G}\hfill & \frac{\kappa -3}{8G}\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{\kappa -3}{8G}}& \frac{1+\kappa }{8G}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1}{2G}\hfill \end{array}\right],$$

(44)

and $f_{\mathrm{in}}$ is the volume fraction of the fluid inclusions, while matrices $A_c$ and $B_c$ represent concentration coefficients defining the average strain ( ${\epsilon }_{\mathrm{ij}}^{\mathrm{in}}$ , i,j = x, y) and stress ( ${\sigma }_{\mathrm{ij}}^{\mathrm{in}}$ , i,j = x, y) within the fluid inclusions as defined by

$$\left[\begin{array}{c}{\epsilon }_{\mathrm{xx}}^{\mathrm{in}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{yy}}^{\mathrm{in}}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{xy}}^{\mathrm{in}}}\end{array}\right]=A_c\left[\begin{array}{c}{\sigma }_{\mathrm{xx}}^{\infty }\hfill \\ \multicolumn{1}{c}{{\sigma }_{\mathrm{yy}}^{\infty }}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xy}}^{\infty }}\end{array}\right],\left[\begin{array}{c}{\sigma }_{\mathrm{xx}}^{\mathrm{in}}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\mathrm{yy}}^{\mathrm{in}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xy}}^{\mathrm{in}}}\end{array}\right]=B_c\left[\begin{array}{c}{\sigma }_{\mathrm{xx}}^{\infty }\hfill \\ \multicolumn{1}{c}{{\sigma }_{\mathrm{yy}}^{\infty }}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xy}}^{\infty }}\end{array}\right].$$

(45)

Specifically, the dilute method and the Mori–Tanaka method are used to evaluate the concentration coefficient matrices. In the dilute method, the interaction between fluid inclusions is ignored, and therefore the concentration coefficient matrices, denoted by $A_c^d$ and $B_c^d$ , may be obtained from the current solution for a single isolated inclusion. In the Mori-Tanaka method, the concentration coefficient matrices, denoted by $A_c^{\mathrm{MT}}$ and $B_c^{\mathrm{MT}}$ , are simply established on its dilute counterparts as [28,35]

$$A_c^{\mathrm{MT}}=A_c^d{[(1-f_{\mathrm{in}})I+f_{\mathrm{in}}{B}_c^d]}^{-1},B_c^{\mathrm{MT}}=B_c^d{[(1-f_{\mathrm{in}})I+f_{\mathrm{in}}{B}_c^d]}^{-1},$$

(46)

where I is a $3\times 3$ identity matrix. Obviously, $A_c^d$ and $B_c^d$ play an essential role in determining the effective properties of the composite material in either of the dilute and Mori–Tanaka methods.

We derive $A_c^d$ and $B_c^d$ as follows. First, since the interaction between the fluid inclusions and the surrounding solid in our current case is represented via not the classical boundary condition but a modified boundary condition, the average strain and stress in an isolated fluid inclusion may be defined equivalently from the following contour integrals as [29]

$${\epsilon }_{\mathrm{xx}}^{\mathrm{in}}=\frac{{\int }_L(u_x{n}_x)\mathrm{ds}}{V},{\epsilon }_{\mathrm{yy}}^{\mathrm{in}}=\frac{{\int }_L(u_y{n}_y)\mathrm{ds}}{V},{\epsilon }_{\mathrm{xy}}^{\mathrm{in}}=\frac{{\int }_L(u_y{n}_x+u_x{n}_y)\mathrm{ds}}{2V},$$

(47)

$${\sigma }_{\mathrm{xx}}^{\mathrm{in}}=\frac{{\int }_L(-p_x)\mathrm{xds}}{V},{\sigma }_{\mathrm{yy}}^{\mathrm{in}}=\frac{{\int }_L(-p_y)\mathrm{yds}}{V},{\sigma }_{\mathrm{xy}}^{\mathrm{in}}=\frac{{\int }_L(-p_x)\mathrm{yds}}{V}=\frac{{\int }_L(-p_y)\mathrm{xds}}{V},$$

(48)

where $n_x$ and $n_y$ are the Cartesian components of unit outward-pointing normal vector of the curve L, and $p_x$ and $p_y$ are the Cartesian components of the interfacial traction imposed on the solid matrix. Then, one may represent equations (47) and (48), in terms of the analytic functions $\varphi (z)$ and $\psi (z)$ , as [28]

$${\epsilon }_{\mathrm{xx}}^{\mathrm{in}}=\frac{\mathrm{Im}(C_1+C_2)}{2V},{\epsilon }_{\mathrm{yy}}^{\mathrm{in}}=\frac{\mathrm{Im}(C_2-C_1)}{2V},{\epsilon }_{\mathrm{xy}}^{\mathrm{in}}=-\frac{\mathrm{Re}(C_1)}{2V},$$

(49)

$${\sigma }_{\mathrm{xx}}^{\mathrm{in}}=\frac{\mathrm{Im}(D_1+D_2)}{2V},{\sigma }_{\mathrm{yy}}^{\mathrm{in}}=\frac{\mathrm{Im}(D_2-D_1)}{2V},{\sigma }_{\mathrm{xy}}^{\mathrm{in}}=-\frac{\mathrm{Re}(D_1)}{2V},$$

(50)

where

$$\begin{array}{c}\multicolumn{1}{c}{C_1={\int }_L\frac{\kappa \varphi (t)-t\overline{{\varphi }^{\prime }(t)}-\overline{\psi (t)}}{2G}\mathrm{dt},C_2={\int }_L\frac{\kappa \overline{\varphi (t)}-\overline{t}{\varphi }^{\prime }(t)-\psi (t)}{2G}\mathrm{dt},}\\ \multicolumn{1}{c}{D_1={\int }_{L}t·d(\varphi (t)+t\overline{{\varphi }^{\prime }(t)}+\overline{\psi (t)}),D_2={\int }_L\overline{t}·d(\varphi (t)+t\overline{{\varphi }^{\prime }(t)}+\overline{\psi (t)}),}\\ \multicolumn{1}{c}{V=\frac{1}{2}\mathrm{Im}\left({\int }_L\overline{t}\mathrm{dt}\right).}\end{array}$$

(51)

Furthermore, substituting the analytic solutions for $\varphi (z)$ and $\psi (z)$ derived in section 3 into equation (51), we identify

$$\begin{array}{c}\multicolumn{1}{c}{C_1=\pi i\frac{[\kappa {\alpha }_{-1}-e^{2i\theta }{K}_{-1}-\overline{{\mathrm{\Gamma }}_2}R{e}^{-i\theta }-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}R{e}^{i\theta }{a}_1]R{e}^{i\theta }}{G},}\\ \multicolumn{1}{c}{C_2=\pi i\frac{\begin{array}{c}-[\kappa \overline{{\alpha }_{-1}}+(\kappa -1){\mathrm{\Gamma }}_{1}R{e}^{-i\theta }\overline{a_1}-e^{-2i\theta }\overline{K_{-1}}-{\mathrm{\Gamma }}_{2}R{e}^{i\theta }]R{e}^{i\theta }{a}_1\hfill \\ \multicolumn{1}{c}{-\sum _{m=2}^{n-2}[\kappa \overline{{\alpha }_{-m}}+(\kappa -1){\mathrm{\Gamma }}_{1}R{e}^{-i\theta }\overline{a_m}-e^{-2i\theta }\overline{K_{-m}}]R{e}^{i\theta }m{a}_m}\\ \multicolumn{1}{c}{-[\kappa \overline{{\alpha }_{-n+1}}+(\kappa -1){\mathrm{\Gamma }}_{1}R{e}^{-i\theta }\overline{a_{n-1}}]R{e}^{i\theta }(n-1)a_{n-1}}\\ \multicolumn{1}{c}{-[\kappa \overline{{\alpha }_{-n}}+(\kappa -1){\mathrm{\Gamma }}_{1}R{e}^{-i\theta }\overline{a_n}]R{e}^{i\theta }n{a}_n+[(\kappa -1){\mathrm{\Gamma }}_1+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}]R^2}\end{array}}{G},}\\ \multicolumn{1}{c}{D_1=-2\pi \mathrm{iR}{e}^{i\theta }\left[\begin{array}{c}(2{\mathrm{\Gamma }}_{1}R{e}^{i\theta }{a}_1+{\alpha }_{-1}+\overline{{\mathrm{\Gamma }}_2}R{e}^{-i\theta }+e^{2i\theta }{K}_{-1})\hfill \\ \multicolumn{1}{c}{-a_{1}R{e}^{i\theta }[2{\mathrm{\Gamma }}_1-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}]}\end{array}\right]}\\ \multicolumn{1}{c}{D_2=-2\pi \mathrm{iR}{e}^{-i\theta }\left[\begin{array}{c}\overline{a_1}(2{\mathrm{\Gamma }}_{1}R{e}^{i\theta }{a}_1+{\alpha }_{-1}+\overline{{\mathrm{\Gamma }}_2}R{e}^{-i\theta }+e^{2i\theta }{K}_{-1})\hfill \\ \multicolumn{1}{c}{+\sum _{m=2}^{n-2}(2{\mathrm{\Gamma }}_{1}R{e}^{i\theta }{a}_m+{\alpha }_{-m}+e^{2i\theta }{K}_{-m})\overline{a_m}m}\\ \multicolumn{1}{c}{+\overline{a_{n-1}}(n-1)(2{\mathrm{\Gamma }}_{1}R{e}^{i\theta }{a}_{n-1}+{\alpha }_{-n+1})}\\ \multicolumn{1}{c}{+\overline{a_n}n(2{\mathrm{\Gamma }}_{1}R{e}^{i\theta }{a}_n+{\alpha }_{-n})-[2{\mathrm{\Gamma }}_1-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}]R{e}^{i\theta }}\end{array}\right],}\\ \multicolumn{1}{c}{V=\pi R^2(1-\sum _{m=1}^{n}m{a}_m\overline{a_m}),}\end{array}$$

(52)

where $\mathrm{\Delta }P$ is introduced from equation (32). Noting equation (12) one immediately finds that $C_1$ , $C_2$ , $D_1$ , and $D_2$ identified in equation (52) rely linearly on the remote loadings ${\sigma }_{\mathrm{xx}}^{\infty }$ , ${\sigma }_{\mathrm{xy}}^{\infty }$ , and ${\sigma }_{\mathrm{yy}}^{\infty }$ . Consequently, it is easy to extract the concentration coefficient matrices $A_c^d$ and $B_c^d$ from the expressions in equation (52), although the lengthy details are simply not recorded here.

Here, since the large number of fluid inclusions in the dilute or Mori–Tanaka method-based homogenization scheme are supposed to have a uniform distribution and identical orientation, the effective properties of the composite material attained would be generally anisotropic (except for the cases in which the fluid inclusions have highly symmetric shapes). However, the realistic scenario for fluid inclusion-filled composites is that the large number of fluid inclusions are oriented randomly in the surrounding solid leading to that the composite structure is overall isotropic. To model such a realistic scenario of isotropic effective properties, one may use the inclination angle-average concentration coefficient matrices $A_c^{d*}$ and $B_c^{d*}$ , instead of the inclination angle-dependent counterparts $A_c^d$ and $B_c^d$ , in corresponding calculations. Specifically, $A_c^{d*}$ and $B_c^{d*}$ are defined by

$$A_c^{d*}=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}{A}_c^d(\theta )d\theta ,B_c^{d*}=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}{B}_c^d(\theta )d\theta ,$$

(53)

where $\theta$ is the inclination angle introduced from equation (16). Substituting $A_c^{d*}$ and $B_c^{d*}$ into equation (46) (replacing $A_c^d$ and $B_c^d$ there), one is able to obtain the angle-average concentration coefficient matrices $A_c^{\mathrm{MT}*}$ and $B_c^{\mathrm{MT}*}$ for applications of the Mori–Tanaka method. Finally, by substituting these inclination angle-average concentration coefficient matrices into equation (43), one may anticipate isotropic effective properties for practical fluid inclusion-filled composite materials with random distributions in both space and orientation. The transverse compliance matrix of the homogenized isotropic composite material is denoted by $H_{\mathrm{eff}}^{*}$ (more specifically it is represented by $H_{\mathrm{eff}}^{*\mathrm{dilute}}$ and $H_{\mathrm{eff}}^{*\mathrm{MT}}$ , respectively, in the dilute and Mori–Tanaka methods), and the effective transverse Young’s modulus and shear modulus are defined by $E_{\mathrm{eff}}=1/H_{\mathrm{eff}}^{*}(1,1)$ and $G_{\mathrm{eff}}=1/(2H_{\mathrm{eff}}^{*}(3,3))$ , respectively.

4.1. Hypotrochoidal and elliptical fluid inclusions

In what follows, we present some specific explicit results for the cases of hypotrochoidal fluid inclusions and elliptical fluid inclusions. For a hypotrochoidal inclusion defined in equation (33), one may simplify equation (52) into

$$\begin{array}{c}\multicolumn{1}{c}{C_1=-\frac{(\kappa E+1)\overline{{\mathrm{\Gamma }}_2}\pi i{R}^2}{G(1-E^2(n-2)a_n^2)},}\\ \multicolumn{1}{c}{C_2=\pi i{R}^2\frac{(\kappa -1){\mathrm{\Gamma }}_1+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}+[\kappa \frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}-(\kappa -1){\mathrm{\Gamma }}_1]{\mathrm{na}}_n^2}{G},}\\ \multicolumn{1}{c}{D_1=\frac{-2\pi i{R}^2(1-E)\overline{{\mathrm{\Gamma }}_2}}{(1-E^2(n-2)a_n^2)},}\\ \multicolumn{1}{c}{D_2=2\pi i{R}^2\{[2{\mathrm{\Gamma }}_1-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}]+a_n^{2}n[-2{\mathrm{\Gamma }}_1+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}]\},}\\ \multicolumn{1}{c}{V=\pi R^2(1-a_n^{2}n),}\end{array}$$

(54)

where $\mathrm{\Delta }P$ is introduced from equation (37). Due to the high symmetry of hypotrochoidal inclusions, the corresponding concentration coefficient matrices ( $A_c^d$ and $B_c^d$ ) are found to be independent of the inclination angle $\theta$ and therefore they are equivalent to the angle-average counterparts ( $A_c^{d*}$ and $B_c^{d*}$ ). Particularly, we identify simple explicit closed-form results for the effective shear modulus (in the case of randomly oriented hypotrochoidal fluid inclusions) as

$$\begin{array}{c}\multicolumn{1}{c}{\frac{G_{\mathrm{eff}}^{\mathrm{dilute}}}{G}=\frac{(B^2(n-2)a_n^2-A^2)(a_n^{2}n-1)}{[(n^2-2n)a_n^4+(2-n)a_n^2]B^2+\mathrm{AB}{f}_{\mathrm{in}}(\kappa +1)-n{A}^2{a}_n^2+A^2},}\end{array}$$

(55)

$$\frac{G_{\mathrm{eff}}^{\mathrm{MT}}}{G}=\frac{(n^2-2n)(f_{\mathrm{in}}-1)B^2{a}_n^4+[(n-2)B^2+A^{2}n](1-f_{\mathrm{in}})a_n^2+f_{\mathrm{in}}\mathrm{AB}-A^2}{(n^2-2n)(f_{\mathrm{in}}-1)B^2{a}_n^4+[(n-2)B^2+A^{2}n](1-f_{\mathrm{in}})a_n^2-f_{\mathrm{in}}\kappa \mathrm{AB}-A^2}.$$

(56)

For an elliptical fluid inclusion, equation (52) reduces to

$$\begin{array}{c}\multicolumn{1}{c}{C_1=-\pi i{R}^2\frac{[\kappa [E\overline{{\mathrm{\Gamma }}_2}+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}{e}^{2i\theta }{a}_1]+\overline{{\mathrm{\Gamma }}_2}+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}{e}^{2i\theta }{a}_1]}{G},}\\ \multicolumn{1}{c}{C_2=\pi i{R}^2\frac{\begin{array}{c}[\kappa E{\mathrm{\Gamma }}_2{e}^{2i\theta }+\frac{\kappa (D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}{a}_1-(\kappa -1){\mathrm{\Gamma }}_1{a}_1+{\mathrm{\Gamma }}_2{e}^{2i\theta }]a_1\hfill \\ \multicolumn{1}{c}{+[(\kappa -1){\mathrm{\Gamma }}_1+\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}]}\end{array}}{G},}\\ \multicolumn{1}{c}{D_1=-2\pi i{R}^2[(1-E)\overline{{\mathrm{\Gamma }}_2}-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}{e}^{2i\theta }{a}_1+a_1{e}^{2i\theta }\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}],}\\ \multicolumn{1}{c}{D_2=-2\pi i{R}^2\left[\begin{array}{c}{a}_1(2{\mathrm{\Gamma }}_1{a}_1+(1-E)\overline{{\mathrm{\Gamma }}_2}{e}^{-2i\theta }-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{A}{a}_1)\hfill \\ \multicolumn{1}{c}{-[2{\mathrm{\Gamma }}_1-\frac{(D{\mathrm{\Gamma }}_1+\mathrm{\Delta }P)}{B}]}\end{array}\right],}\\ \multicolumn{1}{c}{V=\pi \mathrm{ab},}\end{array}$$

(57)

where $\mathrm{\Delta }P$ is introduced from equation (42). Similarly, in terms of the effective shear modulus, we attain also simple explicit closed-form result (in the case of randomly oriented elliptical fluid inclusions) as

$$\frac{G_{\mathrm{eff}}^{\mathrm{dilute}}}{G}=\frac{A({\mathrm{ABGa}}_1^2-\kappa a_1^2\mathrm{BK}-\mathrm{GAB}-\mathrm{AK})(a_1^2-1)}{g_1^{\mathrm{ed}}{A}^2+g_2^{\mathrm{ed}}A-(\kappa +1)0.5\kappa f_{\mathrm{in}}{\mathrm{Ka}}_1^2{B}^2},$$

(58)

$$\frac{G_{\mathrm{eff}}^{\mathrm{MT}}}{G}=\frac{g_1^{\mathrm{eM}}+g_2^{\mathrm{eM}}+0.5B^2{f}_{\mathrm{in}}K\kappa a_1^2}{g_1^{\mathrm{eM}}+g_3^{\mathrm{eM}}-0.5B^2{f}_{\mathrm{in}}K{\kappa }^2{a}_1^2},$$

(59)

where

$$\begin{array}{c}\multicolumn{1}{c}{g_1^{\mathrm{ed}}=(-a_1^4+2a_1^2-1)\mathrm{GB}+{\mathrm{Ka}}_1^2-K,}\\ \multicolumn{1}{c}{g_2^{\mathrm{ed}}=(a_1^2-1)(\kappa +1)f_{\mathrm{in}}G{B}^2+[\kappa a_1^4+((\kappa +1)0.5f_{\mathrm{in}}-\kappa )a_1^2-(\kappa +1)f_{\mathrm{in}}]\mathrm{KB},}\\ \multicolumn{1}{c}{g_1^{\mathrm{eM}}=([(f_{\mathrm{in}}-1)a_1^4+(2-f_{\mathrm{in}})a_1^2-1]\mathrm{GB}+(1-0.5f_{\mathrm{in}}){\mathrm{Ka}}_1^2-K)A^2,}\\ \multicolumn{1}{c}{g_2^{\mathrm{eM}}=((1-a_1^2)f_{\mathrm{in}}G{B}^2+[(1-f_{\mathrm{in}})\kappa {\mathrm{Ka}}_1^4+[0.5f_{\mathrm{in}}(\kappa -1)-\kappa ]{\mathrm{Ka}}_1^2+f_{\mathrm{in}}K]B)A,}\\ \multicolumn{1}{c}{g_3^{\mathrm{eM}}=((a_1^2-1)f_{\mathrm{in}}\kappa G{B}^2+[(1-f_{\mathrm{in}})\kappa {\mathrm{Ka}}_1^4+(f_{\mathrm{in}}-1)\kappa {\mathrm{Ka}}_1^2-f_{\mathrm{in}}K\kappa ]B)A.}\end{array}$$

(60)

5. Numerical examples

We will illustrate our analytical solutions via a few numerical examples of hypotrochoidal liquid inclusions, regular polygonal liquid inclusions, rectangular liquid inclusions, and elliptical liquid inclusions for plane strain deformations. Specifically, we show the incremental circumferential stress distribution around the inclusions, as well as the effective moduli of the composite materials containing randomly distributed and oriented inclusions (in terms of the Mori–Tanaka method).

We specify in the numerical calculations the material constants for the elastic matrix and liquid inclusion as follows:

$$G=1\mathrm{MPa},\nu =0.47,K=2\mathrm{GPa}.$$

(61)

Specially, in the numerical illustration of the incremental circumferential stress, we simply take the remote loading as a uniaxial tension/compression ${\sigma }_{\mathrm{yy}}^{\infty }$ in the y-direction and set the inclination angle $\theta$ of the inclusion to zero.

5.1. Verification

Figure 2 presents the incremental hoop stress, induced by the uniaxial loading (in the y-direction), around hypotrochoidal liquid inclusions for different initial pressure of the inclusions. It is seen from the figure that the presence of the initial pressure consistently reduces the incremental hoop stress concentration around the inclusions, and particularly when the initial pressure vanishes the hoop stress predicted by the current solutions becomes consistent with that determined from the classical solutions for an incompressible liquid inclusion derived in the work by Wang and Schiavone [21]. Here, since the specific liquid inclusion considered here has an extremely large bulk modulus (relative to the shear modulus of the matrix), it undoubtedly may be treated as being incompressible explaining the consistency between the classical solutions in the work by Wang and Schiavone [21] and the current solutions when the initial pressure vanishes.

Figure 2.

Incremental hoop stress around a hypotrochoidal liquid inclusion related to a uniaxial remote loading $(\theta =0)$ .

5.2. Regular polygonal liquid inclusions

The mapping function for a p-sided regular polygonal inclusion may be represented by Schwarz–Christoffel as follows:

$$w(\xi )=R{e}^{i\theta }\underset{1}{\overset{\xi }{\int }}{[1+{(-1)}^p\frac{1}{t^p}]}^{\frac{2}{p}}\mathrm{dt}.$$

(62)

By expanding equation (62) into series in $\xi$ , one may obtain the standard polynomial mapping function for the inclusion (as introduced in equation (16)). For example, the polynomial mapping functions for equilateral triangular, square, regular pentagonal, and regular hexagonal inclusions appear, respectively, to be

$$w(\xi )=R{e}^{i\theta }(\xi +\frac{1}{3}{\xi }^{-2}+\frac{1}{45}{\xi }^{-5}+\frac{1}{162}{\xi }^{-8}+\cdots +a_n{\xi }^{-n}),\mathrm{for}\mathrm{equilateral}\mathrm{triangle},$$

(63)

$$w(\xi )=R{e}^{i\theta }(\xi -\frac{1}{6}{\xi }^{-3}+\frac{1}{56}{\xi }^{-7}-\frac{1}{176}{\xi }^{-11}+\cdots +a_n{\xi }^{-n}),\mathrm{for}\mathrm{square},$$

(64)

$$w(\xi )=R{e}^{i\theta }(\xi +\frac{1}{10}{\xi }^{-4}+\frac{1}{75}{\xi }^{-9}+\frac{4}{875}{\xi }^{-14}+\cdots +a_n{\xi }^{-n}),\mathrm{for}\mathrm{regular}\mathrm{pentagon},$$

(65)

$$w(\xi )=R{e}^{i\theta }(\xi -\frac{1}{15}{\xi }^{-5}+\frac{1}{99}{\xi }^{-11}-\frac{5}{1377}{\xi }^{-17}+\cdots +a_n{\xi }^{-n}),\mathrm{for}\mathrm{regular}\mathrm{hexagon}.$$

(66)

Since the current analytic solutions established in section 3 only allow for a finite-term polynomial mapping function, one may truncate the infinite series of a regular polygonal inclusion by a sufficiently large number n for practical calculations. For example, as shown in Figure 3, the first four terms in the mapping function are sufficient to effectively restore the shapes of some common regular polygons (from a practical engineering point of view).

Figure 3.

Approximately regular polygonal liquid inclusion described by finite-term mapping functions $(\theta =0)$ .

We first calculate the uniaxial remote loading-induced hoop stress distribution around the inclusion for the common regular polygonal inclusion shapes defined in equations (63)–(66) with $\theta =0$ , as shown in Figure 4. We see from Figure 4 that the presence of initial pressure always reduces the circumferential stress around the inclusion, and as the initial pressure increases, the circumferential stress around the inclusion gradually decreases. In addition, we note that as the initial pressure increases, the rate of decrease in the circumferential stress around the inclusion varies, particularly at the vertex, depending on the vertex angle and the direction of the external load relative to the two sides forming the vertex (e.g., typically, if one of the two sides forming the vertex is parallel to the external load, as shown in Figure 4(b) for the square case, the rate of decrease in the circumferential stress near the vertex is relatively small).

Figure 4.

Hoop stress around a regular polygonal liquid inclusion for a uniaxial remote loading $(\theta =0)$ .

We continue to analyze the transverse effective moduli of the composite material containing randomly distributed and randomly oriented cylindrical regular polygonal liquid inclusions. Due to the random distribution and orientation of the inclusions, the composite material exhibits an overall transverse isotropy. The effective shear modulus and Young’s modulus of the composite material relative to the volume fraction of the inclusions are illustrated in Figures 5–8 for practical regular polygonal inclusion shapes shown in Figure 3. We see from these figures that the introduction of macroscale liquid inclusions (without small-scale surface effects) into a solid matrix always weakens the stiffness of the matrix, and for a given volume fraction of inclusions, both the effective Young’s modulus and shear modulus increase when the number of sides of the inclusions increases (although the rate of increase in the effective moduli with increasing number of sides of the inclusions reduces rapidly, see Figures 5 and 6), reflecting that for a given volume fraction of liquid inclusions circular inclusion shape is an optimal shape among all regular polygonal shapes which minimizes such weakening of the matrix. On the other hand, as illustrated in Figures 7 and 8, the presence of initial pressure always enhances the overall tensile/compressive and shear properties of the composite material (for given shape of the liquid inclusions), although such enhancement effects are weakened as the number of sides of the liquid inclusions increase.

Figure 5.

Effective transverse Young’s modulus of a composite material containing randomly distributed and oriented cylindrical regular polygonal liquid inclusions with or without initial pressure relative to the number of sides of the inclusions.

Figure 6.

Effective transverse shear modulus of a composite material containing randomly distributed and oriented cylindrical regular polygonal liquid inclusions with or without initial pressure relative to the number of sides of the inclusions.

Figure 7.

Initial pressure effects on the effective Young’s modulus of a composite material containing randomly distributed and oriented cylindrical regular polygonal liquid inclusions.

Figure 8.

Initial pressure effects on the effective transverse shear modulus of a composite material containing randomly distributed and oriented cylindrical regular polygonal liquid inclusions.

5.3. Slender liquid inclusions

Small-scale liquid inclusions in a soft elastic matrix appear commonly in plump configurations due to surface tension effects. For the currently considered macroscale liquid inclusions, however, they often appear, in terms of slender shapes, for both soft and stiff elastic matrices. We will present some results for the effective moduli of composite materials containing randomly distributed and oriented slender liquid inclusions which are modeled using either elliptical or rectangular shapes.

The mapping function for an elliptical inclusion shape is defined in equation (39), while those for approximately rectangular ones are usually expressed via [36]

$$z=w(\xi )=R{e}^{i\theta }(\xi +\sum _{m=1,3,5,...}{a}_m{\xi }^{-m}).$$

(67)

Here, we are mainly concerned with approximately rectangular liquid inclusions for a few typical aspect ratios, whose specific shapes and corresponding coefficients of the mapping functions are illustrated in Figure 9.

Figure 9.

Approximately rectangular liquid inclusion described by finite-term mapping functions $(\theta =0)$ .

The transverse effective shear modulus and Young’s modulus of the composite material containing randomly distributed and oriented slender cylindrical liquid inclusions (modeled as either rectangular or elliptical shapes) are illustrated in Figures 10 and 11 for different aspect ratios. We see from these figures that for identical aspect ratio and identical volume fraction of the liquid inclusions, compared to the model of elliptical liquid inclusions, that of rectangular liquid inclusions predicts overall smaller transverse effective moduli of the composite material (except for the cases corresponding to very slender liquid inclusions, it would predict slightly larger effective shear modulus). Overall speaking, as demonstrated from Figures 10 and 11, for composite materials with slender liquid inclusions whose aspect ratio exceeds a certain value (e.g., 5), the use of either elliptical or rectangular inclusion shapes in modeling the effective properties is acceptable (the differences between them are negligible). On the other hand, we see from the specific data for relative deviation in Figures 10 and 11 that the presence of initial pressure in the inclusions always enhances the overall tensile/compressive and shear properties of the composite material, and such enhancement effects would be amplified with increasing aspect ratio of the inclusions (as the inclusions become slenderer).

Figure 10.

Effective Young’s modulus of a composite material containing randomly distributed and oriented slender liquid inclusions.

Figure 11.

Effective shear modulus of a composite material containing randomly distributed and oriented slender liquid inclusions.

6. Conclusion

We reinvestigate, in the context of linear elasticity, the plane deformation of a macroscale compressible fluid inclusion of general shape in an elastic solid. In contrast to existing theoretical study of this problem in the literature, we focus on the influence of the pre-existence of the initial pressure in the fluid inclusion on the elastic response of the solid when a constant far-field loading is applied. To this end, we formulate, in the context of the complex variable formalism for plane elasticity, the boundary value problem corresponding to the use of a modified boundary condition accounting for the directional change in the initial pressure-induced traction acting on the solid (during the deformation caused by the far-field loading), and then derive initial pressure-dependent explicit analytic solutions for full incremental elastic field in the solid induced by the far-field loading. Further merging the obtained analytic solutions with some homogenization methods, we attain initial pressure-dependent effective in-plane properties of the solid-based composite containing a group of randomly distributed fluid inclusions occupying a given volume fraction. Numerical examples are given for a soft elastic solid containing, respectively, hypotrochoidal liquid inclusions (for comparative validation only), (approximately) regular polygonal liquid inclusions, (approximately) rectangular liquid inclusions, and elliptical liquid inclusions. The main findings are summarized as follows:

1. The presence of initial pressure in the inclusion reduces the circumferential stress concentration (induced by the far-field loading) around the inclusion and enhances the overall elastic resistance (both tensile/compressive and shear) of the corresponding composite. Such enhancement effects on the overall elastic resistance are weakened as the number of sides of the inclusion increases (in terms of regular polygonal inclusions), while they are amplified as the inclusion becomes slenderer (in terms of rectangular and elliptical inclusions).

2. Slender fluid inclusions may be modeled as either of elliptical and rectangular shapes in predicting the overall resistance of the corresponding fluid–solid composites as long as the aspect ratio of the slender inclusions is around or larger than five.