网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

基于向量式有限元法的磁浮列车磁力耦合系统建模与数值分析  PDF

  • 孙友刚 1,2
  • 徐俊起 2
  • 王素梅 3
  • 袁野 4
  • 倪一清 3
1. 同济大学 铁道与城市轨道交通研究院,上海 201804; 2. 同济大学 国家磁浮交通工程技术研究中心,上海 201804; 3. 香港理工大学 国家轨道交通电气化与自动化工程技术研究中心香港分中心,香港 999077; 4. 江苏大学 电气信息工程学院,江苏 镇江 212036

中图分类号: U237

最近更新:2021-12-14

DOI:10.11908/j.issn.0253-374x.21210

  • 全文
  • 图表
  • 参考文献
  • 作者
  • 出版信息
EN
目录contents

摘要

针对中低速磁浮列车悬浮系统,基于向量式有限元法建立可变刚度的高架轨道梁模型,同时基于牛顿力学方程建立车辆系统模型,并通过可控悬浮电磁力将2个模型耦合。以轨道梁的跨中位移、梁端转角、振动加速度以及悬浮间隙偏差值为重要指标,从所提出的车‒桥磁力耦合模型出发,通过数值仿真得到磁浮列车及轨道线路相应结构构件的振动响应及位移变形响应规律。最后,通过全尺寸磁浮列车现场试验初步验证所提出的磁力耦合模型的有效性。

电磁悬浮(EMS)型磁浮列车是一种新型轨道交通工具,具有无接触、无摩擦、噪声低、易维护等优

1-3。近年来,我国磁浮交通尤其是低速磁浮交通得到迅猛发展,但在大规模商业化前,磁浮交通技术仍存在诸多亟待解决的问题,如磁浮控制系统可靠性、磁浮车辆的承载力、车‒轨耦合作用机理4-5。其中,磁浮列车的车‒轨耦合模型如何更加准确地反映磁浮交通的特征和机理,一直困扰着磁浮技术专家和工程师。磁浮交通的轨道梁多为高架结构,目前大部分磁浮车‒轨耦合建模方法要么忽略或简化非线性项,电磁力线性化后类似弹簧阻尼器(失去本质特征),要么只建立在轨道梁某个位置的静态悬浮(不能反映列车速度影响)。因此,亟需研究由可控电磁力无接触耦合的车‒轨动力学机理,建立反映不同速度下车辆和轨道垂向运行动态特性的计算模型。

在早期磁浮列车动态特性研究中,将磁浮列车简化为轨道上的一个移动

6或者移动质7,但是这样的模型无法反映车辆的动力学行为。翟婉明8和赵春发9同时考虑车辆系统和轨道系统,提出了用弹簧阻尼器替代多电磁铁磁浮的车‒轨耦合模型。然而,将电磁力简化为弹簧阻尼器虽然提高了计算效率,但是不能完全反映电磁力的非线性特征。Yau10考虑比例‒积分(PI)控制器下的非线性电磁力,将磁浮列车简化为一系列考虑二系悬挂系统的二自由振子模型,其中轨道被简化为单跨的欧拉梁。每个振子间没有任何连接,并且不能研究车辆的点头等运动状态。Han11提出了一种基于多体动力学的磁浮车辆动力学模型,该模型考虑了分段开关轨道的灵活性。利用所建立的分析模型,进行了中低速下悬浮稳定性和安全性的仿真研究。Han12基于多体动力学开发了EMS型磁浮列车系统集成模型,分析了列车沿曲线运行时侧风的影响和横向阻尼器的影响。Zhang13设计了低动力作用下的新型中低速磁浮列车悬浮架,并进行了动力学数值仿真验证。对于电磁力控制环节仿真,采用了简化的线性磁力控制模型。Sun14建立了EMS型磁浮列车的车‒轨耦合动力学模型,研究了车‒轨耦合共振的发生机理,但是未考虑动态跑车过程。近些年,学者15-16利用成熟的商业有限元软件进行轨道建模,再联合其他软件模块组成联合仿真模型,大大提高了精确度,但是大量的网格和多个软件成本增大了计算负担和花费,并且无法提供数学模型给控制器设计使用。

