import numpy as np
from echoes import *
import math, random
import matplotlib.pyplot as plt
np.set_printoptions(precision=8, suppress=True)
# to display only 8 significant digits of array components6 Crack compliance
Definition of the crack compliance
Geometrically an elliptical crack can be defined as the limit domain of an ellipsoid for which the smallest aspect ratio tends towards 0. Keeping the definition of the ellipsoid by means of the tensor \(\uu{A}\) in (5.1), the latter tensor becomes 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 recall the second Hill polarization tensor defined as (see (5.6)) \[ \uuuu{Q}=\uuuu{C}-\uuuu{C}:\uuuu{P}:\uuuu{C} \] and in particular the so-called crack compliance (Kachanov, 1992), (Kachanov, 1993), (Sevostianov and Kachanov, 2002), (Barthรฉlรฉmy et al., 2021) \[ \uuuu{H}=\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\). It is also shown in previously cited references that \(\uuuu{H}\) actually derives from a symmetric second-order tensor \(\uu{B}\) as \[ \uuuu{H}= \lim_{\omega\to 0} \omega\,\uuuu{Q}^{-1} =\frac{3}{4}\,\uv{n}\sotimes\uu{B}\sotimes\uv{n} \tag{6.1}\] For more details about this limit, the reader can also refer to (Laws, 1977), (Laws, 1985), (Nemat-Nasser and Hori, 1999), (Barthรฉlรฉmy, 2009), (Kachanov and Sevostianov, 2018), (Barthรฉlรฉmy et al., 2021).
For an arbitrarly anisotropic matrix, an algorithm allowing to estimate the limit (6.1) 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{align} 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{align}\] 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} \]
Numerical implementation
The syntax allowing to calculate the compliance \(\uuuu{H}\) of an elliptical crack embedded in an infinite matrix of elastic stiffness \(\uuuu{C}\) is very close to that corresponding to hill and eshelby in Chapter 5, in other words
H = crack_compliance(shape, C, algo=NUMINT, epsroots=1.e-4, epsabs=1.e-4, epsrel=1.e-4, maxnb=100000)
where the keyword argument algo can be either NUMINT or RESIDUES (default: NUMINT). See (Barthรฉlรฉmy, 2009) for more details.
Warning
As recalled above the elliptical crack is described as a flat ellipsoid. However the calculation of the crack compliance \(\uuuu{H}\) (6.1) relies on an ellipsoidal argument (see Section 5.2) in which even the smallest aspect ratio should be strictly positive for numerical reasons. It is then 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\)).
Isotropic example
C = stiff_Enu(1.,0.2) ; print("C =\n", C)
ฯ = 1.e-4
H = crack_compliance(spheroidal(ฯ), C) ; print("H =\n", H)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 ]
H =
[[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. ]]Anisotropic example
M = np.array([ [2.66011, 1.26432, 0.662772, 1.9402, 1.54905, 1.10384],
[1.26432, 3.6072, 1.78964, 2.0247, 1.2701, 1.1089],
[0.662772, 1.78964, 2.743, 1.3367, 1.2962, 0.897632],
[1.9402, 2.0247 ,1.3367, 4.42684, 2.05632, 1.52686],
[1.54905, 1.2701, 1.2962, 2.05632, 3.54431, 1.3445],
[1.10384, 1.1089, 0.897632, 1.52686, 1.3445, 1.99356] ])
C = tensor(M) ; print(C)
ฯ = 1.e-4
H = crack_compliance(spheroidal(ฯ), C, algo=NUMINT) ; print("H =\n", H)Order 4 ANISO tensor | Param(size=21)=[ 2.66011 1.26432 0.662772 1.9402 1.54905 1.10384 3.6072 1.78964 2.0247 1.2701 1.1089 2.743 1.3367 1.2962 0.897632 4.42684 2.05632 1.52686 3.54431 1.3445 1.99356 ] | Angles(size=0)=[ ]
[ 2.66011 1.26432 0.662772 1.9402 1.54905 1.10384
1.26432 3.6072 1.78964 2.0247 1.2701 1.1089
0.662772 1.78964 2.743 1.3367 1.2962 0.897632
1.9402 2.0247 1.3367 4.42684 2.05632 1.52686
1.54905 1.2701 1.2962 2.05632 3.54431 1.3445
1.10384 1.1089 0.897632 1.52686 1.3445 1.99356 ]
H =
[[ 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0.57167315 -0.08924584 -0.10580477 0. ]
[ 0. 0. -0.08924584 0.2711102 -0.06996412 0. ]
[ 0. 0. -0.10580477 -0.06996412 0.34412187 0. ]
[ 0. 0. 0. 0. 0. 0. ]]H = crack_compliance(spheroidal(ฯ), C, algo=RESIDUES) ; print("H =\n", H)H =
[[ 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0.57167315 -0.08924584 -0.10580477 0. ]
[ 0. 0. -0.08924584 0.2711102 -0.06996412 0. ]
[ 0. 0. -0.10580477 -0.06996412 0.34412187 0. ]
[ 0. 0. 0. 0. 0. 0. ]]This exercise aims at checking the aspect ratio for which the approximation \(\omega\,\uuuu{Q}^{-1}\approx \lim_{\omega\to 0}\omega\,\uuuu{Q}^{-1}\) remains acceptable. For this purpose it is proposed to build the graph of relative distance between these two tensors as function of the aspect ratio \(\omega\) for different reference stiffness tensors \(\uuuu{C}\).
Hints: with hill_dual use preferably algo=NUMINT and decrease the values of epsrel and espabs to avoid roundoff error prior to inversion of \(\uuuu{Q}\)
Solution
plt.figure(figsize=(8,3))
for i in range(4):
A = np.random.rand(6,6)
C = tensor(A.T.dot(A) + np.eye(6)) # generation of an arbitrary positive definite matrix
ฯ = 1.e-4
H = crack_compliance(spheroidal(ฯ), C)
tฯ = np.logspace(-5,1,20)
tabฮด = []
for ฯ in tฯ:
Q = hill_dual(spheroidal(ฯ), C, algo=NUMINT, epsrel=1.e-12, epsabs=1.e-12)
Hฯ = ฯ*np.linalg.inv(Q)
ฮดH = np.linalg.norm(Hฯ-H)/np.linalg.norm(H)
tabฮด.append(ฮดH)
plt.loglog(tฯ,tabฮด,'+-')
plt.xlabel(r"$\omega$")
plt.ylabel(r"$\frac{||\mathbb{H}-\omega\,\mathbb{Q}^{-1}||}{||\mathbb{H}||}$")
plt.grid(True,which='both')
plt.show()
\(\,\)