API — Viscoelasticity

MeanFieldHomogenization.Viscoelasticity Module

MeanFieldHomogenization.Viscoelasticity

Aging linear Volterra (ALV) homogenization: aging linear viscoelasticity at order 4, and its hereditary conduction and diffusion counterpart at order 2. Provides:

• ViscoLaw — relaxation or creep 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 (  for scalar,   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   block matrix.

• hill_kernel — discrete ALV Hill polarization tensor for an ellipsoidal inclusion, isotropic-matrix branch using the time-space decoupling formula (Barthélémy et al., 2016; appendix on the ALV Hill kernel).

• Time-domain homogenization schemes (Voigt, Reuss, Dilute, DiluteDual, Mori-Tanaka, Maxwell, Ponte Castañeda-Willis, self-consistent, asymmetric self-consistent at order 4, differential), ellipsoids, cracks and layered spheres alike, plugged into the existing homogenize dispatcher whenever a phase carries a ViscoLaw property.

All ALV operators are stored as dense Matrix{T} of size   (  for 4-tensor,   for scalar) with explicit zeros above the block diagonal — this is the convention of Sanahuja (2013) and of the C++ library Echoes.

Constitutive laws

MeanFieldHomogenization.Viscoelasticity.AbstractViscoLaw Type

AbstractViscoLaw

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

MeanFieldHomogenization.Viscoelasticity.ViscoLaw Type

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   already in Mandel form.

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

MeanFieldHomogenization.Viscoelasticity.visco_mode Function

visco_mode(law) -> Symbol

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

MeanFieldHomogenization.Viscoelasticity.VALID_VISCO_MODES Constant

VALID_VISCO_MODES :: Tuple{Symbol, Symbol}

Allowed values for ViscoLaw.mode :

• :relaxation — maps strain history to stress.

• :creep — maps stress history to strain.

MeanFieldHomogenization.Viscoelasticity.visco_eval Function

visco_eval(law, t, t_p)

Evaluate the kernel at (t, t_p).

MeanFieldHomogenization.Viscoelasticity.maxwell_relaxation Function

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

Build a generalized Maxwell relaxation kernel

for  , 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 -th branch.

MeanFieldHomogenization.Viscoelasticity.kelvin_creep Function

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

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

for  , otherwise. J_0 () is the instantaneous compliance, J_branches[i] () the -th branch compliance, taus[i] () its retardation time.

MeanFieldHomogenization.Viscoelasticity.maxwell_iso Function

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

Convenience: build an isotropic Maxwell relaxation 4-tensor kernel

with mu (), eta_k (), eta_mu (), the spherical projector   and    the deviatoric projector. The output is a TensND.TensISO{4,3} at every with  .

MeanFieldHomogenization.Viscoelasticity.kelvin_iso Function

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

Convenience: build an isotropic Kelvin creep 4-tensor kernel

with mu_0 (), k_branches[i] (), mu_branches[i] (), taus_k[i] () and taus_mu[i] (). k_branches, mu_branches, taus_k, taus_mu may be empty if no Kelvin branches are required (instantaneous-only compliance).

MeanFieldHomogenization.Viscoelasticity.heaviside_law Function

heaviside_law(C; mode = :relaxation)

Return a ViscoLaw corresponding to a purely elastic kernel    . Useful for testing the elastic limit of the ALV pipeline.

Trapezoidal discretization

MeanFieldHomogenization.Viscoelasticity.trapezoidal_matrix Function

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

Build the discrete   lower-block-triangular matrix representing the Stieltjes integral    on the time grid times using the trapezoidal rule of Sanahuja (2013).   for scalar-valued kernels and   for 4-tensor- or  -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, ...).

MeanFieldHomogenization.Viscoelasticity._trapezoidal_relaxation Function

_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 algebra

MeanFieldHomogenization.Viscoelasticity.volterra_inverse Function

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