综上,EMS型磁浮列车系统的车‒轨动力学耦合机理和模型研究已取得较大进展,但是对不同速度和电磁力非线性共同作用下的动态跑车过程仍然研究不足。就磁浮车辆行驶在桥梁上的运动而言,涉及到2个子系统的相互耦合问题,而且这种耦合行为会随着时间的变化而改变。简化的模型不能完全反映耦合特性,而传统有限元在分析磁浮列车车‒桥的动力耦合行为时计算又非常耗时。向量式有限元法是以向量式力学为基础的新型有限元计算方法。该方法通过数值计算的方式来准确预测和分析结构的真实运动行为,还可以用来分析大变形、大变位、断裂、多体运动等复杂行为。采用点值描述的方式,将连续结构体以一组有限的质点来描述,每个质点可以独立且并行计算,故该方法可以用来分析多个运动变形体的相互作用。基于此,将向量式有限元法作为一种更加简单高效的计算方法进行磁浮列车车‒桥耦合作用分析。

针对上述问题,首先推导出磁浮列车的磁力耦合模型,其中高架轨道模型采用向量式有限元法。然后,以梁的跨中位移、梁端转角、振动加速度以及悬浮间隙偏差为重要指标,对不同车速下的磁浮列车动态通过轨道梁的过程进行仿真分析。最后,利用上海临港中低速磁浮基地的全尺寸试验样车,对该建模方法的准确性和有效性进行初步验证。

1 磁浮列车磁力耦合模型

以向量式有限元法进行车‒桥动力学分析时,可采用点值描述法将车辆和桥梁用一组质点表示,每个质点均满足牛顿第二定律。作用在车辆上的外力以及作用在桥梁上的外力为悬浮力的作用力和反作用力,施加在桥梁上每个质点的内力可以由梁单元的逆向运动得到。

图1所示,一节磁浮列车以匀速v行驶在简支梁上。车体由1个车厢、5个悬浮架组成。每个悬浮架有2个电磁铁模块(左右各1个),每个模块由4个电磁铁线圈构成。左右模块相同位置的每2个电磁铁线圈的电磁力简化为一个集中悬浮力。为了用向量式有限元法模拟车‒桥系统,对于车辆模型采用多体车辆有限元模型,对于桥梁采用向量式有限元梁单元模拟。在建模前,作出如下合理假设:①只考虑二维的车‒桥模型,忽略磁浮车辆的纵向运动(纵向为垂直于重力和轨道中心线的竖直平面方向);②桥梁采用欧拉‒伯努利梁来表示,桥梁的截面均匀。将桥梁离散成一组有质量的点,点与点之间用无质量的单元连接。

图1 磁浮列车通过多跨简支梁

Fig.1 Maglev train passing through multi-span simply supported beam

1.1 车体运动方程

在如图1所示的车‒桥模型中,Fiji=1,,5,j=1,,4)为电磁铁产生的悬浮力,mc为车厢质量,mbi(i=1,,5)为悬浮架质量,zbi为每个悬浮架质心垂向位移,θbi为每个悬浮架转动角度,zc为车厢质心垂向位移,θc为车厢转动角度,lm为2个悬浮架中心距离,lb为悬浮架质心与空气弹簧的距离,lnn=1,,10)为车厢质心到每个空气弹簧的距离,cc为空气弹簧阻尼,kc为空气弹簧刚度,Jc为车厢转动惯量,Jbi为悬浮架转动惯量。由上述模型可以得到车体的运动方程,如下所示:

磁浮车厢垂向沉浮运动时

mcz¨c+ccZ˙+kcZ=0 (1)
Z=10zc-2zb1-2zb2-2zb3-2zb4-2zb5

磁浮车体点头运动时

Jcθ¨c+ccθ¯˙L¯2+ccZ¯˙θL¯c+kcθ¯L¯2+kcZ¯θL¯c=0 (2)
θ¯=(θc,θb1,θb2,θb3,θb4,θb5)
Z¯θ=(zb1,zb2,zb3,zb4,zb5)
L¯2=(2(l12+l22+l32+l42+l52),-lbl1+lbl2,
-lbl3+lbl4,-lbl5-lbl6,-lbl7-lbl8,
-lbl9-lbl10)T
L¯c=(-l1-l2,-l3-l4,-l5+l6,l7+l8,l9+l10)T

悬浮架垂向沉浮运动时

