The Walpole basis end to end

Everything the transversely isotropic machinery rests on, computed rather than quoted: the six as   matrices, the multiplication table, the Gram matrix, the correct decomposition of , the synthetic   product and inverse rules, and the three parametrizations of a transversely isotropic material compared on the same material.

Theory: The Walpole basis and TI parametrizations.

using TensND
using LinearAlgebra
using Printf

n = (0.0, 0.0, 1.0)          # symmetry axis
W = walpole_basis(n)
arr(x) = get_array(x)

arr (generic function with 1 method)

The six tensors, in Kelvin–Mandel form

Built from    and   , in the frame where  .

function show_KM(label, t)
    M = KM(t)
    println(label)
    for r in 1:6
        println("   ", join([abs(M[r, c]) < 1.0e-12 ? "   .   " : @sprintf("%7.4f", M[r, c]) for c in 1:6], " "))
    end
    return println()
end

for i in 1:6
    show_KM("𝕎$i =", W[i])
end

𝕎1 =
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .     1.0000    .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .

𝕎2 =
    0.5000  0.5000    .       .       .       .
    0.5000  0.5000    .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .

𝕎3 =
      .       .       .       .       .       .
      .       .       .       .       .       .
    0.7071  0.7071    .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .

𝕎4 =
      .       .     0.7071    .       .       .
      .       .     0.7071    .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .

𝕎5 =
    0.5000 -0.5000    .       .       .       .
   -0.5000  0.5000    .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .     1.0000

𝕎6 =
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .       .       .       .
      .       .       .     1.0000    .       .
      .       .       .       .     1.0000    .
      .       .       .       .       .       .

and are transposes of each other — the only two that are not major-symmetric individually:

norm(arr(W[3]) - permutedims(arr(W[4]), (3, 4, 1, 2)))

0.0

The multiplication table

The basis is closed under ⊡. Decomposing each product back onto the basis shows the structure: behave exactly as the matrix units of   matrices, while and are orthogonal idempotents.

gram = [sum(arr(W[i]) .* arr(W[j])) for i in 1:6, j in 1:6]
coeffs(P) = [sum(P .* arr(W[k])) / gram[k, k] for k in 1:6]

println("   𝕎ᵢ : 𝕎ⱼ")
println("      ", join([lpad("𝕎$j", 5) for j in 1:6]))
for i in 1:6
    row = String[]
    for j in 1:6
        P = arr(W[i]) ⊡ arr(W[j])
        if norm(P) < 1.0e-12
            push!(row, lpad("0", 5))
        else
            c = coeffs(P)
            k = argmax(abs.(c))
            push!(row, lpad("𝕎$k", 5))
        end
    end
    println("  𝕎$i  ", join(row))
end

   𝕎ᵢ : 𝕎ⱼ
         𝕎1   𝕎2   𝕎3   𝕎4   𝕎5   𝕎6
  𝕎1     𝕎1    0   𝕎3    0    0    0
  𝕎2      0   𝕎2    0   𝕎4    0    0
  𝕎3      0   𝕎3    0   𝕎1    0    0
  𝕎4     𝕎4    0   𝕎2    0    0    0
  𝕎5      0    0    0    0   𝕎5    0
  𝕎6      0    0    0    0    0   𝕎6

The Gram matrix is diagonal

  with  . This is the single fact that makes every TI projection closed-form: the normal equations decouple.

round.(gram, digits = 12)

6×6 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  0.0
 0.0  0.0  0.0  0.0  0.0  2.0

The isotropic tensors on this basis

Easy to get wrong, so worth computing. is not , and is not  .

𝕀, 𝕁, 𝕂 = ISO(Val(3), Val(Float64))
s2 = sqrt(2)

checks = [
    "𝕀 = 𝕎₁+𝕎₂+𝕎₅+𝕎₆" =>
        norm(arr(𝕀) - (arr(W[1]) + arr(W[2]) + arr(W[5]) + arr(W[6]))),
    "𝕁 = (𝕎₁+2𝕎₂+√2𝕎₃+√2𝕎₄)/3" =>
        norm(arr(𝕁) - (arr(W[1]) + 2arr(W[2]) + s2 * arr(W[3]) + s2 * arr(W[4])) / 3),
    "𝕂 = (2𝕎₁+𝕎₂−√2𝕎₃−√2𝕎₄)/3 + 𝕎₅ + 𝕎₆" =>
        norm(arr(𝕂) - ((2arr(W[1]) + arr(W[2]) - s2 * (arr(W[3]) + arr(W[4]))) / 3 + arr(W[5]) + arr(W[6]))),
    "— and the WRONG one: 𝕀 = Σᵢ𝕎ᵢ" =>
        norm(arr(𝕀) - sum(arr(W[i]) for i in 1:6)),
    "— and the WRONG one: 𝕁 = 𝕎₁+𝕎₂" =>
        norm(arr(𝕁) - (arr(W[1]) + arr(W[2]))),
]
for (name, r) in checks
    @printf "%-40s residual = %.3e\n" name r
