Hill polarization tensors in practice
The Hill tensor
The closed forms are collected in Hill polarization tensors; this page is the computational tour. Related scripts go further on individual points: 01_auxiliary_tensors.jl for the geometric tensors 03_hill_conductivity.jl for the 2nd-order transport counterpart, 06_cylinder.jl for the infinite-cylinder limit and 07_hill_ti_coaxial.jl for the transversely-isotropic closed form.
using MeanFieldHomogenization
using TensND
using LinearAlgebra
using Printf§0 A reference matrix
Steel,
const E_ref = 210.0e3
const ν_ref = 0.3
const λ_ref = E_ref * ν_ref / ((1 + ν_ref) * (1 - 2ν_ref))
const μ_ref = E_ref / (2 * (1 + ν_ref))
const k_ref = λ_ref + 2μ_ref / 3
const C_iso = TensISO{3}(3k_ref, 2μ_ref)
@printf "λ = %.2f, μ = %.2f, k = %.2f MPa\n" λ_ref μ_ref k_refλ = 121153.85, μ = 80769.23, k = 175000.00 MPaTwo display helpers — Voigt storage for printing, Mandel for the linear algebra of §5. Double contraction itself is TensND's ⊡, used as P ⊡ C_iso in §4; there is no need to write index loops for it.
const voigt_idx = ((1, 1), (2, 2), (3, 3), (2, 3), (3, 1), (1, 2))
const voigt_lab = ["11", "22", "33", "23", "31", "12"]
function print_voigt(C; label = "P", scale = 1.0e6, unit = "×10⁻⁶ MPa⁻¹")
println(" Voigt[$label] ($unit):")
print(" ")
for l in voigt_lab
@printf "%10s" l
end
println()
for (I, (i, j)) in enumerate(voigt_idx)
@printf " %2s | " voigt_lab[I]
for (k, l) in voigt_idx
@printf "%9.4f " scale * C[i, j, k, l]
end
println()
end
return
end
const mandel_f = (1.0, 1.0, 1.0, √2, √2, √2)
mandel(C) = [
C[i, j, k, l] * mandel_f[I] * mandel_f[J]
for (I, (i, j)) in enumerate(voigt_idx), (J, (k, l)) in enumerate(voigt_idx)
]mandel (generic function with 1 method)§1 Isotropic matrix, four geometries
One call, hill_tensor, for every shape. The algorithm is selected automatically from the symmetry of the matrix and the shape of the inclusion.
The sphere
The only case with a one-line closed form:
let ell = Ellipsoid(1.0)
P = hill_tensor(ell, C_iso)
P1111_theory = 1 / (5 * (λ_ref + 2μ_ref)) + (1 / 3 - 1 / 5) / μ_ref
@printf " P[1,1,1,1] = %12.9e MPa⁻¹\n" P[1, 1, 1, 1]
@printf " theory = %12.9e MPa⁻¹ (err = %.2e)\n" P1111_theory abs(P[1, 1, 1, 1] - P1111_theory)
@printf " isotropy : |P[1111] − P[2222]| = %.2e\n\n" abs(P[1, 1, 1, 1] - P[2, 2, 2, 2])
print_voigt(P)
end P[1,1,1,1] = 2.358276644e-06 MPa⁻¹
theory = 2.358276644e-06 MPa⁻¹ (err = 0.00e+00)
isotropy : |P[1111] − P[2222]| = 0.00e+00
Voigt[P] (×10⁻⁶ MPa⁻¹):
11 22 33 23 31 12
11 | 2.3583 -0.5896 -0.5896 0.0000 0.0000 0.0000
22 | -0.5896 2.3583 -0.5896 0.0000 0.0000 0.0000
33 | -0.5896 -0.5896 2.3583 0.0000 0.0000 0.0000
23 | 0.0000 0.0000 0.0000 1.4739 0.0000 0.0000
31 | 0.0000 0.0000 0.0000 0.0000 1.4739 0.0000
12 | 0.0000 0.0000 0.0000 0.0000 0.0000 1.4739The prolate spheroid — a long fiber
a = 5, b = c = 1. The answer must be transversely isotropic about the long axis, which is a free correctness check.
let ell = Ellipsoid(5.0, 1.0, 1.0)
P = hill_tensor(ell, C_iso)
@printf " P[1,1,1,1] = %12.9e MPa⁻¹ (axial)\n" P[1, 1, 1, 1]
@printf " P[2,2,2,2] = %12.9e MPa⁻¹ (transverse)\n" P[2, 2, 2, 2]
@printf " transverse isotropy: |P[2222] − P[3333]| = %.2e\n" abs(P[2, 2, 2, 2] - P[3, 3, 3, 3])
end P[1,1,1,1] = 5.377298334e-07 MPa⁻¹ (axial)
P[2,2,2,2] = 2.841312301e-06 MPa⁻¹ (transverse)
transverse isotropy: |P[2222] − P[3333]| = 0.00e+00The oblate spheroid — a disc
a = b = 5, c = 1. The component normal to the disc dominates, which is the whole reason flat inclusions are such effective reinforcements — and, in the
let ell = Ellipsoid(5.0, 5.0, 1.0)
P = hill_tensor(ell, C_iso)
@printf " P[1,1,1,1] = %12.9e MPa⁻¹ (in-plane)\n" P[1, 1, 1, 1]
@printf " P[3,3,3,3] = %12.9e MPa⁻¹ (normal — dominant)\n" P[3, 3, 3, 3]
@printf " ratio P[3333]/P[1111] = %.2f\n" P[3, 3, 3, 3] / P[1, 1, 1, 1]
end P[1,1,1,1] = 1.044422018e-06 MPa⁻¹ (in-plane)
P[3,3,3,3] = 3.527507529e-06 MPa⁻¹ (normal — dominant)
ratio P[3333]/P[1111] = 3.38The triaxial ellipsoid
No symmetry left to exploit, so this is where the minor and major symmetries of
let ell = Ellipsoid(4.0, 2.0, 1.0)
P = hill_tensor(ell, C_iso)
print_voigt(P; label = "P triaxial")
err_min = maximum(abs(P[i, j, k, l] - P[j, i, k, l]) for i in 1:3, j in 1:3, k in 1:3, l in 1:3)
err_maj = maximum(abs(P[i, j, k, l] - P[k, l, i, j]) for i in 1:3, j in 1:3, k in 1:3, l in 1:3)
@printf "\n symmetry: minor %.2e, major %.2e\n" err_min err_maj
end Voigt[P triaxial] (×10⁻⁶ MPa⁻¹):
11 22 33 23 31 12
11 | 0.9923 -0.2426 -0.3522 0.0000 0.0000 0.0000
22 | -0.2426 2.0404 -0.7904 0.0000 0.0000 0.0000
33 | -0.3522 -0.7904 3.2752 0.0000 0.0000 0.0000
23 | 0.0000 0.0000 0.0000 1.9571 0.0000 0.0000
31 | 0.0000 0.0000 0.0000 0.0000 1.8616 0.0000
12 | 0.0000 0.0000 0.0000 0.0000 0.0000 0.9866
symmetry: minor 0.00e+00, major 0.00e+00§2 Anisotropic matrix — two algorithms, one answer
With an anisotropic MeanFieldHomogenization offers several independent routes to it — a residue evaluation of the acoustic-tensor poles, and two cubatures (:nestedquadgk, built in, and :decuhr when the DECUHR extension is loaded). Agreement between two of them is the strongest available check, since they share no code.
The matrix here is cubic — C11 = 250, C12 = 100, C44 = 80 GPa — with a Zener ratio
let
C11, C12, C44 = 250.0e3, 100.0e3, 80.0e3
C_arr = zeros(3, 3, 3, 3)
for (I, (i, j)) in enumerate(voigt_idx), (J, (k, l)) in enumerate(voigt_idx)
v = [
C11 C12 C12 0 0 0;
C12 C11 C12 0 0 0;
C12 C12 C11 0 0 0;
0 0 0 C44 0 0;
0 0 0 0 C44 0;
0 0 0 0 0 C44
][I, J]
C_arr[i, j, k, l] = C_arr[j, i, k, l] = C_arr[i, j, l, k] = C_arr[j, i, l, k] = v
end
C_cubic = Tens(C_arr)
@printf "Zener ratio A = 2C44/(C11−C12) = %.4f (isotropic ⇔ 1)\n" 2C44 / (C11 - C12)
let _ell = Ellipsoid(2.0, 1.0, 1.0) # warm-up, so the
hill_tensor(_ell, C_cubic; method = :residues) # timings below do
hill_tensor(_ell, C_cubic; method = :nestedquadgk) # not count compilation
end
ell = Ellipsoid(3.0, 1.0, 1.0)
t_res = @elapsed P_res = hill_tensor(ell, C_cubic; method = :residues, abstol = 1.0e-8, reltol = 1.0e-6)
t_qgk = @elapsed P_qgk = hill_tensor(ell, C_cubic; method = :nestedquadgk, abstol = 1.0e-8, reltol = 1.0e-6)
max_err = maximum(abs(P_res[i, j, k, l] - P_qgk[i, j, k, l]) for i in 1:3, j in 1:3, k in 1:3, l in 1:3)
@printf " :residues %.4f s\n :nestedquadgk %.4f s\n" t_res t_qgk
@printf " max |P_residues − P_nestedquadgk| = %.3e MPa⁻¹\n" max_err
endZener ratio A = 2C44/(C11−C12) = 1.0667 (isotropic ⇔ 1)
:residues 0.0158 s
:nestedquadgk 0.0138 s
max |P_residues − P_nestedquadgk| = 6.036e-10 MPa⁻¹method = :auto does not pick the residues
For an anisotropic 3-D elastic reference the automatic choice is a cubature, not the residue algorithm: the residue acoustic polynomial degenerates when the reference is anisotropic in type and isotropic in value, which is exactly what the self-consistent and differential schemes feed back at their first step. method = :residues remains available explicitly, as above, and :decuhr is preferred over :nestedquadgk whenever import DECUHR, Integrals has been run.
§3 Two dimensions — plane strain
The same entry point, with 2-D inclusions and a 2-D reference.
let
C_iso2 = TensISO{2}(3k_ref, 2μ_ref)
P = hill_tensor(Ellipsoid(1.0; dim = 2), C_iso2)
@printf "circle : P[1111] = %.9e, isotropy err = %.2e\n" P[1, 1, 1, 1] abs(P[1, 1, 1, 1] - P[2, 2, 2, 2])
P = hill_tensor(Ellipsoid(4.0, 1.0), C_iso2)
@printf "ellipse : P[1111] = %.9e (major), P[2222] = %.9e (minor)\n" P[1, 1, 1, 1] P[2, 2, 2, 2]
E1, E2, ν12, G12 = 100.0e3, 200.0e3, 0.3, 40.0e3
ν21 = ν12 * E2 / E1
D = 1 - ν12 * ν21
C2_arr = zeros(2, 2, 2, 2)
C2_arr[1, 1, 1, 1] = E1 / D
C2_arr[2, 2, 2, 2] = E2 / D
for idx in ((1, 1, 2, 2), (2, 2, 1, 1))
C2_arr[idx...] = ν12 * E2 / D
end
for idx in ((1, 2, 1, 2), (1, 2, 2, 1), (2, 1, 1, 2), (2, 1, 2, 1))
C2_arr[idx...] = G12
end
P = hill_tensor(Ellipsoid(4.0, 1.0), Tens(C2_arr); abstol = 1.0e-8, reltol = 1.0e-6)
@printf "ellipse in an orthotropic matrix: P[1111] = %.9e, P[1212] = %.9e\n" P[1, 1, 1, 1] P[1, 2, 1, 2]
endcircle : P[1111] = 2.640056022e-06, isotropy err = 0.00e+00
ellipse : P[1111] = 1.340056022e-06 (major), P[2222] = 3.087955182e-06 (minor)
ellipse in an orthotropic matrix: P[1111] = 3.062568944e-06, P[1212] = 4.818880744e-06§4 The Eshelby tensor
let ell = Ellipsoid(1.0)
S = hill_tensor(ell, C_iso) ⊡ C_iso
S1111_th = (7 - 5ν_ref) / (15 * (1 - ν_ref))
S1122_th = (5ν_ref - 1) / (15 * (1 - ν_ref))
S1212_th = (4 - 5ν_ref) / (15 * (1 - ν_ref))
println(" component computed analytical error")
@printf " S[1,1,1,1] %12.8f %12.8f %.2e\n" S[1, 1, 1, 1] S1111_th abs(S[1, 1, 1, 1] - S1111_th)
@printf " S[1,1,2,2] %12.8f %12.8f %.2e\n" S[1, 1, 2, 2] S1122_th abs(S[1, 1, 2, 2] - S1122_th)
@printf " S[1,2,1,2] %12.8f %12.8f %.2e\n" S[1, 2, 1, 2] S1212_th abs(S[1, 2, 1, 2] - S1212_th)
S2 = hill_tensor(Ellipsoid(5.0, 1.0, 1.0), C_iso) ⊡ C_iso
@printf "\n prolate a/b = 5 : S[1111] = %.7f (axial), S[2222] = %.7f (transverse)\n" S2[1, 1, 1, 1] S2[2, 2, 2, 2]
end component computed analytical error
S[1,1,1,1] 0.52380952 0.52380952 1.11e-16
S[1,1,2,2] 0.04761905 0.04761905 3.47e-17
S[1,2,1,2] 0.23809524 0.23809524 2.78e-17
prolate a/b = 5 : S[1111] = 0.1107873 (axial), S[2222] = 0.6613053 (transverse)§5 From
The dilute estimate is the shortest path from a Hill tensor to a homogenized property:
For voids
let
f = 0.05
P = hill_tensor(Ellipsoid(1.0), C_iso)
M₀ = mandel(C_iso)
Mδ = -M₀
MP = mandel(P)
I6 = Matrix{Float64}(I, 6, 6)
M_A = inv(I6 + MP * Mδ)
M_eff = M₀ + f * Mδ * M_A
S_eff = inv(M_eff)
E1_eff = 1 / S_eff[1, 1]
μ_eff = 1 / (2 * S_eff[6, 6]) # Mandel: the shear block carries the factor 2
k₀, μ₀ = k_ref, μ_ref
μ_eff_th = μ₀ * (1 - 5f * (3k₀ + 4μ₀) / (9k₀ + 8μ₀))
@printf " f = %.2f\n" f
@printf " E_eff = %.1f MPa (E_eff/E₀ = %.4f)\n" E1_eff E1_eff / E_ref
@printf " μ_eff = %.1f MPa (μ_eff/μ₀ = %.4f)\n" μ_eff μ_eff / μ_ref
@printf " analytical dilute μ_eff = %.1f MPa (err = %.2e)\n" μ_eff_th abs(μ_eff - μ_eff_th)
end f = 0.05
E_eff = 188916.7 MPa (E_eff/E₀ = 0.8996)
μ_eff = 73059.4 MPa (μ_eff/μ₀ = 0.9045)
analytical dilute μ_eff = 73059.4 MPa (err = 1.46e-11)In practice none of that is written by hand: an RVE plus Dilute does the same algebra, in any symmetry class and for any inclusion family, and returns a tensor rather than a 6×6 array.
let
f = 0.05
rve = RVE()
add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => C_iso); fraction = :rest)
add_phase!(rve, :PORE, Ellipsoid(1.0), Dict(:C => TensISO{3}(0.0, 0.0)); fraction = f)
k_eff, μ_eff = k_mu(homogenize(rve, Dilute(), :C))
@printf " homogenize(rve, Dilute(), :C): k_eff = %.1f MPa, μ_eff = %.1f MPa\n" k_eff μ_eff
end homogenize(rve, Dilute(), :C): k_eff = 152031.2 MPa, μ_eff = 73059.4 MPaThis page was generated using Literate.jl.