Compute the Volterra inverse of a lower-block-triangular matrix M of size   with = block_size. The result is also lower-block-triangular and satisfies M * volterra_inverse(M)   (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 flops via block forward-substitution; each diagonal block is inverted with the dense inv(...) of stdlib (so   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 materializing the   matrix.

MeanFieldHomogenization.Viscoelasticity.volterra_product Function

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

Discrete Volterra product  : 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.

MeanFieldHomogenization.Viscoelasticity.volterra_divide Function

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

Compute    (Volterra-divide) by direct block forward substitution on the linear system   . Numerically stable where M * volterra_inverse(S) would lose precision through a huge intermediate inverse, e.g. for soft-phase moduli ( ) 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.

MeanFieldHomogenization.Viscoelasticity.volterra_left_divide Function

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

Compute    (Volterra LEFT-divide) by direct forward substitution on the linear system   .

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   ). 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    within the structured class. Auto-promotes mixed inputs (e.g. iso S + TI M) to the more general class before solving.

Symmetry-class conversions

MeanFieldHomogenization.Viscoelasticity.iso_params_from_blocks Function

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

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

where    ( the Mandel-1 unit) augmented with the symmetric off-diagonal part for the upper   block, and   .

and are extracted from the diagonal entry and the (4,4) Mandel-shear entry of each block:

MeanFieldHomogenization.Viscoelasticity.iso_blocks_from_params Function

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

Inverse of iso_params_from_blocks: build a   block matrix whose every block is the iso 4-tensor     in Mandel form.

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

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

MeanFieldHomogenization.Viscoelasticity.ti_params_from_blocks Function

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

Decompose a   block matrix whose every   block is a TI 4-tensor with axis (only   is currently supported) into the six scalar Walpole parameter matrices , each of size  .

MeanFieldHomogenization.Viscoelasticity.ti_blocks_from_params Function

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

Inverse of ti_params_from_blocks: rebuild a   block matrix whose every   block has the TI Mandel structure with axis   and Walpole coefficients .

MeanFieldHomogenization.Viscoelasticity.ortho_params_from_blocks Function

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

Decompose a   block matrix whose every   block is an ortho 4-tensor with material frame (only the canonical frame is currently supported) into the 12 scalar parameter matrices

each of size  , where the   normal block of the Mandel form is stored row-major in entries 1..9 and the 3 shear-diagonal Mandel entries map to .

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

MeanFieldHomogenization.Viscoelasticity.ortho_blocks_from_params Function

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   block matrix whose every   block has the ortho Mandel structure with axes . All off-block-diagonal entries are zero (block 1..3 ↔ block 4..6 coupling is forbidden by orthotropic symmetry).

Hill ALV kernel

MeanFieldHomogenization.Viscoelasticity.hill_kernel Function

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 ).

The output is a   lower-block-triangular Matrix{T} with = 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   block matrix .

2. Extract the iso scalar parameter matrices of ( ,   per Mandel convention).

3. Build the longitudinal    and the shear     matrices: the two eigenvalues of the acoustic tensor, which the Hill tensor inverts. The bulk and shear kernels follow linearly,    and  ; their inverses do not (see the theory page).

4. Take the scalar Volterra inverses   and  .

5. Compute the Mandel forms and of the elastic auxiliary tensors tens_UA(ell) and tens_VA(ell).

6. Assemble block-by-block:      .

MeanFieldHomogenization.Viscoelasticity.second_hill_kernel Function

second_hill_kernel(ell, C0_law::ViscoLaw, times) -> Matrix

Discrete ALV second Hill tensor kernel of the inclusion ell in an isotropic ALV matrix described by C0_law, the counterpart of second_hill_tensor with the double contractions replaced by Volterra products:

being the relaxation kernel of the matrix (a creep law is inverted first) and the hill_kernel built on it. Same   lower-block-triangular Mandel layout as trapezoidal_matrix; for a law constant in time each diagonal block is the elastic .

Schemes — generic (6n × 6n) algebra

MeanFieldHomogenization.Viscoelasticity.voigt_alv Function

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

Voigt (uniform-strain) bound:   . Each C_phases[i] is a   block matrix and fractions[i] is the volume fraction of phase .

MeanFieldHomogenization.Viscoelasticity.reuss_alv Function

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

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

MeanFieldHomogenization.Viscoelasticity.dilute_alv Function

dilute_alv(C_0, contribs, fractions) -> Matrix

Dilute scheme:     where the dilute contributions of each inclusion phase are pre-computed via dilute_contribution_alv (each  ), is the volume fraction of phase , and the matrix is the reference.

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

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv Function

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 linearization around the matrix compliance converges better.

with the compliance contribution kernels      , the Volterra counterpart of the elastic identity       (Barthélémy et al., 2019); contribs_compliance holds them.

Implementation routes through volterra_inverse of C_0 to obtain and back.

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv Function

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.

