API — Viscoelasticity

MeanFieldHom.Viscoelasticity

Ageing linear viscoelastic (ALV) homogenization. Provides:

• ViscoLaw — relaxation R(t,t') or creep J(t,t') kernel, scalar- or 4-tensor-valued, with built-in Maxwell / Kelvin constructors.
• trapezoidal_matrix — discretization of the Stieltjes integral on a time grid into a lower-block-triangular matrix (n×n for scalar, 6n×6n for 4-tensor in Mandel form).
• volterra_inverse — block forward-substitution that takes a discrete relaxation kernel to the corresponding creep kernel (and vice versa).
• iso_params_from_blocks / iso_blocks_from_params and their ti_ / ortho_ counterparts — conversions between symmetry-structured per-component scalar matrices and the full 6n×6n block matrix.
• hill_kernel — discrete ALV Hill polarization tensor for an ellipsoidal inclusion, isotropic-matrix branch using the time-space decoupling formula ([7], App. ALV Hill kernel).
• Time-domain viscoelastic homogenization schemes (Voigt, Reuss, Dilute, DiluteDual, Mori-Tanaka, Maxwell, Self-Consistent), plugged into the existing homogenize dispatcher whenever a phase carries a ViscoLaw property.

All ALV operators are stored as dense Matrix{T} of size (B·n)×(B·n) (B = 6 for 4-tensor, B = 1 for scalar) with explicit zeros above the block diagonal — this is the convention of [56] and the C++ ECHOES reference.

AbstractViscoLaw

Root supertype for viscoelastic kernels. Concrete subtypes provide an eval(law, t, t_p) method returning a scalar or a 4-tensor.

ViscoLaw(eval_fun, mode::Symbol = :relaxation)

Concrete viscoelastic law. Wraps an eval_fun::F callable (t, t_p) -> X together with its mode Symbol (:relaxation or :creep). The output type X may be:

• a Real or Complex scalar (scalar Volterra kernel) ;
• a TensND.AbstractTens{4,3} (4-tensor relaxation / creep tensor) ;
• an AbstractMatrix of size 6×6 already in Mandel form.

Construct via ViscoLaw directly or use the convenience constructors maxwell_relaxation, kelvin_creep, maxwell_iso, kelvin_iso.

visco_mode(law) -> Symbol

Return the mode (:relaxation or :creep) of law.

visco_eval(law, t, t_p)

Evaluate the kernel at (t, t_p).

maxwell_relaxation(C_inf, C_branches, taus; mode = :relaxation)

Build a generalized Maxwell relaxation kernel

R(t, t') = C_inf + Σ_i C_branches[i] · exp(-(t - t')/taus[i])

for t ≥ t', 0 otherwise. C_inf and C_branches[i] may be scalars or 4-tensors (TensND.AbstractTens{4,3}); taus[i] is a positive relaxation time of the i-th branch.

kelvin_creep(J_0, J_branches, taus; mode = :creep)

Build a Kelvin (or Kelvin-Voigt-Generalized) creep kernel

J(t, t') = J_0 + Σ_i J_branches[i] · (1 - exp(-(t - t')/taus[i]))

for t ≥ t', 0 otherwise. J_0 is the instantaneous compliance, J_branches[i] the i-th branch compliance, taus[i] its retardation time.

maxwell_iso(k, mu, eta_k, eta_mu) -> ViscoLaw

Convenience: build an isotropic Maxwell relaxation 4-tensor kernel