M¯buZ¯¨bu+ccZ¯˙m+ccL¯buθ¯˙bu+kcL¯buθ¯bu=F¯bu (3)
M¯bu=mb100000mb200000mb300000mb400000mb5
Z¯¨bu=z¨b100000z¨b200000z¨b300000z¨b400000z¨b5
L¯bu=lb-l100000lb-l300000lb-l500000lb-l700000lb-l9
Z¯˙m=2(z˙b1-z˙c)000002(z˙b2-z˙c)000002(z˙b3-z˙c)000002(z˙b4-z˙c)000002(z˙b5-z˙c)
θ¯bu=θc+θb100000θc+θb200000θc+θb300000θc+θb400000θc+θb5
F¯bu=(-F11-F12-F13-F14,-F21-F22-F23-F24,-F31-F32-F33-F34,
-F41-F42-F43-F44,-F51-F52-F53-F54)T

悬浮架点头运动时

J¯bθ¯¨br+ccL¯brθ¯˙mix+kcL¯brθ¯mix=F¯br (4)
J¯b=Jb100000Jb200000Jb300000Jb400000Jb5
θ¯¨br=θ¨b100000θ¨b200000θ¨b300000θ¨b400000θ¨b5
L¯br=2lb-l1+l2000002lb-l3+l4000002lb-l5+l6000002lb-l7+l8000002lb-l9+l10
θ¯mix=θc+θb100000θc+θb200000θc+θb300000θc+θb400000θc+θb5
F¯br=-F11lb-F12lb2+F13lb2+F14lb,-F21lb-F22lb2+F23lb2+F24lb,-F31lb-F32lb2+F33lb2+F34lb,
-F41lb-F42lb2+F43lb2+F44lb,-F51lb-F52lb2+F53lb2+F54lbT

1.2 可控电磁力模型

对于磁浮列车的悬浮控制,采用双环控制,并取得了成功应

17。因为斩波器的成功应用,所以这里只设计电流环。采用经典的比例‒微分(PD)控制器,悬浮力F计算式如下所示:

F=KeKp(h-h0)+Kdh˙+I0h2 (5)

式中:Ke为电磁力常数;Kp为比例项系数;Kd为微分项系数;h为实时悬浮间隙;h0为理想悬浮间隙;I0为平衡电流。

1.3 轨道梁运动方程

对轨道梁作以下假设:

(1)轨道梁看作是无限长的欧拉‒伯努利梁。

(2)轨道扣件与桥梁之间采用刚性阻尼连接,如图2所示。

图2 桥计算简图

Fig.2 Schematic diagram of rail-bridge calculation

轨道和桥梁模型均采用向量式点值描述,每个点满足牛顿第二定律,而作用在轨道和桥梁上的每个质点的内力可以采用向量式梁单元求得。依据向量式理论,桥梁可采用一组质点进行模拟,相邻质点之间采用向量式梁单元连接,每个质点均满足牛顿第二定律,即:

MαX¨α=Fα,ext+Fα,int+Fα,d (6)

式中:MαX¨α分别为任一质点的质量矩阵和位移矩阵,MαX¨α为任一质点的惯性力;Fα,extFα,intFα,d分别为任一质点受到的内力、外力和阻尼力。

对于作用在桥梁质点上的内力,可采用向量式梁单元理论中的逆向运动进行求解。向量式有限元法的梁理论是采用点值描述的方法来描述梁受力后其位置的变化,以途径单元将梁的变形简化成每个途径单元内的小变形,这样梁的内力就可采用材料力学的公式来计算;为了得到节点位移中连接单元的纯变形,可将当前时刻梁的构型(1t‒2t)做一个虚拟的逆向运动(1d‒2d),如图3所示。此逆向运动包括逆向平动位移(-u1)和逆向转动位移-θ。通过梁的虚拟逆向运动,将刚体运动从全位移中分离,以便得到梁的纯变形,即

图3 梁单元的逆向运动

Fig.3 Reverse motion of beam element

Δe=lt-la,θ1=β1t-β1a-φ,θ2=β2t-β2a-θ (7)

式中:Δeθ1θ2分别为轴向纯变形、梁单元节点1和2的纯转动位移;ltla分别为梁单元的当前时刻

