Thermal expansion and Biot tensor

Before this page

The theory page Levin's theorem, which states and proves the two formulas used below, and the manual section Concentration tensors of the phases, which shows how the concentration tensors of a scheme are obtained and read.

A homogenization scheme gives the effective stiffness of a composite, and with it the concentration tensors of its phases. This page uses the latter for two loadings that the stiffness alone does not cover: a temperature change, which makes each phase expand freely by its own amount, and a fluid pressure in the pores, which acts on the solid as a prestress. Both are answered by Levin's theorem, first numerically on a two-phase composite and a porous solid, then in closed form.

A two-phase composite

The composite is made of a compliant matrix with a large thermal expansion, in which stiff spherical inclusions with a small one are dispersed. The values are illustrative, with moduli in GPa and expansion coefficients in 1/K:

Phase
matrix :M	5	2
inclusions :I	40	30

using MeanFieldHomogenization
using TensND
using Plots
gr()  # headless backend; GKSwstype is set to "100" in make.jl

k_M, μ_M, α_M = 5.0, 2.0, 60.0e-6
k_I, μ_I, α_I = 40.0, 30.0, 5.0e-6
C_M, C_I = iso_stiffness(k_M, μ_M), iso_stiffness(k_I, μ_I)
𝟏 = TensISO{3}(1.0)   # second-order identity

function composite(f; shape = Ellipsoid(1.0))
    rve = RVE()
    add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => C_M); fraction = :rest)
    add_phase!(rve, :I, shape, Dict(:C => C_I); fraction = f)
    return rve
end

A temperature rise gives each phase the free strain    , and Levin's theorem gives the macroscopic one,    , from the stress concentration tensors of the scheme. Taking  , levin_free_gradient returns the effective thermal expansion tensor directly:

ct = concentration_tensors(composite(0.3), MoriTanaka(:M), :C)
α_hom = levin_free_gradient(ct, Dict(:M => α_M * 𝟏, :I => α_I * 𝟏))
α_hom[1, 1]

3.903314917127072e-5

For two phases the theorem needs no concentration tensor at all: the free strain is fixed by the effective compliance alone, and for isotropic phases by the effective bulk modulus,       . The two routes agree:

k_hom, _ = k_mu(ct.effective)
α_exact = α_M + (α_I - α_M) * (1 / k_hom - 1 / k_M) / (1 / k_I - 1 / k_M)
α_hom[1, 1] - α_exact

6.776263578034403e-21

The estimate depends on the scheme

Each scheme has its own concentration tensors, hence its own estimate of the thermal expansion. The bounds give two classical rules: under Reuss every phase carries the macroscopic stress,  , and the expansion coefficients mix by volume; under Voigt every phase undergoes the macroscopic strain, and they mix weighted by the bulk moduli,     .

function α_of(f, scheme)
    ct = concentration_tensors(composite(f), scheme, :C)
    return levin_free_gradient(ct, Dict(:M => α_M * 𝟏, :I => α_I * 𝟏))[1, 1]
end

fs = range(0.0, 0.6; length = 25)
schemes = [
    "Reuss" => Reuss(), "Voigt" => Voigt(), "Mori–Tanaka" => MoriTanaka(:M),
    "self-consistent" => SelfConsistent(), "differential" => DifferentialScheme(:M),
]
plt = plot(;
    xlabel = "volume fraction of inclusions f", ylabel = "α_hom (10⁻⁶/K)",
    title = "k_M = 5, μ_M = 2, α_M = 60e-6; k_I = 40, μ_I = 30, α_I = 5e-6",
    titlefontsize = 9, legend = :topright, framestyle = :box, size = (760, 480),
    left_margin = 5Plots.mm, bottom_margin = 5Plots.mm,
)
for (name, scheme) in schemes
    plot!(plt, fs, [1.0e6 * α_of(f, scheme) for f in fs]; label = name, lw = 2)
end
plt

Reuss and Voigt still frame the estimates, as they frame the bulk modulus from which, for two phases, the expansion follows. Of the estimates, Mori–Tanaka gives the largest expansion, the matrix being the continuous phase that imposes its own.

Aligned platelets: an anisotropic expansion

Flat inclusions aligned with the normal make the composite transversely isotropic, and its expansion with it. The two-phase formula still applies, with the effective compliance in place of the bulk modulus:        .

