mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-17 17:22:10 +00:00
chore(test): drop THM real-physics DOF regression
Remove stale coupling tests superseded by domain-level kernels. - Delete `test/dofs/test_thm_real_physics.jl`.
This commit is contained in:
Jukka Aho
@@ -1,705 +0,0 @@
|
||||
"""
|
||||
🚀 THE ULTIMATE REAL PHYSICS: Thermo-Hydro-Mechanical-Electric Coupling
|
||||
|
||||
This implements COMPLETE REAL PHYSICS for THM-E with ALL coupling terms!
|
||||
|
||||
Field Variables:
|
||||
- T: Temperature (Float64) at VERTICES - continuous H¹ field
|
||||
- u: Displacement (Vec{3}) at VERTICES - continuous H¹ vector field
|
||||
- p: Pore pressure (Float64) at CELLS - discontinuous L² field
|
||||
- φ: Electric potential (Float64) at EDGES - H(curl) field
|
||||
|
||||
═══════════════════════════════════════════════════════════════════════
|
||||
COMPLETE PHYSICS FORMULATION - FULLY COUPLED THM-E SYSTEM
|
||||
═══════════════════════════════════════════════════════════════════════
|
||||
|
||||
1️⃣ THERMAL (Heat Equation with Thermoelastic Coupling):
|
||||
ρcₚ ∂T/∂t - ∇·(κ∇T) = Q + α_T·T₀·E/(1-2ν) ∇·∂u/∂t + β_T·∂p/∂t + S·∇·J
|
||||
|
||||
where:
|
||||
- κ: thermal conductivity [W/(m·K)]
|
||||
- α_T: thermal expansion coefficient [1/K]
|
||||
- T₀: reference temperature [K]
|
||||
- β_T: thermal pressurization coefficient [K/Pa]
|
||||
- S: Seebeck coefficient [V/K]
|
||||
- J: electric current density [A/m²]
|
||||
|
||||
2️⃣ MECHANICAL (Linear Elasticity with Multi-Physics Coupling):
|
||||
ρ ∂²u/∂t² - ∇·σ = f
|
||||
|
||||
where constitutive law includes ALL couplings:
|
||||
σ = C : ε(u) - α_T·(T-T₀)·I - α_p·p·I - e^T·E
|
||||
|
||||
Strain: ε(u) = ½(∇u + ∇uᵀ)
|
||||
Elasticity: C_ijkl = λδ_ij δ_kl + μ(δ_ik δ_jl + δ_il δ_jk)
|
||||
|
||||
Coupling terms:
|
||||
- Thermal stress: α_T·E/(1-2ν)·(T-T₀)·I
|
||||
- Pore pressure: α_p·p·I (Biot coupling)
|
||||
- Piezoelectric: e_kij·E_k (converse piezoelectric effect)
|
||||
|
||||
3️⃣ HYDRAULIC (Darcy Flow with Biot and Thermal Coupling):
|
||||
S_s ∂p/∂t + α_p ∂(∇·u)/∂t + β_T ∂T/∂t - ∇·(k/μ_f ∇p) = q - ζ·∇·J
|
||||
|
||||
where:
|
||||
- k: permeability [m²]
|
||||
- μ_f: fluid viscosity [Pa·s]
|
||||
- S_s: specific storage [1/Pa]
|
||||
- α_p: Biot coefficient [-]
|
||||
- ζ: electro-osmotic coefficient [m²/(V·s)]
|
||||
|
||||
4️⃣ ELECTRIC (Charge Conservation with Multi-Physics Sources):
|
||||
∇·D = ρ_e
|
||||
∇×E = 0 ⟹ E = -∇φ
|
||||
|
||||
where constitutive law:
|
||||
D = ε·E + e:ε(u) - p·∇ζ
|
||||
J = σ_e·E + S·(-κ∇T)
|
||||
|
||||
- D: electric displacement [C/m²]
|
||||
- E: electric field [V/m]
|
||||
- ε: permittivity [F/m]
|
||||
- σ_e: electric conductivity [S/m]
|
||||
- e_kij: piezoelectric tensor (3rd order) [C/m²]
|
||||
|
||||
═══════════════════════════════════════════════════════════════════════
|
||||
COUPLING MATRIX (12 OFF-DIAGONAL BLOCKS):
|
||||
═══════════════════════════════════════════════════════════════════════
|
||||
|
||||
│ T u p φ
|
||||
─────┼──────────────────────────────────────────
|
||||
T │ K_TT K_Tu K_Tp K_Tφ
|
||||
│ (α_T) (β_T) (S)
|
||||
─────┼──────────────────────────────────────────
|
||||
u │ K_uT K_uu K_up K_uφ
|
||||
│ (α_T) (α_p) (e_kij)
|
||||
─────┼──────────────────────────────────────────
|
||||
p │ K_pT K_pu K_pp K_pφ
|
||||
│ (β_T) (α_p) (ζ)
|
||||
─────┼──────────────────────────────────────────
|
||||
φ │ K_φT K_φu K_φp K_φφ
|
||||
│ (S) (e_kij) (ζ)
|
||||
|
||||
Onsager reciprocity: K_ab = K_ba^T for all coupling pairs!
|
||||
|
||||
═══════════════════════════════════════════════════════════════════════
|
||||
APPLICATION DOMAINS:
|
||||
═══════════════════════════════════════════════════════════════════════
|
||||
- Geothermal energy extraction (T-H-M)
|
||||
- Nuclear waste repositories (T-H-M)
|
||||
- CO₂ geological sequestration (H-M)
|
||||
- Electrokinetic soil remediation (E-H-M)
|
||||
- Piezoelectric sensors/actuators (E-M)
|
||||
- Thermoelectric energy harvesting (T-E)
|
||||
- Smart materials (all coupled)
|
||||
|
||||
Use case: PROVING that JuliaFEM handles arbitrarily complex physics elegantly!
|
||||
"""
|
||||
|
||||
using JuliaFEM
|
||||
using Test
|
||||
using Tensors
|
||||
using LinearAlgebra
|
||||
using SparseArrays
|
||||
using Printf
|
||||
|
||||
@testset "🚀 REAL THM-E: Complete Physics on All Entity Types" begin
|
||||
println("\n" * "="^70)
|
||||
println("🚀 REAL THM-E: COMPLETE PHYSICS ON ALL ENTITY TYPES")
|
||||
println("="^70)
|
||||
|
||||
# Create 3D mesh: Two tetrahedra
|
||||
nodes = [
|
||||
Vec{3,Float64}((0.0, 0.0, 0.0)), # Node 1
|
||||
Vec{3,Float64}((1.0, 0.0, 0.0)), # Node 2
|
||||
Vec{3,Float64}((0.5, 1.0, 0.0)), # Node 3
|
||||
Vec{3,Float64}((0.5, 0.5, 1.0)), # Node 4
|
||||
Vec{3,Float64}((1.5, 0.5, 0.5)), # Node 5
|
||||
]
|
||||
connectivity = [
|
||||
(UInt32(1), UInt32(2), UInt32(3), UInt32(4)), # Tet 1
|
||||
(UInt32(2), UInt32(3), UInt32(4), UInt32(5)), # Tet 2
|
||||
]
|
||||
mesh = Mesh{Tetrahedron{4}}(nodes, connectivity)
|
||||
|
||||
println("\n3D Mesh: 2 tetrahedra, 5 nodes")
|
||||
|
||||
# Create ONE element type with ALL FOUR physics fields!
|
||||
println("\nCreating multi-field elements with ALL physics...")
|
||||
|
||||
# Define field spec as a TYPE using @DOFSet (hides NamedTuple implementation)
|
||||
S = @DOFSet{T::DOF{Temperature, Vertex},
|
||||
u::DOF{Displacement{3}, Vertex},
|
||||
p::DOF{Pressure, Cell},
|
||||
φ::DOF{ElectricPotential, Edge}}
|
||||
|
||||
# Step 1: Initialize DOF manager
|
||||
dof_mgr = DOFManager(mesh)
|
||||
|
||||
# Step 2: Register fields and create elements
|
||||
register_fields!(dof_mgr, S)
|
||||
elements = create_elements!(dof_mgr, Element{Tetrahedron{4}, Lagrange{1}, S})
|
||||
|
||||
n_total = dof_mgr.total_dofs
|
||||
|
||||
# Count DOFs by field (from first element structure)
|
||||
elem1 = first(elements)
|
||||
n_T = length(elem1.dof_indices.T)
|
||||
n_u = length(elem1.dof_indices.u)
|
||||
n_p = length(elem1.dof_indices.p)
|
||||
n_φ = length(elem1.dof_indices.φ)
|
||||
|
||||
# Total system DOFs (calculated from DOF manager!)
|
||||
n_T_total = count_field_dofs(dof_mgr, :T)
|
||||
n_u_total = count_field_dofs(dof_mgr, :u)
|
||||
n_p_total = count_field_dofs(dof_mgr, :p)
|
||||
n_φ_total = count_field_dofs(dof_mgr, :φ)
|
||||
|
||||
println(" Temperature: $n_T DOFs per element (total: $n_T_total in system)")
|
||||
println(" Displacement: $n_u DOFs per element (total: $n_u_total in system)")
|
||||
println(" Pressure: $n_p DOFs per element (total: $n_p_total in system)")
|
||||
println(" Electric: $n_φ DOFs per element (total: $n_φ_total in system)")
|
||||
println(" TOTAL SYSTEM DOFs: $n_total")
|
||||
|
||||
@test n_T == 4
|
||||
@test n_u == 12
|
||||
@test n_p == 1
|
||||
@test n_φ == 6
|
||||
# Note: Total may be less than sum due to shared DOFs between elements
|
||||
@test n_total > 0 && n_total ≤ n_T_total + n_u_total + n_p_total + n_φ_total
|
||||
|
||||
# Material parameters (scaled for numerical stability)
|
||||
κ = 1.0 # Thermal conductivity
|
||||
E = 10.0 # Young's modulus (reduced for better conditioning)
|
||||
ν = 0.25 # Poisson's ratio (avoid near-incompressibility)
|
||||
k_perm = 1.0 # Hydraulic permeability
|
||||
σ_e = 1.0 # Electric conductivity
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("ASSEMBLING REAL PHYSICS FROM MULTI-FIELD ELEMENTS (NO MOCKS!)")
|
||||
println("="^70)
|
||||
|
||||
# ONE global system for ALL fields (this is the whole point!)
|
||||
K = spzeros(Float64, n_total, n_total)
|
||||
F = zeros(Float64, n_total)
|
||||
|
||||
println("\n🔥 ONE ELEMENT LOOP - ALL PHYSICS!")
|
||||
println("="^70)
|
||||
|
||||
# ONE LOOP over elements - assemble ALL physics!
|
||||
for (elem_idx, elem) in enumerate(elements)
|
||||
println("\n📦 Element $elem_idx:")
|
||||
|
||||
# Get LOCAL-GLOBAL mapping for coupled assembly
|
||||
n_local = local_dof_count(elem) # Total local DOFs (ALL fields)
|
||||
dof_map = local_to_global_map(elem) # Local → Global mapping
|
||||
|
||||
# Get LOCAL DOF ranges for each field
|
||||
T_local = field_dof_range(elem, :T) # e.g., 1:4
|
||||
u_local = field_dof_range(elem, :u) # e.g., 5:16
|
||||
p_local = field_dof_range(elem, :p) # e.g., 17:17
|
||||
φ_local = field_dof_range(elem, :φ) # e.g., 18:23
|
||||
|
||||
println(" Total local DOFs: $n_local")
|
||||
println(" T local range: $T_local ($(length(T_local)) DOFs)")
|
||||
println(" u local range: $u_local ($(length(u_local)) DOFs)")
|
||||
println(" p local range: $p_local ($(length(p_local)) DOFs)")
|
||||
println(" φ local range: $φ_local ($(length(φ_local)) DOFs)")
|
||||
|
||||
# 🎯 BUILD ONE LOCAL COUPLED MATRIX (THIS IS THE BEEF!)
|
||||
K_local = zeros(n_local, n_local)
|
||||
F_local = zeros(n_local)
|
||||
|
||||
# Get geometry as tuple of Vec{3} (zero-allocation)
|
||||
conn = mesh.connectivity[elem_idx]
|
||||
X_nodes = ntuple(i -> nodes[conn[i]], 4) # NTuple{4, Vec{3}}
|
||||
|
||||
# Get integration points for Tet4 with linear basis (Gauss{1} = 1 point)
|
||||
ips = integration_points(Gauss{1}(), Tetrahedron{4}())
|
||||
ip = ips[1] # Single integration point at centroid
|
||||
ξ = ip.ξ
|
||||
weight = ip.weight
|
||||
|
||||
# Get basis function derivatives w.r.t. parametric coords (returns NTuple{4, Vec{3}})
|
||||
dN_dξ = get_basis_derivatives(Tetrahedron{4}(), Lagrange{1}(), ξ)
|
||||
|
||||
# Compute Jacobian: J = ∑ᵢ Xᵢ ⊗ (∂Nᵢ/∂ξ) - zero allocation with Tensors.jl!
|
||||
J = X_nodes[1] ⊗ dN_dξ[1]
|
||||
@inbounds for i in 2:4
|
||||
J += X_nodes[i] ⊗ dN_dξ[i]
|
||||
end
|
||||
|
||||
# Physical gradients: ∇N = J⁻ᵀ ⋅ (∂N/∂ξ)
|
||||
J_inv_T = transpose(inv(J))
|
||||
∇N = ntuple(i -> J_inv_T ⋅ dN_dξ[i], 4) # NTuple{4, Vec{3}}
|
||||
|
||||
# Volume = det(J) * weight (for reference element)
|
||||
volume = det(J) * weight
|
||||
|
||||
# ================================================================
|
||||
# 1️⃣ THERMAL: Fill thermal block in local matrix
|
||||
# ================================================================
|
||||
|
||||
# Thermal stiffness: K_TT[i,j] = ∫κ(∇Nᵢ·∇Nⱼ) dV
|
||||
@inbounds for i in 1:4, j in 1:4
|
||||
i_local = T_local[i]
|
||||
j_local = T_local[j]
|
||||
K_local[i_local, j_local] += κ * (∇N[i] ⋅ ∇N[j]) * volume
|
||||
end
|
||||
|
||||
# Heat source
|
||||
Q_source = 1.0
|
||||
@inbounds for i in 1:4
|
||||
i_local = T_local[i]
|
||||
F_local[i_local] += Q_source * volume / 4.0
|
||||
end
|
||||
|
||||
# ================================================================
|
||||
# 2️⃣ MECHANICAL: Fill mechanical block in local matrix
|
||||
# ================================================================
|
||||
|
||||
# Lame parameters
|
||||
λ = E * ν / ((1 + ν) * (1 - 2ν))
|
||||
μ = E / (2 * (1 + ν))
|
||||
|
||||
# Build 4th-order elasticity tensor C (isotropic)
|
||||
δ = one(Tensor{2,3})
|
||||
C = λ * δ ⊗ δ + μ * (otimesu(δ, δ) + otimesl(δ, δ))
|
||||
|
||||
# Body force
|
||||
f_body = Vec{3}((0.0, 0.0, -0.1))
|
||||
|
||||
# Fill K_uu and F_u blocks
|
||||
@inbounds for k in 1:4
|
||||
grad_k = ∇N[k]
|
||||
|
||||
# Force vector
|
||||
for α in 1:3
|
||||
i_local = u_local[3*(k-1) + α]
|
||||
F_local[i_local] += (volume / 4.0) * f_body[α]
|
||||
end
|
||||
|
||||
# Stiffness matrix
|
||||
for l in 1:4
|
||||
grad_l = ∇N[l]
|
||||
for α in 1:3, β in 1:3
|
||||
e_α = basevec(Vec{3}, α)
|
||||
e_β = basevec(Vec{3}, β)
|
||||
|
||||
B_k_α = 0.5 * (grad_k ⊗ e_α + e_α ⊗ grad_k)
|
||||
B_l_β = 0.5 * (grad_l ⊗ e_β + e_β ⊗ grad_l)
|
||||
|
||||
k_val = dcontract(B_k_α, dcontract(C, B_l_β)) * volume
|
||||
|
||||
i_local = u_local[3*(k-1) + α]
|
||||
j_local = u_local[3*(l-1) + β]
|
||||
K_local[i_local, j_local] += k_val
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# ================================================================
|
||||
# 3️⃣ HYDRAULIC: Fill pressure block (cell-local)
|
||||
# ================================================================
|
||||
|
||||
K_local[p_local[1], p_local[1]] += k_perm * volume
|
||||
F_local[p_local[1]] += 0.1 * volume
|
||||
|
||||
# ================================================================
|
||||
# 4️⃣ ELECTRIC: Fill electric block (simplified)
|
||||
# ================================================================
|
||||
|
||||
@inbounds for i in 1:6
|
||||
i_local = φ_local[i]
|
||||
K_local[i_local, i_local] += σ_e * volume / 6.0
|
||||
F_local[i_local] += 0.01 * volume / 6.0
|
||||
end
|
||||
|
||||
# ================================================================
|
||||
# 🔗 COUPLING TERMS (This is THE POINT of multi-field elements!)
|
||||
# ================================================================
|
||||
|
||||
# ================================================================
|
||||
# 🔗 COUPLING TERMS - Full Physics Implementation
|
||||
# ================================================================
|
||||
|
||||
# All coupling functions use proper tensor operations - NO simplifications!
|
||||
|
||||
@inline function thermal_expansion_coupling(α_T::Float64, E::Float64, ν::Float64,
|
||||
∇N_T::Vec{3}, ∇N_u::Vec{3},
|
||||
e_α::Vec{3}, vol::Float64)
|
||||
# Full thermo-mechanical coupling: σ = C:ε - α_T·(T-T₀)·I
|
||||
# Linearized: K_Tu = ∫ α_T·E/(1-2ν) · (∇N_T) · (e_α · ∇N_u) dV
|
||||
coupling_strength = α_T * E / (1 - 2*ν)
|
||||
return coupling_strength * (∇N_T ⋅ e_α) * (e_α ⋅ ∇N_u) * vol
|
||||
end
|
||||
|
||||
@inline function biot_coupling(α_p::Float64, ∇N_u::Vec{3}, e_α::Vec{3}, vol::Float64)
|
||||
# Full Biot poroelasticity: σ_eff = σ_total + α_p·p·I
|
||||
# K_up = ∫ α_p · (e_α · ∇N_u) dV (volumetric strain coupling)
|
||||
return α_p * (e_α ⋅ ∇N_u) * vol
|
||||
end
|
||||
|
||||
@inline function thermal_pressurization_coupling(β_T::Float64, ∇N_T::Vec{3}, vol::Float64)
|
||||
# Thermal pressurization in saturated porous media
|
||||
# K_Tp = ∫ β_T · ∇N_T dV (scalar, integrated over volume)
|
||||
# Physically: thermal expansion of pore fluid increases pressure
|
||||
return β_T * norm(∇N_T) * vol
|
||||
end
|
||||
|
||||
@inline function electroosmotic_coupling(ζ::Float64, ∇N_φ::Vec{3}, vol::Float64)
|
||||
# Electro-osmotic flow: fluid flow driven by electric field
|
||||
# K_φp = ∫ ζ · ∇N_φ · ∇N_p dV
|
||||
# Simplified for cell-local pressure (discontinuous)
|
||||
return ζ * norm(∇N_φ) * vol
|
||||
end
|
||||
|
||||
@inline function seebeck_peltier_coupling(S::Float64, ∇N_T::Vec{3}, ∇N_φ::Vec{3}, vol::Float64)
|
||||
# Seebeck effect: J = σ_e·E + S·(-κ∇T)
|
||||
# Peltier effect: Heat flux = Π·J (reciprocal)
|
||||
# K_Tφ = ∫ S · (∇N_T · ∇N_φ) dV
|
||||
return S * (∇N_T ⋅ ∇N_φ) * vol
|
||||
end
|
||||
|
||||
# ----------------------------------------------------------------
|
||||
# Coupling Assembly: Thermo-mechanical (T ↔ u)
|
||||
# Full thermal expansion: σ = C:ε - α_T·E/(1-2ν)·(T-T₀)·I
|
||||
# ----------------------------------------------------------------
|
||||
α_T = 1e-5 # Thermal expansion coefficient [1/K]
|
||||
|
||||
@inbounds for i in 1:4 # Temperature nodes
|
||||
∇N_T = ∇N[i]
|
||||
for k in 1:4 # Displacement nodes
|
||||
∇N_u = ∇N[k]
|
||||
for α in 1:3 # Displacement components (diagonal of I)
|
||||
i_T_local = T_local[i]
|
||||
j_u_local = u_local[3*(k-1) + α]
|
||||
|
||||
e_α = basevec(Vec{3}, α)
|
||||
coupling_val = thermal_expansion_coupling(
|
||||
α_T, E, ν, ∇N_T, ∇N_u, e_α, volume
|
||||
)
|
||||
|
||||
# K_Tu and K_uT blocks (Onsager reciprocity)
|
||||
K_local[i_T_local, j_u_local] += coupling_val
|
||||
K_local[j_u_local, i_T_local] += coupling_val
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# ----------------------------------------------------------------
|
||||
# Coupling Assembly: Hydro-mechanical (u ↔ p)
|
||||
# Full Biot poroelasticity: σ_eff = C:ε - α_p·p·I
|
||||
# ----------------------------------------------------------------
|
||||
α_p = 1e-3 # Biot coefficient (α_p = 1 - K/K_s) [-]
|
||||
|
||||
@inbounds for k in 1:4 # Displacement nodes
|
||||
∇N_u = ∇N[k]
|
||||
for α in 1:3 # Displacement components (trace term)
|
||||
j_u_local = u_local[3*(k-1) + α]
|
||||
p_local_idx = p_local[1]
|
||||
|
||||
e_α = basevec(Vec{3}, α)
|
||||
coupling_val = biot_coupling(α_p, ∇N_u, e_α, volume)
|
||||
|
||||
# K_up and K_pu blocks (Onsager reciprocity)
|
||||
K_local[j_u_local, p_local_idx] += coupling_val
|
||||
K_local[p_local_idx, j_u_local] += coupling_val
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# ----------------------------------------------------------------
|
||||
# Coupling Assembly: Thermal-hydraulic (T ↔ p)
|
||||
# ----------------------------------------------------------------
|
||||
β_T = 1e-6 # Thermal pressurization coefficient
|
||||
|
||||
@inbounds for i in 1:4 # Temperature nodes
|
||||
i_T_local = T_local[i]
|
||||
p_local_idx = p_local[1]
|
||||
|
||||
coupling_val = thermal_pressurization_coupling(β_T, ∇N[i], volume)
|
||||
|
||||
# K_Tp and K_pT blocks (Onsager symmetry)
|
||||
K_local[i_T_local, p_local_idx] += coupling_val
|
||||
K_local[p_local_idx, i_T_local] += coupling_val
|
||||
end
|
||||
|
||||
|
||||
# ----------------------------------------------------------------
|
||||
# Coupling Assembly: Electro-osmotic (p ↔ φ)
|
||||
# Full electrokinetic coupling: v_f = -k/μ_f·∇p + ζ·E
|
||||
# ----------------------------------------------------------------
|
||||
ζ = 1e-7 # Electro-osmotic coefficient [m²/(V·s)]
|
||||
|
||||
@inbounds for i in 1:6 # Electric DOFs on edges
|
||||
i_φ_local = φ_local[i]
|
||||
p_local_idx = p_local[1]
|
||||
|
||||
# Approximate edge gradient (use nodal gradients)
|
||||
node_idx = mod1(i, 4)
|
||||
∇N_φ = ∇N[node_idx]
|
||||
|
||||
coupling_val = electroosmotic_coupling(ζ, ∇N_φ, volume)
|
||||
|
||||
# K_φp and K_pφ blocks (Onsager reciprocity)
|
||||
K_local[i_φ_local, p_local_idx] += coupling_val
|
||||
K_local[p_local_idx, i_φ_local] += coupling_val
|
||||
end
|
||||
|
||||
|
||||
# ----------------------------------------------------------------
|
||||
# Coupling Assembly: Thermo-electric (T ↔ φ)
|
||||
# Full thermoelectric coupling: J = σ_e·E + S·(-κ∇T) (Seebeck)
|
||||
# Q = Π·J (Peltier, where Π = S·T)
|
||||
# ----------------------------------------------------------------
|
||||
S = 1e-6 # Seebeck coefficient [V/K]
|
||||
|
||||
@inbounds for i in 1:4 # Temperature nodes
|
||||
∇N_T = ∇N[i]
|
||||
i_T_local = T_local[i]
|
||||
for j in 1:6 # Electric DOFs on edges
|
||||
j_φ_local = φ_local[j]
|
||||
|
||||
# Approximate edge gradient (use nodal gradients)
|
||||
node_idx = mod1(j, 4)
|
||||
∇N_φ = ∇N[node_idx]
|
||||
|
||||
coupling_val = seebeck_peltier_coupling(S, ∇N_T, ∇N_φ, volume)
|
||||
|
||||
# K_Tφ and K_φT blocks (Onsager reciprocity: Peltier = Seebeck·T)
|
||||
K_local[i_T_local, j_φ_local] += coupling_val
|
||||
K_local[j_φ_local, i_T_local] += coupling_val
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# ----------------------------------------------------------------
|
||||
# Coupling Assembly: Piezoelectric (u ↔ φ) - FULL 3RD ORDER TENSOR!
|
||||
# ----------------------------------------------------------------
|
||||
# Full piezoelectric constitutive laws:
|
||||
# D_k = ε·E_k + e_kij·ε_ij (direct: strain → polarization)
|
||||
# σ_ij = C_ijkl·ε_kl - e_kij·E_k (converse: field → stress)
|
||||
#
|
||||
# Weak form coupling:
|
||||
# K_uφ = ∫ e_kij · (∂N_u^i/∂x_j) · (∂N_φ/∂x_k) dV
|
||||
# K_φu = ∫ e_kij · (∂N_φ/∂x_k) · (∂N_u^i/∂x_j) dV (Onsager reciprocal!)
|
||||
#
|
||||
# For real materials (quartz, PZT, PVDF), e_kij has specific symmetries
|
||||
# Here: simplified diagonal-dominant tensor for demonstration
|
||||
|
||||
# Create 3rd-order piezoelectric tensor e_kij using Tensor{3,3}!
|
||||
# This is THE mathematically correct way - Tensors.jl handles all contractions!
|
||||
e_piezo = Tensor{3,3}((k,i,j) -> k==i==j ? 1e-8 : 0.0)
|
||||
|
||||
# Helper functions for proper tensor contractions
|
||||
@inline function compute_strain_gradient_product(e::Tensor{3,3},
|
||||
∇N_u::Vec{3},
|
||||
∇N_φ::Vec{3},
|
||||
i_comp::Int,
|
||||
vol::Float64)
|
||||
# Contract: e_kij · (∂N_u^i/∂x_j) · (∂N_φ/∂x_k)
|
||||
# This is the FULL piezoelectric coupling integral!
|
||||
result = 0.0
|
||||
for k in 1:3, j in 1:3
|
||||
# e[k,i_comp,j] · (∂N_u/∂x_j) · (∂N_φ/∂x_k)
|
||||
result += e[k,i_comp,j] * ∇N_u[j] * ∇N_φ[k]
|
||||
end
|
||||
return result * vol
|
||||
end
|
||||
|
||||
# Assembly: displacement-electric coupling (FULL tensor contraction!)
|
||||
@inbounds for node_k in 1:4 # Displacement nodes
|
||||
∇N_u = ∇N[node_k]
|
||||
|
||||
for i_comp in 1:3 # Displacement components (stress σ_ij row i)
|
||||
j_u_local = u_local[3*(node_k-1) + i_comp]
|
||||
|
||||
for edge_j in 1:6 # Electric DOFs on edges
|
||||
j_φ_local = φ_local[edge_j]
|
||||
|
||||
# Approximate edge gradient using nodal values
|
||||
node_idx = mod1(edge_j, 4)
|
||||
∇N_φ = ∇N[node_idx]
|
||||
|
||||
# Full tensor contraction: e_kij · (∂u^i/∂x_j) · E_k
|
||||
coupling_val = compute_strain_gradient_product(
|
||||
e_piezo, ∇N_u, ∇N_φ, i_comp, volume
|
||||
)
|
||||
|
||||
# Symmetric (reciprocal) coupling - Onsager reciprocity!
|
||||
# Direct piezoelectric: D = e:ε
|
||||
# Converse piezoelectric: σ = e^T·E (transposed!)
|
||||
K_local[j_u_local, j_φ_local] += coupling_val
|
||||
K_local[j_φ_local, j_u_local] += coupling_val
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# ================================================================
|
||||
# 🚀 SCATTER LOCAL TO GLOBAL (ONE OPERATION!)
|
||||
# ================================================================
|
||||
|
||||
println("\n 📤 Scattering coupled local matrix ($n_local×$n_local) to global")
|
||||
@inbounds for i in 1:n_local
|
||||
I = dof_map[i]
|
||||
F[I] += F_local[i]
|
||||
for j in 1:n_local
|
||||
J = dof_map[j]
|
||||
K[I, J] += K_local[i, j]
|
||||
end
|
||||
end
|
||||
end # End of element loop
|
||||
|
||||
println("\n✓ Assembly complete!")
|
||||
println(" ONE coupled system matrix: $(size(K))")
|
||||
println(" Total non-zeros: $(nnz(K))")
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("APPLYING BOUNDARY CONDITIONS AND SOLVING")
|
||||
println("="^70)
|
||||
|
||||
# Apply BCs - properly constrain ALL fields to avoid singularity!
|
||||
#
|
||||
# Physical interpretation:
|
||||
# - Node 1: Fully grounded (T=0, u=0, reference for all fields)
|
||||
# - Node 2: Prevent rigid motion in x (ux=0)
|
||||
# - Electric: Ground edge 1 to prevent floating potential (φ_edge1=0)
|
||||
|
||||
bc_dofs = [
|
||||
1, # T at node 1 (thermal ground)
|
||||
n_T_total+1, n_T_total+2, n_T_total+3, # u at node 1 (mechanical ground)
|
||||
n_T_total+4, # ux at node 2 (prevent x-rotation)
|
||||
n_T_total+n_u_total+n_p_total+1 # φ at edge 1 (electric ground)
|
||||
]
|
||||
|
||||
for dof in bc_dofs
|
||||
K[dof, :] .= 0.0
|
||||
K[:, dof] .= 0.0
|
||||
K[dof, dof] = 1.0
|
||||
F[dof] = 0.0
|
||||
end
|
||||
|
||||
println("\nBoundary conditions (FULL MULTI-PHYSICS):")
|
||||
println(" Thermal: Node 1 fixed at T=0 K (thermal ground)")
|
||||
println(" Mechanical: Node 1 fully fixed u=(0,0,0) (mechanical ground)")
|
||||
println(" Mechanical: Node 2 ux=0 (prevent rigid rotation)")
|
||||
println(" Electric: Edge 1 fixed at φ=0 V (electric ground)")
|
||||
println(" Hydraulic: Natural BCs (traction-free, no flow prescribed)")
|
||||
|
||||
# Solve
|
||||
println("\n🎯 Solving coupled system...")
|
||||
println(" Matrix size: $(size(K))")
|
||||
println(" Non-zeros: $(nnz(K))")
|
||||
println(" Condition number estimate: checking...")
|
||||
|
||||
# Add small regularization to prevent singularity from weakly coupled terms
|
||||
# This is physically reasonable - represents small stabilization
|
||||
ε_reg = 1e-12
|
||||
for i in 1:n_total
|
||||
K[i,i] += ε_reg
|
||||
end
|
||||
|
||||
println(" Added regularization (ε=$ε_reg) for numerical stability")
|
||||
|
||||
# Solve using robust method
|
||||
sol = try
|
||||
result = K \ F
|
||||
println(" ✓ Solution converged!")
|
||||
result
|
||||
catch e
|
||||
println(" ERROR: System still singular!")
|
||||
println(" This indicates physical model needs more constraints")
|
||||
rethrow(e)
|
||||
end
|
||||
|
||||
# Extract fields
|
||||
T_sol = sol[1:n_T_total]
|
||||
u_sol = sol[n_T_total+1:n_T_total+n_u_total]
|
||||
p_sol = sol[n_T_total+n_u_total+1:n_T_total+n_u_total+n_p_total]
|
||||
φ_sol = sol[n_T_total+n_u_total+n_p_total+1:end]
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("✨ SOLUTION (REAL PHYSICS!)")
|
||||
println("="^70)
|
||||
|
||||
println("\n🌡️ Temperature field:")
|
||||
for i in 1:length(T_sol)
|
||||
println(" Node $i: T = $(@sprintf("%.6f", T_sol[i])) K")
|
||||
end
|
||||
|
||||
println("\n🏗️ Displacement field:")
|
||||
n_disp_nodes = div(length(u_sol), 3)
|
||||
for node_id in 1:n_disp_nodes
|
||||
ux = u_sol[3*(node_id-1)+1]
|
||||
uy = u_sol[3*(node_id-1)+2]
|
||||
uz = u_sol[3*(node_id-1)+3]
|
||||
println(" Node $node_id: u = ($(@sprintf("%.6f", ux)), $(@sprintf("%.6f", uy)), $(@sprintf("%.6f", uz))) m")
|
||||
end
|
||||
|
||||
println("\n💧 Pore pressure field:")
|
||||
for i in 1:length(p_sol)
|
||||
println(" Cell $i: p = $(@sprintf("%.6f", p_sol[i])) Pa")
|
||||
end
|
||||
|
||||
println("\n⚡ Electric potential (edges):")
|
||||
for i in 1:length(φ_sol)
|
||||
println(" Edge $i: φ = $(@sprintf("%.6f", φ_sol[i])) V")
|
||||
end
|
||||
|
||||
# Verification tests
|
||||
@test all(isfinite.(T_sol))
|
||||
@test all(isfinite.(u_sol))
|
||||
@test all(isfinite.(p_sol))
|
||||
@test all(isfinite.(φ_sol))
|
||||
|
||||
@test T_sol[1] ≈ 0.0 atol=1e-10 # BC
|
||||
@test u_sol[1:3] ≈ [0.0, 0.0, 0.0] atol=1e-10 # BC
|
||||
|
||||
# Check non-trivial solution
|
||||
@test maximum(abs.(T_sol[2:end])) > 1e-6
|
||||
@test maximum(abs.(u_sol[4:end])) > 1e-6
|
||||
@test maximum(abs.(p_sol)) > 1e-6
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("🎉 ACHIEVEMENTS UNLOCKED:")
|
||||
println("="^70)
|
||||
println(" ✅ ONE element type with FOUR physics fields!")
|
||||
println(" ✅ ONE local coupled matrix per element (23×23)")
|
||||
println(" ✅ ALL physics assembled together (true coupling!)")
|
||||
println(" ✅ Thermo-mechanical coupling: K_Tu, K_uT (thermal expansion)")
|
||||
println(" ✅ Hydro-mechanical coupling: K_up, K_pu (Biot poroelasticity)")
|
||||
println(" ✅ Thermal-hydraulic coupling: K_Tp, K_pT (thermal pressurization)")
|
||||
println(" ✅ Electro-osmotic coupling: K_φp, K_pφ (electrokinetic flow)")
|
||||
println(" ✅ Thermo-electric coupling: K_Tφ, K_φT (Seebeck/Peltier)")
|
||||
println(" ✅ Piezoelectric coupling: K_uφ, K_φu (Tensor{3,3} elegance!)")
|
||||
println(" ✅ Total: 12 off-diagonal coupling blocks! (ALL physics coupled!)")
|
||||
println(" ✅ Modular coupling functions (inlined for zero overhead)")
|
||||
println(" ✅ 3rd-order tensor formulation (e_kij via Tensor{3,3})")
|
||||
println(" ✅ Onsager reciprocity respected (all couplings symmetric)")
|
||||
println(" ✅ Local-to-global mapping via type system")
|
||||
println(" ✅ REAL thermal diffusion (∫κ∇T·∇T' dV)")
|
||||
println(" ✅ REAL 3D elasticity (∫C:ε:ε dV)")
|
||||
println(" ✅ Zero-allocation Tensors.jl operations")
|
||||
println(" ✅ get_basis_derivatives API (no manual gradients)")
|
||||
println(" ✅ Cell-local pressure DOFs (discontinuous)")
|
||||
println(" ✅ Edge-based electric DOFs")
|
||||
println(" ✅ Full $n_total × $n_total coupled system solved")
|
||||
println(" ✅ Type-safe field access: .T, .u, .p, .φ")
|
||||
println(" ✅ field_dof_range() - extract field blocks (compile-time!)")
|
||||
println(" ✅ local_to_global_map() - scatter operation")
|
||||
println("="^70)
|
||||
|
||||
println("\n💡 THIS IS THE POWER OF MULTI-FIELD ELEMENTS!")
|
||||
println(" ONE element → ONE local matrix → ALL physics coupled!")
|
||||
println(" T ↔ u (thermal expansion), T ↔ p (thermal pressurization)")
|
||||
println(" T ↔ φ (Seebeck/Peltier), u ↔ p (Biot poroelasticity)")
|
||||
println(" u ↔ φ (piezoelectric via Tensor{3,3}!), p ↔ φ (electro-osmotic)")
|
||||
println(" → Complete multi-physics: 4 fields × 6 couplings = 12 blocks!")
|
||||
println(" → Tensors.jl elegance: 3rd-order piezoelectric tensor!")
|
||||
println(" → Modular design: coupling functions inlined for performance!")
|
||||
println(" → Geothermal, nuclear waste, CO2 sequestration, smart materials!")
|
||||
println(" Natural coupling, type-safe, composable, ELEGANT! 🚀")
|
||||
|
||||
end
|
||||
Reference in New Issue
Block a user