Layered sphere — bulk + shear recurrences and imperfect interfaces

LayeredSphere is an n-layer isotropic spherical composite inclusion in an infinite isotropic matrix, with per-layer localization, global contribution tensors and layer / sphere / cumulative averages. The bulk and shear recurrences follow [32] (generalizing the three-phase model of [33]); imperfect interfaces follow [34] — PerfectInterface plus one primal/dual pair per physics:

Elasticity (primal / dual)	Conductivity (primal / dual)
SpringInterface(kn, kt)	KapitzaInterface(ρ)
MembraneInterface(κs, μs) (surface elasticity, [35])	SurfaceConductiveInterface(ks) (Barthélémy-Bignonnet 2020)

Convention

Radii are stored in ascending order from the center,

with layer occupying   . Layer 1 is the core, layer is the outermost shell, and the composite sphere is embedded in an infinite matrix for  .

That is the whole problem the recurrences below solve: a concentric pattern in an unbounded reference medium , loaded by a remote uniform strain. Because the strain is not uniform inside such a pattern, it has no Hill tensor at all; what it does have is a volume-averaged concentration tensor, and that is what every scheme consumes.

The three-phase model of Christensen & Lo [33] is one use of this solution rather than a variant of it — take   and let the reference medium be the unknown effective one, and the fixed point is their result. That is a property of the scheme, not of the pattern, so it lives with the schemes.

An interactive view of the layered geometry is in The inclusion zoo.

Moduli are TensISO{4,3} (elasticity) or TensISO{2,3} (conductivity). Interface conditions at each radius are specified in an NTuple{N, AbstractInterface} (default all PerfectInterface).

Bulk (spherical) recurrence — Hervé-Zaoui 1993

Under a purely hydrostatic remote strain, the displacement in layer is u_r^{(k)}(r) = A_k r + B_k / r^2. To stay regular in the incompressibility limit  , the implementation propagates the field-valued state vector   directly, with the intra-layer transfer

where     and     with   . Every entry stays finite as   (hence also for  ), and the per-layer bulk localization   degenerates smoothly (  for an incompressible core).

The entry-point at   is written in the "pressure amplitude" parameterization  , so that

Interface jump matrices

Each interface type provides a 2×2 (bulk) jump matrix such that   .

SpringInterface and KapitzaInterface encode a primal discontinuity (displacement / temperature jump), while MembraneInterface and SurfaceConductiveInterface encode a dual discontinuity (traction / flux jump). All limit to PerfectInterface when their compliance goes to zero.

Conductivity recurrence (Y₁ harmonic)

Under a remote uniform temperature gradient, the temperature field has a Y₁ dependence. The state vector   (amplitudes projected onto the remote gradient direction) propagates through a 2×2 transfer matrix   with  .

Interface jumps for conductivity are given above (Kapitza primal, SurfaceConductive dual, matching the structural pattern of their elastic analogs). The per-layer gradient localization   reduces, in the single-layer case, to the classical   of Maxwell-type composites.

Type genericity & incompressibility

The recurrence consists of small-size matrix arithmetic over the element type; it is exercised with Float64, BigFloat, ForwardDiff.Dual, SymPy.Sym, and Symbolics.Num. Symbolically,

using SymPy; @syms κ₀ μ₀ κ₁ μ₁
s = LayeredSphere((Sym(1),), (TensISO{3}(3κ₁, 2μ₁),))
simplify(MeanFieldHomogenization.LayeredSpheres._bulk_localization(s, κ₀, μ₀)[1])
# → (3κ₀ + 4μ₀) / (3κ₁ + 4μ₀)

and the derivative with respect to any modulus or radius is obtained by wrapping the computation in ForwardDiff.derivative / ForwardDiff.gradient.

Deviatoric (shear) recurrence — Y₂-harmonic 4×4 state vector

Under a remote pure-deviatoric strain, the displacement field in an isotropic layer has the axisymmetric form   ,   , and the four linearly-independent Navier solutions at   are parametrized by the power-law exponents     with material- dependent ratios derived directly from the Navier characteristic equation (using  ):

Mode	Radial dependence
1		— uniform deviatoric strain
2
3
4

The corresponding traction amplitudes are obtained from Hooke's law      :

The state vector   combines displacement and physical traction amplitudes; this form is continuous across every perfect interface and rational in , so the recurrence is type-generic (supports Float64, BigFloat, ForwardDiff.Dual, SymPy.Sym, Symbolics.Num) and remains regular in the incompressibility limit  .

Interface jumps at :

• Perfect: identity.

• Spring:  ,   (traction continuous).

• Membrane: surface-elastic 2D shell generates a jump in the tractions driven by the surface-stress divergence.

Seeding at uses the two regular modes (   and   ; the two singular amplitudes    are forced by regularity at the origin). Propagating both probes and solving a 2×2 linear system for the matrix-side far-field   yields the per-layer amplitudes.

The layer localization is not the mode-1 amplitude alone. Mode 2 has an displacement profile, so it integrates to a non-zero deviatoric strain over a shell of finite thickness, whereas modes 3 and 4 average to zero pointwise:

Dropping the mode-2 term is invisible on degenerate configurations (vanishing core, core ≡ shell, single layer) and wrong by 1–50 % on a genuine multi-layer stack.

For   the recurrence reduces to the classical Eshelby single- sphere result; for   it reproduces the [33] core-shell effective shear modulus and passes the Eshelby consistency tests (  with core ≡ shell ↔ single-layer of radius , etc.).

Averages (Echoes-style)

Three volume-average flavors are provided:

• layer_strain_average(sphere, C₀, ε∞, k) — mean strain in layer (bulk + deviatoric parts).

• sphere_strain_average(sphere, C₀, ε∞) — mean strain in the whole composite.

• cumulative_strain_average(sphere, C₀, ε∞, r) — mean strain inside the ball of radius .

All three cover the deviatoric part for any   via the shear recurrence above.

Pointwise fields

The recurrences above already carry everything needed to evaluate the field at a point, in any layer and in the matrix; only the reconstruction was missing. Write  ,   ,   .

Spherical part.    with   , hence     .

Deviatoric part. The general — non-axisymmetric —   solution is

which is the  ,    convention above through   ,  . Differentiating with    gives, with    ,

Every generator is transversely isotropic about — the configuration is rotation-invariant about the center — so the pointwise localization tensor , defined by   , is a TensTI{4,T,6}: six Walpole coefficients and an axis, with no major symmetry.

Averaging over directions and over a shell returns exactly the above — mode 1 contributes , modes 3 and 4 contribute nothing pointwise, and mode 2 contributes  , whose shell average is the factor. shell_localization exposes that identity from the same cached amplitudes, so the pointwise and averaged routes cannot drift apart.

The transport problem is the analog:     and        , a second-order transversely isotropic tensor.

Everything is validated pointwise against the C++ reference (Echoes' loc_eE, loc_eS, loc_sE, loc_sS) to , perfect, spring and Gurtin–Murdoch membrane interfaces alike.

See the worked example with figures: n-layer sphere: pointwise fields.