rods/Agents/MATH_SPEC.md

Math specification — literature backbone

Status: Living document. Every implemented term should trace to a citation below or to an explicit heuristic calibration note in code + docs/engineering/validation.md.

Primary references: see references/papers/README.md for citation details and source links.

• Romero, A., and Almeida, A. R. (2014). A Numerical Sucker-Rod Pumping Analysis Tool. SPE Artificial Lift Conference - Latin America and the Caribbean. SPE-169395-MS
• Everitt, T. A., and Jennings, J. W. (1992). An Improved Finite-Difference Calculation of Downhole Dynamometer Cards for Sucker-Rod Pumps. SPE Production Engineering. SPE-18189-PA
• Araujo, O., et al. (SPE-173970). 3D Rod String Dynamics in Deviated Wells (minimum-curvature and coupled dynamics reference used for trajectory coupling).
• Eisner, B., Langbauer, C., and Fruhwirth, R. (2022). A finite element approach for dynamic sucker-rod diagnostics (Newmark + Rayleigh + diagnostic iteration basis).
• Lukasiewicz, H. (as summarized in Eisner et al.) coupled axial/lateral force balance for deviated rod strings.

---

1. Gibbs one-dimensional damped wave (vertical / uniform rod)

Romero & Almeida — Eq. (2) (viscous damped wave, constant (A,E,\rho)):

[ \frac{\partial^2 u}{\partial t^2} = a^2 \frac{\partial^2 u}{\partial x^2} - \varsigma \frac{\partial u}{\partial t} ]

with (a^2 = E/\rho) (wave speed) and (\varsigma) viscous damping coefficient per unit mass (paper uses lumped fluid damping narrative).

Variable cross-section (Everitt & Jennings — Eq. 2 form) after multiplying through by (\rho A); axial stiffness gradient:

[ \frac{\partial}{\partial x}\left( EA \frac{\partial u}{\partial x} \right) = \rho A \frac{\partial^2 u}{\partial t^2} + c \rho A \frac{\partial u}{\partial t} - f_{\text{body}} ]

where (f_{\text{body}}) collects distributed weight and buoyancy along the rod (gravity along tangent; buoyancy from fluid — Lukasiewicz axial force balance in Eisner Eq. (2) discussion / Lukasiewicz coupled model referenced in Eisner).

Implementation note: Code uses discrete (E_i, A_i, \rho_i) per node or segment, harmonic or measured surface (u(0,t)), and bottom boundary coupling (pump / valve / spring-damper per roadmap).

---

2. Damping conventions

• Gibbs dimensionless damping (\nu): related to decay rate; Romero cites (\varsigma) proportional to velocity; Eisner Eq. (1) uses (\frac{\pi a \nu}{2L}\frac{\partial u}{\partial t}) form for vertical damped wave.
• Rayleigh damping (FEA, Eisner): (\mathbf{C} = \alpha \mathbf{M} + \beta \mathbf{K}) for bar elements.
• XML factors: UpStrokeDampingFactor, DownStrokeDampingFactor, NonDimensionalFluidDamping modulate effective (\gamma) or (\alpha,\beta) per phase (mapped in solver; see code comments).

---

3. Explicit finite-difference stencil (diagnostic / deviated extension)

Everitt & Jennings — Eq. (3) (conceptual explicit FD recursion transferring displacements downhole from surface card). The paper gives a five-point relation in ((i,j)) space (space index (i), time (j)) with coefficients involving (a), (c), (\Delta t), (\Delta x), and (EA).

Variable (EA): harmonic mean or segment fluxes:

[ \text{flux}{i+\frac12} = E{i+\frac12} A_{i+\frac12} \left( u_{i+1} - u_i \right) ]

Laplacian-like term formed from (\partial/\partial x (EA \partial u/\partial x)) discretization (matches current diagnostic JS prior to C port).

CFL: explicit wave requires (\Delta t \le \text{CFL} \cdot \Delta x / a_{\max}) with (a_{\max} = \max \sqrt{E/\rho}). Code clamps effective CFL (see verbose.numerics.cflEffective).

---

4. Deviated wells — coupled axial / lateral (reference)

Lukasiewicz (via Eisner Eqs. (2)–(3)) — force balance along rod tangent and normal; includes (\rho g A \cos\phi), viscosity term (\nu), Coulomb (\mu N), curvature (R), and lateral equilibrium.

Current implementation: distributed side load (N_i) and inclination-aware Coulomb friction (F_{f,i} = \mu_{\text{eff}} N_i \operatorname{sgn}(v_i)) are computed per node and injected into FDM/diagnostic/FEA updates. API exposes profiles.sideLoadProfile and profiles.frictionProfile when options.enableProfiles=true.

---

5. Three-dimensional trajectory (SPE-173970)

Araujo et al. — Eqs. (5)–(24): minimum-curvature method between survey stations: unit tangent, curvature angle (\gamma_i), radius (r_{c,i}), binormal, center of curvature, position propagation (R(s_{i+1})), and interpolation of any MD (s) within a segment.

Implementation: solver-c/src/trajectory.c uses minimum-curvature style tangent interpolation and exports node-wise (\kappa(s)), inclination, and azimuth on the rod grid for side-load/friction coupling and trajectory3D output.

---

6. Dynamic 1D bar FEM (Eisner et al.)

• Bar stiffness / mass: consistent element (K_e), (M_e) for linear shape functions.
• Newmark-β ((\beta=\frac14), (\gamma=\frac12)) for transient axial dynamics.
• Bottom BC: spring-damper + friction; diagnostic mode adds iterative plunger load adjustment so that computed top reaction matches measured polished load (Eisner Fig. 4 principle — restart / bisection per time step).

---

7. Boundary conditions

Mode	Surface	Bottom
Predictive	Harmonic (u(0,t)) (Romero Eq. (4) style crank motion) or measured surrogate	Pump: spring-damper + friction + valve-state pressure balance; gas fraction inferred from chamber state
Diagnostic	Measured (u(0,t)) and (F(0,t)) from dynamometer card (Everitt)	Same pump BC; FEA adjusts unknown plunger load and reports valve/gas timeline

Surface card QA (Eisner narrative): ≥ 75 samples recommended for Fourier-style tools; we enforce minimum samples and cycle checks in API (POST /solve/validate-card).

---

8. Outputs (API / solver)

• Polished and downhole cards ((x, F)).
• Pump movement: plunger position series, velocity, stroke, period.
• Optional profiles: stressProfile, sideLoadProfile, trajectory3D (when schemaVersion >= 2).
• Comparison: FDM vs FEA peak deltas + point-wise residuals + optional Fourier analytical baseline (comparison.fourier when enabled).

---

9. Symbols (SI internal)

Symbol	Unit	Meaning
(u)	m	Axial displacement
(F)	N	Axial force (tension positive)
(E)	Pa	Young’s modulus
(A)	m²	Cross-sectional area
(\rho)	kg/m³	Rod density
(\phi)	rad	Inclination from vertical
(\kappa)	1/m	Path curvature

---

10. Implementation quality rules (for next execution)

• Every newly introduced solver term must be tagged as one of:
  • equation-backed (with paper/equation reference),
  • heuristic (with rationale + planned replacement).
• No silent heuristic drift: if a coefficient changes, update validation fixtures and notes.
• Any new field wired into equations must be reflected in docs/engineering/field-traceability.md.
• Any new coupled term must include at least:
  • one unit-scale numerical test,
  • one integration case assertion,
  • one regression guard.

---

End of math spec backbone. Extend with equation numbers from PDFs as features land.