A_duts[i] is the dilute concentration kernel of phase , contribs[i]    , fractions[i] its volume fraction in the RVE, and f_matrix the matrix volume fraction .

MeanFieldHomogenization.Viscoelasticity.maxwell_alv Function

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 . The Volterra-discrete formula is

where    and the dilute contributions are the ones computed by dilute_contribution_alv.

When every inclusion phase has the shape of the distribution and the matrix is the reference, Maxwell reduces to Mori-Tanaka.

Pass H_0 = hill_kernel(distribution_shape, C_0_law, times), the distribution shape being the one declared on the RVE; there is no default.

MeanFieldHomogenization.Viscoelasticity.self_consistent_alv Function

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   block matrix until convergence.

The initial estimate is the discretized matrix kernel . Each iteration rebuilds the per-phase Hill kernels using the current estimate's iso parameters, computes the dilute concentration tensors, and forms .

Returns the converged effective relaxation matrix.

Keyword arguments

• abstol — absolute Frobenius tolerance on  .

• reltol — additive relative tolerance (multiplied by ).

• maxiters — hard iteration cap.

• damping — Picard relaxation   damping (0 = no damping).

• verbose — print residual norms each iteration.

• select_best— return the best iterate seen (rather than the last) when convergence stalls.

MeanFieldHomogenization.Viscoelasticity.asymmetric_self_consistent_alv Function

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

with the matrix kernel and the iteration index, mirroring the reference ASC form. Returns the   effective relaxation matrix once the residual   falls below abstol + reltol   (or after maxiters Picard steps).

MeanFieldHomogenization.Viscoelasticity.pcw_alv Function

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.

MeanFieldHomogenization.Viscoelasticity.differential_alv Function

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 aging linear viscoelasticity, solved as a SciML ODE on the fictitious incorporation time   (Norris, 1985):

with the volume balance      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   instead, through     , 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.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv Function

dilute_concentration_alv(C_E, C_0, P) -> Matrix

Dilute strain concentration kernel     ,    (Barthélémy et al., 2019; eq. 16). All inputs are   block matrices ; the result is also   and lower-block-triangular.

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv Function

dilute_contribution_alv(C_E, C_0, P) -> Matrix

Dilute strain contribution kernel    (Barthélémy et al., 2019; eq. 17). This is the size-independent stiffness contribution of a single inclusion.

Internal guards

MeanFieldHomogenization.Viscoelasticity._alv_diff_keeps_iso Function

_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 — 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.

Schemes — iso fast path

MeanFieldHomogenization.Viscoelasticity.voigt_alv_iso Function

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

Iso-form Voigt bound:   ,   .

MeanFieldHomogenization.Viscoelasticity.reuss_alv_iso Function

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.

MeanFieldHomogenization.Viscoelasticity.dilute_alv_iso Function

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

Iso-form Dilute scheme:    , evaluated on the iso components.

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_iso Function

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

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

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_iso Function

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

Iso-form Mori-Tanaka:

reduced to two scalar Volterra problems on .

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_iso Function

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

Iso-form Maxwell scheme:

with the Hill kernel of the distribution shape, reduced to two scalar Volterra problems on .

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_iso Function

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

Iso-form dilute concentration     ,   , reduced to two scalar Volterra problems on   .

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_iso Function

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

Iso-form dilute contribution   .

Schemes — TI fast path

MeanFieldHomogenization.Viscoelasticity.voigt_alv_ti Function

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

TI-form Voigt bound:    for each Walpole component .

MeanFieldHomogenization.Viscoelasticity.reuss_alv_ti Function

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

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

MeanFieldHomogenization.Viscoelasticity.dilute_alv_ti Function

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

TI-form Dilute scheme:     (component-wise).

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_ti Function

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

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

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_ti Function

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

TI-form Mori-Tanaka:

all in TI Walpole form.

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_ti Function

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

TI-form Maxwell scheme.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_ti Function

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

TI-form dilute concentration     ,   , reduced to a   block-Volterra inverse + two   scalar inverses.

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_ti Function

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

TI-form dilute contribution   .

Schemes — ortho fast path

MeanFieldHomogenization.Viscoelasticity.voigt_alv_ortho Function

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

Ortho-form Voigt bound:    for each component .

MeanFieldHomogenization.Viscoelasticity.reuss_alv_ortho Function

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

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

