Degenerate shapes and extreme aspect ratios

Before this page

Hill polarization tensors for the closed forms and the cubatures, and COD tensors for the cracks.

The closed forms of the Hill, Eshelby and crack opening tensors are written for a generic shape. At the shapes where two semi-axes meet, where a spheroid becomes a sphere, a disc or a needle, or where the two acoustic roots of a transversely isotropic matrix coincide, they reach removable limits. Next to those shapes they are differences of nearly equal terms. Evaluated as written, they lose every digit there. Their derivatives fare worse: the derivative with respect to a semi-axis at the point where it meets another one, or with respect to a modulus at an isotropic point, is a limit as well, and a code that merges the two semi-axes or the two roots drops it.

This page lists those situations, says for each one what goes wrong, and how the package evaluates it instead. Every figure quoted is a relative error, measured against a reference that shares no code with the formula it checks: the nested cubature at tolerances , the same integrals in 256-bit arithmetic, and, for derivatives, ForwardDiff through the cubature and finite differences of the closed form.

Situation	What goes wrong	Treatment	Accuracy
spheroid near the sphere, isotropic matrix	divided difference of two potentials	the potentials as the integrals they stand for
triaxial ellipsoid with two close semi-axes	elliptic integrals and their divided differences	expansion around the spheroid
derivative where two semi-axes meet	the merged shape no longer depends on their gap	shapes merged only when their derivatives coincide
spheroid in a TI matrix: sphere, coincident roots, disc, needle	cancellations, a square root of zero	three regimes in the mean and squared half-difference of the roots
elliptical crack near the circle		power series
a cubature under ForwardDiff	NaN error estimate, dual integration bounds	error norm on the values, plain bounds	—
self-consistent estimate of TI class	estimate rebuilt as a full tensor	projection within the class	faster

The integrals behind the first four rows live in src/Core/spheroid_integrals.jl.

The integrals of a spheroid

With      and   , over   , three families of integrals carry the shape of a spheroid of axial over transverse semi-axis :

With   they give the Newton potentials of the spheroid; with  , a root of the acoustic polynomial of a transversely isotropic matrix, the Hill tensor of a coaxial spheroid in it. is elementary through  , which is below  :

With   the denominator is   , and these are the integrals of (Barthélémy, 2020; Appendix D) in the variables of that article:  ,  , and  , so that the three closed forms above are those of its Table D.1. The moments of that the package evaluates are differences of them,    and    . The application Hill tensors by transformation proves the closed forms and the partial fractions symbolically, and the correspondence between both writings.

These cancel as  , the sphere. Writing     with   , all three families are instead power series in there, with the coefficients       :

They are used for   , where 32 terms leave a truncation below (64 for the moments of high power ). Away from the sphere the moments follow from   and    ; follows from by two recurrences, the integral of and    :

Spheroid near the sphere, isotropic matrix

The closed forms in an isotropic matrix — Hill, Eshelby and concentration tensors, and , and conduction — are built on the Newton potentials of the ellipsoid. For a spheroid of axial semi-axis and transverse semi-axis , the mixed second-order potential is classically the divided difference     of two first-order ones that tend to the same value at the sphere. A Hill tensor built from it is off by at    and by more than its own size at    .

newton_potential_3d evaluates instead, with  ,

then   ,    and   , none of which cancels.

Where the integral forms come from

The potentials are        with  . The change of variable    maps   onto   with    and    . The axial potential becomes

and the mixed one, with one more factor   , the moment divided by .

Against the nested cubature at tight tolerances, the Hill tensor of a spheroid in an isotropic matrix is now within   at every distance from the sphere tried, from to .

Triaxial ellipsoid close to a spheroid

A triaxial ellipsoid goes through Legendre's elliptic integrals, and its mixed potentials through divided differences such as    . When two semi-axes come together the first cancel as   and the second as   : the error grows as with the relative gap , to   at  . Taking the two semi-axes as equal below some tolerance only trades that error for one of the order of the tolerance.

When the closer pair of semi-axes differ by less than 2 % (in squares), the potentials are expanded around the spheroid they are close to. With the third semi-axis,  ,   and the same change of variable as above, every potential takes the form

where   are the exponents of  ,  ,   in the integrand and     . With    and   ,    , and

so that each potential is a combination of the moments . The series converges as    ; eleven terms leave a truncation below at the 2 % threshold, beyond which the elliptic integrals are accurate to . The expansion is exact when  , derivatives included, which the elliptic integrals and their divided differences cannot be.

Against the cubature, a triaxial ellipsoid with two semi-axes to apart — next to a flat, an elongated or a near-spherical spheroid — is within , and within   on the elliptic side of the threshold. Two semi-axes of plain floating-point type within of each other are still taken as equal, which costs an error of that order.

Where two semi-axes meet: derivatives

An ellipsoid whose two semi-axes are equal is classified as a spheroid (or a sphere), and its closed form depends on the shared semi-axis only. Differentiate with respect to one of the two at the point where they meet, and that closed form has nothing to give: merged, the derivative of the Hill tensor of a sphere with respect to one of its semi-axes is zero, and the derivatives along the gap of the shapes next to a spheroid are off by 100 % to 550 %, without a word.

