using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks; using ToolPathParser; namespace ClassicalCuttingForce { internal class ConvexCurvesForce { private const double PI = Math.PI; private const double AngularIncrement = 0.006; // h1 private const double AxialIncrement = 0.006; // da public void CalculateCuttingForce(int direction, ToolPosition position, CuttingForceCoefficients coefficients, ToolParameters toolParams) { // 刀具参数 int N = toolParams.NumberOfTeeth; // 刀齿数 N double r = toolParams.ToolRadius; // 刀具半径 r double beta = toolParams.HelixAngle; // 刀具螺旋角 beta double R_bull = toolParams.ArcRadius; // 环形圆弧半径(mm) double a_p = position.AxialDepth; // 轴向切削深度 a_p double a_e = position.RadialDepth; // 径向切削深度 a_e double R = position.CurvatureRadius; // 目标曲率半径 R double nrot = position.SpindleSpeed; // 主轴转速 nrot double fd = position.FeedRate; // 进给速度 fd // 计算中间变量 double thetaPit = 2.0 * PI / N; // 齿间角 double beta_rad = beta * PI / 180.0; // 螺旋角转换为弧度 double Rr = r - R_bull; // 径向偏移量(mm) double kb = Math.Tan(beta * PI / 180.0) / r; // 刀具常数 double trot = 60.0 / nrot; // 主轴旋转周期 double c = fd / (nrot * N); // 每齿进给量 double e = R + r; // 注意：这里是 R+r，不是 R-r double H = 2.0 * PI * e; double nrev = fd / H; double trev = 60.0 / nrev; double krev = trev / trot * N; double dphi = 2.0 * PI / krev; // 计算切入角 double thetaSt = Math.Asin(-((Math.Pow(a_e, 2) + 2.0 * R * a_e - 2.0 * Math.Pow(r, 2) - 2.0 * R * r) / (2.0 * (R + r))) / r) + 0.5 * PI; double thetaEx = PI; // 切出角 // 计算切厚最大角 double x1 = (2.0 * Math.Pow(R, 2) * Math.Cos(2.0 * dphi) - 2.0 * R * r - 2.0 * R * a_e + Math.Pow(a_e, 2) * Math.Cos(2.0 * dphi) - 2.0 * Math.Pow(R, 2) - Math.Pow(a_e, 2) + Math.Sin(2.0 * dphi) * Math.Sqrt(-a_e * (2.0 * R + a_e) * (a_e - 2.0 * r) * (2.0 * R + a_e + 2.0 * r)) + 2.0 * R * a_e * Math.Cos(2.0 * dphi) + 2.0 * R * r * Math.Cos(2.0 * dphi)) / (4.0 * Math.Sin(dphi) * (R + r)); double thetaMax = Math.Acos(x1 / r) + 0.5 * PI; // 计算相对坐标和距离 double x_rel = position.X - position.Cx; // P_coords - O_coords double y_rel = position.Y - position.Cy; double d = Math.Sqrt(x_rel * x_rel + y_rel * y_rel); // 计算角向和轴向积分数目 int K = (int)Math.Floor(2.0 * PI / AngularIncrement); int L = (int)Math.Floor(a_p / AxialIncrement); // 初始化结果数组 double[] Fx = new double[K]; double[] Fy = new double[K]; double[] Fz = new double[K]; double[] Fx_global = new double[K]; double[] Fy_global = new double[K]; double[] Fz_global = new double[K]; // 环形刀轴向几何参数数组 double[] r_z = new double[L]; // z处的刀具半径 double[] kappa_z = new double[L]; // z处的轴向浸入角(rad) double[] psi_z = new double[L]; // z处的径向滞后角(rad) double[] db_z = new double[L]; // z处的微元切削宽度(mm) double[] ds_z = new double[L]; // z处的微元刃口长度(mm) double[] tool_helixangel_z = new double[L];// z处的螺旋角(rad) int m = 0; int n = 0; int M = L * N * K; double[] h = new double[M]; double[] R1 = new double[N * K]; // 主计算循环 for (int i = 0; i < K; i++) { double R2 = i * AngularIncrement; // 螺旋槽底部切削刃的接触角 for (int k = 0; k < N; k++) { R1[n] = R2 + k * thetaPit; // 第k个齿当前角度的位置 double phi_j = R1[n]; for (int j = 0; j < L; j++) { double alph = (j + 1) * AxialIncrement; double z = (j + 1) * AxialIncrement;// 当前轴向高度(mm) // 计算z处的刀具半径r(z)（分区计算） if (z >= 0 && z < R_bull) // 环形圆弧区(MN区) { r_z[j] = Rr + Math.Sqrt(R_bull * R_bull - Math.Pow(R_bull - z, 2)); } else // 圆柱区(NS区) { r_z[j] = r; } // 计算轴向浸入角、螺旋角、径向滞后角（分区） if (z >= 0 && z < R_bull) // 环形圆弧区 { kappa_z[j] = Math.Acos((R_bull - z) / R_bull); tool_helixangel_z[j] = Math.Atan((r_z[j] - Rr) * Math.Tan(beta_rad) / R_bull); psi_z[j] = (z / R_bull) * Math.Tan(tool_helixangel_z[j]); } else // 圆柱区 { kappa_z[j] = PI / 2.0; tool_helixangel_z[j] = beta_rad; psi_z[j] = (z * Math.Tan(beta_rad)) / r; } // 计算微元切削宽度db(z)，避免除零 double sin_kap = Math.Sin(kappa_z[j]); if (sin_kap < 1e-6) sin_kap = 1e-6; db_z[j] = AxialIncrement / sin_kap; // 计算微元刃口长度ds(z) ds_z[j] = AxialIncrement / (Math.Cos(tool_helixangel_z[j]) * sin_kap); // 修正径向浸入角，角度归一化 double phi_j_z = phi_j - psi_z[j]; while (phi_j_z > 2 * PI) phi_j_z -= 2 * PI; while (phi_j_z < 0) phi_j_z += 2 * PI; if (phi_j_z >= thetaSt && phi_j_z <= thetaEx) { double chipThickness; double phi = phi_j_z - 0.5 * PI; if (phi_j_z <= thetaMax) { // lp计算公式 double sqrt_term = 2.0 * R * a_e - 2.0 * R * r + Math.Pow(R, 2) * Math.Pow(Math.Sin(phi), 2) + Math.Pow(r, 2) * Math.Pow(Math.Sin(phi), 2) + Math.Pow(a_e, 2) - Math.Pow(r, 2) + 2.0 * R * r * Math.Pow(Math.Sin(phi), 2); double lp = R * Math.Sin(phi) - Math.Sqrt(sqrt_term) + r * Math.Sin(phi); chipThickness = r - lp; } else { // lQ计算公式 double lQ = CalculateLQ(R, r, a_e, dphi, phi); chipThickness = r - lQ; } h[m] = chipThickness; // 计算微分切削力 double dft = (coefficients.ktc * h[m]) * db_z[j] + coefficients.kte * ds_z[j]; double dfr = (coefficients.krc * h[m]) * db_z[j] + coefficients.kre * ds_z[j]; double dfu = (coefficients.kuc * h[m]) * db_z[j] + coefficients.kue * ds_z[j]; if (direction == 1) { double[,] T = { { -Math.Cos(phi_j_z), -Math.Sin(kappa_z[j]) * Math.Sin(phi_j_z), -Math.Cos(kappa_z[j]) * Math.Sin(phi_j_z) }, { Math.Sin(phi_j_z), -Math.Sin(kappa_z[j]) * Math.Cos(phi_j_z), -Math.Cos(kappa_z[j]) * Math.Cos(phi_j_z) }, { 0, Math.Cos(kappa_z[j]), -Math.Sin(kappa_z[j]) } }; double fxComponent = T[0, 0] * dft + T[0, 1] * dfr + T[0, 2] * dfu; double fyComponent = T[1, 0] * dft + T[1, 1] * dfr + T[1, 2] * dfu; double fzComponent = T[2, 0] * dft + T[2, 1] * dfr + T[2, 2] * dfu; Fx[i] += fxComponent; Fy[i] += fyComponent; Fz[i] += fzComponent; } else if (direction == 2) { // 坐标变换矩阵 double[,] T = { { -Math.Cos(phi_j_z), -Math.Sin(kappa_z[j]) * Math.Sin(phi_j_z), -Math.Cos(kappa_z[j]) * Math.Sin(phi_j_z) }, { Math.Sin(phi_j_z), -Math.Sin(kappa_z[j]) * Math.Cos(phi_j_z), -Math.Cos(kappa_z[j]) * Math.Cos(phi_j_z) }, { 0, -Math.Cos(kappa_z[j]), Math.Sin(kappa_z[j]) } }; double fxComponent = T[0, 0] * dft + T[0, 1] * dfr + T[0, 2] * dfu; double fyComponent = T[1, 0] * dft + T[1, 1] * dfr + T[1, 2] * dfu; double fzComponent = T[2, 0] * dft + T[2, 1] * dfr + T[2, 2] * dfu; Fx[i] += fxComponent; Fy[i] += fyComponent; Fz[i] += fzComponent; } m++; } } n++; } if (direction == 1) { // 坐标转换到工件坐标系 double Fx_from_Y1 = Fy[i] * (x_rel / d); double Fy_from_Y1 = Fy[i] * (y_rel / d); double Fx_from_X1 = Fx[i] * (y_rel / d); double Fy_from_X1 = Fx[i] * (-x_rel / d); Fx_global[i] = Fx_from_X1 + Fx_from_Y1; Fy_global[i] = Fy_from_X1 + Fy_from_Y1; Fz_global[i] = Fz[i]; } else if (direction == 2) { // 坐标转换到工件坐标系 double Fx_from_Y1 = Fy[i] * (x_rel / d); double Fy_from_Y1 = Fy[i] * (y_rel / d); double Fx_from_X1 = Fx[i] * (y_rel / d); double Fy_from_X1 = Fx[i] * (-x_rel / d); Fx_global[i] = Fx_from_X1 + Fx_from_Y1; Fy_global[i] = Fy_from_X1 + Fy_from_Y1; Fz_global[i] = -Fz[i]; } } // 计算并设置结果 SetResults(position, Fx_global, Fy_global, Fz_global, K); } // lQ计算公式 private double CalculateLQ(double R, double r, double a_e, double dphi, double phi) { double term1 = Math.Pow(R, 2) * Math.Cos(dphi + 2.0 * phi); double term2 = 2.0 * R * r; double term3 = Math.Pow(r, 2) * Math.Cos(dphi + 2.0 * phi); double term4 = (Math.Pow(R, 2) * Math.Cos(2.0 * phi)) / 2.0; double term5 = (Math.Pow(r, 2) * Math.Cos(2.0 * phi)) / 2.0; double term6 = (Math.Pow(R, 2) * Math.Cos(2.0 * dphi + 2.0 * phi)) / 2.0; double term7 = (Math.Pow(r, 2) * Math.Cos(2.0 * dphi + 2.0 * phi)) / 2.0; double term8 = Math.Pow(R, 2); double term9 = Math.Pow(R, 2) * Math.Cos(dphi); double term10 = Math.Pow(r, 2) * Math.Cos(dphi); double term11 = 2.0 * R * r * Math.Cos(dphi + 2.0 * phi); double term12 = R * r * Math.Cos(2.0 * phi); double term13 = R * r * Math.Cos(2.0 * dphi + 2.0 * phi); double term14 = 2.0 * R * r * Math.Cos(dphi); double sqrtContent = term1 - term2 + term3 - term4 - term5 - term6 - term7 - term8 + term9 + term10 + term11 - term12 - term13 + term14; double lQ = Math.Sqrt(sqrtContent) - R * Math.Sin(dphi + phi) - r * Math.Sin(dphi + phi) + R * Math.Sin(phi) + r * Math.Sin(phi); return lQ; } private void SetResults(ToolPosition position, double[] Fx_global, double[] Fy_global, double[] Fz_global, int K) { // 计算平均值 double sumFx = 0, sumFy = 0, sumFz = 0; for (int i = 0; i < K; i++) { sumFx += Fx_global[i]; sumFy += Fy_global[i]; sumFz += Fz_global[i]; } position.AverageFx = sumFx / K; position.AverageFy = sumFy / K; position.AverageFz = sumFz / K; // 计算最大值 double maxAbsFx = 0, maxAbsFy = 0, maxAbsFz = 0; int idxFx = 0, idxFy = 0, idxFz = 0; for (int i = 0; i < K; i++) { if (Math.Abs(Fx_global[i]) > maxAbsFx) { maxAbsFx = Math.Abs(Fx_global[i]); idxFx = i; } if (Math.Abs(Fy_global[i]) > maxAbsFy) { maxAbsFy = Math.Abs(Fy_global[i]); idxFy = i; } if (Math.Abs(Fz_global[i]) > maxAbsFz) { maxAbsFz = Math.Abs(Fz_global[i]); idxFz = i; } } position.MaxFx = Fx_global[idxFx]; position.MaxFy = Fy_global[idxFy]; position.MaxFz = Fz_global[idxFz]; } } }