1.一种气体静压主轴转子运动轨迹多场耦合数值分析方法，包括以下步骤：步骤一：首先建立空气轴承的数学模型，包括径向轴承和止推轴承，其主导方程为雷诺方程，对瞬态雷诺方程进行无量纲处理；采用交替隐式差分格式(ADI)对无量纲瞬态雷诺方程进行离散，采用Thomas追赶法对离散方程进行求解，得到无量纲气膜压力分布，通过气膜压力的积分得到轴承气膜恢复力和气膜力矩；具体包括：建立径向轴承和止推轴承的瞬态雷诺方程，并对其无量纲化，其结果分别为:径向轴承：

其中参数为：

上述式(1)及其中间变量中，P为无量纲气膜压力；p为气膜压力；Pr为参考压力；H为无量纲气膜厚度；h为气膜厚度；Cr为初始气膜厚度；X为无量纲周向长度；x为气膜周向坐标；L为径向轴承特征长度；Z为轴向无量纲长度；z为气膜轴向坐标；τ为无量纲时间；ω为主轴转速；τ为时间；t为时间；Λ为轴承数； 为法向气体流速； 为径向轴承气体流量；δk为克罗内科符号；μ为气体动力粘度；pa为环境压力；ρa为环境空气密度；ρ为气体密度；

止推轴承：

其中参数为：

上述式(2)及其中间变量中，θ为无量纲周向长度；xr为轴承周向坐标；r为轴承径向坐标；R为无量纲径向长度； 为止推轴承气体流量；Rm为止推轴承特征长度；Λt为轴承数；

2

用S＝P代入式(1)，采用交替隐式差分格式(ADI)对式(1)进行离散，分别对Z方向与θ方向进行差分，第n时间步在Z方向进行差分得到：第n+1时间步在X方向进行差分得到：

n+1 n+1 n+1

Ei,jS i,j‑1+Fi,jS i,j+Gi,jS i,j+1＝Ii,j     (8)ΔZ为Z方向网格步长；Δθ为θ方向网格步长；Δτ为无量纲时间步长；在第n时间步中，上标为 表示该物理量为当前时间步的值，上标为n表示该物理量为上一时间步的值；

在第n+1时间步中，上标为n+1表示该物理量为当前时间步的值，上标为 表示该物理量为上一时间步的值；

采用Thomas追赶法对差分后的方程进行求解；Z＝0为大气边界条件，Z＝M为对称边界条件，θ＝0和θ＝2π为周期性边界条件；

2

用S＝P代入式(2)，采用交替隐式差分格式(ADI)对式(2)进行离散，分别对R方向与θ方向进行差分，第n时间步在R方向进行差分得到：第n+1时间步在θ方向进行差分得到：

n+1 n+1 n+1

Eti,jS i,j‑1+Fti,jS i,j+Gti,jS i,j+1＝Iti,j (18)Δθ为θ方向网格步长；ΔR为R方向网格步长；Δτ为无量纲时间步长；在第n时间步中，上标为 表示该物理量为当前时间步的值，上标为n表示该物理量为上一时间步的值；

在第n+1时间步中，上标为n+1表示该物理量为当前时间步的值，上标为 表示该物理量为上一时间步的值；

采用Thomas追赶法对差分后的方程进行求解；R＝Rb为大气边界条件，R＝Ra为对称边界条件，θ＝0和θ＝2π为周期性边界条件；

求解得到无量纲气膜压力分布，对气膜压力进行积分可得到气膜恢复力以及气膜力矩：上述式(24)‑(29)中fbx,fby分别为作用于主轴x,y方向的气膜力；tbjx,tbjy分别为径向轴承作用于主轴x,y轴主轴的气膜力矩；tbtx,tbty分别为止推轴承作用于主轴x,y轴主轴的气膜力矩；θα为主轴偏心方向与x轴所成夹角；θβ为主轴绕z轴偏转角度；Ra为止推轴承外径；Rb为止推轴承内径；sinθt(i)为止推轴承节点沿轴心法向与x轴所成的夹角；Lj为气膜节点到主轴质心的距离；Rt为止推轴承节点到主轴质心的距离；

