diff --git a/AGENTS.md b/AGENTS.md index 2f314b3..2679e26 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -23,6 +23,10 @@ coefficients. This file holds only directives that are specifically for you (the AI assistant) and would be noise in the README. +## Critical: Sign conventions + +If you realize you need to make the decision about a sign convention, first validate what Ardupilot uses and prefer that, then stick with helicopter convention, then wind turbine convention. Once you decide, document the decision in this AGENTS.md so that future code is written using this convention. + ## Rotor rotation direction — CCW from above (American convention) **This project uses the American helicopter convention: the rotor spins @@ -124,6 +128,44 @@ See [PITT_PETERS_DESIGN.md](docs/PITT_PETERS_DESIGN.md) for the full coefficient definitions, the cross-product derivation, and the BladeAD sign differences (our C_L_hub = -BladeAD C_Mx, our C_M_hub = +BladeAD C_My). +## In-plane hub force (H-force) sign convention + +The rotor's net in-plane hub force (classical "H-force") is assembled in +`dynbem_rs/src/bem_common.rs` and validated by +`validation_rs/src/checks/h_force.rs`. The signs below are load-bearing; +do not flip them without re-running that check. + +### Profile / induced-drag H-force (implemented) + +Each blade element's tangential aerodynamic force `dFt = dQ / r` (the +same force whose moment produces `dQ`, recovered by undoing the `* r` +weighting) is projected onto the fixed hub axes with the SAME +`(sin psi, cos psi)` pairing used everywhere else in the sweep +(`v_t_extra`, `Mx_hub`/`My_hub`), then azimuth-averaged: + + Fx_hub += (dQ / r) * sin(psi) + Fy_hub += (dQ / r) * cos(psi) + +**Decision (which convention this project uses):** the pairing above is +the chosen convention, pinned by two validated directional cases: + +- Disk level, hub translating at +X through still air (apparent headwind + from the nose): the resulting force is REARWARD (-X) -- drag opposing + the direction of flight, matching the textbook H-force definition. +- Stationary disk, +Y crosswind: force is +Y (downwind). + +The term is zero for azimuth-independent (pure axial / hover) loading and +grows with edgewise flow, vanishing again as the hub axis re-aligns with +the relative wind. See the `assemble_result` / `SweepCtx::run` doc +comments for the full rationale. + +### Flapping-tilt H-force + +The second classical H-force contribution -- the thrust vector tilting +into the disk plane by the local flap angle beta(psi) -- and its sign +convention are documented in the blade-flapping section below, since it +shares the flap-angle sign convention and the harmonic flap solve. + ## Pitt-Peters inflow model Implementation: `dynbem_rs/src/pitt_peters.rs`. @@ -275,12 +317,14 @@ is not consumed by the Rust aero solvers. Use the `pitch_actuation:` YAML block (with a `servoflap:` sub-block) to enable active servo-flap feathering dynamics. -## Quasi-static blade flapping (hub moment reduction) +## Quasi-static blade flapping (hub moment reduction + flapback H-force) `FlapProperties` (in `dynbem_rs/src/rotor_definition.rs`) models -out-of-plane blade flexibility as an equivalent spring-hinge. The blade -absorbs most aerodynamic pitching/rolling moment via deflection; only a -fraction reaches the airframe (hub). +out-of-plane blade flexibility as an equivalent spring-hinge and solves a +phase-correct 1/rev flap response. The blade absorbs most of the +aerodynamic 1/rev moment via deflection; only the fraction set by the +centrifugal stiffness reaches the airframe (hub moment), and the tilted +disk contributes an in-plane hub force (the flapping-tilt "H-force"). Parameters: @@ -288,23 +332,72 @@ Parameters: - `omega_nr_rad_s` -- non-rotating flap natural frequency [rad/s] (K_beta = I_b * omega_NR^2). 0.0 = freely hinged (no spring). -Physics: +### Flap-angle sign convention (load-bearing) + +Project frame: psi=0 at +X (nose/North), CCW from above (see +rotor-rotation section). The flap angle is + + beta(psi) = beta_0 + beta_1c*cos(psi) + beta_1s*sin(psi) + +with **beta > 0 = blade tip UP** (toward -z, the thrust direction). A +disk that flaps up at the nose (beta_1c > 0) is tilted AFT (thrust vector +leans toward -X), which is the classical flapback. + +### Phase-correct 1/rev solve + +A real rotor running near nu_beta ~ 1.0-1.15 (1/rev forcing near +resonance) responds to the cyclic aero flap moment with a ~90 deg lag +driven by aerodynamic flap damping. That lag is what points the flapback +aft rather than sideways, so the flap solve must be phase-correct (a +magnitude-only scaling of the hub moment cannot produce it). + +`apply_flap_dynamics()` in `dynbem_rs/src/bem_common.rs` implements the +damped 1/rev harmonic balance (' = d/dpsi): + + I_b*Omega^2 * (beta'' + d*beta' + nu^2*beta) = M_beta(psi) + nu_beta^2 = 1 + (omega_NR / Omega)^2 (centrifugal + spring) + d = gamma/8 = 0.5*rho*a*S3 / I_b (aero flap damping) + S3 = integral c(r)*r^3 dr (over the radial grid; Omega cancels) + +The sweep's azimuth-averaged hub moments ARE the 1/rev aero flap-moment +harmonics: `mx_hub = N_b*M1s/2`, `my_hub = N_b*M1c/2`. At 1/rev +(beta'' = -beta) the balance is the linear system + + [[a, d], [-d, a]] [beta_1c; beta_1s] = (1/(I_b*Omega^2)) [M1c; M1s], + a = nu^2 - 1 + +Transmitted hub moment (Johnson): +`M_hub = (N_b/2)*I_b*Omega^2*(nu^2-1)*beta_1`. Flapping-tilt H-force +(hub frame, T = mean disk thrust): + + Fx_flap = -T*beta_1c / 2 + Fy_flap = +T*beta_1s / 2 + +added to the profile-drag H-force before `assemble_result`. - nu_beta^2 = 1 + (omega_NR / Omega)^2 - hub_moment_factor = (nu_beta^2 - 1) / nu_beta^2 +### Decisions this convention pins -- Freely hinged (omega_NR=0): factor=0, no moment transfer. -- Rigid blade (omega_NR >> Omega): factor->1, full moment transfer. -- Typical hingeless rotor: nu_beta ~ 1.05-1.15, factor ~ 0.05-0.15. +- Freely hinged (omega_NR=0, a=0): transmits **zero** hub moment, but + STILL produces a flapback H-force -- damping alone gives beta_1 a pure + 90 deg lag. +- Rigid blade / no damping limit: transmitted moment -> full aero moment. +- **H-force direction (validated, not hand-traced):** forward flight along + +X through still air with real thrust gives a REARWARD (-X) flapping + H-force that adds to (and dominates) the profile-drag term. Pinned by + the `flapback` case in `validation_rs/src/checks/h_force.rs` + (flap F_north ~ -2.86 N vs rigid -1.37 N at Lock ~8; free hinge + transmits ~0 moment). If you ever change the flap-angle sign, flip + `dfx_hub` and `dfy_hub` together (a beta sign flip) and re-run that + check. -Implementation: +Implementation notes: -- `apply_flap_reduction()` in `dynbem_rs/src/bem_common.rs` scales - `mx_hub, my_hub` by the factor before `assemble_result`. -- Applied in all three models (Pitt-Peters, Oye, quasi-static BEM). -- The inflow ODE (Pitt-Peters) still uses full aerodynamic moments - (the wake responds to disk loading, not what the airframe sees). +- `apply_flap_dynamics()` is applied in all three BEM-family models + (Pitt-Peters, Oye, quasi-static BEM) at the result-assembly step only. +- The inflow ODE (Pitt-Peters) still uses the FULL aerodynamic moments + (the wake responds to disk loading, not to what the airframe sees). - Thrust and torque are unchanged by flapping. +- `FlapProperties::nu_beta_sq()` feeds the harmonic solve's stiffness term. YAML schema: @@ -312,10 +405,11 @@ YAML schema: I_blade_flap_kgm2: 0.012 omega_nr_rad_s: 8.0 -Absent `flap:` section = rigid blade (no moment reduction), preserving -backward compatibility. +Absent `flap:` section = rigid blade (no moment reduction, no flapping +H-force), preserving backward compatibility. -Tests: `tests/test_flap_hinge.py`. +Tests: `tests/test_flap_hinge.py`, +`validation_rs/src/checks/h_force.rs` (flapback case). ## Windmill solver: non-axial v_t_extra extension diff --git a/Cargo.toml b/Cargo.toml index 30e668f..83585a9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,7 +4,7 @@ members = ["dynbem", "dynbem_rs", "validation_rs"] # Fields below are inherited by member crates via `field.workspace = true`. [workspace.package] -version = "0.7.0" +version = "0.8.0" edition = "2021" license = "MIT" authors = ["Kristof"] diff --git a/README.md b/README.md index 9718c5c..ab17f24 100644 --- a/README.md +++ b/README.md @@ -181,7 +181,7 @@ enabled in the release profile for symbol resolution). | Document | Content | |---|---| -| [BEM_COMMON.md](docs/BEM_COMMON.md) | Coordinate system, kinematics, force kernel, QS BEM, VRS, servo-flap, output assembly, model comparison | +| [BEM_COMMON.md](docs/BEM_COMMON.md) | Coordinate system, kinematics, force kernel, QS BEM, VRS, servo-flap, blade flapping, in-plane hub force (H-force), output assembly, model comparison | | [PITT_PETERS_DESIGN.md](docs/PITT_PETERS_DESIGN.md) | Pitt-Peters 3-state dynamic inflow (formal math + implementation notes) | | [OYE_DESIGN.md](docs/OYE_DESIGN.md) | Oye 2-stage annular dynamic inflow | | [VPM_DESIGN.md](docs/VPM_DESIGN.md) | VPM free-wake solver | diff --git a/Research/Harris_CR-2008-215370/fig_5_18_pca2_hforce.md b/Research/Harris_CR-2008-215370/fig_5_18_pca2_hforce.md new file mode 100644 index 0000000..0f7a3d6 --- /dev/null +++ b/Research/Harris_CR-2008-215370/fig_5_18_pca2_hforce.md @@ -0,0 +1,98 @@ +# Fig. 5-18 -- PCA-2 Rotor H-Force Coefficient vs Advance Ratio + +**Source:** Harris, F.D. (2008). *Introduction to Autogyros, Helicopters, +and Other V/STOL Aircraft* -- NASA/CR-2008-215370, Fig. 5-18 (p.88). +Re-analysis of the Wheatley & Hood PCA-2 wind-tunnel data (NACA TR 515), +converted to rotor-axis coefficients. + +**Rotor:** Pitcairn PCA-2 autogyro, 4 blades, D = 45 ft (R = 6.86 m), +untwisted, autorotating (no shaft power). Same rotor as +`Research/Wheatley_Hood_NACA515/`. + +## Coefficient definition (confirm against Harris symbol page) + +Harris uses the American (Army/NACA) coefficient convention, NOT the +1/2-rho dynamic-pressure form: + + CH = H / (rho * pi * R^2 * (Omega*R)^2) + +where H is the rotor in-plane (H-) force along the rotor-axis drag +direction, rho air density, R tip radius, Omega rotor speed. Fig. 5-18 +plots the BARE CH (Fig. 5-25 by contrast plots CH/sigma). The exact +reference area/velocity should be verified against Harris's front-matter +symbol list (image-only in the PDF) before using these numbers +dimensionally. To recover a force: H = CH * rho * pi * R^2 * (Omega*R)^2. + +Nominal test RPMs (from the paper): 98.6, 118.7, 137.6, 147.9, giving +Omega*R ~ 232, 280, 324, 348 ft/s (~70.7, 85.3, 98.8, 106.1 m/s). + +## Confidence: LOW (visual digitization) + +These points were read by eye from the rendered scatter plot, NOT from a +tabulated source. Estimated reading tolerance: mu +/- 0.01, +CH +/- 0.00003. Use as a trend/order-of-magnitude validation target. The +higher-fidelity primary data for this same rotor is the Wheatley & Hood +tables (`Research/Wheatley_Hood_NACA515/page_10_table_iii.md`, `_iv.md`), +which give CD/C_y in airplane axes at HIGH confidence; Fig. 5-18 is +Harris's rotor-axis re-projection of that data across all four RPM runs. + +The clear physical trend (the point of the figure): CH rises roughly +linearly with advance ratio, ~0.0002 at mu ~ 0.14 to ~0.001 at mu ~ 0.7, +and is essentially independent of RPM (all four series collapse onto one +curve). + +## Data (grouped by nominal RPM series) + +Column key -- mu: advance ratio | CH: rotor-axis H-force coefficient | +rpm_series: nominal RPM marker series in the figure. + +| mu | CH | rpm_series | +|-------|----------|------------| +| 0.13 | 0.00020 | 98.6 | +| 0.17 | 0.00024 | 98.6 | +| 0.205 | 0.00027 | 98.6 | +| 0.25 | 0.00030 | 98.6 | +| 0.27 | 0.00037 | 98.6 | +| 0.31 | 0.00039 | 98.6 | +| 0.345 | 0.00041 | 98.6 | +| 0.38 | 0.00049 | 98.6 | +| 0.41 | 0.00065 | 98.6 | +| 0.44 | 0.00057 | 98.6 | +| 0.46 | 0.00062 | 98.6 | +| 0.49 | 0.00066 | 98.6 | +| 0.53 | 0.00072 | 98.6 | +| 0.59 | 0.00078 | 98.6 | +| 0.635 | 0.00086 | 98.6 | +| 0.70 | 0.00097 | 98.6 | +| 0.14 | 0.00020 | 118.7 | +| 0.18 | 0.00021 | 118.7 | +| 0.22 | 0.00030 | 118.7 | +| 0.245 | 0.00040 | 118.7 | +| 0.285 | 0.00043 | 118.7 | +| 0.34 | 0.00047 | 118.7 | +| 0.375 | 0.00055 | 118.7 | +| 0.42 | 0.00062 | 118.7 | +| 0.50 | 0.00066 | 118.7 | +| 0.595 | 0.00077 | 118.7 | +| 0.605 | 0.00085 | 118.7 | +| 0.135 | 0.00007 | 137.6 | +| 0.205 | 0.00016 | 137.6 | +| 0.21 | 0.00028 | 137.6 | +| 0.24 | 0.00036 | 137.6 | +| 0.275 | 0.00042 | 137.6 | +| 0.30 | 0.00046 | 137.6 | +| 0.335 | 0.00050 | 137.6 | +| 0.37 | 0.00053 | 137.6 | +| 0.41 | 0.00060 | 137.6 | +| 0.45 | 0.00062 | 137.6 | +| 0.49 | 0.00070 | 137.6 | +| 0.52 | 0.00075 | 137.6 | +| 0.14 | 0.00018 | 147.9 | +| 0.175 | 0.00009 | 147.9 | +| 0.225 | 0.00031 | 147.9 | +| 0.28 | 0.00046 | 147.9 | +| 0.32 | 0.00054 | 147.9 | +| 0.35 | 0.00060 | 147.9 | +| 0.415 | 0.00060 | 147.9 | +| 0.135 | 0.00019 | transition | +| 0.72 | 0.00107 | transition | diff --git a/Research/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs.md b/Research/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs.md new file mode 100644 index 0000000..1c8458b --- /dev/null +++ b/Research/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs.md @@ -0,0 +1,96 @@ +# Harris CR-2008-215370 -- PCA-2 rotor tabulated coefficients (pages 504-505) + +Source: Harris, "Rotor Performance at High Advance Ratio: Theory versus +Test", NASA CR-2008-215370. Printed pages 504-505 (PDF page indices +511-512). Extracted directly from the report's embedded text layer +(not a figure digitization), so these values are HIGH confidence -- they +are the author's own tabulated numbers and are the primary source behind +Fig. 5-18 (Rotor CH vs advance ratio) and the related CT / CY figures. + +This supersedes the earlier eyeballed digitization in +`fig_5_18_pca2_hforce.md` for any quantitative use. + +## Column definitions + +Coefficients use the American (Harris) normalization +`C = force / (rho * pi * R^2 * (Omega*R)^2)` (bare coefficient, NOT +divided by solidity sigma). Rotor CH is the in-plane hub force +coefficient (the H-force), Rotor CT the thrust coefficient, Rotor CY the +side-force coefficient. + +| Column | Meaning | +| --- | --- | +| Advance Ratio | mu = V / (Omega*R) | +| Shaft AoA (deg) | rotor shaft angle of attack | +| Airplane CL | airplane (test article) lift coefficient | +| Airplane CD | airplane drag coefficient | +| L/D | airplane lift-to-drag ratio | +| RPM | rotor rotational speed | +| Airplane Cl_roll | airplane rolling moment coefficient | +| Airplane Cm_pitch | airplane pitching moment coefficient | +| Airplane Cy_lat | airplane lateral force coefficient | +| Rotor CT | rotor thrust coefficient | +| Rotor CH | rotor in-plane (H-force) coefficient | +| Rotor CY | rotor side-force coefficient | + +## Data + +| Advance Ratio | Shaft AoA (deg) | Airplane CL | Airplane CD | L/D | RPM | Airplane Cl_roll | Airplane Cm_pitch | Airplane Cy_lat | Rotor CT | Rotor CH | Rotor CY | +| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | +| 0.1435 | 16.98 | 0.494326 | 0.166101 | 2.986 | 98.63 | 0.000310 | -0.000191 | 0.010333 | 0.005862 | 0.000164 | 0.000116 | +| 0.1966 | 10.56 | 0.279403 | 0.065150 | 4.296 | 98.63 | 0.000136 | 0.000181 | 0.005568 | 0.005729 | 0.000257 | 0.000111 | +| 0.2478 | 7.53 | 0.176877 | 0.034163 | 5.183 | 98.63 | 0.000075 | 0.000212 | 0.003526 | 0.005617 | 0.000335 | 0.000110 | +| 0.2986 | 5.52 | 0.117805 | 0.020489 | 5.745 | 98.63 | 0.000060 | 0.000202 | 0.002249 | 0.005363 | 0.000408 | 0.000101 | +| 0.3490 | 4.41 | 0.085388 | 0.014319 | 5.962 | 98.63 | 0.000029 | 0.000202 | 0.001775 | 0.005279 | 0.000472 | 0.000109 | +| 0.3993 | 3.47 | 0.063016 | 0.010506 | 5.988 | 98.63 | 0.000013 | 0.000156 | 0.001175 | 0.005081 | 0.000533 | 0.000094 | +| 0.4494 | 3.08 | 0.050612 | 0.008575 | 5.904 | 98.63 | 0.000010 | 0.000169 | 0.001229 | 0.005161 | 0.000592 | 0.000124 | +| 0.4995 | 2.66 | 0.040677 | 0.007640 | 5.331 | 98.63 | -0.000003 | 0.000174 | 0.001258 | 0.005119 | 0.000718 | 0.000157 | +| 0.5496 | 2.16 | 0.032842 | 0.006428 | 5.114 | 98.63 | -0.000008 | 0.000160 | 0.001023 | 0.004996 | 0.000784 | 0.000155 | +| 0.5997 | 1.77 | 0.027555 | 0.005879 | 4.688 | 98.63 | -0.000008 | 0.000106 | 0.000886 | 0.004986 | 0.000904 | 0.000159 | +| 0.6498 | 1.29 | 0.022296 | 0.005152 | 4.327 | 98.63 | 0.000033 | 0.000095 | 0.000839 | 0.004729 | 0.000982 | 0.000177 | +| 0.6999 | 1.08 | 0.019634 | 0.004833 | 3.981 | 98.63 | 0.000004 | 0.000113 | 0.000966 | 0.004824 | 0.001093 | 0.000237 | +| 0.1435 | 16.97 | 0.502316 | 0.168932 | 2.976 | 118.7 | 0.000342 | 0.000057 | 0.010706 | 0.005957 | 0.000169 | 0.000120 | +| 0.1967 | 10.46 | 0.282081 | 0.065175 | 4.318 | 118.7 | 0.000163 | 0.000262 | 0.006041 | 0.005782 | 0.000258 | 0.000121 | +| 0.2479 | 7.48 | 0.179419 | 0.034337 | 5.242 | 118.7 | 0.000094 | 0.000279 | 0.003998 | 0.005696 | 0.000334 | 0.000125 | +| 0.2986 | 5.45 | 0.119231 | 0.020489 | 5.821 | 118.7 | 0.000074 | 0.000251 | 0.002492 | 0.005426 | 0.000409 | 0.000112 | +| 0.3490 | 4.31 | 0.086116 | 0.014187 | 6.054 | 118.7 | 0.000057 | 0.000222 | 0.001931 | 0.005322 | 0.000470 | 0.000118 | +| 0.3993 | 3.43 | 0.063047 | 0.010711 | 5.874 | 118.7 | 0.000038 | 0.000187 | 0.001337 | 0.005083 | 0.000554 | 0.000107 | +| 0.4494 | 3.00 | 0.050442 | 0.008411 | 5.983 | 118.7 | 0.000028 | 0.000204 | 0.001380 | 0.005142 | 0.000584 | 0.000140 | +| 0.4995 | 2.61 | 0.041019 | 0.007399 | 5.544 | 118.7 | 0.000029 | 0.000180 | 0.001367 | 0.005160 | 0.000691 | 0.000171 | +| 0.5496 | 2.06 | 0.032368 | 0.006025 | 5.368 | 118.7 | 0.000012 | 0.000171 | 0.001180 | 0.004921 | 0.000735 | 0.000179 | +| 0.5997 | 1.93 | 0.028991 | 0.005651 | 5.124 | 118.7 | -0.000005 | 0.000152 | 0.001191 | 0.005244 | 0.000840 | 0.000214 | +| 0.6497 | 1.63 | 0.024436 | 0.005461 | 4.479 | 118.7 | 0.000005 | 0.000155 | 0.000909 | 0.005187 | 0.001006 | 0.000192 | +| 0.6998 | 1.38 | 0.021369 | 0.005008 | 4.271 | 118.7 | -0.000006 | 0.000110 | 0.001036 | 0.005256 | 0.001100 | 0.000254 | +| 0.1438 | 16.54 | 0.511681 | 0.168126 | 3.05 | 137.6 | 0.000321 | 0.000253 | 0.011200 | 0.006054 | 0.000175 | 0.000126 | +| 0.1967 | 10.44 | 0.287443 | 0.066729 | 4.317 | 137.6 | 0.000181 | 0.000342 | 0.006339 | 0.005892 | 0.000271 | 0.000127 | +| 0.2479 | 7.36 | 0.181626 | 0.034224 | 5.307 | 137.6 | 0.000117 | 0.000317 | 0.004251 | 0.005763 | 0.000334 | 0.000133 | +| 0.2987 | 5.35 | 0.120000 | 0.020784 | 5.77 | 137.6 | 0.000112 | 0.000269 | 0.002590 | 0.005461 | 0.000428 | 0.000117 | +| 0.3490 | 4.26 | 0.086831 | 0.014101 | 6.148 | 137.6 | 0.000077 | 0.000223 | 0.002108 | 0.005365 | 0.000466 | 0.000129 | +| 0.3993 | 3.39 | 0.063909 | 0.010398 | 6.149 | 137.6 | 0.000039 | 0.000193 | 0.001442 | 0.005150 | 0.000528 | 0.000115 | +| 0.4494 | 2.97 | 0.051300 | 0.008533 | 6.025 | 137.6 | 0.000028 | 0.000202 | 0.001399 | 0.005228 | 0.000594 | 0.000142 | +| 0.4995 | 2.57 | 0.041302 | 0.007328 | 5.639 | 137.6 | 0.000026 | 0.000180 | 0.001416 | 0.005194 | 0.000683 | 0.000177 | +| 0.5496 | 2.10 | 0.033197 | 0.006169 | 5.379 | 137.6 | 0.000012 | 0.000181 | 0.001412 | 0.005046 | 0.000749 | 0.000214 | +| 0.5996 | 1.97 | 0.029218 | 0.005870 | 4.983 | 137.6 | 0.000034 | 0.000165 | 0.001203 | 0.005287 | 0.000875 | 0.000216 | +| 0.6497 | 1.81 | 0.025760 | 0.005344 | 4.826 | 137.6 | 0.000033 | 0.000146 | 0.001324 | 0.005467 | 0.000956 | 0.000280 | +| 0.6998 | 1.27 | 0.021027 | 0.004685 | 4.485 | 137.6 | 0.000011 | 0.000126 | 0.000860 | 0.005170 | 0.001033 | 0.000211 | +| 0.1438 | 16.52 | 0.516425 | 0.169777 | 3.046 | 147.9 | 0.000338 | 0.000316 | 0.011512 | 0.006110 | 0.000180 | 0.000130 | +| 0.1967 | 10.38 | 0.288981 | 0.066816 | 4.312 | 147.9 | 0.000191 | 0.000343 | 0.006578 | 0.005922 | 0.000273 | 0.000132 | +| 0.2479 | 7.35 | 0.182513 | 0.034661 | 5.269 | 147.9 | 0.000128 | 0.000339 | 0.004313 | 0.005793 | 0.000345 | 0.000135 | +| 0.2987 | 5.30 | 0.120677 | 0.020534 | 5.876 | 147.9 | 0.000100 | 0.000270 | 0.002611 | 0.005490 | 0.000419 | 0.000118 | +| 0.3490 | 4.23 | 0.087459 | 0.014184 | 6.17 | 147.9 | 0.000079 | 0.000237 | 0.002164 | 0.005403 | 0.000470 | 0.000133 | +| 0.3993 | 3.44 | 0.065334 | 0.010607 | 6.177 | 147.9 | 0.000066 | 0.000220 | 0.001446 | 0.005266 | 0.000533 | 0.000116 | +| 0.4494 | 2.91 | 0.051259 | 0.009104 | 5.636 | 147.9 | 0.000048 | 0.000192 | 0.001201 | 0.005227 | 0.000657 | 0.000122 | +| 0.4995 | 2.57 | 0.041501 | 0.007368 | 5.644 | 147.9 | 0.000038 | 0.000177 | 0.001426 | 0.005220 | 0.000688 | 0.000178 | +| 0.5496 | 2.06 | 0.033192 | 0.006260 | 5.307 | 147.9 | 0.000026 | 0.000160 | 0.001268 | 0.005046 | 0.000766 | 0.000192 | +| 0.5996 | 1.99 | 0.029903 | 0.005895 | 5.072 | 147.9 | 0.000028 | 0.000170 | 0.001173 | 0.005411 | 0.000874 | 0.000211 | +| 0.6497 | 1.71 | 0.025599 | 0.005539 | 4.627 | 147.9 | 0.000003 | 0.000154 | 0.001314 | 0.005431 | 0.001008 | 0.000278 | +| 0.6998 | 1.45 | 0.022245 | 0.005113 | 4.349 | 147.9 | -0.000005 | 0.000110 | 0.001209 | 0.005470 | 0.001115 | 0.000296 | + +## Notes + +- 48 rows = 4 RPM series (98.63, 118.7, 137.6, 147.9) x 12 advance-ratio + sweep points each (mu ~ 0.14 to 0.70). +- Rotor CH increases monotonically with mu (~0.00016 at mu~0.14 to + ~0.0011 at mu~0.70) and is essentially RPM-independent -- consistent + with Fig. 5-18. +- RPM labels are as printed in the report (98.63, 118.7, 137.6, 147.9). diff --git a/Research/csv/Harris_CR-2008-215370/fig_5_18_pca2_hforce.csv b/Research/csv/Harris_CR-2008-215370/fig_5_18_pca2_hforce.csv new file mode 100644 index 0000000..a3539f7 --- /dev/null +++ b/Research/csv/Harris_CR-2008-215370/fig_5_18_pca2_hforce.csv @@ -0,0 +1,49 @@ +mu,CH,rpm_series +0.13,0.00020,98.6 +0.17,0.00024,98.6 +0.205,0.00027,98.6 +0.25,0.00030,98.6 +0.27,0.00037,98.6 +0.31,0.00039,98.6 +0.345,0.00041,98.6 +0.38,0.00049,98.6 +0.41,0.00065,98.6 +0.44,0.00057,98.6 +0.46,0.00062,98.6 +0.49,0.00066,98.6 +0.53,0.00072,98.6 +0.59,0.00078,98.6 +0.635,0.00086,98.6 +0.70,0.00097,98.6 +0.14,0.00020,118.7 +0.18,0.00021,118.7 +0.22,0.00030,118.7 +0.245,0.00040,118.7 +0.285,0.00043,118.7 +0.34,0.00047,118.7 +0.375,0.00055,118.7 +0.42,0.00062,118.7 +0.50,0.00066,118.7 +0.595,0.00077,118.7 +0.605,0.00085,118.7 +0.135,0.00007,137.6 +0.205,0.00016,137.6 +0.21,0.00028,137.6 +0.24,0.00036,137.6 +0.275,0.00042,137.6 +0.30,0.00046,137.6 +0.335,0.00050,137.6 +0.37,0.00053,137.6 +0.41,0.00060,137.6 +0.45,0.00062,137.6 +0.49,0.00070,137.6 +0.52,0.00075,137.6 +0.14,0.00018,147.9 +0.175,0.00009,147.9 +0.225,0.00031,147.9 +0.28,0.00046,147.9 +0.32,0.00054,147.9 +0.35,0.00060,147.9 +0.415,0.00060,147.9 +0.135,0.00019,transition +0.72,0.00107,transition diff --git a/Research/csv/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs__column_definitions.csv b/Research/csv/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs__column_definitions.csv new file mode 100644 index 0000000..23eb8a5 --- /dev/null +++ b/Research/csv/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs__column_definitions.csv @@ -0,0 +1,13 @@ +Column,Meaning +Advance Ratio,mu = V / (Omega*R) +Shaft AoA (deg),rotor shaft angle of attack +Airplane CL,airplane (test article) lift coefficient +Airplane CD,airplane drag coefficient +L/D,airplane lift-to-drag ratio +RPM,rotor rotational speed +Airplane Cl_roll,airplane rolling moment coefficient +Airplane Cm_pitch,airplane pitching moment coefficient +Airplane Cy_lat,airplane lateral force coefficient +Rotor CT,rotor thrust coefficient +Rotor CH,rotor in-plane (H-force) coefficient +Rotor CY,rotor side-force coefficient diff --git a/Research/csv/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs__data.csv b/Research/csv/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs__data.csv new file mode 100644 index 0000000..0ff8e52 --- /dev/null +++ b/Research/csv/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs__data.csv @@ -0,0 +1,49 @@ +Advance Ratio,Shaft AoA (deg),Airplane CL,Airplane CD,L/D,RPM,Airplane Cl_roll,Airplane Cm_pitch,Airplane Cy_lat,Rotor CT,Rotor CH,Rotor CY +0.1435,16.98,0.494326,0.166101,2.986,98.63,0.000310,-0.000191,0.010333,0.005862,0.000164,0.000116 +0.1966,10.56,0.279403,0.065150,4.296,98.63,0.000136,0.000181,0.005568,0.005729,0.000257,0.000111 +0.2478,7.53,0.176877,0.034163,5.183,98.63,0.000075,0.000212,0.003526,0.005617,0.000335,0.000110 +0.2986,5.52,0.117805,0.020489,5.745,98.63,0.000060,0.000202,0.002249,0.005363,0.000408,0.000101 +0.3490,4.41,0.085388,0.014319,5.962,98.63,0.000029,0.000202,0.001775,0.005279,0.000472,0.000109 +0.3993,3.47,0.063016,0.010506,5.988,98.63,0.000013,0.000156,0.001175,0.005081,0.000533,0.000094 +0.4494,3.08,0.050612,0.008575,5.904,98.63,0.000010,0.000169,0.001229,0.005161,0.000592,0.000124 +0.4995,2.66,0.040677,0.007640,5.331,98.63,-0.000003,0.000174,0.001258,0.005119,0.000718,0.000157 +0.5496,2.16,0.032842,0.006428,5.114,98.63,-0.000008,0.000160,0.001023,0.004996,0.000784,0.000155 +0.5997,1.77,0.027555,0.005879,4.688,98.63,-0.000008,0.000106,0.000886,0.004986,0.000904,0.000159 +0.6498,1.29,0.022296,0.005152,4.327,98.63,0.000033,0.000095,0.000839,0.004729,0.000982,0.000177 +0.6999,1.08,0.019634,0.004833,3.981,98.63,0.000004,0.000113,0.000966,0.004824,0.001093,0.000237 +0.1435,16.97,0.502316,0.168932,2.976,118.7,0.000342,0.000057,0.010706,0.005957,0.000169,0.000120 +0.1967,10.46,0.282081,0.065175,4.318,118.7,0.000163,0.000262,0.006041,0.005782,0.000258,0.000121 +0.2479,7.48,0.179419,0.034337,5.242,118.7,0.000094,0.000279,0.003998,0.005696,0.000334,0.000125 +0.2986,5.45,0.119231,0.020489,5.821,118.7,0.000074,0.000251,0.002492,0.005426,0.000409,0.000112 +0.3490,4.31,0.086116,0.014187,6.054,118.7,0.000057,0.000222,0.001931,0.005322,0.000470,0.000118 +0.3993,3.43,0.063047,0.010711,5.874,118.7,0.000038,0.000187,0.001337,0.005083,0.000554,0.000107 +0.4494,3.00,0.050442,0.008411,5.983,118.7,0.000028,0.000204,0.001380,0.005142,0.000584,0.000140 +0.4995,2.61,0.041019,0.007399,5.544,118.7,0.000029,0.000180,0.001367,0.005160,0.000691,0.000171 +0.5496,2.06,0.032368,0.006025,5.368,118.7,0.000012,0.000171,0.001180,0.004921,0.000735,0.000179 +0.5997,1.93,0.028991,0.005651,5.124,118.7,-0.000005,0.000152,0.001191,0.005244,0.000840,0.000214 +0.6497,1.63,0.024436,0.005461,4.479,118.7,0.000005,0.000155,0.000909,0.005187,0.001006,0.000192 +0.6998,1.38,0.021369,0.005008,4.271,118.7,-0.000006,0.000110,0.001036,0.005256,0.001100,0.000254 +0.1438,16.54,0.511681,0.168126,3.05,137.6,0.000321,0.000253,0.011200,0.006054,0.000175,0.000126 +0.1967,10.44,0.287443,0.066729,4.317,137.6,0.000181,0.000342,0.006339,0.005892,0.000271,0.000127 +0.2479,7.36,0.181626,0.034224,5.307,137.6,0.000117,0.000317,0.004251,0.005763,0.000334,0.000133 +0.2987,5.35,0.120000,0.020784,5.77,137.6,0.000112,0.000269,0.002590,0.005461,0.000428,0.000117 +0.3490,4.26,0.086831,0.014101,6.148,137.6,0.000077,0.000223,0.002108,0.005365,0.000466,0.000129 +0.3993,3.39,0.063909,0.010398,6.149,137.6,0.000039,0.000193,0.001442,0.005150,0.000528,0.000115 +0.4494,2.97,0.051300,0.008533,6.025,137.6,0.000028,0.000202,0.001399,0.005228,0.000594,0.000142 +0.4995,2.57,0.041302,0.007328,5.639,137.6,0.000026,0.000180,0.001416,0.005194,0.000683,0.000177 +0.5496,2.10,0.033197,0.006169,5.379,137.6,0.000012,0.000181,0.001412,0.005046,0.000749,0.000214 +0.5996,1.97,0.029218,0.005870,4.983,137.6,0.000034,0.000165,0.001203,0.005287,0.000875,0.000216 +0.6497,1.81,0.025760,0.005344,4.826,137.6,0.000033,0.000146,0.001324,0.005467,0.000956,0.000280 +0.6998,1.27,0.021027,0.004685,4.485,137.6,0.000011,0.000126,0.000860,0.005170,0.001033,0.000211 +0.1438,16.52,0.516425,0.169777,3.046,147.9,0.000338,0.000316,0.011512,0.006110,0.000180,0.000130 +0.1967,10.38,0.288981,0.066816,4.312,147.9,0.000191,0.000343,0.006578,0.005922,0.000273,0.000132 +0.2479,7.35,0.182513,0.034661,5.269,147.9,0.000128,0.000339,0.004313,0.005793,0.000345,0.000135 +0.2987,5.30,0.120677,0.020534,5.876,147.9,0.000100,0.000270,0.002611,0.005490,0.000419,0.000118 +0.3490,4.23,0.087459,0.014184,6.17,147.9,0.000079,0.000237,0.002164,0.005403,0.000470,0.000133 +0.3993,3.44,0.065334,0.010607,6.177,147.9,0.000066,0.000220,0.001446,0.005266,0.000533,0.000116 +0.4494,2.91,0.051259,0.009104,5.636,147.9,0.000048,0.000192,0.001201,0.005227,0.000657,0.000122 +0.4995,2.57,0.041501,0.007368,5.644,147.9,0.000038,0.000177,0.001426,0.005220,0.000688,0.000178 +0.5496,2.06,0.033192,0.006260,5.307,147.9,0.000026,0.000160,0.001268,0.005046,0.000766,0.000192 +0.5996,1.99,0.029903,0.005895,5.072,147.9,0.000028,0.000170,0.001173,0.005411,0.000874,0.000211 +0.6497,1.71,0.025599,0.005539,4.627,147.9,0.000003,0.000154,0.001314,0.005431,0.001008,0.000278 +0.6998,1.45,0.022245,0.005113,4.349,147.9,-0.000005,0.000110,0.001209,0.005470,0.001115,0.000296 diff --git a/docs/API_FIELD_MAPPING.md b/docs/API_FIELD_MAPPING.md index 27325cc..1d86d12 100644 --- a/docs/API_FIELD_MAPPING.md +++ b/docs/API_FIELD_MAPPING.md @@ -167,7 +167,6 @@ RotorInputs(collective_rad, tilt_lon, tilt_lat, R_hub, |------------|-----------------|--------|------|-------| | `I_blade_flap_kgm2` | `.I_blade_flap_kgm2` | Y | f64 | Flap inertia [kg⋅m²] | | `omega_nr_rad_s` | `.omega_nr_rad_s` | Y | f64 | Natural frequency [rad/s] | -| **Method:** | `.hub_moment_factor(omega_rad_s)` | Y | f64 | Compute reduction factor | **PyClass Name:** `_dynbem.FlapProperties` diff --git a/docs/BEM_COMMON.md b/docs/BEM_COMMON.md index 9af490f..3794e23 100644 --- a/docs/BEM_COMMON.md +++ b/docs/BEM_COMMON.md @@ -561,27 +561,29 @@ hub-frame aero moment rotated to world frame (hub rolling/pitching moments); `M_spin` is the reaction torque about the hub axis; `Q_spin` is the scalar shaft torque. -When blade flapping is configured (Section 14a), the hub moments are reduced -before rotation to world frame: +When blade flapping is configured (Section 14a), the phase-correct 1/rev +flap solve sets the hub moments transmitted to the airframe and adds a +flapping-tilt contribution to the in-plane hub force before rotation to +world frame: -$$M_{x,\text{out}} = f_\text{hub}\,M_{x,\text{hub}}, \qquad M_{y,\text{out}} = f_\text{hub}\,M_{y,\text{hub}}$$ +$$[M_{x,\text{out}},\;M_{y,\text{out}}] = \tfrac{N_b}{2}\,I_b\,\Omega^2\,(\nu_\beta^2-1)\,[\beta_{1s},\;\beta_{1c}]$$ -Otherwise $M_{x,\text{out}} = M_{x,\text{hub}}$ (rigid blade, backward -compatible). +Otherwise $M_{x,\text{out}} = M_{x,\text{hub}}$ (rigid blade, no flap +solve). --- -## 14a. Quasi-Static Blade Flapping (Hub Moment Reduction) +## 14a. Quasi-Static Blade Flapping (Hub Moment + Flapback H-force) Source: `FlapProperties` in `dynbem_rs/src/rotor_definition.rs`, -`apply_flap_reduction` in `dynbem_rs/src/bem_common.rs`. +`apply_flap_dynamics` in `dynbem_rs/src/bem_common.rs`. -A flexible blade (or one with a flap hinge) absorbs most of the -aerodynamic pitching/rolling moment via out-of-plane deflection. Only the -fraction determined by the blade's flap frequency ratio reaches the -airframe hub. This is modelled as a quasi-static (algebraic) moment -reduction applied after the psi-loop integration and before output -assembly. +A flexible blade (or one with a flap hinge) responds to the 1/rev cyclic +aerodynamic moment with an out-of-plane flap `beta(psi)`. Only the +fraction set by the blade's centrifugal flap stiffness reaches the +airframe as a hub moment; the tilted disk also contributes an in-plane +hub force (the flapping-tilt "H-force"). Both are solved phase-correctly +after the psi-loop integration and before output assembly. ### Equivalent Spring-Hinge Model @@ -595,40 +597,63 @@ natural frequency (from root bending stiffness $K_\beta = I_b\,\omega_\text{NR}^ and $\Omega$ is the rotor speed. The "1" comes from centrifugal stiffening at the rotating speed. -### Hub Moment Reduction Factor +### Phase-Correct 1/rev Harmonic Solve -$$f_\text{hub} = \frac{\nu_\beta^2 - 1}{\nu_\beta^2}$$ +Near $\nu_\beta \approx 1$ the 1/rev forcing is near resonance, so the +flap response is set by aerodynamic damping and lags the forcing by +~90 deg. A magnitude-only scaling cannot reproduce this, so the flap DOF +is solved as a damped 1/rev harmonic balance ($' = d/d\psi$): -| Configuration | $`\omega_\text{NR}`$ | $`\nu_\beta`$ | $`f_\text{hub}`$ | Physical meaning | +$$I_b\,\Omega^2\,(\beta'' + d\,\beta' + \nu_\beta^2\,\beta) = M_\beta(\psi)$$ + +with aerodynamic flap damping (Lock number over 8), integrated over the +radial grid to allow taper: + +$$d = \frac{\gamma}{8} = \frac{\tfrac{1}{2}\,\rho\,a\,S_3}{I_b}, \qquad S_3 = \int c(r)\,r^3\,dr$$ + +($\Omega$ cancels). The azimuth-averaged hub moments from the sweep ARE +the 1/rev aero flap-moment harmonics, `mx_hub = N_b*M1s/2`, +`my_hub = N_b*M1c/2`, and at 1/rev ($\beta'' = -\beta$) the balance is the +linear system + +$$\begin{bmatrix} a & d \\ -d & a \end{bmatrix} \begin{bmatrix} \beta_{1c} \\ \beta_{1s} \end{bmatrix} = \frac{1}{I_b\,\Omega^2} \begin{bmatrix} M_{1c} \\ M_{1s} \end{bmatrix}, \qquad a = \nu_\beta^2 - 1$$ + +### Transmitted Hub Moment and Flapback H-force + +Transmitted hub moment (Johnson): + +$$[M_{x,\text{out}},\;M_{y,\text{out}}] = \tfrac{N_b}{2}\,I_b\,\Omega^2\,(\nu_\beta^2-1)\,[\beta_{1s},\;\beta_{1c}]$$ + +Flapping-tilt H-force (hub frame, $T$ = mean disk thrust), added to the +profile-drag H-force before assembly: + +$$F_{x,\text{flap}} = -\tfrac{1}{2}\,T\,\beta_{1c}, \qquad F_{y,\text{flap}} = +\tfrac{1}{2}\,T\,\beta_{1s}$$ + +| Configuration | $`\omega_\text{NR}`$ | $`\nu_\beta`$ | Transmitted moment | Flapback H-force | |---|---|---|---|---| -| Freely hinged (teetering) | 0 | 1.0 | 0 | No moment transfer to airframe | -| Typical hingeless | $`0.2{-}0.4\,\Omega`$ | 1.02-1.08 | 0.04-0.14 | Blade absorbs 86-96% of moment | -| Very stiff / rigid | $`\gg \Omega`$ | $`\gg 1`$ | $`\to 1`$ | Full moment transfer (legacy behaviour) | +| Freely hinged (teetering) | 0 | 1.0 | 0 (a=0) | present (pure 90 deg lag) | +| Typical hingeless | $`0.2{-}0.4\,\Omega`$ | 1.02-1.08 | small fraction of aero moment | present | +| Very stiff / rigid | $`\gg \Omega`$ | $`\gg 1`$ | $`\to`$ full aero moment | $`\to 0`$ | + +Note the flapback H-force is nonzero even for a freely hinged blade that +transmits **zero** hub moment: aerodynamic damping alone gives $\beta_1$ a +pure 90 deg lag. Its direction is validated (forward flight along +X gives +a rearward, flow-opposing force) by the `flapback` case in +`validation_rs/src/checks/h_force.rs`. -### Where the Reduction is Applied +### Where the Solve is Applied -The reduction acts **only on the output to the airframe**. The Pitt-Peters +The flap solve acts **only on the output to the airframe**. The Pitt-Peters inflow ODE (Section 10) still uses the *full* aerodynamic moments $`(M_{x,\text{hub}}, M_{y,\text{hub}})`$ because the wake responds to disk loading, not to what the airframe sees. Thrust and torque are unaffected. +The profile-drag H-force is computed in the sweep regardless of flapping; +absent a `flap:` section only the flapping-tilt term is omitted. -This is physically correct: the blade flaps to a new equilibrium angle -that redistributes the aerodynamic reaction between blade inertia -(centrifugal restoring) and hub structure (root bending), but the total -aerodynamic moment on the disk -- which drives the wake -- is unchanged. - -### Relation to Attitude Damping - -On a rigid-blade rotor, when the airframe pitches at rate $q$, the hub -moment is proportional to the rate (aerodynamic damping from the -asymmetric thrust distribution). With flapping, the blade absorbs most -of this moment -- the effective aerodynamic damping derivative seen by the -airframe is multiplied by $f_\text{hub}$. However, for an articulated or -flexible rotor the *cyclic inflow feedback* (Pitt-Peters $\lambda_c$, -$\lambda_s$) provides additional dynamic damping that is not captured -by this quasi-static factor. The two effects are complementary: the -flap reduces the instantaneous moment transfer while the inflow lag -provides the dynamic (rate-dependent) component. +This is physically correct: the blade flaps to redistribute the +aerodynamic reaction between blade inertia (centrifugal restoring) and +hub structure, but the total aerodynamic moment on the disk -- which +drives the wake -- is unchanged. --- diff --git a/docs/EMPIRICAL_VALIDATION.md b/docs/EMPIRICAL_VALIDATION.md index 28d4ae7..2460ef9 100644 --- a/docs/EMPIRICAL_VALIDATION.md +++ b/docs/EMPIRICAL_VALIDATION.md @@ -264,6 +264,32 @@ the TR 515 wind-tunnel data into rotor-axis (CT, CH, CY) figures (§5.3 pp. 87–89), which are easier to compare with BEM output than the airplane-axes CL/CD form in the original tables. +**Primary tabulated coefficients (H-force).** Beyond the §5.3 figures, +the report tabulates the rotor-axis coefficients directly. Printed pages +504-505 give a 48-row appendix (Advance Ratio, Shaft AoA, airplane +CL/CD/L-over-D, RPM, airplane roll/pitch/lateral coeffs, and Rotor CT / +CH / CY) for four RPM series (98.63, 118.7, 137.6, 147.9) across +mu ~ 0.14-0.70. These were extracted straight from the PDF text layer +(not a figure digitization, so HIGH confidence) into +[`Research/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs.md`](../Research/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs.md). +Rotor CH is the in-plane hub-force coefficient and is the primary +empirical anchor for the model's H-force. + +**H-force validation check.** +[`validation_rs/src/checks/harris_hforce_empirical.rs`](../validation_rs/src/checks/harris_hforce_empirical.rs) +(module `harris_hforce_empirical`) runs the quasi-static BEM against the +137.6 RPM series (12 points). Since the PCA-2 autorotates and blade pitch +is not tabulated, each point trims collective so the BEM matches the +tabulated Rotor CT, then compares the resulting Rotor CH. The blade is +modelled freely hinged so the flapping-tilt H-force is included. CT trims +exactly (<0.01%); CH agreement is best in the mid-range (mu 0.30-0.45, +0.5-18%) and degrades at the extremes (low-mu over-predicts ~80%, high-mu +under-predicts ~40-49%) -- the expected simple-BEM bias at high advance +ratio (reverse flow, radial flow, real polars). Per-case ceilings live in +`harris_hforce_empirical.csv`; re-baseline with +`REWRITE_EMPIRICAL_CSV=1 cargo run --release -p validation_rs -- harris_hforce`. + + **Rotor (Pitcairn PCA-2).** 4 blades, untwisted, mixed airfoils (symmetric outer 22-3/4 in chord + cambered inner 14-25/32 in). Diameter 45 ft. Tested at four nominal RPMs (98.6, 118.7, 137.6, 147.9), three pitch diff --git a/docs/MASTER_FIELD_MAPPING.md b/docs/MASTER_FIELD_MAPPING.md index b97b1b8..310b335 100644 --- a/docs/MASTER_FIELD_MAPPING.md +++ b/docs/MASTER_FIELD_MAPPING.md @@ -58,7 +58,6 @@ Comprehensive single-reference table of all Rust struct fields and their Python | | flap | .flap | ServoFlapGeometry | ✓ | r/o; nested | | **FlapProperties** | I_blade_flap_kgm2 | .I_blade_flap_kgm2 | f64 | ✓ | r/o | | | omega_nr_rad_s | .omega_nr_rad_s | f64 | ✓ | r/o | -| | *(method)* | .hub_moment_factor(ω) | f64 | ✓ | computed | | **RotorDefinition** | blade | .blade | BladeGeometry | ✓ | r/o; nested | | | airfoil | .airfoil | LinearPolarParameters | ✓ | r/o; nested | | | control | .control | Option | ✓ | r/o | diff --git a/dynbem/pyproject.toml b/dynbem/pyproject.toml index 7d9c539..aeb7b54 100644 --- a/dynbem/pyproject.toml +++ b/dynbem/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "dynbem" -version = "0.7.0" +version = "0.8.0" description = "Rotor aerodynamics: BEM, Pitt-Peters, Oye dynamic inflow, and VPM free-wake -- all in one library" requires-python = ">=3.10" license = "MIT" diff --git a/dynbem/python/dynbem/rotor_definition.py b/dynbem/python/dynbem/rotor_definition.py index 9a62806..4a38495 100644 --- a/dynbem/python/dynbem/rotor_definition.py +++ b/dynbem/python/dynbem/rotor_definition.py @@ -206,10 +206,6 @@ def _to_rust(self): omega_nr_rad_s=float(self.omega_nr_rad_s), ) - def hub_moment_factor(self, omega_rad_s): - """Fraction of aero hub moment that passes to airframe at given omega.""" - return self._to_rust().hub_moment_factor(omega_rad_s) - class InertiaProperties: """Rotor inertia data (Python-only, not used in Rust math).""" diff --git a/dynbem/src/wrappers.rs b/dynbem/src/wrappers.rs index a90230f..f11d609 100644 --- a/dynbem/src/wrappers.rs +++ b/dynbem/src/wrappers.rs @@ -470,11 +470,6 @@ impl PyFlapProperties { fn omega_nr_rad_s(&self) -> f64 { self.0.omega_nr_rad_s } - - /// Compute the hub moment reduction factor at a given rotor speed. - fn hub_moment_factor(&self, omega_rad_s: f64) -> f64 { - self.0.hub_moment_factor(omega_rad_s) - } } #[pyclass(name = "RotorDefinition", module = "dynbem._dynbem")] diff --git a/dynbem_rs/README.md b/dynbem_rs/README.md index c7c3f72..a3f457f 100644 --- a/dynbem_rs/README.md +++ b/dynbem_rs/README.md @@ -199,8 +199,8 @@ All quantities are in SI units in the world (NED) frame. | Field | Type | Description | |---|---|---| -| `F_world` | `Vec3` | Net aerodynamic force on the rotor hub, N, in NED world frame. In hover with the rotor level, thrust appears as a negative Z component (i.e. `F_world.0[2] < 0` lifts the aircraft). | -| `M_hub_world` | `Vec3` | Aerodynamic hub moment (pitching/rolling) about the hub centre, N·m, in NED world frame. Corresponds to `M_x`/`M_y` in standard helicopter notation. This is what the airframe sees after blade flap reduction (if `FlapProperties` is set). In NED, positive Mx is roll-right, positive My is pitch-up. | +| `F_world` | `Vec3` | Net aerodynamic force on the rotor hub, N, in NED world frame. In hover with the rotor level, thrust appears as a negative Z component (i.e. `F_world.0[2] < 0` lifts the aircraft). In edgewise/forward flight it also carries the in-plane hub force ("H-force") -- the profile/induced-drag term plus, when `FlapProperties` is set, the flapping-tilt term. | +| `M_hub_world` | `Vec3` | Aerodynamic hub moment (pitching/rolling) about the hub centre, N·m, in NED world frame. Corresponds to `M_x`/`M_y` in standard helicopter notation. This is what the airframe sees after the phase-correct 1/rev blade-flap solve (if `FlapProperties` is set), which reduces and rotates the transmitted moment. In NED, positive Mx is roll-right, positive My is pitch-up. | | `Q_spin` | `f64` | Aerodynamic reaction torque opposing rotation, N·m. Always positive for a powered rotor (torque opposes spin). Pass straight through to `crate::mechanical::step_omega` -- no sign flip needed. | | `M_spin` | `Vec3` | `Q_spin` expressed as a vector along the rotor spin axis in world frame, N·m. Direction is `R_hub * [0,0,1]` (the hub +Z axis in world coordinates). Convenience for rigid-body integrators that accept a 3-vector torque input. | @@ -297,7 +297,7 @@ are equivalent. +-- aero_model.rs AeroModel trait + RotorStateExt trait + IntegrationMethod +-- bem_common.rs RadialGrid, PolarTable, kinematics(), element_force(), | run_psi_loop, assemble_result(), - | apply_flap_reduction() + | apply_flap_dynamics() +-- common.rs numerical floors (EPS_*), vrs_lambda1 VRS polynomial +-- cyclic.rs swashplate -> theta_1c, theta_1s mapping +-- mechanical.rs omega_derivative(), step_omega() -- single canonical diff --git a/dynbem_rs/src/bem_common.rs b/dynbem_rs/src/bem_common.rs index 11f11d6..1ffb9ca 100644 --- a/dynbem_rs/src/bem_common.rs +++ b/dynbem_rs/src/bem_common.rs @@ -232,36 +232,161 @@ pub fn vrs_regime(t_total: f64, v_climb: f64, rho: f64, area: f64) -> VrsRegime // world-frame outputs for every model. // --------------------------------------------------------------------------- -/// Apply quasi-static flap reduction to hub moments. +/// Result of the 1/rev blade-flapping solve. /// -/// Returns (mx_hub_reduced, my_hub_reduced). If `flap` is None, moments -/// pass through unchanged (rigid-blade assumption). +/// `mx_out`/`my_out` are the hub moments that reach the airframe (the flap +/// deflection redistributes the aerodynamic reaction between blade +/// centrifugal restoring and hub structure). `dfx_hub`/`dfy_hub` are the +/// additional in-plane hub force (the flapping-tilt contribution to +/// H-force): as the blade flaps by beta(psi) the thrust vector tilts into +/// the disk plane. Both are in the hub frame and are added to the +/// profile-drag H-force before `assemble_result`. +pub struct FlapOutput { + pub mx_out: f64, + pub my_out: f64, + pub dfx_hub: f64, + pub dfy_hub: f64, +} + +/// Phase-correct quasi-static blade-flapping solve (1/rev harmonic balance +/// with aerodynamic flap damping). +/// +/// The flap DOF responds to the RAW aerodynamic hub moment (the caller must +/// pass the un-reduced `mx_hub`/`my_hub`; the inflow ODE still uses those +/// full moments because the wake responds to disk loading, not to what the +/// airframe sees). The solve returns both the reduced hub moments and the +/// flapping-tilt H-force. +/// +/// ## Convention (see AGENTS.md "blade flapping" section) +/// +/// - Azimuth psi=0 at +X, CCW from above; blade flap `beta > 0` = tip up +/// (toward -z, the thrust direction). `beta(psi) = beta_0 + beta_1c cos +/// psi + beta_1s sin psi`. +/// - The sweep's azimuth-averaged hub moments are the 1/rev aerodynamic +/// flap-moment harmonics: `mx_hub = N_b M1s / 2`, `my_hub = N_b M1c / 2`. +/// - Flap equation (per rev, ' = d/dpsi): +/// `I_b Omega^2 (beta'' + d beta' + nu^2 beta) = M_beta(psi)` with +/// aerodynamic damping `d = gamma/8`, Lock number `gamma = rho a c R^4 / +/// I_b`. Computed here from the actual radial geometry to allow taper: +/// `d = 0.5 rho a S3 / I_b`, `S3 = sum c(r) r^3 dr` (Omega cancels). +/// - 1/rev balance gives the damped transfer +/// `[[a, d],[-d, a]] [beta_1c; beta_1s] = (1/(I_b Omega^2)) [M1c; M1s]`, +/// `a = nu^2 - 1`. At `nu ~ 1` the damping dominates and produces the +/// classical ~90 deg flap lag (flapback) rather than an in-phase response. +/// - Transmitted hub moment (Johnson): `M_hub = (N_b/2) I_b Omega^2 (nu^2-1) +/// beta_1`. With no damping this reduces to the full aero moment. +/// - Flapping-tilt H-force: `Fx = - = -T beta_1c/2`, +/// `Fy = + = +T beta_1s/2` (T = mean disk thrust). +/// +/// The H-force sign is pinned by `validation_rs/src/checks/h_force.rs` +/// (forward flight with collective must give a rearward, flow-opposing +/// force that adds to the profile-drag term). #[inline] -pub fn apply_flap_reduction( +pub fn apply_flap_dynamics( + t_total: f64, mx_hub: f64, my_hub: f64, flap: Option<&crate::rotor_definition::FlapProperties>, + grid: &RadialGrid, + cl_alpha: f64, + rho: f64, + n_b: usize, omega: f64, -) -> (f64, f64) { - match flap { - Some(fp) => { - let f = fp.hub_moment_factor(omega); - (mx_hub * f, my_hub * f) +) -> FlapOutput { + let fp = match flap { + Some(fp) => fp, + None => { + return FlapOutput { + mx_out: mx_hub, + my_out: my_hub, + dfx_hub: 0.0, + dfy_hub: 0.0, + } } - None => (mx_hub, my_hub), + }; + + let i_b = fp.I_blade_flap_kgm2; + if !(i_b > 0.0) || omega.abs() < 1e-6 { + // Degenerate: no inertia or not spinning -> no meaningful flap solve. + return FlapOutput { + mx_out: mx_hub, + my_out: my_hub, + dfx_hub: 0.0, + dfy_hub: 0.0, + }; + } + + let nu2 = fp.nu_beta_sq(omega); + let a_stiff = nu2 - 1.0; + + // Aerodynamic flap damping (nondimensional, = gamma/8). Integrated over + // the radial grid so taper is handled: C_beta = 0.5 rho a Omega * S3, + // S3 = sum c(r) r^3 dr, and gamma/8 = C_beta/(I_b Omega) = 0.5 rho a + // S3 / I_b (Omega cancels). + let mut s3 = 0.0; + for i in 0..grid.n_elements { + let r = grid.r_mid[i]; + s3 += grid.chord[i] * r * r * r * grid.dr; + } + let d_damp = 0.5 * rho * cl_alpha * s3 / i_b; + + // 1/rev aerodynamic flap-moment harmonics (per blade). + let n_bf = n_b as f64; + let m1c = 2.0 * my_hub / n_bf; + let m1s = 2.0 * mx_hub / n_bf; + + // Damped 1/rev flap response: + // [[a, d],[-d, a]] [b1c; b1s] = (1/(I_b Omega^2)) [M1c; M1s] + // => [b1c; b1s] = (1/(I_b Omega^2 D)) [[a, -d],[d, a]] [M1c; M1s] + let om2 = omega * omega; + let det = a_stiff * a_stiff + d_damp * d_damp; + let inv = 1.0 / (i_b * om2 * det); + let b1c = inv * (a_stiff * m1c - d_damp * m1s); + let b1s = inv * (d_damp * m1c + a_stiff * m1s); + + // Transmitted hub moment (Johnson): M_hub = (N_b/2) I_b Omega^2 (nu^2-1) beta1. + let k = 0.5 * n_bf * i_b * om2 * a_stiff; + let my_out = k * b1c; + let mx_out = k * b1s; + + // Flapping-tilt H-force from the thrust vector tilting with beta(psi). + let dfx_hub = -0.5 * t_total * b1c; + let dfy_hub = 0.5 * t_total * b1s; + + FlapOutput { + mx_out, + my_out, + dfx_hub, + dfy_hub, } } +/// +/// `fx_hub, fy_hub` are the net in-plane (H-force) components in the hub +/// frame -- the world-frame reaction of the blades' tangential aerodynamic +/// force, summed over the disk (see `SweepCtx::run`'s `fx_acc`/`fy_acc`). +/// This is zero whenever the tangential loading is azimuth-independent +/// (pure axial flow / hover) and grows with edgewise flow (`v_edge`), +/// vanishing again as the hub axis re-aligns with the relative wind. +/// Pass `0.0, 0.0` for callers that don't yet compute it (e.g. VPM). +/// +/// The caller is expected to have already folded in the flapping-tilt +/// contribution to H-force (`FlapOutput::dfx_hub`/`dfy_hub` from +/// `apply_flap_dynamics`) when a `FlapProperties` is configured; the +/// profile/induced-drag part comes from `SweepCtx::run`. #[inline] pub fn assemble_result( t_total: f64, q_total: f64, mx_hub: f64, my_hub: f64, + fx_hub: f64, + fy_hub: f64, hub_axis: Vec3, r_hub: &Mat3, ) -> AeroResult { - let f_world = hub_axis * (-t_total); + let f_hub = Vec3::new(fx_hub, fy_hub, -t_total); + let f_world = *r_hub * f_hub; let mxyz_hub = Vec3::new(mx_hub, my_hub, 0.0); let m_orbital = *r_hub * mxyz_hub; let m_spin = hub_axis * q_total; @@ -361,17 +486,30 @@ pub struct SweepCtx<'a, P: Polar> { impl<'a, P: Polar> SweepCtx<'a, P> { /// Run one full psi x radial sweep with the given kernel. Returns the - /// azimuth-averaged (T, Q, Mx_hub, My_hub) over the rotor disk. + /// azimuth-averaged (T, Q, Mx_hub, My_hub, Fx_hub, Fy_hub) over the rotor + /// disk. + /// + /// `Fx_hub`/`Fy_hub` are the net in-plane hub force (H-force): each + /// element's tangential aerodynamic force is `dFt = dQ / r` (the same + /// force whose moment produces `dQ`; recovered by undoing the `* r` + /// weighting rather than recomputing it), projected onto the fixed hub + /// x/y axes with the same `(sin psi, cos psi)` pairing used everywhere + /// else in this sweep (`v_t_extra = v_in_hub_x*sin psi + v_in_hub_y*cos + /// psi`, `Mx_hub`/`My_hub` below). This is zero for azimuth-independent + /// (pure axial) loading and grows with edgewise flow -- see + /// `assemble_result`. /// /// `self.omega > 0` is assumed (caller filters out the not-spinning case /// before invoking). Reverse-flow region (`v_t <= 0`) is skipped /// per-element. #[inline(always)] - pub fn run(&self, kernel: &mut K) -> (f64, f64, f64, f64) { + pub fn run(&self, kernel: &mut K) -> (f64, f64, f64, f64, f64, f64) { let mut t_acc = 0.0; let mut q_acc = 0.0; let mut mx_acc = 0.0; let mut my_acc = 0.0; + let mut fx_acc = 0.0; + let mut fy_acc = 0.0; let inv_n_psi = self.n_psi_inv; let grid = self.grid; let n_r = grid.n_elements; @@ -410,6 +548,9 @@ impl<'a, P: Polar> SweepCtx<'a, P> { t_acc += dt; q_acc += dq; rdt_sum += r * dt; + let d_ft = dq / r; + fx_acc += d_ft * sin_psi; + fy_acc += d_ft * cos_psi; } mx_acc += rdt_sum * sin_psi; my_acc += rdt_sum * cos_psi; @@ -419,6 +560,8 @@ impl<'a, P: Polar> SweepCtx<'a, P> { q_acc * inv_n_psi, mx_acc * inv_n_psi, my_acc * inv_n_psi, + fx_acc * inv_n_psi, + fy_acc * inv_n_psi, ) } } diff --git a/dynbem_rs/src/oye.rs b/dynbem_rs/src/oye.rs index 18da1d3..5a78932 100644 --- a/dynbem_rs/src/oye.rs +++ b/dynbem_rs/src/oye.rs @@ -6,7 +6,7 @@ use std::f64::consts::PI; use crate::aero_io::{AeroResult, RotorInputs}; use crate::aero_model::{AeroModel, RotorStateExt}; use crate::bem_common::{ - apply_flap_reduction, assemble_result, build_psi_trig_table, element_force, kinematics, + apply_flap_dynamics, assemble_result, build_psi_trig_table, element_force, kinematics, v_mass_flow_disk, vrs_regime, ElementCtx, PsiKernel, RadialGrid, SweepCtx, }; use crate::common::{ @@ -281,32 +281,33 @@ impl AeroModel for OyeBEMModel