MeanFieldHomogenization.Viscoelasticity.dilute_alv_ortho Function

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

Ortho-form Dilute scheme:    .

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_ortho Function

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

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

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_ortho Function

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

Ortho-form Mori-Tanaka:

all in ortho form.

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_ortho Function

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

Ortho-form Maxwell scheme.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_ortho Function

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

Ortho-form dilute concentration     ,   .

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_ortho Function

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

Ortho-form dilute contribution   .

Cracks

The stiffness contribution of a crack is a method of property_contribution_alv, documented with the layered spheres below.

MeanFieldHomogenization.Viscoelasticity.cod_kernel_alv Function

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

Discrete ALV COD-tensor data for a penny crack   in an isotropic ALV matrix. Returns the named tuple (B_n = …, B_t = …) of two   scalar Volterra matrices ( = length(times)). Each approximates the COD coefficient at the time pair .

Interface stiffness (Sevostianov-style spring-like interface)

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

(Sevostianov and Kachanov, 2007; Barthélémy et al., 2019), where is the semi-minor axis (semi_minor) of the elliptic crack. Limits :

• Rn / Rt = nothing (default) → traction-free penny limit, recovers the existing B_n, B_t.

•   (rigid bonding) →   (no opening).

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

MeanFieldHomogenization.Viscoelasticity.inverse_property_contribution_alv Function

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

Discrete   size-independent compliance contribution of a penny crack in an isotropic ALV matrix. Computed via the time-space decoupling formula

where , are the canonical Walpole basis tensors of the crack normal axis. When   the result is in TI form and routes through the existing TI ALV fast path; arbitrary orientation requires a   Mandel rotation per block (not yet implemented).

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

MeanFieldHomogenization.Viscoelasticity.delta_inverse_property_alv Function

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

Apply the Budiansky-O'Connell crack density factor to the size-independent compliance contribution produced by inverse_property_contribution_alv, giving the fractional compliance correction     (penny / elliptic geometry,    ) — same pre-factor as the elastic case.

MeanFieldHomogenization.Viscoelasticity.property_contribution_alv_at Function

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

Layered-sphere counterpart of property_contribution_alv_at for cracks: the same size-independent ALV contribution, but against a pre-discretized   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.

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

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

MeanFieldHomogenization.Viscoelasticity.delta_property_alv Function

delta_property_alv(crack, Ñ, ε) -> Matrix

Apply the Budiansky-O'Connell crack density factor to the size-independent stiffness contribution produced by property_contribution_alv, giving the dilute stiffness correction     (penny / elliptic) or    (ribbon). Same pre-factors as the elastic case.

Layered spheres

One function per notion serves elasticity and conduction, the order of the reference law selecting the method: gradient_gradient_concentration_alv and property_contribution_alv return a   matrix for a fourth-order reference and a   one for a second-order reference.

MeanFieldHomogenization.Viscoelasticity.gradient_gradient_concentration_alv Function

gradient_gradient_concentration_alv(sphere, ref, times; external = true) -> Matrix

Volume-averaged gradient concentration of a whole LayeredSphere in an isotropic aging reference ref (a ViscoLaw, or the Volterra matrix of a running estimate), under a unit Volterra far-field gradient. The order of ref selects the physics.

In elasticity, the   block matrix of the strain concentration,      with     and     (Volterra products), , being the displacement jumps of the spring interfaces. In conduction, the   gradient concentration   , adding the temperature jumps     of the Kapitza interfaces. external decides whether the outer interface belongs to the sphere, as in the elastic gradient_gradient_concentration.

This is what the ALV dilute / Mori–Tanaka / Maxwell schemes use when an inclusion phase is a LayeredSphere. Explicit name in elasticity: strain_strain_concentration_alv.

MeanFieldHomogenization.Viscoelasticity.property_contribution_alv Function

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

Discrete   crack stiffness contribution     , mirror of the elastic [property_contribution(crack, C₀)] formula, with the discretized C_ref. C_ref may be a ViscoLaw (the matrix law — relaxation auto-built through _trapezoidal_relaxation) or a pre-discretized   reference matrix (used by SC iterations against the running estimate C_n).

property_contribution_alv(sphere, ref, times; external = true) -> Matrix

Size-independent ALV contribution of an inclusion to the property of an aging reference ref: the dilute effective property at volume fraction is   .

For a LayeredSphere the order of ref selects the physics. In elasticity the   stiffness contribution has the iso parameters