Two semi-axes carried by ForwardDiff.Dual numbers are therefore merged only when their derivatives coincide too; otherwise the ellipsoid is triaxial, and the expansion above, exact at  , carries the derivative along the gap. The same rule classifies ellipses, cylinders and cracks. A spheroid parametrized by its aspect ratio, , keeps its two equal semi-axes constant, so it stays a spheroid under differentiation. For conduction, the derivative of the Hill tensor of a sphere of radius 1 with respect to its third semi-axis is

which the package returns to rounding. Against ForwardDiff through the cubature these derivatives agree to  , and against central finite differences of the closed form to the accuracy of the differences, to  .

Spheroid in a transversely isotropic matrix

The Walpole coefficients of a spheroid coaxial with a transversely isotropic matrix, given in Hill polarization tensors, involve   of the two acoustic roots   and of  . The package evaluates them with

as

These are the coefficients of that page, in which appear in differences of terms or larger than the result for a needle; the combinations    and     are computed directly instead.

The are symmetric in the two roots: they depend on them only through the mean and the squared half-difference ,

rational in the moduli and free of their scale. Four situations would otherwise break the evaluation:

• near the sphere,  , the closed forms cancel;

• coincident roots — an isotropic matrix stored as transversely isotropic, or one within of isotropy — make the partial fractions divide by  , and the square root of has an infinite derivative; merging the roots there instead makes the derivatives along every modulus that breaks isotropy wrong by 24 % to a factor 14, at an isotropic point included;

• a flat spheroid,  , makes the partial fractions of cancel as ;

• a needle,  , makes those of and cancel as .

They are evaluated in three regimes:

1. near the sphere,   : the power series    , where        needs only     and    ;

2. nearly coincident roots,  : with  ,  , truncated at  , so that  ;

3. elsewhere, partial fractions,   with      , arranged so as not to cancel. For a flat spheroid ( ), and come from by    . For a needle,    turns the single-root integrals into a part rational in , summed exactly, and a remainder in :

   whose partial fractions cancel as only.

The sign of the partial fractions

For   and  , gives      , as the Python script appended to Barthélémy (2020) does. The last row of its Table D.2 prints the numerator with the opposite sign, which gives ; a quadrature of the definition settles it, as the application Hill tensors by transformation shows.

No square root of is taken in the first two regimes, so derivatives are exact where the roots coincide. In derivatives, the split   balances the truncation of the second regime against the cancellation of the third, which grows as .

Against the same integrals in 256-bit arithmetic, on generic, strongly anisotropic and nearly isotropic matrices ( to ), isotropic values included, and on from to , the Walpole coefficients are within . Their derivatives with respect to the moduli are within up to  , and   at  .

Elliptical crack near the circle

The crack opening tensor of an elliptical crack of aspect ratio   in an isotropic matrix involves    and   ,   , which cancel as   with a relative error near . Below   they are replaced by the series of the complete integrals,

and   : within   of the same closed forms in 256-bit arithmetic, from the circle to  .

A kink at the circle

depends on the ratio of the smaller semi-axis to the larger, so the roles of and swap at the circle, and has a kink there as a function of either: is smallest at the circle and grows on both sides. Its derivative at the circle is one-sided. ForwardDiff orders two dual numbers of equal value by their partials, so the semi-axis being increased is sorted as the larger one, and the derivative returned is the one in the direction of the increase. Next to the circle, it agrees with the cubature to  .

Cubatures under ForwardDiff

Two traps stop a cubature under differentiation.

• The error estimate. QuadGK compares the Kronrod and Gauss estimates of each segment through a norm. LinearAlgebra.norm rescales by the largest entry and divides by it unless that entry iszero, and a dual number of value zero with non-zero partials is not iszero. On a segment where the two rules agree to the last bit, which a smooth integrand produces, the estimate comes out NaN and the integration stops with a DomainError — the derivative of the anisotropic conduction tensor of a sphere, for one. Every cubature of the package hands QuadGK the norm of the values: the refinement makes the decisions an evaluation without derivatives makes, so the value of a differentiated evaluation is that evaluation itself, and the derivatives ride on the same nodes.

• The integration bounds. QuadGK has no rule for a dual interval: bounds taken in the number type of the crack make the derivative with respect to a crack semi-axis in an anisotropic matrix a MethodError. The bounds are plain numbers; the derivatives live in the integrand.

Self-consistent estimates of transversely isotropic class

The self-consistent scheme projects each iterate on the major-symmetric tensors. The products of transversely isotropic tensors that build an iterate land in the six- or eight-parameter Walpole form. Projected by rebuilding them as generic tensors, the class is lost, and the Hill tensor of a coaxial spheroid in the next iterate goes through a cubature. In the Walpole basis the major transpose swaps and and changes the sign of the major-antisymmetric , , so the projection is the five-parameter tensor

exact, and the estimate keeps its class. On a 3:1 spheroid the scheme then uses the closed form throughout: 46 µs instead of 4.2 ms, and the fixed point agrees with the one reached through the cubature at tight tolerances to  .

Where to go next

The next chapter leaves the ellipsoid for patterns whose inner strain is not uniform: Layered sphere solves the generalized Eshelby problem of a concentric multilayer. The closed forms treated here are those of Hill polarization tensors and COD tensors, and the self-consistent iteration that relies on them is in Self-consistent schemes.