{ 0.0 }; let mut dt_avg = vec![0.0; n]; - let (t_total, q_total, mx_hub, my_hub) = if omega_r > EPS_OMEGA_R && omega > 1.0 { - let mut kernel = OyeKernel { - lambda_climb, - w: state.w_slice(), - dt_avg: &mut dt_avg, + let (t_total, q_total, mx_hub, my_hub, fx_hub, fy_hub) = + if omega_r > EPS_OMEGA_R && omega > 1.0 { + let mut kernel = OyeKernel { + lambda_climb, + w: state.w_slice(), + dt_avg: &mut dt_avg, + }; + let sweep = SweepCtx { + grid: &self.grid, + polar: &self.polar, + col: loop_collective, + omega, + omega_r, + rho, + n_b: blade.n_blades, + n_psi: self.n_psi_elements, + n_psi_inv: 1.0 / (self.n_psi_elements as f64), + psi_trig: &self.psi_trig, + v_in_hub_x: v_inplane_hub[0], + v_in_hub_y: v_inplane_hub[1], + theta_1c: loop_theta_1c, + theta_1s: loop_theta_1s, + }; + sweep.run(&mut kernel) + } else { + (0.0, 0.0, 0.0, 0.0, 0.0, 0.0) }; - let sweep = SweepCtx { - grid: &self.grid, - polar: &self.polar, - col: loop_collective, - omega, - omega_r, - rho, - n_b: blade.n_blades, - n_psi: self.n_psi_elements, - n_psi_inv: 1.0 / (self.n_psi_elements as f64), - psi_trig: &self.psi_trig, - v_in_hub_x: v_inplane_hub[0], - v_in_hub_y: v_inplane_hub[1], - theta_1c: loop_theta_1c, - theta_1s: loop_theta_1s, - }; - sweep.run(&mut kernel) - } else { - (0.0, 0.0, 0.0, 0.0) - }; let w_mean: f64 = if n > 0 { state.w_slice().iter().sum::() / (n as f64) @@ -374,9 +375,27 @@ impl AeroModel for OyeBEMModel

{ n, ); - let (mx_out, my_out) = - apply_flap_reduction(mx_hub, my_hub, self.defn.flap.as_ref(), inputs.omega_rad_s); - let result = assemble_result(t_total, q_total, mx_out, my_out, hub_axis, &inputs.R_hub); + let flap = apply_flap_dynamics( + t_total, + mx_hub, + my_hub, + self.defn.flap.as_ref(), + &self.grid, + self.defn.airfoil.CL_alpha_per_rad, + rho, + self.defn.blade.n_blades, + omega, + ); + let result = assemble_result( + t_total, + q_total, + flap.mx_out, + flap.my_out, + fx_hub + flap.dfx_hub, + fy_hub + flap.dfy_hub, + hub_axis, + &inputs.R_hub, + ); let derivative = OyeRotorState { n_elements: n, W_int: d_w_int, diff --git a/dynbem_rs/src/pitt_peters.rs b/dynbem_rs/src/pitt_peters.rs index 51da476..04bdc78 100644 --- a/dynbem_rs/src/pitt_peters.rs +++ b/dynbem_rs/src/pitt_peters.rs @@ -74,10 +74,10 @@ // the MU_T_FLOOR clamp). Recomputing v_mf from mu_t_eff would re-impose // the MU_T_FLOOR clamp -- inflating v_mf whenever v_mf/omega_r < // MU_T_FLOOR -- and corrupt the coefficients (see the inline note). -// F. Quasi-static blade flap reduction (`apply_flap_reduction`) scales the -// hub moments the AIRFRAME sees, but the inflow ODE uses the FULL -// aerodynamic moments (the wake responds to disk loading, not to what -// the flapping blade passes to the hub). +// F. Quasi-static blade flap dynamics (`apply_flap_dynamics`) set the +// hub moments the AIRFRAME sees (and the flapping-tilt H-force), but +// the inflow ODE uses the FULL aerodynamic moments (the wake responds +// to disk loading, not to what the flapping blade passes to the hub). // G. Servo-flap feathering pre-pass (Kaman path): in ServoFlap mode the // swashplate collective AND cyclic are reinterpreted as flap commands // and the feathering response replaces the direct blade-pitch path. @@ -87,7 +87,7 @@ use std::f64::consts::PI; use crate::aero_io::{AeroResult, RotorInputs}; use crate::aero_model::{AeroModel, RotorStateExt}; use crate::bem_common::{ - apply_flap_reduction, assemble_result, build_psi_trig_table, element_force, kinematics, + apply_flap_dynamics, assemble_result, build_psi_trig_table, element_force, kinematics, v_mass_flow_disk, vrs_regime, ElementCtx, PsiKernel, RadialGrid, SweepCtx, }; use crate::common::{vrs_lambda1, EPS_DENOM, EPS_OMEGA_R, MU_T_FLOOR}; @@ -251,7 +251,8 @@ impl AeroModel for PittPetersModel

{ // ------------------------------------------------------------------ // Blade element forces // ------------------------------------------------------------------ - let (mut t_total, mut q_total, mut mx_hub, mut my_hub) = (0.0, 0.0, 0.0, 0.0); + let (mut t_total, mut q_total, mut mx_hub, mut my_hub, mut fx_hub, mut fy_hub) = + (0.0, 0.0, 0.0, 0.0, 0.0, 0.0); if omega_r > EPS_OMEGA_R && omega > 1.0 { let mut kernel = PpKernel { lambda_total, @@ -275,11 +276,13 @@ impl AeroModel for PittPetersModel

{ theta_1c: loop_theta_1c, theta_1s: loop_theta_1s, }; - let (t, q, mx, my) = sweep.run(&mut kernel); + let (t, q, mx, my, fx, fy) = sweep.run(&mut kernel); t_total = t; q_total = q; mx_hub = mx; my_hub = my; + fx_hub = fx; + fy_hub = fy; } // ------------------------------------------------------------------ @@ -396,13 +399,33 @@ impl AeroModel for PittPetersModel

{ let d_lam_c = beta_c * d_lam_c_wind - beta_s * d_lam_s_wind; let d_lam_s = beta_s * d_lam_c_wind + beta_c * d_lam_s_wind; - // Outputs -- NON-STANDARD (flag F): apply quasi-static flap reduction - // to the hub moments before assembling the world-frame result. The + // Outputs -- NON-STANDARD (flag F): apply the 1/rev blade-flapping + // solve to the hub moments before assembling the world-frame result, + // and fold its flapping-tilt H-force into the in-plane hub force. The // inflow ODE above uses the FULL aerodynamic moments (they drive the // wake), but the airframe only sees the fraction that passes through - // the blade's flap stiffness. - let (mx_out, my_out) = apply_flap_reduction(mx_hub, my_hub, self.defn.flap.as_ref(), omega); - let result = assemble_result(t_total, q_total, mx_out, my_out, hub_axis, &inputs.R_hub); + // the blade's flap dynamics. + let flap = apply_flap_dynamics( + t_total, + mx_hub, + my_hub, + self.defn.flap.as_ref(), + &self.grid, + self.defn.airfoil.CL_alpha_per_rad, + rho, + blade.n_blades, + omega, + ); + let result = assemble_result( + t_total, + q_total, + flap.mx_out, + flap.my_out, + fx_hub + flap.dfx_hub, + fy_hub + flap.dfy_hub, + hub_axis, + &inputs.R_hub, + ); let derivative = PittPetersRotorState { lambda_0: d_lam0, lambda_c: d_lam_c, diff --git a/dynbem_rs/src/quasi_static_bem.rs b/dynbem_rs/src/quasi_static_bem.rs index 507ec81..fbb0ea0 100644 --- a/dynbem_rs/src/quasi_static_bem.rs +++ b/dynbem_rs/src/quasi_static_bem.rs @@ -6,7 +6,7 @@ use std::f64::consts::PI; use crate::aero_io::{AeroResult, RotorInputs}; use crate::aero_model::{AeroModel, RotorStateExt}; use crate::bem_common::{ - apply_flap_reduction, assemble_result, build_psi_trig_table, kinematics, ElementCtx, PsiKernel, + apply_flap_dynamics, assemble_result, build_psi_trig_table, kinematics, ElementCtx, PsiKernel, RadialGrid, SweepCtx, }; use crate::common::{EPS_DENOM, EPS_OMEGA_R, MIN_LOSS_FACTOR}; @@ -738,6 +738,8 @@ impl AeroModel for QuasiStaticBEM

{ let mut q_total: f64 = 0.0; let mut mx_hub: f64 = 0.0; let mut my_hub: f64 = 0.0; + let mut fx_hub: f64 = 0.0; + let mut fy_hub: f64 = 0.0; if (mu > 0.01 || has_cyclic || has_feathering_cyclic) && omega > 1.0 { let mut kernel = BemKernel { @@ -762,14 +764,18 @@ impl AeroModel for QuasiStaticBEM

{ theta_1c: loop_theta_1c, theta_1s: loop_theta_1s, }; - let (t, q, mx, my) = sweep.run(&mut kernel); + let (t, q, mx, my, fx, fy) = sweep.run(&mut kernel); t_total = t; q_total = q; mx_hub = mx; my_hub = my; + fx_hub = fx; + fy_hub = fy; } else { // Axial: try wind-turbine windmill solver first when v_climb < 0, // fall back to helicopter quadratic per element. + // (Fx_hub/Fy_hub stay 0: azimuth-symmetric loading has no net + // in-plane force.) for i_r in 0..n { let r = r_arr[i_r]; let geom = BEMElementGeometry::new( @@ -796,9 +802,27 @@ impl AeroModel for QuasiStaticBEM

{ } } - let (mx_out, my_out) = - apply_flap_reduction(mx_hub, my_hub, self.defn.flap.as_ref(), inputs.omega_rad_s); - let result = assemble_result(t_total, q_total, mx_out, my_out, hub_axis, &inputs.R_hub); + let flap = apply_flap_dynamics( + t_total, + mx_hub, + my_hub, + self.defn.flap.as_ref(), + grid, + self.defn.airfoil.CL_alpha_per_rad, + rho, + n_blades, + omega, + ); + let result = assemble_result( + t_total, + q_total, + flap.mx_out, + flap.my_out, + fx_hub + flap.dfx_hub, + fy_hub + flap.dfy_hub, + hub_axis, + &inputs.R_hub, + ); let derivative = QuasiStaticRotorState; // suppress unused warning (use_tip_loss read inside windmill helper). diff --git a/dynbem_rs/src/rotor_definition.rs b/dynbem_rs/src/rotor_definition.rs index 12ec269..5135a4a 100644 --- a/dynbem_rs/src/rotor_definition.rs +++ b/dynbem_rs/src/rotor_definition.rs @@ -233,19 +233,6 @@ impl FlapProperties { } 1.0 + (self.omega_nr_rad_s / omega).powi(2) } - - /// Fraction of the aerodynamic hub moment that passes to the airframe. - /// - /// factor = (nu_beta^2 - 1) / nu_beta^2 - /// - /// - Freely hinged (omega_NR=0): nu^2=1, factor=0 (no moment transfer) - /// - Rigid blade (omega_NR>>Omega): factor->1 (full moment transfer) - /// - Typical hingeless: factor ~ 0.05-0.15 - #[inline] - pub fn hub_moment_factor(&self, omega: f64) -> f64 { - let nu2 = self.nu_beta_sq(omega); - (nu2 - 1.0) / nu2 - } } /// How blade pitch is actuated. diff --git a/dynbem_rs/src/vpm/mod.rs b/dynbem_rs/src/vpm/mod.rs index e0ad182..78a263f 100644 --- a/dynbem_rs/src/vpm/mod.rs +++ b/dynbem_rs/src/vpm/mod.rs @@ -546,11 +546,16 @@ impl AeroModel for VpmRotor

{ let (fc, kin) = self.flight_condition(inputs); // One sub-step per call; dt IS the sub-step duration. let (res, out_state) = self.march_window(&fc, Some(state), state.psi, dt, 1, 1); + // TODO: VPM's free-wake march doesn't yet accumulate the in-plane + // hub force (H-force) the BEM-family models compute in + // bem_common::SweepCtx::run -- pass 0.0, 0.0 until that's added here. let result = assemble_result( res.thrust, res.torque, res.mx_hub, res.my_hub, + 0.0, + 0.0, kin.hub_axis, &inputs.R_hub, ); diff --git a/tests/test_flap_hinge.py b/tests/test_flap_hinge.py index 8712602..7e8f63a 100644 --- a/tests/test_flap_hinge.py +++ b/tests/test_flap_hinge.py @@ -1,13 +1,12 @@ """Tests for quasi-static blade flapping (hub moment reduction). Covers: - 1. FlapProperties hub_moment_factor correctness at known omega. - 2. With flap configured, hub moments are reduced vs rigid blade. - 3. Thrust is unchanged (flapping only affects moments, not axial force). - 4. Zero omega_nr (freely hinged) gives near-zero hub moment. - 5. Large omega_nr (stiff blade) gives nearly unchanged hub moment. - 6. YAML round-trip for flap section. - 7. All three models (pitt_peters, oye, bem) support flap. + 1. With flap configured, hub moments are reduced vs rigid blade. + 2. Thrust is unchanged (flapping only affects moments, not axial force). + 3. Zero omega_nr (freely hinged) gives near-zero hub moment. + 4. Large omega_nr (stiff blade) gives nearly unchanged hub moment. + 5. YAML round-trip for flap section. + 6. All three models (pitt_peters, oye, bem) support flap. """ import math @@ -53,33 +52,6 @@ def thrust_N(result): return -float(result.F_world[2]) -# --------------------------------------------------------------------------- -# 1. FlapProperties.hub_moment_factor correctness -# --------------------------------------------------------------------------- - -def test_hub_moment_factor_math(): - """Verify the reduction factor formula directly.""" - fp = _RustFlapProperties(I_blade_flap_kgm2=0.01, omega_nr_rad_s=10.0) - omega = 30.0 - # nu_beta^2 = 1 + (10/30)^2 = 1 + 1/9 = 10/9 - # factor = (10/9 - 1) / (10/9) = (1/9) / (10/9) = 1/10 - expected = 1.0 / 10.0 - assert fp.hub_moment_factor(omega) == pytest.approx(expected, rel=1e-10) - - -def test_hub_moment_factor_zero_stiffness(): - """Freely-hinged blade: factor should be 0 (no moment transfer).""" - fp = _RustFlapProperties(I_blade_flap_kgm2=0.01, omega_nr_rad_s=0.0) - assert fp.hub_moment_factor(30.0) == pytest.approx(0.0, abs=1e-15) - - -def test_hub_moment_factor_stiff(): - """Very stiff blade: factor approaches 1.""" - fp = _RustFlapProperties(I_blade_flap_kgm2=0.01, omega_nr_rad_s=1000.0) - # nu^2 = 1 + (1000/30)^2 ~ 1112; factor = 1111/1112 ~ 0.9991 - f = fp.hub_moment_factor(30.0) - assert f > 0.99 - assert f < 1.0 # --------------------------------------------------------------------------- @@ -245,12 +217,19 @@ def test_yaml_no_flap_section(): # --------------------------------------------------------------------------- def test_moment_reduction_matches_factor(): - """The ratio of flap/rigid moments should match hub_moment_factor.""" + """Phase-correct 1/rev flap solve reduces and rotates the hub moment. + + The model no longer scales the transmitted hub moment by the scalar + hub_moment_factor = (nu^2-1)/nu^2. The phase-correct harmonic solve adds + aerodynamic flap damping, which (a) reduces the transmitted-moment + magnitude below the rigid-blade value and (b) rotates the moment between + axes (the ~90 deg flap lag) -- something a scalar factor cannot do. So + the flap moment is NOT a simple positive multiple of the rigid moment. + """ + import math + omega = 25.0 omega_nr = 8.0 - # factor = (omega_nr/omega)^2 / (1 + (omega_nr/omega)^2) - nu2 = 1.0 + (omega_nr / omega) ** 2 - expected_factor = (nu2 - 1.0) / nu2 flap = FlapProperties(I_blade_flap_kgm2=0.01, omega_nr_rad_s=omega_nr) defn_rigid = make_defn(flap=None) @@ -265,8 +244,19 @@ def test_moment_reduction_matches_factor(): res_flap, _ = model_flap.compute_forces(inputs, model_flap.initial_rotor_state()) mx_rigid = float(res_rigid.m_hub_world[0]) + my_rigid = float(res_rigid.m_hub_world[1]) mx_flap = float(res_flap.m_hub_world[0]) + my_flap = float(res_flap.m_hub_world[1]) + + mag_rigid = math.hypot(mx_rigid, my_rigid) + mag_flap = math.hypot(mx_flap, my_flap) + + # Flap must strictly reduce the transmitted-moment magnitude. + assert mag_rigid > 1e-10 + assert mag_flap < mag_rigid - if abs(mx_rigid) > 1e-10: - actual_ratio = mx_flap / mx_rigid - assert actual_ratio == pytest.approx(expected_factor, rel=1e-6) + # The transmitted moment is rotated relative to the rigid one (phase lag), + # so the two vectors are not parallel: the cross product is non-negligible + # relative to the product of magnitudes. + cross = mx_rigid * my_flap - my_rigid * mx_flap + assert abs(cross) > 0.05 * mag_rigid * mag_flap diff --git a/validation_rs/src/checks/h_force.rs b/validation_rs/src/checks/h_force.rs new file mode 100644 index 0000000..e28796d --- /dev/null +++ b/validation_rs/src/checks/h_force.rs @@ -0,0 +1,307 @@ +// Directional check: the rotor's net in-plane hub force ("H-force") under +// a pure crosswind, with the disk held level and collective at ~0 (no +// thrust). See AGENTS.md "Rotor rotation direction" for the hub-frame +// convention and dynbem_rs/src/bem_common.rs::assemble_result for how +// Fx_hub/Fy_hub are assembled into F_world. +// +// Physical claim under test: a horizontal rotor disk is not "invisible" to +// a crosswind just because collective is zero. Blade profile drag varies +// with azimuth (advancing/retreating asymmetry from the edgewise flow), and +// that asymmetric tangential loading sums to a net in-plane force that +// points along the wind -- the same mechanism that would make a freely +// flying rotor drift downwind. This is separate from (and in addition to) +// the hub *moment* (Mx_hub/My_hub) already checked in cyclic_sign.rs. + +use crate::helpers::*; +use crate::report::Report; +use dynbem_rs::aero_io::{Mat3, RotorInputs, Vec3}; +use dynbem_rs::aero_model::AeroModel; +use dynbem_rs::quasi_static_bem::QuasiStaticBEM; + +pub fn check_h_force(r: &mut Report) { + r.begin_module( + "h_force", + "Directional: level-disk crosswind produces an in-plane hub force along the wind", + ); + + let defn = theory_rotor(8, 0.0); + let polar = theory_polar(); + let rotor = QuasiStaticBEM::build(defn, 36, polar); + + // Reasonable operating RPM (same as the rest of the theory suite), disk + // perfectly level (R_hub = identity), collective ~0 so there is no + // vertical thrust -- isolates the in-plane term. + let omega = OMEGA; + let wind_speed = 8.0; // m/s crosswind + let wind = Vec3::new(0.0, wind_speed, 0.0); // pure +Y (East) crosswind + + let inputs = RotorInputs { + collective_rad: 0.0, + tilt_lon: 0.0, + tilt_lat: 0.0, + R_hub: Mat3::eye(), + v_hub_world: Vec3::zero(), + wind_world: wind, + rho_kg_m3: RHO, + omega_rad_s: omega, + }; + + let (result, _) = rotor.compute_forces(&inputs, &rotor.initial_state()); + let f = result.F_world.0; + + r.info("level_disk_crosswind", "thrust_z_N", f[2], f64::NAN); + r.info("level_disk_crosswind", "f_north_N", f[0], f64::NAN); + r.info("level_disk_crosswind", "f_east_N", f[1], f64::NAN); + + // Collective ~0 -> no vertical thrust. + r.assert_bool( + "level_disk_crosswind", + "no_vertical_thrust", + f[2], + 0.0, + f[2].abs() < 1.0, + &format!( + "collective=0 should give ~zero vertical force, got F_z={:.3} N", + f[2] + ), + ); + + // The wind is purely +Y (East): the horizontal force should point the + // same way (downwind), not across or against it. + r.assert_bool( + "level_disk_crosswind", + "force_along_wind", + f[1], + 0.0, + f[1] > 0.05, + &format!( + "crosswind +Y should push the disk +Y (downwind), got F_east={:.3} N", + f[1] + ), + ); + r.assert_bool( + "level_disk_crosswind", + "no_cross_axis_force", + f[0], + 0.0, + f[0].abs() < 0.1 * f[1].abs().max(1.0), + &format!( + "pure +Y crosswind should not push North/South, got F_north={:.3} N (F_east={:.3} N)", + f[0], f[1] + ), + ); + + // Reversing the wind should reverse the force. + let inputs_rev = RotorInputs { + wind_world: Vec3::new(0.0, -wind_speed, 0.0), + ..inputs.clone() + }; + let (result_rev, _) = rotor.compute_forces(&inputs_rev, &rotor.initial_state()); + let f_rev = result_rev.F_world.0; + r.info("reversed_crosswind", "f_east_N", f_rev[1], f64::NAN); + r.assert_bool( + "reversed_crosswind", + "force_flips_with_wind", + f_rev[1], + 0.0, + f_rev[1] < -0.05, + &format!( + "reversing wind to -Y should reverse the force to F_east<0, got {:.3} N", + f_rev[1] + ), + ); + + // No wind at all (hover, zero collective) -> no in-plane force either. + let inputs_calm = RotorInputs { + wind_world: Vec3::zero(), + ..inputs.clone() + }; + let (result_calm, _) = rotor.compute_forces(&inputs_calm, &rotor.initial_state()); + let f_calm = result_calm.F_world.0; + r.info("no_wind", "f_east_N", f_calm[1], f64::NAN); + r.assert_bool( + "no_wind", + "no_force_without_wind", + f_calm[1], + 0.0, + f_calm[1].abs() < 0.5, + &format!( + "zero wind + zero collective should give ~zero in-plane force, got F_east={:.3} N", + f_calm[1] + ), + ); + + // --- disk tilting toward the wind should shrink the H-force to zero --- + // + // Roll the disk about the North (X) axis by phi: hub_axis rotates from + // [0,0,1] (level) toward [0,1,0] (pointing straight into the +Y wind). + // As phi -> 90 deg the relative wind becomes purely axial (v_edge -> 0), + // and the in-plane hub force -- which only exists because of the + // azimuthal (edgewise-flow) asymmetry -- should vanish with it. + fn hub_force_mag(rotor: &QuasiStaticBEM, phi: f64) -> f64 { + let (s, c) = phi.sin_cos(); + let r_hub = Mat3([[1.0, 0.0, 0.0], [0.0, c, s], [0.0, -s, c]]); + let inputs = RotorInputs { + collective_rad: 0.0, + tilt_lon: 0.0, + tilt_lat: 0.0, + R_hub: r_hub, + v_hub_world: Vec3::zero(), + wind_world: Vec3::new(0.0, 8.0, 0.0), + rho_kg_m3: RHO, + omega_rad_s: OMEGA, + }; + let (result, _) = rotor.compute_forces(&inputs, &rotor.initial_state()); + // Undo the rotation to read the in-plane hub-frame force back out. + let f_hub = r_hub.transpose() * result.F_world; + f_hub.0[0].hypot(f_hub.0[1]) + } + + let mag_level = hub_force_mag(&rotor, 0.0); + let mag_mid = hub_force_mag(&rotor, 45.0_f64.to_radians()); + let mag_aligned = hub_force_mag(&rotor, 89.0_f64.to_radians()); + + r.info("tilt_toward_wind", "h_force_level_N", mag_level, f64::NAN); + r.info("tilt_toward_wind", "h_force_45deg_N", mag_mid, f64::NAN); + r.info("tilt_toward_wind", "h_force_89deg_N", mag_aligned, f64::NAN); + r.assert_bool( + "tilt_toward_wind", + "shrinks_monotonically", + mag_mid, + mag_level, + mag_mid < mag_level && mag_aligned < mag_mid, + &format!( + "H-force should shrink as the disk tilts into the wind: level={:.4} N, 45deg={:.4} N, 89deg={:.4} N", + mag_level, mag_mid, mag_aligned + ), + ); + r.assert_bool( + "tilt_toward_wind", + "vanishes_when_aligned", + mag_aligned, + 0.0, + mag_aligned < 0.02, + &format!( + "H-force should be ~0 once the disk axis aligns with the wind, got {:.4} N", + mag_aligned + ), + ); + + // --- textbook H-force: drag opposing the direction of flight --- + // + // Classical rotor theory defines H-force as a drag-like force that + // opposes the aircraft's forward flight, not "whatever direction the + // wind happens to blow a stationary rotor". Model that directly: the + // hub itself moves at +X through still air (wind_world = 0), so the + // relative wind seen by the disk is v_rel = wind - v_hub = -X (apparent + // headwind from the nose). The in-plane force should point rearward + // (-X, opposing the direction of flight), matching the standard + // definition of rotor profile-drag H-force. + let inputs_forward_flight = RotorInputs { + collective_rad: 0.0, + tilt_lon: 0.0, + tilt_lat: 0.0, + R_hub: Mat3::eye(), + v_hub_world: Vec3::new(8.0, 0.0, 0.0), + wind_world: Vec3::zero(), + rho_kg_m3: RHO, + omega_rad_s: OMEGA, + }; + let (result_ff, _) = rotor.compute_forces(&inputs_forward_flight, &rotor.initial_state()); + let f_ff = result_ff.F_world.0; + r.info("forward_flight_still_air", "f_north_N", f_ff[0], f64::NAN); + r.assert_bool( + "forward_flight_still_air", + "drag_opposes_flight_direction", + f_ff[0], + 0.0, + f_ff[0] < -0.05, + &format!( + "flying forward (+X) through still air should give a rearward (-X) H-force, got F_north={:.3} N", + f_ff[0] + ), + ); + + // --- flapping-tilt H-force: flapback drag --- + // + // Classical rotor theory's second (usually larger) H-force term comes + // from blade flapping: in forward flight the advancing/retreating lift + // asymmetry drives a 1/rev flap response that, with the ~90 deg + // aerodynamic-damping phase lag, tilts the disk AFT (flapback). The + // thrust vector tilting rearward is an extra rearward in-plane force. + // + // Model it directly with a freely-hinged flap blade (omega_NR=0, so no + // hub moment is transmitted -- this isolates the flapping-tilt H-force + // from the transmitted-moment path) at a realistic Lock number, in + // forward flight with real thrust (nonzero collective). Compare against + // the identical rigid-blade rotor: enabling flap must make the H-force + // MORE rearward, and the flapping contribution should dominate the + // small profile-drag term. + let i_beta = 0.03; // -> Lock number ~8 (see info below) + let defn_flap = theory_rotor_flap(8, i_beta); + let rotor_flap = QuasiStaticBEM::build(defn_flap, 36, theory_polar()); + let rotor_rigid = QuasiStaticBEM::build(theory_rotor(8, 0.0), 36, theory_polar()); + + let collective = 8.0_f64.to_radians(); + let ff = |c: f64| RotorInputs { + collective_rad: c, + tilt_lon: 0.0, + tilt_lat: 0.0, + R_hub: Mat3::eye(), + v_hub_world: Vec3::new(8.0, 0.0, 0.0), // forward flight +X, still air + wind_world: Vec3::zero(), + rho_kg_m3: RHO, + omega_rad_s: OMEGA, + }; + + let (res_flap, _) = rotor_flap.compute_forces(&ff(collective), &rotor_flap.initial_state()); + let (res_rigid, _) = rotor_rigid.compute_forces(&ff(collective), &rotor_rigid.initial_state()); + let f_flap = res_flap.F_world.0; + let f_rigid = res_rigid.F_world.0; + + r.info("flapback", "lock_number", lock_number(i_beta), f64::NAN); + r.info("flapback", "thrust_z_N", -f_flap[2], f64::NAN); + r.info("flapback", "f_north_flap_N", f_flap[0], f64::NAN); + r.info("flapback", "f_north_rigid_N", f_rigid[0], f64::NAN); + + // Flapback must be rearward (-X). + r.assert_bool( + "flapback", + "flapping_h_force_is_rearward", + f_flap[0], + 0.0, + f_flap[0] < -0.05, + &format!( + "flapping in forward flight (+X) should tilt the disk aft -> rearward (-X) H-force, got F_north={:.3} N", + f_flap[0] + ), + ); + // Flapping adds to (dominates) the rigid-blade profile-drag H-force. + r.assert_bool( + "flapback", + "flapping_adds_rearward_drag", + f_flap[0], + f_rigid[0], + f_flap[0] < f_rigid[0] - 0.05, + &format!( + "enabling flap should make the H-force more rearward: flap F_north={:.3} N vs rigid {:.3} N", + f_flap[0], f_rigid[0] + ), + ); + // A freely-hinged blade transmits no hub moment: the flapping shows up + // in the H-force, not as an airframe roll/pitch moment. + let m_flap = res_flap.M_hub_world.0; + r.info("flapback", "hub_moment_x_N_m", m_flap[0], f64::NAN); + r.info("flapback", "hub_moment_y_N_m", m_flap[1], f64::NAN); + r.assert_bool( + "flapback", + "free_hinge_transmits_no_moment", + m_flap[0].hypot(m_flap[1]), + 0.0, + m_flap[0].hypot(m_flap[1]) < 1e-6, + &format!( + "freely-hinged flap should transmit ~zero hub moment, got |M|={:.3e} N*m", + m_flap[0].hypot(m_flap[1]) + ), + ); +} diff --git a/validation_rs/src/checks/harris_hforce_empirical.csv b/validation_rs/src/checks/harris_hforce_empirical.csv new file mode 100644 index 0000000..b1518ce --- /dev/null +++ b/validation_rs/src/checks/harris_hforce_empirical.csv @@ -0,0 +1,13 @@ +label,mu,alpha_deg,n_rpm,ct_meas,ch_meas,ch_max_err +mu014,0.1438,16.54,137.6,0.006054,0.000175,0.838207100363872 +mu020,0.1967,10.44,137.6,0.005892,0.000271,0.4418713393979635 +mu025,0.2479,7.36,137.6,0.005763,0.000334,0.3134231423489616 +mu030,0.2987,5.35,137.6,0.005461,0.000428,0.05114374548820241 +mu035,0.349,4.26,137.6,0.005365,0.000466,0.034536489325208586 +mu040,0.3993,3.39,137.6,0.00515,0.000528,0.15475382764989576 +mu045,0.4494,2.97,137.6,0.005228,0.000594,0.20672673195595606 +mu050,0.4995,2.57,137.6,0.005194,0.000683,0.2920380088660227 +mu055,0.5496,2.1,137.6,0.005046,0.000749,0.3631780001194528 +mu060,0.5996,1.97,137.6,0.005287,0.000875,0.42144772248736306 +mu065,0.6497,1.81,137.6,0.005467,0.000956,0.44362762441920345 +mu070,0.6998,1.27,137.6,0.00517,0.001033,0.5190008105470585 diff --git a/validation_rs/src/checks/harris_hforce_empirical.rs b/validation_rs/src/checks/harris_hforce_empirical.rs new file mode 100644 index 0000000..8563a33 --- /dev/null +++ b/validation_rs/src/checks/harris_hforce_empirical.rs @@ -0,0 +1,136 @@ +// Harris NASA CR-2008-215370 PCA-2 rotor in-plane hub force (H-force) vs +// advance ratio. Primary tabulated data (report pages 504-505, extracted from +// the text layer) lives in +// Research/Harris_CR-2008-215370/pages_504_505_appendix_pca2_coeffs.md. +// +// Method: the PCA-2 rotor autorotates, so blade pitch is not tabulated. For +// each operating point (mu, shaft alpha, RPM) the collective is trimmed so the +// quasi-static BEM matches the tabulated Rotor CT, then the resulting Rotor CH +// is compared against the tabulated value. Trimming to CT isolates the H-force +// model from the unknown pitch (and, since flapping does not change thrust in +// our model, the trim is independent of the flap DOF). +// +// The blade is modelled freely hinged (omega_NR = 0) so the flapping-tilt +// H-force is captured; for a free hinge that force is independent of blade +// inertia (see pca2_rotor_flap docs). Coefficients use the Harris/American +// normalization C = force / (rho * pi * R^2 * (Omega*R)^2). +// +// Per-case error ceilings live in harris_hforce_empirical.csv. To re-baseline +// after a model change: +// REWRITE_EMPIRICAL_CSV=1 cargo run --release -p validation_rs -- harris_hforce + +use crate::helpers::*; +use crate::report::Report; +use dynbem_rs::aero_model::AeroModel; +use dynbem_rs::polar::LinearPolar; +use dynbem_rs::quasi_static_bem::QuasiStaticBEM; +use serde::{Deserialize, Serialize}; + +#[derive(Deserialize, Serialize)] +struct Row { + label: String, + mu: f64, + alpha_deg: f64, + n_rpm: f64, + ct_meas: f64, + ch_meas: f64, + ch_max_err: f64, +} + +/// Compute (CT, CH) for the PCA-2 rotor at a given collective pitch and +/// operating point. CH is the in-plane hub-force coefficient (magnitude of the +/// disk-plane force), CT the thrust coefficient, both Harris-normalized. +fn ct_ch_at(rotor: &QuasiStaticBEM, pitch_deg: f64, row: &Row) -> (f64, f64) { + let inp = harris_pca2_inputs(pitch_deg, row.mu, row.alpha_deg, row.n_rpm); + let (res, _) = rotor.compute_forces(&inp, &rotor.initial_state()); + let f = res.F_world.0; + let norm = pca2_coeff_norm(row.n_rpm); + let ct = -f[2] / norm; + let ch = f[0].hypot(f[1]) / norm; + (ct, ch) +} + +/// Bisection trim: find the collective (deg) whose BEM CT matches ct_target. +/// CT is monotone-increasing in pitch over the bracket. Returns the trimmed +/// pitch and the (CT, CH) achieved there. +fn trim_to_ct(rotor: &QuasiStaticBEM, row: &Row) -> (f64, f64, f64) { + let mut lo = -10.0_f64; + let mut hi = 20.0_f64; + let (ct_lo, ch_lo) = ct_ch_at(rotor, lo, row); + let (ct_hi, ch_hi) = ct_ch_at(rotor, hi, row); + // Guard: if the target is outside the bracket, clamp to the nearer edge. + if row.ct_meas <= ct_lo { + return (lo, ct_lo, ch_lo); + } + if row.ct_meas >= ct_hi { + return (hi, ct_hi, ch_hi); + } + let mut mid = 0.5 * (lo + hi); + let (mut ct_mid, mut ch_mid) = ct_ch_at(rotor, mid, row); + for _ in 0..50 { + if (ct_mid - row.ct_meas).abs() < 1e-6 { + break; + } + if ct_mid < row.ct_meas { + lo = mid; + } else { + hi = mid; + } + mid = 0.5 * (lo + hi); + let pair = ct_ch_at(rotor, mid, row); + ct_mid = pair.0; + ch_mid = pair.1; + } + (mid, ct_mid, ch_mid) +} + +pub fn check_harris_hforce(r: &mut Report) { + r.begin_module( + "harris_hforce_empirical", + "QS BEM Rotor CH vs Harris CR-2008-215370 PCA-2 (trim-to-CT, per-case sweep)", + ); + + let csv_data = include_str!("harris_hforce_empirical.csv"); + let rewrite = std::env::var("REWRITE_EMPIRICAL_CSV").is_ok(); + let mut new_rows: Vec = Vec::new(); + + // Freely-hinged PCA-2 blade; i_beta value is irrelevant for a free hinge + // (see pca2_rotor_flap docs), Lock ~5 chosen only as a plausible number. + let i_beta = 1700.0; + let defn = pca2_rotor_flap(16, i_beta); + let polar = polar_for(&defn.airfoil); + let rotor = QuasiStaticBEM::build(defn, 48, polar); + + for rec in csv_rows(csv_data) { + let row: Row = rec + .deserialize(None) + .expect("harris_hforce_empirical.csv parse"); + let (pitch, ct, ch) = trim_to_ct(&rotor, &row); + let err_pct = (ch - row.ch_meas).abs() / row.ch_meas * 100.0; + + let case = format!("{} (mu={:.3})", row.label, row.mu); + r.info(&case, "trim_pitch_deg", pitch, f64::NAN); + r.info(&case, "CT_trimmed", ct, row.ct_meas); + r.info(&case, "CH_bem", ch, row.ch_meas); + r.check(&case, "CH_err_pct", err_pct, 0.0, row.ch_max_err * 100.0); + + if rewrite { + let new_tol = (err_pct / 100.0 + 0.03).max(0.03); + new_rows.push(Row { + ch_max_err: new_tol, + ..row + }); + } + } + + if rewrite { + let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR")) + .join("src/checks/harris_hforce_empirical.csv"); + let mut wtr = csv::Writer::from_path(&path).expect("open CSV for rewrite"); + for row in &new_rows { + wtr.serialize(row).expect("write row"); + } + wtr.flush().expect("flush"); + eprintln!("Rewrote {}", path.display()); + } +} diff --git a/validation_rs/src/checks/mod.rs b/validation_rs/src/checks/mod.rs index e34d92d..51fc31a 100644 --- a/validation_rs/src/checks/mod.rs +++ b/validation_rs/src/checks/mod.rs @@ -9,6 +9,8 @@ mod descent_empirical; mod flap_directional; mod flapping_harmonics; mod glauert_forward_inflow; +mod h_force; +mod harris_hforce_empirical; mod hover_castles_gray; mod hover_cq_empirical; mod hover_empirical; @@ -29,6 +31,8 @@ pub use descent_empirical::check_descent_empirical; pub use flap_directional::check_flap_directional; pub use flapping_harmonics::check_flapping_harmonics; pub use glauert_forward_inflow::check_glauert_forward_inflow; +pub use h_force::check_h_force; +pub use harris_hforce_empirical::check_harris_hforce; pub use hover_castles_gray::check_hover_castles_gray; pub use hover_cq_empirical::check_hover_cq_empirical; pub use hover_empirical::check_hover_empirical; @@ -42,7 +46,7 @@ pub use wake_skew::check_wake_skew; use crate::report::Report; -/// Run all 17 validation checks in order. +/// Run all 18 validation checks in order. pub fn run_all_checks(report: &mut Report) { run_filtered_checks(report, None); } @@ -67,6 +71,8 @@ pub fn run_filtered_checks(report: &mut Report, filter: Option<&str>) { maybe!("autorotation", check_autorotation); maybe!("measured_companions", check_measured_companions); maybe!("cyclic_sign", check_cyclic_sign); + maybe!("h_force", check_h_force); + maybe!("harris_hforce", check_harris_hforce); maybe!("servo_flap", check_servo_flap); maybe!("cyclic_phase_servo", check_cyclic_phase_servo); maybe!("flap_directional", check_flap_directional); diff --git a/validation_rs/src/helpers.rs b/validation_rs/src/helpers.rs index 36a752d..6354e8a 100644 --- a/validation_rs/src/helpers.rs +++ b/validation_rs/src/helpers.rs @@ -153,6 +153,26 @@ pub fn pca2_rotor(n_elements: usize) -> RotorDefinition { } } +/// PCA-2 rotor with a freely-hinged flap DOF (omega_NR = 0), used by the +/// Harris H-force check. `i_beta` is the blade flap inertia [kg*m^2]. +/// +/// NOTE: for a freely-hinged blade (omega_NR = 0, so nu_beta = 1 exactly), +/// the 1/rev flap response is a damping-limited resonance whose amplitude is +/// forcing/damping -- INDEPENDENT of the blade inertia (see +/// bem_common::apply_flap_dynamics: i_beta cancels out of b1c/b1s when +/// a_stiff = 0). So the flapping-tilt H-force this produces does not depend on +/// the exact `i_beta` value; it is set here to a physically plausible number +/// (Lock ~5 for the 6.85 m PCA-2 blade) only because the field is required. +pub fn pca2_rotor_flap(n_elements: usize, i_beta: f64) -> RotorDefinition { + let mut defn = pca2_rotor(n_elements); + defn.flap = Some(FlapProperties { + I_blade_flap_kgm2: i_beta, + omega_nr_rad_s: 0.0, + }); + defn.name = "pca2_harris_flap".to_string(); + defn +} + // --------------------------------------------------------------------------- // VPM rotor creation // --------------------------------------------------------------------------- @@ -404,6 +424,34 @@ pub fn bem_cq(q_spin: f64, rpm: f64) -> f64 { q_spin / (RHO * PI * R_TIP * R_TIP * (omega * R_TIP).powi(2) * R_TIP) } +/// PCA-2 forward-flight RotorInputs at a given shaft angle of attack, mirroring +/// `wheatley_fc` but for the quasi-static BEM (RotorInputs) path. Disk held +/// level (R_hub = identity); the freestream is applied as the hub velocity so +/// the relative wind carries both the edgewise (mu) and axial (shaft-tilt) +/// components. Used by the Harris CH check. +pub fn harris_pca2_inputs(pitch_deg: f64, mu: f64, alpha_deg: f64, n_rpm: f64) -> RotorInputs { + let omega = omega_from_rpm(n_rpm); + let a = alpha_deg.to_radians(); + let v = omega * PCA2_R * mu / a.cos(); + RotorInputs { + collective_rad: pitch_deg.to_radians(), + tilt_lon: 0.0, + tilt_lat: 0.0, + R_hub: Mat3::eye(), + v_hub_world: Vec3::new(v * a.cos(), 0.0, -v * a.sin()), + wind_world: Vec3::zero(), + rho_kg_m3: RHO, + omega_rad_s: omega, + } +} + +/// Harris/American normalization divisor for PCA-2 coefficients: +/// `C = force / (rho * pi * R^2 * (Omega*R)^2)`. +pub fn pca2_coeff_norm(n_rpm: f64) -> f64 { + let omega = omega_from_rpm(n_rpm); + RHO * PI * PCA2_R * PCA2_R * (omega * PCA2_R).powi(2) +} + /// Run a BEM model for n_steps with ExplicitEuler and return the final CT. pub fn run_ct(model: &M, inp: &RotorInputs, rpm: f64, n_steps: usize) -> f64 { let mut state = model.initial_state();