diff --git a/ROADMAP.md b/ROADMAP.md index 50d663a..9ac48ea 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -26,7 +26,7 @@ hex meshes today and it gives the best elements where stresses are read. | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | | M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, connectors; amplitudes, multiple steps; nonlinear static steps in decks | planned | -| M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; heat transfer | planned | +| M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started | | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned | | M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `modelcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started | | M11 Equations of state and shells | Polynomial, ideal-gas, Mie–Grüneisen, JWL with programmed burn, tabulated; energy coupling; shock verification (Rankine–Hugoniot, Sod, Noh, Sedov); Belytschko–Lin–Tsay shells with self-contact on quad meshes | planned | diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 03710c9..77ba472 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -294,6 +294,7 @@ fn check_expressible(model: &Model) -> Result<()> { )); } } + bad.extend(crate::thermal_items(model)); for (name, c) in model.couplings() { if c.kind == CouplingKind::Distributing && !(1..=3).all(|d| c.dofs.contains(d)) { bad.insert(format!("coupling `{name}`: an Abaqus distributing coupling always couples translations 1-3")); @@ -911,6 +912,8 @@ impl Deck { }; self.line(&format!("{target}, {s}")); } + // Refused up front (check_expressible). + BoundaryCondition::Temperature { .. } => {} BoundaryCondition::Displacement { values, .. } => { for (dof, v) in (1u8..).zip(values) { if let Some(v) = v { @@ -1068,8 +1071,8 @@ impl Deck { self.comment("loads have no effect on natural frequencies and are not written"); } } - Procedure::ExplicitDynamic(_) => { - return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: explicit dynamic steps", step.name) }); + Procedure::ExplicitDynamic(_) | Procedure::HeatTransfer { .. } | Procedure::CoupledTemperatureDisplacement {} => { + return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: {} steps", step.name, step.procedure.name()) }); } } self.output(step); diff --git a/crates/decks/src/lib.rs b/crates/decks/src/lib.rs index 973aaa9..35ba863 100644 --- a/crates/decks/src/lib.rs +++ b/crates/decks/src/lib.rs @@ -115,6 +115,23 @@ pub(crate) fn unsupported(format: &'static str, items: BTreeSet) -> Resu if items.is_empty() { Ok(()) } else { Err(Error::Unsupported { format, what: items.into_iter().collect::>().join("; ") }) } } +/// Heat transfer in a model, which no writer expresses yet: heat transfer and coupled +/// temperature–displacement steps, temperature boundary conditions, heat loads and temperature +/// fields. Thermal material properties alone are left out of the decks (nothing uses them). +pub(crate) fn thermal_items(model: &Model) -> BTreeSet { + let mut out = BTreeSet::new(); + for step in model.steps().iter().filter(|s| s.procedure.is_thermal()) { + out.insert(format!("step `{}`: {} steps (not written to decks yet)", step.name, step.procedure.name())); + } + for name in model.boundary_conditions().iter().filter(|(_, b)| b.is_thermal()).map(|(n, _)| n) { + out.insert(format!("boundary condition `{name}`: prescribed temperatures (not written to decks yet)")); + } + for (name, load) in model.loads().iter().filter(|(_, l)| l.is_thermal() || matches!(l, Load::Temperature { .. })) { + out.insert(format!("load `{name}`: {} loads (not written to decks yet)", load.kind())); + } + out +} + /// The concentrated loads of a set of active loads: the resolved nodal values of every force, /// traction (degrees of freedom 1–3) and moment (4–6), summed per node and degree of freedom so /// loads sharing a node add up; zero sums are left out. diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index c011f97..28610ae 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -202,7 +202,7 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { let used: BTreeSet<&str> = model.sections().values().map(|s| s.material.as_str()).collect(); let materials: BTreeMap<&str, u32> = used.into_iter().enumerate().map(|(i, m)| Ok((m, id(i)?))).collect::>()?; // Everything LS-DYNA cannot express, reported at once. - let mut refused = BTreeSet::new(); + let mut refused = crate::thermal_items(model); let mut parts = Vec::new(); let mut hourglasses: BTreeMap = BTreeMap::new(); for (name, sec) in model.sections() { @@ -311,6 +311,8 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { Load::BodyForce { .. } => { refused.insert(format!("load `{l}`: a body force per unit volume (LS-DYNA's body loads are accelerations; use gravity)")); } + // In thermal_items. + Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. } | Load::Temperature { .. } => {} } } } @@ -818,8 +820,8 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste d.card([Field::Blank, Field::Blank, Field::Blank, Field::Blank, Field::from(1)]); } } - Procedure::ExplicitDynamic(_) => { - return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: explicit dynamic steps", step.name) }); + Procedure::ExplicitDynamic(_) | Procedure::HeatTransfer { .. } | Procedure::CoupledTemperatureDisplacement {} => { + return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: {} steps", step.name, step.procedure.name()) }); } } d.line("*CONTROL_TERMINATION"); @@ -872,6 +874,9 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste } Load::BodyForce { .. } => return Err(broken(format!("load `{name}`: a body force"))), Load::Centrifugal { .. } => return Err(broken(format!("load `{name}`: a centrifugal load"))), + Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. } | Load::Temperature { .. } => { + return Err(broken(format!("load `{name}`: a {} load", load.kind()))); + } } } // The other concentrated loads node by node, summed where they share a node. diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index 3455abd..96c0a51 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -327,6 +327,7 @@ impl<'m> Deck<'m> { self.plan_discrete()?; self.plan_constraints()?; self.plan_loads()?; + self.problems.extend(crate::thermal_items(self.model)); for sub in &mut self.subcases { if matches!(sub.step.procedure, Procedure::Frequency { .. }) { sub.method = Some(next_id(&mut self.next_sid)?); @@ -705,6 +706,8 @@ impl<'m> Deck<'m> { OutputVariable::ReactionForce => "SPCFORCES", OutputVariable::Stress => "STRESS", OutputVariable::Strain => "STRAIN", + // Heat transfer steps are refused (thermal_items). + OutputVariable::Temperature | OutputVariable::HeatFlux => continue, }; let _ = writeln!(out, " {request} = ALL"); } @@ -928,7 +931,7 @@ impl<'m> Deck<'m> { let a = angular_velocity / (2.0 * std::f64::consts::PI); entry(Card::new("RFORCE").int(*sid).blank().blank().real(a).reals(*axis), out, &mut first)? } - Load::BodyForce { .. } => {} + Load::BodyForce { .. } | Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. } | Load::Temperature { .. } => {} } } for (sid, loads) in &self.combinations { diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index d405ebf..94f7c1e 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -418,7 +418,14 @@ pub static CATALOGUE: &[Capability] = &[ deck("deck.write.nastran", "Nastran Deck", "nastran", "Write Nastran bulk data (SOL 101, 103)"), deck("deck.write.lsdyna", "LS-DYNA Deck", "lsdyna", "Write an LS-DYNA keyword file"), soon("deck.read", "Read Deck", Solve, "Decks", "M9", "Read a deck, keeping keywords ModelCraft does not know"), - ready("solve.static", "Static", Solve, "Implicit", "Linear static solve of a step with ModelCraft's own solver"), + ready_more( + "solve.static", + "Static", + Solve, + "Implicit", + "Linear static (or coupled temperature–displacement) solve of a step with ModelCraft's own solver; heat transfer in the drop-down", + &["solve.heat"], + ), ready("solve.frequency", "Frequency", Solve, "Implicit", "Natural frequencies and mode shapes of a frequency step"), ready( "solve.explicit", diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index efd4da3..4db1f74 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -90,6 +90,7 @@ pub static ALL: &[CommandSpec] = &[ solve::RUN, solve::STATIC, solve::FREQUENCY, + solve::HEAT, solve::EXPLICIT, solve::HISTORY, solve::FLUID, diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 151ebbd..62079b6 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -5,7 +5,9 @@ use modelcraft_explicit::{ExplicitResult, solve_explicit}; use modelcraft_lbm::{FluidResult, Obstacle}; use modelcraft_mesh::surface::{BoundarySurface, SurfaceOptions}; use modelcraft_model::FluidDomain; -use modelcraft_solver::{FrequencyOptions, Limits, Location, MassKind, StaticResult, StepResult as ImplicitResult, solve_frequency, solve_static}; +use modelcraft_solver::{ + FrequencyOptions, Limits, Location, MassKind, StaticResult, StepResult as ImplicitResult, solve_frequency, solve_heat, solve_static, +}; use serde_json::{Value, json}; use super::{labels, opt, points, req}; @@ -22,11 +24,19 @@ pub const RUN: CommandSpec = CommandSpec { pub const STATIC: CommandSpec = CommandSpec { id: "solve.static", label: "Solve Static Step", - params: r#"{"step": "Load", "max_dofs": 6000000, "max_factor_entries": 400000000} — solves a linear static step (the only static step when `step` is left out) and keeps the result in the session; the caps may only be lowered"#, + params: r#"{"step": "Load", "max_dofs": 6000000, "max_factor_entries": 400000000} — solves a linear static or coupled temperature–displacement step (the only one when `step` is left out) and keeps the result in the session; the caps may only be lowered"#, kind: CommandKind::Model, run: run_static, }; +pub const HEAT: CommandSpec = CommandSpec { + id: "solve.heat", + label: "Solve Heat Transfer Step", + params: r#"{"step": "Heat", "max_dofs": 6000000, "max_factor_entries": 400000000} — steady or transient heat conduction of a heat transfer step (the only one when `step` is left out): temperatures, heat flux and the heat flow at prescribed temperatures, kept in the session; the caps may only be lowered"#, + kind: CommandKind::Model, + run: run_heat, +}; + pub const FREQUENCY: CommandSpec = CommandSpec { id: "solve.frequency", label: "Solve Frequency Step", @@ -97,8 +107,9 @@ fn limits(p: &Value) -> Result { }) } -fn largest(values: impl Iterator) -> Value { - match values.fold(None, |best: Option<(u32, f64)>, (l, v)| if best.is_none_or(|(_, b)| v > b) { Some((l, v)) } else { best }) { +/// The node with the value that `better` prefers, as `{"node", "value"}`. +fn extreme(values: impl Iterator, better: fn(f64, f64) -> bool) -> Value { + match values.fold(None, |best: Option<(u32, f64)>, (l, v)| if best.is_none_or(|(_, b)| better(v, b)) { Some((l, v)) } else { best }) { Some((node, value)) => json!({"node": node, "value": value}), None => Value::Null, } @@ -113,7 +124,8 @@ fn run_all(s: &mut Session, _: &Value) -> Result { .iter() .map(|st| { let solve: Solve = match st.procedure { - Procedure::Static {} => run_static, + Procedure::Static {} | Procedure::CoupledTemperatureDisplacement {} => run_static, + Procedure::HeatTransfer { .. } => run_heat, Procedure::Frequency { .. } => run_frequency, Procedure::ExplicitDynamic(_) => run_explicit, Procedure::Fluid { .. } => run_fluid, @@ -132,6 +144,14 @@ fn run_all(s: &mut Session, _: &Value) -> Result { Ok(json!({"steps": done})) } +fn largest(values: impl Iterator) -> Value { + extreme(values, |v, b| v > b) +} + +fn smallest(values: impl Iterator) -> Value { + extreme(values, |v, b| v < b) +} + fn run_static(s: &mut Session, p: &Value) -> Result { let limits = limits(p)?; let step = opt::(p, "step")?; @@ -152,12 +172,35 @@ fn run_static(s: &mut Session, p: &Value) -> Result { "displacement": p.displacement, "reaction": p.reaction, })).collect::>(), + "max_temperature": r.temperature.as_ref().map(|t| largest(r.nodes.iter().copied().zip(t.iter().copied()))), "warnings": r.warnings, }); s.store_result(StepResult::Implicit(ImplicitResult::Static(r))); Ok(summary) } +fn run_heat(s: &mut Session, p: &Value) -> Result { + let limits = limits(p)?; + let step = opt::(p, "step")?; + let r = solve_heat(&s.model, step.as_deref(), &limits)?; + let t = r.temperature(); + let summary = json!({ + "step": r.step, + "dofs": r.stats.dofs, + "equations": r.stats.equations, + "factor_entries": r.stats.factor_entries, + "frames": r.frames.len(), + "end_time": r.times.last(), + "min_temperature": smallest(r.nodes.iter().copied().zip(t.iter().copied())), + "max_temperature": largest(r.nodes.iter().copied().zip(t.iter().copied())), + "heat_supplied": r.applied_total, + "heat_removed": r.reaction_total, + "warnings": r.warnings, + }); + s.store_result(StepResult::Implicit(ImplicitResult::Heat(r))); + Ok(summary) +} + fn run_frequency(s: &mut Session, p: &Value) -> Result { let limits = limits(p)?; let step = opt::(p, "step")?; @@ -452,6 +495,7 @@ fn probe_implicit(result: &ImplicitResult, step: &str, nodes: &[u32], p: &Value) "reaction_moment": r.reaction_moment.get(i), "stress": r.stress.get(i), "von_mises": r.von_mises.get(i), + "temperature": r.temperature.as_ref().and_then(|t| t.get(i)), }), }) }) @@ -477,6 +521,17 @@ fn probe_implicit(result: &ImplicitResult, step: &str, nodes: &[u32], p: &Value) }) .collect::>()? } + ImplicitResult::Heat(r) => { + let (time, t) = crate::results::heat_frame(r, crate::results::Pick { mode: 1, frame: opt::(p, "frame")? })?; + let out = nodes + .iter() + .map(|&l| { + let i = r.nodes.binary_search(&l).map_err(|_| EngineError::Params(format!("node {l} has no temperature (it belongs to no solid element)")))?; + Ok(json!({"node": l, "temperature": t.get(i), "heat_flux": r.flux.get(i), "reaction": r.reaction.get(i)})) + }) + .collect::>>()?; + return Ok(json!({"step": step, "time": time, "nodes": out})); + } }; let mut response = json!({"step": step, "nodes": out}); if let ImplicitResult::Static(r) = result @@ -1172,4 +1227,40 @@ pub(crate) mod tests { assert!(s.execute("solve.probe", &p).is_err(), "{p}"); } } + + #[test] + fn solve_heat_and_probe_temperatures() { + // The cantilever bar as a conductor: 0 at the root, 100 at the tip, linear between. + let mut s = Session::default(); + cantilever(&mut s); + let script = json!([ + {"command": "model.material.edit", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 1000, "poisson": 0}, "density": 1, "thermal": {"conductivity": 2, "specific_heat": 3, "expansion": 1e-5}}}, + {"command": "model.boundary_condition.create", "params": {"name": "cold", "type": "temperature", "target": {"set": "root"}, "value": 0}}, + {"command": "model.boundary_condition.create", "params": {"name": "hot", "type": "temperature", "target": {"set": "tip"}, "value": 100}}, + {"command": "model.step.create", "params": {"name": "Heat", "procedure": {"type": "heat_transfer"}, "boundary_conditions": ["cold", "hot"]}}, + {"command": "model.step.create", "params": {"name": "Warm", "procedure": {"type": "heat_transfer", "transient": {"end_time": 10, "increment": 1}}, "boundary_conditions": ["cold", "hot"]}}, + ]); + s.run_script(&script).unwrap(); + let r = s.execute("solve.heat", &json!({"step": "Heat"})).unwrap(); + assert!((r["max_temperature"]["value"].as_f64().unwrap() - 100.0).abs() < 1e-9, "{r}"); + assert!(r["min_temperature"]["value"].as_f64().unwrap().abs() < 1e-9, "{r}"); + let probe = s.execute("solve.probe", &json!({"step": "Heat", "nodes": [6]})).unwrap(); + assert!((probe["nodes"][0]["temperature"].as_f64().unwrap() - 50.0).abs() < 1e-9, "{probe}"); + assert!((probe["nodes"][0]["heat_flux"][0].as_f64().unwrap() + 20.0).abs() < 1e-9, "{probe}"); + let w = s.execute("solve.heat", &json!({"step": "Warm"})).unwrap(); + assert_eq!(w["frames"], 11, "{w}"); + let early = s.execute("solve.probe", &json!({"step": "Warm", "nodes": [6], "frame": 1})).unwrap(); + assert_eq!(early["time"], 1.0); + assert!(s.execute("solve.probe", &json!({"step": "Warm", "nodes": [6], "frame": 99})).is_err()); + assert!(s.execute("solve.probe", &json!({"step": "Warm", "nodes": [999]})).is_err()); + assert!(s.execute("solve.heat", &json!({"step": "Bend"})).is_err(), "not a heat step"); + assert!(s.execute("solve.heat", &json!({"step": "Heat", "max_dofs": 0})).is_err()); + // The results display shows temperatures, undeformed. + let view = crate::results::ResultsView::for_result(&s.latest_result("Heat").unwrap().0.clone()); + assert_eq!((view.contour, view.deformed), (Some(crate::results::Field::Temperature), false)); + let (heat, _) = s.latest_result("Heat").unwrap(); + let fields: Vec<_> = crate::results::Field::available(heat).collect(); + assert_eq!(fields, vec![crate::results::Field::Temperature, crate::results::Field::HeatFlux]); + assert!(crate::results::displacements(heat, crate::results::Pick::default()).is_err()); + } } diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 934ffd1..8e8f103 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -428,6 +428,22 @@ mod tests { json!({"name": "A", "procedure": {"type": "fluid", "domain": "A", "end_time": -1, "steps": 3}}), json!({"name": "Bend", "procedure": {"type": "fluid", "domain": "nope", "end_time": 1}}), json!({"step": "Bend", "max_cells": 1e300, "max_steps": -3}), + // Heat transfer. + json!({"name": "A", "thermal": {"conductivity": -1}}), + json!({"name": "A", "thermal": {"conductivity": [1e308, f64::MIN_POSITIVE, 0]}, "density": 1}), + json!({"name": "A", "thermal": {"conductivity": "hot", "specific_heat": 0, "expansion": [1, 2]}}), + json!({"name": "A", "thermal": {}}), + json!({"name": "A", "thermal": {"expansion": 1e308, "reference_temperature": -1e308}}), + json!({"name": "A", "type": "temperature", "target": {"reference_point": "RP"}, "value": 1}), + json!({"name": "A", "type": "temperature", "target": {"set": "root"}, "value": "warm"}), + json!({"name": "A", "type": "surface_flux", "set": "end", "flux": 0}), + json!({"name": "A", "type": "convection", "set": "end", "film": -1, "sink": 1e308}), + json!({"name": "A", "type": "heat_source", "set": "all", "power": -1e308}), + json!({"name": "A", "procedure": {"type": "heat_transfer", "transient": {"end_time": 1e308, "increment": 1e-308}}}), + json!({"name": "A", "procedure": {"type": "heat_transfer", "transient": {"end_time": 1, "increment": 0.1, "theta": 0.1}}}), + json!({"name": "A", "procedure": {"type": "heat_transfer", "transient": {"end_time": -1, "increment": 0}}}), + json!({"name": "A", "procedure": {"type": "coupled_temperature_displacement", "extra": 1}}), + json!({"step": "Heat", "frame": u64::MAX, "nodes": [1]}), json!({"geometry": {"body": "Bar"}, "elements": [1, 99], "faces": [{"element": 1, "face": 200}]}), json!({"geometry": {"body": "Bar", "vertex": "v"}, "nodes": [1], "edges": [{"element": 1, "edge": 0}]}), ]; diff --git a/crates/engine/src/results.rs b/crates/engine/src/results.rs index 1212844..f38f43a 100644 --- a/crates/engine/src/results.rs +++ b/crates/engine/src/results.rs @@ -9,7 +9,7 @@ use modelcraft_explicit::{ExplicitResult, Frame}; use modelcraft_model::Model; -use modelcraft_solver::{FrequencyResult, StaticResult}; +use modelcraft_solver::{FrequencyResult, HeatResult, StaticResult}; use serde::{Deserialize, Serialize}; use serde_json::{Value, json}; @@ -22,6 +22,7 @@ pub enum Kind<'a> { Explicit(&'a ExplicitResult), /// Per cell, not per node: no nodal fields to display yet. Fluid(&'a modelcraft_lbm::FluidResult), + Heat(&'a HeatResult), } impl StepResult { @@ -29,6 +30,7 @@ impl StepResult { match self { StepResult::Implicit(modelcraft_solver::StepResult::Static(r)) => Kind::Static(r), StepResult::Implicit(modelcraft_solver::StepResult::Frequency(r)) => Kind::Frequency(r), + StepResult::Implicit(modelcraft_solver::StepResult::Heat(r)) => Kind::Heat(r), StepResult::Explicit(r) => Kind::Explicit(r), StepResult::Fluid(r) => Kind::Fluid(r), } @@ -43,6 +45,16 @@ pub struct Pick { pub frame: Option, } +/// The temperature frame of a heat result `pick` asks for. +pub fn heat_frame(r: &HeatResult, pick: Pick) -> Result<(f64, &[f64])> { + let last = r.frames.len().checked_sub(1).ok_or_else(|| EngineError::Params(format!("step `{}` has no temperatures", r.step)))?; + let k = pick.frame.unwrap_or(last); + match (r.times.get(k), r.frames.get(k)) { + (Some(t), Some(f)) => Ok((*t, f)), + _ => Err(EngineError::Params(format!("step `{}` has frames 0 to {last}; there is no frame {k}", r.step))), + } +} + /// The field frame of an explicit result `pick` asks for. pub fn frame(r: &ExplicitResult, pick: Pick) -> Result<&Frame> { let last = r.frames.len().checked_sub(1).ok_or_else(|| EngineError::Params(format!("step `{}` has no field frames", r.step)))?; @@ -71,10 +83,14 @@ pub enum Field { StressXy, StressXz, StressYz, + /// Temperature (heat transfer and coupled steps). + Temperature, + /// Heat flux magnitude (heat transfer steps, at the end of the step). + HeatFlux, } impl Field { - pub const ALL: [Field; 15] = [ + pub const ALL: [Field; 17] = [ Field::Displacement, Field::DisplacementX, Field::DisplacementY, @@ -90,6 +106,8 @@ impl Field { Field::StressXy, Field::StressXz, Field::StressYz, + Field::Temperature, + Field::HeatFlux, ]; /// The id commands take (`von_mises`…). @@ -110,6 +128,8 @@ impl Field { Field::StressXy => "stress_xy", Field::StressXz => "stress_xz", Field::StressYz => "stress_yz", + Field::Temperature => "temperature", + Field::HeatFlux => "heat_flux", } } @@ -130,6 +150,8 @@ impl Field { Field::StressXy => "Stress XY", Field::StressXz => "Stress XZ", Field::StressYz => "Stress YZ", + Field::Temperature => "Temperature", + Field::HeatFlux => "Heat flux magnitude", } } @@ -139,20 +161,46 @@ impl Field { /// Whether it is a stress. pub fn is_stress(self) -> bool { - !matches!(self, Field::Displacement | Field::DisplacementX | Field::DisplacementY | Field::DisplacementZ | Field::Velocity) + !matches!( + self, + Field::Displacement + | Field::DisplacementX + | Field::DisplacementY + | Field::DisplacementZ + | Field::Velocity + | Field::Temperature + | Field::HeatFlux + ) + } + + /// Whether it is a heat transfer quantity. + pub fn is_thermal(self) -> bool { + matches!(self, Field::Temperature | Field::HeatFlux) } /// The fields a result has: displacements everywhere, stresses for static steps and explicit /// frames that kept them (mode shapes have no stress scale), velocity for explicit frames. pub fn available(r: &StepResult) -> impl Iterator + '_ { - let (stress, velocity) = match r.kind() { - Kind::Static(_) => (true, false), - Kind::Frequency(_) => (false, false), - Kind::Explicit(e) => e.frames.last().map_or((false, false), |f| (!f.stress.is_empty(), !f.velocity.is_empty())), + let (stress, velocity, temperature, flux) = match r.kind() { + Kind::Static(s) => (true, false, s.temperature.is_some(), false), + Kind::Frequency(_) => (false, false, false, false), + Kind::Explicit(e) => { + e.frames.last().map_or((false, false, false, false), |f| (!f.stress.is_empty(), !f.velocity.is_empty(), false, false)) + } // Cell fields have their own display (M12 streamlines); no nodal field. - Kind::Fluid(_) => return None.into_iter().flatten(), + Kind::Fluid(_) => return Vec::new().into_iter(), + Kind::Heat(_) => return Field::ALL.into_iter().filter(|f| f.is_thermal()).collect::>().into_iter(), }; - Some(Field::ALL.into_iter().filter(move |f| if *f == Field::Velocity { velocity } else { stress || !f.is_stress() })).into_iter().flatten() + Field::ALL + .into_iter() + .filter(move |f| match f { + Field::Velocity => velocity, + Field::Temperature => temperature, + Field::HeatFlux => flux, + f => stress || !f.is_stress(), + }) + .collect::>() + .into_iter() } } @@ -190,11 +238,18 @@ impl ResultsView { /// (and an explicit step that kept stresses), the first mode with displacement contours for /// a frequency step; an explicit step at its last frame, to true scale. pub fn for_result(r: &StepResult) -> ResultsView { - let contour = if Field::available(r).any(|f| f == Field::VonMises) { Field::VonMises } else { Field::Displacement }; + let heat = matches!(r.kind(), Kind::Heat(_)); + let contour = if heat { + Field::Temperature + } else if Field::available(r).any(|f| f == Field::VonMises) { + Field::VonMises + } else { + Field::Displacement + }; let scale = matches!(r.kind(), Kind::Explicit(_)).then_some(1.0); ResultsView { step: r.step().to_string(), - deformed: true, + deformed: !heat, scale, contour: Some(contour), bands: DEFAULT_BANDS, @@ -228,6 +283,7 @@ pub fn displacements(r: &StepResult, pick: Pick) -> Result<&[[f64; 3]]> { Ok(&f.displacement) } Kind::Fluid(f) => Err(EngineError::Params(format!("step `{}` is a fluid step: its fields are per cell (solve.probe)", f.step))), + Kind::Heat(h) => Err(EngineError::Params(format!("step `{}` is heat transfer and has no displacements", h.step))), } } @@ -239,7 +295,7 @@ pub fn field_values(model: &Model, r: &StepResult, field: Field, pick: Pick) -> let stresses = || -> Result> { match r.kind() { Kind::Static(s) => Ok(std::borrow::Cow::Borrowed(s.stress.as_slice())), - Kind::Frequency(_) | Kind::Fluid(_) => Err(unavailable()), + Kind::Frequency(_) | Kind::Fluid(_) | Kind::Heat(_) => Err(unavailable()), Kind::Explicit(e) => { let f = frame(e, pick)?; if f.stress.is_empty() { @@ -261,7 +317,16 @@ pub fn field_values(model: &Model, r: &StepResult, field: Field, pick: Pick) -> let f = frame(e, pick)?; if f.velocity.is_empty() { Err(unavailable()) } else { Ok(f.velocity.iter().map(|v| magnitude(*v)).collect()) } } - Kind::Static(_) | Kind::Frequency(_) | Kind::Fluid(_) => Err(unavailable()), + Kind::Static(_) | Kind::Frequency(_) | Kind::Fluid(_) | Kind::Heat(_) => Err(unavailable()), + }, + Field::Temperature => match r.kind() { + Kind::Heat(h) => Ok(heat_frame(h, pick)?.1.to_vec()), + Kind::Static(StaticResult { temperature: Some(t), .. }) => Ok(t.clone()), + _ => Err(unavailable()), + }, + Field::HeatFlux => match r.kind() { + Kind::Heat(h) => Ok(h.flux.iter().map(|v| magnitude(*v)).collect()), + _ => Err(unavailable()), }, Field::VonMises => each(&modelcraft_fem::von_mises), Field::MaxPrincipal => principal(0), @@ -347,6 +412,12 @@ impl Session { Kind::Fluid(f) => { out["frames"] = f.frames.iter().map(|f| f.time).collect(); } + Kind::Heat(h) => { + out["frames"] = json!(h.times); + if let Ok((t, _)) = heat_frame(h, v.pick()) { + out["time"] = json!(t); + } + } Kind::Static(_) => {} } out @@ -390,7 +461,7 @@ mod tests { } let (modes, _) = s.latest_result("Modes").unwrap(); assert_eq!(Field::available(modes).count(), 4); - assert_eq!(Field::available(bend).count(), Field::ALL.len() - 1, "all but velocity"); + assert_eq!(Field::available(bend).count(), Field::ALL.len() - 3, "all but velocity and the thermal fields"); assert!(field_values(&s.model, modes, Field::VonMises, Pick { mode: 1, frame: None }).is_err()); assert!(field_values(&s.model, modes, Field::Displacement, Pick { mode: 9, frame: None }).is_err()); assert_eq!(field_values(&s.model, modes, Field::Displacement, Pick { mode: 2, frame: None }).unwrap().len(), modes.nodes().len()); diff --git a/crates/explicit/src/elements.rs b/crates/explicit/src/elements.rs index c2e3b2f..e483b51 100644 --- a/crates/explicit/src/elements.rs +++ b/crates/explicit/src/elements.rs @@ -311,6 +311,17 @@ fn update(mat: &Material, es: &mut ElementState, q: &Mat3, d: &[f64; 6], vol: f6 } } +/// The mechanical part of a strain increment: less the thermal strain increment `thermal` in +/// every normal direction (isotropic expansion, the same in any frame). +fn mechanical(mut d: [f64; 6], thermal: f64) -> [f64; 6] { + if thermal != 0.0 { + for x in d.iter_mut().take(3) { + *x -= thermal; + } + } + d +} + /// The density and dilatational wave speed of an element of rest volume `rest`, now of volume /// `volume`. fn current(mat: &Material, rest: f64, volume: f64) -> (f64, f64) { @@ -386,6 +397,7 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor let g = mid.gradient(&e.du); let q = rotation_of(sys, mat, es, &g, h.label, m.time)?; let (d, _) = rotation::split(&g); + let d = mechanical(d, h.thermal_rate * m.dt); let (rho, c) = current(mat, h.kinematics.volume, now.volume); update(mat, es, &q, &d, d[0] + d[1] + d[2], m.dt, rho, mid.volume); if let Some(failure) = es.failed { @@ -466,8 +478,10 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor let averaged = t.nodes.iter().map(|n| ws.theta.get(*n).copied().unwrap_or(0.0)).sum::() / 4.0; let q = rotation_of(sys, mat, es, &inc.g, t.label, m.time)?; let (d, _) = rotation::split(&inc.g); + let thermal = t.thermal_rate * m.dt; + let d = mechanical(d, thermal); let (rho, c) = current(mat, t.kinematics.volume, inc.now.volume); - update(mat, es, &q, &d, averaged, m.dt, rho, inc.mid.volume); + update(mat, es, &q, &d, averaged - 3.0 * thermal, m.dt, rho, inc.mid.volume); if let Some(failure) = es.failed { state.eroded.push((t.index, failure)); continue; diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index b2b5428..cac6fd7 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -39,6 +39,9 @@ pub struct Material { pub j2: Option, /// The equivalent plastic strain at which elements erode. pub erosion: Option, + /// The (isotropic) expansion coefficient and reference temperature, when the material + /// expands. + pub expansion: Option<(f64, f64)>, } /// A one-point hexahedron. @@ -58,6 +61,8 @@ pub struct Hex { /// The element's stable length: its critical increment without bulk viscosity is /// `length / wave_speed`. pub length: f64, + /// The rate of its thermal strain, α ΔT / end time (every normal direction). + pub thermal_rate: f64, } /// A nodal-pressure tetrahedron. @@ -72,6 +77,8 @@ pub struct Tet { /// The frequency bound ω² of the undeformed element. pub omega2: f64, pub length: f64, + /// The rate of its thermal strain, α ΔT / end time (every normal direction). + pub thermal_rate: f64, } /// A fixed, frictionless rigid plane. @@ -145,8 +152,14 @@ fn material(model: &Model, name: &str) -> Result { } Elastic::EngineeringConstants { .. } => None, }; + let expansion = match m.thermal.and_then(|t| t.expansion.map(|a| (a.axes(), t.reference_temperature))) { + None => None, + Some(([a, b, c], t_ref)) if a == b && b == c => Some((a, t_ref)), + Some(_) => return Err(Error::NotSupported(format!("orthotropic thermal expansion (material `{name}`) in explicit dynamics"))), + }; Ok(Material { name: name.to_string(), + expansion, elasticity: d, density, wave_speed: (stiffest / density).sqrt(), @@ -278,7 +291,7 @@ impl System { let rest = hex_scales(&k, &gamma, &mass, mat, rho, settings.hourglass); sys.lump(&nodes, &mass); let length = 2.0 * c / rest.omega2.sqrt(); - sys.hexes.push(Hex { label, index, nodes, material: mi, kinematics: k, gamma, mass, rest, length }); + sys.hexes.push(Hex { label, index, nodes, material: mi, kinematics: k, gamma, mass, rest, length, thermal_rate: 0.0 }); } ElementKind::Tet4 => { let nodes: [usize; 4] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 4 nodes")))?; @@ -287,7 +300,17 @@ impl System { let mass = k.nodal_volume.map(|v| rho * v); let omega2 = frequency_bound(&k, k.gradient_norm(), mat.modulus, &mass); sys.lump(&nodes, &mass); - sys.tets.push(Tet { label, index, nodes, material: mi, kinematics: k, mass, omega2, length: 2.0 * c / omega2.sqrt() }); + sys.tets.push(Tet { + label, + index, + nodes, + material: mi, + kinematics: k, + mass, + omega2, + length: 2.0 * c / omega2.sqrt(), + thermal_rate: 0.0, + }); } ElementKind::Rod2 | ElementKind::Rod3 => { let area = area.ok_or_else(|| Error::Internal("truss area".into()))?; @@ -319,6 +342,7 @@ impl System { sys.discrete_increment = sys.discrete.critical_increment(&sys.dof_mass(), sys.nonlinear); sys.prescribe(model, resolved, step, settings.end_time)?; sys.load(model, resolved, step)?; + sys.temperatures(model, resolved, step, settings.end_time)?; sys.conditions(resolved, settings)?; Ok(sys) } @@ -355,6 +379,9 @@ impl System { if let Target::ReferencePoint(rp) = bc.target() { return Err(Error::NotSupported(format!("boundary conditions on reference points (`{name}` on `{rp}`)"))); } + if bc.is_thermal() { + return Err(Error::NotSupported(format!("temperature boundary conditions (`{name}`) in explicit dynamics"))); + } for &label in resolved.target_nodes(bc.target()) { let Some(node) = self.node(label) else { continue }; for kind in self.layout.kinds(node).iter() { @@ -380,6 +407,51 @@ impl System { Ok(()) } + /// The thermal strain rates of the step's temperature loads: each element's mean nodal + /// temperature change from its material's reference (nodes no load reaches stay at it) is + /// reached linearly over the step, like a prescribed displacement, so the thermal strain + /// α ΔT grows at the rate α ΔT / end time (no thermal shock at the start). + fn temperatures(&mut self, model: &Model, resolved: &Resolved, step: &Step, end_time: f64) -> Result<()> { + let mut field: Vec> = vec![None; self.labels.len()]; + let mut any = false; + for name in &step.loads { + let Some(Load::Temperature { target, value }) = model.load(name) else { continue }; + any = true; + for &label in resolved.target_nodes(target) { + let Some(slot) = self.node(label).and_then(|n| field.get_mut(n)) else { continue }; + match *slot { + Some(old) if old != *value => { + return Err(invalid(format!("node {label} is given both {old} and {value} (temperature load `{name}` and another)"))); + } + _ => *slot = Some(*value), + } + } + } + if !any { + return Ok(()); + } + if !self.rods.is_empty() { + self.warnings + .push("trusses take no thermal strain in explicit dynamics yet; the temperature loads strain the solid elements only".into()); + } + let materials = &self.materials; + let rate = |nodes: &[usize], material: usize| -> f64 { + let Some((alpha, reference)) = materials.get(material).and_then(|m| m.expansion) else { return 0.0 }; + let dt: f64 = nodes.iter().map(|n| field.get(*n).copied().flatten().map_or(0.0, |t| t - reference)).sum::() / nodes.len() as f64; + alpha * dt / end_time + }; + for h in &mut self.hexes { + h.thermal_rate = rate(&h.nodes, h.material); + } + for t in &mut self.tets { + t.thermal_rate = rate(&t.nodes, t.material); + } + if self.hexes.iter().map(|h| h.thermal_rate).chain(self.tets.iter().map(|t| t.thermal_rate)).any(|r| !r.is_finite()) { + return Err(invalid("the thermal strain rate is not finite")); + } + Ok(()) + } + /// The constant load vector: forces and tractions as the model distributes them, consistent /// pressure loads, gravity and body forces with the lumped (row-sum) weights. fn load(&mut self, model: &Model, resolved: &Resolved, step: &Step) -> Result<()> { @@ -394,6 +466,11 @@ impl System { Load::Moment { .. } => { return Err(Error::NotSupported(format!("moments (`{name}`): solid nodes have no rotations"))); } + Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. } => { + return Err(Error::NotSupported(format!("{} loads (`{name}`) in explicit dynamics", load.kind()))); + } + // Thermal strains, set by `temperatures`. + Load::Temperature { .. } => {} Load::Force { .. } | Load::Traction { .. } => { for (label, v) in resolved.nodal_loads.get(name).map(Vec::as_slice).unwrap_or_default() { let node = self.node(*label).ok_or_else(|| off_mesh(*label))?; diff --git a/crates/explicit/tests/thermal.rs b/crates/explicit/tests/thermal.rs new file mode 100644 index 0000000..cc3e4f8 --- /dev/null +++ b/crates/explicit/tests/thermal.rs @@ -0,0 +1,140 @@ +//! Thermal strain in explicit dynamics: a temperature load reached slowly over the step expands +//! a free bar by α ΔT L without stress; a bar held at both ends carries σ = −E α ΔT, and with +//! perfectly plastic J2 the stress stops at the yield stress. + +use modelcraft_explicit::{Limits, solve_explicit}; +use modelcraft_model::*; + +const YOUNG: f64 = 1.0e4; +const ALPHA: f64 = 1.0e-4; +const DT: f64 = 20.0; +const LENGTH: f64 = 4.0; + +/// A 4 × 1 × 1 bar of hex8 (one-point) or Kuhn-split tet4 cells, ν = 0, heated by DT over a +/// step much longer than the wave transit, held at x = 0 (and at x = L when `clamped`) along x +/// and on the planes y = 0, z = 0. +fn bar(kind: ElementKind, clamped: bool, plastic: Option) -> Model { + let mut m = Model::default(); + let n = 4usize; + let label = |i: usize, j: usize, k: usize| (1 + i + (n + 1) * (j + 2 * k)) as u32; + for k in 0..2 { + for j in 0..2 { + for i in 0..=n { + m.set_node(label(i, j, k), [i as f64, j as f64, k as f64]).unwrap(); + } + } + } + let corners = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]; + let mut next = 1; + for i in 0..n { + let at = |c: [usize; 3]| label(i + c[0], c[1], c[2]); + if kind == ElementKind::Hex8 { + m.set_element(next, kind, corners.iter().map(|c| at(*c)).collect()).unwrap(); + next += 1; + continue; + } + for perm in [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]] { + let mut v = [[0usize; 3]; 4]; + for s in 0..3 { + v[s + 1] = v[s]; + v[s + 1][perm[s]] = 1; + } + let x = v.map(|c| c.map(|k| k as f64)); + let e = |a: usize| [0, 1, 2].map(|k| x[a][k] - x[0][k]); + let (a, b, c) = (e(1), e(2), e(3)); + let vol = a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]); + if vol < 0.0 { + v.swap(1, 2); + } + m.set_element(next, kind, v.iter().map(|c| at(*c)).collect()).unwrap(); + next += 1; + } + } + let all: Vec = m.elements().keys().copied().collect(); + let nodes: Vec = m.nodes().keys().copied().collect(); + let pick = |f: &dyn Fn([f64; 3]) -> bool| -> Vec { nodes.iter().copied().filter(|l| f(m.node(*l).unwrap().xyz)).collect() }; + let (x0, x1, y0, z0) = (pick(&|p| p[0] == 0.0), pick(&|p| p[0] == LENGTH), pick(&|p| p[1] == 0.0), pick(&|p| p[2] == 0.0)); + m.add_set("all", Set { elements: all, nodes: nodes.clone(), ..Set::default() }).unwrap(); + for (name, set) in [("x0", x0), ("x1", x1), ("y0", y0), ("z0", z0)] { + m.add_set(name, Set { nodes: set, ..Set::default() }).unwrap(); + } + let thermal = Thermal { expansion: Some(Directional::Isotropic(ALPHA)), ..Thermal::default() }; + m.add_material( + "Mat", + Material { + elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), + density: Some(1.0), + thermal: Some(thermal), + plastic, + ..Material::default() + }, + ) + .unwrap(); + let integration = if kind == ElementKind::Hex8 { Integration::Reduced } else { Integration::Full }; + let formulation = Formulation { integration, ..Formulation::default() }; + m.add_section("bar", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap(); + let mut step = Step::new("Heat", Procedure::ExplicitDynamic(ExplicitDynamic::new(20.0))); + let mut bcs = vec![("x0", Symmetry::X), ("y0", Symmetry::Y), ("z0", Symmetry::Z)]; + if clamped { + bcs.push(("x1", Symmetry::X)); + } + for (set, plane) in bcs { + let name = format!("hold_{set}"); + m.add_boundary_condition(&name, BoundaryCondition::Symmetry { target: Target::Set(set.into()), plane }).unwrap(); + step.boundary_conditions.insert(name); + } + m.add_load("warm", Load::Temperature { target: Target::Set("all".into()), value: DT }).unwrap(); + step.loads.insert("warm".into()); + m.add_step(step, None).unwrap(); + m +} + +#[test] +fn a_free_bar_expands_without_stress() { + for kind in [ElementKind::Hex8, ElementKind::Tet4] { + let r = solve_explicit(&bar(kind, false, None), None, &Limits::default()).unwrap(); + let last = r.frames.last().unwrap(); + let tip = r.nodes.iter().position(|l| *l == 5).unwrap(); + let want = ALPHA * DT * LENGTH; + assert!((last.displacement[tip][0] - want).abs() < 0.02 * want, "{kind:?}: {:?} vs {want}", last.displacement[tip]); + let worst = last.stress.iter().flatten().fold(0.0f64, |m, s| m.max(s.abs())); + assert!(worst < 0.02 * YOUNG * ALPHA * DT, "{kind:?}: {worst}"); + } +} + +#[test] +fn a_clamped_bar_carries_minus_e_alpha_dt_up_to_yield() { + let sigma = -YOUNG * ALPHA * DT; + for kind in [ElementKind::Hex8, ElementKind::Tet4] { + let r = solve_explicit(&bar(kind, true, None), None, &Limits::default()).unwrap(); + for s in &r.frames.last().unwrap().stress { + assert!((s[0] - sigma).abs() < 0.02 * sigma.abs(), "{kind:?}: {s:?}"); + } + // Perfectly plastic at half that stress: the bar yields and the stress stays at −σ_y. + let plastic = Plastic::Bilinear { yield_stress: 0.5 * sigma.abs(), hardening_modulus: 0.0, hardening: Hardening::Isotropic }; + let r = solve_explicit(&bar(kind, true, Some(plastic)), None, &Limits::default()).unwrap(); + let last = r.frames.last().unwrap(); + for (s, ep) in last.stress.iter().zip(&last.plastic_strain) { + let mises = modelcraft_fem::von_mises(s); + assert!((mises - 0.5 * sigma.abs()).abs() < 0.02 * sigma.abs(), "{kind:?}: {s:?}"); + assert!(*ep > 0.0, "{kind:?}: plastic strain {ep}"); + } + } +} + +#[test] +fn orthotropic_expansion_and_heat_loads_are_refused() { + let mut m = bar(ElementKind::Hex8, false, None); + let thermal = Thermal { expansion: Some(Directional::Orthotropic([1e-5, 2e-5, 1e-5])), ..Thermal::default() }; + m.replace_material( + "Mat", + Material { + elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), + density: Some(1.0), + thermal: Some(thermal), + ..Material::default() + }, + ) + .unwrap(); + assert!(solve_explicit(&m, None, &Limits::default()).unwrap_err().to_string().contains("orthotropic thermal expansion")); +} diff --git a/crates/fem/src/element.rs b/crates/fem/src/element.rs index 9b3bf4b..790c192 100644 --- a/crates/fem/src/element.rs +++ b/crates/fem/src/element.rs @@ -195,6 +195,26 @@ struct Gradients { det: f64, } +/// A thermal strain field over an element: the expansion coefficients along X, Y, Z and the +/// nodal temperature changes T − T_ref. +#[derive(Debug, Clone, PartialEq)] +pub struct ThermalStrain { + pub alpha: [f64; 3], + pub delta: Vec, +} + +impl ThermalStrain { + pub(crate) fn check(&self, nodes: usize) -> Result<()> { + if self.delta.len() != nodes { + return Err(Error::Invalid(format!("{} temperatures for a {nodes}-node element", self.delta.len()))); + } + if self.alpha.iter().chain(&self.delta).any(|v| !v.is_finite()) { + return Err(Error::Invalid("temperatures and expansion coefficients must be finite".into())); + } + Ok(()) + } +} + /// A solid element: its type, its node coordinates and its formulation. #[derive(Debug, Clone, Copy)] pub struct Solid<'a> { @@ -423,16 +443,83 @@ impl<'a> Solid<'a> { Ok(f) } + /// The thermal strain α ΔT at a point with shape values `n` (xx, yy, zz; no shear). + fn thermal_strain(thermal: &ThermalStrain, n: &[f64]) -> [f64; 6] { + let dt: f64 = n.iter().zip(&thermal.delta).map(|(a, b)| a * b).sum(); + [thermal.alpha[0] * dt, thermal.alpha[1] * dt, thermal.alpha[2] * dt, 0.0, 0.0, 0.0] + } + + /// The equivalent nodal forces of a thermal strain, ∫ Bᵀ D ε_th dV, with the stiffness rule + /// (the initial-strain load of Zienkiewicz, Taylor & Zhu 2013, §6.2.3); incompatible modes + /// are condensed out as in the stiffness, f = f_u − K_ua K_aa⁻¹ f_a. + pub fn thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> Result> { + thermal.check(self.ty.nodes)?; + if self.formulation == Formulation::IncompatibleModes { + let IncompatibleParts { kua, kaa, j0, .. } = self.incompatible_parts(d, false)?; + let (mut f, mut fa) = self.thermal_forces(d, thermal, Some(&j0))?; + Cholesky::new(&kaa)?.solve(&mut fa); + for (i, fi) in f.iter_mut().flatten().enumerate() { + *fi -= (0..INCOMPATIBLE_DOFS).map(|a| kua.get(i * INCOMPATIBLE_DOFS + a).copied().unwrap_or(0.0) * fa[a]).sum::(); + } + return Ok(f); + } + Ok(self.thermal_forces(d, thermal, None)?.0) + } + + /// ∫ Bᵀ D ε_th dV on the nodes and, with the centre Jacobian `j0`, on the incompatible modes. + fn thermal_forces( + &self, + d: &Elasticity, + thermal: &ThermalStrain, + j0: Option<&([[f64; 3]; 3], f64)>, + ) -> Result<(Vec<[f64; 3]>, [f64; INCOMPATIBLE_DOFS])> { + let mut f = vec![[0.0; 3]; self.ty.nodes]; + let mut fa = [0.0; INCOMPATIBLE_DOFS]; + let add = |out: &mut [f64], g: &[f64; 3], sig: &[f64; 6], w: f64| { + let b = b_block(g); + for (c, o) in out.iter_mut().enumerate() { + *o += w * (0..6).map(|s| b[s][c] * sig[s]).sum::(); + } + }; + for p in self.rule().points() { + let g = self.gradients(p.xi)?; + let w = g.det * p.weight; + let sig = stress(d, &Self::thermal_strain(thermal, g.shape.values())); + for (fi, gi) in f.iter_mut().zip(g.dx.get(..self.ty.nodes).unwrap_or_default()) { + add(fi, gi, &sig, w); + } + if let Some(j0) = j0 { + for (k, gk) in self.incompatible_gradients(p.xi, g.det, j0).iter().enumerate() { + add(fa.get_mut(3 * k..3 * k + 3).unwrap_or_default(), gk, &sig, w); + } + } + } + Ok((f, fa)) + } + /// Strain and stress at every integration point for the nodal displacements `u`. pub fn stresses(&self, d: &Elasticity, u: &[[f64; 3]]) -> Result> { + self.stresses_with(d, u, None) + } + + /// Strain and stress at every integration point for the nodal displacements `u` and, if + /// given, a thermal strain: the stress is D (ε − ε_th), the strain reported the total ε. + pub fn stresses_with(&self, d: &Elasticity, u: &[[f64; 3]], thermal: Option<&ThermalStrain>) -> Result> { if u.len() != self.ty.nodes { return Err(Error::Invalid(format!("{} displacements for a {}-node element", u.len(), self.ty.nodes))); } - // Incompatible modes: recover the internal amplitudes α = −K_aa⁻¹ K_au u. + if let Some(t) = thermal { + t.check(self.ty.nodes)?; + } + // Incompatible modes: recover the internal amplitudes from K_aa α + K_au u = f_a (the + // thermal load on the modes; zero without one). let alpha = if self.formulation == Formulation::IncompatibleModes { let IncompatibleParts { kua, kaa, j0, .. } = self.incompatible_parts(d, false)?; let flat: Vec = u.iter().flatten().copied().collect(); - let mut rhs = [0.0; INCOMPATIBLE_DOFS]; + let mut rhs = match thermal { + Some(t) => self.thermal_forces(d, t, Some(&j0))?.1, + None => [0.0; INCOMPATIBLE_DOFS], + }; for (i, ui) in flat.iter().enumerate() { for (a, r) in rhs.iter_mut().enumerate() { *r -= kua.get(i * INCOMPATIBLE_DOFS + a).copied().unwrap_or(0.0) * ui; @@ -471,7 +558,14 @@ impl<'a> Solid<'a> { xyz[k] += n * x[k]; } } - out.push(PointStress { xyz, strain, stress: stress(d, &strain) }); + let sigma = match thermal { + Some(t) => { + let eth = Self::thermal_strain(t, g.shape.values()); + stress(d, &[0, 1, 2, 3, 4, 5].map(|k| strain[k] - eth[k])) + } + None => stress(d, &strain), + }; + out.push(PointStress { xyz, strain, stress: sigma }); } Ok(out) } diff --git a/crates/fem/src/heat.rs b/crates/fem/src/heat.rs new file mode 100644 index 0000000..53b17f2 --- /dev/null +++ b/crates/fem/src/heat.rs @@ -0,0 +1,233 @@ +//! Heat conduction kernels: conductivity, heat capacity, internal heat sources, surface fluxes and +//! convection, and the heat flux at the integration points. +//! +//! The Galerkin form of transient conduction ρc Ṫ = ∇·(k∇T) + Q with flux and convection on the +//! boundary (Zienkiewicz, Taylor & Zhu 2013, §3.2 and ch. 7 of vol. 1; Hughes 2000, ch. 1, §8.1) +//! gives, element by element, with one temperature per node ([`ElementDofs`] `TEMPERATURE`): +//! - conductivity K = ∫ ∇Nᵀ k ∇N dV, k the diagonal conductivity in the global axes (the +//! material axes until orientations land), integrated with the element's full rule; +//! - capacity C = ∫ ρc N Nᵀ dV with the mass rule (consistent); +//! - source f = ∫ N Q dV; +//! - on faces: a flux q into the body gives f = ∫ N q dA; convection with film coefficient h to a +//! sink at T∞ adds h ∫ N Nᵀ dA to K and h T∞ ∫ N dA to f. +//! +//! Faces are integrated on their own (possibly curved) geometry with the model's face rule +//! (the one pressures use), exact for N Nᵀ on flat faces. + +use modelcraft_model::ElementKind; + +use crate::basis::Basis; +use crate::dense::Matrix; +use crate::dof::{DofSet, ElementDofs}; +use crate::element::{ElementType, element_type}; +use crate::quadrature::Rule; +use crate::{Error, Result}; + +/// What a heat conduction element carries: one temperature per node. +pub const HEAT_DOFS: ElementDofs = ElementDofs::new(DofSet::TEMPERATURE); + +/// The heat flux −k∇T at one integration point. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct PointFlux { + pub xyz: [f64; 3], + pub flux: [f64; 3], +} + +/// A solid element conducting heat. +#[derive(Debug, Clone, Copy)] +pub struct Conduction<'a> { + pub ty: &'static ElementType, + basis: Basis, + xyz: &'a [[f64; 3]], +} + +impl<'a> Conduction<'a> { + pub fn new(kind: ElementKind, xyz: &'a [[f64; 3]]) -> Result> { + let ty = element_type(kind)?; + if xyz.len() != ty.nodes { + return Err(Error::Invalid(format!("a {} element has {} nodes, got {}", kind.name(), ty.nodes, xyz.len()))); + } + modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?; + Ok(Conduction { ty, basis: Basis::of(kind), xyz }) + } + + /// Visit the points of `rule`: shape values, physical gradients and the weight w det J. + fn each(&self, rule: Rule, mut visit: impl FnMut(&[f64], &[[f64; 3]], f64)) -> Result<()> { + for p in rule.points() { + let s = self.basis.eval(p.xi); + let g = self.basis.gradients(&s, self.xyz)?; + visit(s.values(), g.gradients(), g.det * p.weight); + } + Ok(()) + } + + /// K = ∫ ∇Nᵀ k ∇N dV, n × n. + pub fn conductivity(&self, k: [f64; 3]) -> Result { + let mut m = Matrix::zeros(self.ty.nodes); + self.each(self.ty.full, |_, g, w| { + for (i, gi) in g.iter().enumerate() { + for (j, gj) in g.iter().enumerate() { + m.add(i, j, w * (k[0] * gi[0] * gj[0] + k[1] * gi[1] * gj[1] + k[2] * gi[2] * gj[2])); + } + } + })?; + Ok(m) + } + + /// C = ∫ ρc N Nᵀ dV, n × n; `capacity` is ρc (energy per volume per temperature). + pub fn capacity(&self, capacity: f64) -> Result { + let mut m = Matrix::zeros(self.ty.nodes); + self.each(self.ty.mass, |n, _, w| { + for (i, a) in n.iter().enumerate() { + for (j, b) in n.iter().enumerate() { + m.add(i, j, w * capacity * a * b); + } + } + })?; + Ok(m) + } + + /// f = ∫ N Q dV for a uniform heat source Q per volume. + pub fn source(&self, q: f64) -> Result> { + let mut f = vec![0.0; self.ty.nodes]; + self.each(self.ty.mass, |n, _, w| { + for (fi, ni) in f.iter_mut().zip(n) { + *fi += w * q * ni; + } + })?; + Ok(f) + } + + /// The heat flux −k∇T at the full rule's points for nodal temperatures `t`. + pub fn flux(&self, k: [f64; 3], t: &[f64]) -> Result> { + if t.len() != self.ty.nodes { + return Err(Error::Invalid(format!("{} temperatures for a {}-node element", t.len(), self.ty.nodes))); + } + let mut out = Vec::new(); + self.each(self.ty.full, |n, g, _| { + let (mut grad, mut xyz) = ([0.0; 3], [0.0; 3]); + for ((gi, ti), (ni, xi)) in g.iter().zip(t).zip(n.iter().zip(self.xyz)) { + for c in 0..3 { + grad[c] += gi[c] * ti; + xyz[c] += ni * xi[c]; + } + } + out.push(PointFlux { xyz, flux: [0, 1, 2].map(|c| -k[c] * grad[c]) }); + })?; + Ok(out) + } +} + +/// Visit the points of a face (node coordinates in the face's node order): shape values and +/// the area weight, with the model's face rule ([`modelcraft_model::nodal::integrate_face`], the +/// one pressures and tractions use). +fn each_face_point(xyz: &[[f64; 3]], mut visit: impl FnMut(&[f64], f64)) -> Result<()> { + let mut collapsed = false; + modelcraft_model::nodal::integrate_face(xyz, |n, nda| { + let area = (nda[0] * nda[0] + nda[1] * nda[1] + nda[2] * nda[2]).sqrt(); + collapsed |= !(area > 0.0 && area.is_finite()); + visit(n, area); + }) + .map_err(|e| Error::Invalid(e.to_string()))?; + if collapsed { + return Err(Error::Degenerate("a face is collapsed (zero area)".into())); + } + Ok(()) +} + +/// ∫ N dA for every node of a face: the nodal shares of a uniform flux. +pub fn face_weights(xyz: &[[f64; 3]]) -> Result> { + let mut f = vec![0.0; xyz.len()]; + each_face_point(xyz, |n, w| { + for (fi, ni) in f.iter_mut().zip(n) { + *fi += w * ni; + } + })?; + Ok(f) +} + +/// h ∫ N Nᵀ dA for a face: the conductivity convection adds. +pub fn face_film(xyz: &[[f64; 3]], film: f64) -> Result { + let mut m = Matrix::zeros(xyz.len()); + each_face_point(xyz, |n, w| { + for (i, a) in n.iter().enumerate() { + for (j, b) in n.iter().enumerate() { + m.add(i, j, w * film * a * b); + } + } + })?; + Ok(m) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::shape::{HEX_CORNERS, TET_CORNERS}; + + fn box_hex(a: [f64; 3]) -> Vec<[f64; 3]> { + HEX_CORNERS.iter().map(|c| [0, 1, 2].map(|k| (c[k] + 1.0) / 2.0 * a[k])).collect() + } + + #[test] + fn conductivity_rows_sum_to_zero_and_reproduce_a_linear_field() { + // A uniform gradient G in a 2 × 1 × 3 box: K T = boundary flux, zero at interior + // balance; here check K·1 = 0 and Tᵀ K T = k |G|² V. + let xyz = box_hex([2.0, 1.0, 3.0]); + for kind in [ElementKind::Hex8] { + let e = Conduction::new(kind, &xyz).unwrap(); + let k = [2.0, 3.0, 5.0]; + let m = e.conductivity(k).unwrap(); + let ones = vec![1.0; 8]; + assert!(m.mul_vec(&ones).iter().all(|v| v.abs() < 1e-12)); + let g = [0.3, -0.7, 0.2]; + let t: Vec = xyz.iter().map(|x| g[0] * x[0] + g[1] * x[1] + g[2] * x[2]).collect(); + let kt = m.mul_vec(&t); + let energy: f64 = kt.iter().zip(&t).map(|(a, b)| a * b).sum(); + let exact = (k[0] * g[0] * g[0] + k[1] * g[1] * g[1] + k[2] * g[2] * g[2]) * 6.0; + assert!((energy - exact).abs() < 1e-12 * exact, "{energy} vs {exact}"); + let fl = e.flux(k, &t).unwrap(); + for p in fl { + for c in 0..3 { + assert!((p.flux[c] + k[c] * g[c]).abs() < 1e-13); + } + } + } + } + + #[test] + fn capacity_and_source_integrate_the_volume() { + let tet: Vec<[f64; 3]> = TET_CORNERS.iter().map(|c| [2.0 * c[0], 3.0 * c[1], c[2]]).collect(); + let e = Conduction::new(ElementKind::Tet4, &tet).unwrap(); + let v = 1.0; + let c = e.capacity(4.0).unwrap(); + assert!((c.data.iter().sum::() - 4.0 * v).abs() < 1e-13); + assert!((c.get(0, 0) - 4.0 * v / 10.0).abs() < 1e-13 && (c.get(0, 1) - 4.0 * v / 20.0).abs() < 1e-13); + let f = e.source(3.0).unwrap(); + assert!(f.iter().all(|x| (x - 0.75).abs() < 1e-13), "{f:?}"); + assert!(Conduction::new(ElementKind::Tet4, &tet[..3]).is_err()); + assert!(e.flux([1.0; 3], &[0.0; 3]).is_err()); + } + + #[test] + fn face_weights_and_film() { + let quad = [[0., 0., 0.], [2., 0., 0.], [2., 3., 0.], [0., 3., 0.]]; + let w = face_weights(&quad).unwrap(); + assert!(w.iter().all(|x| (x - 1.5).abs() < 1e-14)); + let m = face_film(&quad, 2.0).unwrap(); + // h A / 9 on the diagonal of a bilinear rectangle, sum h A. + assert!((m.data.iter().sum::() - 12.0).abs() < 1e-12); + assert!((m.get(0, 0) - 2.0 * 6.0 / 9.0).abs() < 1e-13); + let mut t6 = vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]]; + for (a, b) in [(0, 1), (1, 2), (2, 0)] { + t6.push([0, 1, 2].map(|k| (t6[a][k] + t6[b][k]) / 2.0)); + } + let w = face_weights(&t6).unwrap(); + for (i, x) in w.iter().enumerate() { + let want = if i < 3 { 0.0 } else { 0.5 / 3.0 }; + assert!((x - want).abs() < 1e-14, "{i}: {x}"); + } + assert!(face_weights(&quad[..2]).is_err()); + assert!(face_weights(&[[0.0; 3]; 4]).is_err()); + assert!(face_film(&[[f64::NAN, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], 1.0).is_err()); + } +} diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index d5352e3..f0ddcdc 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -22,6 +22,8 @@ //! - [`beam`]: two-node Timoshenko beams with profile section properties; //! - [`frame`]: local axes of shells and beams and the rotation of their matrices; //! - [`face`]: consistent nodal forces of pressure on element faces; +//! - [`heat`]: heat conduction (conductivity, capacity, sources, surface flux and convection, +//! heat flux); thermal strain loads and stresses are in [`element`]; //! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass //! shape vectors and element frequency bounds for explicit dynamics; //! - [`rotation`]: incrementally objective stress updates (Hughes–Winget Jaumann, Flanagan–Taylor @@ -47,6 +49,7 @@ pub mod dof; pub mod element; pub mod face; pub mod frame; +pub mod heat; pub mod jet; pub mod material; pub mod plane; @@ -60,7 +63,7 @@ pub mod uniform; pub use basis::Basis; pub use dof::{DofKind, DofLayout, DofSet, ElementDofs}; -pub use element::{ELEMENT_TYPES, ElementType, Family, Formulation, PointStress, Solid, element_type}; +pub use element::{ELEMENT_TYPES, ElementType, Family, Formulation, PointStress, Solid, ThermalStrain, element_type}; pub use material::{Elasticity, elasticity, principal, von_mises}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] diff --git a/crates/fem/src/plane.rs b/crates/fem/src/plane.rs index e925b2b..df23620 100644 --- a/crates/fem/src/plane.rs +++ b/crates/fem/src/plane.rs @@ -44,7 +44,7 @@ use modelcraft_model::{ElementKind, Extent, SectionKind}; use crate::basis::Basis; use crate::dense::{Cholesky, Matrix}; use crate::dof::{DofSet, ElementDofs, Internal}; -use crate::element::{ElementType, Family, Formulation, PointStress, extrapolation, lump, symmetrize}; +use crate::element::{ElementType, Family, Formulation, PointStress, ThermalStrain, extrapolation, lump, symmetrize}; use crate::material::Elasticity; use crate::quadrature::{Point, Rule}; use crate::shape::{MAX_NODES, Shape}; @@ -471,20 +471,98 @@ impl<'a> Plane<'a> { Ok(f) } + /// The thermal strain α ΔT at a point with shape values `n` (xx, yy, zz, xy; zz is the hoop + /// direction of an axisymmetric element). + fn thermal_strain(thermal: &ThermalStrain, n: &[f64]) -> [f64; 4] { + let dt: f64 = n.iter().zip(&thermal.delta).map(|(a, b)| a * b).sum(); + [thermal.alpha[0] * dt, thermal.alpha[1] * dt, thermal.alpha[2] * dt, 0.0] + } + + /// ∫ Bᵀ D ε_th dV on the nodes' (u, v) and on the incompatible modes. + fn thermal_forces( + d2: &Elasticity2, + thermal: &ThermalStrain, + points: &[(Point, At)], + modes: Option<&Vec<[Column; 2]>>, + nodes: usize, + ) -> (Vec, Vec) { + let mut fu = vec![0.0; 2 * nodes]; + let mut fa = vec![0.0; INCOMPATIBLE_DOFS]; + for (i, (_, a)) in points.iter().enumerate() { + let eps = Self::thermal_strain(thermal, a.shape.values()); + let s: [f64; 4] = std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * eps[c]).sum()); + // Bᵀ σ for one strain column g: u takes εxx, εzz and γxy; v takes εyy and γxy. + let add = |out: &mut [f64], k: usize, g: &Column| { + if let Some(x) = out.get_mut(2 * k) { + *x += a.dv * (g[0] * s[0] + g[2] * s[2] + g[1] * s[3]); + } + if let Some(x) = out.get_mut(2 * k + 1) { + *x += a.dv * (g[1] * s[1] + g[0] * s[3]); + } + }; + for (k, g) in a.g.iter().take(nodes).enumerate() { + add(&mut fu, k, g); + } + if let Some(cols) = modes.and_then(|m| m.get(i)) { + for (k, g) in cols.iter().enumerate() { + add(&mut fa, k, g); + } + } + } + (fu, fa) + } + + /// The equivalent nodal forces (u, v) of a thermal strain, ∫ Bᵀ D ε_th dV with the + /// theory's elasticity (plane stress leaves εzz free, plane strain holds it, axisymmetric + /// strains the hoop direction too); incompatible modes are condensed as in the stiffness. + pub fn thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> Result> { + thermal.check(self.ty.nodes)?; + let d2 = self.theory.elasticity(d)?; + let n = 2 * self.ty.nodes; + let mut fu = if self.formulation == Formulation::IncompatibleModes { + let (full, points, modes) = self.full_stiffness(&d2)?; + let (mut fu, mut fa) = Self::thermal_forces(&d2, thermal, &points, modes.as_ref(), self.ty.nodes); + Cholesky::new(&sub_block(&full, n, n, INCOMPATIBLE_DOFS))?.solve(&mut fa); + for (i, x) in fu.iter_mut().enumerate() { + *x -= (0..INCOMPATIBLE_DOFS).map(|a| full.get(i, n + a) * fa.get(a).copied().unwrap_or(0.0)).sum::(); + } + fu + } else { + Self::thermal_forces(&d2, thermal, &self.points()?, None, self.ty.nodes).0 + }; + fu.resize(n, 0.0); + Ok(fu.as_chunks::<2>().0.to_vec()) + } + /// Strain and stress at every integration point for the nodal displacements `u` (u, v). pub fn stresses(&self, d: &Elasticity, u: &[[f64; 2]]) -> Result> { + self.stresses_with(d, u, None) + } + + /// Strain and stress at every integration point for the nodal displacements `u` (u, v) and, + /// if given, a thermal strain: the stress is D (ε − ε_th), the strain reported the total. + pub fn stresses_with(&self, d: &Elasticity, u: &[[f64; 2]], thermal: Option<&ThermalStrain>) -> Result> { if u.len() != self.ty.nodes { return Err(Error::Invalid(format!("{} displacements for a {}-node element", u.len(), self.ty.nodes))); } + if let Some(t) = thermal { + t.check(self.ty.nodes)?; + } let d2 = self.theory.elasticity(d)?; let (full, points, modes) = if self.formulation == Formulation::IncompatibleModes { self.full_stiffness(&d2)? } else { (Matrix::zeros(0), self.points()?, None) }; - // Incompatible modes: α = −K_aa⁻¹ K_au u. + // Incompatible modes: K_aa α = f_a − K_au u (f_a the thermal load on the modes). let alpha = if modes.is_some() { let n = 2 * self.ty.nodes; let flat: Vec = u.iter().flatten().copied().collect(); let mut rhs: Vec = (0..INCOMPATIBLE_DOFS).map(|a| -flat.iter().enumerate().map(|(i, x)| full.get(n + a, i) * x).sum::()).collect(); + if let Some(t) = thermal { + let (_, fa) = Self::thermal_forces(&d2, t, &points, modes.as_ref(), self.ty.nodes); + for (r, f) in rhs.iter_mut().zip(fa) { + *r += f; + } + } Cholesky::new(&sub_block(&full, n, n, INCOMPATIBLE_DOFS))?.solve(&mut rhs); rhs } else { @@ -507,11 +585,19 @@ impl<'a> Plane<'a> { add(g, [alpha.get(2 * k).copied().unwrap_or(0.0), alpha.get(2 * k + 1).copied().unwrap_or(0.0)]); } } + let eth = thermal.map(|t| Self::thermal_strain(t, a.shape.values())); if let Theory::PlaneStress { .. } = self.theory { - e[2] = -(d[2][0] * e[0] + d[2][1] * e[1] + d[2][3] * e[3]) / d[2][2]; + e[2] = match eth { + // σzz = 0 with the thermal strain: εzz = εth,zz − D₃·(ε − εth)/D₃₃ in plane. + Some(t) => t[2] - (d[2][0] * (e[0] - t[0]) + d[2][1] * (e[1] - t[1]) + d[2][3] * e[3]) / d[2][2], + None => -(d[2][0] * e[0] + d[2][1] * e[1] + d[2][3] * e[3]) / d[2][2], + }; } // Plane stress: D̃ has no zz row or column, so the εzz just recovered adds nothing. - let s: [f64; 4] = std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * e[c]).sum()); + let s: [f64; 4] = match eth { + Some(t) => std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * (e[c] - t[c])).sum()), + None => std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * e[c]).sum()), + }; out.push(PointStress { xyz: a.x, strain: [e[0], e[1], e[2], e[3], 0.0, 0.0], stress: [s[0], s[1], s[2], s[3], 0.0, 0.0] }); } Ok(out) diff --git a/crates/fem/src/tests.rs b/crates/fem/src/tests.rs index f54c260..ac7cf9a 100644 --- a/crates/fem/src/tests.rs +++ b/crates/fem/src/tests.rs @@ -208,3 +208,33 @@ fn section_formulations() { assert!(Formulation::from_section(ElementKind::Tet4, &hybrid).is_err()); assert!(Formulation::from_section(ElementKind::Quad4, &full).is_err()); } + +#[test] +fn uniform_thermal_strain_is_free_expansion_and_its_load_balances_it() { + // A uniform ΔT with orthotropic α: the free expansion u = diag(α ΔT) x is stress-free, and + // the thermal load equals K u for it (the element is in equilibrium with no other load). + let d = orthotropic([1000.0, 800.0, 600.0], [0.3, 0.25, 0.2], [400.0, 300.0, 350.0]).unwrap(); + let alpha = [1.2e-5, 2.0e-5, 0.7e-5]; + for (kind, f) in CASES { + let xyz = distorted(kind, true); + let solid = Solid::new(kind, &xyz, f).unwrap(); + let th = crate::ThermalStrain { alpha, delta: vec![80.0; xyz.len()] }; + let u: Vec<[f64; 3]> = xyz.iter().map(|x| [0, 1, 2].map(|k| alpha[k] * 80.0 * x[k])).collect(); + for p in solid.stresses_with(&d, &u, Some(&th)).unwrap() { + assert!(p.stress.iter().all(|s| s.abs() < 1e-9), "{kind:?} {f:?}: {:?}", p.stress); + for k in 0..3 { + assert!((p.strain[k] - alpha[k] * 80.0).abs() < 1e-12); + } + } + let load = solid.thermal_load(&d, &th).unwrap(); + let ku = solid.stiffness(&d).unwrap().mul_vec(&u.iter().flatten().copied().collect::>()); + let scale = ku.iter().fold(0.0f64, |m, v| m.max(v.abs())); + for (a, b) in load.iter().flatten().zip(&ku) { + assert!((a - b).abs() < 1e-9 * scale, "{kind:?} {f:?}: {a} vs {b}"); + } + // Without a thermal strain the two stress paths agree bit for bit. + assert_eq!(solid.stresses_with(&d, &u, None).unwrap(), solid.stresses(&d, &u).unwrap()); + assert!(solid.thermal_load(&d, &crate::ThermalStrain { alpha, delta: vec![1.0; 2] }).is_err()); + assert!(solid.thermal_load(&d, &crate::ThermalStrain { alpha: [f64::NAN; 3], delta: vec![1.0; xyz.len()] }).is_err()); + } +} diff --git a/crates/fem/src/truss.rs b/crates/fem/src/truss.rs index 0c11b60..647d4b6 100644 --- a/crates/fem/src/truss.rs +++ b/crates/fem/src/truss.rs @@ -29,7 +29,7 @@ use modelcraft_model::{Elastic, ElementKind}; use crate::basis::Basis; use crate::dense::Matrix; use crate::dof::{DofSet, ElementDofs}; -use crate::element::PointStress; +use crate::element::{PointStress, ThermalStrain}; use crate::quadrature::gauss_legendre; use crate::{Error, Result}; @@ -196,12 +196,51 @@ impl<'a> Truss<'a> { /// Strain and stress at the integration points as global tensors (see the module docs). pub fn stresses(&self, young: f64, u: &[[f64; 3]]) -> Result> { - let strains = self.axial_strains(u)?; + self.stresses_with(young, u, None) + } + + /// The axial thermal strain at each integration point: the thermal strain tensor + /// diag(α) ΔT along the tangent, Σ αₖ tₖ² ΔT, with ΔT interpolated from the nodes. + fn thermal_strains(&self, thermal: &ThermalStrain) -> Result> { + thermal.check(self.nodes())?; Ok(self .points()? .iter() - .zip(strains) - .map(|(p, e)| PointStress { xyz: p.xyz, strain: uniaxial(e, p.tangent, 2.0), stress: uniaxial(young * e, p.tangent, 1.0) }) + .map(|p| { + let dt: f64 = p.n.iter().zip(&thermal.delta).map(|(a, b)| a * b).sum(); + (0..3).map(|k| thermal.alpha[k] * p.tangent[k] * p.tangent[k]).sum::() * dt + }) + .collect()) + } + + /// The equivalent nodal forces of a thermal strain, ∫ Bᵀ E A ε_th ds. + pub fn thermal_load(&self, young: f64, area: f64, thermal: &ThermalStrain) -> Result> { + let strains = self.thermal_strains(thermal)?; + let mut f = vec![[0.0; 3]; self.nodes()]; + for (p, e) in self.points()?.iter().zip(strains) { + let force = young * area * e * p.weight * p.jacobian; + for (fi, ds) in f.iter_mut().zip(p.ds) { + for k in 0..3 { + fi[k] += force * ds * p.tangent[k]; + } + } + } + Ok(f) + } + + /// Strain and stress at the integration points with an optional thermal strain: the + /// stress is E (ε − ε_th), the strain reported the total. + pub fn stresses_with(&self, young: f64, u: &[[f64; 3]], thermal: Option<&ThermalStrain>) -> Result> { + let strains = self.axial_strains(u)?; + let thermal = match thermal { + Some(t) => self.thermal_strains(t)?, + None => vec![0.0; strains.len()], + }; + Ok(self + .points()? + .iter() + .zip(strains.into_iter().zip(thermal)) + .map(|(p, (e, th))| PointStress { xyz: p.xyz, strain: uniaxial(e, p.tangent, 2.0), stress: uniaxial(young * (e - th), p.tangent, 1.0) }) .collect()) } @@ -401,4 +440,26 @@ mod tests { assert!(err < 0.01 * span.powi(4), "{span}: {err}"); } } + + #[test] + fn thermal_strain_follows_the_axis() { + // A skew rod3 with orthotropic α: the free expansion along the member, u = ε_t (x − x₀) + // with ε_t = Σ αₖ tₖ² ΔT, is stress-free and balances the thermal load (K u = f_th). + let t = [2.0f64, 1.0, 2.0].map(|c| c / 3.0); + let xyz = [[0.0; 3], t.map(|c| 3.0 * c), t.map(|c| 1.5 * c)]; + let truss = Truss::new(ElementKind::Rod3, &xyz).unwrap(); + let alpha = [1.0e-5, 2.0e-5, 3.0e-5]; + let th = ThermalStrain { alpha, delta: vec![40.0; 3] }; + let eps: f64 = (0..3).map(|k| alpha[k] * t[k] * t[k]).sum::() * 40.0; + let u: Vec<[f64; 3]> = xyz.iter().map(|x| x.map(|c| eps * c)).collect(); + for p in truss.stresses_with(2e5, &u, Some(&th)).unwrap() { + assert!(p.stress.iter().all(|s| s.abs() < 1e-9), "{:?}", p.stress); + } + let f = truss.thermal_load(2e5, 2.0, &th).unwrap(); + let ku = truss.stiffness(2e5, 2.0).unwrap().mul_vec(&u.iter().flatten().copied().collect::>()); + for (a, b) in f.iter().flatten().zip(&ku) { + assert!((a - b).abs() < 1e-9 * 2e5 * 2.0 * eps, "{a} vs {b}"); + } + assert!(truss.thermal_load(2e5, 2.0, &ThermalStrain { alpha, delta: vec![1.0] }).is_err()); + } } diff --git a/crates/model/src/bcs.rs b/crates/model/src/bcs.rs index 88920ca..e4babb4 100644 --- a/crates/model/src/bcs.rs +++ b/crates/model/src/bcs.rs @@ -142,20 +142,28 @@ pub enum BoundaryCondition { Displacement { target: Target, values: [Option; 6] }, /// Mirror symmetry about a plane through the nodes. Symmetry { target: Target, plane: Symmetry }, + /// A prescribed temperature (degree of freedom [`TEMPERATURE_DOF`]) at the target's nodes, + /// in heat transfer and coupled temperature–displacement steps. + Temperature { target: Target, value: f64 }, } +/// The model's number of the temperature degree of freedom (the deck formats' convention). +pub const TEMPERATURE_DOF: u8 = 11; + impl BoundaryCondition { pub fn target(&self) -> &Target { match self { - BoundaryCondition::Fixed { target } | BoundaryCondition::Displacement { target, .. } | BoundaryCondition::Symmetry { target, .. } => { - target - } + BoundaryCondition::Fixed { target } + | BoundaryCondition::Displacement { target, .. } + | BoundaryCondition::Symmetry { target, .. } + | BoundaryCondition::Temperature { target, .. } => target, } } - /// The constrained degrees of freedom. + /// The constrained structural degrees of freedom (1–6); none for a temperature. pub fn dofs(&self) -> Dofs { match self { + BoundaryCondition::Temperature { .. } => Dofs::none(), BoundaryCondition::Fixed { .. } => Dofs::all(), BoundaryCondition::Displacement { values, .. } => { Dofs(values.iter().enumerate().filter(|(_, v)| v.is_some()).fold(0, |m, (i, _)| m | (1 << i))) @@ -164,15 +172,28 @@ impl BoundaryCondition { } } - /// The prescribed value of a constrained degree of freedom (1–6); `None` when it is free. + /// Whether it prescribes temperatures (rather than motion). + pub fn is_thermal(&self) -> bool { + matches!(self, BoundaryCondition::Temperature { .. }) + } + + /// The prescribed value of a constrained degree of freedom (1–6, or [`TEMPERATURE_DOF`]); + /// `None` when it is free. pub fn value(&self, dof: u8) -> Option { match self { + BoundaryCondition::Temperature { value, .. } => (dof == TEMPERATURE_DOF).then_some(*value), BoundaryCondition::Displacement { values, .. } => values.get(usize::from(dof).checked_sub(1)?).copied().flatten(), _ => self.dofs().contains(dof).then_some(0.0), } } pub fn validate(&self) -> Result<()> { + if let BoundaryCondition::Temperature { target, value } = self { + finite("temperature", *value)?; + if let Target::ReferencePoint(r) = target { + return Err(Error::Invalid(format!("reference points have no temperature (`{r}`); prescribe it on a set"))); + } + } if let BoundaryCondition::Displacement { values, .. } = self { for (i, v) in values.iter().enumerate() { if let Some(v) = v { @@ -241,6 +262,11 @@ mod tests { assert_eq!(BoundaryCondition::Fixed { target: t.clone() }.value(6), Some(0.0)); assert!(BoundaryCondition::Displacement { target: t.clone(), values: [None; 6] }.validate().is_err()); assert!(BoundaryCondition::Displacement { target: t, values: [Some(f64::NAN), None, None, None, None, None] }.validate().is_err()); + let hot = BoundaryCondition::Temperature { target: Target::Set("s".into()), value: 80.0 }; + assert_eq!((hot.value(TEMPERATURE_DOF), hot.value(1), hot.dofs().list()), (Some(80.0), None, vec![])); + assert!(hot.is_thermal()); + assert!(BoundaryCondition::Temperature { target: Target::Set("s".into()), value: f64::INFINITY }.validate().is_err()); + assert!(BoundaryCondition::Temperature { target: Target::ReferencePoint("RP".into()), value: 1.0 }.validate().is_err()); let json = serde_json::json!({"type": "symmetry", "target": {"reference_point": "RP"}, "plane": "z"}); let parsed: BoundaryCondition = serde_json::from_value(json.clone()).unwrap(); assert_eq!(serde_json::to_value(&parsed).unwrap(), json); diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index 76b881b..82d9fa1 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -169,6 +169,8 @@ impl Model { let mut d = Diagnostics::default(); let mut resolved = Resolved::default(); let units = self.units; + // Whether a step needs elasticity (every step but pure heat transfer). + let structural = self.steps.is_empty() || self.steps.iter().any(|s| s.procedure.is_structural()); // A model of fluid steps alone needs no mesh: its domains carry their own grid; one of // point masses, springs and dashpots alone is its own mesh. @@ -289,7 +291,7 @@ impl Model { for (other, labels) in taken { d.error(&subject, format!("elements {} are also in section `{other}`", first_labels(&labels))); } - if self.materials.get(&sec.material).is_some_and(|m| m.elastic.is_none()) { + if structural && self.materials.get(&sec.material).is_some_and(|m| m.elastic.is_none()) { d.error(&subject, format!("material `{}` has no elastic constants", sec.material)); } } @@ -307,7 +309,13 @@ impl Model { let why = match s.procedure { Procedure::Frequency { .. } => "extracts frequencies", Procedure::ExplicitDynamic(_) => "is dynamic", - Procedure::Static {} | Procedure::Fluid { .. } => return None, + Procedure::HeatTransfer { transient: Some(_) } => "is transient heat transfer", + Procedure::Static {} + | Procedure::Fluid { .. } + | Procedure::HeatTransfer { transient: None } + | Procedure::CoupledTemperatureDisplacement {} => { + return None; + } }; Some(format!("step `{}` {why}", s.name)) }) @@ -327,6 +335,33 @@ impl Model { } } + // Heat transfer needs conductivities (and heat capacities when transient); thermal strain + // needs expansion coefficients. + let used: BTreeSet<&str> = self.sections.values().map(|s| s.material.as_str()).collect(); + let thermal = |m: &str| self.materials.get(m).and_then(|m| m.thermal); + for s in self.steps.iter().filter(|s| s.procedure.is_thermal()) { + for &m in &used { + let t = thermal(m); + if t.and_then(|t| t.conductivity).is_none() { + d.error(format!("material `{m}`"), format!("has no conductivity, and step `{}` conducts heat", s.name)); + } + if matches!(s.procedure, Procedure::HeatTransfer { transient: Some(_) }) && t.and_then(|t| t.specific_heat).is_none() { + d.error(format!("material `{m}`"), format!("has no specific heat, and step `{}` is transient heat transfer", s.name)); + } + } + } + let strains_thermally = self.steps.iter().any(|s| { + matches!(s.procedure, Procedure::CoupledTemperatureDisplacement {}) + || s.loads.iter().any(|l| matches!(self.loads.get(l), Some(Load::Temperature { .. }))) + }); + if strains_thermally { + for &m in &used { + if thermal(m).and_then(|t| t.expansion).is_none() { + d.warning(format!("material `{m}`"), "has no expansion coefficient, so temperatures do not strain it"); + } + } + } + // Reference points and couplings. resolved.reference_point_nodes = self.reference_point_nodes(); for (name, label) in &resolved.reference_point_nodes { @@ -447,9 +482,24 @@ impl Model { } // The fluid domain holds its own conditions; nothing structural applies. Procedure::Fluid { .. } => continue, + Procedure::HeatTransfer { .. } | Procedure::CoupledTemperatureDisplacement {} => { + let bcs = s.boundary_conditions.iter().filter_map(|b| self.boundary_conditions.get(b)); + let loads: Vec<&Load> = s.loads.iter().filter_map(|l| self.loads.get(l)).collect(); + let held = bcs.clone().any(BoundaryCondition::is_thermal) || loads.iter().any(|l| matches!(l, Load::Convection { .. })); + if !held && matches!(s.procedure, Procedure::HeatTransfer { transient: None } | Procedure::CoupledTemperatureDisplacement {}) { + d.warning(&subject, "no temperature boundary condition or convection is active, so the temperature level is undetermined"); + } + } } - // A free body is a natural case of transient dynamics (flight, impact). - if s.boundary_conditions.is_empty() && !s.procedure.is_dynamic() { + for l in s.loads.iter().filter(|l| self.loads.get(*l).is_some_and(|l| l.is_thermal() && !s.procedure.is_thermal())) { + d.warning(&subject, format!("heat load `{l}` has no effect in a {} step", s.procedure.name())); + } + // A free body is a natural case of transient dynamics (flight, impact); heat transfer + // has no rigid motion. + if s.procedure.is_structural() + && s.boundary_conditions.iter().filter_map(|b| self.boundary_conditions.get(b)).all(BoundaryCondition::is_thermal) + && !s.procedure.is_dynamic() + { d.warning(&subject, "no boundary condition is active, so the model is free to move as a rigid body"); } } @@ -708,6 +758,8 @@ impl Model { Load::Force { vector, .. } | Load::Traction { vector, .. } | Load::BodyForce { vector, .. } => vector[2] == 0.0, Load::Gravity { acceleration, .. } => acceleration[2] == 0.0, Load::Moment { .. } | Load::Pressure { .. } => true, + // Scalars. + Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. } | Load::Temperature { .. } => true, Load::Centrifugal { point, axis, .. } => { let ok = if extent == Some(Extent::Revolved) { axis[0] == 0.0 && axis[2] == 0.0 && point[0] == 0.0 && point[2] == 0.0 diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index 5dec6cc..85fb301 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -50,18 +50,20 @@ use std::collections::BTreeMap; use serde::{Deserialize, Serialize}; -pub use bcs::{BoundaryCondition, Dofs, Symmetry, Target}; +pub use bcs::{BoundaryCondition, Dofs, Symmetry, TEMPERATURE_DOF, Target}; pub use check::{Diagnostic, Resolved, ResolvedSet, Severity}; pub use constraints::{Coupling, CouplingKind, ReferencePoint}; pub use discrete::{Direction, Mass, PointMass, ResolvedDirection, ResolvedMass, ResolvedSpring, Spring}; pub use fluids::{BoxFace, BoxFaces, FaceCondition, FluidBoundary, FluidDomain, Profile, SurfaceBoundary, SurfaceCondition}; pub use geometry::{Classification, Entity, GeometryRef}; pub use loads::{Distribution, Load, centrifugal_acceleration}; -pub use materials::{Damage, Elastic, Erosion, Fluid, Hardening, HardeningPoint, Material, Plastic}; +pub use materials::{Damage, Directional, Elastic, Erosion, Fluid, Hardening, HardeningPoint, Material, Plastic, Thermal}; pub use mesh::{Element, ElementEdge, ElementFace, ElementKind, Node, check_coordinates}; pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, SectionKind}; pub use sets::Set; -pub use steps::{BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, OutputVariable, Procedure, RigidWall, Step, StressRate}; +pub use steps::{ + BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, OutputVariable, Procedure, RigidWall, Step, StressRate, Transient, +}; pub use units::UnitSystem; /// The largest node or element label. Nastran's GRID and element ids stop at 99 999 999 (eight diff --git a/crates/model/src/loads.rs b/crates/model/src/loads.rs index d14da26..24dfedd 100644 --- a/crates/model/src/loads.rs +++ b/crates/model/src/loads.rs @@ -10,7 +10,7 @@ use std::collections::BTreeMap; use serde::{Deserialize, Serialize}; use crate::bcs::Target; -use crate::{Error, Model, Result, context, finite, insert_new, not_in_use, remove_existing, replace_existing, require}; +use crate::{Error, Model, Result, context, finite, insert_new, not_in_use, positive, remove_existing, replace_existing, require}; /// How a force or moment on a set is shared among the set's nodes. #[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] @@ -70,6 +70,17 @@ pub enum Load { #[serde(default, skip_serializing_if = "Option::is_none")] set: Option, }, + /// Heat flux per area into the set's element faces (positive heats the body). + SurfaceFlux { set: String, flux: f64 }, + /// Convection from the set's element faces to a surrounding medium at `sink` temperature: + /// heat flux h (T_sink − T) into the body, with film coefficient `film` (h). + Convection { set: String, film: f64, sink: f64 }, + /// Heat generated per volume in the set's elements. + HeatSource { set: String, power: f64 }, + /// A temperature field at the target's nodes in a stress step: it strains the material + /// thermally (ε = α (T − T_ref)); nodes it does not reach stay at their material's reference + /// temperature. + Temperature { target: Target, value: f64 }, } /// The centrifugal acceleration ω² r at `x` of a rotation at `angular_velocity` about the axis @@ -94,14 +105,29 @@ impl Load { Load::Gravity { .. } => "gravity", Load::BodyForce { .. } => "body_force", Load::Centrifugal { .. } => "centrifugal", + Load::SurfaceFlux { .. } => "surface_flux", + Load::Convection { .. } => "convection", + Load::HeatSource { .. } => "heat_source", + Load::Temperature { .. } => "temperature", } } + /// Whether it is a heat transfer load (flux, convection, heat source), active in heat + /// transfer and coupled temperature–displacement steps. + pub fn is_thermal(&self) -> bool { + matches!(self, Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. }) + } + /// The set the load acts on, if any. pub fn set(&self) -> Option<&str> { match self { - Load::Force { target, .. } | Load::Moment { target, .. } => target.set(), - Load::Pressure { set, .. } | Load::BodyForce { set, .. } | Load::Traction { set, .. } => Some(set), + Load::Force { target, .. } | Load::Moment { target, .. } | Load::Temperature { target, .. } => target.set(), + Load::Pressure { set, .. } + | Load::BodyForce { set, .. } + | Load::Traction { set, .. } + | Load::SurfaceFlux { set, .. } + | Load::Convection { set, .. } + | Load::HeatSource { set, .. } => Some(set), Load::Gravity { set, .. } | Load::Centrifugal { set, .. } => set.as_deref(), } } @@ -109,7 +135,7 @@ impl Load { /// The reference point the load acts on, if any. pub fn reference_point(&self) -> Option<&str> { match self { - Load::Force { target, .. } | Load::Moment { target, .. } => target.reference_point(), + Load::Force { target, .. } | Load::Moment { target, .. } | Load::Temperature { target, .. } => target.reference_point(), _ => None, } } @@ -137,6 +163,24 @@ impl Load { } return Ok(()); } + Load::SurfaceFlux { flux: v, .. } | Load::HeatSource { power: v, .. } => { + if finite(self.kind(), *v)? == 0.0 { + return Err(Error::Invalid(format!("a {} must not be zero", self.kind()))); + } + return Ok(()); + } + Load::Convection { film, sink, .. } => { + positive("film coefficient", *film)?; + finite("sink temperature", *sink)?; + return Ok(()); + } + Load::Temperature { target, value } => { + finite("temperature", *value)?; + if let Target::ReferencePoint(r) = target { + return Err(Error::Invalid(format!("reference points have no temperature (`{r}`); apply it to a set"))); + } + return Ok(()); + } }; for c in vector { finite(self.kind(), *c)?; @@ -185,3 +229,23 @@ impl Model { remove_existing(&mut self.loads, "load", name).map(|_| ()) } } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn heat_loads_validate() { + let set = || "s".to_string(); + assert!(Load::SurfaceFlux { set: set(), flux: 2.0 }.validate().is_ok()); + assert!(Load::SurfaceFlux { set: set(), flux: 0.0 }.validate().is_err()); + assert!(Load::HeatSource { set: set(), power: f64::NAN }.validate().is_err()); + assert!(Load::Convection { set: set(), film: 0.0, sink: 1.0 }.validate().is_err()); + assert!(Load::Convection { set: set(), film: 1.0, sink: f64::INFINITY }.validate().is_err()); + assert!(Load::Temperature { target: Target::ReferencePoint("RP".into()), value: 1.0 }.validate().is_err()); + assert!(Load::Temperature { target: Target::Set(set()), value: -40.0 }.validate().is_ok()); + let l = Load::Convection { set: set(), film: 10.0, sink: 20.0 }; + assert!(l.is_thermal() && l.set() == Some("s") && l.kind() == "convection"); + assert!(!Load::Temperature { target: Target::Set(set()), value: 1.0 }.is_thermal()); + } +} diff --git a/crates/model/src/materials.rs b/crates/model/src/materials.rs index c3d7ef6..0571bf1 100644 --- a/crates/model/src/materials.rs +++ b/crates/model/src/materials.rs @@ -1,4 +1,5 @@ -//! Materials: linear elasticity, mass density, fluid behaviour, and the inelastic behaviours of +//! Materials: linear elasticity, mass density, fluid behaviour, thermal properties (conductivity, +//! specific heat, expansion), and the inelastic behaviours of //! explicit dynamics: von Mises (J2) plasticity with a bilinear, tabulated or Johnson–Cook flow stress, //! Johnson–Cook damage, and element erosion criteria. //! @@ -237,6 +238,67 @@ impl Plastic { } } +/// A property given per material axis: one value for every direction, or three along the +/// material axes 1, 2, 3 (the global X, Y, Z until material orientations land). +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(untagged)] +pub enum Directional { + Isotropic(f64), + Orthotropic([f64; 3]), +} + +impl Directional { + /// The values along the three axes. + pub fn axes(self) -> [f64; 3] { + match self { + Directional::Isotropic(v) => [v; 3], + Directional::Orthotropic(v) => v, + } + } +} + +/// Heat conduction and thermal expansion. +#[derive(Debug, Clone, Copy, Default, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Thermal { + /// Conductivity k (power per length per temperature): q = −k ∇T. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub conductivity: Option, + /// Specific heat c (energy per mass per temperature); with the density, the heat capacity + /// of transient heat transfer. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub specific_heat: Option, + /// Linear expansion coefficient α (strain per temperature): ε_th = α (T − T_ref). + #[serde(default, skip_serializing_if = "Option::is_none")] + pub expansion: Option, + /// The temperature at which thermal strain is zero. + #[serde(default)] + pub reference_temperature: f64, +} + +impl Thermal { + pub fn validate(&self) -> Result<()> { + if let Some(k) = self.conductivity { + for v in k.axes() { + positive("conductivity", v)?; + } + } + if let Some(c) = self.specific_heat { + positive("specific heat", c)?; + } + if let Some(a) = self.expansion { + for v in a.axes() { + finite("expansion coefficient", v)?; + } + } + finite("reference temperature", self.reference_temperature)?; + if self.conductivity.is_none() && self.specific_heat.is_none() && self.expansion.is_none() { + return Err(Error::Invalid("thermal properties need a conductivity, a specific heat or an expansion coefficient".into())); + } + Ok(()) + } +} + /// Damage that accumulates with plastic strain; a material point fails when it reaches 1. #[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] #[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] @@ -270,7 +332,7 @@ pub struct Erosion { pub struct Material { #[serde(default, skip_serializing_if = "Option::is_none")] pub elastic: Option, - /// Mass per volume, needed by frequency steps and gravity. + /// Mass per volume, needed by frequency steps, gravity and transient heat transfer. #[serde(default, skip_serializing_if = "Option::is_none")] pub density: Option, /// How it flows, for fluid domains (needs a density). @@ -286,12 +348,14 @@ pub struct Material { /// Erosion criteria (explicit dynamics). #[serde(default, skip_serializing_if = "Option::is_none")] pub erosion: Option, + #[serde(default, skip_serializing_if = "Option::is_none")] + pub thermal: Option, } impl Material { pub fn validate(&self) -> Result<()> { - if self.elastic.is_none() && self.density.is_none() && self.fluid.is_none() { - return Err(Error::Invalid("a material needs elastic constants, a density or fluid behaviour".into())); + if self.elastic.is_none() && self.density.is_none() && self.fluid.is_none() && self.thermal.is_none() { + return Err(Error::Invalid("a material needs elastic constants, a density, fluid behaviour or thermal properties".into())); } if let Some(f) = &self.fluid { positive("viscosity", f.viscosity())?; @@ -332,6 +396,9 @@ impl Material { return Err(Error::Invalid("a plastic-strain erosion criterion needs plasticity".into())); } } + if let Some(t) = &self.thermal { + t.validate()?; + } Ok(()) } } @@ -476,4 +543,32 @@ mod tests { let text = serde_json::to_string(&Material { plastic: Some(jc()), damage: Some(damage), ..steel() }).unwrap(); assert_eq!(serde_json::from_str::(&text).unwrap(), Material { plastic: Some(jc()), damage: Some(damage), ..steel() }); } + + #[test] + fn thermal_properties() { + let m: Material = serde_json::from_value(serde_json::json!({ + "thermal": {"conductivity": [1, 2, 3], "specific_heat": 460, "expansion": 1.2e-5, "reference_temperature": 20} + })) + .unwrap(); + let t = m.thermal.unwrap(); + assert_eq!(t.conductivity.unwrap().axes(), [1.0, 2.0, 3.0]); + assert_eq!(t.expansion.unwrap().axes(), [1.2e-5; 3]); + assert!(m.validate().is_ok(), "a heat-only material needs no elasticity"); + let back: Material = serde_json::from_value(serde_json::to_value(&m).unwrap()).unwrap(); + assert_eq!(back, m); + for bad in [ + Thermal::default(), + Thermal { conductivity: Some(Directional::Isotropic(0.0)), ..Thermal::default() }, + Thermal { conductivity: Some(Directional::Orthotropic([1.0, -1.0, 1.0])), ..Thermal::default() }, + Thermal { specific_heat: Some(f64::NAN), ..Thermal::default() }, + Thermal { expansion: Some(Directional::Isotropic(f64::INFINITY)), ..Thermal::default() }, + Thermal { expansion: Some(Directional::Isotropic(1e-5)), reference_temperature: f64::NAN, ..Thermal::default() }, + ] { + assert!(bad.validate().is_err(), "{bad:?}"); + } + assert!(serde_json::from_value::(serde_json::json!({"thermal": {"conductivity": [1, 2]}})).is_err()); + assert!(serde_json::from_value::(serde_json::json!({"thermal": {"heat": 1}})).is_err()); + // Negative expansion (invar-like behaviour of some composites) is allowed. + assert!(Thermal { expansion: Some(Directional::Isotropic(-1e-6)), ..Thermal::default() }.validate().is_ok()); + } } diff --git a/crates/model/src/nodal.rs b/crates/model/src/nodal.rs index 399b50b..e3dfa1c 100644 --- a/crates/model/src/nodal.rs +++ b/crates/model/src/nodal.rs @@ -97,8 +97,9 @@ fn rule(corners: usize) -> Vec<(f64, f64, f64)> { /// Visit the quadrature points of a face (node coordinates in the face's node order): the shape /// function values and the weighted normal (∂x/∂ξ × ∂x/∂η) w, whose length is the point's share of -/// the area and whose direction is the face's right-hand normal. -fn integrate_face(xyz: &[[f64; 3]], mut visit: impl FnMut(&[f64], [f64; 3])) -> Result<()> { +/// the area and whose direction is the face's right-hand normal. Every surface load (pressure, +/// traction, heat flux, convection) integrates with it, so they share one rule. +pub fn integrate_face(xyz: &[[f64; 3]], mut visit: impl FnMut(&[f64], [f64; 3])) -> Result<()> { if xyz.iter().flatten().any(|c| !c.is_finite()) { return Err(Error::Invalid("face coordinates must be finite".into())); } diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index 72a260a..198d70e 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -76,6 +76,62 @@ pub enum Procedure { #[serde(default, skip_serializing_if = "Option::is_none")] magic: Option, }, + /// Heat conduction: the steady temperature field, or its history over time when `transient` + /// is given. Temperature boundary conditions and heat loads (flux, convection, sources) + /// act; mechanical ones are ignored. + HeatTransfer { + #[serde(default, skip_serializing_if = "Option::is_none")] + transient: Option, + }, + /// Steady heat conduction and the static equilibrium of the thermal strains it causes, + /// together: every boundary condition and load of the step acts. + CoupledTemperatureDisplacement {}, +} + +/// Most time increments of a transient heat transfer step. +pub const MAX_HEAT_INCREMENTS: f64 = 100_000.0; + +/// Time integration of a transient heat transfer step: the generalized trapezoidal (θ) method +/// with fixed increments, from a uniform initial temperature. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Transient { + pub end_time: f64, + pub increment: f64, + /// 1 is backward Euler (the default), ½ Crank–Nicolson; θ ≥ ½ is unconditionally stable. + #[serde(default = "default_theta")] + pub theta: f64, + #[serde(default)] + pub initial_temperature: f64, +} + +fn default_theta() -> f64 { + 1.0 +} + +impl Transient { + /// Backward Euler from 0. + pub fn new(end_time: f64, increment: f64) -> Transient { + Transient { end_time, increment, theta: default_theta(), initial_temperature: 0.0 } + } + + /// The number of increments (the last one may be shorter). + pub fn increments(&self) -> usize { + (self.end_time / self.increment - 1e-9).ceil().max(1.0) as usize + } + + pub fn validate(&self) -> Result<()> { + let end = positive("end time", self.end_time)?; + let dt = positive("time increment", self.increment)?; + if end / dt > MAX_HEAT_INCREMENTS { + return Err(Error::Invalid(format!("a transient heat step takes at most {MAX_HEAT_INCREMENTS} increments, got {}", end / dt))); + } + if !(0.5..=1.0).contains(&self.theta) { + return Err(Error::Invalid(format!("θ is in [0.5, 1] (unconditionally stable), got {}", self.theta))); + } + finite("initial temperature", self.initial_temperature)?; + Ok(()) + } } /// The settings of an explicit dynamic step. @@ -336,6 +392,8 @@ impl Procedure { Procedure::Static {} => "static", Procedure::Frequency { .. } => "frequency", Procedure::ExplicitDynamic(_) => "explicit_dynamic", + Procedure::HeatTransfer { .. } => "heat_transfer", + Procedure::CoupledTemperatureDisplacement {} => "coupled_temperature_displacement", Procedure::Fluid { .. } => "fluid", } } @@ -344,10 +402,21 @@ impl Procedure { pub fn fluid_domain(&self) -> Option<&str> { match self { Procedure::Fluid { domain, .. } => Some(domain), - Procedure::Static {} | Procedure::Frequency { .. } | Procedure::ExplicitDynamic(_) => None, + _ => None, } } + /// Whether the procedure solves for the motion of the solids (not pure heat transfer or + /// fluid flow). + pub fn is_structural(&self) -> bool { + !matches!(self, Procedure::HeatTransfer { .. } | Procedure::Fluid { .. }) + } + + /// Whether the procedure solves for temperatures. + pub fn is_thermal(&self) -> bool { + matches!(self, Procedure::HeatTransfer { .. } | Procedure::CoupledTemperatureDisplacement {}) + } + /// Whether the procedure is transient dynamics (it needs mass; a free body is legitimate). pub fn is_dynamic(&self) -> bool { matches!(self, Procedure::ExplicitDynamic(_)) @@ -357,7 +426,7 @@ impl Procedure { pub fn sets(&self) -> impl Iterator { let explicit = match self { Procedure::ExplicitDynamic(e) => Some(e), - Procedure::Static {} | Procedure::Frequency { .. } | Procedure::Fluid { .. } => None, + _ => None, }; explicit.into_iter().flat_map(ExplicitDynamic::sets) } @@ -387,6 +456,9 @@ impl Procedure { return Err(Error::Invalid(format!("the two-relaxation-time parameter Λ must be above 0 and at most 1, got {m}"))); } } + if let Procedure::HeatTransfer { transient: Some(t) } = self { + return t.validate(); + } if let Procedure::Frequency { modes, min_hz, max_hz } = *self { if !(1..=MAX_MODES).contains(&modes) { return Err(Error::Invalid(format!("a frequency step extracts 1 to {MAX_MODES} modes, got {modes}"))); @@ -413,6 +485,10 @@ impl Procedure { Procedure::ExplicitDynamic(_) => [OutputVariable::Displacement, OutputVariable::Stress].into_iter().collect(), // Velocity, pressure and wall forces, always. Procedure::Fluid { .. } => BTreeSet::new(), + Procedure::HeatTransfer { .. } => [OutputVariable::Temperature, OutputVariable::HeatFlux].into_iter().collect(), + Procedure::CoupledTemperatureDisplacement {} => { + OutputVariable::ALL.into_iter().chain([OutputVariable::Temperature, OutputVariable::HeatFlux]).collect() + } } } } @@ -427,9 +503,14 @@ pub enum OutputVariable { ReactionForce, Stress, Strain, + /// Nodal temperatures (heat transfer). + Temperature, + /// Heat flux at the integration points (heat transfer). + HeatFlux, } impl OutputVariable { + /// The structural variables. pub const ALL: [OutputVariable; 4] = [OutputVariable::Displacement, OutputVariable::ReactionForce, OutputVariable::Stress, OutputVariable::Strain]; } @@ -671,4 +752,30 @@ mod tests { } assert!(serde_json::from_value::(serde_json::json!({"type": "explicit_dynamic", "end_time": 1, "dt": 2})).is_err()); } + + #[test] + fn heat_transfer_procedures() { + let steady: Procedure = serde_json::from_value(serde_json::json!({"type": "heat_transfer"})).unwrap(); + assert_eq!(steady, Procedure::HeatTransfer { transient: None }); + assert!(steady.is_thermal() && steady.validate().is_ok()); + assert_eq!(steady.default_output(), [OutputVariable::Temperature, OutputVariable::HeatFlux].into_iter().collect()); + let t: Procedure = + serde_json::from_value(serde_json::json!({"type": "heat_transfer", "transient": {"end_time": 10, "increment": 0.3}})).unwrap(); + let Procedure::HeatTransfer { transient: Some(tr) } = t else { panic!() }; + assert_eq!((tr.theta, tr.initial_temperature, tr.increments()), (1.0, 0.0, 34)); + assert_eq!(Transient::new(1.0, 0.25).increments(), 4); + for bad in [ + Transient { theta: 0.4, ..Transient::new(1.0, 0.1) }, + Transient { theta: f64::NAN, ..Transient::new(1.0, 0.1) }, + Transient::new(1.0, 0.0), + Transient::new(-1.0, 0.1), + Transient::new(1e300, 1e-300), + Transient { initial_temperature: f64::INFINITY, ..Transient::new(1.0, 0.1) }, + ] { + assert!(Procedure::HeatTransfer { transient: Some(bad) }.validate().is_err(), "{bad:?}"); + } + let coupled: Procedure = serde_json::from_value(serde_json::json!({"type": "coupled_temperature_displacement"})).unwrap(); + assert!(coupled.is_thermal() && coupled.default_output().contains(&OutputVariable::Stress)); + assert_eq!(coupled.name(), "coupled_temperature_displacement"); + } } diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 6359c71..d4408c9 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -14,7 +14,10 @@ use crate::constraints::{Budget, Reduction, coupling_rows}; use crate::eigen::{Operator, Wanted, lanczos}; use crate::ordering::nested_dissection; use crate::sparse::{Csr, NodeGraph, add_block, add_diagonal, add_scalar}; -use crate::system::{Mesh, Structural, StructuralKernel, inelastic_warnings, loads, prescribed, structural_body_forces, translation_block}; +use crate::system::{ + Mesh, Structural, StructuralKernel, add_thermal_loads, inelastic_warnings, loads, prescribed, structural_body_forces, temperature_loads, + translation_block, +}; use crate::{Error, Limits, Result}; /// Size and effort of a solution. @@ -99,6 +102,8 @@ pub struct StaticResult { pub reaction_total: [f64; 3], pub stats: Stats, pub warnings: Vec, + /// The temperatures of a coupled temperature–displacement step, per node. + pub temperature: Option>, } /// One natural mode. @@ -293,17 +298,43 @@ fn budget(limits: &Limits, mesh: &Mesh, graph: &NodeGraph, matrices: usize) -> B } /// Solve a linear static step: the step named `step`, or the model's only static step. +/// +/// A coupled temperature–displacement step is solved here too: its steady temperatures first +/// ([`crate::heat`]), then the equilibrium with their thermal strains. The coupling is one-way +/// (temperature strains the solid; deformation does not heat it), which is the whole of linear +/// quasi-static thermoelasticity without thermoelastic dissipation, so solving the two in turn +/// is the monolithic solution exactly. pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Result { let resolved = model.resolve()?; - let step = model.select_step(step, |p| matches!(p, Procedure::Static {}), "static")?; + let step = model.select_step( + step, + |p| matches!(p, Procedure::Static {} | Procedure::CoupledTemperatureDisplacement {}), + "static or coupled temperature–displacement", + )?; let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect(); warnings.extend(inelastic_warnings(model)); + // The temperature field: the coupled step's own heat solution, or the step's temperature + // loads. + let (field, temperature) = match step.procedure { + Procedure::CoupledTemperatureDisplacement {} => { + let (labels, t) = crate::heat::temperatures(model, &resolved, step, limits)?; + let at = |l: &u32| labels.binary_search(l).ok().and_then(|i| t.get(i)).copied(); + (Some(mesh.labels.iter().map(at).collect::>()), Some(mesh.labels.iter().map(|l| at(l).unwrap_or(0.0)).collect())) + } + _ => (temperature_loads(model, &resolved, &mesh, step)?, None), + }; + if field.is_some() && !mesh.structural.is_empty() { + warnings.push("shells and beams take no thermal strain yet; the temperatures strain the solid, 2D and truss elements only".into()); + } let (graph, order) = pattern(&mesh); ignored_dashpots(&mesh, &mut warnings); let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; let red = reduction(&setup, true, &order, budget(limits, &mesh, &graph, 1), &mut warnings)?; - let f = loads(model, &resolved, &mesh, step, &red.active)?; + let mut f = loads(model, &resolved, &mesh, step, &red.active)?; + if let Some(field) = &field { + add_thermal_loads(&mut f, &mesh, field)?; + } let k = assemble(&mesh, &graph, None)?.k; if red.neq() == 0 { return Err(Error::Invalid("every degree of freedom is prescribed; there is nothing to solve".into())); @@ -343,7 +374,8 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul let kernel = e.kernel()?; let at = |err| Error::element(e.label, err); let ue: Vec<[f64; 3]> = e.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); - let points = kernel.stresses(&e.d, &ue).map_err(at)?; + let thermal = field.as_deref().and_then(|t| e.thermal_strain(t)); + let points = kernel.stresses_with(&e.d, &ue, thermal.as_ref()).map_err(at)?; let values: Vec<[f64; 6]> = points.iter().map(|p| p.stress).collect(); let (ty, rule) = (kernel.ty(), kernel.rule()); let cached = extrapolations.iter().position(|(t, r, _)| std::ptr::eq(*t, ty) && *r == rule); @@ -369,7 +401,8 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul let truss = t.truss()?; let at = |err| Error::element(t.label, err); let ue: Vec<[f64; 3]> = t.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); - let points = truss.stresses(t.young, &ue).map_err(at)?; + let thermal = field.as_deref().and_then(|f| t.thermal_strain(f)); + let points = truss.stresses_with(t.young, &ue, thermal.as_ref()).map_err(at)?; let nodal = truss.extrapolate(&points.iter().map(|p| p.stress).collect::>()).map_err(at)?; for (n, s) in t.nodes.iter().zip(&nodal) { if let (Some(acc), Some(c)) = (sum.get_mut(*n), count.get_mut(*n)) { @@ -432,6 +465,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul reaction_total, stats: Stats { dofs: mesh.dofs(), equations: red.neq(), factor_entries: chol.entries(), lanczos_vectors: 0 }, warnings, + temperature, }) } @@ -585,6 +619,7 @@ pub enum Location { pub enum StepResult { Static(StaticResult), Frequency(FrequencyResult), + Heat(crate::heat::HeatResult), } impl StepResult { @@ -592,6 +627,7 @@ impl StepResult { match self { StepResult::Static(r) => &r.step, StepResult::Frequency(r) => &r.step, + StepResult::Heat(r) => &r.step, } } @@ -599,6 +635,7 @@ impl StepResult { match self { StepResult::Static(r) => &r.nodes, StepResult::Frequency(r) => &r.nodes, + StepResult::Heat(r) => &r.nodes, } } @@ -611,6 +648,7 @@ impl StepResult { match self { StepResult::Static(r) => r.reference_points.iter().position(|p| p.label == label), StepResult::Frequency(r) => r.reference_points.iter().position(|p| p.1 == label), + StepResult::Heat(_) => None, } .map(Location::Point) } diff --git a/crates/solver/src/heat.rs b/crates/solver/src/heat.rs new file mode 100644 index 0000000..a5c5e7c --- /dev/null +++ b/crates/solver/src/heat.rs @@ -0,0 +1,491 @@ +//! Heat transfer: steady and transient conduction with prescribed temperatures, surface fluxes, +//! convection and internal heat sources. +//! +//! Every node of the solid elements carries one temperature ([`modelcraft_fem::heat::HEAT_DOFS`]). +//! The conductivity K (with the convection films) and capacity C are assembled on the node +//! graph in element label order ([`modelcraft_fem::heat`]); prescribed temperatures are +//! eliminated, the free ones numbered by nested dissection and solved by sparse Cholesky. +//! +//! - **Steady**: K T = f. +//! - **Transient**, the generalized trapezoidal (θ) method with n equal increments Δt over the +//! step (T. J. R. Hughes, *The Finite Element Method*, Dover 2000, §8.2; unconditionally stable +//! for θ ≥ ½, θ = 1 backward Euler, θ = ½ Crank–Nicolson): +//! (C + θΔt K) Tₙ₊₁ = (C − (1 − θ)Δt K) Tₙ + Δt f, from a uniform initial temperature, with +//! the prescribed temperatures held from the first increment on. The matrix is factored once. +//! +//! The heat flux −k∇T is reported at the integration points and averaged at the nodes; the heat +//! flow at prescribed temperatures is the reaction K T − f (C Ṫ included in transient steps). + +use modelcraft_fem::element::{Extrapolation, extrapolation}; +use modelcraft_fem::heat::{Conduction, PointFlux, face_film, face_weights}; +use modelcraft_fem::{DofLayout, DofSet}; +use modelcraft_model::{ElementKind, Load, Model, Procedure, Resolved, SectionKind, Step, TEMPERATURE_DOF}; + +use crate::cholesky::{Factor, FactorError, UpperColumns}; +use crate::ordering::nested_dissection; +use crate::sparse::{Csr, NodeGraph, add_block}; +use crate::system::hold; +use crate::{Error, Limits, Result, Stats}; + +/// Most temperature values a transient result keeps (frames × nodes); frames are thinned +/// evenly to stay below it, the last increment always kept. +pub const MAX_TEMPERATURE_VALUES: usize = 50_000_000; + +/// A solid element conducting heat. +#[derive(Debug, Clone)] +pub struct HeatElem { + pub label: u32, + pub kind: ElementKind, + /// Indices into [`HeatMesh::labels`]. + pub nodes: Vec, + pub xyz: Vec<[f64; 3]>, + /// Conductivity along X, Y, Z. + pub conductivity: [f64; 3], + /// ρc, when the material has a density and a specific heat. + pub capacity: Option, +} + +impl HeatElem { + fn kernel(&self) -> Result> { + Conduction::new(self.kind, &self.xyz).map_err(|e| Error::element(self.label, e)) + } +} + +/// The conducting mesh: the nodes of the solid elements by label, and the elements by label. +#[derive(Debug, Clone)] +pub struct HeatMesh { + pub labels: Vec, + pub xyz: Vec<[f64; 3]>, + pub elements: Vec, + pub layout: DofLayout, +} + +impl HeatMesh { + pub fn new(model: &Model, resolved: &Resolved, max_dofs: usize) -> Result { + let mut elements = Vec::new(); + for (name, section) in model.sections() { + if !matches!(section.kind, SectionKind::Solid {}) { + return Err(Error::NotSupported(format!( + "heat transfer in {} sections (section `{name}`); it has solid elements so far", + section.kind.name() + ))); + } + let material = model + .material(§ion.material) + .ok_or_else(|| Error::Invalid(format!("section `{name}`: material `{}` is missing", section.material)))?; + let thermal = material.thermal.unwrap_or_default(); + let conductivity = + thermal.conductivity.ok_or_else(|| Error::Invalid(format!("material `{}` has no conductivity", section.material)))?.axes(); + let capacity = material.density.zip(thermal.specific_heat).map(|(rho, c)| rho * c); + let set = resolved.sets.get(§ion.set).ok_or_else(|| Error::Invalid(format!("set `{}` is not resolved", section.set)))?; + for &label in &set.elements { + let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; + elements.push((label, el.kind, el.nodes.clone(), conductivity, capacity)); + } + } + elements.sort_by_key(|e| e.0); + let layout = DofLayout::uniform(elements.iter().flat_map(|e| e.2.iter().copied()), DofSet::TEMPERATURE); + if layout.dofs() > max_dofs { + return Err(Error::Limit(format!("the model has {} temperatures; the limit is {max_dofs}", layout.dofs()))); + } + let labels: Vec = (0..layout.nodes()).filter_map(|s| layout.label(s)).collect(); + let xyz: Vec<[f64; 3]> = labels + .iter() + .map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("node {l} is missing")))) + .collect::>()?; + let mut mesh = HeatMesh { labels, xyz, elements: Vec::with_capacity(elements.len()), layout }; + for (label, kind, nodes, conductivity, capacity) in elements { + let idx: Vec = nodes.iter().filter_map(|n| mesh.layout.node(*n)).collect(); + let xyz = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); + mesh.elements.push(HeatElem { label, kind, nodes: idx, xyz, conductivity, capacity }); + } + Ok(mesh) + } + + fn node_of(&self, what: &str, label: u32) -> Result { + self.layout.node(label).ok_or_else(|| Error::Invalid(format!("{what}: node {label} belongs to no solid element"))) + } + + fn element(&self, label: u32) -> Option<&HeatElem> { + self.elements.binary_search_by_key(&label, |e| e.label).ok().and_then(|i| self.elements.get(i)) + } +} + +/// The heat flux in one element. +#[derive(Debug, Clone, PartialEq)] +pub struct ElementFlux { + pub label: u32, + pub kind: ElementKind, + pub points: Vec, +} + +/// The result of a heat transfer step. +#[derive(Debug, Clone, PartialEq)] +pub struct HeatResult { + pub step: String, + /// Node labels with a temperature, sorted; per-node vectors follow this order. + pub nodes: Vec, + /// The time of each frame (one frame at time 0 for a steady step). + pub times: Vec, + /// Nodal temperatures per frame; the last frame is the end of the step. + pub frames: Vec>, + /// Heat flux averaged at the nodes, at the end of the step. + pub flux: Vec<[f64; 3]>, + /// Heat flux at the integration points, at the end of the step. + pub elements: Vec, + /// The heat each prescribed temperature supplies to the body at the end of the step + /// (negative where it takes heat away; zero at free nodes). + pub reaction: Vec, + /// Heat supplied by fluxes, convection and sources (at the end temperatures), and the total + /// reaction heat flow; in a steady state they cancel. + pub applied_total: f64, + pub reaction_total: f64, + pub stats: Stats, + pub warnings: Vec, +} + +impl HeatResult { + /// The temperatures at the end of the step. + pub fn temperature(&self) -> &[f64] { + self.frames.last().map(Vec::as_slice).unwrap_or_default() + } +} + +/// The assembled heat system of a step: conductivity (with films), capacity when transient, +/// and the heat supplied. +struct HeatSystem { + k: Csr, + c: Option, + f: Vec, +} + +fn assemble(model: &Model, resolved: &Resolved, mesh: &HeatMesh, graph: &NodeGraph, step: &Step, transient: bool) -> Result { + let layout = &mesh.layout; + let mut k = graph.zero_matrix(layout); + let mut c = transient.then(|| graph.zero_matrix(layout)); + for e in &mesh.elements { + let kernel = e.kernel()?; + let at = |err| Error::element(e.label, err); + add_block(&mut k, layout, &e.nodes, DofSet::TEMPERATURE, &kernel.conductivity(e.conductivity).map_err(at)?.data)?; + if let Some(c) = c.as_mut() { + let rho_c = e.capacity.ok_or_else(|| Error::Invalid(format!("element {}: its material needs a density and a specific heat", e.label)))?; + add_block(c, layout, &e.nodes, DofSet::TEMPERATURE, &kernel.capacity(rho_c).map_err(at)?.data)?; + } + } + let mut f = vec![0.0; layout.dofs()]; + let add = |f: &mut [f64], nodes: &[usize], v: &[f64]| { + for (n, x) in nodes.iter().zip(v) { + if let Some(y) = f.get_mut(*n) { + *y += x; + } + } + }; + for name in &step.loads { + let load = model.load(name).ok_or_else(|| Error::Invalid(format!("load `{name}` is missing")))?; + let faces = || -> Result, Vec<[f64; 3]>)>> { + let set = load.set().unwrap_or_default(); + resolved + .sets + .get(set) + .map(|s| s.faces.as_slice()) + .unwrap_or_default() + .iter() + .map(|fc| { + let labels = model + .face_nodes(*fc) + .ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} face {} is missing", fc.element, fc.face)))?; + let idx: Vec = labels.iter().map(|l| mesh.node_of(&format!("load `{name}`"), *l)).collect::>()?; + let xyz = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); + Ok((fc.element, idx, xyz)) + }) + .collect() + }; + match load { + Load::SurfaceFlux { flux, .. } => { + for (label, idx, xyz) in faces()? { + let w = face_weights(&xyz).map_err(|e| Error::element(label, e))?; + add(&mut f, &idx, &w.iter().map(|w| w * flux).collect::>()); + } + } + Load::Convection { film, sink, .. } => { + for (label, idx, xyz) in faces()? { + let at = |e| Error::element(label, e); + add_block(&mut k, layout, &idx, DofSet::TEMPERATURE, &face_film(&xyz, *film).map_err(at)?.data)?; + let w = face_weights(&xyz).map_err(at)?; + add(&mut f, &idx, &w.iter().map(|w| w * film * sink).collect::>()); + } + } + Load::HeatSource { set, power } => { + for label in resolved.sets.get(set).map(|s| s.elements.as_slice()).unwrap_or_default() { + let e = mesh.element(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {label} is not a solid element")))?; + let q = e.kernel()?.source(*power).map_err(|err| Error::element(e.label, err))?; + add(&mut f, &e.nodes, &q); + } + } + // Mechanical loads and temperature fields do not act on heat conduction. + _ => {} + } + } + Ok(HeatSystem { k, c, f }) +} + +/// The prescribed temperature of every node (`None` = free) under the step's boundary +/// conditions. +fn prescribed(model: &Model, resolved: &Resolved, mesh: &HeatMesh, step: &Step) -> Result>> { + let mut out = vec![None; mesh.layout.dofs()]; + for name in &step.boundary_conditions { + let bc = model.boundary_condition(name).ok_or_else(|| Error::Invalid(format!("boundary condition `{name}` is missing")))?; + let Some(v) = bc.value(TEMPERATURE_DOF) else { continue }; + for &label in resolved.target_nodes(bc.target()) { + let Some(n) = mesh.layout.node(label) else { continue }; + hold(out.get_mut(n), v, |old| format!("node {label} is held at both {old} and {v} (boundary condition `{name}` and another)"))?; + } + } + Ok(out) +} + +/// The free temperatures numbered in nested-dissection order. +struct Equations { + /// Per node, its equation or `usize::MAX` when prescribed. + eq: Vec, + /// Per equation, its node. + node: Vec, +} + +impl Equations { + fn new(mesh: &HeatMesh, graph: &NodeGraph, fixed: &[Option]) -> Equations { + let mut eq = vec![usize::MAX; mesh.labels.len()]; + let mut node = Vec::new(); + for n in nested_dissection(&mesh.xyz, &graph.adjacency) { + if fixed.get(n).is_some_and(Option::is_none) + && let Some(e) = eq.get_mut(n) + { + *e = node.len(); + node.push(n); + } + } + Equations { eq, node } + } + + /// The free-free block of a node-level matrix, rows and columns in equation order. + fn reduce(&self, a: &Csr) -> Csr { + let mut ptr = vec![0]; + let (mut idx, mut val) = (Vec::new(), Vec::new()); + for &n in &self.node { + let (cols, vals) = a.row(n); + let mut row: Vec<(usize, f64)> = + cols.iter().zip(vals).filter_map(|(c, v)| self.eq.get(*c).filter(|e| **e != usize::MAX).map(|e| (*e, *v))).collect(); + row.sort_unstable_by_key(|t| t.0); + for (c, v) in row { + idx.push(c); + val.push(v); + } + ptr.push(idx.len()); + } + Csr { n: self.node.len(), ptr, idx, val } + } + + /// Solve for the free temperatures with right-hand side `rhs` (over every node; its free + /// rows are used) and write them into `t`, whose prescribed entries are kept. + fn solve(&self, chol: &Factor, rhs: &[f64], q: &mut Vec, t: &mut [f64]) -> Result<()> { + q.clear(); + q.extend(self.node.iter().map(|n| rhs.get(*n).copied().unwrap_or(0.0))); + chol.solve(q)?; + for (n, v) in self.node.iter().zip(q.iter()) { + if let Some(x) = t.get_mut(*n) { + *x = *v; + } + } + if t.iter().any(|x| !x.is_finite()) { + return Err(Error::Singular("the temperatures are not finite; the model is ill-conditioned".into())); + } + Ok(()) + } +} + +fn factor(a: &Csr, limits: &Limits, mesh: &HeatMesh, eqs: &Equations) -> Result { + Factor::new(&UpperColumns::from_symmetric(a), limits.max_factor_entries).map_err(|e| match e { + FactorError::TooLarge { entries, limit } => { + Error::Limit(format!("the factorization needs {entries} entries; the limit is {limit} (refine less or raise the limit)")) + } + FactorError::NotPositiveDefinite { column } => { + let label = eqs.node.get(column).and_then(|n| mesh.labels.get(*n)).copied().unwrap_or(0); + Error::Singular(format!( + "the temperature is undetermined at node {label}: no temperature boundary condition or convection reaches its part of the model" + )) + } + }) +} + +/// y = A x. +fn mul(a: &Csr, x: &[f64]) -> Vec { + let mut y = vec![0.0; a.n]; + a.mul(x, &mut y); + y +} + +/// The temperature history of a heat transfer or coupled step and what produced it. +struct Solution { + mesh: HeatMesh, + eqs: Equations, + k: Csr, + c: Option, + f: Vec, + fixed: Vec>, + times: Vec, + frames: Vec>, + /// Ṫ over the last increment (zero when steady). + rate: Vec, + factor_entries: usize, +} + +impl Solution { + fn temperature(&self) -> &[f64] { + self.frames.last().map(Vec::as_slice).unwrap_or_default() + } +} + +/// Solve the step's temperatures (steady unless the step is transient). +fn run(model: &Model, resolved: &Resolved, step: &Step, limits: &Limits) -> Result { + let transient = match step.procedure { + Procedure::HeatTransfer { transient } => transient, + _ => None, + }; + let mesh = HeatMesh::new(model, resolved, limits.max_dofs)?; + let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); + let fixed = prescribed(model, resolved, &mesh, step)?; + let eqs = Equations::new(&mesh, &graph, &fixed); + if eqs.node.is_empty() { + return Err(Error::Invalid("every temperature is prescribed; there is nothing to solve".into())); + } + let HeatSystem { k, c, f } = assemble(model, resolved, &mesh, &graph, step, transient.is_some())?; + let n = mesh.labels.len(); + let held: Vec = fixed.iter().map(|v| v.unwrap_or(0.0)).collect(); + let mut q = Vec::with_capacity(eqs.node.len()); + let (times, frames, rate, factor_entries) = match transient.zip(c.as_ref()) { + Some((tr, c)) => { + let steps = tr.increments(); + let dt = tr.end_time / steps as f64; + // (C + θΔt K) Tₙ₊₁ = (C − (1 − θ)Δt K) Tₙ + Δt f, the prescribed values held from the + // first increment on (their part A_fp T_p moves to the right-hand side). + let a = c.add_scaled(tr.theta * dt, &k)?; + let b = c.add_scaled(-(1.0 - tr.theta) * dt, &k)?; + let chol = factor(&eqs.reduce(&a), limits, &mesh, &eqs)?; + let a_held = mul(&a, &held); + let every = (steps.saturating_mul(n) / MAX_TEMPERATURE_VALUES.max(1)).saturating_add(1); + let mut t = vec![tr.initial_temperature; n]; + let (mut frames, mut times) = (vec![t.clone()], vec![0.0]); + let mut rhs = vec![0.0; n]; + let mut next = held.clone(); + for i in 1..=steps { + next.copy_from_slice(&held); + b.mul(&t, &mut rhs); + for ((r, fj), aj) in rhs.iter_mut().zip(&f).zip(&a_held) { + *r += dt * fj - aj; + } + eqs.solve(&chol, &rhs, &mut q, &mut next)?; + std::mem::swap(&mut t, &mut next); + if i % every == 0 || i == steps { + frames.push(t.clone()); + times.push(if i == steps { tr.end_time } else { i as f64 * dt }); + } + } + // After the swap `next` holds the state before the last increment. + let rate = t.iter().zip(&next).map(|(a, b)| (a - b) / dt).collect(); + (times, frames, rate, chol.entries()) + } + None => { + let rhs: Vec = f.iter().zip(mul(&k, &held)).map(|(a, b)| a - b).collect(); + let chol = factor(&eqs.reduce(&k), limits, &mesh, &eqs)?; + let mut t = held; + eqs.solve(&chol, &rhs, &mut q, &mut t)?; + (vec![0.0], vec![t], vec![0.0; n], chol.entries()) + } + }; + Ok(Solution { mesh, eqs, k, c, f, fixed, times, frames, rate, factor_entries }) +} + +/// The steady temperatures of a coupled temperature–displacement step, by node label (sorted). +pub(crate) fn temperatures(model: &Model, resolved: &Resolved, step: &Step, limits: &Limits) -> Result<(Vec, Vec)> { + let mut s = run(model, resolved, step, limits)?; + let t = s.frames.pop().unwrap_or_default(); + Ok((s.mesh.labels, t)) +} + +/// Solve a heat transfer step: the step named `step`, or the model's only one. +pub fn solve_heat(model: &Model, step: Option<&str>, limits: &Limits) -> Result { + let resolved = model.resolve()?; + let step = model.select_step(step, |p| matches!(p, Procedure::HeatTransfer { .. }), "heat transfer")?; + let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); + let s = run(model, &resolved, step, limits)?; + let t = s.temperature(); + // Reactions: the heat each prescribed node supplies, K T + C Ṫ − f. + let kt = mul(&s.k, t); + let ct = s.c.as_ref().map(|c| mul(c, &s.rate)); + let reaction: Vec = kt + .iter() + .zip(&s.f) + .zip(&s.fixed) + .enumerate() + .map(|(i, ((kt, f), fixed))| if fixed.is_some() { kt + ct.as_ref().and_then(|c| c.get(i)).copied().unwrap_or(0.0) - f } else { 0.0 }) + .collect(); + // Heat supplied: fluxes and sources, and convection at the end temperatures, Σ (f − K_film T) + // = Σ (f − K T) since conduction alone has zero row sums. + let applied_total: f64 = kt.iter().zip(&s.f).map(|(a, b)| b - a).sum(); + let reaction_total: f64 = reaction.iter().sum(); + let (flux, elements) = fluxes(&s.mesh, t)?; + // Measured against the gross heat flows, which cancel in the net sums. + let scale = reaction.iter().chain(&s.f).map(|v| v.abs()).sum::(); + if s.c.is_none() && scale > 0.0 && (applied_total + reaction_total).abs() > 1e-6 * scale { + warnings.push(format!("heat supplied and removed differ by {:.3e}", (applied_total + reaction_total).abs())); + } + let stats = Stats { dofs: s.mesh.layout.dofs(), equations: s.eqs.node.len(), factor_entries: s.factor_entries, lanczos_vectors: 0 }; + Ok(HeatResult { + step: step.name.clone(), + nodes: s.mesh.labels, + times: s.times, + frames: s.frames, + flux, + elements, + reaction, + applied_total, + reaction_total, + stats, + warnings, + }) +} + +/// The heat flux at the integration points and averaged at the nodes (one extrapolation per +/// element kind). +fn fluxes(mesh: &HeatMesh, t: &[f64]) -> Result<(Vec<[f64; 3]>, Vec)> { + let mut sum = vec![[0.0; 3]; mesh.labels.len()]; + let mut count = vec![0u32; mesh.labels.len()]; + let mut out = Vec::with_capacity(mesh.elements.len()); + let mut extrapolations: Vec<(ElementKind, Extrapolation)> = Vec::new(); + for e in &mesh.elements { + let kernel = e.kernel()?; + let at = |err| Error::element(e.label, err); + let te: Vec = e.nodes.iter().map(|n| t.get(*n).copied().unwrap_or(0.0)).collect(); + let points = kernel.flux(e.conductivity, &te).map_err(at)?; + let padded: Vec<[f64; 6]> = points.iter().map(|p| [p.flux[0], p.flux[1], p.flux[2], 0.0, 0.0, 0.0]).collect(); + let index = match extrapolations.iter().position(|(k, _)| *k == e.kind) { + Some(i) => i, + None => { + extrapolations.push((e.kind, extrapolation(kernel.ty, kernel.ty.full).map_err(at)?)); + extrapolations.len() - 1 + } + }; + let nodal = extrapolations.get(index).map(|(_, x)| x.apply(&padded)).unwrap_or_default(); + for (n, v) in e.nodes.iter().zip(&nodal) { + if let (Some(s), Some(c)) = (sum.get_mut(*n), count.get_mut(*n)) { + for k in 0..3 { + s[k] += v[k]; + } + *c += 1; + } + } + out.push(ElementFlux { label: e.label, kind: e.kind, points }); + } + let flux = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); + Ok((flux, out)) +} diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 50db5b9..756038d 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -27,6 +27,7 @@ pub mod cholesky; pub mod constraints; pub mod discrete; pub mod eigen; +pub mod heat; pub mod ordering; pub mod sparse; pub mod system; @@ -35,6 +36,7 @@ pub use analysis::{ ElementStress, FrequencyOptions, FrequencyResult, Location, MassKind, Mode, PointResult, StaticResult, Stats, StepResult, solve_frequency, solve_static, }; +pub use heat::{HeatResult, solve_heat}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 48e3962..8b5a9b1 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -9,7 +9,7 @@ use modelcraft_fem::beam::{Beam, BeamSection}; use modelcraft_fem::dense::Matrix; -use modelcraft_fem::element::{Formulation, Solid}; +use modelcraft_fem::element::{Formulation, Solid, ThermalStrain}; use modelcraft_fem::plane::{Plane, Theory}; use modelcraft_fem::shell::{Shell, ShellSection}; use modelcraft_fem::truss::{Truss, young_modulus}; @@ -35,6 +35,9 @@ pub struct Elem { pub density: Option, /// The 2D theory of an element of a 2D section; `None` for solids. pub theory: Option, + /// The expansion coefficients along X, Y, Z and the reference temperature, when the + /// material expands. + pub expansion: Option<([f64; 3], f64)>, } /// An element's kernel: a solid or a 2D element. @@ -90,11 +93,20 @@ impl Kernel<'_> { } } - /// Integration-point stresses for nodal displacements given as 3-vectors. - pub fn stresses(&self, d: &Elasticity, u: &[[f64; 3]]) -> modelcraft_fem::Result> { + /// Integration-point stresses for nodal displacements given as 3-vectors, with a thermal + /// strain when given. + pub fn stresses_with(&self, d: &Elasticity, u: &[[f64; 3]], thermal: Option<&ThermalStrain>) -> modelcraft_fem::Result> { match self { - Kernel::Solid(s) => s.stresses(d, u), - Kernel::Plane(p) => p.stresses(d, &u.iter().map(|v| [v[0], v[1]]).collect::>()), + Kernel::Solid(s) => s.stresses_with(d, u, thermal), + Kernel::Plane(p) => p.stresses_with(d, &u.iter().map(|v| [v[0], v[1]]).collect::>(), thermal), + } + } + + /// The equivalent nodal forces of a thermal strain, as 3-vectors (Z zero for 2D). + pub fn thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> modelcraft_fem::Result> { + match self { + Kernel::Solid(s) => s.thermal_load(d, thermal), + Kernel::Plane(p) => Ok(p.thermal_load(d, thermal)?.into_iter().map(|v| [v[0], v[1], 0.0]).collect()), } } @@ -121,6 +133,12 @@ impl Elem { }; k.map_err(|e| Error::element(self.label, e)) } + + /// The element's thermal strain in a temperature field (nodes the field leaves out are at + /// the reference temperature); `None` when the material does not expand or nothing changes. + pub fn thermal_strain(&self, field: &[Option]) -> Option { + thermal_strain(self.expansion, &self.nodes, field) + } } /// A truss element (rod2, rod3) ready for the kernel. @@ -134,12 +152,29 @@ pub struct Bar { pub young: f64, pub area: f64, pub density: Option, + /// The expansion coefficients along X, Y, Z and the reference temperature, when the + /// material expands. + pub expansion: Option<([f64; 3], f64)>, } impl Bar { pub fn truss(&self) -> Result> { Truss::new(self.kind, &self.xyz).map_err(|e| Error::element(self.label, e)) } + + /// The member's thermal strain in a temperature field, as for [`Elem::thermal_strain`]. + pub fn thermal_strain(&self, field: &[Option]) -> Option { + thermal_strain(self.expansion, &self.nodes, field) + } +} + +/// The thermal strain of an element of nodes `nodes` in a temperature field (nodes the field +/// leaves out are at the reference temperature); `None` when the material does not expand or +/// nothing changes. +fn thermal_strain(expansion: Option<([f64; 3], f64)>, nodes: &[usize], field: &[Option]) -> Option { + let (alpha, reference) = expansion?; + let delta: Vec = nodes.iter().map(|n| field.get(*n).copied().flatten().map_or(0.0, |t| t - reference)).collect(); + delta.iter().any(|d| *d != 0.0).then_some(ThermalStrain { alpha, delta }) } /// What a shell or beam element is made of. @@ -394,7 +429,8 @@ impl Mesh { young_modulus(elastic).map_err(|e| Error::Invalid(format!("material `{}` of section `{name}`: {e}", section.material)))?; for &label in &set.elements { let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; - trusses.push((label, el.kind, el.nodes.clone(), young, area, material.density)); + let expansion = material.thermal.and_then(|t| t.expansion.map(|a| (a.axes(), t.reference_temperature))); + trusses.push((label, el.kind, el.nodes.clone(), young, area, material.density, expansion)); } continue; } @@ -406,7 +442,8 @@ impl Mesh { Some(_) => plane::formulation(el.kind, §ion.formulation).map(|f| (f, ElementDofs::new(DofSet::PLANE))), } .map_err(|e| Error::element(label, e))?; - elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density, theory, dofs)); + let expansion = material.thermal.and_then(|t| t.expansion.map(|a| (a.axes(), t.reference_temperature))); + elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density, theory, dofs, expansion)); } } if let (Some(_), Some(name)) = (planar, model.couplings().keys().next()) { @@ -472,15 +509,15 @@ impl Mesh { warnings, discrete: crate::discrete::Discrete::default(), }; - for (label, kind, nodes, formulation, d, density, theory, _) in elements { + for (label, kind, nodes, formulation, d, density, theory, _, expansion) in elements { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); - mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory }); + mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory, expansion }); } - for (label, kind, nodes, young, area, density) in trusses { + for (label, kind, nodes, young, area, density, expansion) in trusses { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); - mesh.trusses.push(Bar { label, kind, nodes: idx, xyz: exyz, young, area, density }); + mesh.trusses.push(Bar { label, kind, nodes: idx, xyz: exyz, young, area, density, expansion }); } for (label, kind, nodes, density, section) in structural { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); @@ -713,6 +750,9 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ mesh.discrete.add_gravity(&mut f, *acceleration, mesh); } } + // Heat loads act in heat transfer; temperature fields through their thermal strain + // (add_thermal_loads). + Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. } | Load::Temperature { .. } => {} Load::BodyForce { set, vector } => { for label in resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default() { if let Some(t) = mesh.truss(*label) { @@ -779,6 +819,52 @@ fn structural_body(f: &mut [f64], mesh: &Mesh, e: &Structural, b: [f64; 3]) -> R Ok(()) } +/// The temperature field of a static step's temperature loads, per mesh node; `None` when the +/// step has none. +pub fn temperature_loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step) -> Result>>> { + let mut field: Option>> = None; + for name in &step.loads { + let Some(Load::Temperature { target, value }) = model.load(name) else { continue }; + let field = field.get_or_insert_with(|| vec![None; mesh.labels.len()]); + for &label in resolved.target_nodes(target) { + let Some(n) = mesh.node(label) else { continue }; + hold(field.get_mut(n), *value, |old| format!("node {label} is given both {old} and {value} (temperature load `{name}` and another)"))?; + } + } + Ok(field) +} + +/// Set a node's value, refusing a different one already there (`conflict` names both). +pub(crate) fn hold(slot: Option<&mut Option>, value: f64, conflict: impl FnOnce(f64) -> String) -> Result<()> { + let slot = slot.ok_or_else(|| Error::Internal("node index".into()))?; + match *slot { + Some(old) if old != value => Err(Error::Invalid(conflict(old))), + _ => { + *slot = Some(value); + Ok(()) + } + } +} + +/// Add the equivalent nodal forces of the thermal strains of a temperature field. +pub fn add_thermal_loads(f: &mut [f64], mesh: &Mesh, field: &[Option]) -> Result<()> { + for e in &mesh.elements { + let Some(th) = e.thermal_strain(field) else { continue }; + let forces = e.kernel()?.thermal_load(&e.d, &th).map_err(|err| Error::element(e.label, err))?; + for (node, v) in e.nodes.iter().zip(forces) { + add(f, &translations(mesh, *node)?, v); + } + } + for t in &mesh.trusses { + let Some(th) = t.thermal_strain(field) else { continue }; + let forces = t.truss()?.thermal_load(t.young, t.area, &th).map_err(|err| Error::element(t.label, err))?; + for (node, v) in t.nodes.iter().zip(forces) { + add(f, &translations(mesh, *node)?, v); + } + } + Ok(()) +} + /// The uniform force per volume a step's gravity and body forces put on each shell and beam /// element ([`Mesh::structural`] order): beam end forces are recovered with it. pub fn structural_body_forces(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step) -> Result> { diff --git a/crates/solver/tests/plane.rs b/crates/solver/tests/plane.rs index d50f2fb..5fda0b5 100644 --- a/crates/solver/tests/plane.rs +++ b/crates/solver/tests/plane.rs @@ -567,3 +567,76 @@ fn two_dimensional_results_are_deterministic_and_errors_name_the_problem() { let e = solve_static(&g.model, None, &Limits::default()).unwrap_err().to_string(); assert!(e.contains("singular") && e.contains("node"), "{e}"); } + +/// Give the material an expansion coefficient and heat every node by `dt` above the reference. +fn heat(g: &mut Grid, alpha: f64, dt: f64) { + let mut m = g.model.material("Mat").unwrap().clone(); + m.thermal = Some(modelcraft_model::Thermal { + expansion: Some(modelcraft_model::Directional::Isotropic(alpha)), + reference_temperature: 20.0, + ..Default::default() + }); + g.model.replace_material("Mat", m).unwrap(); + let nodes: Vec = g.param.keys().copied().collect(); + g.model.add_set("hot", Set { nodes, ..Set::default() }).unwrap(); + g.load("heat", Load::Temperature { target: Target::Set("hot".into()), value: 20.0 + dt }); +} + +#[test] +fn thermal_expansion_in_every_theory() { + // Uniform ΔT, held only by symmetry: plane stress and axisymmetric expand freely + // (u = αΔT x, no stress; εzz or the hoop strain αΔT), plane strain holds εzz = 0, so + // σzz = −EαΔT and the in-plane strain is (1 + ν)αΔT. + let (e, nu, alpha, dt) = (2e5, 0.3, 1.2e-5, 80.0); + for (kind, f, name) in variants() { + for section in [SectionKind::PlaneStress { thickness: 0.5 }, SectionKind::PlaneStrain { thickness: 0.5 }, SectionKind::Axisymmetric {}] { + let axisymmetric = matches!(section, SectionKind::Axisymmetric {}); + let mut g = grid([2, 2], kind, |p| [0.5 + 2.0 * p[0] + 0.2 * p[1], p[1]]); + g.define(section, e, nu, None, f, Procedure::Static {}); + heat(&mut g, alpha, dt); + g.symmetry("bottom", Symmetry::Y, |p| near(p[1], 0.0) && near(p[0], 0.0)); + if !axisymmetric { + g.symmetry("left", Symmetry::X, |p| near(p[0], 0.0) && near(p[1], 0.0)); + let corner = g.nodes_where(|p| near(p[0], 1.0) && near(p[1], 0.0)); + g.model.add_set("roller", Set { nodes: corner, ..Set::default() }).unwrap(); + g.bc("roller", BoundaryCondition::Symmetry { target: Target::Set("roller".into()), plane: Symmetry::Y }); + } + let r = g.solve(); + let plane_strain = matches!(section, SectionKind::PlaneStrain { .. }); + let strain = if plane_strain { (1.0 + nu) * alpha * dt } else { alpha * dt }; + for (l, i) in r.nodes.iter().zip(0..) { + let x = g.model.node(*l).unwrap().xyz; + let x0 = g.model.node(g.nodes_where(|p| near(p[0], 0.0) && near(p[1], 0.0))[0]).unwrap().xyz; + // Axisymmetric: the radial displacement is αΔT r (no rigid radial motion). + let ux = if axisymmetric { strain * x[0] } else { strain * (x[0] - x0[0]) }; + assert!((r.displacement[i][0] - ux).abs() < 1e-9, "{name} {section:?} node {l}: {:?}", r.displacement[i]); + let uy = if axisymmetric { alpha * dt * x[1] } else { strain * x[1] }; + assert!((r.displacement[i][1] - uy).abs() < 1e-9, "{name} {section:?} node {l}: {:?}", r.displacement[i]); + } + for s in &r.stress { + let zz = if plane_strain { -e * alpha * dt } else { 0.0 }; + assert!((s[0].abs() + s[1].abs() + s[3].abs()) < 1e-6 && (s[2] - zz).abs() < 1e-6, "{name} {section:?}: {s:?}"); + } + } + } +} + +#[test] +fn a_plate_held_on_every_edge_carries_the_biaxial_thermal_stress() { + // Plane stress, normal displacements held on all four edges: σxx = σyy = −EαΔT / (1 − ν). + let (e, nu, alpha, dt) = (2e5, 0.3, 1.2e-5, 50.0); + let want = -e * alpha * dt / (1.0 - nu); + for (kind, f, name) in variants() { + let mut g = grid([2, 2], kind, |p| [p[0], p[1]]); + g.define(SectionKind::PlaneStress { thickness: 1.0 }, e, nu, None, f, Procedure::Static {}); + heat(&mut g, alpha, dt); + g.symmetry("x0", Symmetry::X, |p| near(p[0], 0.0)); + g.symmetry("x1", Symmetry::X, |p| near(p[0], 1.0)); + g.symmetry("y0", Symmetry::Y, |p| near(p[1], 0.0)); + g.symmetry("y1", Symmetry::Y, |p| near(p[1], 1.0)); + let r = g.solve(); + for s in &r.stress { + assert!((s[0] - want).abs() < 1e-8 * want.abs() && (s[1] - want).abs() < 1e-8 * want.abs() && s[3].abs() < 1e-6, "{name}: {s:?}"); + } + } +} diff --git a/crates/solver/tests/thermal.rs b/crates/solver/tests/thermal.rs new file mode 100644 index 0000000..710f702 --- /dev/null +++ b/crates/solver/tests/thermal.rs @@ -0,0 +1,440 @@ +//! Analytic verification of heat transfer and thermal stress (docs/fem.md, docs/elements.md): +//! steady conduction through a slab (prescribed temperatures, a heat source, convection), +//! radial conduction through a thick cylinder, transient conduction into a semi-infinite solid +//! (erfc), free and constrained thermal expansion, a bimetal strip (Timoshenko 1925) and the +//! coupled temperature–displacement step. + +mod common; + +use common::*; +use modelcraft_model::*; +use modelcraft_solver::{Limits, solve_heat, solve_static}; + +const K: f64 = 45.0; + +/// Replace the grid's material with `elastic` plus thermal properties. +fn thermal(g: &mut Grid, elastic: Elastic, density: Option, t: Thermal) { + g.model.replace_material("Mat", Material { elastic: Some(elastic), density, thermal: Some(t), ..Material::default() }).unwrap(); +} + +fn conductor(k: f64) -> Thermal { + Thermal { conductivity: Some(Directional::Isotropic(k)), ..Thermal::default() } +} + +fn hold(g: &mut Grid, name: &str, value: f64, f: impl Fn([f64; 3]) -> bool) { + g.node_set(name, f); + g.bc(name, BoundaryCondition::Temperature { target: Target::Set(name.into()), value }); +} + +fn temperature_at(r: &modelcraft_solver::HeatResult, label: u32) -> f64 { + r.temperature()[r.nodes.binary_search(&label).unwrap()] +} + +/// Abramowitz & Stegun 7.1.26 is too coarse here; erfc by its continued fraction for large x and +/// the Taylor series for small x, both to ~1e-15. +fn erfc(x: f64) -> f64 { + if x < 2.0 { + // erf x = 2/√π Σ (−1)ⁿ x^(2n+1) / (n! (2n + 1)). + let mut sum = 0.0; + let mut term = x; + for n in 0..200 { + let t = term / (2 * n + 1) as f64; + sum += t; + term *= -x * x / (n + 1) as f64; + if t.abs() < 1e-17 { + break; + } + } + 1.0 - 2.0 / std::f64::consts::PI.sqrt() * sum + } else { + // Lentz evaluation of erfc x = e^(−x²)/√π · 1/(x + ½/(x + 1/(x + 3/2/(x + …)))). + let mut f = x; + for n in (1..200).rev() { + f = x + n as f64 / 2.0 / f; + } + (-x * x).exp() / std::f64::consts::PI.sqrt() / f + } +} + +#[test] +fn erfc_matches_reference_values() { + for (x, want) in [(0.0, 1.0), (0.5, 0.479_500_122_186_953_5), (1.0, 0.157_299_207_050_285_1), (2.5, 4.069_520_174_449_59e-4)] { + assert!((erfc(x) - want).abs() < 1e-13 * want.max(1e-3), "{x}: {} vs {want}", erfc(x)); + } +} + +#[test] +fn slab_steady_conduction_is_linear_with_the_fourier_flux() { + // A 4 × 1 × 1 slab, 20 at x = 0 and 100 at x = 4: T linear, q_x = −k ΔT / L everywhere. + for (kind, f, name) in variants() { + let mut g = grid([4, 1, 1], kind, |p| [4.0 * p[0], p[1], p[2]]); + g.define(iso(2e5, 0.3), None, f, Procedure::HeatTransfer { transient: None }); + thermal(&mut g, iso(2e5, 0.3), None, conductor(K)); + hold(&mut g, "cold", 20.0, |p| near(p[0], 0.0)); + hold(&mut g, "hot", 100.0, |p| near(p[0], 1.0)); + let r = solve_heat(&g.model, None, &Limits::default()).unwrap(); + for (l, p) in &g.param { + let t = temperature_at(&r, *l); + assert!((t - (20.0 + 80.0 * p[0])).abs() < 1e-9, "{name} node {l}: {t}"); + } + for e in &r.elements { + for p in &e.points { + assert!((p.flux[0] + K * 20.0).abs() < 1e-8 && p.flux[1].abs() < 1e-8 && p.flux[2].abs() < 1e-8, "{name}: {:?}", p.flux); + } + } + // The heat entering at the hot face leaves at the cold one: q A = k ΔT / L · 1. + let hot: f64 = g.nodes_where(|p| near(p[0], 1.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()]).sum(); + assert!((hot - K * 20.0).abs() < 1e-8, "{name}: {hot}"); + assert!(r.reaction_total.abs() < 1e-8 && r.warnings.iter().all(|w| !w.contains("differ")), "{name}"); + } +} + +#[test] +fn slab_with_a_heat_source_and_convection() { + // Uniform source Q in a slab of length L held at 0 at x = 0, convecting to T∞ at x = L: + // T(x) = −Q x²/(2k) + a x with k a − Q L = h (T∞ − T(L)) from the energy balance at x = L, + // a = (Q L + h T∞ + h Q L²/(2k)) / (k + h L). Quadratic in x: exact at the nodes of the + // quadratic elements. + let (q, h, sink, l) = (300.0, 25.0, 60.0, 2.0); + let a = (q * l + h * sink + h * q * l * l / (2.0 * K)) / (K + h * l); + let exact = |x: f64| -q * x * x / (2.0 * K) + a * x; + for (kind, f, name) in variants().into_iter().filter(|v| matches!(v.0, ElementKind::Hex20 | ElementKind::Tet10)) { + let mut g = grid([4, 1, 1], kind, |p| [l * p[0], 0.5 * p[1], 0.5 * p[2]]); + g.define(iso(2e5, 0.3), None, f, Procedure::HeatTransfer { transient: None }); + thermal(&mut g, iso(2e5, 0.3), None, conductor(K)); + hold(&mut g, "base", 0.0, |p| near(p[0], 0.0)); + g.face_set("end", |p| near(p[0], 1.0)); + g.load("film", Load::Convection { set: "end".into(), film: h, sink }); + g.load("source", Load::HeatSource { set: "all".into(), power: q }); + let r = solve_heat(&g.model, None, &Limits::default()).unwrap(); + for (lb, p) in &g.param { + let t = temperature_at(&r, *lb); + assert!((t - exact(l * p[0])).abs() < 1e-9 * exact(l).abs(), "{name} x = {}: {t} vs {}", l * p[0], exact(l * p[0])); + } + // Steady balance: what the source and the film supply leaves through the held face, + // (Q L + h (T∞ − T(L))) A. + let supplied = (q * l + h * (sink - exact(l))) * 0.25; + assert!((r.reaction_total + supplied).abs() < 1e-9 * supplied, "{name}: {} vs {supplied}", -r.reaction_total); + assert!((r.applied_total - supplied).abs() < 1e-9 * supplied, "{name}"); + } + // A flux on a face instead of convection: T = −Q x²/(2k) + (Q L + q″) x / k. + let qs = 700.0; + let mut g = grid([4, 1, 1], ElementKind::Hex20, |p| [l * p[0], 0.5 * p[1], 0.5 * p[2]]); + g.define(iso(2e5, 0.3), None, full(), Procedure::HeatTransfer { transient: None }); + thermal(&mut g, iso(2e5, 0.3), None, conductor(K)); + hold(&mut g, "base", 0.0, |p| near(p[0], 0.0)); + g.face_set("end", |p| near(p[0], 1.0)); + g.load("flux", Load::SurfaceFlux { set: "end".into(), flux: qs }); + g.load("source", Load::HeatSource { set: "all".into(), power: q }); + let r = solve_heat(&g.model, None, &Limits::default()).unwrap(); + for (lb, p) in &g.param { + let x = l * p[0]; + let want = -q * x * x / (2.0 * K) + (q * l + qs) * x / K; + assert!((temperature_at(&r, *lb) - want).abs() < 1e-9 * want.abs().max(1.0)); + } +} + +#[test] +fn thick_cylinder_radial_conduction() { + // A quarter of a hollow cylinder r = 1..2 held at 100 inside and 0 outside: + // T = 100 (1 − ln r / ln 2). The cut faces are insulated by symmetry. + for (kind, n, tol) in [(ElementKind::Hex20, [6, 6, 1], 2e-4), (ElementKind::Hex8, [16, 16, 1], 2e-3)] { + let mut g = grid(n, kind, |p| { + let (r, th) = (1.0 + p[0], std::f64::consts::FRAC_PI_2 * p[1]); + [r * th.cos(), r * th.sin(), 0.3 * p[2]] + }); + g.define(iso(2e5, 0.3), None, full(), Procedure::HeatTransfer { transient: None }); + thermal(&mut g, iso(2e5, 0.3), None, conductor(K)); + hold(&mut g, "inner", 100.0, |p| near(p[0], 0.0)); + hold(&mut g, "outer", 0.0, |p| near(p[0], 1.0)); + let r = solve_heat(&g.model, None, &Limits::default()).unwrap(); + let mut worst: f64 = 0.0; + for (l, p) in &g.param { + let rad = 1.0 + p[0]; + worst = worst.max((temperature_at(&r, *l) - 100.0 * (1.0 - rad.ln() / 2f64.ln())).abs()); + } + assert!(worst < tol * 100.0, "{kind:?}: {worst}"); + // The heat flow through the quarter: Q = (π/2) k L ΔT / ln 2. + let q_exact = std::f64::consts::FRAC_PI_2 * K * 0.3 * 100.0 / 2f64.ln(); + let q_in: f64 = g.nodes_where(|p| near(p[0], 0.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()]).sum::().abs(); + assert!((q_in - q_exact).abs() < 10.0 * tol * q_exact, "{kind:?}: {q_in} vs {q_exact}"); + } +} + +#[test] +fn transient_conduction_into_a_semi_infinite_solid() { + // A long bar at 0 whose end x = 0 is held at 1 from t = 0: T = erfc(x / (2√(κt))), κ = k/(ρc), + // until the front reaches the far end (here 2√(κt) ≤ L/4). + let (rho, c, length): (f64, f64, f64) = (7.8, 0.46, 10.0); + let kappa = K / (rho * c); + let end = (length / 8.0).powi(2) / kappa; + for (theta, increments) in [(0.5, 200usize), (1.0, 800)] { + let mut g = grid([100, 1, 1], ElementKind::Hex8, |p| [length * p[0], 0.1 * p[1], 0.1 * p[2]]); + let tr = Transient { end_time: end, increment: end / increments as f64, theta, initial_temperature: 0.0 }; + g.define(iso(2e5, 0.3), Some(rho), full(), Procedure::HeatTransfer { transient: Some(tr) }); + thermal(&mut g, iso(2e5, 0.3), Some(rho), Thermal { specific_heat: Some(c), ..conductor(K) }); + hold(&mut g, "end", 1.0, |p| near(p[0], 0.0)); + let r = solve_heat(&g.model, None, &Limits::default()).unwrap(); + assert_eq!(*r.times.last().unwrap(), end); + let mut worst: f64 = 0.0; + for (l, p) in &g.param { + let x = length * p[0]; + worst = worst.max((temperature_at(&r, *l) - erfc(x / (2.0 * (kappa * end).sqrt()))).abs()); + } + assert!(worst < 5e-3, "θ = {theta}: {worst}"); + // Frames start at the initial temperature. + assert!(r.frames[0].iter().all(|t| *t == 0.0) && r.times[0] == 0.0); + } +} + +/// A 4 × 1 × 1 bar of `kind` heated by `dt` through a temperature load, held as `hold` says. +fn heated_bar(kind: ElementKind, f: Formulation, alpha: Directional, dt: f64, poisson: f64) -> Grid { + let mut g = grid([4, 1, 1], kind, |p| [4.0 * p[0], p[1], p[2]]); + g.define(iso(2e5, poisson), None, f, Procedure::Static {}); + thermal(&mut g, iso(2e5, poisson), None, Thermal { expansion: Some(alpha), reference_temperature: 20.0, ..Thermal::default() }); + g.node_set("everything", |_| true); + g.load("heat", Load::Temperature { target: Target::Set("everything".into()), value: 20.0 + dt }); + g +} + +#[test] +fn free_thermal_expansion_is_stress_free() { + // Held only against rigid motion (3-2-1 by symmetry planes through a corner): ΔL = α L ΔT + // along each axis, no stress. + let alpha = [1.2e-5, 2.0e-5, 0.8e-5]; + for (kind, f, name) in variants() { + let mut g = heated_bar(kind, f, Directional::Orthotropic(alpha), 80.0, 0.3); + g.symmetry("x0", Symmetry::X, |p| near(p[0], 0.0)); + g.symmetry("y0", Symmetry::Y, |p| near(p[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |p| near(p[2], 0.0)); + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + for (i, l) in r.nodes.iter().enumerate() { + let x = g.model.node(*l).unwrap().xyz; + for k in 0..3 { + assert!((r.displacement[i][k] - alpha[k] * 80.0 * x[k]).abs() < 1e-12, "{name} node {l}"); + } + assert!(r.stress[i].iter().all(|s| s.abs() < 1e-6), "{name}: {:?}", r.stress[i]); + } + assert!(r.reaction.iter().flatten().all(|v| v.abs() < 1e-6), "{name}"); + } +} + +#[test] +fn constrained_bar_carries_minus_e_alpha_dt() { + // Both ends held axially, the sides free: σ_xx = −E α ΔT, the others zero, for any ν. + for (kind, f, name) in variants() { + let mut g = heated_bar(kind, f, Directional::Isotropic(1.2e-5), 50.0, 0.3); + g.symmetry("x0", Symmetry::X, |p| near(p[0], 0.0)); + g.symmetry("x1", Symmetry::X, |p| near(p[0], 1.0)); + g.symmetry("y0", Symmetry::Y, |p| near(p[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |p| near(p[2], 0.0)); + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + let want = -2e5 * 1.2e-5 * 50.0; + for s in &r.stress { + assert!((s[0] - want).abs() < 1e-8 * want.abs(), "{name}: {s:?}"); + assert!(s[1..].iter().all(|v| v.abs() < 1e-8 * want.abs()), "{name}: {s:?}"); + } + // The ends push back with σ A. + let rx: f64 = g.nodes_where(|p| near(p[0], 1.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()][0]).sum(); + assert!((rx - want).abs() < 1e-8 * want.abs(), "{name}: {rx}"); + } +} + +#[test] +fn bimetal_strip_bends_as_timoshenko_says() { + // Two equal layers (thickness t each, same E, ν = 0) with α₁, α₂, heated by ΔT, clamped at + // x = 0: curvature κ = 3 Δα ΔT / (2 h) with h = 2t (Timoshenko 1925, m = n = 1), tip rise + // κ L² / 2. + let (len, t, dt, a1, a2) = (20.0, 0.5, 100.0, 1.0e-5, 2.0e-5); + let mut g = grid([20, 1, 4], ElementKind::Hex20, |p| [len * p[0], 1.0 * p[1], 2.0 * t * p[2]]); + let all: Vec = g.model.elements().keys().copied().collect(); + let (bottom, top): (Vec, Vec) = + all.into_iter().partition(|e| g.model.element(*e).unwrap().nodes.iter().all(|n| g.param[n][2] <= 0.5 + 1e-12)); + g.model.add_set("bottom", Set { elements: bottom, ..Set::default() }).unwrap(); + g.model.add_set("top", Set { elements: top, ..Set::default() }).unwrap(); + for (name, alpha) in [("Low", a1), ("High", a2)] { + g.model + .add_material( + name, + Material { + elastic: Some(iso(2e5, 0.0)), + density: None, + thermal: Some(Thermal { expansion: Some(Directional::Isotropic(alpha)), ..Thermal::default() }), + ..Material::default() + }, + ) + .unwrap(); + } + g.model.add_section("low", Section { set: "bottom".into(), material: "Low".into(), kind: SectionKind::Solid {}, formulation: full() }).unwrap(); + g.model.add_section("high", Section { set: "top".into(), material: "High".into(), kind: SectionKind::Solid {}, formulation: full() }).unwrap(); + g.model.add_step(Step::new("Step", Procedure::Static {}), None).unwrap(); + g.fixed("root", |p| near(p[0], 0.0)); + g.node_set("everything", |_| true); + g.load("heat", Load::Temperature { target: Target::Set("everything".into()), value: dt }); + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + let kappa = 1.5 * (a2 - a1) * dt / (2.0 * t); + let want = -kappa * len * len / 2.0; + // The tip's mid-height nodes (the high-α layer on top lengthens, so the strip curls down). + let tip: Vec = g + .param + .iter() + .filter(|(_, p)| near(p[0], 1.0) && near(p[2], 0.5)) + .map(|(l, _)| r.displacement[r.nodes.binary_search(l).unwrap()][2]) + .collect(); + assert!(!tip.is_empty()); + for w in tip { + assert!((w - want).abs() < 0.02 * want.abs(), "{w} vs {want}"); + } +} + +#[test] +fn coupled_step_matches_heat_then_static() { + // A slab held at 20 and 100 at its ends, free to expand (3-2-1): the coupled step gives the + // heat step's temperatures and the static step's displacements for that field. + let (alpha, tref) = (1.2e-5, 20.0); + let build = |procedure: Procedure| { + let mut g = grid([4, 1, 1], ElementKind::Hex20, |p| [4.0 * p[0], p[1], p[2]]); + g.define(iso(2e5, 0.3), None, full(), procedure); + let t = Thermal { expansion: Some(Directional::Isotropic(alpha)), reference_temperature: tref, ..conductor(K) }; + thermal(&mut g, iso(2e5, 0.3), None, t); + g + }; + let mut c = build(Procedure::CoupledTemperatureDisplacement {}); + hold(&mut c, "cold", 20.0, |p| near(p[0], 0.0)); + hold(&mut c, "hot", 100.0, |p| near(p[0], 1.0)); + c.symmetry("x0", Symmetry::X, |p| near(p[0], 0.0)); + c.symmetry("y0", Symmetry::Y, |p| near(p[1], 0.0)); + c.symmetry("z0", Symmetry::Z, |p| near(p[2], 0.0)); + let r = solve_static(&c.model, None, &Limits::default()).unwrap(); + let temps = r.temperature.clone().unwrap(); + for (i, l) in r.nodes.iter().enumerate() { + let x = c.model.node(*l).unwrap().xyz; + assert!((temps[i] - (20.0 + 20.0 * x[0])).abs() < 1e-9); + } + // The heat step alone gives the same temperatures. + let mut h = build(Procedure::HeatTransfer { transient: None }); + hold(&mut h, "cold", 20.0, |p| near(p[0], 0.0)); + hold(&mut h, "hot", 100.0, |p| near(p[0], 1.0)); + let hr = solve_heat(&h.model, None, &Limits::default()).unwrap(); + assert_eq!(hr.temperature(), temps.as_slice()); + // A static step with the same field as temperature loads (one per plane x = const) gives the + // same displacements and stresses. + let mut s = build(Procedure::Static {}); + s.symmetry("x0", Symmetry::X, |p| near(p[0], 0.0)); + s.symmetry("y0", Symmetry::Y, |p| near(p[1], 0.0)); + s.symmetry("z0", Symmetry::Z, |p| near(p[2], 0.0)); + for k in 0..=8 { + let name = format!("plane{k}"); + s.node_set(&name, |p| near(p[0], k as f64 / 8.0)); + s.load(&name, Load::Temperature { target: Target::Set(name.clone()), value: 20.0 + 20.0 * 4.0 * k as f64 / 8.0 }); + } + let sr = solve_static(&s.model, None, &Limits::default()).unwrap(); + let scale = sr.displacement.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs())); + assert!(scale > 0.0); + for (a, b) in r.displacement.iter().flatten().zip(sr.displacement.iter().flatten()) { + assert!((a - b).abs() < 1e-12 * scale, "{a} vs {b}"); + } + for (a, b) in r.stress.iter().flatten().zip(sr.stress.iter().flatten()) { + assert!((a - b).abs() < 1e-6, "{a} vs {b}"); + } + assert!(sr.temperature.is_none()); +} + +#[test] +fn heat_steps_refuse_what_they_cannot_solve() { + // Nothing holds the temperature level. + let mut g = grid([2, 1, 1], ElementKind::Hex8, |p| p); + g.define(iso(2e5, 0.3), None, full(), Procedure::HeatTransfer { transient: None }); + thermal(&mut g, iso(2e5, 0.3), None, conductor(K)); + g.face_set("end", |p| near(p[0], 1.0)); + g.load("flux", Load::SurfaceFlux { set: "end".into(), flux: 10.0 }); + let e = solve_heat(&g.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("undetermined"), "{e}"); + // No conductivity: the model check says so. + let mut g = grid([2, 1, 1], ElementKind::Hex8, |p| p); + g.define(iso(2e5, 0.3), None, full(), Procedure::HeatTransfer { transient: None }); + hold(&mut g, "end", 1.0, |p| near(p[0], 0.0)); + let e = solve_heat(&g.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("conductivity"), "{e}"); + // Transient without a specific heat. + let mut g = grid([2, 1, 1], ElementKind::Hex8, |p| p); + g.define(iso(2e5, 0.3), Some(1.0), full(), Procedure::HeatTransfer { transient: Some(Transient::new(1.0, 0.1)) }); + thermal(&mut g, iso(2e5, 0.3), Some(1.0), conductor(K)); + hold(&mut g, "end", 1.0, |p| near(p[0], 0.0)); + let e = solve_heat(&g.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("specific heat"), "{e}"); + // Conflicting prescribed temperatures. + let mut g = grid([2, 1, 1], ElementKind::Hex8, |p| p); + g.define(iso(2e5, 0.3), None, full(), Procedure::HeatTransfer { transient: None }); + thermal(&mut g, iso(2e5, 0.3), None, conductor(K)); + hold(&mut g, "a", 1.0, |p| near(p[0], 0.0)); + hold(&mut g, "b", 2.0, |p| near(p[0], 0.0)); + assert!(solve_heat(&g.model, None, &Limits::default()).is_err()); + // A hard cap on the unknowns. + let mut g = grid([2, 1, 1], ElementKind::Hex8, |p| p); + g.define(iso(2e5, 0.3), None, full(), Procedure::HeatTransfer { transient: None }); + thermal(&mut g, iso(2e5, 0.3), None, conductor(K)); + hold(&mut g, "a", 1.0, |p| near(p[0], 0.0)); + let tight = Limits { max_dofs: 3, ..Limits::default() }; + assert!(solve_heat(&g.model, None, &tight).unwrap_err().to_string().contains("limit")); + // A static step is not a heat step. + let g = heated_bar(ElementKind::Hex8, full(), Directional::Isotropic(1e-5), 10.0, 0.3); + assert!(solve_heat(&g.model, None, &Limits::default()).is_err()); +} + +/// One rod3 or two rod2 along X over [0, `len`], heated by `dt`, fixed at x = 0, the other +/// nodes held in Y and Z, and the end held along X too when `held`. +fn heated_rod(kind: ElementKind, len: f64, alpha: [f64; 3], dt: f64, held: bool) -> modelcraft_solver::StaticResult { + let mut m = Model::default(); + for (label, x) in [(1, 0.0), (2, len), (3, 0.5 * len)] { + m.set_node(label, [x, 0.0, 0.0]).unwrap(); + } + // One rod3, or two rod2 in a row. + let elements = if kind == ElementKind::Rod3 { + m.set_element(1, kind, vec![1, 2, 3]).unwrap(); + vec![1] + } else { + m.set_element(1, kind, vec![1, 3]).unwrap(); + m.set_element(2, kind, vec![3, 2]).unwrap(); + vec![1, 2] + }; + m.add_set("rod", Set { elements, nodes: vec![1, 2, 3], ..Set::default() }).unwrap(); + m.add_set("a", Set { nodes: vec![1], ..Set::default() }).unwrap(); + m.add_set("b", Set { nodes: vec![2], ..Set::default() }).unwrap(); + m.add_set("c", Set { nodes: vec![3], ..Set::default() }).unwrap(); + let thermal = Thermal { expansion: Some(Directional::Orthotropic(alpha)), ..Thermal::default() }; + m.add_material("Mat", Material { elastic: Some(iso(2e5, 0.3)), thermal: Some(thermal), ..Material::default() }).unwrap(); + m.add_section("rod", Section { set: "rod".into(), material: "Mat".into(), kind: SectionKind::Rod { area: 4.0 }, formulation: full() }).unwrap(); + let mut step = Step::new("Step", Procedure::Static {}); + let lateral = [None, Some(0.0), Some(0.0)]; + for (name, v) in [("a", [Some(0.0); 3]), ("b", if held { [Some(0.0); 3] } else { lateral }), ("c", lateral)] { + let values = [v[0], v[1], v[2], None, None, None]; + m.add_boundary_condition(name, BoundaryCondition::Displacement { target: Target::Set(name.into()), values }).unwrap(); + step.boundary_conditions.insert(name.into()); + } + m.add_load("warm", Load::Temperature { target: Target::Set("rod".into()), value: dt }).unwrap(); + step.loads.insert("warm".into()); + m.add_step(step, None).unwrap(); + solve_static(&m, None, &Limits::default()).unwrap() +} + +#[test] +fn trusses_take_thermal_strain_along_their_axis() { + // Along X the member's thermal strain is αₓ ΔT (the projection Σ αₖ tₖ² is checked on the + // kernel): free it lengthens by αₓ ΔT L stress-free; held at both ends it carries −E αₓ ΔT. + let alpha = [1.0e-5, 2.0e-5, 3.0e-5]; + let (e, area, dt, len) = (2e5, 4.0, 50.0, 3.0); + let sigma = -e * alpha[0] * dt; + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let r = heated_rod(kind, len, alpha, dt, true); + for p in r.elements.iter().flat_map(|el| &el.points) { + assert!((p.stress[0] - sigma).abs() < 1e-9 * sigma.abs(), "{kind:?}: {:?}", p.stress); + } + let reaction = r.reaction[r.nodes.binary_search(&2).unwrap()][0]; + assert!((reaction - sigma * area).abs() < 1e-9 * (sigma * area).abs(), "{kind:?}: {reaction}"); + let r = heated_rod(kind, len, alpha, dt, false); + let u = r.displacement[r.nodes.binary_search(&2).unwrap()][0]; + assert!((u - alpha[0] * dt * len).abs() < 1e-12, "{kind:?}: {u}"); + assert!(r.elements.iter().flat_map(|el| &el.points).all(|p| p.stress[0].abs() < 1e-9)); + } +} diff --git a/crates/ui-egui/src/results.rs b/crates/ui-egui/src/results.rs index 0469818..0a8315c 100644 --- a/crates/ui-egui/src/results.rs +++ b/crates/ui-egui/src/results.rs @@ -69,16 +69,19 @@ impl Plot { pub fn mode(&self) -> Option<(usize, f64)> { match self.result.kind() { Kind::Frequency(f) => f.modes.get(self.pick.mode.checked_sub(1)?).map(|m| (self.pick.mode, m.frequency)), - Kind::Static(_) | Kind::Explicit(_) | Kind::Fluid(_) => None, + Kind::Static(_) | Kind::Explicit(_) | Kind::Fluid(_) | Kind::Heat(_) => None, } } - /// An explicit step's field frame shown, its time and the number of frames. + /// An explicit or transient heat step's frame shown, its time and the number of frames. pub fn frame(&self) -> Option<(usize, f64, usize)> { - let Kind::Explicit(e) = self.result.kind() else { return None }; - let n = e.frames.len(); + let (n, time): (usize, &dyn Fn(usize) -> Option) = match self.result.kind() { + Kind::Explicit(e) => (e.frames.len(), &|k| e.frames.get(k).map(|f| f.time)), + Kind::Heat(h) if h.times.len() > 1 => (h.times.len(), &|k| h.times.get(k).copied()), + _ => return None, + }; let k = self.pick.frame.unwrap_or(n.checked_sub(1)?); - Some((k, e.frames.get(k)?.time, n)) + Some((k, time(k)?, n)) } /// The node's position in the plot, if it is drawn. @@ -107,7 +110,7 @@ pub fn swing(result: &StepResult, frame: Option) -> f64 { match result.kind() { Kind::Frequency(_) => a.sin(), Kind::Static(_) => 0.5 * (1.0 - a.cos()), - Kind::Explicit(_) | Kind::Fluid(_) => 1.0, + Kind::Explicit(_) | Kind::Fluid(_) | Kind::Heat(_) => 1.0, } } diff --git a/docs/elements.md b/docs/elements.md index f729ba0..4eee6b6 100644 --- a/docs/elements.md +++ b/docs/elements.md @@ -83,6 +83,6 @@ hex8/20/27 with first and second derivatives). | Cohesive (2D, 3D) | Adhesive layers, delamination | planned | | Gasket | Through-thickness gasket behaviour | planned | | Infinite elements | Unbounded domains | planned | -| Heat transfer (all shapes), coupled temperature–displacement | Thermal and thermal-stress analyses | in progress (FEM foundation, phase B); temperature DOF kind done | +| Heat transfer (all shapes), coupled temperature–displacement | Thermal and thermal-stress analyses | done for solids (tet4/10, wedge6, hex8/20: steady and transient conduction, flux, convection, sources; coupled temperature–displacement); thermal strain in solids, trusses, 2D elements (plane stress, plane strain, axisymmetric) and explicit dynamics (one-point hex8, nodal-pressure tet4, with plasticity); 2D/shell/beam conduction and shell/beam thermal strain planned | | Acoustic | Fluid acoustics, structural–acoustic coupling | planned | | Rigid elements | Rigid bodies from meshes | planned (M7) | diff --git a/docs/features.md b/docs/features.md index 0db067a..520da15 100644 --- a/docs/features.md +++ b/docs/features.md @@ -121,7 +121,7 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md). | Decks | Nastran Deck — Write Nastran bulk data (SOL 101, 103) | ready | `deck.write {"format":"nastran"}` | | Decks | LS-DYNA Deck — Write an LS-DYNA keyword file | ready | `deck.write {"format":"lsdyna"}` | | Decks | Read Deck — Read a deck, keeping keywords ModelCraft does not know | coming soon (M9) | | -| Implicit | Static — Linear static solve of a step with ModelCraft's own solver | ready | `solve.static` | +| Implicit | Static — Linear static (or coupled temperature–displacement) solve of a step with ModelCraft's own solver; heat transfer in the drop-down | ready | `solve.static`, `solve.heat` | | Implicit | Frequency — Natural frequencies and mode shapes of a frequency step | ready | `solve.frequency` | | Explicit | Explicit — Transient dynamics of an explicit step by central differences (contact and erosion coming in M10) | ready | `solve.explicit` | | Explicit | Blast & Impact — EOS-driven fluids with fluid–structure contact (ALE or SPH) | coming soon (M14) | | diff --git a/docs/fem.md b/docs/fem.md index 3e2fb41..175925f 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -101,3 +101,35 @@ of the family's degree; on a curved map, linear fields exact with zero Hessian, physical Hessian of a quadratic interpolant consistent with finite differences of its physical gradient. The solver and explicit migration to `DofLayout` was checked bit for bit: every result of the solver, explicit, engine and io test suites hashes identically before and after. + +## Heat transfer and thermal strain + +A heat conduction element carries one temperature per node (`heat::HEAT_DOFS`, kind `T`); +`heat::Conduction` gives the conductivity ∫ ∇Nᵀ k ∇N dV (full rule), the capacity ∫ ρc N Nᵀ dV +and a source ∫ N Q dV (mass rule), and the heat flux −k∇T at the integration points; +`heat::face_weights` and `face_film` integrate surface flux and convection on tri3/tri6/quad4/quad8 +faces. The solver (`modelcraft_solver::heat`) numbers the temperatures with +`DofLayout::uniform(…, DofSet::TEMPERATURE)` and assembles with the same `sparse::add_block`. + +Thermal strain ε_th = α (T − T_ref) along the material axes enters a solid as the load +`Solid::thermal_load(d, &ThermalStrain { alpha, delta })` = ∫ Bᵀ D ε_th dV and its stress as +`Solid::stresses_with(d, u, Some(&thermal))` = D (ε − ε_th); `delta` is the nodal T − T_ref. +2D elements have the same pair (`Plane::thermal_load`, `Plane::stresses_with`): plane stress +leaves εzz free (σzz = 0), plane strain holds it (σzz = −D₃·ε_th), axisymmetric elements strain +the hoop direction too. Trusses (rod2, rod3) take the axial part Σ αₖ tₖ² ΔT along their tangent (`Truss::thermal_load`, +`Truss::stresses_with`). Shells and beams will provide the same two functions for their strains; +springs, dashpots and point masses have no thermal strain. + +Explicit dynamics takes a step's temperature loads as a thermal strain rate per element, +α ΔT / end time (the change reached linearly over the step, like a prescribed displacement, so +there is no thermal shock), subtracted from each increment's strain increment before the +constitutive update, so it feeds J2 and Johnson–Cook plasticity alike; isotropic expansion only +(the rate is the same in any frame). Explicit trusses take no thermal strain yet (a warning +says so). + +Model side: `Material::thermal` (`conductivity`, `specific_heat`, `expansion`, each isotropic +or per axis, and `reference_temperature`), `BoundaryCondition::Temperature`, loads +`SurfaceFlux`, `Convection`, `HeatSource` (heat transfer) and `Temperature` (a temperature field +in a static step), procedures `HeatTransfer { transient }` and `CoupledTemperatureDisplacement`. +Commands: `solve.heat`, and `solve.static` for coupled steps. The deck writers refuse heat +transfer by name for now. diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 476b189..d75ddd9 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1119,3 +1119,25 @@ update are written from (cited again in the module docs): and undeformed length of a flat-ended projectile). Published experimental Taylor-bar lengths (Wilkins & Guinan 1973) are not used: that needs the owner's sign-off for a tagged benchmark suite. + +## Part K: Heat transfer and thermal stress, as implemented + +- **Conduction.** Galerkin form of ρc Ṫ = ∇·(k∇T) + Q with surface flux and convection: + O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, *The Finite Element Method: Its Basis and + Fundamentals*, 7th ed., Butterworth-Heinemann 2013, §3.2; T. J. R. Hughes, *The Finite Element + Method*, Dover 2000, ch. 1. Kernels in `modelcraft_fem::heat`, solver in `modelcraft_solver::heat`. +- **Time integration.** The generalized trapezoidal (θ) family, unconditionally stable for + θ ≥ ½ (Hughes 2000, §8.2): backward Euler (θ = 1, the default) and Crank–Nicolson (θ = ½). +- **Thermal strain.** Initial-strain load ∫ Bᵀ D ε_th dV and stress D(ε − ε_th) + (Zienkiewicz, Taylor & Zhu 2013, §6.2.3); incompatible modes condensed with the same load. +- **Coupled temperature–displacement.** One-way coupling (temperature strains the solid) solved + in turn, which is the exact monolithic solution of linear quasi-static thermoelasticity + without thermoelastic dissipation (B. A. Boley, J. H. Weiner, *Theory of Thermal Stresses*, + Wiley 1960, ch. 1). +- **Verification** (`crates/solver/tests/thermal.rs`, analytic): steady slab conduction, a slab + with a heat source and convection or flux (quadratic profile), radial conduction through a + thick cylinder (T ∝ ln r), the semi-infinite solid T = erfc(x / 2√(κt)) (H. S. Carslaw, + J. C. Jaeger, *Conduction of Heat in Solids*, 2nd ed., Oxford 1959, §2.4), free and + constrained thermal expansion (σ = −EαΔT), the bimetal strip (S. Timoshenko, *Analysis of + bi-metal thermostats*, J. Opt. Soc. Am. 11(3):233–255, 1925), and the coupled step against + the heat step followed by the static step.