Add prescribed kinematics for immersed boundaries

docs/documentation/case.md

@@ -366,6 +366,17 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | `moving_ibm`         | Integer | Sets the method used for IB movement. |
 | `vel(i)`             | Real    | Initial velocity of the moving IB in the i-th direction. |
 | `angular_vel(i)`     | Real    | Initial angular velocity of the moving IB in the i-th direction. |
+| `kin_model`          | Integer | Prescribed kinematics (requires `moving_ibm = 1`, 3D): [0] off; [1] hinged flapping (roll + pitch); [2] smoothed pitch ramp and hold. |
+| `kin_hinge(i)`       | Real    | Hinge point, i-th component. |
+| `kin_offset(i)`      | Real    | Body-frame vector from the hinge to the patch centroid, i-th component. |
+| `kin_phi0`, `kin_theta0` | Real | Roll and pitch amplitudes (rad). |
+| `kin_theta_mean`     | Real    | Mean pitch angle (rad), held before onset and superposed after. |
+| `kin_freq`           | Real    | Flapping frequency (cycles per unit time). |
+| `kin_phase`          | Real    | Pitch phase lead relative to roll (rad); `pi/2` makes pitch lead by a quarter cycle. |
+| `kin_t0`             | Real    | Onset time of flapping. |
+| `kin_ramp`           | Real    | Duration of the raised-cosine amplitude ramp after onset (0 = instantaneous). |
+| `kin_pitch_rate`     | Real    | `kin_model = 2`: nominal pitch rate \f$\Omega\f$ (rad per unit time); the pitch time is `kin_theta0`/\f$\Omega\f$. |
+| `kin_smooth`         | Real    | `kin_model = 2`: smoothing parameter \f$a\f$ of the Eldredge log-cosh ramp (11 in the AIAA canonical cases). |
 | `coefficient_of_restitution`     | Real    | A number 0 to 1 describing how elastic IB collisions are |
 | `collision_model`     | Integer    | Integer to select the collision model being used for IB collisions. |
 | `collision_time`     | Real    | Amount of simulation time used to resolve collisions |
