diff --git a/ROADMAP.md b/ROADMAP.md index bb965ca..f9d3a74 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -41,7 +41,7 @@ pressure) is a degree of freedom of the existing DOF framework, not a separate s | M12 CFD: lattice Boltzmann | D3Q19/D3Q27 on grids voxelised from the face-tagged surface; TRT, regularised and cumulant collision; interpolated walls, inlets and outlets; LES (Smagorinsky, WALE) with wall models; thermal LBM and conjugate heat transfer with the structure's thermal solver; non-Newtonian viscosity (power-law, Carreau–Yasuda, Cross, Herschel–Bulkley, Bingham) and temperature-dependent viscosity; porous media (Brinkman, Darcy–Forchheimer); multiphase and free surface (M13); validation by the lid-driven cavity, backward-facing step, cylinder wake (Strouhal) and turbulent channel; octree refinement shared with the hex mesher; wgpu backend (deterministic per device, checked against the CPU) (CPU core: done — `fieldcraft-lbm` with D3Q19 TRT and Guo forcing, halfway bounce-back and moving walls, slip and periodic faces, Zou–He velocity and pressure faces, exact voxelisation of obstacle bodies, momentum-exchange forces per face; fluid materials, fluid domains and fluid steps in the model; `model.fluid_domain.*`, `solve.fluid`; verified by Poiseuille channel and pipe, Couette, Taylor–Green decay, sphere-array Stokes drag and mass conservation; Schäfer–Turek 2D-1 as a benchmark) | started | | M13 Fluid–structure interaction | Immersed-boundary coupling of the lattice Boltzmann fluid with the explicit or implicit structure; interface force transfer; sub-cycling; Aitken then interface quasi-Newton coupling; free-surface flow | planned | | M14 Blast and impact fluids | EOS-driven fluids in the explicit solver: ALE multi-material or SPH, with fluid–structure contact | planned | -| M15 Nonlinear implicit core | Newton–Raphson with line search, automatic increments and cutbacks, convergence on force, displacement and energy; geometric nonlinearity (total and updated Lagrangian, large rotations for shells and beams, follower pressure); every `fem` material behind one interface with its consistent tangent (plasticity, damage and thermal strain from the explicit solver reused implicitly); arc-length (Riks) for snap-through; linear buckling (eigen) with preload stress stiffening, and imperfections from buckling modes; nonlinear stabilisation; restart. Verified by the elastica, Lee's frame, the snap-through arch, the hinged cylindrical shell and plastic collapse loads (solids core: done — `nonlinear_static` steps with automatic increments and cutbacks, Newton–Raphson with the consistent tangent, force, displacement and energy convergence and an optional line search; a material point interface (Kirchhoff stress and spatial tangent from F, with history) with St Venant–Kirchhoff and compressible neo-Hookean; total Lagrangian tet4, tet10, wedge6, hex8, hex20 (full and reduced), hybrid hex8/hex20/tet10 by the three-field mixed form, plane strain and axisymmetric elements; follower pressure with its symmetrised load stiffness; Cauchy stress and logarithmic strain, increment history and frames, `solve.nonlinear`; Abaqus and LS-DYNA decks; verified by the patch and objectivity tests, uniaxial stretch, the elastica and tube inflation by a follower pressure (finite-strain Cook's membrane pending its transcribed reference). Next: J2 plasticity and Mooney–Rivlin on the interface, thermal strain, incompatible modes (EAS), shells and beams, arc-length) | started | +| M15 Nonlinear implicit core | Newton–Raphson with line search, automatic increments and cutbacks, convergence on force, displacement and energy; geometric nonlinearity (total and updated Lagrangian, large rotations for shells and beams, follower pressure); every `fem` material behind one interface with its consistent tangent (plasticity, damage and thermal strain from the explicit solver reused implicitly); arc-length (Riks) for snap-through; linear buckling (eigen) with preload stress stiffening, and imperfections from buckling modes; nonlinear stabilisation; restart. Verified by the elastica, Lee's frame, the snap-through arch, the hinged cylindrical shell and plastic collapse loads (solids core: done — `nonlinear_static` steps with automatic increments and cutbacks, Newton–Raphson with the consistent tangent, force, displacement and energy convergence and an optional line search; a material point interface (Kirchhoff stress and spatial tangent from F, with history) with St Venant–Kirchhoff and compressible neo-Hookean; total Lagrangian tet4, tet10, wedge6, hex8, hex20 (full and reduced), hybrid hex8/hex20/tet10 by the three-field mixed form, plane strain and axisymmetric elements; follower pressure with its symmetrised load stiffness; Cauchy stress and logarithmic strain, increment history and frames, `solve.nonlinear`; Abaqus and LS-DYNA decks; verified by the patch and objectivity tests, uniaxial stretch, the elastica and tube inflation by a follower pressure (finite-strain Cook's membrane pending its transcribed reference; J2 plasticity on the interface: Simo's exponential return map at finite strain and the radial return at small strain, both with consistent tangents and the return shared with explicit dynamics, bilinear isotropic and kinematic (small strain), tabulated and Johnson–Cook (rate-independent) hardening, hybrid elements included, load paths for cyclic loading, a tangent predictor for path-dependent materials, verified by uniaxial tension to large strain, Hill's thick cylinder (plastic zone and collapse pressure), a cantilever's collapse and the Bauschinger effect, with quadratic Newton convergence. Next: Mooney–Rivlin and Yeoh, thermal strain, incompatible modes (EAS), shells and beams, arc-length) | started | | M16 Structural dynamics | Implicit transient (Newmark, HHT-α, generalised-α); modal superposition transient; steady-state harmonic (direct and modal), with Rayleigh, modal, structural and material damping; random vibration from PSDs and response spectrum (SRSS, CQC); complex eigenvalues with damping; rotating machinery (Campbell diagrams, gyroscopic terms); substructures (Craig–Bampton) | planned | | M17 Constitutive library and material library | Elastic: transversely isotropic, fully anisotropic (21 constants), temperature-dependent. Elastomers: hyperelastic neo-Hookean, Mooney–Rivlin, Yeoh, Ogden, Arruda–Boyce, Gent, Marlow, with Mullins softening and finite-strain viscoelasticity (Prony series; parallel networks after Bergström–Boyce), foams. Viscoelasticity: small- and finite-strain Prony series, time–temperature superposition (WLF, Arrhenius), storage and loss moduli for harmonic steps. Viscoplasticity and creep: Norton power law, time and strain hardening, Perzyna, Anand, Chaboche combined hardening. Plasticity: Hill anisotropic, Drucker–Prager, Mohr–Coulomb, Cam-clay, crushable foam, concrete damaged plasticity, Gurson. Fluids: Newtonian and non-Newtonian viscosity shared with M12. Material library: curated, cited data sets (steels, aluminium and titanium alloys, polymers, elastomers, composites, tissues) with units and temperature dependence; calibration from test data (uniaxial, biaxial, shear, relaxation, DMA) by least squares with Drucker stability checks; user materials as a safe Rust trait | planned | | M18 Composites and anisotropy | Material orientations (rectangular, cylindrical, spherical, from geometry edges and faces, per-element fields, draping of fibres over curved parts); layups on shells and continuum shells, ply stacks on swept hex; classical lamination theory checks; ply stresses and strains through the thickness; failure indices (Tsai–Wu, Tsai–Hill, Hashin, Puck, LaRC) and progressive damage; homogenisation (rule of mixtures, Mori–Tanaka, RVEs with periodic boundary conditions); delamination through M20's cohesive zones. Verified by laminate closed forms and the NAFEMS composite benchmarks | planned | @@ -76,7 +76,7 @@ What exists and what is missing, by area; each gap is in a milestone above. | Area | Today | Missing | |---|---|---| -| Solution procedures | Linear static, nonlinear static (solids with geometric nonlinearity), frequency (Lanczos), steady and transient heat, coupled temperature–displacement, explicit dynamics, LBM flow | Nonlinear shells and beams, plasticity in implicit steps, arc-length, buckling (M15); implicit, harmonic, random and spectrum dynamics (M16) | +| Solution procedures | Linear static, nonlinear static (solids with geometric and material nonlinearity: hyperelasticity, J2 plasticity), frequency (Lanczos), steady and transient heat, coupled temperature–displacement, explicit dynamics, LBM flow | Nonlinear shells and beams, arc-length, buckling (M15); implicit, harmonic, random and spectrum dynamics (M16) | | Contact, friction, adhesion | Explicit planar rigid walls (frictionless); couplings, springs, dashpots | Implicit and explicit deformable contact, tie, friction laws, thermal contact, wear (M7); cohesive adhesion and debonding (M20) | | Materials | Isotropic and orthotropic elasticity, density, thermal, J2 plasticity (bilinear, table, Johnson–Cook), Johnson–Cook damage and erosion, compressible neo-Hookean hyperelasticity, Newtonian fluids | Elastomers beyond neo-Hookean (Mooney–Rivlin and the rest), viscoelasticity, viscoplasticity and creep, anisotropic and pressure-dependent plasticity, non-Newtonian viscosity, a cited material library and calibration (M17) | | Composites and anisotropy | Orthotropic constants | Orientations, layups, ply results, failure criteria, homogenisation (M18) | diff --git a/docs/fem.md b/docs/fem.md index 4295d5a..7e3d10f 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -172,6 +172,22 @@ converges, so a cut-back increment leaves no trace. The materials: | Elastic (isotropic or orthotropic), no option | St Venant–Kirchhoff: S = D E, E = ½(FᵀF − I). Right for large rotations with small strains; it softens without bound in large compression. | | Isotropic elastic with `"hyperelastic": {"type": "neo_hookean"}` | Compressible neo-Hookean with the isochoric–volumetric split, W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)², μ and κ from E and ν (C₁₀ = μ/2, D₁ = 2/κ). Its small-strain limit is the linear material, which linear steps use. | +| Isotropic elastic with `"plastic"` (bilinear isotropic or kinematic, table, Johnson–Cook) | J2 plasticity (`fieldcraft_fem::elastoplastic`): at finite strain Simo's multiplicative model (CMAME 99, 1992) with a Hencky elastic energy, τ = K ln J I + 2G dev ½ ln bₑ, the history Cₚ⁻¹ and the equivalent plastic strain, integrated by the exponential return map on the principal logarithmic strains; at small strain the additive radial return with the back stress. Both call `J2::radial_return`, the update of explicit dynamics, and its consistent tangent `J2::tangent`. | + +**Plasticity's tangent.** With Cₚ⁻¹ frozen over an iteration, L_v τ = 𝔻ᵉᵖ : ½ (∂ ln b/∂b) : (d b + b d) − (d τ + τ d), 𝔻ᵉᵖ the algorithmic tangent of the return +(K I⊗I + 2Gθ ℙ − 2Gθ̄ n⊗n, Simo & Hughes 1998 Box 3.2) and the derivative of the logarithm the +Daleckii–Krein divided differences θ_AB = artanh(r)/r, r = (λ²_A − λ²_B)/(λ²_A + λ²_B), in the +principal frame of the trial bₑ: smooth through equal stretches, so uniaxial and undeformed +states need no special case. The tangent is symmetric (isotropic hardening) and checked against +central differences of the stress and of every element's forces from a plastic history. Hybrid +elements take the deviatoric return and the volume K/2 (ln J)². Refused by name in nonlinear +steps: kinematic hardening at finite strain (a back stress is not coaxial with bₑ, which the +principal-frame return needs; it works with small strain), the Johnson–Cook strain-rate term (a +static step has no time; C = 0 is the rate-independent law), and incompatible-mode elements of a +plastic material. Damage, erosion and adiabatic heating act in explicit dynamics only (a warning +says so). A hardening slope below −3G leaves the point without a tangent: the increment cuts +back and the step stops at the minimum increment, naming the element. + "Isotropic elastic" keeps meaning the same thing everywhere: the neo-Hookean model is chosen only by the explicit option, and explicit dynamics refuses a hyperelastic material by name. @@ -190,9 +206,11 @@ differences of its forces. | plane stress | refused (σ_zz = 0 needs a local iteration for the thickness) | | trusses, shells, beams | refused (large rotations of structural elements come later in M15) | -Without nonlinear geometry (`"nonlinear_geometry": false`) the elements are the linear ones: -every material is elastic today, so such a step is the linear solution reached by Newton (the -path the next materials, plasticity first, take with small strains). +Without nonlinear geometry (`"nonlinear_geometry": false`) an elastic element is the linear one +(K u), and an element of a plastic material takes the same loops with small strain +(`Kinematics::Small`: ε = sym ∇u, the reference gradients, f = ∫ Bᵀσ dV, K = ∫ Bᵀ𝔻B dV with the +consistent tangent, no geometric stiffness; hybrid elements project 1 + tr ε and K tr ε), so an +elastic model is still the linear solution reached by Newton. **Follower pressure** (`fieldcraft_fem::face::follower_pressure`, `follower_edge_pressure`): fᵢ = −p ∫ Nᵢ (x,ξ × x,η) dξ dη on the current face, with its load stiffness @@ -207,6 +225,23 @@ the rate of convergence of open loaded surfaces; it is revisited with the symmet (LDLᵀ) factorization that arc-length and Taylor–Hood need. Without nonlinear geometry a pressure is a dead load. +**Load path.** A step's `load_path`, `[[fraction, factor], …]` from `[0, 0]` to fraction 1, makes +the load factor piecewise linear in the step fraction (without one it equals the fraction), so +loads and prescribed displacements can come down or reverse within a step: cyclic loading. +Increments end at every point of the path, and the extrapolated predictor scales with the change +of the load factor, so it turns with the path. LS-DYNA writes the path as the step's curve; +Abaqus refuses it by name for now (it would be an amplitude on every load and boundary +condition). + +**Predictors.** An increment starts from the linear extrapolation of the last one; with +path-dependent materials an increment without one (the first, after a turn of the path to zero +change) starts from the tangent predictor: one solve with the converged state's tangent for the +change of the loads and of the prescribed displacements, so a prescribed motion is spread into the +interior before the material sees it (started from the boundary's jump alone, the elements next +to it would flow plastically far beyond the solution). Its work is then the energy criterion's +reference. A first tangent that is not positive definite is a cutback when the undeformed +model's tangent is (the strain of the increment, not a model that can move freely). + **Newton–Raphson** with the consistent tangent refactored every iteration, an optional line search (off by default: on large elastic rotations the full step converges in fewer iterations), and convergence when the force residual (∞-norm over the largest force seen), the @@ -219,8 +254,10 @@ tangent is not positive definite, or the iterations run out; the step stops with and the reason at the minimum increment, the increment cap, or `Limits::max_newton_iterations`. **Results**: the final displacements and reactions, the **Cauchy stress** and the -**logarithmic strain** (ln V, engineering shears) at the integration points and averaged at the -nodes, every increment's load factor, iterations and residuals, and a displacement frame per +**logarithmic strain** (ln V, engineering shears; the small strain without nonlinear geometry) at +the integration points and averaged at the nodes, the **equivalent plastic strain** at the +integration points of plastic elements (`solve.probe` with `elements`, the largest in the step +summary), every increment's load factor, iterations and residuals, and a displacement frame per converged increment (`results.frame`, the slider in the results panel). **Couplings** stay the linear constraints of the linear solver (small rotations of the @@ -229,9 +266,11 @@ centrifugal loads are refused by name in nonlinear steps. **Decks**: Abaqus `*STEP, NLGEOM=YES|NO, INC=` with `*STATIC` (initial, period 1, minimum, maximum) and `*CONTROLS, PARAMETERS=LINE SEARCH`; the neo-Hookean option as `*HYPERELASTIC, NEO -HOOKE`. LS-DYNA: the ramped implicit run with full Newton (`ilimit` 1, `maxref`) and +HOOKE`; plasticity as for explicit dynamics (`*PLASTIC`); a load path, and a plastic deck of +several general steps (Abaqus carries the plastic state from one to the next, every model step +starts undeformed), are refused. LS-DYNA: the ramped implicit run with full Newton (`ilimit` 1, `maxref`) and `*CONTROL_IMPLICIT_AUTO`; hyperelastic materials in nonlinear steps and steps without nonlinear -geometry are refused. Nastran refuses nonlinear steps (SOL 400) by name. Iteration caps and +geometry are refused; a load path is the step's curve. Nastran refuses nonlinear steps (SOL 400) by name. Iteration caps and tolerances are left to each solver's own controls, whose norms differ from ours. **Verification** (`crates/solver/tests/nonlinear.rs`, `crates/fem/src/{constitutive,finite,face}.rs`, @@ -242,4 +281,16 @@ test on a distorted mesh (all solids, hybrid ones included); a 90° rigid rotati the elastica (Bisshopp–Drucker) at P L²/(E I) = 1, 2, 5, 10 with hex20 (0.2–0.5 %); the inflation of a thick neo-Hookean tube by a follower pressure (Rivlin's closed form); Newton's quadratic convergence; the linear limit; a finite-strain Cook's membrane is pending its transcribed reference; and the failures (an element inverted by -crushing stops at the minimum increment, naming it; the caps; what is refused). +crushing stops at the minimum increment, naming it; the caps; what is refused). Plasticity +(`crates/solver/tests/plasticity.rs`, `crates/fem/src/elastoplastic.rs`, `tests/oracle/` +`nl-j2-uniaxial`, `j2-bauschinger`, `plastic-thick-cylinder`, `plastic-hinge-cantilever`): the +consistent tangents against central differences from a plastic history (material point and +every element, finite and small strain, full and hybrid); uniaxial tension to a logarithmic +strain of 0.2 in every element kind, hybrid ones included (Kirchhoff stress exact to 10⁻⁷, and +0.13 % from explicit dynamics' hypoelastic Cauchy stress, the elastic volume change between the +two laws); the Bauschinger effect on a load path back to zero strain (exact); Hill's thick +cylinder, the bore displacement at a plastic radius of 1.5a (0.009 %) and the collapse pressure +(2·10⁻⁸); a cantilever's collapse 3.1 % above beam theory's hinge, converged under refinement; +Newton's order with the consistent tangent 1.9 (small strain) and 2.3 (finite strain) in a +plastic cantilever; and the refusals and the softening material that stops at the minimum +increment. diff --git a/docs/oracle.md b/docs/oracle.md index 79288f9..35a43b9 100644 --- a/docs/oracle.md +++ b/docs/oracle.md @@ -2,7 +2,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the policy is [verification.md](verification.md). Every case is our own model, built with FieldCraft commands from the problem's published description, and compared with the reference value its `expected.json` cites, within the stated tolerance and for the stated reason. Accepted results are frozen in each case's `fixture.json`: a later change that moves one beyond its band (10⁻⁹ of its size) fails CI even while the reference still passes. No reference value was produced with a commercial product. -**142 of 148 quantities pass** (6 known failure(s)) in 40 active cases: 31 theory, 9 published benchmark, 0 experiment; 10 pending cases wait for capabilities still being built. +**149 of 155 quantities pass** (6 known failure(s)) in 44 active cases: 35 theory, 9 published benchmark, 0 experiment; 9 pending cases wait for capabilities still being built. ## Theory: closed-form solutions, patch tests and mesh invariants @@ -31,6 +31,8 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the | [free-beam-frequencies](../tests/oracle/free-beam-frequencies/) | sixth frequency (the last rigid-body mode) | 0 | ≈ 0 | ≈ 0 | ± 8.1264e-5 | pass | matches | | [free-beam-frequencies](../tests/oracle/free-beam-frequencies/) | first free–free bending frequency (mode 7) | 0.081264 | 0.080618 | -0.795% | 1.5% | pass | matches | | [free-beam-frequencies](../tests/oracle/free-beam-frequencies/) | total mass | 20 | 20 | < 1e-7% | 1e-12 rel. | pass | matches | +| [j2-bauschinger](../tests/oracle/j2-bauschinger/) | axial stress back at zero strain after a pull to ε_a = 0.01 | -247.52 | -247.52 | < 1e-7% | 1e-6 rel. | pass | blessed | +| [j2-bauschinger](../tests/oracle/j2-bauschinger/) | accumulated equivalent plastic strain | 0.016089 | 0.016089 | < 1e-7% | 1e-6 rel. | pass | blessed | | [j2-uniaxial-explicit](../tests/oracle/j2-uniaxial-explicit/) | axial Cauchy stress at a logarithmic strain of 0.2 | 643.56 | 643.57 | 0.000174% | 1% | pass | matches | | [j2-uniaxial-explicit](../tests/oracle/j2-uniaxial-explicit/) | equivalent plastic strain | 0.19678 | 0.19678 | -0.00000113% | 1% | pass | matches | | [kirsch-plate](../tests/oracle/kirsch-plate/) | σ_xx at the top of the hole (0, 1, 0) over the remote stress | 3 | 3.0121 | 0.402% | 2% | pass | matches | @@ -57,6 +59,8 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip shortening δ_h/L at P L²/(E I) = 5 | 0.38763 | 0.38955 | 0.495% | 2% | pass | matches | | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip deflection δ_v/L at P L²/(E I) = 10 | 0.81061 | 0.81421 | 0.444% | 1% | pass | matches | | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip shortening δ_h/L at P L²/(E I) = 10 | 0.555 | 0.558 | 0.541% | 2% | pass | matches | +| [nl-j2-uniaxial](../tests/oracle/nl-j2-uniaxial/) | axial Cauchy stress at a logarithmic strain of 0.2 | 642.74 | 642.74 | < 1e-7% | 1e-6 rel. | pass | blessed | +| [nl-j2-uniaxial](../tests/oracle/nl-j2-uniaxial/) | equivalent plastic strain | 0.19678 | 0.19678 | < 1e-7% | 1e-6 rel. | pass | blessed | | [nl-neo-hookean-uniaxial](../tests/oracle/nl-neo-hookean-uniaxial/) | axial Cauchy stress at the corner (1, 1, 1), stretch λ = 2 | 8.0601 | 8.0601 | < 1e-7% | 1e-6 rel. | pass | matches | | [nl-neo-hookean-uniaxial](../tests/oracle/nl-neo-hookean-uniaxial/) | lateral displacement u_y at the corner (1, 1, 1), l − 1 | -0.18686 | -0.18686 | < 1e-7% | 1e-6 rel. | pass | matches | | [nl-patch-hex8](../tests/oracle/nl-patch-hex8/) | u_x at inner node 22 | 0.13583 | 0.13583 | < 1e-7% | 1e-8 rel. | pass | matches | @@ -121,6 +125,9 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the | [patch-wedge6](../tests/oracle/patch-wedge6/) | σ_yy at the inner node (0.46, 0.55, 0.53) | 2 | 2 | < 1e-7% | 1e-9 rel. | pass | matches | | [patch-wedge6](../tests/oracle/patch-wedge6/) | σ_zz at the inner node (0.46, 0.55, 0.53) | 3 | 3 | < 1e-7% | 1e-9 rel. | pass | matches | | [patch-wedge6](../tests/oracle/patch-wedge6/) | σ_xy at the inner node | 0 | ≈ 0 | ≈ 0 | ± 1.0000e-9 | pass | matches | +| [plastic-hinge-cantilever](../tests/oracle/plastic-hinge-cantilever/) | tip force at a tip deflection of 15 (4.5 times the first yield's) | -312.5 | -322.61 | -3.23% | 5% | pass | blessed | +| [plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/) | bore displacement at the pressure that takes the plastic zone to c = 15 (a = 10, b = 20, σ_y = 250) | 0.024355 | 0.024358 | 0.0089% | 0.2% | pass | blessed | +| [plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/) | radial force on the bore (per radian times 2π, height 1) at a bore displacement of 0.3, 24 times the first yield's: the collapse pressure 2k ln 2 = 200.094 times 2πa | 12572 | 12572 | 0.00000178% | 0.1% | pass | blessed | | [pure-bending-coupling](../tests/oracle/pure-bending-coupling/) | bending stress σ_xx at mid-span, top fibre (10, 0.5, 1) | 15 | 15 | -0.0000103% | 0.5% | pass | matches | | [pure-bending-coupling](../tests/oracle/pure-bending-coupling/) | rotation of the reference point about y | 0.3 | 0.30355 | 1.18% | 1% | known failure: 1.2% high: the distributing coupling weights every node of the end face equally, corner and mid-side nodes of hex20 alike, which is not the consistent linear traction; the end face deforms locally and the fitted rotation includes it. Area-consistent weights would remove it | — | | [rotating-cylinder-axisymmetric](../tests/oracle/rotating-cylinder-axisymmetric/) | radial stress σ_r on the axis (r = 0) | 0.38571 | 0.38608 | 0.0958% | 1% | pass | matches | @@ -173,6 +180,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the - **[cdt-cylinder](../tests/oracle/cdt-cylinder/)** (theory): Constrained Delaunay mesh of a cylinder: volume within the chordal bound, valid elements. Source: mesh invariants: an inscribed polygon of sagitta s loses at most 2s/R of a circle's area (elementary geometry); constrained Delaunay tetrahedralization after Diazzi et al., ACM TOG 42(6) 2023. Tolerances — volume: the tessellation inscribes the cylinder within the chordal tolerance s = 0.05: the polygon section loses between 0 and 2s/R = 0.5% of the area; the reference is the middle of that band; invalid: every element of a valid mesh has positive volume - **[cdt-stepped-block](../tests/oracle/cdt-stepped-block/)** (theory): Constrained Delaunay mesh of a stepped block: exact volume, valid elements. Source: mesh invariants: a tetrahedralization of a polyhedron fills its volume exactly (Si, TetGen paper ACM TOMS 41(2) 2015 for the definition of a CDT; Diazzi et al., ACM TOG 42(6) 2023); the volume is 60·40·10 + 40·40·10 + 20·40·10 − 10·10·10. Tolerances — volume: a planar boundary is tessellated exactly and the constrained Delaunay mesh fills it exactly: only rounding; invalid: every element of a valid mesh has positive volume; undetermined: a linear tetrahedron's sign is exact - **[free-beam-frequencies](../tests/oracle/free-beam-frequencies/)** (theory): Natural frequencies of a free beam: rigid-body modes and the first bending mode. Source: Euler–Bernoulli free–free beam, β₁L = 4.7300; Timoshenko, Young and Weaver, Vibration Problems in Engineering, 4th ed., Wiley 1974, §5.3. Tolerances — rigid_6: a free body has six zero frequencies; numerical rigid-body frequencies below 10⁻³ of the first elastic one (eigenvalues below 10⁻⁶ of it) are rounding of the shifted solve; f7: Euler–Bernoulli neglects shear and rotary inertia, about 0.8% at L/h = 20 for the free–free mode; mass: the mass matrix integrates the density exactly +- **[j2-bauschinger](../tests/oracle/j2-bauschinger/)** (theory): Bauschinger effect of linear kinematic hardening: tension, then back to zero strain. Source: J2 plasticity with Prager's linear kinematic hardening in uniaxial stress (W. Prager, Ing.-Arch. 1955; Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3): from σ_a = σ_y + E_t (ε_a − σ_y/E) the material unloads elastically by 2σ_y and then hardens with E_t, so at zero strain σ = σ_a − 2σ_y − E_t (ε_a − 2σ_y/E). Tolerances — residual_stress: a homogeneous small-strain state the element represents exactly, with increments that end at the load path's turn; isotropic hardening would give −281.835 (it yields again only at −σ_a); plastic_strain: (σ_a − σ_y)/H forward plus E/(E + H) (ε_a − 2σ_y/E) in reverse - **[j2-uniaxial-explicit](../tests/oracle/j2-uniaxial-explicit/)** (theory): Uniaxial tension of a J2 bilinear-hardening cube to large strain (explicit). Source: von Mises plasticity with linear isotropic hardening under uniaxial stress (Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3); logarithmic strain for a hypoelastic–plastic law under coaxial stretching. Tolerances — sigma: uniaxial J2 with linear hardening: σ = σ_y + E_t (ε − σ_y/E) in logarithmic strain; the stretch takes 10⁴ wave transits, so inertia and the ringing of the start (ρcv ≈ 0.1% of σ) stay well under 1%; plastic_strain: ε_p = ε − σ/E, as for the stress - **[kirsch-plate](../tests/oracle/kirsch-plate/)** (theory): Plate with a circular hole in tension (Kirsch). Source: Kirsch, Zeitschrift des VDI 42 (1898) 797–807; Timoshenko and Goodier, Theory of Elasticity, 3rd ed., 1970, §35. Tolerances — kt: the reference is the infinite plate; the finite plate (W = 60a, Heywood's width correction about +0.1%, and the loaded ends 30a away) raises the peak by a few tenths of a percent, and the stress extrapolated to the hole surface adds up to about 1%; side: infinite-plate value −S; finite width and extrapolation as for Kt, on a smaller value - **[kirsch-plate-plane-stress](../tests/oracle/kirsch-plate-plane-stress/)** (theory): Plate with a circular hole in tension, 2D plane stress (Kirsch). Source: Kirsch, Zeitschrift des VDI 42 (1898) 797–807; Timoshenko and Goodier, Theory of Elasticity, 3rd ed., 1970, §35. Tolerances — kt: the reference is the infinite plate; the finite plate (W = 60a, loaded ends 30a away) raises the peak by a few tenths of a percent, and the stress extrapolated to the hole surface adds up to about 1%; side: infinite-plate value −S; finite width and extrapolation as for Kt, on a smaller value @@ -189,6 +197,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the - **[nafems-le10](../tests/oracle/nafems-le10/)** (benchmark): NAFEMS LE10 thick plate under pressure (hex20). Source: NAFEMS, The Standard NAFEMS Benchmarks, TNSB Rev. 3, October 1990, test LE10 (target σ_yy = −5.38 MPa at D); as reproduced in the Abaqus 2017 Benchmarks Guide, 'LE10: Thick plate under pressure'. Tolerances — syy_D: the target is published to three figures; at D the stress is extrapolated to a corner where the hole edge meets the loaded face, and published converged solutions scatter by 1–2% (Abaqus Benchmarks Guide fine meshes: −5.51 to −5.89) - **[nl-cylinder-inflation](../tests/oracle/nl-cylinder-inflation/)** (theory): Inflation of a thick neo-Hookean cylinder by a follower pressure (plane strain). Source: Finite inflation of an incompressible neo-Hookean thick-walled tube in plane strain: r² = R² + a² − A², σ_θ − σ_r = μ (λ² − λ⁻²) with λ = r/R, integrated through the wall: p = (μ/2)[(ln λ_a² − λ_a⁻²) − (ln λ_b² − λ_b⁻²)] (R. S. Rivlin, Large elastic deformations of isotropic materials VI, Phil. Trans. R. Soc. A 242:173–195, 1949; R. W. Ogden, Non-Linear Elastic Deformations, Dover 1997, §5.3.2). Tolerances — ua: the closed form is incompressible and our material nearly so (κ/μ = 1000, a volume change of order p/κ ≈ 4·10⁻⁴); hybrid hex20, 6 × 16, follow the wall's strain gradient well inside; ub: as at the bore - **[nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/)** (theory): Large-deflection cantilever under a dead tip load against the elastica (solids). Source: The elastica of a cantilever under a vertical end load: K. E. Bisshopp, D. C. Drucker, Large deflection of cantilever beams, Quarterly of Applied Mathematics 3(3):272–275, 1945 (elliptic-integral solution); values computed by integrating θ″ = −(P/EI) cos θ by shooting (our own RK4, 20000 steps), agreeing with the table of K. Mattiasson, Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals, IJNME 17(1):145–153, 1981, to all six digits. Tolerances — v1: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u1: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v2: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u2: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v5: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u5: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v10: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u10: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's +- **[nl-j2-uniaxial](../tests/oracle/nl-j2-uniaxial/)** (theory): Uniaxial tension of a J2 bilinear-hardening cube to large strain (implicit, finite strain). Source: J2 plasticity with linear isotropic hardening in uniaxial stress (Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3) in the logarithmic strain, which Simo's exponential return map (CMAME 99:61–112, 1992) integrates exactly on a proportional path: Kirchhoff stress τ = σ_y + E_t (ε − σ_y/E), Cauchy stress σ = τ/J with ln J = τ/(3K) (Hencky's volume law). Tolerances — sigma: a homogeneous deformation the element represents exactly, reached on a proportional path: the tolerance only admits the Newton tolerance and rounding; plastic_strain: ε_p = (τ − σ_y)/H, as for the stress - **[nl-neo-hookean-uniaxial](../tests/oracle/nl-neo-hookean-uniaxial/)** (theory): Uniaxial stretch of a neo-Hookean block to twice its length. Source: Closed form of the compressible neo-Hookean model W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)² in uniaxial stress (G. A. Holzapfel, Nonlinear Solid Mechanics, Wiley 2000, §6.5, eq. 6.148 with the volumetric term of §6.4): σ_yy = μ J^(−5/3)(l² − Ī/3) + κ(J − 1) = 0 fixes the lateral stretch l, then σ_xx = μ J^(−5/3)(λ² − l²). Tolerances — sigma: a homogeneous deformation every element represents exactly: the tolerance only admits rounding and the Newton tolerance; lateral: a homogeneous deformation every element represents exactly: the tolerance only admits rounding and the Newton tolerance - **[nl-patch-hex8](../tests/oracle/nl-patch-hex8/)** (theory): Nonlinear patch test, hex8: a homogeneous finite deformation of a distorted mesh. Source: The patch test (B. M. Irons, A. Razzaque, in A. K. Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, 557–587) at finite strain: a boundary moved by x = F X makes F the solution everywhere (J. Bonet, R. D. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis, 2nd ed., Cambridge 2008, §8); the neo-Hookean Cauchy stress σ = μ J^(−5/3) dev b + κ (J − 1) I of that F. Tolerances — ux: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; uy: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; uz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; sxx: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; syy: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; szz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; sxy: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; sxz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; syz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance - **[nl-rigid-rotation](../tests/oracle/nl-rigid-rotation/)** (theory): A quarter turn strains and stresses nothing (objectivity). Source: Objectivity of the total Lagrangian formulation and of the hyperelastic stress: a rigid motion x = R X + c gives C = I and σ = 0 whatever the rotation (Bonet & Wood 2008, §3.4 and §5.2; T. Belytschko, W. K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, 2nd ed., Wiley 2014, §5.5). Tolerances — ux: the inner nodes follow the rigid rotation exactly; uy: the inner nodes follow the rigid rotation exactly; uz: the inner nodes follow the rigid rotation exactly; von_mises: zero up to rounding of stresses of order E = 1000: a 90° turn leaves no strain @@ -199,6 +208,8 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the - **[patch-tet10](../tests/oracle/patch-tet10/)** (theory): Constant-stress patch test, tet10. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1) - **[patch-tet4](../tests/oracle/patch-tet4/)** (theory): Constant-stress patch test, tet4. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1) - **[patch-wedge6](../tests/oracle/patch-wedge6/)** (theory): Constant-stress patch test, wedge6. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1) +- **[plastic-hinge-cantilever](../tests/oracle/plastic-hinge-cantilever/)** (theory): Plastic collapse of a cantilever: the fully plastic moment at the root. Source: Limit analysis (B. G. Neal, The Plastic Methods of Structural Analysis, 3rd ed., Chapman & Hall 1977, ch. 1–2): beam theory's plastic hinge at the root, the fully plastic moment M_p = σ_y b h²/4 of a rectangular section, collapses a cantilever of length L at P_c = M_p/L = 312.5 (b = h = 10, L = 200, σ_y = 250); a bound for the solid, whose collapse lies above it (W. Prager, P. G. Hodge, Theory of Perfectly Plastic Solids, Wiley 1951: a continuum mechanism allows no normal velocity jump, which the beam's hinge rotation is). Tolerances — collapse_force: the solid's collapse lies above beam theory's hinge: refined through the depth (4 to 8 elements) and along the length (40 to 160) the force settles 3.1 % above M_p/L (322.3), so the band is that bound's gap; the root is a symmetry plane (no rotation, free lateral contraction) held on average through a coupling, the tip pushed through another; reduced hex20, 20 × 2 × 4 +- **[plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/)** (theory): Elastic–perfectly plastic thick cylinder: plastic zone and collapse pressure. Source: Hill, The Mathematical Theory of Plasticity, Oxford 1950, ch. V §2: a von Mises cylinder in plane strain, elastically incompressible, with the plastic zone a ≤ r ≤ c under the pressure p = k (2 ln(c/a) + 1 − c²/b²), k = σ_y/√3, and the bore displacement u_a = k c²/(2 G a); collapse when c = b at p = 2k ln(b/a). Tolerances — u_bore: 20 reduced quad8 elements through the wall resolve the logarithmic plastic field and ν = 0.4999 stands in for Hill's incompressible elasticity (0.009 % apart); u_a ∝ c², so 0.2 % on u_a is 0.1 % on the plastic radius; collapse_force: once the wall is plastic through, perfect plasticity holds the force at the collapse value (σ_θ − σ_r = 2k everywhere, σ_z the mean of the others as plastic flow dominates); the band admits the quadrature of the logarithmic field - **[pure-bending-coupling](../tests/oracle/pure-bending-coupling/)** (theory): Pure bending through a distributing coupling. Source: Euler–Bernoulli pure bending, σ = M_y z/I, θ = ML/(EI) (exact in elasticity for ν = 0); Timoshenko and Goodier, Theory of Elasticity, 3rd ed., 1970, §102. Tolerances — sigma_top: with ν = 0 pure bending is a quadratic displacement field that hex20 represents exactly; the coupling's nodal force pattern differs from the linear traction only within a depth or two of the end (Saint-Venant), five depths from the probe; rotation: the reference point follows the best-fit rigid motion of the end section, which the end zone's local deformation disturbs by well under 1% - **[rotating-cylinder-axisymmetric](../tests/oracle/rotating-cylinder-axisymmetric/)** (theory): Long rotating solid cylinder, axisymmetric. Source: Timoshenko and Goodier, Theory of Elasticity, 3rd ed., McGraw-Hill 1970, §32 (long rotating cylinder, ε_z = 0). Tolerances — radial_axis: stress extrapolated to the axis, where σ_r = σ_θ: under 1% with sixteen quadratic elements across; hoop_rim: stress extrapolated to the surface: under 1% - **[rotating-disc-plane-stress](../tests/oracle/rotating-disc-plane-stress/)** (theory): Rotating disc with a central hole, plane stress. Source: Timoshenko and Goodier, Theory of Elasticity, 3rd ed., McGraw-Hill 1970, §32 (rotating disc of uniform thickness with a hole). Tolerances — hoop_bore: stress extrapolated to the bore surface, where the a²b²/r² term is steepest: 16 radial quadratic elements graded towards the bore leave about 0.2% (8 uniform: −2.4%; 32 graded: −0.05%); radial_max: the nearest node to √(ab) is within half an element of it, where σ_r is stationary (an error of order (h/r)², under 0.5%), and nodal averaging adds a few tenths of a percent; hoop_rim: as at the bore @@ -218,6 +229,5 @@ Reference values already cited; each case gets its recipe when the capability it | [nl-cook-membrane](../tests/oracle/nl-cook-membrane/) | benchmark | Cook's membrane, nearly incompressible neo-Hookean, large load (plane strain) | the published finite-strain reference value, transcribed from the paper with its load and material constants | vertical displacement of the top corner C (48, 60) = to be transcribed | R. D. Cook, Improved two-dimensional finite element, J. Struct. Div. ASCE 100(ST9):1851–1863, 1974 (the tapered panel (0, 0), (48, 44), (48, 60), (0, 44), clamped on x = 0, shear on x = 48); finite-strain versions with a neo-Hookean material are published by J. C. Simo, F. Armero, IJNME 33:1413–1449, 1992; the values to compare with are to be transcribed from the paper | | [pinched-cylinder](../tests/oracle/pinched-cylinder/) | benchmark | Pinched cylinder with rigid diaphragms | shell elements (M8) | radial displacement under the load = -1.8248e-5 | Belytschko, Stolarski, Liu, Carpenter and Ong, 'Stress projection for membrane and shear locking in shell finite elements', Computer Methods in Applied Mechanics and Engineering 51 (1985) 221–258 (the shell obstacle course); reference 1.8248 × 10⁻⁵ under the load | | [pinched-hemisphere](../tests/oracle/pinched-hemisphere/) | benchmark | Hemispherical shell with an 18° hole under alternating point loads | shell elements (M8), or mapped blocks that follow doubly curved faces | displacement under the load along the load = 0.094 | MacNeal and Harder, 'A proposed standard set of problems to test finite element accuracy', Finite Elements in Analysis and Design 1 (1985) 3–20, Fig. 9 (p. 9), Table 5(b) (p. 11): 0.0940 (lower bound 0.0924 from their [7]) | -| [plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/) | theory | Elastic–perfectly plastic thick cylinder: collapse pressure | implicit (static) J2 plasticity: explicit J2 is in, a collapse load needs the static solver | collapse pressure for a = 10, b = 20, σ_y = 250 = 200.09 | Hill, The Mathematical Theory of Plasticity, Oxford 1950, ch. V §2: plane strain von Mises collapse pressure p = (2/√3)σ_y ln(b/a) | | [scordelis-lo-roof-shell](../tests/oracle/scordelis-lo-roof-shell/) | benchmark | Scordelis–Lo roof with shell elements | shell elements (MITC4, M8) | vertical displacement at the midpoint of the free edge = -0.3024 | MacNeal and Harder, 'A proposed standard set of problems to test finite element accuracy', Finite Elements in Analysis and Design 1 (1985) 3–20, Fig. 8 (p. 8), Table 5(a) (p. 11); Belytschko, Stolarski, Liu, Carpenter and Ong, 'Stress projection for membrane and shear locking in shell finite elements', Computer Methods in Applied Mechanics and Engineering 51 (1985) 221–258 (the shell obstacle course) | | [taylor-bar](../tests/oracle/taylor-bar/) | experiment | Taylor bar impact (copper) | the measured final shapes, to be transcribed from the paper (explicit Johnson–Cook plasticity is in) | final over initial length = to be transcribed | Wilkins and Guinan, 'Impact of cylinders on a rigid boundary', Journal of Applied Physics 44 (1973) 1200–1206 (final length and mushroom diameter vs impact speed) | diff --git a/docs/parity.md b/docs/parity.md index e90ea57..66a739d 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -44,7 +44,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | Mesh | Mesh editing (move, merge, split, collapse, remesh region), imported orphan meshes | M6 | yes | started (node/element add/remove) | | Model | Sets and surfaces bound to geometry names | M3 | yes | done | | Model | Materials: elastic (isotropic, orthotropic), density, unit-system check | M3 | | done | -| Model | Materials: plasticity, hyperelasticity, damage, thermal | M8 | | started (von Mises plasticity: bilinear, table, Johnson–Cook; Johnson–Cook damage, erosion; written to Abaqus `*PLASTIC`, Nastran `MATS1`, LS-DYNA plastic materials where expressible; compressible neo-Hookean for nonlinear steps, Abaqus `*HYPERELASTIC, NEO HOOKE`) | +| Model | Materials: plasticity, hyperelasticity, damage, thermal | M8 | | started (von Mises plasticity: bilinear, table, Johnson–Cook; Johnson–Cook damage, erosion; plasticity in explicit and nonlinear static steps; written to Abaqus `*PLASTIC`, Nastran `MATS1`, LS-DYNA plastic materials where expressible; compressible neo-Hookean for nonlinear steps, Abaqus `*HYPERELASTIC, NEO HOOKE`) | | Model | Sections: solid, shell, beam profiles; assignment by set | M3 | | done | | Model | 2D sections: plane stress, plane strain, axisymmetric (XY plane, Y the axis); element edges in sets; centrifugal load | M4 | | done | | Model | Composite layups, material orientations, continuum shells | M8 | yes | planned | @@ -53,7 +53,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | Model | Loads: force and moment, pressure, gravity, body force | M3 | | done | | Model | Reference points, kinematic and distributing coupling | M3 | | done | | Model | Steps: linear static, frequency; output presets | M3 | | done | -| Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8, M15 | | started (nonlinear static with nonlinear geometry, increments and tolerances, written to Abaqus `*STEP, NLGEOM` and LS-DYNA implicit, refused by Nastran; explicit dynamics; heat transfer; amplitudes planned) | +| Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8, M15 | | started (nonlinear static with or without nonlinear geometry, increments, tolerances and a load path (cyclic loading), written to Abaqus `*STEP, NLGEOM` and LS-DYNA implicit (the load path as LS-DYNA's curve; refused by Abaqus for now), refused by Nastran; explicit dynamics; heat transfer; amplitudes per load planned) | | Interaction | Tie between non-matching meshes | M7 | | planned | | Interaction | Contact (general and pairs), friction, automatic pair detection | M7 | | planned | | Interaction | Springs and dashpots (ground or two-point; DOF, vector or line), point masses with rotary inertia | M4 | | done (implicit; explicit translational; all three decks with what each cannot take reported) | @@ -63,7 +63,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) | | Decks | Readers keeping unknown keywords | M9 | | planned | | Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); quad4/tri3 MITC shells and Timoshenko beams with profiles, mixed with solids; reference points, kinematic and distributing couplings, moments; centrifugal loads) | -| Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff and neo-Hookean; verified by the elastica, tube inflation, patch and objectivity tests; plasticity, shells, contact, arc-length next) | +| Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff, neo-Hookean and J2 plasticity (Simo's exponential return map at finite strain, radial return at small strain, consistent tangents, the return shared with explicit dynamics); verified by the elastica, tube inflation, patch and objectivity tests, uniaxial tension to large strain, Hill's thick cylinder, a cantilever's collapse and the Bauschinger effect; Mooney–Rivlin, thermal strain, shells, contact, arc-length next) | | Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; contact planned) | | Results | Deformed shape, contours, mode shapes, probes, reaction sums | M4 | | started (results in the session, `solve.probe`, reaction sums; display planned) | | Results | XY and path plots, cuts, display groups, free-body forces, animation | M9 | | planned | diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 6e494e8..f003348 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1186,3 +1186,28 @@ update are written from (cited again in the module docs): transcribed finite-strain values), R. D. Cook, J. Struct. Div. ASCE 100(ST9):1851–1863, 1974, and J. C. Simo, F. Armero, IJNME 33:1413–1449, 1992; tube inflation, R. S. Rivlin, Phil. Trans. R. Soc. A 242:173–195, 1949, and R. W. Ogden, *Non-Linear Elastic Deformations*, Dover 1997. + +## Part N: Implicit plasticity (finite and small strain), as implemented + +- Multiplicative J2 plasticity with a Hencky elastic energy and the exponential return map in + the principal logarithmic strains, which keeps the classical radial return of the + infinitesimal theory: J. C. Simo, *Algorithms for static and dynamic multiplicative plasticity + that preserve the classical return mapping schemes of the infinitesimal theory*, CMAME + 99(1):61–112, 1992, DOI [10.1016/0045-7825(92)90123-2](https://doi.org/10.1016/0045-7825(92)90123-2); + J. C. Simo, T. J. R. Hughes, *Computational Inelasticity*, Springer 1998, ch. 9. +- The radial return (shared with explicit dynamics) and its consistent tangent + K I⊗I + 2Gθ ℙ − 2Gθ̄ n⊗n: Simo & Hughes 1998, §3.3, Box 3.2 (A4.1); J. C. Simo, R. L. Taylor, + *Consistent tangent operators for rate-independent elastoplasticity*, CMAME 48(1):101–118, + 1985, DOI [10.1016/0045-7825(85)90070-2](https://doi.org/10.1016/0045-7825(85)90070-2). +- The derivative of the matrix logarithm by divided differences (Daleckii–Krein), smooth through + equal stretches: R. Bhatia, *Matrix Analysis*, Springer 1997, §V.3. +- Kinematic hardening at finite strain (refused for now; the route when it comes): additive + plasticity in the Lagrangian logarithmic strain space, C. Miehe, N. Apel, M. Lambrecht, CMAME + 191(47–48):5383–5425, 2002, DOI [10.1016/S0045-7825(02)00438-3](https://doi.org/10.1016/S0045-7825(02)00438-3). +- Prager's linear kinematic hardening: W. Prager, *The theory of plasticity: a survey of recent + achievements*, Proc. Inst. Mech. Eng. 169:41–57, 1955. +- Verification references: the thick cylinder's plastic zone and collapse pressure, R. Hill, + *The Mathematical Theory of Plasticity*, Oxford 1950, ch. V; the cantilever's plastic hinge, + B. G. Neal, *The Plastic Methods of Structural Analysis*, 3rd ed., Chapman & Hall 1977, and the + continuum's velocity fields, W. Prager, P. G. Hodge, *Theory of Perfectly Plastic Solids*, + Wiley 1951. diff --git a/tests/oracle/j2-bauschinger/README.md b/tests/oracle/j2-bauschinger/README.md new file mode 100644 index 0000000..6d4448e --- /dev/null +++ b/tests/oracle/j2-bauschinger/README.md @@ -0,0 +1,12 @@ +# Bauschinger effect of linear kinematic hardening + +One hex8 cube (E = 200000, ν = 0.3, σ_y = 250, H = 2000, Prager's linear kinematic hardening) +held by symmetry on three faces, its fourth face pulled to a strain of 0.01 and brought back to +zero within one nonlinear static step without nonlinear geometry, by the load path +[[0, 0], [0.5, 1], [1, 0]] in increments of 0.05. + +At ε_a the stress is σ_a = σ_y + E_t (ε_a − σ_y/E) = 267.327. The yield surface keeps its size +and moves with the back stress, so unloading is elastic only until σ_a − 2σ_y and then hardens +with E_t: back at zero strain σ = σ_a − 2σ_y − E_t (ε_a − 2σ_y/E) = −247.525. Isotropic hardening +would yield again only at −σ_a and end at −281.835. The equivalent plastic strain accumulates +(σ_a − σ_y)/H forward and E/(E + H)(ε_a − 2σ_y/E) in reverse, 0.0160891. diff --git a/tests/oracle/j2-bauschinger/expected.json b/tests/oracle/j2-bauschinger/expected.json new file mode 100644 index 0000000..dae2a49 --- /dev/null +++ b/tests/oracle/j2-bauschinger/expected.json @@ -0,0 +1,45 @@ +{ + "title": "Bauschinger effect of linear kinematic hardening: tension, then back to zero strain", + "kind": "theory", + "source": "J2 plasticity with Prager's linear kinematic hardening in uniaxial stress (W. Prager, Ing.-Arch. 1955; Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3): from σ_a = σ_y + E_t (ε_a − σ_y/E) the material unloads elastically by 2σ_y and then hardens with E_t, so at zero strain σ = σ_a − 2σ_y − E_t (ε_a − 2σ_y/E)", + "quantities": [ + { + "id": "residual_stress", + "quantity": "axial stress back at zero strain after a pull to ε_a = 0.01", + "measure": { + "command": "solve.probe", + "params": { + "step": "Pull", + "elements": [ + 1 + ] + }, + "pointer": "/elements/0/stress/0/0" + }, + "value": -247.52475247524754, + "tolerance": { + "relative": 1e-06 + }, + "reason": "a homogeneous small-strain state the element represents exactly, with increments that end at the load path's turn; isotropic hardening would give −281.835 (it yields again only at −σ_a)" + }, + { + "id": "plastic_strain", + "quantity": "accumulated equivalent plastic strain", + "measure": { + "command": "solve.probe", + "params": { + "step": "Pull", + "elements": [ + 1 + ] + }, + "pointer": "/elements/0/plastic_strain/0" + }, + "value": 0.016089108910891083, + "tolerance": { + "relative": 1e-06 + }, + "reason": "(σ_a − σ_y)/H forward plus E/(E + H) (ε_a − 2σ_y/E) in reverse" + } + ] +} \ No newline at end of file diff --git a/tests/oracle/j2-bauschinger/fixture.json b/tests/oracle/j2-bauschinger/fixture.json new file mode 100644 index 0000000..fc46941 --- /dev/null +++ b/tests/oracle/j2-bauschinger/fixture.json @@ -0,0 +1,13 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "plastic_strain": { + "band": 1.6089108910891082e-11, + "value": 0.01608910891089074 + }, + "residual_stress": { + "band": 2.4752475247524754e-7, + "value": -247.52475247484173 + } + } +} diff --git a/tests/oracle/j2-bauschinger/recipe.json b/tests/oracle/j2-bauschinger/recipe.json new file mode 100644 index 0000000..1df0a23 --- /dev/null +++ b/tests/oracle/j2-bauschinger/recipe.json @@ -0,0 +1,16 @@ +[ + {"command": "mesh.block", "params": {"body": "Cube", "kind": "hex8", "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]], "divisions": [1, 1, 1], "faces": {"6": "x0", "3": "y0", "1": "z0", "4": "x1"}}}, + {"command": "model.set.create", "params": {"name": "all", "geometry": [{"body": "Cube"}]}}, + {"command": "model.set.create", "params": {"name": "x0", "geometry": [{"body": "Cube", "face": "x0"}]}}, + {"command": "model.set.create", "params": {"name": "y0", "geometry": [{"body": "Cube", "face": "y0"}]}}, + {"command": "model.set.create", "params": {"name": "z0", "geometry": [{"body": "Cube", "face": "z0"}]}}, + {"command": "model.set.create", "params": {"name": "x1", "geometry": [{"body": "Cube", "face": "x1"}]}}, + {"command": "model.material.create", "params": {"name": "Steel", "elastic": {"type": "isotropic", "young": 200000.0, "poisson": 0.3}, "plastic": {"type": "bilinear", "yield_stress": 250.0, "hardening_modulus": 2000.0, "hardening": "kinematic"}}}, + {"command": "model.section.create", "params": {"name": "cube", "set": "all", "material": "Steel", "kind": {"type": "solid"}}}, + {"command": "model.boundary_condition.create", "params": {"name": "sx", "type": "symmetry", "target": {"set": "x0"}, "plane": "x"}}, + {"command": "model.boundary_condition.create", "params": {"name": "sy", "type": "symmetry", "target": {"set": "y0"}, "plane": "y"}}, + {"command": "model.boundary_condition.create", "params": {"name": "sz", "type": "symmetry", "target": {"set": "z0"}, "plane": "z"}}, + {"command": "model.boundary_condition.create", "params": {"name": "pull", "type": "displacement", "target": {"set": "x1"}, "values": [0.01, null, null, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Pull", "procedure": {"type": "nonlinear_static", "nonlinear_geometry": false, "initial_increment": 0.05, "max_increment": 0.05, "load_path": [[0, 0], [0.5, 1], [1, 0]]}, "boundary_conditions": ["sx", "sy", "sz", "pull"], "loads": []}}, + {"command": "solve.nonlinear", "params": {"step": "Pull"}, "id": "solve"} +] diff --git a/tests/oracle/nl-j2-uniaxial/README.md b/tests/oracle/nl-j2-uniaxial/README.md new file mode 100644 index 0000000..ca90bfc --- /dev/null +++ b/tests/oracle/nl-j2-uniaxial/README.md @@ -0,0 +1,16 @@ +# Uniaxial tension of a J2 bilinear-hardening cube to large strain (implicit) + +The cube of `j2-uniaxial-explicit` solved by the implicit nonlinear solver: one fully integrated +hex8, E = 200000, ν = 0.3, σ_y = 250, H = 2000 (isotropic), held by symmetry on three faces and +stretched on the fourth to λ = e^0.2 in a nonlinear static step, lateral faces free. + +Simo's exponential return map integrates this proportional path exactly in the logarithmic +strain, so the Kirchhoff stress is the small-strain curve τ = σ_y + E_t (ε − σ_y/E) = 643.564 at +ε = 0.2 and the equivalent plastic strain (τ − σ_y)/H = 0.196782. The Cauchy stress is +σ = τ/J with Hencky's volume law ln J = τ/(3K): 642.737. + +Cross-check with the explicit solver: `j2-uniaxial-explicit` reaches 643.57 for the Cauchy stress +with the same plastic strain. Its hypoelastic update takes the logarithmic strain rate against +the Cauchy stress, while the implicit model is hyperelastic in the Kirchhoff stress: the two +laws differ by the factor J = 1.0013 (the elastic volume change), which is the 0.13 % between +them; both meet the closed form of their own law. diff --git a/tests/oracle/nl-j2-uniaxial/expected.json b/tests/oracle/nl-j2-uniaxial/expected.json new file mode 100644 index 0000000..48a1532 --- /dev/null +++ b/tests/oracle/nl-j2-uniaxial/expected.json @@ -0,0 +1,45 @@ +{ + "title": "Uniaxial tension of a J2 bilinear-hardening cube to large strain (implicit, finite strain)", + "kind": "theory", + "source": "J2 plasticity with linear isotropic hardening in uniaxial stress (Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3) in the logarithmic strain, which Simo's exponential return map (CMAME 99:61–112, 1992) integrates exactly on a proportional path: Kirchhoff stress τ = σ_y + E_t (ε − σ_y/E), Cauchy stress σ = τ/J with ln J = τ/(3K) (Hencky's volume law)", + "quantities": [ + { + "id": "sigma", + "quantity": "axial Cauchy stress at a logarithmic strain of 0.2", + "measure": { + "command": "solve.probe", + "params": { + "step": "Pull", + "elements": [ + 1 + ] + }, + "pointer": "/elements/0/stress/0/0" + }, + "value": 642.7365391418857, + "tolerance": { + "relative": 1e-06 + }, + "reason": "a homogeneous deformation the element represents exactly, reached on a proportional path: the tolerance only admits the Newton tolerance and rounding" + }, + { + "id": "plastic_strain", + "quantity": "equivalent plastic strain", + "measure": { + "command": "solve.probe", + "params": { + "step": "Pull", + "elements": [ + 1 + ] + }, + "pointer": "/elements/0/plastic_strain/0" + }, + "value": 0.19678217821782182, + "tolerance": { + "relative": 1e-06 + }, + "reason": "ε_p = (τ − σ_y)/H, as for the stress" + } + ] +} \ No newline at end of file diff --git a/tests/oracle/nl-j2-uniaxial/fixture.json b/tests/oracle/nl-j2-uniaxial/fixture.json new file mode 100644 index 0000000..08df0ad --- /dev/null +++ b/tests/oracle/nl-j2-uniaxial/fixture.json @@ -0,0 +1,13 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "plastic_strain": { + "band": 1.9678217821795454e-10, + "value": 0.19678217821795452 + }, + "sigma": { + "band": 6.427365391418857e-7, + "value": 642.7365390754569 + } + } +} diff --git a/tests/oracle/nl-j2-uniaxial/recipe.json b/tests/oracle/nl-j2-uniaxial/recipe.json new file mode 100644 index 0000000..2316eb4 --- /dev/null +++ b/tests/oracle/nl-j2-uniaxial/recipe.json @@ -0,0 +1,16 @@ +[ + {"command": "mesh.block", "params": {"body": "Cube", "kind": "hex8", "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]], "divisions": [1, 1, 1], "faces": {"6": "x0", "3": "y0", "1": "z0", "4": "x1"}}}, + {"command": "model.set.create", "params": {"name": "all", "geometry": [{"body": "Cube"}]}}, + {"command": "model.set.create", "params": {"name": "x0", "geometry": [{"body": "Cube", "face": "x0"}]}}, + {"command": "model.set.create", "params": {"name": "y0", "geometry": [{"body": "Cube", "face": "y0"}]}}, + {"command": "model.set.create", "params": {"name": "z0", "geometry": [{"body": "Cube", "face": "z0"}]}}, + {"command": "model.set.create", "params": {"name": "x1", "geometry": [{"body": "Cube", "face": "x1"}]}}, + {"command": "model.material.create", "params": {"name": "Steel", "elastic": {"type": "isotropic", "young": 200000.0, "poisson": 0.3}, "plastic": {"type": "bilinear", "yield_stress": 250.0, "hardening_modulus": 2000.0, "hardening": "isotropic"}}}, + {"command": "model.section.create", "params": {"name": "cube", "set": "all", "material": "Steel", "kind": {"type": "solid"}}}, + {"command": "model.boundary_condition.create", "params": {"name": "sx", "type": "symmetry", "target": {"set": "x0"}, "plane": "x"}}, + {"command": "model.boundary_condition.create", "params": {"name": "sy", "type": "symmetry", "target": {"set": "y0"}, "plane": "y"}}, + {"command": "model.boundary_condition.create", "params": {"name": "sz", "type": "symmetry", "target": {"set": "z0"}, "plane": "z"}}, + {"command": "model.boundary_condition.create", "params": {"name": "pull", "type": "displacement", "target": {"set": "x1"}, "values": [0.22140275816016985, null, null, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Pull", "procedure": {"type": "nonlinear_static"}, "boundary_conditions": ["sx", "sy", "sz", "pull"], "loads": []}}, + {"command": "solve.nonlinear", "params": {"step": "Pull"}, "id": "solve"} +] diff --git a/tests/oracle/plastic-hinge-cantilever/README.md b/tests/oracle/plastic-hinge-cantilever/README.md new file mode 100644 index 0000000..02b9799 --- /dev/null +++ b/tests/oracle/plastic-hinge-cantilever/README.md @@ -0,0 +1,13 @@ +# Plastic collapse of a cantilever + +A steel bar 200 long with a 10 × 10 section (E = 200000, ν = 0.3, elastic–perfectly plastic at +σ_y = 250), reduced hex20 20 × 2 × 4, in a nonlinear static step without nonlinear geometry. The +root x = 0 is a plane of symmetry (no rotation, the section free to contract laterally), held in +Y, Z and against twist on average through a distributing coupling; a coupling at the tip pushes +the end down by 15, 4.5 times the deflection at first yield. + +Beam theory's collapse is a plastic hinge at the root: M_p = σ_y b h²/4 = 62500 and +P_c = M_p/L = 312.5. The solid's collapse lies above it, since a continuum cannot rotate about a +hinge (that needs a jump of the normal velocity across the section): refined to 8 elements +through the depth or 160 along the length, the force settles at 322.3, 3.1 % above. The case +checks that the implicit plasticity reaches this plateau and holds it. diff --git a/tests/oracle/plastic-hinge-cantilever/expected.json b/tests/oracle/plastic-hinge-cantilever/expected.json new file mode 100644 index 0000000..64a7b36 --- /dev/null +++ b/tests/oracle/plastic-hinge-cantilever/expected.json @@ -0,0 +1,20 @@ +{ + "title": "Plastic collapse of a cantilever: the fully plastic moment at the root", + "kind": "theory", + "source": "Limit analysis (B. G. Neal, The Plastic Methods of Structural Analysis, 3rd ed., Chapman & Hall 1977, ch. 1–2): beam theory's plastic hinge at the root, the fully plastic moment M_p = σ_y b h²/4 of a rectangular section, collapses a cantilever of length L at P_c = M_p/L = 312.5 (b = h = 10, L = 200, σ_y = 250); a bound for the solid, whose collapse lies above it (W. Prager, P. G. Hodge, Theory of Perfectly Plastic Solids, Wiley 1951: a continuum mechanism allows no normal velocity jump, which the beam's hinge rotation is)", + "quantities": [ + { + "id": "collapse_force", + "quantity": "tip force at a tip deflection of 15 (4.5 times the first yield's)", + "measure": { + "from": "solve", + "pointer": "/reference_points/1/reaction/2" + }, + "value": -312.5, + "tolerance": { + "relative": 0.05 + }, + "reason": "the solid's collapse lies above beam theory's hinge: refined through the depth (4 to 8 elements) and along the length (40 to 160) the force settles 3.1 % above M_p/L (322.3), so the band is that bound's gap; the root is a symmetry plane (no rotation, free lateral contraction) held on average through a coupling, the tip pushed through another; reduced hex20, 20 × 2 × 4" + } + ] +} \ No newline at end of file diff --git a/tests/oracle/plastic-hinge-cantilever/fixture.json b/tests/oracle/plastic-hinge-cantilever/fixture.json new file mode 100644 index 0000000..1d752f6 --- /dev/null +++ b/tests/oracle/plastic-hinge-cantilever/fixture.json @@ -0,0 +1,9 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "collapse_force": { + "band": 3.2260685847425373e-7, + "value": -322.6068584742537 + } + } +} diff --git a/tests/oracle/plastic-hinge-cantilever/recipe.json b/tests/oracle/plastic-hinge-cantilever/recipe.json new file mode 100644 index 0000000..e0374d6 --- /dev/null +++ b/tests/oracle/plastic-hinge-cantilever/recipe.json @@ -0,0 +1,17 @@ +[ + {"command": "mesh.block", "params": {"body": "Beam", "kind": "hex20", "corners": [[0, 0, -5.0], [200.0, 0, -5.0], [200.0, 10.0, -5.0], [0, 10.0, -5.0], [0, 0, 5.0], [200.0, 0, 5.0], [200.0, 10.0, 5.0], [0, 10.0, 5.0]], "divisions": [20, 2, 4], "faces": {"6": "root", "4": "tip"}}}, + {"command": "model.set.create", "params": {"name": "all", "geometry": [{"body": "Beam"}]}}, + {"command": "model.set.create", "params": {"name": "root", "geometry": [{"body": "Beam", "face": "root"}]}}, + {"command": "model.set.create", "params": {"name": "tip", "geometry": [{"body": "Beam", "face": "tip"}]}}, + {"command": "model.material.create", "params": {"name": "Steel", "elastic": {"type": "isotropic", "young": 200000.0, "poisson": 0.3}, "plastic": {"type": "bilinear", "yield_stress": 250.0}}}, + {"command": "model.section.create", "params": {"name": "beam", "set": "all", "material": "Steel", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "model.boundary_condition.create", "params": {"name": "plane", "type": "symmetry", "target": {"set": "root"}, "plane": "x"}}, + {"command": "model.reference_point.create", "params": {"name": "Root", "xyz": [0, 5.0, 0]}}, + {"command": "model.coupling.create", "params": {"name": "support", "reference_point": "Root", "set": "root", "kind": "distributing"}}, + {"command": "model.boundary_condition.create", "params": {"name": "hold", "type": "displacement", "target": {"reference_point": "Root"}, "values": [null, 0, 0, 0, null, null]}}, + {"command": "model.reference_point.create", "params": {"name": "Tip", "xyz": [200.0, 5.0, 0]}}, + {"command": "model.coupling.create", "params": {"name": "end", "reference_point": "Tip", "set": "tip", "kind": "distributing"}}, + {"command": "model.boundary_condition.create", "params": {"name": "push", "type": "displacement", "target": {"reference_point": "Tip"}, "values": [null, null, -15.0, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Collapse", "procedure": {"type": "nonlinear_static", "nonlinear_geometry": false}, "boundary_conditions": ["plane", "hold", "push"], "loads": []}}, + {"command": "solve.nonlinear", "params": {"step": "Collapse"}, "id": "solve"} +] diff --git a/tests/oracle/plastic-thick-cylinder/README.md b/tests/oracle/plastic-thick-cylinder/README.md index daa8eeb..8962c28 100644 --- a/tests/oracle/plastic-thick-cylinder/README.md +++ b/tests/oracle/plastic-thick-cylinder/README.md @@ -1,3 +1,14 @@ -# Elastic–perfectly plastic thick cylinder: collapse pressure +# Elastic–perfectly plastic thick cylinder: plastic zone and collapse pressure -The Lamé cylinder of `lame-cylinder`, elastic–perfectly plastic, loaded to collapse. +The Lamé cylinder of `lame-cylinder-axisymmetric` (a = 10, b = 20, an axisymmetric slice held +axially on both faces, so plane strain), now elastic–perfectly plastic: E = 200000, +ν = 0.4999 (Hill's solution takes the elasticity incompressible), σ_y = 250, one row of 20 +reduced quad8 elements through the wall, in two nonlinear static steps without nonlinear +geometry: + +- **Pressure**: the bore pressure p = k (2 ln(c/a) + 1 − c²/b²) = 180.195 (k = σ_y/√3) that + takes the plastic zone to c = 15. The bore then moves u_a = k c²/(2 G a) = 0.024355: the plastic + zone moves incompressibly (u r constant) and the elastic ring outside yields at its inner edge. +- **Collapse**: the bore pushed out by 0.3 (24 times its first-yield displacement), past the + pressure at which the wall is plastic through; the radial force on the bore stays at the + collapse pressure 2k ln(b/a) = 200.094 times 2πa per unit height. diff --git a/tests/oracle/plastic-thick-cylinder/expected.json b/tests/oracle/plastic-thick-cylinder/expected.json index 6d74b98..cbbff6d 100644 --- a/tests/oracle/plastic-thick-cylinder/expected.json +++ b/tests/oracle/plastic-thick-cylinder/expected.json @@ -1,14 +1,43 @@ { - "title": "Elastic–perfectly plastic thick cylinder: collapse pressure", + "title": "Elastic–perfectly plastic thick cylinder: plastic zone and collapse pressure", "kind": "theory", - "source": "Hill, The Mathematical Theory of Plasticity, Oxford 1950, ch. V §2: plane strain von Mises collapse pressure p = (2/√3)σ_y ln(b/a)", - "status": "pending", - "requires": "implicit (static) J2 plasticity: explicit J2 is in, a collapse load needs the static solver", + "source": "Hill, The Mathematical Theory of Plasticity, Oxford 1950, ch. V §2: a von Mises cylinder in plane strain, elastically incompressible, with the plastic zone a ≤ r ≤ c under the pressure p = k (2 ln(c/a) + 1 − c²/b²), k = σ_y/√3, and the bore displacement u_a = k c²/(2 G a); collapse when c = b at p = 2k ln(b/a)", "quantities": [ { - "id": "p_collapse", - "quantity": "collapse pressure for a = 10, b = 20, σ_y = 250", - "value": 200.0944 + "id": "u_bore", + "quantity": "bore displacement at the pressure that takes the plastic zone to c = 15 (a = 10, b = 20, σ_y = 250)", + "measure": { + "command": "solve.probe", + "params": { + "step": "Pressure", + "at": [ + [ + 10.0, + 0, + 0 + ] + ] + }, + "pointer": "/nodes/0/displacement/0" + }, + "value": 0.02435534068380524, + "tolerance": { + "relative": 0.002 + }, + "reason": "20 reduced quad8 elements through the wall resolve the logarithmic plastic field and ν = 0.4999 stands in for Hill's incompressible elasticity (0.009 % apart); u_a ∝ c², so 0.2 % on u_a is 0.1 % on the plastic radius" + }, + { + "id": "collapse_force", + "quantity": "radial force on the bore (per radian times 2π, height 1) at a bore displacement of 0.3, 24 times the first yield's: the collapse pressure 2k ln 2 = 200.094 times 2πa", + "measure": { + "from": "collapse", + "pointer": "/reaction_total/0" + }, + "value": 12572.299154203694, + "tolerance": { + "relative": 0.001 + }, + "reason": "once the wall is plastic through, perfect plasticity holds the force at the collapse value (σ_θ − σ_r = 2k everywhere, σ_z the mean of the others as plastic flow dominates); the band admits the quadrature of the logarithmic field" } ] -} +} \ No newline at end of file diff --git a/tests/oracle/plastic-thick-cylinder/fixture.json b/tests/oracle/plastic-thick-cylinder/fixture.json new file mode 100644 index 0000000..26a446c --- /dev/null +++ b/tests/oracle/plastic-thick-cylinder/fixture.json @@ -0,0 +1,13 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "collapse_force": { + "band": 0.000012572299377396248, + "value": 12572.299377396248 + }, + "u_bore": { + "band": 2.4357507851220506e-11, + "value": 0.024357507851220504 + } + } +} diff --git a/tests/oracle/plastic-thick-cylinder/recipe.json b/tests/oracle/plastic-thick-cylinder/recipe.json new file mode 100644 index 0000000..9059e22 --- /dev/null +++ b/tests/oracle/plastic-thick-cylinder/recipe.json @@ -0,0 +1,15 @@ +[ + {"command": "mesh.block", "params": {"body": "Pipe", "kind": "quad8", "corners": [[10.0, 0, 0], [20.0, 0, 0], [20.0, 1, 0], [10.0, 1, 0]], "divisions": [20, 1], "bias": [1, 1], "edges": [{"corners": [1, 4], "name": "bore"}, {"corners": [1, 2], "name": "bottom"}, {"corners": [3, 4], "name": "top"}]}}, + {"command": "model.set.create", "params": {"name": "all", "geometry": [{"body": "Pipe", "face": "area"}]}}, + {"command": "model.set.create", "params": {"name": "bore", "geometry": [{"body": "Pipe", "edge": "bore"}]}}, + {"command": "model.set.create", "params": {"name": "ends", "geometry": [{"body": "Pipe", "edge": "bottom"}, {"body": "Pipe", "edge": "top"}]}}, + {"command": "model.material.create", "params": {"name": "Steel", "elastic": {"type": "isotropic", "young": 200000.0, "poisson": 0.4999}, "plastic": {"type": "bilinear", "yield_stress": 250.0}}}, + {"command": "model.section.create", "params": {"name": "pipe", "set": "all", "material": "Steel", "kind": {"type": "axisymmetric"}, "formulation": {"integration": "reduced"}}}, + {"command": "model.boundary_condition.create", "params": {"name": "ends", "type": "symmetry", "target": {"set": "ends"}, "plane": "y"}}, + {"command": "model.load.create", "params": {"name": "p", "type": "pressure", "set": "bore", "magnitude": 180.19538034924005}}, + {"command": "model.boundary_condition.create", "params": {"name": "expand", "type": "displacement", "target": {"set": "bore"}, "values": [0.3, null, null, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Pressure", "procedure": {"type": "nonlinear_static", "nonlinear_geometry": false}, "boundary_conditions": ["ends"], "loads": ["p"]}}, + {"command": "model.step.create", "params": {"name": "Collapse", "procedure": {"type": "nonlinear_static", "nonlinear_geometry": false}, "boundary_conditions": ["ends", "expand"], "loads": []}}, + {"command": "solve.nonlinear", "params": {"step": "Pressure"}, "id": "pressure"}, + {"command": "solve.nonlinear", "params": {"step": "Collapse"}, "id": "collapse"} +]