(1t‒2t)和初始时刻(1a‒2a)的单元长度;(β1tβ2t)和(β1aβ2a)分别为节点1和2在当前时刻和初始时刻的转动位移。

在得到梁单元的节点变形后,可采用一组内插函数来计算梁单元上任一点的变形,此插值函数即变形函数。在求得梁单元上每一点的变形后,根据节点内力平衡和节点变形所产生的虚功与梁单元的变形虚功相等的条件,即可求得梁单元的内力。

轨道梁上列车范围内的外力P仅包括轨道与悬浮架之间的悬浮力,即:

P=F (8)

式中:F为悬浮力。可以利用中心差分法求解车辆和桥梁的微分方

18-19,表达式如下所示:

Xα,-1=Xα,0-Δt1-a0Δt2MαX˙α,0+Fα,0,ext+Fα,0,int+a1Fα,0,intΔtΔt22Mα (9)
Xα,n+1=4C1Xα,n+C2C1Xα,n-1+2h2Fα,n,ext+Fα,n,int-a1Fα,n,int-Fα,n-1,intΔtC1Mα (10)

式中:a0a1为阻尼因子;C1=a0h+2C2=a0h-2分别为与阻尼因子a0有关的阻尼系数;Xα,-1Xα,0Xα,nXα,n-1Xα,n+1分别是-1时刻、初始时刻、n时刻、(n-1)时刻和(n+1)时刻的位移矩阵;Fα,0,extFα,n,ext分别是初始时刻和n时刻的外力矩阵;Fα,0,intFα,n-1,intFα,n,int分别是初始时刻、(n-1)时刻和n时刻的内力矩阵。

2 数值仿真与分析

根据所提出的磁浮列车磁力耦合模型,通过数值仿真来研究该方法的有效性和相关规律。工况为一节磁浮列车分别以60 km·h-1、80 km·h-1和100 km·h-1的速度通过五跨的简支轨道梁。每跨轨道梁之间有简支梁跨端结构缝。磁浮系统相关仿真参数的取值如表1表2所示。

表 1 欧拉伯努利梁参数
Tab.1 Parameters of Euler-Bernolli beam
跨长L/m

抗弯刚度EI/