, being the per-layer concentration matrices, , the per-layer Volterra moduli and , the displacement jumps of the spring interfaces, plus the surface stress of the membranes. In conduction the   conductivity contribution is       plus the surface flux of the surface-conductive interfaces. ref may also be the Volterra matrix of a running estimate, which is how the self-consistent and differential schemes call it; external has the meaning it has in gradient_gradient_concentration_alv.

The crack method, documented next, is the elastic ALV stiffness contribution of a crack family. Explicit names: stiffness_contribution_alv and conductivity_contribution_alv.

The per-layer kernels behind them, of the bulk and shear modes in elasticity and of the single mode of conduction:

MeanFieldHomogenization.Viscoelasticity.bulk_concentration_alv Function

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

Per-layer bulk concentration matrices (  each), such that    in the Volterra sense.

For the bulk harmonic, the volume-averaged volumetric strain in layer is exactly (the mode-2 amplitude has   ⇒ traceless contribution). Therefore    with 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" ; sanahuja2013.

MeanFieldHomogenization.Viscoelasticity.bulk_state_seq_alv Function

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 , reconstruct the   block from the per-layer amplitudes via    and   .

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

MeanFieldHomogenization.Viscoelasticity.bulk_amplitude_seq_alv Function

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

Forward-propagate the bulk mode amplitudes (each an   Volterra matrix) layer by layer using closed-form interface transitions. Core regularity sets   and the amplitude seed is  . 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 . The only inversions performed are volterra_divide by  , which remain stable for any non-degenerate modulus.

MeanFieldHomogenization.Viscoelasticity.shear_concentration_alv Function

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

Per-layer deviatoric ALV concentration matrices (  each), defined by    in the Volterra sense — the   Volterra matrix that maps a unit deviatoric remote strain to the volume-averaged deviatoric strain in layer .

The -harmonic recurrence uses a   time-major fundamental matrix per layer, two probe states with seed   and propagated outward, and a final   Volterra solve that picks the linear combination matching unit far-field  . Per-layer 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 (Sanahuja, 2013).

MeanFieldHomogenization.Viscoelasticity.gradient_concentration_alv Function

gradient_concentration_alv(sphere, K0_law, times) -> NTuple{N, Matrix}

Per-layer gradient concentrations of a LayeredSphere in an isotropic aging reference: the   Volterra matrices with   , averaged over the material of layer (the conduction twin of bulk_concentration_alv). For a non-aging problem they are the elastic ones of gradient_gradient_concentration at every time step.

The names of elasticity and of conduction

Each generic function also answers to the name it has in elasticity, and the contribution to the name it has in conduction. Such a name forwards to the generic function, accepts one order only, and refuses the other with an error that names the generic function. The crack functions are elastic only, and their names in elasticity are plain synonyms.

Generic function	Elasticity (order 4)	Conduction (order 2)
gradient_gradient_concentration_alv	strain_strain_concentration_alv
property_contribution_alv	stiffness_contribution_alv	conductivity_contribution_alv
property_contribution_alv_at	stiffness_contribution_alv_at
inverse_property_contribution_alv	compliance_contribution_alv
delta_property_alv	delta_stiffness_alv
delta_inverse_property_alv	delta_compliance_alv

MeanFieldHomogenization.Viscoelasticity.strain_strain_concentration_alv Function

strain_strain_concentration_alv(sphere, C0_law, times; external = true) -> Matrix

Name in elasticity of gradient_gradient_concentration_alv, for a fourth-order reference only.

MeanFieldHomogenization.Viscoelasticity.stiffness_contribution_alv Function

stiffness_contribution_alv(sphere_or_crack, C0_law, times; kw...) -> Matrix

Name in elasticity of property_contribution_alv, for a fourth-order reference only.

MeanFieldHomogenization.Viscoelasticity.conductivity_contribution_alv Function

conductivity_contribution_alv(sphere, K0_law, times; external = true) -> Matrix

Name in conduction of property_contribution_alv, for a second-order reference only.

MeanFieldHomogenization.Viscoelasticity.stiffness_contribution_alv_at Function

stiffness_contribution_alv_at(sphere_or_crack, C_ref, args...; kw...) -> Matrix

Name in elasticity of property_contribution_alv_at, whose methods are all elastic.

MeanFieldHomogenization.Viscoelasticity.compliance_contribution_alv Function