步骤二：建立气体静压主轴中的永磁同步电机(PMSM)二维模型，采用有限元方法求解，在得到径向电磁力后，通过曲线拟合的方法对其进行简化，得到简化后的数学模型，便于后续计算；利用电机二维切片近似等效三维电机，从而降低电机数值分析难度，计算得到三维永磁同步电机的径向电磁力及力矩；具体包括：建立气体静压主轴中的永磁同步电机(PMSM)模型，可采用有限元方法进行求解，永磁同步电机的电磁力可通过对磁通密度进行积分得到：fmx＝frcosθα+fθsinθα     (32)

fmy＝frsinθα+fθcosθα      (33)

tmx＝fmyLm      (34)

tmy＝fmxLm       (35)

上述式(30)‑(35)中，fr,fθ分别为径向、切向电磁力；Br,Bθ分别为径向、切向磁通密度；

μ0为空气磁导率；fmx,fmy分别为作用于主轴x,y方向的电磁力；tmx,tmy分别为作用于主轴x,y轴的电磁力矩；Lm为电机节点到主轴质心的距离；

对积分得到的电磁力进行曲线拟合，得到其简化的数学表达形式如式(36)(37)所示：上述式(36)(37)中，ε为主轴偏心率；θr为主轴径向电磁力与x轴所成角度；

在主轴实际工作中倾斜角很小的条件下，忽略轴向电磁力，将二维电机切片在轴向阵列，则可近似模拟三维电机，降低电机计算难度；

步骤三：建立气体静压主轴转子的5自由度动力学方程，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程，进行多场耦合下的转子动力学数值分析，求解可得转子运动轨迹；具体包括：将主轴转子视为刚体，考虑主轴的5自由度运动；对速度、加速度、角速度及角加速度作线性化处理，建立气体静压主轴转子的动力学方程，如式(38)‑(52)所示，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程：Vx＝Vx0+AxΔτ (43)

Vy＝Vy0+AyΔτ (44)

VZ＝VZ0+AZΔτ (45)

ΩX＝ΩX0+AXΔτ (46)ΩY＝ΩY0+AYΔτ (47)其中的各个参数可由以下方式得到：

上述式(38)‑(52)中，Ax,Ay,Az分别为当前时间步主轴在x,y,z方向的无量纲加速度；X,Y,Z分别为当前时间步主轴在x,y,z方向的无量纲位移；Fex,Fey,Fez分别为用于主轴x,y,z方向的无量纲外载荷；Fmx,Fmy,Fmz分别为作用于主轴x,y,z方向的无量纲电磁力；Fbx,Fby,Fbz分别为作用于主轴x,y,z方向的无量纲气膜力；AX,AY分别为当前时间步主轴绕x,y方向的角加速度；ΘX,ΘY分别为当前时间步主轴绕x,y方向的角度；Tex,Tey分别为作用于主轴x,y轴的无量外力矩；Tmx,Tmy为作用于主轴x,y轴的无量纲电磁力矩；Tbx,Tby分别为作用于主轴x,y轴的无量纲气膜力矩；M为主轴无量纲质量；IX,IY,IZ分别为主轴对x,y,z轴无量纲转动惯量；

Vx,Vy,Vz分别为当前时间步主轴在x,y,z方向无量纲速度；ΩX,ΩY分别为当前时间步主轴绕x,y轴角速度；X0,Y0,Z0分别为上一时间步主轴在x,y,z方向的无量纲位移；ΘX0,ΘY0分别为上一时间步主轴绕x,y方向的角度；Vx0,Vy0,Vz0分别为上一时间步主轴在x,y,z方向无量纲速度；ΩX0,ΩY0分别为上一时间步主轴绕x,y轴角速度；Δτ为无量纲时间间隔；fex,fey,fez分别为用于主轴x,y,z方向的外载荷；tex,tey分别为作用于主轴x,y轴的外力矩；ix,iy,iz分别为主轴对x,y,z轴转动惯量；

Ls为主轴长度。