src/cpp/wheelrailcontact.cpp - master - Loylik/fastsim

src/cpp/wheelrailcontact.cpp
122 строки4 KB
GAlekseyV
unify API between C++ and Python versions
28 янв 2026, 13:30
28 янв 2026, 13:30075d3ca
// Реализации (header-only, простые оценки)
#include "wheelrailcontact.hpp"

#include <iostream>
#include <numbers>

namespace fastsim {
ContactPatch calculateContactPatch(const ContactGeometry& geom) {
  double K1 = (1 - geom.poissonRatio * geom.poissonRatio) /
              (std::numbers::pi * geom.youngModulus);
  double K2 = (1 - geom.poissonRatio * geom.poissonRatio) /
              (2 * std::numbers::pi * geom.youngModulus);
  double K3 = (1.0 / 2.0) * (1 / geom.wheelRadius + 1 / geom.railRadius);
  double K4 =
      (1.0 / 2.0) * ((1 / geom.wheelRadius) - (1 / geom.railRadius));

  double theta = std::acos(K4 / K3);
  double m = 1.018;
  double n = 0.9835;

  ContactPatch result;
  result.a = m * std::cbrt(3 * geom.normalForce * std::numbers::pi *
                           (K1 + K2) / (4 * K3));
  result.b = n * std::cbrt(3 * geom.normalForce * std::numbers::pi *
                           (K1 + K2) / (4 * K3));

  result.contactArea = std::numbers::pi * result.a * result.b;
  result.maxPressure = 3 * geom.normalForce / (2 * result.contactArea);

  return result;
}

void calculateLCoefficients(const ContactKinematics& kin, const Material& mat,
                            double& L1, double& L2, double& L3) {
  if (kin.a > kin.b) {
    double beta = kin.b / kin.a;
    L1 = (-0.119 - 0.486 * beta + 1.25 * sqrt(beta)) * kin.a / mat.G;
    L2 = (-0.142 - 0.337 * beta + 1.207 * sqrt(beta)) * kin.a / mat.G;
    L3 = (-0.122 - 0.383 * beta + 1.039 * sqrt(beta)) * kin.a / mat.G;
  } else {
    double beta = kin.a / kin.b;
    L1 = (0.818 - 0.196 * beta + 0.024 * sqrt(beta)) * kin.a / mat.G;
    L2 = (1.112 - 0.519 * beta + 0.133 * beta * beta) * kin.a / mat.G;
    L3 = (0.527 + 0.428 * beta - 0.719 * beta * beta +
          0.298 * beta * beta * beta) *
         kin.a / mat.G;
  }
}

ContactOutput computeContact(const ContactKinematics& kin, const Material& mat) {
  double L1 = 0;
  double L2 = 0;
  double L3 = 0;

  int m = kin.m;
  int n = kin.n;

  fastsim::calculateLCoefficients(kin, mat, L1, L2, L3);

  double w_x = 0;
  double w_y = 0;
  double x_L = 0;
  double x = 0;
  double y = 0;

  double dx = 0;
  double p_z = 0;
  double p = 0;

  double c = sqrt(kin.a * kin.b);

  double L = (fabs(kin.xi_x) * L1 + fabs(kin.xi_y) * L2 +
              c * fabs(kin.fi) * L3) /
             sqrt(kin.xi_x * kin.xi_x + kin.xi_y * kin.xi_y +
                  c * c * kin.fi * kin.fi);

  double dy = 2.0 * kin.b / m;
  double Fx = 0.0;
  double Fy = 0.0;
  double maxPz = 0.0;

  for (int j = 1; j <= m; j++) {
    y = -kin.b + (j - 1.0 / 2.0) * dy;
    x_L = kin.a * sqrt(1.0 - y * y / (kin.b * kin.b));
    dx = 2.0 * x_L / n;
    double p_y = 0.0;
    double p_x = 0.0;
    bool Z1 = false;
    bool Z2 = true;
    bool Z3 = true;

    for (int k = 1; k <= n; k++) {
      x = x_L - (k - 1.0 / 2.0) * dx;
      double f = kin.f_kin;
      if (Z1) {
        f = kin.f_stat;
      }
      p_z = 2.0 * kin.Fz * (x_L * x_L - x * x) /
            (std::numbers::pi * kin.a * kin.a * kin.a * kin.b);
      maxPz = std::max(p_z, maxPz);
      w_x = kin.xi_x - kin.fi * y;
      w_y = kin.xi_y + kin.fi * (x + 1.0 / 2.0 * dx);
      p_x = p_x - w_x * dx / L;
      p_y = p_y - w_y * dx / L;
      p = sqrt(p_x * p_x + p_y * p_y);

      if (p > f * p_z) {
        p_x = p_x * kin.f_kin * p_z / p;
        p_y = p_y * kin.f_kin * p_z / p;
        Z2 = false;
        Z3 = false;
      } else {
        Z2 = true;
      }
      Fx = Fx + p_x * dx * dy;
      Fy = Fy + p_y * dx * dy;
    }
  }
  return ContactOutput{Fx, Fy, maxPz};
}

}  // namespace fastsim