(N·m2

每米质量m/(kg·m-1阻尼比ξ
24 3.15×1010 4.6×103 0.035
表 2 磁浮车辆参数
Tab.2 Parameters of maglev vehicle
物理量数值物理量数值
mc/kg 1.04×104 Jc/(kg·m4 3.85×105
mbi/kg 2.23×103 Jbi/(kg·m4 1.15×103
kc/(N·m-1 1.6×105 cc/(N·s·m-2 1.0×104
I0/A 25.74 lm/m 1.150
lb/m 0.575 l1/m 7.450
l2/m 5.150 l3/m 4.300
l4/m 2.000 l5/m 1.150
Kp 65 000 Kd 1 500

对60、80、100 km·h-1车速下的磁浮列车动态通过轨道梁的过程进行仿真,梁的跨中垂向位移、梁端转角、梁的中点垂向加速度、车厢垂向加速度、悬浮间隙偏差和电磁力的仿真结果如图4~9所示。

图4 第3根轨道梁的跨中垂向位移

Fig.4 Mid-span vertical displacement of the third track beam

图5 第3根轨道梁端转角

Fig.5 Angle at the end of the third track beam

图6 第3根轨道梁中点垂向加速度

Fig.6 Vertical acceleration at the midpoint of the third track beam

图7 车厢垂向加速度

Fig.7 Vertical acceleration of the carriage

图8 悬浮间隙变化

Fig.8 Change of airgap

图9 电磁力变化

Fig.9 Change of electromagnetic force

从图4~9可以看出,60、80、100 km·h-1车速下,梁的跨中垂向位移分别为2.012、2.014、2.015 mm,而梁端转角分别为0.268 4‰、0.269 4‰和 0.274 1‰。图89的周期波动是由车辆通过多跨简支梁时产生的,尖峰跳跃是由梁端边界条件产生的。由图4~9可以看出,在低速下轨道梁跨中垂向位移、梁端转角和车厢垂向加速度随速度变化的趋势并不明显,而轨道梁中点垂向加速度、悬浮间隙偏差以及电磁力的振动响应有随速度增大而增加的趋势。

3 试验验证

在同济大学国家磁浮交通工程技术研究中心的上海临港磁浮试验基地进行相关试验,对模型有效性进行初步验证。如图10a所示,上海临港中低速磁浮试验列车为五悬浮架结

20。每个悬浮架系统主要由悬浮架、空簧、悬浮电磁铁和悬浮控制器组成。控制算法在悬浮控制器内编程实现,可由上位机实时采集悬浮间隙和悬浮电流,在车厢里可以布置加速度传感器。如图10b所示,在磁浮轨道梁的中点位置布置传感器(压电式加速度传感器、激光位移传感器、拉线式位移传感器、光纤光栅应变传感器),采样频率为5 000 Hz。磁浮列车分别以20、40、60 km·h-1的速度通过高架轨道梁,并通过电压式数据采集仪以及光纤光栅数据解调仪采集相关试验数据。

图10 试验现场布置

Fig.10 Test on site

选取工况为:磁浮列车以60 km·h-1的速度通过轨道梁时的实测数据作对比分析。轨道梁跨中垂向位移实测值和仿真值的对比如图11所示。

图11 跨中垂向位移的实测值和仿真值对比

Fig.11 Comparison of mid-span displacement between measured data and simulation data

实测的最大跨中位移为2.013 mm,仿真的最大跨中位移为2.011 mm。实测数据和仿真数据的曲线重合度较高。平均电流、悬浮间隙和轨道梁振动一阶频率如表3所示。通过数据对比,初步验证了所提出的向量式有限元磁浮列车磁力耦合系统建模方法的有效性。

表3 各物理量实测数据和仿真数据对比(60 km·h-1
Tab.3 Comparison of various physical quantities between measured data and simulation data (60 km·h-1
物理量仿真值实测值
电流平均值/A 26.24 25.81
悬浮间隙平均值/mm 8.53 8.82

轨道梁中点垂向振动

一阶频率/Hz

6.49 6.85

4 结语

提出了一种基于向量式有限元法的磁力耦合计算模型。分别对60、80、100 km·h-1车速下的磁浮列车动态通过轨道梁的过程进行仿真,获得了梁的跨中垂向位移、梁端转角、轨道梁中点垂向加速度、车厢垂向加速度、悬浮间隙偏差和电磁力的仿真结果。可以看出,轨道梁中点垂向加速度、悬浮间隙偏差以及电磁力对速度的敏感性高于轨道梁跨中垂向位移、梁端转角和车厢垂向加速度。因此,在控制器设计时可考虑引入轨道垂向加速度作为状态反馈,为中低速磁浮的提速设计新的控制算法。最后,在同济大学国家磁浮交通工程技术研究中心的上海临港磁浮试验基地进行整车试验,验证了所提出的基于向量式有限元法的磁浮列车磁力耦合模型的有效性。

作者贡献声明

孙友刚:总体方案设计和仿真。

徐俊起:悬浮控制算法设计。

王素梅:向量式有限元建模。

袁 野:试验数据采集。

倪一清:试验方案设计。

参考文献

1

LEE H WKIM K CLEE J. Review of maglev train technologies [J]. IEEE Transactions on Magnetics2006427): 1917. [百度学术

2

KUSAGAWA SBABA JSHUTOH Ket al. Multipurpose design optimization of EMS-type magnetically levitated vehicle based on genetic algorithm [J]. IEEE Transactions on Applied Superconductivity2004142): 1922. [百度学术

3

ZHANG LHUANG JHUANG Let al. Stability and bifurcation analysis in a maglev system with multiple delays[J]. International Journal of Bifurcation & Chaos2015255):1. [百度学术

4

周丹峰李杰余佩倡. 磁浮交通轨排耦合自激振动分析及自适应控制方法[J].自动化学报20194512): 2328. [百度学术

ZHOU DanfengLI JieYU Peichanget al. Analysis and adaptive control of the track induced self-excited vibration for the maglev transport[J]. Acta Automatica Sinica20194512): 2328. [百度学术

5

WU HZENG XYU Y. Motion stability of high-speed maglev systems in consideration of aerodynamic effects: a study of a single magnetic suspension system [J]. Acta Mechanica Sinica2017336): 1084. [百度学术

6

CAI YCHEN S SROTE D Met al. Vehicle/guideway interaction for high speed vehicles on a flexible guideway [J]. Journal of Sound and Vibration19941755): 625. [百度学术

7

周又和武建军郑晓静.磁浮列车的动力稳定性分析与Liapunov指数[J].力学学报2000321):42. [百度学术

ZHOU YouheWU JianjunZHENG Xiaojinget al. Analysis of dynamic stability for magnetic levitation vehicles by Liapunov characteristic number[J]. Acta Mechanica Sinica2000321): 42. [百度学术

8

翟婉明赵春发.磁浮车辆/轨道系统动力学(Ⅰ): 磁/轨相互作用及稳定性[J].机械工程学报2005417):1. [百度学术

ZHAI WanmingZHAO Chunfa. Dynamics of maglev vehicle/guideway systems (Ⅰ):magnet/rail interaction and system stability[J]. Chinese Journal of Mechanical Engineering2005417): 1. [百度学术

9

赵春发翟婉明.磁浮车辆/轨道系统动力学(Ⅱ): 建模与仿真[J].机械工程学报2005418):163. [百度学术

ZHAO ChunfaZHAI Wanming. Dynamics of maglev vehicle/guideway system (Ⅱ):modeling and simulation[J]. Chinese Journal of Mechanical Engineering2005418): 163. [百度学术

10

YAU J D. Vibration control of maglev vehicles traveling over a flexible guideway [J]. Journal of Sound and Vibration20093211): 184. [百度学术

11

HAN J BHAN H SLEE J Met al. Dynamic modeling and simulation of EMS maglev vehicle to evaluate the levitation stability and operational safety over an elastic segmented switch track[J]. Journal of Mechanical Science & Technology2018327):2987. [百度学术

12

HAN J BKIM K J. Characteristics of vibration in magnetically levitated trains subjected to crosswind [J]. Proceedings of the Institution of Mechanical Engineers, Part F: Journal of Rail and Rapid Transit20182325): 1347. [百度学术

13

ZHANG MLUO SGAO Cet al. Research on the mechanism of a newly developed levitation frame with mid-set air spring[J]. Vehicle System Dynamics20185612): 1797. [百度学术

14

SUN YXU JQIANG Het al. Hopf bifurcation analysis ofmaglev vehicle-guideway interaction vibration system and stability control based on fuzzy adaptive theory[J]. Computers in Industry2019108197. [百度学术

15

刘德军李小珍洪沁烨.中低速磁浮列车‒大跨度连续梁耦合振动研究[J].铁道工程学报2017349):53. [百度学术

LIU DejunLI XiaozhenHONG Qinyeet al. Coupling vibration study of medium-low speed maglev train-long span continuous bridge system[J]. Journal of Railway Engineering Society2017349):53. [百度学术

16

向湘林龙志强.基于虚拟样机的中速磁浮列车单悬浮架耦合振动特性研究[J].黑龙江大学工程学报2020114):61. [百度学术

XIANG XianglinLONG Zhiqiang. Coupling vibration characteristics study of maglev single bogie based on virtual prototype [J]. Journal of Engineering of Heilongjiang University2020114):61. [百度学术

17

刘恒坤常文森.磁悬浮列车的双环控制[J].控制工程2007142):198. [百度学术

LIU HengkunCHANG Wensen. Double-loop control of maglev train [J]. Control Engineering of China2007142): 198. [百度学术

18

王素梅. 基于向量式有限元的风‒车‒轨‒桥耦合系统动力分析[D]. 杭州浙江大学2018. [百度学术

WANG Sumei. Vector form intrinsic finite element method-based dynamic analysis of wind-train-rail-bridge interaction system [D]. HangzhouZhejiang University2018. [百度学术

19

WANG S MYAU J DDUAN Y Fet al. Prediction of crosswind-induced derailment of train-rail-bridge system by vector mechanics [J]. Journal of Engineering Mechanics202014612): 04020132. [百度学术

20

孙友刚李万莉林国斌.低速磁浮列车悬浮系统动力学建模及非线性控制[J].同济大学学报(自然科学版)2017455):741. [百度学术

SUN YougangLI WanliLIN Guobinet al. Dynamic modeling and nonlinear control research on magnetic suspension systems of low-speed maglev train[J]. Journal of Tongji University(Natural Science)2017455):741. [百度学术