embedded3 S2-3c: the 2D solver on the flag wake's kinematics (the reference for the 3D full-span run)

Co-Authored-By: Claude Fable 5.1 <[email protected]>
2026-09-17 23:14:42 -05:00
parent 6d2872e7d0
commit f2cdda4691
1 changed files with 240 additions and 0 deletions
@@ -0,0 +1,240 @@
+//! S2-3c: the 2D solver on the flag wake's exact kinematics — the same
+//! cylinder + capsule flag, the same first-mode motion (tip 84 mm at
+//! 1.930 Hz), the same inflow mean (Ū 1.0, the 2D parabola), TVD, on the
+//! ny 62 spacing (378 × 62, h 6.6 mm). Its loads are the 2D answer for
+//! this kinematics; the 3D full-span run must reproduce them. Loads by
+//! the traction sampler over the flag's and the cylinder's surface
+//! samples per unit span, every 10 steps; the last period's mean drag and
+//! median-filtered lift swing.
+//!
+//! `cargo test --release -p rtx-cfd --test embedded3_flag_reference_2d -- --ignored --nocapture`
+use rtx_cfd::solvers::incompressible::{
+    AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
+    FlowField, MgPrecision, MgSmoother, PoissonSolverKind, SideBoundary,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+
+const H: f64 = 0.41;
+const L: f64 = 2.5;
+const CX: f64 = 0.2;
+const CY: f64 = 0.2;
+const R_CYL: f64 = 0.05;
+const FLAG_X0: f64 = 0.6;
+const FLAG_LEN: f64 = 0.35;
+const FLAG_HALF: f64 = 0.01;
+const AMP: f64 = 0.084;
+const FREQ: f64 = 1.930;
+const U_MEAN: f64 = 1.0;
+const RHO: f64 = 1000.0;
+const NU: f64 = 1e-3;
+const BETA_L: f64 = 1.875_104_069;
+
+fn mode(s: f64) -> f64 {
+    let b = BETA_L;
+    let sigma = (b.sinh() - b.sin()) / (b.cosh() + b.cos());
+    let phi = |s: f64| (b * s).cosh() - (b * s).cos() - sigma * ((b * s).sinh() - (b * s).sin());
+    phi(s) / phi(1.0)
+}
+
+fn deflection(s: f64, t: f64) -> (f64, f64) {
+    let w = 2.0 * std::f64::consts::PI * FREQ;
+    (
+        AMP * mode(s) * (w * t).sin(),
+        AMP * mode(s) * w * (w * t).cos(),
+    )
+}
+
+fn centreline(m: usize, n: usize, t: f64) -> (f64, f64, f64) {
+    let s = m as f64 / n as f64;
+    let (d, v) = deflection(s, t);
+    (FLAG_X0 + s * FLAG_LEN, CY + d, v)
+}
+
+/// Distance to the capsule flag and the centreline velocity at the foot.
+fn flag_sdf(x: f64, y: f64, t: f64) -> (f64, f64) {
+    let n = 40;
+    let mut best = f64::INFINITY;
+    let mut v_best = 0.0;
+    for m in 0..n {
+        let (ax, ay, av) = centreline(m, n, t);
+        let (bx, by, bv) = centreline(m + 1, n, t);
+        let (ex, ey) = (bx - ax, by - ay);
+        let u = (((x - ax) * ex + (y - ay) * ey) / (ex * ex + ey * ey)).clamp(0.0, 1.0);
+        let d = ((x - ax - u * ex).powi(2) + (y - ay - u * ey).powi(2)).sqrt();
+        if d < best {
+            best = d;
+            v_best = av + u * (bv - av);
+        }
+    }
+    (best - FLAG_HALF, v_best)
+}
+
+fn cyl_sdf(x: f64, y: f64) -> f64 {
+    ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL
+}
+
+/// Surface samples `(x, y, nx, ny, ds)` at `t`: the flag's two sides and
+/// tip along the deflected centreline, the cylinder's circle; samples
+/// inside the other body are dropped.
+fn samples(t: f64, ds: f64) -> Vec<(f64, f64, f64, f64, f64)> {
+    let mut out = Vec::new();
+    let n = ((FLAG_LEN / ds).ceil() as usize).max(8);
+    for m in 0..n {
+        let (ax, ay, _) = centreline(m, n, t);
+        let (bx, by, _) = centreline(m + 1, n, t);
+        let (ex, ey) = (bx - ax, by - ay);
+        let len = (ex * ex + ey * ey).sqrt();
+        let (tx, ty) = (ex / len, ey / len);
+        let (nx, ny) = (-ty, tx);
+        let (mx, my) = (0.5 * (ax + bx), 0.5 * (ay + by));
+        for sign in [1.0, -1.0] {
+            let (px, py) = (mx + sign * FLAG_HALF * nx, my + sign * FLAG_HALF * ny);
+            if cyl_sdf(px, py) > 0.0 {
+                out.push((px, py, sign * nx, sign * ny, len));
+            }
+        }
+    }
+    // The tip: a semicircle around the last centreline point.
+    let (tx0, ty0, _) = centreline(n, n, t);
+    let (px, py, _) = centreline(n - 1, n, t);
+    let ang0 = (ty0 - py).atan2(tx0 - px);
+    let n_arc = ((std::f64::consts::PI * FLAG_HALF / ds).ceil() as usize).max(4);
+    for k in 0..n_arc {
+        let a = ang0 - std::f64::consts::FRAC_PI_2
+            + (k as f64 + 0.5) / n_arc as f64 * std::f64::consts::PI;
+        out.push((
+            tx0 + FLAG_HALF * a.cos(),
+            ty0 + FLAG_HALF * a.sin(),
+            a.cos(),
+            a.sin(),
+            std::f64::consts::PI * FLAG_HALF / n_arc as f64,
+        ));
+    }
+    let n_c = ((2.0 * std::f64::consts::PI * R_CYL / ds).ceil() as usize).max(16);
+    for k in 0..n_c {
+        let a = (k as f64 + 0.5) / n_c as f64 * 2.0 * std::f64::consts::PI;
+        let (px, py) = (CX + R_CYL * a.cos(), CY + R_CYL * a.sin());
+        if flag_sdf(px, py, t).0 > 0.0 {
+            out.push((
+                px,
+                py,
+                a.cos(),
+                a.sin(),
+                2.0 * std::f64::consts::PI * R_CYL / n_c as f64,
+            ));
+        }
+    }
+    out
+}
+
+fn median(v: &[f64]) -> f64 {
+    let mut s = v.to_vec();
+    s.sort_by(|a, b| a.partial_cmp(b).unwrap());
+    s[s.len() / 2]
+}
+
+#[tokio::test]
+#[ignore = "S2-3c: the 2D reference for the flag wake's kinematics (minutes on the host)"]
+async fn flag_wake_2d_reference() -> CfdResult<()> {
+    let ny = 62;
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let dt = 8.817e-4;
+    let period = 1.0 / FREQ;
+    let t_end = 2.0 * period;
+    let config = CfdConfig::new()
+        .with_density(RHO)
+        .with_viscosity(RHO * NU)
+        .with_reference_velocity(U_MEAN)
+        .with_reference_length(2.0 * R_CYL);
+    let mut solver = EmbeddedPisoSolver::new(
+        config,
+        EmbeddedParameters {
+            corrector_steps: 3,
+            tolerance: 1e-8,
+            boundaries: AleBoundaries {
+                right: SideBoundary::PressureOutlet,
+                ..AleBoundaries::default()
+            },
+            poisson_solver: PoissonSolverKind::Multigrid,
+            poisson_precision: MgPrecision::F64,
+            poisson_smoother: MgSmoother::Lexicographic,
+            convection_scheme: ConvectionScheme::TvdVanAlbada,
+        },
+    )?;
+    solver.set_boundary_velocity(|x, y, _t| {
+        if x <= 0.0 {
+            (1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2), 0.0)
+        } else {
+            (0.0, 0.0)
+        }
+    });
+    let flag = EmbeddedBody::from_sdf(|x, y, t| flag_sdf(x, y, t).0)
+        .with_surface_velocity(|x, y, t| (0.0, flag_sdf(x, y, t).1));
+    solver.set_moving_body(EmbeddedBody::union(
+        EmbeddedBody::circle(CX, CY, R_CYL),
+        flag,
+    ));
+    let mut field = FlowField::new(nx, ny, h, h)?;
+    for j in 0..ny {
+        let u0 =
+            1.5 * U_MEAN * ((j as f64 + 0.5) * h) * (H - (j as f64 + 0.5) * h) / (0.5 * H).powi(2);
+        for i in 0..=nx {
+            field.u[(j, i)] = u0;
+        }
+    }
+    solver.initialize(&mut field)?;
+    let steps = (t_end / dt).ceil() as usize;
+    println!("  2D reference: {nx}×{ny}, h {h:.4e}, dt {dt:.3e}, {steps} steps");
+    let mu = RHO * NU;
+    let start = std::time::Instant::now();
+    let mut drag = Vec::new();
+    let mut lift = Vec::new();
+    let mut worst = 0.0_f64;
+    for step in 0..steps {
+        let r = solver.advance(&mut field, dt).await?;
+        worst = worst.max(r.solver_result.final_residual);
+        let t = (step + 1) as f64 * dt;
+        if (step + 1) % 10 == 0 {
+            let mask = solver.mask().expect("mask");
+            let body = solver.body().expect("body");
+            let (mut fx, mut fy, mut skipped) = (0.0, 0.0, 0usize);
+            for (x, y, nx_, ny_, ds) in samples(t, 0.5 * h) {
+                match mask.traction_at(body, &field.u, &field.v, &field.p, mu, t, x, y, nx_, ny_) {
+                    Some((tx, ty)) => {
+                        fx += tx * ds;
+                        fy += ty * ds;
+                    }
+                    None => skipped += 1,
+                }
+            }
+            if t >= t_end - period {
+                drag.push(fx);
+                lift.push(fy);
+            }
+            if (step + 1) % 100 == 0 {
+                println!(
+                    "  t {t:6.3} tip {:+.4}: drag {fx:7.1} lift {fy:+8.1} N/m (skipped {skipped}); residual {:.1e}; [{:.0} s]",
+                    deflection(1.0, t).0,
+                    r.solver_result.final_residual,
+                    start.elapsed().as_secs_f64()
+                );
+            }
+        }
+    }
+    let filt: Vec<f64> = (0..lift.len())
+        .map(|i| median(&lift[i.saturating_sub(5)..(i + 6).min(lift.len())]))
+        .collect();
+    let lo = filt.iter().cloned().fold(f64::INFINITY, f64::min);
+    let hi = filt.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
+    println!(
+        "  FINAL 2D reference: last period drag mean {:.1} N/m, lift swing {:+.1} … {:+.1} (raw {:+.1} … {:+.1}); worst residual {worst:.1e}; {:.0} s",
+        drag.iter().sum::<f64>() / drag.len() as f64,
+        lo,
+        hi,
+        lift.iter().cloned().fold(f64::INFINITY, f64::min),
+        lift.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
+        start.elapsed().as_secs_f64()
+    );
+    Ok(())
+}