compliance_contribution_alv(crack, C_M_law, times; Rn, Rt) -> Matrix

Name in elasticity of inverse_property_contribution_alv, whose methods are all elastic.

MeanFieldHomogenization.Viscoelasticity.delta_stiffness_alv Function

delta_stiffness_alv(crack, Ñ, ε) -> Matrix

Name in elasticity of delta_property_alv.

MeanFieldHomogenization.Viscoelasticity.delta_compliance_alv Function

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

Name in elasticity of delta_inverse_property_alv.

Order-2 ALV (conductivity / diffusion)

The kernels below act on Volterra matrices, which do not carry their order in their type: a   matrix with even also has rows. They keep the order in their names, as the iso, TI and ortho fast paths keep their symmetry class; homogenize_alv chooses among them from the property.

MeanFieldHomogenization.Viscoelasticity.hill_kernel_order2 Function

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

Build the discrete Hill kernel (size  ) for an ellipsoidal inclusion in an isotropic ALV matrix. Uses the time-space decoupling formula   , with    the discretized isotropic reference.

MeanFieldHomogenization.Viscoelasticity.second_hill_kernel_order2 Function

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

Discrete ALV second Hill tensor kernel (size  ) of an ellipsoidal inclusion in an isotropic ALV matrix of conductivity law K_0_law:     , with the hill_kernel_order2; the order-2 counterpart of second_hill_kernel.

MeanFieldHomogenization.Viscoelasticity.voigt_alv_order2 Function

voigt_alv_order2(matrices, fractions) -> Matrix

Order-2 Voigt bound:   .

MeanFieldHomogenization.Viscoelasticity.reuss_alv_order2 Function

reuss_alv_order2(matrices, fractions) -> Matrix

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

MeanFieldHomogenization.Viscoelasticity.dilute_alv_order2 Function

dilute_alv_order2(K_0, contribs, fractions) -> Matrix

Order-2 Dilute scheme:    .

MeanFieldHomogenization.Viscoelasticity.dilute_dual_alv_order2 Function

dilute_dual_alv_order2(K_0, contribs_compliance, fractions) -> Matrix

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

MeanFieldHomogenization.Viscoelasticity.mori_tanaka_alv_order2 Function

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

Order-2 Mori-Tanaka:

MeanFieldHomogenization.Viscoelasticity.maxwell_alv_order2 Function

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.

MeanFieldHomogenization.Viscoelasticity.self_consistent_alv_order2 Function

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

Order-2 (conduction, diffusion) self-consistent estimate in the aging linear Volterra setting: the fixed point of      , every concentration being taken in the running estimate , by Picard iteration from the matrix (the phase matrix, or the one with fraction = :rest). A LayeredSphere phase enters through its recurrence. Every phase must keep the estimate isotropic, as the order-2 Hill kernel against it requires: a sphere of isotropic conductivity, a layered sphere, or symmetrize = :iso. Reached through homogenize_alv(rve, SelfConsistent(), :K; times).

MeanFieldHomogenization.Viscoelasticity.homogenize_alv_order2 Function

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.

MeanFieldHomogenization.Viscoelasticity.dilute_concentration_alv_order2 Function

dilute_concentration_alv_order2(K_E, K_0, P) -> Matrix

Order-2 dilute concentration     ,   , all matrices  .

MeanFieldHomogenization.Viscoelasticity.dilute_contribution_alv_order2 Function

dilute_contribution_alv_order2(K_E, K_0, P) -> Matrix

Order-2 dilute contribution   .

MeanFieldHomogenization.Viscoelasticity.iso_order2_params_from_blocks Function

iso_order2_params_from_blocks(M) -> α::Matrix

Decompose a   block matrix whose every   block is   (iso 2-tensor) into the scalar   Volterra matrix .

MeanFieldHomogenization.Viscoelasticity.iso_order2_blocks_from_params Function

iso_order2_blocks_from_params(α::AbstractMatrix) -> Matrix

Inverse of iso_order2_params_from_blocks: build a   block matrix   per block.

Structured ALV kernel types

MeanFieldHomogenization.Viscoelasticity.AbstractALVKernel Type

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   AbstractMatrix view via lazy getindex.

MeanFieldHomogenization.Viscoelasticity.ALVKernelISO Type

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

Iso ALV kernel: stores the two   Volterra parameter matrices α and β (with   and   in Mandel form). Materializes as a   block matrix on demand via Matrix(K) or getindex.

