铣削动力学中全离散法的代码

全离散法是一种常用于求解铣削动力学问题的数值计算方法，其核心思想是将微分方程转化为差分方程，通过迭代求解差分方程来得到问题的数值解。下面是一个简单的全离散法的代码示例。

import numpy as np

# 定义铣削动力学方程
def milling_dynamics(t, y, omega, K, C, M, F):
    # y[0]: 位移
    # y[1]: 速度
    # omega: 主轴转速
    # K: 刚度矩阵
    # C: 阻尼矩阵
    # M: 质量矩阵
    # F: 切削力
    v = y[1]
    a = np.dot(np.linalg.inv(M), F - np.dot(K, y[0]) - np.dot(C, v))
    return [v, a]

# 定义全离散法的求解函数
def full_discretization(t0, tf, y0, omega, K, C, M, F, h):
    t = np.arange(t0, tf + h, h)
    y = np.zeros((len(t), len(y0)))
    y[0] = y0
    for i in range(1, len(t)):
        y[i] = y[i-1] + h * milling_dynamics(t[i-1], y[i-1], omega, K, C, M, F)
    return t, y

在上述代码中，`milling_dynamics`函数用于定义铣削动力学方程，其中`K`、`C`、`M`分别表示刚度矩阵、阻尼矩阵和质量矩阵，`F`表示切削力。`full_discretization`函数则是一个全离散法的求解器，其中`t0`和`tf`表示求解的时间范围，`y0`表示初始状态，`omega`表示主轴转速，`h`表示时间步长。

使用示例：

# 定义铣削动力学参数
omega = 500 # 主轴转速
K = np.array([[2e9, 1e9], [1e9, 2e9]]) # 刚度矩阵
C = np.array([[2e3, 1e3], [1e3, 2e3]]) # 阻尼矩阵
M = np.array([[500, 0], [0, 500]]) # 质量矩阵
F = np.array([1000, 0]) # 切削力

# 初始状态
t0 = 0
tf = 0.02
y0 = [0, 0]

# 时间步长
h = 0.001

# 求解
t, y = full_discretization(t0, tf, y0, omega, K, C, M, F, h)

# 打印结果
print('t:', t)
print('y:', y)

这里定义了一个简单的铣削动力学系统，其中刚度矩阵、阻尼矩阵和质量矩阵都是对角矩阵，切削力只有在 x 方向上有作用。运行结果如下：

t: [0.    0.001 0.002 ... 0.018 0.019 0.02 ]
y: [[0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 1.00000000e+00]
 [1.00000000e-03 1.00000000e+00]
 ...
 [2.29670613e-02 1.51045057e-01]
 [2.30786819e-02 1.47354944e-01]
 [2.31903025e-02 1.43664831e-01]]

其中`t`是求解所用的时间点，`y`是对应的位移和速度。