@@ -411,6 +422,10 @@ Additional details on this specification can be found in [NACA airfoil](https://
 - `angular_vel(i)` is the initial angular velocity of the IB about the x, y, z axes for i=1, 2, 3 in radians per second. When `moving_ibm` equals 2, this rotation rate is just the starting rate of the object, which will then change due to external torques. If `moving_ibm` equals 1, then this is constant if it is a number, or can be described analytically with an expression.
 
   Moving-IB analytic expressions use the same Python syntax and error-reporting as IC patch expressions (see the "Analytical Definition of Primitive Variables" section above).
+
+- `kin_model = 1` prescribes hinged flapping kinematics at run time (no analytic expressions, so the binary is shared across parameter values): roll \f$\phi\f$ about the lab \f$x\f$ axis through `kin_hinge` and pitch \f$\theta\f$ about the body spanwise (\f$y\f$) axis through the hinge, composed as \f$R = R_x(\phi) R_y(\theta)\f$. With \f$\tau = t - t_0\f$ and amplitude envelope \f$A(\tau)\f$ (0 before onset, raised cosine over `kin_ramp`, then 1): \f$\phi = A \phi_0 \sin(2\pi f \tau)\f$, \f$\theta = \theta_m + A \theta_0 \sin(2\pi f \tau + \psi)\f$. The centroid follows \f$x_c = x_h + R\,\mathbf{r}_\mathrm{off}\f$ and the ghost-cell velocities use the lab-frame angular velocity \f$\dot\phi \mathbf{e}_x + \dot\theta R_x(\phi)\mathbf{e}_y\f$. Set the initial `x[y,z]_centroid` and `angles` consistently with \f$t = 0\f$ so pre-process marks the body in the right place.
+
+- `kin_model = 2` is the smoothed linear pitch-ramp-and-hold of the AIAA low-Reynolds-number canonical cases (Eldredge et al. 2009, Ol et al. 2010) about the hinge, with no roll: \f$\theta(t) = \theta_m + \frac{\theta_0}{2}\left[1 + \frac{1}{a t_p}\log\frac{\cosh(a\tau)}{\cosh(a(\tau - t_p))}\right]\f$, \f$\tau = t - t_0\f$, \f$t_p = \theta_0/\Omega\f$, so the angle rises from `kin_theta_mean` by `kin_theta0` at nominal rate `kin_pitch_rate` starting at `kin_t0`, smoothed by `kin_smooth`. The same hinge, offset and centroid conventions as `kin_model = 1` apply.
   Available variables: `x` (`x_cc(i)`), `y` (`y_cc(j)`), `z` (`z_cc(k)`), `t` (current simulation time), and `r` (the IB patch radius).
   The same intrinsic functions and `pi` constant apply; bare `e` is not available.
 

examples/3D_ibm_flapping_plate/case.py

@@ -0,0 +1,136 @@
+#!/usr/bin/env python3
+# Flapping flat plate: a rigid wing in steady glide begins to flap, driven by the runtime prescribed-kinematics
+# option (patch_ib%kin_model = 1) rather than by analytic expressions in the case file. See readme.md.
+import json
+import math
+
+U, rho, gamma, Ma, Re = 1.0, 1.0, 1.4, 0.2, 1000.0
+P = rho * U**2 / (gamma * Ma**2)
+cs = math.sqrt(gamma * P / rho)
+
+c = 1.0  # chord
+semi_span = 2.0 * c  # aspect ratio 4 overall, half of it simulated
+phi0 = math.radians(30.0)  # roll amplitude
+theta0 = math.radians(20.0)  # pitch amplitude
+theta_glide = math.radians(5.0)  # incidence held during the glide
+St = 0.3  # Strouhal number based on the peak-to-peak tip excursion
+A_tip = 2.0 * semi_span * math.sin(phi0)
+freq = St * U / A_tip
+T = 1.0 / freq
+t0 = 2.0  # flapping begins here; before it the wing glides
+t_end = t0 + 0.75 * T
+
+dx = 0.06 * c
+thick = 4 * dx  # at least four cells across the section: see readme.md
+x0, x1 = -2.0 * c, 5.0 * c
+y0, y1 = 0.0, 3.0 * c  # y = 0 is the symmetry plane at the wing root
+z0, z1 = -2.5 * c, 2.5 * c
+m, n, p = int((x1 - x0) / dx) - 1, int((y1 - y0) / dx) - 1, int((z1 - z0) / dx) - 1
+dt = 0.4 * dx / (U + cs)
+nt = int(t_end / dt)
+
+case = {
+    "run_time_info": "T",
+    "parallel_io": "T",
+    "prim_vars_wrt": "T",
+    "ib_state_wrt": "T",
+    "format": "silo",
+    "precision": "double",
+    "x_domain%beg": x0,
+    "x_domain%end": x1,
+    "y_domain%beg": y0,
+    "y_domain%end": y1,
+    "z_domain%beg": z0,
+    "z_domain%end": z1,
+    "m": m,
+    "n": n,
+    "p": p,
+    "cyl_coord": "F",
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": nt,
+    "t_step_save": max(1, nt // 20),
+    "num_patches": 1,
+    "num_fluids": 1,
+    "model_eqns": "5eq",
+    "alt_soundspeed": "F",
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": "rk3",
+    "weno_order": 5,
+    "weno_eps": 1.0e-10,
+    "weno_Re_flux": "T",
+    "weno_avg": "T",
+    "avg_state": "arithmetic",
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "F",
+    "riemann_solver": "hllc",
+    "low_Mach": 2,
+    "wave_speeds": "direct",
+    "viscous": "T",
+    "fd_order": 4,
+    "patch_icpp(1)%geometry": 9,
+    "patch_icpp(1)%x_centroid": 0.5 * (x0 + x1),
+    "patch_icpp(1)%y_centroid": 0.5 * (y0 + y1),
+    "patch_icpp(1)%z_centroid": 0.5 * (z0 + z1),
+    "patch_icpp(1)%length_x": x1 - x0,
+    "patch_icpp(1)%length_y": y1 - y0,
+    "patch_icpp(1)%length_z": z1 - z0,
+    "patch_icpp(1)%vel(1)": U,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%vel(3)": 0.0,
+    "patch_icpp(1)%pres": P,
+    "patch_icpp(1)%alpha_rho(1)": rho,
+    "patch_icpp(1)%alpha(1)": 1.0,
+    "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+    "fluid_pp(1)%eos": "ideal_gas",
+    "fluid_pp(1)%Re(1)": Re,
+    "bc_x%beg": -7,
+    "bc_x%grcbc_in": "T",
+    "bc_x%vel_in(1)": U,
+    "bc_x%vel_in(2)": 0.0,
+    "bc_x%vel_in(3)": 0.0,
+    "bc_x%pres_in": P,
+    "bc_x%alpha_rho_in(1)": rho,
+    "bc_x%alpha_in(1)": 1.0,
+    "bc_x%end": -8,
+    "bc_x%grcbc_out": "T",
+    "bc_x%pres_out": P,
+    "bc_y%beg": -2,
+    "bc_y%end": -9,
+    "bc_z%beg": -9,
+    "bc_z%end": -9,
+    "ib": "T",
+    "num_ibs": 1,
+    "patch_ib(1)%geometry": 9,
+    # the plate reaches past the symmetry plane so its root is a continuation of the wing, not a wall
+    "patch_ib(1)%x_centroid": 0.25 * c,
+    "patch_ib(1)%y_centroid": 0.75 * c,
+    "patch_ib(1)%z_centroid": 0.0,
+    "patch_ib(1)%length_x": c,
+    "patch_ib(1)%length_y": semi_span + 0.5 * c,
+    "patch_ib(1)%length_z": thick,
+    "patch_ib(1)%slip": "F",
+    "patch_ib(1)%moving_ibm": 1,
+    "patch_ib(1)%angles(1)": 0.0,
+    "patch_ib(1)%angles(2)": theta_glide,
+    "patch_ib(1)%angles(3)": 0.0,
+    # prescribed kinematics: roll about the lab x axis through the hinge, pitch about the body spanwise axis
+    # through the same hinge, pitch leading roll by a quarter cycle
+    "patch_ib(1)%kin_model": 1,
+    "patch_ib(1)%kin_hinge(1)": 0.0,
+    "patch_ib(1)%kin_hinge(2)": 0.0,
+    "patch_ib(1)%kin_hinge(3)": 0.0,
+    "patch_ib(1)%kin_offset(1)": 0.25 * c,
+    "patch_ib(1)%kin_offset(2)": 0.75 * c,
+    "patch_ib(1)%kin_offset(3)": 0.0,
+    "patch_ib(1)%kin_phi0": phi0,
+    "patch_ib(1)%kin_theta0": theta0,
+    "patch_ib(1)%kin_theta_mean": theta_glide,
+    "patch_ib(1)%kin_freq": freq,
+    "patch_ib(1)%kin_phase": math.pi / 2,
+    "patch_ib(1)%kin_t0": t0,
+    "patch_ib(1)%kin_ramp": 0.5 * T,
+}
+print(json.dumps(case, indent=4))

examples/3D_ibm_flapping_plate/readme.md

@@ -0,0 +1,87 @@
+# 3D Flapping Flat Plate (prescribed immersed-boundary kinematics)
+
+## Case Set Up
+
+A rigid rectangular plate of aspect ratio 4 glides at a fixed incidence and then begins to flap. Roll is about
+the streamwise axis through a hinge at the wing root; pitch is about the body spanwise axis through the same
+hinge, leading roll by a quarter cycle so the wing feathers through stroke reversal. Half the span is simulated
+with a symmetry plane at the root and a free tip.
+
+| | |
+|---|---|
+| Chord Reynolds number | 1 000 |
+| Mach number | 0.2 |
+| Roll amplitude | 30 degrees |
+| Pitch amplitude | 20 degrees, mean 5 degrees (the glide incidence) |
+| Strouhal number | 0.3, on the peak-to-peak tip excursion |
+| Onset | flapping ramps up over half a cycle after 2 convective times of glide |
+
+The motion is set by `patch_ib(1)%kin_model = 1` — the body state is evaluated from the current time at each
+Runge-Kutta stage, rather than integrated from analytic velocity expressions. Two consequences worth knowing:
+one case-optimized binary serves every Strouhal number and every ensemble member, because nothing about the
+motion is compiled in; and restarts reproduce the trajectory exactly, because nothing is integrated.
+
+Note the plate is defined to reach 0.5 c **past** the symmetry plane. An immersed boundary marks only cells
+inside the body, so a plate ending exactly on the plane would have fluid on the far side of its root and the root
+would be reconstructed as a no-slip face rather than the middle of a continuous wing.
+
+## Numerics
+
+WENO5 with the HLLC Riemann solver and the Thornber low-Mach velocity correction; `mp_weno` off, since the flow
+is smooth and the monotonicity clipping only costs accuracy here. Characteristic inflow and outflow with the
+generalized relaxation treatment. Fourth-order finite differences for the immersed boundary.
+
+The section is set to four cells thick deliberately. An immersed body reconstructs its interior cells from the
+fluid outside, and with fewer than about four cells across there is no real interior left: on a 2D pitching plate
+the lift through a ramp was 30 to 40 percent high at two cells and converged by four. Keep that ratio if you
+change `dx`.
+
+## Verifying it
+
+The prescribed motion has a closed form, so this case can be checked without any reference data. Read
+`D/ib1_forces.dat`, which carries one line per time step, and compare columns `ax, ay` (roll and pitch angles),
+`xc, yc, zc` (centroid), `vx, vy, vz` (centroid velocity) and `wx, wy, wz` (lab-frame angular velocity) against
+
+```
+phi(t)   = A(tau) phi_0 sin(2 pi f tau),        tau = t - t_0
+theta(t) = theta_m + A(tau) theta_0 sin(2 pi f tau + pi/2)
+centroid = hinge + Rx(phi) Ry(theta) offset
+omega    = phi' e_x + theta' Rx(phi) e_y
+```
+
+with `A` the raised-cosine onset envelope. All of these agree to round-off (1e-15) at every step.
+
+The angular velocity is the interesting one. It is **not** the vector of Euler-angle rates: with the rotation
+composed as `R = Rx(phi) Ry(theta)`, the lab-frame angular velocity is `phi' e_x + theta' Rx(phi) e_y`. Using
+`(phi', theta', 0)` instead — which is what the analytic-expression path effectively does, since it advances the
+angles componentwise and then uses the same array as a vector — is wrong once both angles are moving. Sampling
+the velocity carried by the body cells adjacent to the fluid and comparing with the rigid-body velocity
+distinguishes them: 8e-13 with the correct vector, and 4e-3 to 1e-1 with the Euler rates, at roll angles near 30
+degrees.
+
+Two physical checks are worth running as well, because they need no reference either:
+
+- During the glide, before flapping starts, the wing should carry the lift of a finite wing at its incidence.
+  Lifting-line theory gives `C_L = 2 pi alpha / (1 + 2/AR)` = 0.37 at 5 degrees and AR 4; the case measures 0.33.
+- Drag during the glide must be positive.
+
+## A caution on the forces
+
+The lift from `D/ib1_forces.dat` behaves sensibly here, but **the streamwise force on a thin inclined plate is
+currently unreliable** — see issue #1849. The transverse component is unaffected. If you need drag or thrust
+from a case like this one, cross-check it against a control-volume momentum balance rather than taking the
+volume integral at face value.
+
+## Running
+
+```bash
+./mfc.sh run examples/3D_ibm_flapping_plate/case.py -n 8 --gpu
+```
+
+0.48 M cells and 1 750 steps, which is three quarters of a flapping cycle: enough to see the leading-edge and
+tip vortices form after onset, and to check the kinematics. Raise `t_end` for the approach to a periodic state.
+
+## References
+
+Prescribed kinematics of this form, and the Strouhal range, follow the animal-cruise literature; the canonical
+pitch-ramp companion case is `examples/3D_ibm_pitchup_plate`.

examples/3D_ibm_pitchup_plate/case.py

@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+# Canonical pitch-up: a flat plate pitches about its leading edge from 0 to 45 degrees on the smoothed linear
+# ramp of the AIAA low-Reynolds-number canonical cases, driven by patch_ib%kin_model = 2. See readme.md.
+import json
+import math
+
+U, rho, gamma, Ma, Re = 1.0, 1.0, 1.4, 0.2, 300.0
+P = rho * U**2 / (gamma * Ma**2)
+cs = math.sqrt(gamma * P / rho)
+
+c = 1.0  # chord
+semi_span = 2.0 * c  # aspect ratio 4 overall, half of it simulated
+th_max = math.radians(45.0)  # the canonical maneuver pitches from 0 to 45 degrees
+K = math.pi / 8  # reduced pitch rate: the plate pitches over one chord of travel (case C1)
+Omega = 2 * K * U  # nominal pitch rate, rad per c/U
+t_p = th_max / Omega  # pitch duration
+a_smooth = 21.0  # Eldredge smoothing for C1
+t0 = 3.0  # the ramp starts here; before it the plate sits at zero incidence
+t_end = t0 + t_p + 4.0
+
+dx = 0.0125 * c
+thick = 4 * dx  # at least four cells across the section: see readme.md
+x0, x1 = -2.0 * c, 4.0 * c
+y0, y1 = 0.0, 2.8 * c  # y = 0 is the symmetry plane at the wing root
+z0, z1 = -2.0 * c, 2.0 * c
+m, n, p = int((x1 - x0) / dx) - 1, int((y1 - y0) / dx) - 1, int((z1 - z0) / dx) - 1
+dt = 0.4 * dx / (U + cs)
+nt = int(t_end / dt)
+
+case = {
+    "run_time_info": "T",
+    "parallel_io": "T",
+    "prim_vars_wrt": "T",
+    "ib_state_wrt": "T",
+    "format": "silo",
+    "precision": "double",
+    "x_domain%beg": x0,
+    "x_domain%end": x1,
+    "y_domain%beg": y0,
+    "y_domain%end": y1,
+    "z_domain%beg": z0,
+    "z_domain%end": z1,
+    "m": m,
+    "n": n,
+    "p": p,
+    "cyl_coord": "F",
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": nt,
+    "t_step_save": max(1, nt // 20),
+    "num_patches": 1,
+    "num_fluids": 1,
+    "model_eqns": "5eq",
+    "alt_soundspeed": "F",
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": "rk3",
+    "weno_order": 5,
+    "weno_eps": 1.0e-10,
+    "weno_Re_flux": "T",
+    "weno_avg": "T",
+    "avg_state": "arithmetic",
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "F",
+    "riemann_solver": "hllc",
+    "low_Mach": 2,
+    "wave_speeds": "direct",
+    "viscous": "T",
+    "fd_order": 4,
+    "patch_icpp(1)%geometry": 9,
+    "patch_icpp(1)%x_centroid": 0.5 * (x0 + x1),
+    "patch_icpp(1)%y_centroid": 0.5 * (y0 + y1),
+    "patch_icpp(1)%z_centroid": 0.5 * (z0 + z1),
+    "patch_icpp(1)%length_x": x1 - x0,
+    "patch_icpp(1)%length_y": y1 - y0,
+    "patch_icpp(1)%length_z": z1 - z0,
+    "patch_icpp(1)%vel(1)": U,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%vel(3)": 0.0,
+    "patch_icpp(1)%pres": P,
+    "patch_icpp(1)%alpha_rho(1)": rho,
+    "patch_icpp(1)%alpha(1)": 1.0,
+    "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+    "fluid_pp(1)%eos": "ideal_gas",
+    "fluid_pp(1)%Re(1)": Re,
+    "bc_x%beg": -7,
+    "bc_x%grcbc_in": "T",
+    "bc_x%vel_in(1)": U,
+    "bc_x%vel_in(2)": 0.0,
+    "bc_x%vel_in(3)": 0.0,
+    "bc_x%pres_in": P,
+    "bc_x%alpha_rho_in(1)": rho,
+    "bc_x%alpha_in(1)": 1.0,
+    "bc_x%end": -8,
+    "bc_x%grcbc_out": "T",
+    "bc_x%pres_out": P,
+    "bc_y%beg": -2,
+    "bc_y%end": -9,
+    "bc_z%beg": -9,
+    "bc_z%end": -9,
+    "ib": "T",
+    "num_ibs": 1,
+    "patch_ib(1)%geometry": 9,
+    # the plate reaches past the symmetry plane so its root is a continuation of the wing, not a wall
+    "patch_ib(1)%x_centroid": 0.5 * c,
+    "patch_ib(1)%y_centroid": 0.75 * c,
+    "patch_ib(1)%z_centroid": 0.0,
+    "patch_ib(1)%length_x": c,
+    "patch_ib(1)%length_y": semi_span + 0.5 * c,
+    "patch_ib(1)%length_z": thick,
+    "ib_neighborhood_radius": 3,
+    "patch_ib(1)%slip": "F",
+    "patch_ib(1)%moving_ibm": 1,
+    "patch_ib(1)%angles(1)": 0.0,
+    "patch_ib(1)%angles(2)": 0.0,
+    "patch_ib(1)%angles(3)": 0.0,
+    # prescribed kinematics: the Eldredge smoothed linear pitch ramp and hold, about the leading edge
+    "patch_ib(1)%kin_model": 2,
+    "patch_ib(1)%kin_hinge(1)": 0.0,
+    "patch_ib(1)%kin_hinge(2)": 0.0,
+    "patch_ib(1)%kin_hinge(3)": 0.0,
+    "patch_ib(1)%kin_offset(1)": 0.5 * c,
+    "patch_ib(1)%kin_offset(2)": 0.75 * c,
+    "patch_ib(1)%kin_offset(3)": 0.0,
+    "patch_ib(1)%kin_theta0": th_max,
+    "patch_ib(1)%kin_theta_mean": 0.0,
+    "patch_ib(1)%kin_pitch_rate": Omega,
+    "patch_ib(1)%kin_smooth": a_smooth,
+    "patch_ib(1)%kin_t0": t0,
+}
+print(json.dumps(case, indent=4))