mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
test(solvers): add modal analysis new API visionary test
New 560-line test file for modal analysis NEW API (test-driven development): - Visionary tests for ModalSolver and EigenvalueProblem - Tests single Tet4 element modal analysis - Tests with and without geometric stiffness - Tests heat transfer modal analysis - Extensive documentation of intended API design - Tests currently skipped (@test_skip) until implementation complete - Documents eigenvalue problem K φ = λ M φ - Documents mass matrix assembly and mode shape normalization This test file serves as both test suite and API design documentation for the new modal analysis solver interface.
This commit is contained in:
Jukka Aho
@@ -0,0 +1,560 @@
|
||||
"""
|
||||
# Modal Analysis - NEW API (Test-Driven Development)
|
||||
|
||||
**What:** Shows how modal (eigenvalue) analysis SHOULD work with the NEW API
|
||||
|
||||
**Why:**
|
||||
- **Eigenvalue problems** - Different solver type (not Newton-Krylov)
|
||||
- **Natural frequencies** - Validates dynamic behavior
|
||||
- **Geometric stiffness** - Prestress effects on frequencies
|
||||
- **Multiple physics** - Works for elasticity AND heat transfer
|
||||
- **Mass matrix required** - First time we need M, not just K
|
||||
|
||||
**NEW API Concepts:**
|
||||
1. **ModalSolver** - Eigenvalue solver (not LinearSolver or NewtonSolver)
|
||||
2. **EigenvalueProblem** - K φ = λ M φ (generalized eigenvalue problem)
|
||||
3. **Geometric stiffness** - K_g from initial stress/displacement
|
||||
4. **Modal shapes** - φ_i (eigenvectors = mode shapes)
|
||||
5. **Natural frequencies** - ω_i = √λ_i
|
||||
|
||||
**Test Problems:**
|
||||
|
||||
## Test 1: Single Tet4 Element
|
||||
- 4-node tetrahedron with 3 nodes fixed
|
||||
- Validates basic eigenvalue computation
|
||||
- Tests with and without geometric stiffness
|
||||
|
||||
## Test 2: Heat Transfer Modal Analysis
|
||||
- Two Quad4 elements (plane heat)
|
||||
- Eigenvalues of thermal diffusion operator
|
||||
- Validates modal analysis works for temperature field
|
||||
|
||||
## Test 3: Vibrating Beam (Future)
|
||||
- Fixed-fixed beam
|
||||
- Compare natural frequencies to analytical solution
|
||||
- ω_i = λ_i² √(EI/ρA)
|
||||
|
||||
**Expected Behavior (when implemented):**
|
||||
✅ Eigenvalues computed correctly
|
||||
✅ Mode shapes orthogonal (φ_i^T M φ_j = δ_ij)
|
||||
✅ Geometric stiffness affects frequencies
|
||||
✅ Works for both elasticity and heat transfer
|
||||
✅ Sparse eigenvalue solver (only need first k modes)
|
||||
|
||||
**Status:** 🚧 VISIONARY TEST - Implementation in progress
|
||||
"""
|
||||
|
||||
using Test
|
||||
using JuliaFEM
|
||||
using Tensors
|
||||
using LinearAlgebra
|
||||
using Statistics
|
||||
|
||||
@testset "Modal Analysis - NEW API (TDD)" begin
|
||||
|
||||
# =============================================================================
|
||||
# SINGLE ELEMENT: TET4 WITHOUT GEOMETRIC STIFFNESS
|
||||
# =============================================================================
|
||||
|
||||
@testset "Single Tet4 Modal Analysis (Visionary)" begin
|
||||
@test_skip begin # Skip until implemented
|
||||
|
||||
# Geometry
|
||||
X = Dict(
|
||||
1 => Vec(2.0, 3.0, 4.0),
|
||||
2 => Vec(6.0, 3.0, 2.0),
|
||||
3 => Vec(2.0, 5.0, 1.0),
|
||||
4 => Vec(4.0, 3.0, 6.0)
|
||||
)
|
||||
|
||||
# Initial displacement (for geometric stiffness later)
|
||||
u0 = Dict(
|
||||
1 => Vec(0.0, 0.0, 0.0),
|
||||
2 => Vec(0.0, 0.0, 0.0),
|
||||
3 => Vec(0.0, 0.0, 0.0),
|
||||
4 => Vec(0.25, 0.25, 0.25)
|
||||
)
|
||||
|
||||
# Material
|
||||
material = LinearElastic(E=96.0, ν=1 / 3, ρ=420.0)
|
||||
|
||||
# Physics
|
||||
elastic_physics = ContinuumPhysics{Displacement}(
|
||||
material=material,
|
||||
formulation=FullThreeD(),
|
||||
finite_strain=false,
|
||||
geometric_stiffness=false # No K_g initially
|
||||
)
|
||||
|
||||
# Create element
|
||||
tet_element = create_element(
|
||||
Tet4, (1, 2, 3, 4),
|
||||
geometry=X,
|
||||
displacement=u0
|
||||
)
|
||||
|
||||
# Domain
|
||||
domain = Domain(
|
||||
name="TET",
|
||||
elements=[tet_element],
|
||||
physics=elastic_physics
|
||||
)
|
||||
|
||||
# Boundary conditions (3 nodes fixed)
|
||||
bc_fixed = DirichletBC(
|
||||
name="FIXED_FACE",
|
||||
nodes=[1, 2, 3],
|
||||
dof=:displacement,
|
||||
values=[0.0, 0.0, 0.0]
|
||||
)
|
||||
|
||||
# NEW: EigenvalueProblem (not static or transient)
|
||||
problem = EigenvalueProblem(
|
||||
domains=[domain],
|
||||
boundary_conditions=[bc_fixed]
|
||||
)
|
||||
|
||||
# NEW: ModalSolver
|
||||
solver = ModalSolver(
|
||||
n_modes=2, # Number of modes to compute
|
||||
which=:SM, # Smallest Magnitude (or :LM for largest)
|
||||
method=:arpack # ARPACK, KrylovKit, or direct
|
||||
)
|
||||
|
||||
# Solve: K φ = λ M φ
|
||||
solution = solve!(problem, solver)
|
||||
|
||||
# Extract eigenvalues and eigenvectors
|
||||
λ = solution.eigenvalues # [λ_1, λ_2]
|
||||
φ = solution.eigenvectors # [φ_1, φ_2]
|
||||
|
||||
# Validate eigenvalues (without geometric stiffness)
|
||||
@test isapprox(λ, [4 / 3, 1 / 3], rtol=1e-6)
|
||||
|
||||
# Validate orthogonality: φ_i^T M φ_j = δ_ij
|
||||
for i in 1:2
|
||||
for j in 1:2
|
||||
orthogonality = φ[i]' * solution.mass_matrix * φ[j]
|
||||
expected = (i == j) ? 1.0 : 0.0
|
||||
@test isapprox(orthogonality, expected, atol=1e-6)
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# GEOMETRIC STIFFNESS EFFECT
|
||||
# =============================================================================
|
||||
|
||||
@testset "Tet4 with Geometric Stiffness (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
# Same setup as before, but enable geometric stiffness
|
||||
elastic_physics = ContinuumPhysics{Displacement}(
|
||||
material=LinearElastic(E=96.0, ν=1 / 3, ρ=420.0),
|
||||
formulation=FullThreeD(),
|
||||
finite_strain=false,
|
||||
geometric_stiffness=true # NEW: Enable K_g
|
||||
)
|
||||
|
||||
problem = EigenvalueProblem(
|
||||
domains=[domain],
|
||||
boundary_conditions=[bc_fixed],
|
||||
initial_displacement=u0 # Required for K_g
|
||||
)
|
||||
|
||||
solver = ModalSolver(n_modes=2, which=:SM)
|
||||
solution = solve!(problem, solver)
|
||||
|
||||
λ = solution.eigenvalues
|
||||
|
||||
# With geometric stiffness: K_eff = K + K_g(u0)
|
||||
# Eigenvalues should be different!
|
||||
@test isapprox(λ, [5 / 3, 2 / 3], rtol=1e-6)
|
||||
|
||||
# Geometric stiffness changes natural frequencies
|
||||
@test λ[1] > 4 / 3 # Stiffened by prestress
|
||||
@test λ[2] > 1 / 3
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# HEAT TRANSFER MODAL ANALYSIS
|
||||
# =============================================================================
|
||||
|
||||
@testset "Heat Transfer Modal (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
# Geometry: Two quad elements stacked vertically
|
||||
X = Dict(
|
||||
1 => Vec(0.0, 0.0),
|
||||
2 => Vec(1.0, 0.0),
|
||||
3 => Vec(1.0, 3.0),
|
||||
4 => Vec(0.0, 3.0),
|
||||
5 => Vec(0.0, 3.0),
|
||||
6 => Vec(1.0, 3.0),
|
||||
7 => Vec(1.0, 9.0),
|
||||
8 => Vec(0.0, 9.0)
|
||||
)
|
||||
|
||||
# Material (thermal)
|
||||
thermal_material = ThermalMaterial(
|
||||
conductivity=36.0,
|
||||
density=6.0,
|
||||
specific_heat=1.0
|
||||
)
|
||||
|
||||
# Physics
|
||||
thermal_physics = ContinuumPhysics{Temperature}(
|
||||
material=thermal_material,
|
||||
formulation=PlaneHeat()
|
||||
)
|
||||
|
||||
# Elements
|
||||
elem1 = create_element(Quad4, (1, 2, 3, 4), geometry=X)
|
||||
elem2 = create_element(Quad4, (4, 3, 7, 8), geometry=X)
|
||||
|
||||
# Domain
|
||||
domain = Domain(
|
||||
name="THERMAL",
|
||||
elements=[elem1, elem2],
|
||||
physics=thermal_physics
|
||||
)
|
||||
|
||||
# Boundary conditions (fixed temperature at ends)
|
||||
bc_bottom = DirichletBC(
|
||||
nodes=[1, 2],
|
||||
dof=:temperature,
|
||||
value=0.0
|
||||
)
|
||||
|
||||
bc_top = DirichletBC(
|
||||
nodes=[7, 8],
|
||||
dof=:temperature,
|
||||
value=0.0
|
||||
)
|
||||
|
||||
# Eigenvalue problem for heat equation
|
||||
# M ∂T/∂t + K T = 0
|
||||
# Modal: K φ = λ M φ
|
||||
problem = EigenvalueProblem(
|
||||
domains=[domain],
|
||||
boundary_conditions=[bc_bottom, bc_top]
|
||||
)
|
||||
|
||||
solver = ModalSolver(n_modes=1, which=:SM)
|
||||
solution = solve!(problem, solver)
|
||||
|
||||
λ = solution.eigenvalues
|
||||
|
||||
# First eigenvalue should be 1.0 (analytical)
|
||||
@test isapprox(λ[1], 1.0, rtol=1e-6)
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# NATURAL FREQUENCIES FROM EIGENVALUES
|
||||
# =============================================================================
|
||||
|
||||
@testset "Natural Frequency Computation (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
# Solve modal problem
|
||||
solution = solve_modal_problem()
|
||||
|
||||
# Extract eigenvalues
|
||||
λ = solution.eigenvalues # [rad²/s²]
|
||||
|
||||
# Natural frequencies: ω = √λ [rad/s]
|
||||
ω = sqrt.(λ)
|
||||
|
||||
# Convert to Hz: f = ω / (2π)
|
||||
f = ω ./ (2π)
|
||||
|
||||
# Validate units and values
|
||||
@test all(λ .≥ 0) # Eigenvalues non-negative
|
||||
@test all(ω .≥ 0) # Frequencies non-negative
|
||||
@test all(f .≥ 0)
|
||||
|
||||
# For simple beam: f_1 ≈ 1-10 Hz (typical)
|
||||
@test 1.0 < f[1] < 100.0
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# MODE SHAPE VISUALIZATION (FUTURE)
|
||||
# =============================================================================
|
||||
|
||||
@testset "Mode Shape Properties (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
solution = solve_modal_problem()
|
||||
|
||||
φ = solution.eigenvectors # Mode shapes
|
||||
M = solution.mass_matrix
|
||||
K = solution.stiffness_matrix
|
||||
λ = solution.eigenvalues
|
||||
|
||||
# Property 1: Orthogonality w.r.t. mass matrix
|
||||
# φ_i^T M φ_j = δ_ij
|
||||
for i in 1:length(φ)
|
||||
for j in 1:length(φ)
|
||||
orth_M = φ[i]' * M * φ[j]
|
||||
expected = (i == j) ? 1.0 : 0.0
|
||||
@test isapprox(orth_M, expected, atol=1e-6)
|
||||
end
|
||||
end
|
||||
|
||||
# Property 2: Eigenvalue equation
|
||||
# K φ_i = λ_i M φ_i
|
||||
for i in 1:length(φ)
|
||||
Kφ = K * φ[i]
|
||||
λMφ = λ[i] * (M * φ[i])
|
||||
@test Kφ ≈ λMφ rtol = 1e-6
|
||||
end
|
||||
|
||||
# Property 3: Mode shapes normalized
|
||||
# φ_i^T M φ_i = 1
|
||||
for i in 1:length(φ)
|
||||
norm_M = φ[i]' * M * φ[i]
|
||||
@test isapprox(norm_M, 1.0, rtol=1e-6)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# MASS MATRIX ASSEMBLY
|
||||
# =============================================================================
|
||||
|
||||
@testset "Mass Matrix Assembly (Visionary)" begin
|
||||
# Pseudo-code showing mass matrix assembly
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("MASS MATRIX ASSEMBLY (NODAL APPROACH)")
|
||||
println("="^70)
|
||||
|
||||
mass_assembly_pseudo = """
|
||||
# Similar to stiffness, but uses density ρ and N^T N
|
||||
for node_i in nodes
|
||||
for elem in node_to_elements[node_i]
|
||||
for node_j in elem.nodes
|
||||
# Mass matrix block (3×3 for displacement)
|
||||
M_ij = ∫_Ω ρ N_i N_j dΩ
|
||||
|
||||
# For displacement (3D):
|
||||
# M_ij = m_ij * I_3×3 (often lumped)
|
||||
|
||||
# Consistent mass (full integration)
|
||||
M_block = compute_mass_block(elem, node_i, node_j, ρ)
|
||||
|
||||
# Or lumped mass (diagonal only)
|
||||
if lumped
|
||||
M_block = (i == j) ? diag(M_block) : 0
|
||||
end
|
||||
|
||||
M_nodal[node_i, node_j] += M_block
|
||||
end
|
||||
end
|
||||
end
|
||||
"""
|
||||
|
||||
println(mass_assembly_pseudo)
|
||||
println("="^70)
|
||||
println("✓ Mass matrix M = ∫_Ω ρ N^T N dΩ")
|
||||
println("✓ Consistent mass: Full integration (more accurate)")
|
||||
println("✓ Lumped mass: Diagonal only (faster, explicit dynamics)")
|
||||
println("✓ Nodal assembly works same as stiffness!")
|
||||
println("="^70)
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# KEY ARCHITECTURAL INSIGHTS
|
||||
# =============================================================================
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("MODAL ANALYSIS ARCHITECTURE INSIGHTS (NEW API)")
|
||||
println("="^70)
|
||||
println("✓ Modal analysis is EigenvalueProblem (K φ = λ M φ)")
|
||||
println("✓ Requires both stiffness K AND mass M matrices")
|
||||
println("✓ ModalSolver uses sparse eigenvalue methods (ARPACK)")
|
||||
println("✓ Geometric stiffness K_g affects natural frequencies")
|
||||
println("✓ Works for ANY physics (displacement, temperature, etc.)")
|
||||
println("✓ Mode shapes φ orthogonal w.r.t. mass matrix")
|
||||
println("✓ Natural frequencies ω = √λ, in Hz: f = ω/(2π)")
|
||||
println("✓ Same nodal assembly pattern for M as for K!")
|
||||
println("="^70)
|
||||
|
||||
end
|
||||
|
||||
"""
|
||||
# IMPLEMENTATION NOTES
|
||||
|
||||
## Generalized Eigenvalue Problem
|
||||
|
||||
**Mathematical formulation:**
|
||||
|
||||
K φ = λ M φ
|
||||
|
||||
where:
|
||||
- K = stiffness matrix (from elastic/thermal energy)
|
||||
- M = mass matrix (from kinetic/thermal energy)
|
||||
- φ = eigenvector (mode shape)
|
||||
- λ = eigenvalue (ω² for vibrations, α for heat)
|
||||
|
||||
**For structural vibrations:**
|
||||
- ω = √λ (natural frequency in rad/s)
|
||||
- f = ω/(2π) (natural frequency in Hz)
|
||||
- φ = displacement mode shape
|
||||
|
||||
**For heat diffusion:**
|
||||
- α = λ (thermal diffusivity eigenvalue)
|
||||
- φ = temperature mode shape
|
||||
|
||||
## Stiffness Matrix
|
||||
|
||||
**Elasticity:**
|
||||
|
||||
K = ∫_Ω B^T C B dΩ
|
||||
|
||||
where:
|
||||
- B = strain-displacement matrix
|
||||
- C = elasticity tensor
|
||||
|
||||
**Heat transfer:**
|
||||
|
||||
K = ∫_Ω k ∇N^T ∇N dΩ
|
||||
|
||||
where:
|
||||
- k = thermal conductivity
|
||||
|
||||
## Mass Matrix
|
||||
|
||||
**Elasticity (consistent):**
|
||||
|
||||
M = ∫_Ω ρ N^T N dΩ
|
||||
|
||||
where ρ = density.
|
||||
|
||||
**Elasticity (lumped):**
|
||||
|
||||
M_ii = ∑_elements ∫_Ω_e ρ N_i dΩ_e
|
||||
|
||||
(Diagonal only, faster for explicit dynamics)
|
||||
|
||||
**Heat transfer:**
|
||||
|
||||
M = ∫_Ω ρc N^T N dΩ
|
||||
|
||||
where:
|
||||
- ρ = density
|
||||
- c = specific heat
|
||||
|
||||
## Geometric Stiffness
|
||||
|
||||
**Definition:** Stiffness contribution from initial stress state.
|
||||
|
||||
K_g = ∫_Ω G^T σ₀ G dΩ
|
||||
|
||||
where:
|
||||
- G = geometric matrix (relates δε to ∇(δu))
|
||||
- σ₀ = initial stress
|
||||
|
||||
**Effect:** Prestress changes natural frequencies:
|
||||
- Tension → higher frequencies (stiffening)
|
||||
- Compression → lower frequencies (softening)
|
||||
|
||||
**Total effective stiffness:**
|
||||
|
||||
K_eff = K + K_g
|
||||
|
||||
Then solve: K_eff φ = λ M φ
|
||||
|
||||
## Sparse Eigenvalue Solvers
|
||||
|
||||
**Problem:** For large systems, computing ALL eigenvalues is expensive.
|
||||
|
||||
**Solution:** Sparse eigenvalue methods (compute only k << n eigenvalues).
|
||||
|
||||
**Methods:**
|
||||
1. **ARPACK** - Arnoldi iteration (standard in Julia)
|
||||
2. **KrylovKit.jl** - Modern Krylov methods
|
||||
3. **Lanczos** - For symmetric problems (K and M symmetric)
|
||||
|
||||
**JuliaFEM approach:**
|
||||
|
||||
```julia
|
||||
using Arpack
|
||||
|
||||
# Solve K φ = λ M φ for k smallest eigenvalues
|
||||
λ, φ = eigs(K, M, nev=k, which=:SM)
|
||||
```
|
||||
|
||||
## Mode Shape Normalization
|
||||
|
||||
**Goal:** Normalize eigenvectors for convenience.
|
||||
|
||||
**Mass normalization (standard):**
|
||||
|
||||
φ_i^T M φ_i = 1
|
||||
|
||||
**Advantages:**
|
||||
- Orthogonality: φ_i^T M φ_j = δ_ij
|
||||
- Energy interpretation clear
|
||||
- Modal damping ratios well-defined
|
||||
|
||||
**Implementation:**
|
||||
|
||||
```julia
|
||||
function normalize_modes!(φ, M)
|
||||
for i in 1:length(φ)
|
||||
# Compute φ_i^T M φ_i
|
||||
norm_M = φ[i]' * M * φ[i]
|
||||
|
||||
# Normalize
|
||||
φ[i] ./= sqrt(norm_M)
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
## Nodal Assembly for Mass Matrix
|
||||
|
||||
```julia
|
||||
function assemble_mass_nodal!(M_nodal, ρ, elements, node_to_elements)
|
||||
Threads.@threads for node_i in 1:n_nodes
|
||||
for elem_idx in node_to_elements[node_i]
|
||||
elem = elements[elem_idx]
|
||||
|
||||
for node_j in elem.nodes
|
||||
# Mass matrix block (3×3 for displacement)
|
||||
M_ij = compute_mass_block(elem, node_i, node_j, ρ)
|
||||
|
||||
# Add to global (no atomic on diagonal if i = j)
|
||||
if node_i == node_j
|
||||
M_nodal[node_i, node_i] += M_ij # Direct write
|
||||
else
|
||||
atomic_add!(M_nodal[node_i, node_j], M_ij)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
**Key:** SAME PATTERN as stiffness assembly!
|
||||
|
||||
## Next Steps
|
||||
|
||||
1. Implement `EigenvalueProblem` type
|
||||
2. Implement `ModalSolver` with ARPACK integration
|
||||
3. Implement mass matrix assembly (nodal)
|
||||
4. Implement geometric stiffness computation
|
||||
5. Add mode shape normalization
|
||||
6. Validate against analytical solutions
|
||||
7. Performance benchmarks (large systems)
|
||||
|
||||
"""
|
||||
Reference in New Issue
Block a user