R(t, t') = 3 k · exp(-(t - t')/eta_k) · 𝕁
         + 2 mu · exp(-(t - t')/eta_mu) · 𝕂

with 𝕁 the spherical projector (1/3) 𝟙 ⊗ 𝟙 and 𝕂 = 𝕀 - 𝕁 the deviatoric projector. The output is a TensND.TensISO{4,3} at every (t, t') with t ≥ t'.

kelvin_iso(k_0, mu_0, k_branches, mu_branches, taus_k, taus_mu) -> ViscoLaw

Convenience: build an isotropic Kelvin creep 4-tensor kernel

J(t, t') = (1/(3 k_0)) 𝕁 + (1/(2 mu_0)) 𝕂
         + Σ_i (1/(3 k_branches[i])) (1 - exp(-(t-t')/taus_k[i])) 𝕁
         + Σ_i (1/(2 mu_branches[i])) (1 - exp(-(t-t')/taus_mu[i])) 𝕂

k_branches, mu_branches, taus_k, taus_mu may be empty if no Kelvin branches are required (instantaneous-only compliance).

heaviside_law(C; mode = :relaxation)

Return a ViscoLaw corresponding to a purely elastic kernel R(t,t') = C · H(t-t'). Useful for testing the elastic limit of the ALV pipeline.

trapezoidal_matrix(law::ViscoLaw, times::AbstractVector{<:Real}) -> Matrix

Build the discrete (B·n) × (B·n) lower-block-triangular matrix representing the Stieltjes integral $y(t_i) = ∫_{t_0}^{t_i} f(t_i, τ) \, dx(τ)$ on the time grid times = (t_0, …, t_{n-1}) using the trapezoidal rule of [56]. B = 1 for scalar-valued kernels and B = 6 for 4-tensor- or 6×6-matrix-valued kernels (Mandel form).

The exact block layout is described at the top of trapezoidal.jl. The output is type-stable in the element type returned by visco_eval(law, ...).

_trapezoidal_relaxation(law::ViscoLaw, times, B) -> Matrix

Build the discrete relaxation block matrix from a ViscoLaw, regardless of whether the law is in :relaxation or :creep mode. When the law is :creep, the trapezoidal compliance matrix is inverted (block forward substitution at block_size = B) to obtain the corresponding relaxation matrix — the convention every ALV scheme assumes internally.

B is the block size (6 for order-4 4-tensor / Mandel, 3 for order-2 vector-tensor).

volterra_inverse(M::AbstractMatrix; block_size::Int = 6) -> Matrix

Compute the Volterra inverse of a lower-block-triangular matrix M of size (B·n) × (B·n) with B = block_size. The result is also lower-block-triangular and satisfies M * volterra_inverse(M) = H 𝟙 (block-diagonal identity).

block_size must be a positive divisor of size(M, 1) and size(M, 2), typically 1 (scalar Volterra kernel) or 6 (4-tensor in Mandel form).

The cost is O(B³ n²) flops via block forward-substitution; each diagonal block is inverted with the dense inv(...) of stdlib (so B = 6 costs only a few hundred flops per block).

volterra_inverse(K::AbstractALVKernel) -> AbstractALVKernel

Volterra inverse of a structured ALV kernel. Stays in the same symmetry class : iso ↦ iso, TI ↦ TI, ortho ↦ ortho. Avoids materialising the (6n × 6n) matrix.

volterra_product(A::AbstractMatrix, B::AbstractMatrix) -> Matrix

Discrete Volterra product A ∘ B: a regular matrix multiplication of two lower-block-triangular matrices, returned as a fresh Matrix. The result is again lower-block-triangular by construction.

Equivalent to A * B, exposed as a named function to make the viscoelastic algebra explicit at call sites.

volterra_divide(M, S; block_size = 1) -> Matrix

Compute T = M ∘ S^{-vol} (Volterra-divide) by direct block forward substitution on the linear system T · S = M. Numerically stable where M * volterra_inverse(S) would lose precision through a huge intermediate inverse, e.g. for soft-phase moduli (κ, μ → 0) or step-activated ViscoLaw kernels in multi-layer ALV recurrences.

block_size must divide size(M, 1); typical values are 1 (scalar Volterra) and 6 (4-tensor Mandel). The result is lower-block-triangular if both M and S are.

volterra_left_divide(S, M; block_size = 1) -> Matrix

Compute T = S^{-vol} ∘ M (Volterra LEFT-divide) by direct forward substitution on the linear system S · T = M.

Use this rather than volterra_divide when the closed-form algebra requires the inverse on the LEFT of the numerator (e.g. the Hervé–Zaoui bulk/shear transition matrices, where T = M_b^{-1} · M_a). On a non-uniform time grid Volterra trapezoidal matrices do not form a commutative algebra, so the order matters and right-vs-left divides give different results.

block_size must divide size(M, 1); typical values are 1 (scalar Volterra) and 6 (4-tensor Mandel). Both arguments must be lower- block-triangular of the same size; the result is lower-block-triangular.

volterra_left_divide(S::AbstractALVKernel, M::AbstractALVKernel)
    -> AbstractALVKernel

Volterra left-divide T = S^{-vol} ∘ M within the structured class. Auto-promotes mixed inputs (e.g. iso S + TI M) to the more general class before solving.

iso_params_from_blocks(M) -> (α::Matrix, β::Matrix)

Decompose a 6n×6n block matrix whose every 6×6 block is an isotropic 4-tensor in Mandel form into the two scalar parameter matrices α and β, both of size n × n. The block layout is

  M_block(i,j) = α[i,j] · 𝕁_Mandel + β[i,j] · 𝕂_Mandel

where 𝕁_Mandel = (1/3) e₁ e₁ᵀ (e₁ the Mandel-1 unit) augmented with the symmetric off-diagonal part (1/3) for the upper 3×3 block, and 𝕂_Mandel = 𝕀_Mandel - 𝕁_Mandel.

α[i,j] and β[i,j] are extracted from the diagonal entry M_block[1,1] and the (4,4) Mandel-shear entry of each block: α[i,j] = (M_block[1,1] + 2 M_block[1,2]) β[i,j] = M_block[4,4].

iso_blocks_from_params(α::Matrix, β::Matrix) -> Matrix

Inverse of iso_params_from_blocks: build a 6n×6n block matrix whose every block is the iso 4-tensor α[i,j] · 𝕁 + β[i,j] · 𝕂 in Mandel form.

Both α and β must be n × n and have a common element type.

Note: a kron(α, 𝕁_M) + kron(β, 𝕂_M) formulation was experimented (P3.3.1) and reverted — the two extra (6n × 6n) intermediate allocations made it ~3× slower than this single-pass scalar loop, which writes each entry exactly once.

ti_params_from_blocks(M; axis = (0, 0, 1))
    -> NTuple{6, Matrix{T}}

Decompose a 6n×6n block matrix whose every 6×6 block is a TI 4-tensor with axis n (only n = e₃ is currently supported) into the six scalar Walpole parameter matrices (ℓ₁, ℓ₂, ℓ₃, ℓ₄, ℓ₅, ℓ₆), each of size n × n.

ti_blocks_from_params(ℓ::NTuple{6, AbstractMatrix}; axis = (0, 0, 1))
    -> Matrix{T}

Inverse of ti_params_from_blocks: rebuild a 6n×6n block matrix whose every 6×6 block has the TI Mandel structure with axis n = e₃ and Walpole coefficients (ℓ₁[i,j], …, ℓ₆[i,j]).

ortho_params_from_blocks(M; axes = ((1,0,0),(0,1,0),(0,0,1)))
    -> NTuple{12, Matrix{T}}

Decompose a 6n × 6n block matrix whose every 6×6 block is an ortho 4-tensor with material frame (e₁, e₂, e₃) (only the canonical frame is currently supported) into the 12 scalar parameter matrices

(o₁₁, o₁₂, o₁₃, o₂₁, o₂₂, o₂₃, o₃₁, o₃₂, o₃₃, o₄, o₅, o₆),

each of size n × n, where the 3×3 normal block of the Mandel form is stored row-major in entries 1..9 and the 3 shear-diagonal Mandel entries (M[4,4], M[5,5], M[6,6]) map to (o₄, o₅, o₆).

Note that closure under Volterra product does not preserve major symmetry (o₁₂ ≠ o₂₁ in general), so the full 9-entry normal block is needed even for major-symmetric inputs.

ortho_blocks_from_params(o::NTuple{12, AbstractMatrix};
                          axes = ((1,0,0),(0,1,0),(0,0,1)))
    -> Matrix{T}

Inverse of ortho_params_from_blocks: rebuild a 6n × 6n block matrix whose every 6×6 block has the ortho Mandel structure with axes (e₁, e₂, e₃). All off-block-diagonal entries are zero (block 1..3 ↔ block 4..6 coupling is forbidden by orthotropic symmetry).

hill_kernel(ell::AbstractEllipsoidalInclusion,
            C0_law::ViscoLaw,
            times::AbstractVector{<:Real}) -> Matrix

Discrete ALV Hill polarization tensor kernel for the inclusion ell in an isotropic ALV matrix described by C0_law (a ViscoLaw returning a TensND.TensISO{4,3} 4-tensor at each (t, t')).

The output is a (6n × 6n) lower-block-triangular Matrix{T} with n = length(times), in the same Mandel block convention as trapezoidal_matrix.

Implementation follows the time-space decoupling formula of the ECHOES manual appendix viscoelastic_hill_kernel.qmd:

1. Discretize the matrix kernel C0_law on times to a 6n × 6n block matrix R̃_M.
2. Extract the iso scalar parameter matrices (α, β) of R̃_M (α = 3K, β = 2μ per Mandel convention).
3. Build the longitudinal M_long = (α + 2β)/3 and the shear M_shear = β/2 n×n matrices.
4. Take the scalar Volterra inverses J_long = M_long^{-vol} and J_shear = M_shear^{-vol}.
5. Compute the Mandel forms U^A_M and V^A_M of the elastic auxiliary tensors tens_UA(ell) and tens_VA(ell).
6. Assemble block-by-block: P̃_block(i,j) = J_long[i,j] · U^A_M + J_shear[i,j] · (V^A_M - U^A_M).

voigt_alv(C_phases::AbstractVector, fractions::AbstractVector) -> Matrix

Voigt (uniform-strain) bound: C_eff = Σ_r f_r · C̃^r. Each C_phases[r] is a (6n × 6n) block matrix and fractions[r] is the volume fraction of phase r.

reuss_alv(C_phases::AbstractVector, fractions::AbstractVector) -> Matrix

Reuss (uniform-stress) bound: invert each C̃^r, do the volume average of compliances, then invert the result.

J̃eff = Σr fr · J̃^r , C̃eff = (J̃_eff)^{-vol}

dilute_alv(C_0, contribs, fractions) -> Matrix

Dilute scheme: C̃_eff = C̃^0 + Σ_r f_r · Ñ^{r,dil} where the dilute contributions Ñ^{r,dil} of each inclusion phase are pre-computed via dilute_contribution_alv (each (6n × 6n)), f_r is the volume fraction of phase r, and the matrix C̃^0 is the reference.

Note that for the dilute scheme the matrix volume fraction f_0 does not appear: the inclusions are treated as if they were embedded in an infinite matrix.

dilute_dual_alv(C_0, contribs, fractions) -> Matrix

DiluteDual scheme: the dilute contribution is computed in compliance space. Equivalent in form to dilute_alv but with the matrix acting on the compliance kernel. Useful when the inclusions are weak (e.g. cracks) — the linearisation around the matrix compliance converges better.

J̃eff = J̃^0 + Σr fr · H̃^{r,dil} C̃eff = (J̃_eff)^{-vol}

with H̃^{r,dil} = (J^r - J^0) ∘ B̃^{r,dil} and B̃^{r,dil} = (𝟙 + (𝟙 - Ã^{r,dil}))^{-vol} … see Sevostianov 2008 for the dual formulation in the elastic case.

Implementation routes through volterra_inverse of C_0 to obtain J^0 and back.

mori_tanaka_alv(C_0, A_duts, contribs, fractions, f_matrix) -> Matrix

Mori-Tanaka scheme: the average matrix-strain is taken as the reference "infinite-matrix" strain, hence each inclusion experiences a perturbed strain proportional to the inverse of the volume-weighted concentration average.

C̃eff = C̃^0 + (Σr fr Ñ^{r,dil}) ∘ (f0 · 𝟙 + Σs fs Ã^{s,dil})^{-vol}

A_duts[r] is the dilute concentration kernel of phase r, contribs[r] = ΔC^r ∘ A_duts[r], fractions[r] its volume fraction in the RVE, and f_matrix the matrix volume fraction f_0.

maxwell_alv(C_0, contribs, fractions; H_0) -> Matrix

Maxwell scheme: the reference is replaced by a "host medium" with a prescribed distribution shape, whose Hill kernel is H_0 = P̃_d. The Volterra-discrete formula is

C̃eff = C̃^0 + Σ̃ ∘ (𝟙 - P̃d ∘ Σ̃)^{-vol}

where Σ̃ = Σ_r f_r Ñ^{r,dil} and the dilute contributions are the ones computed by dilute_contribution_alv.

When the distribution shape coincides with the inclusion shape and the matrix is the reference, Maxwell reduces to Mori-Tanaka.

Pass H_0 = hill_kernel(distribution_shape, C_0_law, times) for the default uniform-sphere distribution.

self_consistent_alv(rve, prop; times,
                    abstol = 1e-10, reltol = 1e-8, maxiters = 200,
                    damping = 0.0, verbose = false,
                    select_best = false) -> Matrix

Self-consistent ALV homogenization. Iterates the symmetric Picard fixed point on the (6n × 6n) block matrix until convergence.

The initial estimate is the discretized matrix kernel C̃^0. Each iteration rebuilds the per-phase Hill kernels using the current estimate's iso parameters, computes the dilute concentration tensors, and forms C̃_{m+1}.

Returns the converged effective relaxation matrix.

Keyword arguments

• abstol — absolute Frobenius tolerance on ‖C̃_{m+1} − C̃_m‖.
• reltol — additive relative tolerance (multiplied by ‖C̃_m‖).
• maxiters — hard iteration cap.
• damping — Picard relaxation 0 ≤ damping < 1 (0 = no damping).
• verbose — print residual norms each iteration.
• select_best— return the best iterate seen (rather than the last) when convergence stalls.

asymmetric_self_consistent_alv(rve::RVE, prop::Symbol; times,
                                abstol = 1e-10, reltol = 1e-8,
                                maxiters = 200, damping = 0.0,
                                verbose = false, select_best = false)
    -> Matrix{T}

Asymmetric self-consistent viscoelastic homogenization. The iteration update reads

`C^{n+1} = C_M + Σ_i f_i (C_i − C_M) ∘ A^{dil,i}(C^n)`,

mirroring the reference ASC form. Returns the (6n × 6n) effective relaxation matrix once the residual ‖C^{n+1} − C^n‖_F falls below abstol + reltol · ‖C^n‖_F (or after maxiters Picard steps).

pcw_alv(C_0, contribs, fractions; H_dist) -> Matrix

Ponte-Castañeda & Willis (1995) viscoelastic homogenization in single-distribution-shape form. The formula is algebraically identical to maxwell_alv ; the only difference is that the Hill kernel H_dist is computed against the distribution shape stored in rve.distribution_shape, not against any individual phase.

Use homogenize_alv with scheme = PonteCastanedaWillis() — the dispatcher reads the RVE's distribution shape and forwards here.

differential_alv(rve::RVE, prop::Symbol; times,
                  nsteps = 100, trajectory = nothing,
                  abstol = 1e-8, reltol = 1e-6, alg = nothing,
                  formulation = :stiffness) -> Matrix{T}

Differential homogenization in ageing linear viscoelasticity, solved as a SciML ODE on the fictitious incorporation time τ ∈ [0, 1] (Norris 1985; user's hand-written DEM note) :

\[\frac{\mathrm d \tilde{\mathbb C}}{\mathrm d \tau} = \sum_\alpha \frac{\mathrm d \varphi_\alpha}{\mathrm d \tau} (\tilde{\mathbb C}_\alpha - \tilde{\mathbb C}) \circ \tilde{\mathbb A}_\alpha^{dil}(\tilde{\mathbb C}) + \sum_c \frac{\mathrm d \varepsilon_c}{\mathrm d \tau} \Delta\tilde{\mathbb C}^{crack}_c(\tilde{\mathbb C})\]

with the volume balance df = (𝟙 − f ⊗ 𝐔)·dφ inverted by Sherman- Morrison for solid phases (cracks contribute their density derivative directly). nsteps is the density of save points along τ ; the integration step is controlled by abstol / reltol.

formulation = :compliance integrates the dual ODE on the creep function J̃ = C̃^{-vol} instead, through H̃_α = −J̃ ∘ Ñ_α ∘ J̃, and inverts the result — the same choice as the elastic DifferentialScheme.

Supported inclusion geometries: ellipsoids / spheroids (through the ALV Hill kernel), LayeredSphere (through the bulk + shear ALV recurrences) and flat cracks (CrackDensity). Because the reference of the ODE is the running medium and the ALV Hill kernel is isotropic only, every phase must keep that medium isotropic — see _alv_diff_keeps_iso; the ODE throws an explicit ArgumentError otherwise instead of returning a wrong answer.

dilute_concentration_alv(C_E, C_0, P) -> Matrix

Dilute strain concentration kernel Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol} ([55] eq. 16). All inputs are (6n × 6n) block matrices ; the result is also (6n × 6n) and lower-block-triangular.

dilute_contribution_alv(C_E, C_0, P) -> Matrix

Dilute strain contribution kernel Ñ = ΔC̃ ∘ Ã^dil ([55] eq. 17). This is the size-independent stiffness contribution of a single inclusion.

_alv_diff_keeps_iso(geom, sym, C_r) -> Bool

Whether a solid phase leaves the running effective medium of the differential ALV ODE inside the isotropic class.

This matters only for the differential scheme. The ALV Hill kernel hill_kernel exists for an isotropic reference only, and where Mori-Tanaka or the dilute scheme evaluate it against the (fixed, isotropic) matrix, the differential scheme evaluates it against the running estimate C̃(τ) — which an aligned, non-spherical inclusion progressively takes out of the iso class. Rather than silently reading (α, β) off a matrix that is no longer isotropic, the ODE refuses to start (see differential_alv).

An isotropic orientation average (symmetrize = :iso) restores the property for any shape, which is what randomly oriented inclusions or cracks mean physically.

voigt_alv_iso(αβ_phases, fractions) -> (α_eff, β_eff)

Iso-form Voigt bound: α_eff = Σ_r f_r α_r, β_eff = Σ_r f_r β_r.

reuss_alv_iso(αβ_phases, fractions) -> (α_eff, β_eff)

Iso-form Reuss bound: invert per-phase compliances, average, invert back — done independently on the two scalar Volterra components.

dilute_alv_iso(αβ_0, contribs_iso, fractions) -> (α_eff, β_eff)

Iso-form Dilute scheme: αβ_eff = αβ_0 + Σ_r f_r · (Ñ_r in iso form).

dilute_dual_alv_iso(αβ_0, contribs_compliance_iso, fractions) -> (α_eff, β_eff)

Iso-form DiluteDual: invert the matrix to compliance space, average, invert back.

mori_tanaka_alv_iso(αβ_0, A_duts_iso, contribs_iso, fractions, f_M)
    -> (α_eff, β_eff)

Iso-form Mori-Tanaka: C̃_eff = C̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}, reduced to two scalar Volterra problems on (α, β).

maxwell_alv_iso(αβ_0, contribs_iso, fractions, αβ_H_0) -> (α_eff, β_eff)

Iso-form Maxwell scheme: C̃_eff = C̃_0 + Σ̃ ∘ (𝟙 - P̃_d ∘ Σ̃)^{-vol}, reduced to two scalar Volterra problems on (α, β).

dilute_concentration_alv_iso(αβ_E, αβ_0, αβ_P) -> (α_A, β_A)

Iso-form dilute concentration Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol} reduced to two scalar Volterra problems on (α_E − α_0, β_E − β_0).

dilute_contribution_alv_iso(αβ_E, αβ_0, αβ_P) -> (α_N, β_N)

Iso-form dilute contribution Ñ = ΔC̃ ∘ Ã^dil.

voigt_alv_ti(ℓ_phases, fractions) -> NTuple{6, Matrix}

TI-form Voigt bound: ℓᵢ_eff = Σ_r f_r ℓᵢ_r for each Walpole component.

reuss_alv_ti(ℓ_phases, fractions) -> NTuple{6, Matrix}

TI-form Reuss bound: invert per-phase compliances in TI Walpole form, average, invert back.

dilute_alv_ti(ℓ_0, contribs_ti, fractions) -> NTuple{6, Matrix}

TI-form Dilute scheme: ℓ_eff = ℓ_0 + Σ_r f_r · Ñ_r (component-wise).

dilute_dual_alv_ti(ℓ_0, contribs_compliance_ti, fractions)
    -> NTuple{6, Matrix}

TI-form DiluteDual: invert to compliance Walpole form, average, invert back.

mori_tanaka_alv_ti(ℓ_0, A_duts_ti, contribs_ti, fractions, f_M)
    -> NTuple{6, Matrix}

TI-form Mori-Tanaka: C̃_eff = C̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}, all in TI Walpole form.

maxwell_alv_ti(ℓ_0, contribs_ti, fractions, ℓ_H_0) -> NTuple{6, Matrix}

TI-form Maxwell scheme.

dilute_concentration_alv_ti(ℓ_E, ℓ_0, ℓ_P) -> NTuple{6, Matrix}

TI-form dilute concentration Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol} reduced to a (2n)×(2n) block-Volterra inverse + two n×n scalar inverses.

dilute_contribution_alv_ti(ℓ_E, ℓ_0, ℓ_P) -> NTuple{6, Matrix}

TI-form dilute contribution Ñ = ΔC̃ ∘ Ã^dil.

voigt_alv_ortho(o_phases, fractions) -> NTuple{12, Matrix}

Ortho-form Voigt bound: oₖ_eff = Σ_r f_r oₖ_r for each component.

reuss_alv_ortho(o_phases, fractions) -> NTuple{12, Matrix}

Ortho-form Reuss bound: invert per-phase compliances in ortho form, average, invert back.

dilute_alv_ortho(o_0, contribs_ortho, fractions) -> NTuple{12, Matrix}

Ortho-form Dilute scheme: o_eff = o_0 + Σ_r f_r · Ñ_r.

dilute_dual_alv_ortho(o_0, contribs_compliance_ortho, fractions)
    -> NTuple{12, Matrix}

Ortho-form DiluteDual: invert to compliance ortho form, average, invert back.

mori_tanaka_alv_ortho(o_0, A_duts_ortho, contribs_ortho, fractions, f_M)
    -> NTuple{12, Matrix}

Ortho-form Mori-Tanaka: C̃_eff = C̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}, all in ortho form.

maxwell_alv_ortho(o_0, contribs_ortho, fractions, o_H_0) -> NTuple{12, Matrix}

Ortho-form Maxwell scheme.

dilute_concentration_alv_ortho(o_E, o_0, o_P) -> NTuple{12, Matrix}

Ortho-form dilute concentration Ã^dil = (𝟙 + P̃ ∘ ΔC̃)^{-vol}.

dilute_contribution_alv_ortho(o_E, o_0, o_P) -> NTuple{12, Matrix}

Ortho-form dilute contribution Ñ = ΔC̃ ∘ Ã^dil.

cod_kernel_alv(crack::EllipticCrack, C_M_law::ViscoLaw, times;
               Rn = nothing, Rt = nothing) -> NamedTuple

Discrete ALV COD-tensor data for a penny crack η = 1 in an isotropic ALV matrix. Returns the named tuple (B_n = …, B_t = …) of two n × n scalar Volterra matrices (n = length(times)). Each B̃[i,j] approximates the COD coefficient at the time pair (t_i, t_j).

Interface stiffness (Sevostianov-style spring-like interface)

When the crack carries finite interface stiffness kernels Rn(t,t') (normal) and Rt(t,t') (tangential), pass them as scalar ViscoLaws through the Rn / Rt keyword arguments. The traction-free COD matrices B̃_n, B̃_t are then post-corrected via the algebraic identity

B̃_eff = (b · K + B̃^{-1})^{-vol} = B̃ ∘ (𝟙 + b · K ∘ B̃)^{-vol}

(see [58], [55]) where b = semi_minor of the elliptic crack. Limits :

• Rn / Rt = nothing (default) → traction-free penny limit, recovers the existing B_n, B_t.
• Rn, Rt → ∞ (rigid bonding) → B̃_eff_n, B̃_eff_t → 0 (no opening).

Throws if the matrix law is not iso or the crack is not a penny.

compliance_contribution_alv(crack, C_M_law::ViscoLaw, times) -> Matrix{T}

Discrete (6n × 6n) size-independent compliance contribution H̃ of a penny crack in an isotropic ALV matrix. Computed via the time-space decoupling formula

H̃ = (3/4) · B̃_n · W₁(n̂) + (3/8) · B̃_t · W₆(n̂)

where W₁, W₆ are the canonical Walpole basis tensors of the crack normal axis. When n̂ = e_3 the result is in TI form and routes through the existing TI ALV fast path; arbitrary orientation requires a 6×6 Mandel rotation per (i, j) block (not yet implemented).

Convention: same as the elastic compliance_contribution — the Budiansky-O'Connell density factor is applied separately via delta_compliance_alv.

delta_compliance_alv(crack, H̃, ε) -> Matrix

Apply the Budiansky-O'Connell crack density factor to the size-independent compliance contribution H̃ produced by compliance_contribution_alv, giving the fractional compliance correction ΔJ̃ = (4π/3) ε³ᵈ · H̃ (penny / elliptic geometry, ε³ᵈ = N a b²) — same pre-factor as the elastic case.

stiffness_contribution_alv(sphere, C0_law, times) -> Matrix{T}

(6n × 6n) size-independent ALV stiffness contribution of a layered sphere relative to its iso ALV matrix C0_law. Iso parameters (α-, β-blocks of the assembled matrix) are α = 3 Σ_k f_k (M_κ_k − M_κ_0) ∘ α_k, β = 2 Σ_k f_k (M_μ_k − M_μ_0) ∘ β_k, where α_k, β_k are the per-layer localization matrices and M_κ_k, M_μ_k the per-layer Volterra moduli.

The dilute-scheme effective stiffness with this inclusion at volume fraction f is C̃_eff = C̃_0 + f · stiffness_contribution_alv(sphere, …).

stiffness_contribution_alv(crack, C_ref, times) -> Matrix{T}

Discrete (6n × 6n) crack stiffness contribution Ñ = − C̃_ref ∘ H̃ ∘ C̃_ref, mirror of the elastic [stiffness_contribution(crack, C₀)] formula. C_ref may be a ViscoLaw (the matrix law — relaxation auto-built through _trapezoidal_relaxation) or a pre-discretized (6n × 6n) reference matrix (used by SC iterations against the running estimate C_n).

stiffness_contribution_alv_at(sphere::LayeredSphere, C_ref::AbstractMatrix, times)

Layered-sphere counterpart of stiffness_contribution_alv_at for cracks: the same size-independent ALV contribution, but against a pre-discretized (6n × 6n) reference — the running effective medium of the differential (or self-consistent) ODE rather than the matrix law.

The reference must be isotropic, as everywhere in the layered-sphere ALV recurrences.

stiffness_contribution_alv_at(crack, C_ref::AbstractMatrix) -> Matrix

Variant that takes a pre-discretized (6n × 6n) reference matrix. The compliance contribution is recomputed from the iso parameters of C_ref (only iso ALV matrices are currently supported by compliance_contribution_alv).

delta_stiffness_alv(crack, Ñ, ε) -> Matrix

Apply the Budiansky-O'Connell crack density factor to the size-independent stiffness contribution Ñ produced by stiffness_contribution_alv, giving the dilute stiffness correction ΔC̃ = (4π/3) ε³ᵈ · Ñ (penny / elliptic) or π ε²ᵈ · Ñ (ribbon). Same pre-factors as the elastic case.

bulk_localization_alv(sphere, C0_law, times) -> NTuple{N, Matrix}

Per-layer bulk localization matrices α_k(t,t') (n × n each), such that <ε_v>_layer_k = ⟨α_k⟩ ∘ ε_v_∞ in the Volterra sense.

For the bulk Y₀ harmonic, the volume-averaged volumetric strain in layer k is exactly 3 A_k (the mode-2 amplitude B_k has u_r = B_k/r² ⇒ traceless contribution). Therefore α_k = A_k · A_M^{-1} with A_M the matrix-side mode-1 amplitude.

Implementation : amplitude-space recurrence with closed-form interface transitions (bulk_amplitude_seq_alv). Stable for arbitrary modulus contrasts (pores, step-activated layers …).

Reference : Hervé-Zaoui 1993 (elastic) ; ECHOES manual ch07 §"n-layer ALV bulk recurrence" ; Sanahuja IJSS 2013.

bulk_state_seq_alv(sphere, C0_law, times)
    -> (inside_states::NTuple{N, Matrix}, s_matrix::Matrix)

Backwards-compatible wrapper that exposes the σ-state form of the bulk recurrence: at each layer k, reconstruct the 2n × n block (u_r ; σ_rr) from the per-layer amplitudes via u_r = A_k r_k + B_k / r_k² and σ_rr = 3 κ_k A_k - 4 μ_k B_k / r_k³.

Most users should rely on bulk_localization_alv directly. This helper is retained for the existing test that asserts the matrix-side state matches the analytical Hervé-Zaoui boundary condition.

bulk_amplitude_seq_alv(sphere, C0_law, times)
    -> (NTuple{N, (A_k, B_k)}, A_M, B_M)

Forward-propagate the bulk mode amplitudes (A_k, B_k) (each an n × n Volterra matrix) layer by layer using closed-form interface transitions. Core regularity sets B_1 = 0 and the amplitude seed is A_1 = I_n. Returns the inner-layer amplitudes plus the matrix-side (A_M, B_M).

Supports PerfectInterface, SpringInterface (primal) and MembraneInterface (dual) — same set as the σ-state recurrence — but with no matrix inversion of a fundamental M(r; κ, μ). The only inversions performed are volterra_divide(_, 3κ + 4μ) which remain stable for any non-degenerate modulus.

shear_localization_alv(sphere, C0_law, times) -> NTuple{N, Matrix}

Per-layer deviatoric ALV localization matrices β_k(t,t') (n × n each), defined by <ε_d>_layer_k = ⟨β_k⟩ ∘ ε_d_∞ in the Volterra sense — the n × n Volterra matrix that maps a unit deviatoric remote strain to the volume-averaged deviatoric strain in layer k.

The Y₂-harmonic recurrence uses a (4n × 4n) time-major fundamental matrix per layer, two probe states with seed (a, b) = (I, 0) and (0, I) propagated outward, and a final (2n × 2n) Volterra solve that picks the linear combination matching unit far-field (a_{N+1}, b_{N+1}) = (I, 0). Per-layer β_k is the mode-1 amplitude block extracted from the combined inside state.

Reference : ECHOES manual ch07 §"n-layer ALV shear recurrence" ; Hervé-Zaoui 1993 generalized to ALV via [56].

strain_strain_loc_alv(sphere, C0_law, times) -> Matrix{T}

(6n × 6n) block matrix describing the volume-averaged strain-strain localization across the entire layered sphere under a unit Volterra far-field strain. In iso form this is A_avg = ⟨α⟩ 𝕁 + ⟨β⟩ 𝕂 with ⟨α⟩ = Σ_k f_k α_k(t,t') and ⟨β⟩ = Σ_k f_k β_k(t,t') (Volterra products).

This is the analog used by the ALV dilute / MT / Maxwell schemes when the inclusion phase is a LayeredSphere.

hill_kernel_order2(ell, K_0_law::ViscoLaw, times) -> Matrix

Build the discrete Hill kernel P̃ (size (3n × 3n)) for an ellipsoidal inclusion in an isotropic ALV matrix. Uses the time-space decoupling formula P̃[block(i,j)] = α₀^{-vol}[i,j] · 𝐈^A.

voigt_alv_order2(matrices, fractions) -> Matrix

Order-2 Voigt bound: K̃_eff = Σ_r f_r K̃_r.

reuss_alv_order2(matrices, fractions) -> Matrix

Order-2 Reuss bound: invert each compliance, average, invert back.

dilute_alv_order2(K_0, contribs, fractions) -> Matrix

Order-2 Dilute scheme: K̃_eff = K̃_0 + Σ_r f_r Ñ_r.

dilute_dual_alv_order2(K_0, contribs_compliance, fractions) -> Matrix

Order-2 DiluteDual: invert to compliance space, average, invert back.

mori_tanaka_alv_order2(K_0, A_duts, contribs, fractions, f_M) -> Matrix

Order-2 Mori-Tanaka: K̃_eff = K̃_0 + (Σ_r f_r Ñ_r) ∘ (f_0 𝟙 + Σ_s f_s Ã_s)^{-vol}.

maxwell_alv_order2(K_0, contribs, fractions; H_0) -> Matrix

Order-2 Maxwell scheme. H_0 is the Hill kernel of the (matrix-only) distribution shape — defaults to a sphere when not specified by homogenize_alv_order2.

homogenize_alv_order2(rve, scheme, prop; times)

Backwards-compatible alias for homogenize_alv when the matrix property is order-2. New code should call homogenize_alv directly — the dispatch on order-2 vs order-4 is automatic from the law sample.

AbstractALVKernel{T} <: AbstractMatrix{T}

Abstract supertype for structured ALV kernel wrappers. Concrete subtypes (ALVKernelISO, ALVKernelTI, ALVKernelOrtho) store the symmetry-class parameters compactly and present a (6n × 6n) AbstractMatrix view via lazy getindex.

ALVKernelISO{T}(α::Matrix{T}, β::Matrix{T})
ALVKernelISO(M::AbstractMatrix)

Iso ALV kernel: stores the two n × n Volterra parameter matrices α and β (with α = 3K(t,t') and β = 2μ(t,t') in Mandel form). Materialises as a (6n × 6n) block matrix on demand via Matrix(K) or getindex.

ALVKernelTI{T}(ℓ::NTuple{6, Matrix{T}}; axis = (0, 0, 1))
ALVKernelTI(M::AbstractMatrix; axis = (0, 0, 1))

TI ALV kernel: stores the six n × n Walpole parameter matrices (ℓ₁, ℓ₂, ℓ₃, ℓ₄, ℓ₅, ℓ₆) with the canonical axis (currently only e₃ supported). Materialises as a (6n × 6n) block matrix on demand.

ALVKernelOrtho{T}(o::NTuple{12, Matrix{T}}; axes = canonical)
ALVKernelOrtho(M::AbstractMatrix; axes = canonical)

Ortho ALV kernel: stores the twelve n × n parameter matrices of the ortho closure (9 entries of the full unsymmetric 3×3 normal block in Mandel form + 3 shears) with the canonical material frame (e₁, e₂, e₃). Materialises as a (6n × 6n) block matrix on demand.

homogenize_alv(rve, scheme, prop::Symbol; times) -> Matrix

ALV pipeline: build the discrete block matrices of every phase, compute the ALV Hill kernel and dilute concentration tensors, and dispatch on scheme to the corresponding _alv scheme function.

The function dispatches on the order of the matrix property (read once from the matrix ViscoLaw sample type):

• order-4 (4-tensor / 6×6 Mandel kernel) → returns (6n × 6n) relaxation matrix following the standard Hill-kernel + dilute-concentration pipeline.
• order-2 (2-tensor / 3×3 kernel; conductivity / diffusion / permittivity) → returns (3n × 3n) matrix via the order-2 ALV pipeline (hill_kernel_order2, time-space decoupling).

Supports two inclusion-geometry families (order-4 only):

• Single-shape ellipsoidal (Ellipsoid, Spheroid): the standard Hill-kernel + dilute-concentration pipeline.
• LayeredSphere: bulk + shear ALV recurrences (see bulk_localization_alv and shear_localization_alv) feed stiffness_contribution_alv and strain_strain_loc_alv directly — no Hill kernel is needed. In this case phase_property(rve, name, :C) is ignored (the per-layer moduli stored in the geometry are used instead); pass any ViscoLaw (e.g. heaviside_law(C_0)) as a placeholder.

has_visco_property(rve, prop::Symbol = :C) -> Bool

Return true if any phase of rve carries a ViscoLaw under the key prop.