end

𝕀 = 𝕎₁+𝕎₂+𝕎₅+𝕎₆                          residual = 0.000e+00
𝕁 = (𝕎₁+2𝕎₂+√2𝕎₃+√2𝕎₄)/3                 residual = 0.000e+00
𝕂 = (2𝕎₁+𝕎₂−√2𝕎₃−√2𝕎₄)/3 + 𝕎₅ + 𝕎₆       residual = 2.077e-16
— and the WRONG one: 𝕀 = Σᵢ𝕎ᵢ            residual = 1.414e+00
— and the WRONG one: 𝕁 = 𝕎₁+𝕎₂           residual = 1.000e+00

The last two are off by and respectively — not round-off.

Reading the Walpole coefficients of the three directly:

for (name, t) in (("𝕀", 𝕀), ("𝕁", 𝕁), ("𝕂", 𝕂))
    B, _, _ = proj_tens(:TI, arr(t), [0.0, 0.0, 1.0])
    println(rpad(name, 3), " → (ℓ₁,…,ℓ₆) = ", round.(get_ℓ(B), digits = 6))
end

𝕀   → (ℓ₁,…,ℓ₆) = (1.0, 1.0, 0.0, 0.0, 1.0, 1.0)
𝕁   → (ℓ₁,…,ℓ₆) = (0.333333, 0.666667, 0.471405, 0.471405, 0.0, 0.0)
𝕂   → (ℓ₁,…,ℓ₆) = (0.666667, 0.333333, -0.471405, -0.471405, 1.0, 1.0)

The synthetic algebra

With  , the double contraction is   and the inverse is .

A = TensTI{4}(2.0, 1.0, 0.5, 0.9, 1.3, n)      # major-symmetric, N=5
B = TensTI{4}(3.0, 1.5, 0.2, 1.1, 0.4, n)
P = A ⊡ B

println("typeof(A)   = ", typeof(A))
println("typeof(A⊡B) = ", typeof(P))

typeof(A)   = TensTI{4, Float64, 5}
typeof(A⊡B) = TensTI{4, Float64, 6}

The product of two major-symmetric TI tensors is generally not major-symmetric —     unless the blocks commute — so the storage widens from N=5 to N=6:

ℓA, ℓB, ℓP = get_ℓ(A), get_ℓ(B), get_ℓ(P)
LA = [ℓA[1] ℓA[3]; ℓA[4] ℓA[2]]
LB = [ℓB[1] ℓB[3]; ℓB[4] ℓB[2]]

println("L_A · L_B   = ", round.(LA * LB, digits = 10))
println("block of P  = ", round.([ℓP[1] ℓP[3]; ℓP[4] ℓP[2]], digits = 10))
println("ℓ₅ rule     : ", ℓA[5] * ℓB[5], " vs ", ℓP[5])
println("ℓ₆ rule     : ", ℓA[6] * ℓB[6], " vs ", ℓP[6])
println("ℓ₃ == ℓ₄ ?  ", isapprox(ℓP[3], ℓP[4]))

L_A · L_B   = [6.1 1.15; 1.7 1.6]
block of P  = [6.1 1.15; 1.7 1.6]
ℓ₅ rule     : 0.9900000000000001 vs 0.9900000000000001
ℓ₆ rule     : 0.52 vs 0.52
ℓ₃ == ℓ₄ ?  false

Against the dense 81-component product:

norm(arr(P) - arr(A) ⊡ arr(B))

3.597533769998862e-16

Inversion, likewise closed-form:

norm(arr(A ⊡ inv(A)) - arr(𝕀))

0.0

The symmetrized basis

Merging the pair,   , spans exactly the five-dimensional major-symmetric subspace, with Gram :