MeanFieldHomogenization.Viscoelasticity.ALVKernelTI Type

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

TI ALV kernel: stores the six   Walpole parameter matrices with the canonical axis (currently only supported). Materializes as a   block matrix on demand.

MeanFieldHomogenization.Viscoelasticity.ALVKernelOrtho Type

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

Ortho ALV kernel: stores the twelve   parameter matrices of the ortho closure (9 entries of the full unsymmetric   normal block in Mandel form + 3 shears) with the canonical material frame . Materializes as a   block matrix on demand.

Orientation symmetrization

MeanFieldHomogenization.Viscoelasticity._maybe_symmetrize_alv Function

_maybe_symmetrize_alv(M, sym) -> Matrix

Apply the orientation-averaging projection corresponding to sym to a   ALV block matrix : NoSymmetrize (passthrough), IsoSymmetrize (block-by-block exact average) and TISymmetrize(axis) (block-by-block exact azimuthal average about the axis).

MeanFieldHomogenization.Viscoelasticity._maybe_symmetrize_alv2 Function

_maybe_symmetrize_alv2(M, sym) -> Matrix

Order-2 counterpart of _maybe_symmetrize_alv, for   ALV matrices whose blocks are 2-tensors.

Public dispatcher

MeanFieldHomogenization.Viscoelasticity.homogenize_alv Function

homogenize_alv(lam::Laminate, ::Laminated, prop; times, kw...)

Aging-viscoelastic effective operator of a periodic multilayer, dispatching on the order of the layer laws exactly as the elastic Laminated scheme dispatches on the order of the layer tensors.

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 /   Mandel kernel) → returns   relaxation matrix following the standard Hill-kernel + dilute-concentration pipeline.

• order-2 (2-tensor /   kernel; conductivity / diffusion / permittivity) → returns   matrix via the order-2 ALV pipeline (hill_kernel_order2, time-space decoupling).

Supports two inclusion-geometry families, in both orders:

• Single-shape ellipsoidal (Ellipsoid, Spheroid): the standard Hill-kernel + dilute-concentration pipeline.

• LayeredSphere: the ALV recurrences feed its concentration and contribution directly — no Hill kernel is needed. At order 4 these are gradient_gradient_concentration_alv and property_contribution_alv (bulk and shear recurrences, bulk_concentration_alv, shear_concentration_alv); at order 2, gradient_gradient_concentration_alv and property_contribution_alv. The phase property is ignored (the layers carry their own kernels); pass any ViscoLaw (e.g. heaviside_law(C_0)) as a placeholder.

symmetrize is honored in both orders, with the projector of the matching tensor order (  Mandel blocks for order 4, 2-tensor blocks for order 2 — see _maybe_symmetrize_alv and _maybe_symmetrize_alv2). Voigt and Reuss average the phase relaxation matrices (Voigt) or their Volterra inverses (Reuss) directly, each over its orientation distribution symmetrize; the geometry does not enter them.

Reference-updating schemes need an isotropic running medium

Every ALV Hill kernel, order 2 and order 4 alike, is built for an isotropic reference: that is what decouples the time and space parts. Voigt, Reuss, Dilute, DiluteDual, MoriTanaka, Maxwell and PCW evaluate it against the fixed, isotropic matrix and accept any inclusion shape or orientation. SelfConsistent and DifferentialScheme evaluate it against their running estimate, so every inclusion phase must keep that estimate isotropic — either a spherical inclusion with an isotropic phase law, or an isotropic orientation average symmetrize = :iso, which is also what randomly oriented inclusions and cracks mean physically. A LayeredSphere, whose contribution is isotropic by construction, is always accepted.

An RVE satisfying neither raises an ArgumentError naming the phase, rather than reading iso parameters off a matrix that is not isotropic. A non-isotropic ALV matrix is refused for the same reason. The elastic pipeline has no such restriction — its Hill tensor is available for anisotropic references.

MeanFieldHomogenization.Viscoelasticity.has_visco_property Function

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

Return true if any member (phase or layer) of cell carries a ViscoLaw under the key prop.

Reads through the raw accessor cell_container_property: this is a type inspection, and going through the resolving accessor would run a whole inner homogenization on every declaratively nested cell just to look at a type. A nested cell is instead recursed into, so a ViscoLaw buried one scale down is still reported.