ct = concentration_tensors(composite(0.3; shape = Ellipsoid(1.0, 1.0, 0.1)), MoriTanaka(:M), :C)
α_hom = levin_free_gradient(ct, Dict(:M => α_M * 𝟏, :I => α_I * 𝟏))
α_exact = α_M * 𝟏 + (α_I - α_M) / (1 / k_I - 1 / k_M) * (3 * (𝟏 ⊡ ct.compliance) - 𝟏 / k_M)
(α₁₁ = α_hom[1, 1], α₃₃ = α_hom[3, 3], difference = maximum(abs, get_array(α_hom - α_exact)))

(α₁₁ = 2.3902698704529678e-5, α₃₃ = 6.128663515384499e-5, difference = 1.3552527156068805e-20)

In the plane of the platelets the stiff inclusions restrain the matrix, while across them the composite expands even more than the matrix alone: held back in the plane, the matrix expands further across it by the Poisson effect.

A saturated porous solid: the Biot tensor

A fluid at pressure filling the pores acts on them as the prestress  , and Levin's theorem in its second form,    , gives the macroscopic prestress  , where is the Biot tensor. The solid is the matrix above, with a porosity of 10 % in spherical pores and a family of penny-shaped cracks of normal and density 0.1, filled with the same fluid:

rve = RVE()
add_phase!(rve, :S, Ellipsoid(1.0), Dict(:C => C_M); fraction = :rest)
add_phase!(rve, :P, Ellipsoid(1.0), Dict(:C => zero(C_M)); fraction = 0.1)
add_phase!(rve, :c, PennyCrack(1.0), Dict(:C => C_M); density = 0.1)

ct = concentration_tensors(rve, MoriTanaka(:S), :C)
B = -levin_preflux(ct, Dict(:P => -𝟏, :c => -𝟏))
(B₁₁ = B[1, 1], B₃₃ = B[3, 3])

(B₁₁ = 0.3696059571184822, B₃₃ = 0.5265333599152829)

A crack family has no volume: it enters with its density as weight and its concentration tensor per unit density, and its openings across make the larger coefficient. With a single solid phase the theorem reduces to the classical      of biot_tensor:

maximum(abs, get_array(B - biot_tensor(ct.effective, C_M)))

2.7755575615628914e-16

With several solid phases that shortcut is lost and only the theorem applies; Poromechanics gives an example.

In closed form

Mori–Tanaka on spheres runs on symbolic moduli, and so do the concentration tensors and Levin's theorem: the two results above come out as formulas. With SymPy:

using MeanFieldHomogenization
using TensND
using SymPyPythonCall

@syms k0::positive μ0::positive k1::positive μ1::positive f::positive
@syms α0::positive α1::positive
𝟏 = TensISO{3}(Sym(1))

rve = RVE()
add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => iso_stiffness(k0, μ0)); fraction = :rest)
add_phase!(rve, :I, Ellipsoid(1.0), Dict(:C => iso_stiffness(k1, μ1)); fraction = f)
ct = concentration_tensors(rve, MoriTanaka(:M), :C)
α_MT = levin_free_gradient(ct, Dict(:M => α0 * 𝟏, :I => α1 * 𝟏))[1, 1]

The thermal expansion of the Mori–Tanaka estimate is

which does not depend on the shear modulus of the inclusions, and returns at   and at  :

simplify(α_MT - (α0 + (α1 - α0) * f * k1 * (3k0 + 4μ0) / (k0 * (3k1 + 4μ0) + 4μ0 * f * (k1 - k0))))

Turning the inclusions into pores filled with a fluid gives the Biot coefficient of the Mori–Tanaka estimate for spherical pores of porosity ,

which tends to for an incompressible solid, the limit of Terzaghi's effective stress, and to the porosity as  :

rve_p = RVE()
add_phase!(rve_p, :S, Ellipsoid(1.0), Dict(:C => iso_stiffness(k0, μ0)); fraction = :rest)
add_phase!(rve_p, :P, Ellipsoid(1.0), Dict(:C => zero(iso_stiffness(k0, μ0))); fraction = f)
ct_p = concentration_tensors(rve_p, MoriTanaka(:S), :C)
b_MT = -levin_preflux(ct_p, Dict(:P => -𝟏))[1, 1]
simplify(b_MT - f * (3k0 + 4μ0) / (3f * k0 + 4μ0))

Where to go next

The same schemes apply at tensor order 2, to a conductivity or a diffusivity: Transport properties opens that extension. The theorem and its proof are in Levin's theorem, and the concentration tensors of each scheme at the end of its theory page, from Voigt and Reuss bounds to The equivalent inclusion method. The tutorial Symbolic spheres derives the Mori–Tanaka moduli used above in closed form, and Poromechanics closes the poroelastic law with the Biot modulus.