CN112364452A - 一种气体静压主轴转子运动轨迹多场耦合数值分析方法 - Google Patents

一种气体静压主轴转子运动轨迹多场耦合数值分析方法 Download PDF

Info

Publication number
CN112364452A
CN112364452A CN202011236278.0A CN202011236278A CN112364452A CN 112364452 A CN112364452 A CN 112364452A CN 202011236278 A CN202011236278 A CN 202011236278A CN 112364452 A CN112364452 A CN 112364452A
Authority
CN
China
Prior art keywords
main shaft
dimensionless
time step
axis
theta
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011236278.0A
Other languages
English (en)
Other versions
CN112364452B (zh
Inventor
陈国达
陈燚杰
卢奇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University of Technology ZJUT
Original Assignee
Zhejiang University of Technology ZJUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University of Technology ZJUT filed Critical Zhejiang University of Technology ZJUT
Priority to CN202011236278.0A priority Critical patent/CN112364452B/zh
Publication of CN112364452A publication Critical patent/CN112364452A/zh
Application granted granted Critical
Publication of CN112364452B publication Critical patent/CN112364452B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/12—Timing analysis or timing optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Magnetic Bearings And Hydrostatic Bearings (AREA)
  • Connection Of Motors, Electrical Generators, Mechanical Devices, And The Like (AREA)

Abstract

一种气体静压主轴转子运动轨迹多场耦合数值分析方法，包括：(1)：首先建立空气轴承的数学模型，采用交替隐式差分格式(ADI)对无量纲瞬态雷诺方程进行离散，采用Thomas追赶法对离散方程进行求解，得到无量纲气膜压力分布，通过气膜压力的积分得到轴承气膜恢复力和气膜力矩；(2)建立气体静压主轴中的永磁同步电机(PMSM)二维模型，利用电机二维切片近似等效三维电机，从而降低电机数值分析难度，计算得到三维永磁同步电机的径向电磁力及力矩；(3)建立气体静压主轴转子的5自由度动力学方程，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程，进行多场耦合下的转子动力学数值分析，求解可得转子运动轨迹。

Description

一种气体静压主轴转子运动轨迹多场耦合数值分析方法
技术领域
本发明涉及气体静压主轴转子运动仿真分析领域，尤其涉及一种气体静压主轴转子运动轨迹多场耦合数值分析方法。
背景技术
气体静压主轴通常由电机直接驱动，由气体轴承支承。由于其结构紧凑、效率高、无污染、精度高等优点，因此，它在超精密加工机床、精密仪器、光电装备等领域有着广泛的应用。
动力系数法和恢复力法是研究气体轴承-转子系统非线性动力学行为的两种主要的分析方法。最常用的是摄动法，该方法适应参数的线性化形式，通过求解摄动雷诺方程得到动力系数(即刚度和阻尼)。此外，还采用了相应阻抗的有理函数逼近和递推最小二乘方法来求得空气轴承的动力系数。采用动态系数法研究气体静压主轴的非线性动态行为时，通常将动态系数视为常数，这确实提高了效率，但忽略了动态系数变化所带来的影响。动力系数法是基于动力系数的获取，在考虑时变动态系数时，需要做大量的额外工作。而恢复力法可以通过求解瞬态雷诺方程直接得到，相比之下，恢复力法似乎更有效。瞬态雷诺方程可以采用有限元法(FEM)或有限差分法(FDM)求解，由于网格划分和迭代较少，有限差分法更容易获得瞬态恢复力，特别是在考虑多个物理场耦合作用的情况。目前，已知电磁因素对系统稳定性的影响不容忽视，而电磁因素对气体静压主轴动态特性的影响还有待进一步研究与探索，但已有的转子运动轨迹数值分析方法尚未实现对电磁因素的耦合。实现气体静压主轴转子运动轨迹的多场耦合数值分析是研究主轴运动误差形成机理的重要手段，对于进一步提高气体静压主轴运动精度有重要意义。
发明内容
本发明要克服现有技术的上述不足，提出一种气体静压主轴转子运动轨迹多场耦合数值分析方法，能够有效地将电磁因素耦合至转子动力学方程，实现多物理场作用下转子运动轨迹的数值分析。
本发明的一种气体静压主轴转子运动轨迹多场耦合数值分析方法，包括如下步骤：
步骤一：首先建立空气轴承的数学模型，包括径向轴承和止推轴承，其主导方程为雷诺方程，对瞬态雷诺方程进行无量纲处理。采用交替隐式差分格式(ADI)对无量纲瞬态雷诺方程进行离散，采用Thomas追赶法对离散方程进行求解，得到无量纲气膜压力分布，通过气膜压力的积分得到轴承气膜恢复力和气膜力矩。
步骤二：建立气体静压主轴中的永磁同步电机(PMSM)二维模型，采用有限元方法求解，在得到径向电磁力后，通过曲线拟合的方法对其进行简化，得到简化后的数学模型，便于后续计算。利用电机二维切片近似等效三维电机，从而降低电机数值分析难度，计算得到三维永磁同步电机的径向电磁力及力矩。
步骤三：建立气体静压主轴转子的5自由度动力学方程，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程，进行多场耦合下的转子动力学数值分析，求解可得转子运动轨迹。
与现有方法相比，本发明的有益效果是：
1)考虑了气体静压主轴的多场耦合因素，特别是考虑了电磁力因素，能够更为精确地进行转子运动轨迹的数值分析。
2)采用有限差分法结合恢复力法能够大大提高计算效率。
3)可以为气体静压主轴的多目标优化设计提供理论基础。
附图说明
图1是本发明的径向气膜网格划分及边界条件示意图；
图2是本发明的止推气膜网格划分及边界条件示意图；
图3是本发明的无量纲气膜压力求解实例图；
图4是本发明中永磁同步电机实例示意图；
图5是本发明中永磁同步电机电磁力实例仿真图；
图6是本发明中永磁同步电机电磁力实例曲线拟合图；
图7是本发明的求解流程图；
图8是本发明中案例求解结果主轴轴端Z向位移；
图9是本发明中案例求解结果主轴轴端在X-Y面内轨迹；
图10是本发明中案例求解结果主轴轴端空间轨迹；
具体实施方式
下面结合附图与具体实施案例对本发明作进一步详细说明。
本发明所述方法可由安装在计算机上的软件程序实现。
本发明的一种气体静压主轴转子运动轨迹多场耦合数值分析方法，包括以下步骤：
步骤一：首先建立空气轴承的数学模型，包括径向轴承和止推轴承，其主导方程为雷诺方程，对瞬态雷诺方程进行无量纲处理。采用交替隐式差分格式(ADI)对无量纲瞬态雷诺方程进行离散，采用Thomas追赶法对离散方程进行求解，得到无量纲气膜压力分布，通过气膜压力的积分得到轴承气膜恢复力和气膜力矩。
步骤一的具体实施内容如下：
建立径向轴承和止推轴承的瞬态雷诺方程如，并对其无量纲化，其结果分别为：
径向轴承：
Figure BDA0002766806480000041
其中参数为：
Figure BDA0002766806480000042
τ＝ωt，
Figure BDA0002766806480000043
上述式(1)及其中间变量中，P为无量纲气膜压力；p为气膜压力；Pr为参考压力；H为无量纲气膜厚度；h为气膜厚度；Cr为初始气膜厚度；X为无量纲周向长度；x为气膜周向坐标；L为径向轴承特征长度；Z为轴向无量纲长度；z为气膜轴向坐标；τ为无量纲时间；ω为主轴转速；τ为时间；t为时间；Λ为轴承数；
Figure BDA0002766806480000044
为法向气体流速；
Figure BDA0002766806480000045
为径向轴承气体流量；δk为克罗内科符号；μ为气体动力粘度；pa为环境压力；ρa为环境空气密度；ρ为气体密度；
止推轴承：
Figure BDA0002766806480000046
其中参数为：
Figure BDA0002766806480000047
τ＝ωt，
Figure BDA0002766806480000048
上述式(2)及其中间变量中，θ为无量纲周向长度；xr为轴承周向坐标；r为轴承径向坐标；R为无量纲径向长度；
Figure BDA0002766806480000051
为止推轴承气体流量；Rm为止推轴承特征长度；Λt为轴承数；
用S＝P2代入式(1)，采用交替隐式差分格式(ADI)对式(1)进行离散，分别对Z方向与θ方向进行差分，第n时间步在Z方向进行差分得到：
Figure BDA0002766806480000052
Figure BDA0002766806480000053
Figure BDA0002766806480000054
Figure BDA0002766806480000055
Figure BDA0002766806480000056
第n+1时间步在X方向进行差分得到：
Ei，jSn+1 i，j-1+Fi，jSn+1 i，j+Gi，jSn+1 i，j+1＝Ii，j (8)
Figure BDA0002766806480000057
Figure BDA0002766806480000058
Figure BDA0002766806480000059
Figure BDA0002766806480000061
式(3)-(11)中的变量Ai，j Bi，j Ci，j Di，j Ei，j Fi，j Gi，j Ii，j仅是为了方便书写，不具有特殊物理含义；ΔZ为Z方向网格步长；Δθ为θ方向网格步长；Δτ为无量纲时间步长；在第n时间步中，上标为
Figure BDA0002766806480000062
表示该物理量为当前时间步的值，上标为n表示该物理量为上一时间步的值；在第n+1时间步中，上标为n+1表示该物理量为当前时间步的值，上标为
Figure BDA0002766806480000063
表示该物理量为上一时间步的值；
采用Thomas追赶法对差分后的方程进行求解。其边界条件附图1所示，Z＝0为大气边界条件，Z＝M为对称边界条件，θ＝0和θ＝2π为周期性边界条件。
用S＝P2代入式(2)，采用交替隐式差分格式(ADI)对式(2)进行离散，分别对R方向与θ方向进行差分，第n时间步在R方向进行差分得到：
Figure BDA0002766806480000064
Figure BDA0002766806480000065
Figure BDA0002766806480000066
Figure BDA0002766806480000067
Figure BDA0002766806480000071
第n+1时间步在θ方向进行差分得到：
Eti，jSn+1 i，j-1+Fti，jSn+1 i，j+Gti，jSn+1 i，j+1＝Iti，j (18)
Figure BDA0002766806480000072
Figure BDA0002766806480000073
Figure BDA0002766806480000074
Figure BDA0002766806480000075
式(13)-(22)中的变量Ati，j Bti，j Cti，j Dti，j Eti，j Fti，j Gti，j Iti，j仅是为了方便书写，不具有特殊物理含义；Δθ为θ方向网格步长；ΔR为R方向网格步长；Δτ为无量纲时间步长；在第n时间步中，上标为
Figure BDA0002766806480000076
表示该物理量为当前时间步的值，上标为n表示该物理量为上一时间步的值；在第n+1时间步中，上标为n+1表示该物理量为当前时间步的值，上标为
Figure BDA0002766806480000081
表示该物理量为上一时间步的值；
采用Thomas追赶法对差分后的方程进行求解。其边界条件附图2所示，R＝Rb为大气边界条件，R＝Ra为对称边界条件，θ＝0和θ＝2π为周期性边界条件。
求解得到无量纲气膜压力分布，如附图3所示为某一具体实例的气膜压力分布，对气膜压力进行积分可得到气膜恢复力以及气膜力矩：
Figure BDA0002766806480000082
Figure BDA0002766806480000083
Figure BDA0002766806480000084
Figure BDA0002766806480000085
Figure BDA0002766806480000086
Figure BDA0002766806480000087
Figure BDA0002766806480000088
上述式(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)模型，可采用有限元方法进行求解，永磁同步电机的电磁力可通过对磁通密度进行积分得到：
Figure BDA0002766806480000091
Figure BDA0002766806480000092
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为电机节点到主轴质心的距离；
附图4所示为永磁同步电机结构示例，对某一具体永磁同步电机实例进行仿真计算得到电磁力如附图5所示，附图6为该实例中对积分得到的电磁力进行曲线拟合，对电磁力进行简化，得到其简化的数学表达形式如式(36)(37)所示：
Figure BDA0002766806480000101
Figure BDA0002766806480000102
上述式(36)(37)中，ε为主轴偏心率；θr为主轴径向电磁力与x轴所成角度；
在主轴实际工作中倾斜角很小的条件下，忽略轴向电磁力，将二维电机切片在轴向阵列，则可近似模拟三维电机，降低电机计算难度。
步骤三：建立气体静压主轴转子的5自由度动力学方程，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程，进行多场耦合下的转子动力学数值分析，求解可得转子运动轨迹。
步骤三的具体实施内容如下：
将主轴转子视为刚体，考虑主轴的5自由度运动。对速度、加速度、角速度及角加速度作线性化处理，建立气体静压主轴转子的动力学方程，如式(38)-(52)所示，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程：
Figure BDA0002766806480000103
Figure BDA0002766806480000104
Figure BDA0002766806480000105
Figure BDA0002766806480000106
Figure BDA0002766806480000107
Vx＝Vx0+AxΔτ (43)
Vy＝Vy0+AyΔτ (44)
VZ＝VZ0+AzΔτ (45)
ΩX＝ΩX0+AXΔτ (46)
ΩY＝ΩY0+AYΔτ (47)
Figure BDA0002766806480000111
Figure BDA0002766806480000112
Figure BDA0002766806480000113
Figure BDA0002766806480000114
Figure BDA0002766806480000115
其中的各个参数可由以下方式得到：
Figure BDA0002766806480000116
Figure BDA0002766806480000117
Figure BDA0002766806480000118
Figure BDA0002766806480000119
上述式(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为主轴长度；参照附图7所示的求解流程，对某求解案例进行求解：图8为主轴轴端Z向位移，图9为主轴轴端在X-Y面内轨迹，图10为主轴轴端空间轨迹。
附图所示结果可有效证明该方法的正确性，该求解方法效率高，能够为气体静压主轴的多目标优化设计、提高气体静压主轴的运动精度提供重要的理论基础。
本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也包括本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (1)

1.一种气体静压主轴转子运动轨迹多场耦合数值分析方法，包括以下步骤：
步骤一：首先建立空气轴承的数学模型，包括径向轴承和止推轴承，其主导方程为雷诺方程，对瞬态雷诺方程进行无量纲处理；采用交替隐式差分格式(ADI)对无量纲瞬态雷诺方程进行离散，采用Thomas追赶法对离散方程进行求解，得到无量纲气膜压力分布，通过气膜压力的积分得到轴承气膜恢复力和气膜力矩；具体包括：
建立径向轴承和止推轴承的瞬态雷诺方程如，并对其无量纲化，其结果分别为：
径向轴承：
Figure FDA0002766806470000011
其中参数为：
Figure FDA0002766806470000012
τ＝ωt，
Figure FDA0002766806470000013
上述式(1)及其中间变量中，P为无量纲气膜压力；p为气膜压力；Pr为参考压力；H为无量纲气膜厚度；h为气膜厚度；Cr为初始气膜厚度；X为无量纲周向长度；x为气膜周向坐标；L为径向轴承特征长度；Z为轴向无量纲长度；z为气膜轴向坐标；τ为无量纲时间；ω为主轴转速；τ为时间；t为时间；Λ为轴承数；
Figure FDA0002766806470000014
为法向气体流速；
Figure FDA0002766806470000015
为径向轴承气体流量；δk为克罗内科符号；μ为气体动力粘度；pa为环境压力；ρa为环境空气密度；ρ为气体密度；
止推轴承：
Figure FDA0002766806470000016
其中参数为：
Figure FDA0002766806470000021
τ＝ωt，
Figure FDA0002766806470000022
上述式(2)及其中间变量中，θ为无量纲周向长度；xr为轴承周向坐标；r为轴承径向坐标；R为无量纲径向长度；
Figure FDA0002766806470000023
为止推轴承气体流量；Rm为止推轴承特征长度；Λt为轴承数；
用S＝P2代入式(1)，采用交替隐式差分格式(ADI)对式(1)进行离散，分别对Z方向与θ方向进行差分，第n时间步在Z方向进行差分得到：
Figure FDA0002766806470000024
Figure FDA0002766806470000025
Figure FDA0002766806470000026
Figure FDA0002766806470000027
Figure FDA0002766806470000028
第n+1时间步在X方向进行差分得到：
Ei，jSn+1 i，j-1+Fi，jSn+1 i，j+Gi，jSn+1 i，j+1＝Ii，j (8)
Figure FDA0002766806470000029
Figure FDA00027668064700000210
Figure FDA00027668064700000211
Figure FDA0002766806470000031
式(3)-(11)中的变量Ai，j Bi，j Ci，j Di，j Ei，j Fi，j Gi，j Ii，j仅是为了方便书写，不具有特殊物理含义；ΔZ为Z方向网格步长；Δθ为θ方向网格步长；Δτ为无量纲时间步长；在第n时间步中，上标为
Figure FDA0002766806470000032
表示该物理量为当前时间步的值，上标为n表示该物理量为上一时间步的值；在第n+1时间步中，上标为n+1表示该物理量为当前时间步的值，上标为
Figure FDA0002766806470000033
表示该物理量为上一时间步的值；
采用Thomas追赶法对差分后的方程进行求解；其边界条件附图1所示，Z＝0为大气边界条件，Z＝M为对称边界条件，θ＝0和θ＝2π为周期性边界条件；
用S＝P2代入式(2)，采用交替隐式差分格式(ADI)对式(2)进行离散，分别对R方向与θ方向进行差分，第n时间步在R方向进行差分得到：
Figure FDA0002766806470000034
Figure FDA0002766806470000035
Figure FDA0002766806470000036
Figure FDA0002766806470000037
Figure FDA0002766806470000041
第n+1时间步在θ方向进行差分得到：
Eti，jSn+1 i，j-1+Fti，jSn+1 i，j+Gti，jSn+1 i，j+1＝Iti，j (18)
Figure FDA0002766806470000042
Figure FDA0002766806470000043
Figure FDA0002766806470000044
Figure FDA0002766806470000045
式(13)-(22)中的变量Ati，j Bti，j Cti，j Dti，j Eti，j Fti，j Gti，j Iti，j仅是为了方便书写，不具有特殊物理含义；Δθ为θ方向网格步长；ΔR为R方向网格步长；Δτ为无量纲时间步长；在第n时间步中，上标为
Figure FDA0002766806470000046
表示该物理量为当前时间步的值，上标为n表示该物理量为上一时间步的值；在第n+1时间步中，上标为n+1表示该物理量为当前时间步的值，上标为
Figure FDA0002766806470000051
表示该物理量为上一时间步的值；
采用Thomas追赶法对差分后的方程进行求解；其边界条件附图2所示，R＝Rb为大气边界条件，R＝Ra为对称边界条件，θ＝0和θ＝2π为周期性边界条件；
求解得到无量纲气膜压力分布，如附图3所示为某一具体实例的气膜压力分布，对气膜压力进行积分可得到气膜恢复力以及气膜力矩：
Figure FDA0002766806470000052
Figure FDA0002766806470000053
Figure FDA0002766806470000054
Figure FDA0002766806470000055
Figure FDA0002766806470000056
Figure FDA0002766806470000057
Figure FDA0002766806470000058
上述式(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)模型，可采用有限元方法进行求解，永磁同步电机的电磁力可通过对磁通密度进行积分得到：
Figure FDA0002766806470000061
Figure FDA0002766806470000062
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)所示：
Figure FDA0002766806470000063
Figure FDA0002766806470000064
上述式(36)(37)中，ε为主轴偏心率；θr为主轴径向电磁力与x轴所成角度；
在主轴实际工作中倾斜角很小的条件下，忽略轴向电磁力，将二维电机切片在轴向阵列，则可近似模拟三维电机，降低电机计算难度；
步骤三：建立气体静压主轴转子的5自由度动力学方程，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程，进行多场耦合下的转子动力学数值分析，求解可得转子运动轨迹；具体包括：
将主轴转子视为刚体，考虑主轴的5自由度运动；对速度、加速度、角速度及角加速度作线性化处理，建立气体静压主轴转子的动力学方程，如式(38)-(52)所示，将得到的气膜恢复力、永磁同步电机的径向电磁力耦合至动力学方程：
Figure FDA0002766806470000071
Figure FDA0002766806470000072
Figure FDA0002766806470000073
Figure FDA0002766806470000074
Figure FDA0002766806470000075
Vx＝Vx0+AxΔτ (43)
Vy＝Vy0+AyΔτ (44)
VZ＝VZ0+AZΔτ (45)
ΩX＝ΩX0+AXΔτ (46)
ΩY＝ΩY0+AYΔτ (47)
Figure FDA0002766806470000081
Figure FDA0002766806470000082
Figure FDA0002766806470000083
Figure FDA0002766806470000084
Figure FDA0002766806470000085
其中的各个参数可由以下方式得到：
Figure FDA0002766806470000086
Figure FDA0002766806470000087
Figure FDA0002766806470000088
Figure FDA0002766806470000089
上述式(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为主轴长度。
CN202011236278.0A 2020-11-09 2020-11-09 一种气体静压主轴转子运动轨迹多场耦合数值分析方法 Active CN112364452B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011236278.0A CN112364452B (zh) 2020-11-09 2020-11-09 一种气体静压主轴转子运动轨迹多场耦合数值分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011236278.0A CN112364452B (zh) 2020-11-09 2020-11-09 一种气体静压主轴转子运动轨迹多场耦合数值分析方法

Publications (2)

Publication Number Publication Date
CN112364452A true CN112364452A (zh) 2021-02-12
CN112364452B CN112364452B (zh) 2024-05-03

Family

ID=74510318

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011236278.0A Active CN112364452B (zh) 2020-11-09 2020-11-09 一种气体静压主轴转子运动轨迹多场耦合数值分析方法

Country Status (1)

Country Link
CN (1) CN112364452B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114676577A (zh) * 2022-03-29 2022-06-28 西安交通大学 基于多物理场数字孪生模型的空气静压主轴模拟方法及系统
CN115062874A (zh) * 2022-08-16 2022-09-16 成都禀证科技有限责任公司 一种水体污染物监测预测分析方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104776998A (zh) * 2015-03-26 2015-07-15 北京工业大学 一种基于动态刚度系数和阻尼系数的转子轴心轨迹求解方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104776998A (zh) * 2015-03-26 2015-07-15 北京工业大学 一种基于动态刚度系数和阻尼系数的转子轴心轨迹求解方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114676577A (zh) * 2022-03-29 2022-06-28 西安交通大学 基于多物理场数字孪生模型的空气静压主轴模拟方法及系统
CN114676577B (zh) * 2022-03-29 2024-04-02 西安交通大学 基于多物理场数字孪生模型的空气静压主轴模拟方法及系统
CN115062874A (zh) * 2022-08-16 2022-09-16 成都禀证科技有限责任公司 一种水体污染物监测预测分析方法及系统
CN115062874B (zh) * 2022-08-16 2022-11-18 成都禀证科技有限责任公司 一种水体污染物监测预测分析方法及系统

Also Published As

Publication number Publication date
CN112364452B (zh) 2024-05-03

Similar Documents

Publication Publication Date Title
Chen Vibration modelling and verifications for whole aero-engine
CN104331565B (zh) 轴类磁悬浮刚性转子系统的动力学建模方法及控制方法
Liu et al. Field dynamic balancing for rigid rotor-AMB system in a magnetically suspended flywheel
Sun et al. A novel 4-DOF hybrid magnetic bearing for DGMSCMG
CN112364452A (zh) 一种气体静压主轴转子运动轨迹多场耦合数值分析方法
Cao et al. A new dynamic model of ball-bearing rotor systems based on rigid body element
CN113268901B (zh) 基于格子Boltzmann动压气体轴承间隙微流动仿真方法
Kessentini et al. Modeling and dynamics of a horizontal axis wind turbine
CN110096784B (zh) 一种具有轴向压差的径向滑动轴承的快速计算与设计方法
CN109359318A (zh) 空气轴承电主轴5自由度耦合刚性转子系统动态设计方法
Zhou et al. An assumed mode method and finite element method investigation of the coupled vibration in a flexible-disk rotor system with lacing wires
Meng et al. Frequency and stability analysis method of asymmetric anisotropic rotor-bearing system based on three-dimensional solid finite element method
CN102830242A (zh) 一种基于磁悬浮惯性执行机构的姿态角速度测量方法
CN105095583B (zh) 一种微尺度下的静压主轴模态分析方法
CN104965991A (zh) 基于传递函数的机翼颤振速度确定方法
CN117521244A (zh) 机动飞行状态下弹性支承结构振动响应分析方法及系统
Wu et al. Modeling and experimental study on the micro-vibration transmission of a control moment gyro
CN108984936A (zh) 高速双联滚动轴承电主轴转子系统动态设计方法
Su et al. Analysis of dynamic characteristic for misalignment-spline gear shaft based on whole transfer matrix method
CN111209639A (zh) 一种叶轮-轴承-转子系统的高效定量建模方法
Han et al. Experimental and computational analysis of microscale shrouded coaxial rotor in hover
Yim et al. Effect of tangential torque on the dynamics of flexible rotors
Acree Jr et al. High-speed wind tunnel tests of a full-scale proprotor on the tiltrotor test rig
Rao et al. A methodology for dynamic coefficients and nonlinear response of multi-lobe journal bearings
CN114676577B (zh) 基于多物理场数字孪生模型的空气静压主轴模拟方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant