The cluster model on cubic arrays
The whole page as one script
Every block of this page, in order, is scripts/91_cluster_cubic_arrays.jl. From a clone of the repository, julia --project=docs -e 'using Pkg; Pkg.instantiate()' once, then julia scripts/91_cluster_cubic_arrays.jl, runs it end to end.
Before this page
Two-inclusion interaction tensors, on which the scheme rests, The cluster model, where it is derived, and Particle assemblies and N-body schemes for the assemblies it acts on.
The cluster model of Molinari and El Mouden (1996) solves for the mean strain of every inclusion, accounting for the pairwise interaction with each neighbor inside a cluster of radius R_c. Unlike the one-site schemes it therefore sees where the inclusions are, which is why it acts on a ParticleAssembly rather than on an RVE.
This page reproduces the reference results of their §4 on simple-cubic arrays of stiff spheres: the convergence in cluster radius (their Fig. 3), the comparison with the other schemes of the package (their Fig. 5) and with Mori-Tanaka in particular (their Fig. 6).
The scheme rests on the two-inclusion interaction tensor
whose self counterpart is
using MeanFieldHomogenization
using TensND
using LinearAlgebra
using Printf
using Plots
gr()
default(; left_margin = 5Plots.mm, bottom_margin = 5Plots.mm)Matrix: k = 1, μ = 0.4 (hence ν = 0.3, the value used throughout the paper). "Rigid" spheres are four orders of magnitude stiffer, which is numerically indistinguishable from the rigid limit at these contrasts.
k_m, μ_m = 1.0, 0.4
C_m = TensISO{3}(3k_m, 2μ_m)
C_rigid = TensISO{3}(3 * 1.0e4, 2 * 1.0e4)
C_void = TensISO{3}(3 * 1.0e-9, 2 * 1.0e-9)
𝟏 = TensISO{3}(1.0)
μ_of(C) = canonical_components(C)[1, 2, 1, 2]
κ_of(C) = (𝟏 ⊡ C ⊡ 𝟏) / 9 # 𝟏 : ℂ : 𝟏 / 9κ_of (generic function with 1 method)The one-site schemes read the same assembly: RVE(asm) forgets the positions and keeps the derived fractions, so no parallel cell has to be built by hand for the comparisons below — and none can silently drift out of step with the assembly it is meant to mirror.
§1 Convergence in the cluster radius
Their Fig. 3. Below one period the cluster contains no neighbor at all and the scheme degenerates exactly onto Mori-Tanaka — the identity proved in their Appendix C. Beyond R_c ≈ 2 periods the estimate is flat; the residual wobble is the discrete shell structure of a spherical cutoff, visible in their own figure too.
cutoffs = 0.0:0.25:5.0
conv = Dict{Float64, Vector{Float64}}()
for f in (0.2, 0.5)
conv[f] = [
μ_of(
homogenize(
cubic_lattice(
:sc, Dict(:C => C_m), Dict(:C => C_rigid);
fraction = f, cutoff = c
), ClusterModel(), :C
)
) / μ_m for c in cutoffs
]
end
for f in (0.2, 0.5)
@printf "f = %.1f : Mori-Tanaka %.4f cluster(R_c=0) %.4f cluster(R_c=3) %.4f\n" f (μ_of(homogenize(cubic_lattice(:sc, Dict(:C => C_m), Dict(:C => C_rigid); fraction = f), MoriTanaka(), :C)) / μ_m) conv[f][1] conv[f][13]
end
p1 = plot(;
xlabel = "R_c / L", ylabel = "μ_eff / μ_m",
title = "Convergence in cluster radius (SC, rigid spheres)",
legend = :right
)
for f in (0.2, 0.5)
plot!(p1, cutoffs, conv[f]; marker = :circle, ms = 3, label = "cluster, f = $f")
hline!(
p1, [μ_of(homogenize(cubic_lattice(:sc, Dict(:C => C_m), Dict(:C => C_rigid); fraction = f), MoriTanaka(), :C)) / μ_m];
ls = :dash, label = "Mori-Tanaka, f = $f"
)
end
p1
§2 Against the other schemes of the package
Their Fig. 5 and Fig. 6. The one-site self-consistent scheme is unsuited to a matrix-inclusion morphology — it does not see the connectivity of the matrix — while the three-phase and differential schemes, which do, sit much closer to the cluster estimate. Mori-Tanaka overestimates the shear modulus here: for a simple-cubic (i.e. non-isotropic) distribution of stiff spheres it is not a bound.
fracs = 0.05:0.05:0.45
cluster = Float64[]
mt = Float64[]
sc = Float64[]
diff_ = Float64[]
for f in fracs
asm = cubic_lattice(
:sc, Dict(:C => C_m), Dict(:C => C_rigid); fraction = f, cutoff = 3.0
)
push!(cluster, μ_of(homogenize(asm, ClusterModel(), :C)) / μ_m)
push!(mt, μ_of(homogenize(asm, MoriTanaka(), :C)) / μ_m)
push!(sc, μ_of(homogenize(asm, SelfConsistent(), :C)) / μ_m)
push!(diff_, μ_of(homogenize(asm, DifferentialScheme(; nsteps = 200), :C)) / μ_m)
end
p2 = plot(
fracs, cluster;
marker = :circle, ms = 3, label = "cluster (R_c = 3L)",
xlabel = "f", ylabel = "μ_eff / μ_m",
title = "Simple-cubic array of rigid spheres", legend = :topleft
)
plot!(p2, fracs, mt; marker = :square, ms = 3, ls = :dash, label = "Mori-Tanaka")
plot!(p2, fracs, sc; marker = :diamond, ms = 3, ls = :dot, label = "self-consistent")
plot!(p2, fracs, diff_; marker = :utriangle, ms = 3, ls = :dashdot, label = "differential")
p2
§3 What the cluster changes, and what it cannot
The interaction tensor has a strictly vanishing isotropic part, so summing it over a cluster cannot alter the spherical part of the problem. The effective bulk modulus of a cubic array is therefore exactly the Mori-Tanaka one, at every volume fraction — only the shear response sees the arrangement.
println("\n f κ_cluster / κ_MT μ_cluster / μ_MT")
for f in (0.1, 0.2, 0.3, 0.4)
asm = cubic_lattice(:sc, Dict(:C => C_m), Dict(:C => C_rigid); fraction = f, cutoff = 3.0)
Ccl = homogenize(asm, ClusterModel(), :C)
Cmt = homogenize(asm, MoriTanaka(), :C) # same cell, positions ignored
@printf " %.1f %.14f %.6f\n" f (κ_of(Ccl) / κ_of(Cmt)) (μ_of(Ccl) / μ_of(Cmt))
end
f κ_cluster / κ_MT μ_cluster / μ_MT
0.1 1.00000000000000 0.984727
0.2 1.00000000000000 0.954691
0.3 1.00000000000000 0.923106
0.4 1.00000000000000 0.898740A cubic array is cubic, not isotropic: the two shear constants differ, and the gap is the anisotropy the one-site schemes cannot express.
println("\n f C_1212/μ_m (C_1111-C_1122)/2μ_m anisotropy ratio")
for f in (0.1, 0.2, 0.3, 0.4)
asm = cubic_lattice(:sc, Dict(:C => C_m), Dict(:C => C_rigid); fraction = f, cutoff = 3.0)
A = canonical_components(homogenize(asm, ClusterModel(), :C))
μ1 = A[1, 2, 1, 2] / μ_m
μ2 = (A[1, 1, 1, 1] - A[1, 1, 2, 2]) / (2μ_m)
@printf " %.1f %.5f %.5f %.4f\n" f μ1 μ2 (μ1 / μ2)
end
f C_1212/μ_m (C_1111-C_1122)/2μ_m anisotropy ratio
0.1 1.21772 1.27200 0.9573
0.2 1.46292 1.68688 0.8672
0.3 1.76552 2.28203 0.7737
0.4 2.17455 3.06598 0.7093§4 Porous arrays: the spatial distribution matters
Their Fig. 16. At equal porosity the three cubic arrangements give different shear moduli, the simple-cubic one being the softest — an effect no scheme driven by volume fractions alone can produce.
p3 = plot(;
xlabel = "porosity f", ylabel = "μ_eff / μ_m",
title = "Spherical voids: effect of the arrangement", legend = :topright
)
for (kind, mk) in ((:sc, :circle), (:bcc, :square), (:fcc, :diamond))
fmax = min(0.45, max_packing_fraction(kind) - 0.02)
fs = 0.05:0.05:fmax
ys = [
μ_of(
homogenize(
cubic_lattice(
kind, Dict(:C => C_m), Dict(:C => C_void); fraction = f, cutoff = 3.0
), ClusterModel(), :C
)
) / μ_m for f in fs
]
plot!(p3, fs, ys; marker = mk, ms = 3, label = uppercase(string(kind)))
end
fs = 0.05:0.05:0.45
plot!(
p3, fs, [
μ_of(
homogenize(
cubic_lattice(:sc, Dict(:C => C_m), Dict(:C => C_void); fraction = f),
MoriTanaka(), :C
)
) / μ_m for f in fs
];
ls = :dash, color = :black, label = "Mori-Tanaka"
)
p3
§5 The equivalent inclusion method gives the same answer
On a periodic assembly the two N-body schemes of the package are the same linear system — the identity Brisard et al. (2014) state in their §3.1, their order-zero influence pseudotensors being the interaction tensors of El Mouden. They agree to machine precision, which is the acceptance gate of both implementations.
println("\n cutoff cluster EIM ‖difference‖")
for c in (0.0, 1.5, 3.0)
asm = cubic_lattice(:sc, Dict(:C => C_m), Dict(:C => C_rigid); fraction = 0.3, cutoff = c)
A = homogenize(asm, ClusterModel(), :C)
B = homogenize(asm, EquivalentInclusion(), :C)
@printf " %.1f %.10f %.10f %.2e\n" c μ_of(A) μ_of(B) norm(A - B)
end
p_full = plot(p1, p2, p3; layout = (1, 3), left_margin = 8Plots.mm, bottom_margin = 8Plots.mm,
size = (1500, 520), titlefontsize = 9)
p_full
The last figure gathers the panels of §1, §2 and §4 in one row. Read with the bulk-modulus table of §3, they show the cluster estimate departing from Mori-Tanaka through the shear modulus alone, by an amount that depends on the arrangement of the inclusions.
Where to go next
The identity of §5 is the starting point of The equivalent inclusion method, against a published table, which leaves the periodic array for a random assembly of pores inside a circular statistical volume element and reproduces a published table of the method. The cubic estimate obtained here can itself feed a second scale, which requires the interaction tensor in an anisotropic reference: that is the subject of the tutorial Chaining scales through an N-body scheme.
This page was generated using Literate.jl.