Ws = walpole_basis_sym(n)
[
    norm(arr(Ws[1]) - arr(W[1])),
    norm(arr(Ws[2]) - arr(W[2])),
    norm(arr(Ws[3]) - (arr(W[3]) + arr(W[4]))),
    norm(arr(Ws[4]) - arr(W[5])),
    norm(arr(Ws[5]) - arr(W[6])),
]

5-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0

[sum(arr(Ws[i]) .* arr(Ws[i])) for i in 1:5]

5-element Vector{Float64}:
 1.0
 1.0
 1.9999999999999996
 2.0
 2.0

Three parametrizations, one material

A unidirectional carbon/epoxy ply — a standard transversely isotropic composite, the fiber direction being the symmetry axis. Engineering constants first (moduli in GPa).

nv = [0.0, 0.0, 1.0]
E₁, E₃, ν₁₂, ν₃₁, G₃₁ = 9.0, 140.0, 0.40, 0.30, 4.6

𝕊 = tens_TI_eng(E₁, E₃, ν₁₂, ν₃₁, G₃₁, nv)     # compliance
ℂ = inv(𝕊)                                       # stiffness

println("engineering → back      : ", round.(arg_TI_eng(𝕊), digits = 8))
println("‖ℂ:𝕊 − 𝕀‖               : ", norm(arr(ℂ ⊡ 𝕊) - arr(𝕀)))

engineering → back      : (9.0, 140.0, 0.4, 0.3, 4.6)
‖ℂ:𝕊 − 𝕀‖               : 5.551115123125782e-17

The same compliance read in component form, :

round.(arg_TI(𝕊), sigdigits = 6)

(0.111111, -0.0444444, -0.00214286, 0.00714286, 0.0543478)

and in Hoenig form, the dimensionless ratios:

Eh, ν₁, ν₂, H, Γ = arg_TI_Hoenig(𝕊)
@printf "E  = %8.4f GPa\nν₁ = %8.4f\nν₂ = %8.4f\nH  = %8.4f   (axial/transverse modulus ratio)\nΓ  = %8.4f   (shear anisotropy)\n" Eh ν₁ ν₂ H Γ

E  =   9.0000 GPa
ν₁ =   0.4000
ν₂ =   0.0193
H  =  15.5556   (axial/transverse modulus ratio)
Γ  =   1.4311   (shear anisotropy)

and quantify the anisotropy directly and dimensionlessly: here  , so the fiber direction is some fifteen times stiffer than the transverse plane, while   says the axial shear modulus departs from its in-plane counterpart by rather less. Rebuilding from the Hoenig parameters returns the same tensor:

norm(arr(tens_TI_Hoenig(Eh, ν₁, ν₂, H, Γ, nv)) - arr(𝕊))

1.962615573354719e-17

An isotropic material sits at the point    with   — the Hoenig parameters measure departure from isotropy:

𝕊iso = inv(3 * 20.0 * 𝕁 + 2 * 8.0 * 𝕂)
Eᵢ, ν₁ᵢ, ν₂ᵢ, Hᵢ, Γᵢ = arg_TI_Hoenig(fromISO(𝕊iso, nv))
@printf "E = %.6f   ν₁ = %.6f   ν₂ = %.6f   H = %.10f   Γ = %.10f\n" Eᵢ ν₁ᵢ ν₂ᵢ Hᵢ Γᵢ

E = 21.176471   ν₁ = 0.323529   ν₂ = 0.323529   H = 1.0000000000   Γ = 1.0000000000

The Walpole coefficients of the ply

Stiffness and compliance are inverse in the synthetic algebra, so their   blocks are inverse matrices and their are reciprocal:

ℓC, ℓS = get_ℓ(ℂ), get_ℓ(𝕊)
LC = [ℓC[1] ℓC[3]; ℓC[4] ℓC[2]]
LS = [ℓS[1] ℓS[3]; ℓS[4] ℓS[2]]

println("L_ℂ · L_𝕊 = ", round.(LC * LS, digits = 12))
println("ℓ₅ᶜ·ℓ₅ˢ   = ", round(ℓC[5] * ℓS[5], digits = 12))
println("ℓ₆ᶜ·ℓ₆ˢ   = ", round(ℓC[6] * ℓS[6], digits = 12))

L_ℂ · L_𝕊 = [1.0 0.0; 0.0 1.0]
ℓ₅ᶜ·ℓ₅ˢ   = 1.0
ℓ₆ᶜ·ℓ₆ˢ   = 1.0

---

This page was generated using Literate.jl.