overset added-mass pin (P5-3): a no-slip cylinder on the moving patch oscillating in a closed box of still fluid, the force fitted over the last two periods — C_m = 1.069 (n = 32) / 1.097 (n = 64) against Stokes's 1.064 plus box blockage: the composite's unsteady pressure response is right, the FSI2 lift excess is not the unsteady term's scale

Co-Authored-By: Claude Fable 5.1 <[email protected]> Claude-Session: https://claude.ai/code/session_0116sg1Qz1gMv9hdcKP1XUam
2026-09-09 18:01:53 -07:00
parent aeb70eb507
commit 33534ee688
1 changed files with 206 additions and 0 deletions
@@ -0,0 +1,206 @@
+//! P5-3's added-mass pin (`docs/overset_metal_campaign.md`, the replay-to-
+//! replay reading): the composite's UNSTEADY pressure response on the
+//! moving patch against potential theory.
+//!
+//! A no-slip circular cylinder of radius `r` on an O-grid patch translates
+//! through a closed box of still fluid, `x_c(t) = A sin ωt`, the patch
+//! carried with it every step (`set_patch_mesh`, the A-P2 moving path)
+//! and the wall carrying the exact velocity. The in-phase force is the
+//! added-mass reaction, `F = −m_a ẍ` with `m_a = ρ π r²` (unbounded
+//! potential flow), plus Stokes's viscous correction `4 / √(π β)` with
+//! `β = r² ω / ν` and a small blockage term for the box. The FSI2
+//! ladder's motion-fixed lift response is ~1.4× the reference's at every
+//! h and dt (2026-09-09); nothing rigid (CFD1–3, all within 1 %) can see
+//! the unsteady scale, this can. Registered before the run: `C_m` in
+//! [1.0, 1.25] at n = 64 means the unsteady scale is right and the
+//! discrepancy is the wake's; `C_m ≥ 1.4` locates it here.
+//!
+//! `RTX_OVERSET_N` (64), `RTX_OVERSET_AM_PERIODS` (4).
+
+use rtx_cfd::mesh::patch_gen::annulus_skewed;
+use rtx_cfd::mesh::{PatchMesh, PatchSide};
+use rtx_cfd::solvers::incompressible::{
+    CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters, EmbeddedPisoSolver,
+    FlowField, NormalDiffusion, OversetField, OversetParameters, OversetPisoSolver, PatchField,
+    PoissonSolverKind,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+use std::f64::consts::PI;
+
+const RHO: f64 = 1.0;
+const NU: f64 = 5e-5;
+const R: f64 = 0.1;
+const R_OUT: f64 = 0.25;
+const CX: f64 = 0.5;
+const CY: f64 = 0.5;
+const AMP: f64 = 0.005;
+const OMEGA: f64 = 2.0 * PI;
+
+fn patch_at(n: usize, x: f64) -> CfdResult<PatchMesh> {
+    annulus_skewed([CX + x, CY], R, R_OUT, 9 * n / 4, n / 4, 0.0, 3.0)
+}
+
+fn env_usize(k: &str, d: usize) -> usize {
+    std::env::var(k)
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(d)
+}
+
+#[tokio::test]
+async fn oscillating_cylinder_added_mass_on_the_moving_patch() -> CfdResult<()> {
+    let n = env_usize("RTX_OVERSET_N", 64);
+    let periods = env_usize("RTX_OVERSET_AM_PERIODS", 4);
+    let h = 1.0 / n as f64;
+    let config = CfdConfig::new()
+        .with_density(RHO)
+        .with_viscosity(RHO * NU)
+        .with_reference_velocity(AMP * OMEGA)
+        .with_reference_length(2.0 * R);
+    let mut background = EmbeddedPisoSolver::new(
+        config.clone(),
+        EmbeddedParameters {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            poisson_solver: PoissonSolverKind::Multigrid,
+            ..EmbeddedParameters::default()
+        },
+    )?;
+    background.set_boundary_velocity(|_, _, _| (0.0, 0.0));
+    let mut patch = CurvilinearPisoSolver::new(
+        config,
+        CurvilinearParameters {
+            tolerance: 1e-6,
+            normal_diffusion: NormalDiffusion::LineImplicit,
+            ..CurvilinearParameters::default()
+        },
+        patch_at(n, 0.0)?,
+    )?;
+    // The wall's exact velocity: the patch's clock is the end of the step.
+    patch.set_side_velocity(PatchSide::Inner, |_, _, t| {
+        (AMP * OMEGA * (OMEGA * t).cos(), 0.0)
+    });
+    let mut patch_field = PatchField::new(patch.mesh());
+    patch.initialize(&mut patch_field, |_, _| (0.0, 0.0));
+    let bg_field = FlowField::new(n, n, h, h)?;
+    let params = OversetParameters {
+        stall_rounds: 0,
+        ..OversetParameters::default()
+    };
+    let mut solver = OversetPisoSolver::new(background, patch, (n, n, h, h), params)?;
+    let mut field = OversetField {
+        background: bg_field,
+        patch: patch_field,
+    };
+    solver.initialize(&mut field)?;
+
+    // Step: the patch's smallest edge against the explicit limits, and the
+    // wall's own displacement per step well under a cell.
+    let mut hp = f64::INFINITY;
+    for c in 0..solver.patch().mesh().cell_count() {
+        for (f, _) in solver.patch().mesh().cell_faces(c) {
+            let d = solver.patch().mesh().faces()[f].d;
+            hp = hp.min((d[0] * d[0] + d[1] * d[1]).sqrt());
+        }
+    }
+    let period = 2.0 * PI / OMEGA;
+    let dt_raw = 0.4 * (hp * hp / (4.0 * NU)).min(h).min(0.2 * hp / (AMP * OMEGA));
+    let steps_per_period = (period / dt_raw).ceil() as usize;
+    let dt = period / steps_per_period as f64;
+    let beta = R * R * OMEGA / NU;
+    println!(
+        "  n = {n}, h = {h:.4}, patch inner edge {hp:.4}, dt = {dt:.3e} ({steps_per_period} per period), A/r = {:.3}, KC = {:.3}, β = {beta:.0}: Stokes C_m ≈ {:.3}",
+        AMP / R,
+        2.0 * PI * AMP / R,
+        1.0 + 4.0 / (PI * beta).sqrt()
+    );
+
+    let mut samples: Vec<(f64, f64, f64, f64)> = Vec::new();
+    let mut rounds_total = 0usize;
+    let mut reclass_total = 0usize;
+    let start = std::time::Instant::now();
+    for step in 0..periods * steps_per_period {
+        let t_new = (step + 1) as f64 * dt;
+        solver.set_patch_mesh(patch_at(n, AMP * (OMEGA * t_new).sin())?)?;
+        let r = solver.advance(&mut field, dt).await?;
+        rounds_total += r.rounds.iter().sum::<usize>();
+        reclass_total += r.reclassified_cells;
+        let load = solver
+            .patch()
+            .surface_force(&field.patch, PatchSide::Inner, solver.time());
+        let f = load.total();
+        samples.push((t_new, f[0], load.pressure[0], load.viscous[0]));
+        if (step + 1) % steps_per_period == 0 {
+            println!(
+                "    period {}: F_x at the end {:+.4e} (pressure {:+.4e}, viscous {:+.4e}), rounds {:.2}/step, reclassified {:.1}/step, {:.0} s",
+                (step + 1) / steps_per_period,
+                f[0],
+                load.pressure[0],
+                load.viscous[0],
+                rounds_total as f64 / (step + 1) as f64,
+                reclass_total as f64 / (step + 1) as f64,
+                start.elapsed().as_secs_f64()
+            );
+        }
+    }
+    // Least squares over the last two periods: F = a sin ωt + b cos ωt (+ c).
+    // ẍ = −A ω² sin ωt, so the added-mass reaction −m_a ẍ = m_a A ω² sin ωt:
+    // C_m = a / (ρ π r² A ω²). ẋ = A ω cos ωt, so b = −c_d A ω (damping).
+    let fit = |col: usize| -> (f64, f64, f64) {
+        let last: Vec<&(f64, f64, f64, f64)> = samples
+            .iter()
+            .filter(|s| s.0 > (periods as f64 - 2.0) * period - 1e-12)
+            .collect();
+        let g = |s: &(f64, f64, f64, f64)| match col {
+            1 => s.1,
+            2 => s.2,
+            _ => s.3,
+        };
+        // Normal equations for [sin, cos, 1].
+        let mut m = [[0.0f64; 3]; 3];
+        let mut rhs = [0.0f64; 3];
+        for s in &last {
+            let b = [(OMEGA * s.0).sin(), (OMEGA * s.0).cos(), 1.0];
+            for i in 0..3 {
+                rhs[i] += b[i] * g(s);
+                for j in 0..3 {
+                    m[i][j] += b[i] * b[j];
+                }
+            }
+        }
+        // 3×3 solve by Cramer's rule.
+        let det = |a: [[f64; 3]; 3]| {
+            a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
+                - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0])
+                + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])
+        };
+        let d = det(m);
+        let mut sol = [0.0; 3];
+        for k in 0..3 {
+            let mut mk = m;
+            for i in 0..3 {
+                mk[i][k] = rhs[i];
+            }
+            sol[k] = det(mk) / d;
+        }
+        (sol[0], sol[1], sol[2])
+    };
+    let (a, b, c) = fit(1);
+    let (ap, bp, _) = fit(2);
+    let (av, bv, _) = fit(3);
+    let m_a = RHO * PI * R * R;
+    let cm = a / (m_a * AMP * OMEGA * OMEGA);
+    let cm_p = ap / (m_a * AMP * OMEGA * OMEGA);
+    let cm_v = av / (m_a * AMP * OMEGA * OMEGA);
+    let damping = -b / (AMP * OMEGA);
+    println!(
+        "  ADDED MASS n = {n}: C_m = {cm:.4} (pressure {cm_p:.4} + viscous {cm_v:.4}) against 1 + Stokes {:.3}; damping coefficient {damping:.4e} (pressure {:+.3e}, viscous {:+.3e}) [N·s/m per m]; mean force {c:+.3e}; {} steps in {:.0} s",
+        1.0 + 4.0 / (PI * beta).sqrt(),
+        -bp / (AMP * OMEGA),
+        -bv / (AMP * OMEGA),
+        samples.len(),
+        start.elapsed().as_secs_f64()
+    );
+    assert!(cm.is_finite(), "non-finite added-mass coefficient");
+    Ok(())
+}