This section recalls some results about the calculation of the Hill polarization tensors related to a matrix of stiffness \(\mathbb{C}\) and an ellipsoid \(\mathcal{E}_{\uu{A}}\) of equation \[ \uv{x}\in\mathcal{E}_{\uu{A}} \quad\Leftrightarrow\quad \uv{x}\cdot(\trans{\uu{A}}\cdot\uu{A})^{-1}\cdot\uv{x}\leq 1 \] where \(\uu{A}\) is an invertible second-order tensor so that \(\trans{\uu{A}}\cdot\uu{A}\) is a positive definite symmetric tensor associated to 3 radii (eigenvalues \(a\geq b \geq c\) possibly written \(\rho_1 \geq \rho_2 \geq \rho_3\) for convenience) and 3 angles (orientation of the frame of eigenvectors \(\uv{e}_1, \uv{e}_2, \uv{e}_3\)) \[ \trans{\uu{A}}\cdot\uu{A}=a^2 \uv{e}_1\otimes\uv{e}_1 + b^2 \uv{e}_2\otimes\uv{e}_2 + c^2 \uv{e}_3\otimes\uv{e}_3 = \sum_{i=1}^3 \rho_i \uv{e}_i\otimes\uv{e}_i \tag{B.1}\]
A general expression of the elastic polarization tensor is derived in (Willis, 1977) (see also (Mura, 1987)) \[\begin{aligned} \uuuu{P}(\uu{A},\uuuu{C})&=\frac{1}{4\pi} \int_{\norm{\uv{\zeta}}=1} (\uu{A}^{-1}\cdot\uv{\zeta})\sotimes \Big((\uu{A}^{-1}\cdot\uv{\zeta})\cdot\uuuu{C} \cdot(\uu{A}^{-1}\cdot\uv{\zeta})\Big)^{-1} \sotimes(\uu{A}^{-1}\cdot\uv{\zeta}) \ud S_\zeta\\ &= \frac{\det{\uu{A}}}{4\pi} \int_{\norm{\uv{\xi}}=1} \frac{\uv{\xi}\sotimes (\uv{\xi}\cdot\uuuu{C} \cdot\uv{\xi})^{-1} \sotimes\uv{\xi}}{\norm{\uu{A}\cdot\uv{\xi}}^3} \ud S_{\xi} \end{aligned} \tag{B.2}\]
When \(\uuuu{C}\) is arbitrarily anisotropic, it is necessary to resort to numerical cubature to estimate \(\uuuu{P}\) as proposed in (Ghahremani, 1977), (Gavazzi and Lagoudas, 1990) or (Masson, 2008). However in some cases of anisotropy, analytical solutions are available ((Withers, 1989), (Barthélémy, 2020)). The case of isotropic matrix is particularly developed in the next section.
In this section, the matrix is assumed isotropic so that its stiffness tensor writes by means of a bulk \(k\) and shear \(\mu\) or Lamé \(\lambda\) and \(\mu\) moduli or even Young modulus \(E\) and Poisson ratio \(\nu\) with \(k=\frac{E}{3(1-2\nu)}\) and \(\mu=\frac{E}{2(1+\nu)}\). \[\begin{aligned} \uuuu{C} =& {} 3k\uuuu{J}+2\mu\uuuu{K} = 3\lambda\uuuu{I}+2\mu\uuuu{K}\\ & \quad\textrm{with}\quad J_{ijkl}=\frac{\delta_{ij}\delta_{kl}}{3}, I_{ijkl}=\frac{\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}}{2} \textrm{ and } \uuuu{K}=\uuuu{I}-\uuuu{J} \end{aligned} \tag{B.3}\]
Introducing (B.3) in (B.2) leads to after some algebra \[ \uuuu{P}= \frac{1}{\lambda+2\,\mu} \uuuu{U} +\frac{1}{\mu}(\uuuu{V}-\uuuu{U}) \] where the tensors \(\uuuu{U}\) and \(\uuuu{V}\), depending only on the ellipsoidal tensor \(\uu{A}\) of (B.1), are given by (see (Barthélémy, 2020)) \[\begin{aligned} \uuuu{U} &= \frac{\det{\uu{A}}}{4\pi} \int_{\norm{\uv{\xi}}=1} \frac{\uv{\xi}\otimes\uv{\xi}\otimes\uv{\xi}\otimes\uv{\xi}} {\norm{\uu{A}\cdot\uv{\xi}}^3}\ud S_{\xi}\\ &= \frac{1}{4\pi} \int_{\norm{\uv{\zeta}}=1} \frac{% (\uu{A}^{-1}\cdot\uv{\zeta}) \otimes (\uu{A}^{-1}\cdot\uv{\zeta}) \otimes (\uu{A}^{-1}\cdot\uv{\zeta}) \otimes (\uu{A}^{-1}\cdot\uv{\zeta}) }{\norm{\uu{A}^{-1}\cdot\uv{\zeta}}^4} \ud S_{\zeta} \end{aligned}\] and \[\begin{aligned} \uuuu{V} &= \frac{\det{\uu{A}}}{4\pi} \int_{\norm{\uv{\xi}}=1} \frac{\uv{\xi}\sotimes\uu{1}\sotimes\uv{\xi}} {\norm{\uu{A}\cdot\uv{\xi}}^3}\ud S_{\xi}\\ &= \frac{1}{4\pi} \int_{\norm{\uv{\zeta}}=1} \frac{% (\uu{A}^{-1}\cdot\uv{\zeta}) \sotimes \uu{1} \sotimes (\uu{A}^{-1}\cdot\uv{\zeta}) }{\norm{\uu{A}^{-1}\cdot\uv{\zeta}}^2} \ud S_{\zeta} \end{aligned}\] For an arbitrary ellipsoid defined by (B.1), the components of \(\uuuu{U}\) and \(\uuuu{V}\) write \[\begin{aligned} U_{iiii}&=\frac{3(I_i-\rho_i^2I_{ii})}{2} \quad\forall\, i\in\{1,2,3\}\\ U_{iijj}=U_{ijij}=U_{ijji}&=\frac{I_j-\rho_i^2I_{ij}}{2} =\frac{I_i-\rho_j^2I_{ij}}{2} \quad\forall\, i\neq j\in\{1,2,3\} \end{aligned}\] and \[\begin{aligned} V_{iiii}&=I_i\quad\forall\, i\in\{1,2,3\}\\ V_{ijij}=V_{ijji}&=\frac{I_i+I_j}{4} \quad\forall\, i\neq j\in\{1,2,3\} \end{aligned}\] where the coefficients \(I_i\) and \(I_{ij}\) are given by (note that \(I_i\) and \(I_{ij}\) are adapted from those provided in (Kellogg, 1929) and (Eshelby, 1957): they differ by a factor of \(4\pi/3\) for \(I_{ij}\) with \(i\neq j\) and by \(4\pi\) for the others)
In this last case of spherical inclusion (\(\uu{A}=\uu{1}\)), \(\uuuu{U}\) and \(\uuuu{V}\) are simply decomposed as \[ \uuuu{U}=\frac{1}{3}\uuuu{J}+\frac{2}{15}\uuuu{K} \quad\textrm{ and }\quad \uuuu{V}=\frac{1}{3}\uuuu{I} \]
The case of cracks corresponds to ellipsoids for which the smallest radius is very small compared to the two others, in other words the characteristic tensor \(\uu{A}\) (B.1) can be written here \[ \uu{A}= \uv{\ell}\otimes\uv{\ell}+ \eta\,\uv{m}\otimes\uv{m}+ \omega\,\uv{n}\otimes\uv{n} \quad\textrm{with}\quad \eta=\frac{b}{a} \quad\textrm{and}\quad \omega=\frac{c}{a} \]
In the case of cracks, it is useful to introduce the second Hill polarization tensor defined as \[ \uuuu{Q}=\uuuu{C}-\uuuu{C}:\uuuu{P}:\uuuu{C} \] and in particular \(\lim_{\omega\to 0}\omega\,\uuuu{Q}^{-1}\) in which it is recalled that \(\uuuu{P}\) and thus \(\uuuu{Q}\) depend on \(\omega\) such that the components \(\left(\uuuu{Q}^{-1}\right)_{nijk}\) (with \(n\) corresponding to the crack normal) behave as \(1/\omega\) when \(\omega\) tends towards \(0\). The analytical expressions of this limit are fully detailed in (Barthélémy et al., 2021) which recalls in particular that \(\uuuu{L}\) actually derives from a symmetric second-order tensor \(\uu{B}\) as \[ \uuuu{L}= \lim_{\omega\to 0} \omega\,\uuuu{Q}^{-1} =\frac{3}{4}\,\uv{n}\sotimes\uu{B}\sotimes\uv{n} \tag{B.4}\]
For an arbitrarly anisotropic matrix, an algorithm allowing to estimate the limit (B.4) is proposed in (Barthélémy, 2009) whereas in the isotropic case \(\uu{B}\) writes \[ \uu{B}= B_{nn}\,\uv{n}\otimes\uv{n} + B_{mm}\,\uv{m}\otimes\uv{m} + B_{\ell\ell}\,\uv{\ell}\otimes\uv{\ell} \] with \[\begin{aligned} B_{nn}&=\frac{8\,\eta\,(1-\nu^2)}{3\,E}\, \frac{1}{\mathcal{E}_\eta}\label{eq:Bnn}\\ B_{mm}&=\frac{8\,\eta\,(1-\nu^2)}{3\,E}\, \frac{1-\eta^2}{\left(1-(1-\nu)\,\eta^2\right) \,\mathcal{E}_\eta-\nu\,\eta^2\,\mathcal{K}_\eta}\\ B_{\ell\ell}&=\frac{8\,\eta\,(1-\nu^2)}{3\,E}\, \frac{1-\eta^2}{(1-\nu-\eta^2)\,\mathcal{E}_\eta+\nu\,\eta^2\,\mathcal{K}_\eta} \end{aligned}\] where \(\mathcal{K}_\eta=\mathcal{K}(\sqrt{1-\eta^2})\) and \(\mathcal{E}_\eta=\mathcal{E}(\sqrt{1-\eta^2})\) are the complete elliptic integrals of respectively the first and second kind (see (Abramowitz and Stegun, 1972)). If the crack is circular, the components of \(\uu{B}\) become \[ B_{nn}=\frac{16\,(1-\nu^2)}{3\,\pi\,E} \quad\textrm{;}\quad B_{mm}=B_{\ell\ell}=\frac{B_{nn}}{1-\nu/2} \]
C = stiff_Enu(1.,0.2) ; print(C)Order 4 ISO tensor | Param(size=2)=[ 1.66667 0.833333 ] | Angles(size=0)=[ ]
[ 1.11111 0.277778 0.277778 0 0 0
0.277778 1.11111 0.277778 0 0 0
0.277778 0.277778 1.11111 0 0 0
0 0 0 0.833333 0 0
0 0 0 0 0.833333 0
0 0 0 0 0 0.833333 ]
Note that in Echoes it is necessary to provide an aspect ratio \(\omega\) for the crack even if the crack compliance is actually calculated as a limit (not depending on \(\omega\))
ω = 1.e-4
L = crack_compliance(spheroidal(ω), C) ; print(L)[[0. 0. 0. 0. 0. 0. ]
[0. 0. 0. 0. 0. 0. ]
[0. 0. 1.22230996 0. 0. 0. ]
[0. 0. 0. 0.67906109 0. 0. ]
[0. 0. 0. 0. 0.67906109 0. ]
[0. 0. 0. 0. 0. 0. ]]tω = np.logspace(-5,1,20)
tabδ = []
for ω in tω:
Q = hill_dual(spheroidal(ω), C)
Lω = ω*np.linalg.inv(Q)
δL = np.linalg.norm(Lω-L)/np.linalg.norm(L)
tabδ.append(δL)
plt.figure(figsize=(8,3))
plt.loglog(tω,tabδ,'+-')
plt.xlabel(r"$\omega$")
plt.ylabel(r"$\frac{||\mathbb{L}-\omega\,\mathbb{Q}^{-1}||}{||\mathbb{L}||}$")
plt.grid(True,which='both')
plt.show()