CN111368478B - 一种基于滚动直线导轨可动结合部的模态参数识别方法 - Google Patents
一种基于滚动直线导轨可动结合部的模态参数识别方法 Download PDFInfo
- Publication number
- CN111368478B CN111368478B CN202010149334.0A CN202010149334A CN111368478B CN 111368478 B CN111368478 B CN 111368478B CN 202010149334 A CN202010149334 A CN 202010149334A CN 111368478 B CN111368478 B CN 111368478B
- Authority
- CN
- China
- Prior art keywords
- guide rail
- test
- matrix
- node
- model
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 78
- 238000005096 rolling process Methods 0.000 title claims abstract description 50
- 238000012360 testing method Methods 0.000 claims abstract description 152
- 238000004458 analytical method Methods 0.000 claims abstract description 25
- 238000006073 displacement reaction Methods 0.000 claims abstract description 16
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 238000012795 verification Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 115
- 230000005284 excitation Effects 0.000 claims description 46
- 238000005316 response function Methods 0.000 claims description 37
- 238000013016 damping Methods 0.000 claims description 22
- 239000013598 vector Substances 0.000 claims description 18
- 230000004044 response Effects 0.000 claims description 13
- 238000012545 processing Methods 0.000 claims description 12
- 238000001228 spectrum Methods 0.000 claims description 8
- 230000001133 acceleration Effects 0.000 claims description 7
- 238000004088 simulation Methods 0.000 claims description 7
- 238000010998 test method Methods 0.000 claims description 7
- 229910000831 Steel Inorganic materials 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- 239000010959 steel Substances 0.000 claims description 5
- 239000000463 material Substances 0.000 claims description 4
- XAGFODPZIPBFFR-UHFFFAOYSA-N aluminium Chemical compound [Al] XAGFODPZIPBFFR-UHFFFAOYSA-N 0.000 claims description 3
- 229910052782 aluminium Inorganic materials 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 230000008676 import Effects 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000035945 sensitivity Effects 0.000 claims description 3
- 239000000725 suspension Substances 0.000 claims description 3
- 241000764238 Isis Species 0.000 claims 1
- 230000033001 locomotion Effects 0.000 abstract description 27
- 238000002474 experimental method Methods 0.000 description 11
- 230000008878 coupling Effects 0.000 description 5
- 238000010168 coupling process Methods 0.000 description 5
- 238000005859 coupling reaction Methods 0.000 description 5
- 238000013519 translation Methods 0.000 description 5
- 230000003595 spectral effect Effects 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000010354 integration Effects 0.000 description 3
- 238000005452 bending Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000000695 excitation spectrum Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 230000000903 blocking effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- 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
Abstract
本发明涉及机床整机动态特性分析的技术领域,公开了一种基于滚动直线导轨可动结合部的模态参数识别方法及其验证方法与测试系统,模态参数识别方法包括:步骤1、构建模拟直线滚动导轨的计算模型,计算模型包括滑块部、导轨部和导轨与滑块之间的结合部,结合部为三维空间八节点六面体单元的有限元;步骤2、建立节点位移与节点受力之间的关系,构建力学模型;步骤3、利用MATLAB中的lsqnonlin对力学模型的结合部的动力学参数求解。对应该方法还有相应有效性的验证方法及测试系统。本发明提供了适用于导轨结合部的八节点六面体模型,该模型将每个节点沿着导轨运动方向的自由度释放,所建立的模型与导轨的实际运动及受力情况相吻合。
Description
技术领域
本发明属于机床整机动态特性分析的技术领域,尤其是一种基于滚动直线导轨可动结合部的模态参数识别方法
背景技术
如中国专利公告号CN 103267621公开了一种基于滚动导轨系统虚拟材料层参数的识别方法,它将滚动导轨结合部看做是一种虚拟材料层,并采用模态试验和CAE建模分析实验相结合的方法,优化得到虚拟材料层的弹性模量、泊松比、剪切模量、虚拟材料层厚度和面积等参数。其并未提出适用于虚拟材料层,即可动结合部的力学模型,合适的力学模型使有限元分析更为准确。
发明内容
本发明的目的在于提供一种基于滚动直线导轨可动结合部的模态参数识别方法。
为解决上述问题而采用的一种基于滚动直线导轨可动结合部的模态参数识别方法,包括以下步骤:
步骤1、构建模拟直线滚动导轨的计算模型,计算模型包括滑块部、导轨部和导轨与滑块之间的结合部,滑块部下面为结合部,结合部下面为导轨部,结合部为三维空间八节点六面体单元的有限元,沿导轨方向视为约束,其中,每个节点具有2个自由度,一个三维空间八节点六面体单元共计16个自由度;由有限单元理论可知,三维实体单元,每节点只要三个自由度,则该有限单元模型,就可完全反映所研究三维结构的所有运动形式(平动、扭转、弯曲)。因此用一个仿照有限元“八节点六面体”的单元,来综合模拟直线滚动导轨可动结合部。
步骤2、建立节点位移与节点受力之间的关系,构建力学模型,力学模型采用基于柔度影响系数的方法完成。这种建模方法将一维柔度影响系数推广到三维柔度影响系数,不仅具有弹性-阻尼模型的优点而且还考虑了结合部单元各节点之间的相互耦合,具有相对较好的精度。
步骤3、利用MATLAB中的lsqnonlin对力学模型的结合部的动力学参数求解。
进一步地,步骤1中结合部为没有质量属性只有刚度和阻尼属性的单元。
进一步地,通过步骤1中节点的位移与其受力之间的关系,构建结合部的力学模型,力学模型构建采用基于柔度影响系数的方法完成。
进一步地,节点的位移用八节点六面体单元上相应的两结点位移之差描述,设各节点位移为xij,i=1,2,…,7,8;j=1,2,其中j=1代表X方向的平动,j=2代表Z方向的平动。
进一步地,节点之间的受力用柔性弹簧简单描述,设各节点所受的力为fij。
进一步地,推导结合部有限单元的刚度矩阵参数
为刚度影响系数,表示在节点m与节点(m+4)的n方向产生单位相对位移,所需要在i节点j方向施加的力,其中/>表示节点号,/>表示方向;
结合部的运动特性通过节点1与节点5、节点2与节点6、节点3与节点7、节点4与节点8之间的相对运动表现出来,而这些节点之间的相对运动可表示为:x1j-x5j,x2j-x6j,x3j-x7j,x4j-x8j,j=1,2;
平衡条件下相对应的两节点所受的力大小相等方向相反,因此f1j=-f5j,f2j=-f6j,f3j=-f7j,f4j=-f8j;
用矩阵描述节点力:{F}=(f11,f12,f21,f22.....f72,f81,f82)T (1)
同样定义{X}=(x11,x12,x21,x22.....x72,x81,x82)T,式(1)可描述为矩阵形式:
[K]×{X}={F} (2)
其中,[K]16×16为结合部单元的刚度矩阵,并且具有以下两种性质:
①对称性;
②可分块性;
[k']为8×8矩阵;
[k']8×8矩阵的具体表达式为:
确定[k']8×8里面的45个独立变量的值,以确定结合部(III)单元的刚度矩阵。在建立的滚动直线导轨副可动结合部单元的动力学模型中,结合部单元刚度矩阵的参数是未知的,需要寻求一种合适的且精度高的方法对其进行参数识别。对[k']8×8里面的参数进行研究发现里面共有45个独立变量需要识别,只要确定了这45个值,就可以确定结合部的单元矩阵。
进一步地,确定[k']8×8里面的45个独立变量的值的方法为非线性最小二乘的迭代算法:
有阻尼多自由度线性振动系统的运动微分方程为:
其傅氏变换为:
-w2[M]{X(w)}+jw[C]{X(w)}+[K]{X(w)}={F(w)} (5)
为计算方便,定义[Z(w)]=[K]+jw[C]-w2[M],则式(5)可转换为:
[Z(w)]{X(w)}={F(w)} (6)
对于一个多自由度系统,在第j个自由度上进行激励,而在第i个自由度上测量响应,便可得到频响函数Hij(w);
当i、j遍历1,2,……n时,共可得到n2个频响函数,他们构成一个n×n的矩阵,即频响函数矩阵
对系统各个自由度上的激励谱和响应谱分别记为Fi(w)、Xi(w),则它们分别构成激励谱向量{F(w)}和响应谱向量{X(w)},其中
那么系统的频响函数H(w)可以表示为
或{X(w)}=[H(w)]{F(w)} (7)
将(7)代入(6)得,
[Z(w)][H(w)]=I (8)
固定激励点j,移动响应点i,采用单点激振、多点拾振(SIMO)方法进行锤击实验测试时,可以测得各响应点i关于固定激励点j,即j固定,i=1,2,……n时的频响函数Hij(w),即为得到频响函数矩阵的第j列;
定义列矩阵
则式(8)又可写为:
[Z(w)][Hj(w)]=Ij (9)
其中,Ij为单位矩阵的第j列,j恒定展开得到
通过锤击实验,以外界激励频率w为参变量获得频响函数列向量,进而求出整个系统m,c,k中的未知参数。从各矩阵的表达式可以看出,矩阵[Z(w)]与m,c,k有关,因此频响函数[H(w)]也与m,c,k相关。频响函数是反映系统固有特性的量,是以外界激励频率w为参变量的非参数模型,可以借助锤击实验获得,这是锤击法进行模态实验的理论依据。
进一步地,整个系统包括结合部的模拟单元,进一步得到含有可动结合部的整个系统的机械位移阻抗;
整个系统包括模型滑块,模型结合部,模型导轨,模型滑块与模型导轨之间为模型结合部;
设模型滑块和模型导轨子结构的动力学方程为:
包含未知动力学参数的模型结合部的动力学方程为:
其中,[Ms],[Cs],[Ks]分别为两个子结构模型滑块、模型导轨的整体质量、阻尼、刚度矩阵;[Mj],[Cj],[Kj]分别为反映模型结合部动力学特性的质量、阻尼、刚度矩阵的未知量,[Cj]按粘性比例阻尼处理,[Cj]=α[Mj]+β[Kj];
通过有限元理论整体刚度矩阵的集成规则,可对(11)和(12)两个方程进行组装,得到整体结构的动力学方程:
子结构模型滑块和模型导轨本身的阻尼相对于模型结合部的阻尼可以忽略不计,即Cs<<Cj,并且模型结合面(2)质量亦可忽略,则[Cj]=β[Kj],式(13)可以进一步简化为:
式(14)进行傅氏变换:
-w2[Ms]{X(w)}+jwβ[Kj]{X(w)}+[Ks+Kj]{X(w)}={F(w)} (15)
引入频响函数,进一步得到:
(-w2[Ms]+jwβ[Kj]+[Ks+Kj]){H(w)}=[I] (16)
其中,机械结构的质量矩阵[MS]和刚度矩阵[KS]都可通过有限元模态分析获得,求出模型结合部单元的刚度矩阵[Kj]。
进一步地,通过有限元理论整体刚度矩阵的集成规则进行矩阵组装及参数求解的方法为:
将结合部的八结点六面体模拟单元的八个结点的编号放在整个有限元模型的最前面,即Node 1~Node 8;
这八个结点每个结点只有两个自由度,而机械结构的每一个结点都有三个自由度,先把16×16的矩阵,扩充为24×24的矩阵,其中被约束方向的自由度用0表示,然后再把此24×24的矩阵扩充为与子结构刚度矩阵[KS]同型,进行整个结构矩阵的组装;
结合部矩阵扩充为矩阵KJ:
即:
其中,n=24为整个结构自由度总数,KJ为结合面刚度刚度矩阵Kj扩充后的矩阵;
当采用单点激振,多点拾振方法时,由式(16)可知
(-w2[Ms]+jwβ[Kj]+[Ks+Kj]){Hj(w)}=[Ij] (17)
其中,{Hj(w)}、[Ij]分别为{H(w)}、[I]的第j列;
展开、整理得到:
式(18)中,左侧为未知量,右侧为常数列向量,定义
因结合面刚度矩阵的对称性原理其中K1为8×8的对称矩阵;
取式(19)的第一行进行分析,并转为标准形式:
此方程中共有9个未知参数,包括一个比例阻尼系数β和8个刚度参数,求解出这9个参数的值。
进一步地,求解式(20)中的9个未知参数值的方法为:
取不同的w值,得到不同的频响函数矩阵,从而得到不同的方程,可以表示为:
其中p为w取值的个数,也是方程的个数;
得到了p个方程之后,其中p≥9,就可以求解出结合部(II)单元刚度矩阵的第一行的8个参数;
然后再取出式(19)的第二行,按照第一行参数相同的求解方法进行求解,只是由于刚度矩阵是对称阵,第二行的未知参数就只有7个,以此类推就可以求出刚度矩阵中的所有未知参数,从而得到结合部的刚度矩阵。
进一步地,对结合部的刚度矩阵中的所求解的未知参数进行优化的方法为:
根据(20)式,定义方程fp(k):
由于阻尼系数β的存在,使得fp(k)为非线性函数,采用非线性最小二乘法(lsqnonlin)对参数进行优化,进而求解最优动力学参数;
进一步地,求解最优动力学参数的方法为:
通过迭代算法,即从被估计参数的某一初值开始,产生一系列的参数估计值,如果这个估计参数序列使目标函数收敛到极小,则认为寻找到最优解。
进一步地,寻找最优解的具体方法为:
通过MATLAB Optimization Toolbox中的非线性最小二乘法来进行参数的优化估计;
由于模态试验得到的频率响应函数是复数,所以在参数求解时非线性函数fp(k)取模(abs),目标函数为:
f(k)=[abs(f1(k))abs(f2(k))abs(f3(k))…abs(fp(k))]T (23)
在MATLAB软件中使用最小二乘法求解f(k)时,目标函数表达为:
其中k为未知列向量,f(k)为函数列向量;
利用MATLAB中的lsqnonlin对结合部动力学参数求解;
进一步地,利用MATLAB中的lsqnonlin对结合部动力学参数求解的流程步骤为:
步骤1:提取滚动导轨可动结合部的频率响应函数的列矩阵;
步骤2:选择属于可动结合部的多个列向量以及提取轨道和滑块的质量矩阵和刚度矩阵和计算出可动结合部刚度矩阵初值导入算式;
步骤3:利用常数矩阵,系数矩阵,目标函数表达式,刚度矩阵的边界条件得到参数;
步骤4:确定参数值;
步骤5:将参数值插入有限元模型;
步骤6:进行有限元模型的模态试验;
步骤7:得到限元模型的分析结果。
关于用于验证基于滚动直线导轨可动结合部的模态参数识别方法有效性的方法,包括以下步骤:
步骤1、采用单点激励多点响应的模态试验方法来获取其整体结构的频率响应函数矩阵中的列向量。
步骤2、提取整体结构,即包含结合部的可靠、高精度的频响函数矩阵,采集结合部单元刚度矩阵的动力学参数的原始数据。
步骤3、得到反映结合部特性的整体结构的前几阶固有频率和相应的模态振型,进行参数识别结果的有限元仿真比较验证。
进一步地,用于步骤1中的滚动直线导轨副的模态试验所使用的导轨副试件包括试件配重、导轨副滑块、导轨副的结合部,导轨副导轨和试件底座,试件配重置于导轨副滑块上,导轨副滑块置于导轨副导轨之上,导轨副滑块与导轨副导轨之间有一导轨副的结合部,导轨副滑块与导轨副导轨滑动连接,导轨副导轨下面连接试件底座,试件底座的大小、材料都与试件配重一样。样的结构不仅能充分的反映结合部的特征,而且试件结构简单,便于试验。
进一步地,导轨副试件被悬挂放置,导轨副的两边不与导轨副试件接触,以模拟导轨副试件的自由几何边界条件。
进一步地,导轨副试件被钢丝绳悬挂,一般钢丝绳在滚动直线导轨副试件上悬挂的位置选择在振幅较小的位置,最佳悬挂点应该是某阶振型的节点处。
进一步地,用于步骤1中模态试验方法的测试模型包括多个测试点,测试点的多少与分布应该根据实体模型和试验需求决定。
进一步地,测试点的布置方法为:
把测试点在试件上按照一定程度地均匀分布,同时可以在关键区域布置更多的测试点,测试点不靠近节点。试验模型测点越多,越能充分表征结构的固有模态属性。不过试验测点的过度增多对表征结构的固有模态属性并没有太大的改善,反而增加了测试的困难和任务量,同时测试点不应该靠近节点,这样的测试点得到的信息才会有更高的信噪比。
进一步地,用于步骤1中模态试验方法的激励点选择为正好落在某阶模态的反节点或者附近。
进一步地,步骤1中模态试验的具体方法为:试验采用铝制锤头,灵敏度设置为0.25mV/N;带宽fmax取2048Hz,采样频率即为fs=2fmax取4096Hz,谱线数L取2048,此时频率分辨率Δf=fmax/L;
试验采取在锤击点的两个方向进行激振,即-Z和-X方向。
用力锤锤击试验试件获取激励,通过传感器连接试验试件将激励信息传导进入LMS系统,再将LMS系统接入试验电脑,通过试验电脑处理并导出数据信息。
关于基于滚动直线导轨可动结合部的模态参数识别方法有效性的验证方法的模态测试系统,包括激振系统、拾振系统、数据采集与处理系统,三者电性连接,数据采集与处理系统可以接收并处理激振系统和拾振系统传递的信息,数据采集与处理系统使用的是LMS Test.Lab模态试验测试与分析系统,激振系统为力锤激振系统,力锤激振系统使用的是脉冲力锤,拾振系统使用的是ICP三向加速度传感器。
进一步地,力锤激振系统包括力锤传感器和第一电荷放大器,拾振系统包括加速度传感器和第二电荷放大器,第一电荷放大器电性连接在力锤传感器和频谱分析仪之间,第二电荷放大器电性连接在加速度传感器和频谱分析仪之间,频谱分析仪与数据采集和数据处理系统电性连接。
本发明提供了适用于导轨结合部的八节点六面体模型,该模型在充分考虑结合部各个节点之间的耦合关系的同时,将每个节点沿着导轨运动方向的自由度释放,所建立的模型与导轨的实际运动及受力情况相吻合。并根据所建立的模型,提出了相应的参数识别方法。通过搭建导轨结合部的实验平台,进行参数识别并验证了模型的有效性。最后,采用所提出的方法,可以建立不同型号的滚动导轨结合部动力学模型参数库,为滚动导轨结合部动力学模型在机床动力学建模中的应用奠定基础。
附图说明
图1为导轨副可能的受力状态示意图。
图2为八节点六面体结合部单元示意图。
图3为结合部单元的等效动力学模型示意图。
图4为导轨可动结合部示意图。
图5为参数识别流程图。
图6为模态试验测试系统试验图。
图7为滚动直线导轨副模态试验试件截面图。
图8为LMS模态测试模型。
图9为实验模态与计算模态的比较示意图表。
附图标记:1I、滑块部;1II、导轨部;1III、结合部;1、模型滑块;2、模型结合部;3、模型导轨;1a、试件配重;2a、导轨副滑块;3a、导轨副的结合部;4a、导轨副导轨;5a、试件底座。
具体实施方式
一种基于滚动直线导轨可动结合部的模态参数识别方法,如图1所示,零部件在空间的六种受力状态分别是沿X、Y、Z轴线方向的力以及绕三个轴线方向的力矩。由于滑块在沿导轨方向是工作方向,因此在受力分析时不在此方向上加力,视为在此方向上固定约束。所以滚动直线导轨副可以有五种受力状态。
在这五种受力状态下滚动直线导轨副会产生相应的形变运动,分别为俯仰运动、侧翻运动、偏航运动、上下运动和左右运动。当然其可动结合部也跟着整个导轨副产生相应的变形运动,而且变形量更大,因为可动结合滚珠与沟槽接触处为点接触,其刚度远小于导轨和滑块材料自身的刚度。由有限单元理论可知,三维实体单元,每节点只要三个自由度,则该有限单元模型,就可完全反映所研究三维结构的所有运动形式(平动、扭转、弯曲)。因此用一个仿照有限元“八节点六面体”的单元,来综合模拟直线滚动导轨可动结合部。
如图2所示,滑块1I在导轨1II的上面,导轨1II与滑块1I之间为结合部1III单元,而且这种结合部1III单元没有质量属性只有刚度和阻尼属性。每一个结合部1III单元有8个节点,每个节点具有2个自由度(因为在沿导轨方向视为约束),因此一个单元共计16个自由度。结合部单元属性可以通过节点1与5、节点2与6、节点3与7、节点4与8之间的相对运动表现出来。因此只要能准确建立这些节点位移与节点受力之间的关系,就等于建立了其力学模型。
结合部结点位移可用两结合部单元上相应的结点位移之差描述,其模型构建采用基于柔度影响系数的方法完成。这种建模方法将一维柔度影响系数推广到三维柔度影响系数,不仅具有弹性-阻尼模型的优点而且还考虑了结合部单元各节点之间的相互耦合,具有相对较好的精度。
如图3所示,该建模方法在滚动直线导轨可动结合部上的应用,以及各节点之间的受力用柔性弹簧简单描述,结合部单元1III只有刚度和阻尼属性无质量属性。设各节点位移为xij,各节点所受的力为fij,i=1,2,…,7,8;j=1,2,其中j=1代表X方向的平动,j=2代表Z方向的平动。
首先推导结合部有限单元的刚度矩阵。如前所述,结合部的运动特性通过节点1与节点5、节点2与节点6、节点3与节点7、节点4与节点8之间的相对运动表现出来,而这些节点之间的相对运动可表示为:x1j-x5j,x2j-x6j,x3j-x7j,x4j-x8j,j=1,2。依据刚度影响系数法的物理意义,则有:
其中,为刚度影响系数,其中/>表示节点号,/>表示方向。的物理意义:仅在节点m与节点(m+4)的n方向产生单位相对位移,而需要在i节点j方向所需施加的力。fij表示节点力,平衡条件下相对应的两节点所受的力大小相等方向相反,因此f1j=-f5j,f2j=-f6j,f3j=-f7j,f4j=-f8j,用矩阵描述节点力:{F}=(f11,f12,f21,f22.....f72,f81,f82)T (1)
同样定义:{X}=(x11,x12,x21,x22.....x72,x81,x82)T,式(1)可描述为矩阵形式:
[K]×{X}={F} (2)
其中[K]16×16为结合部单元的刚度矩阵,并且具有以下两种性质:
①对称性。
②可分块性,[k']为8×8矩阵。
因此只需要得到[k']8×8的矩阵就能得到刚度矩阵[K]16×16,而[k']8×8矩阵的具体表达式为:
在建立的滚动直线导轨副可动结合部单元的动力学模型中,结合部单元刚度矩阵的参数是未知的,需要寻求一种合适的且精度高的方法对其进行参数识别。对[k']8×8里面的参数进行研究发现里面共有45个独立变量需要识别,只要确定了这45个值,就可以确定结合部的单元矩阵。
有关固定机械结合部动力学参数识别的方法有很多,其中一种基于系统动力学矩阵和频率响应函数的参数识别方法,它综合考虑了结合部的多种特性,能够准确的反映结合部的相关特性,并且该方法已经在机床螺栓固定结合部的参数识别中取得了很好的识别效果,基于上述基础,本发明结合实验模态分析所得的实验值(频响函数列)与有限元分析的理论值(实体子结构的质量、刚度矩阵),以未知的结合面参数(结合面刚度矩阵、阻尼矩阵)为设计变量,以位移阻抗矩阵和位移频响函数之积与单位矩阵之差最小为优化设计目标,通过非线性最小二乘的迭代算法求解结合面动力学参数。
有阻尼多自由度线性振动系统的运动微分方程为:
其傅氏变换为:
-w2[M]{X(w)}+jw[C]{X(w)}+[K]{X(w)}={F(w)} (5)
为计算方便,定义[Z(w)]=[K]+jw[C]-w2[M],则式(5)可转换为:
[Z(w)]{X(w)}={F(w)} (6)
对于一个多自由度系统,仅在第j个自由度上进行激励,而在第i个自由度上测量响应,便可得到频响函数Hij(w);当i、j遍历1,2,……n时,共可得到n2个频响函数,他们构成一个n×n的矩阵,即频响函数矩阵
如果对系统各个自由度上的激励谱和响应谱分别记为Fi(w)、Xi(w),则它们分别构成激励谱向量{F(w)}和响应谱向量{X(w)},其中
那么系统的频响函数H(w)可以表示为
或{X(w)}=[H(w)]{F(w)} (7)
将(7)代入(6)得,
[Z(w)][H(w)]=I (8)
根据上述频响函数的物理定义可知:当我们固定激励点j,移动响应点i,采用单点激振、多点拾振(SIMO)方法进行锤击实验测试时,可以测得各响应点i关于固定激励点j的频响函数Hij(w)(j固定,i=1,2,……n),也即可以得到频响函数矩阵的第j列。不妨定义列矩阵
则式(8)又可写为:
[Z(w)][Hj(w)]=Ij (9)
其中,Ij为单位矩阵的第j列,j恒定展开得到
从各矩阵的表达式可以看出,矩阵[Z(w)]与m,c,k有关,因此频响函数[H(w)]也与m,c,k相关。频响函数是反映系统固有特性的量,是以外界激励频率w为参变量的非参数模型,可以借助锤击实验获得,这是锤击法进行模态实验的理论依据。
根据上一节的参数识别方法推导可知,通过锤击试验得到整个结构系统的频响函数列向量,进而求出整个系统m,c,k中的一些未知参数。而整个系统也包括了结合部的模拟单元,因此需要先做的是要得到含有可动结合部的整个系统的机械位移阻抗,如图4所示,模型滑块1与模型导轨3之间为可动模型结合部2。
设模型滑块1和模型导轨3子结构的动力学方程为:
包含未知动力学参数的结合部2的动力学方程为:
其中,[Ms],[Cs],[Ks]分别为两个子结构模型滑块1和模型导轨3的整体质量、阻尼、刚度矩阵;[Mj],[Cj],[Kj]分别为反映结合面动力学特性的质量、阻尼、刚度矩阵(未知量),[Cj]按粘性比例阻尼处理,[Cj]=α[Mj]+β[Kj]。
按照有限元理论整体刚度矩阵的集成规则,可对(11)和(12)两个方程进行组装,得到整体结构的动力学方程:
子结构模型滑块1、模型导轨3本身的阻尼相对于模型结合部2的阻尼可以忽略不计,即Cs<<Cj,并且模型结合部2质量亦可忽略,则[Cj]=β[Kj],式(13)可以进一步简化为:
式(14)进行傅氏变换:
-w2[Ms]{X(w)}+jwβ[Kj]{X(w)}+[Ks+Kj]{X(w)}={F(w)} (15)
引入频响函数,进一步得到:
(-w2[Ms]+jwβ[Kj]+[Ks+Kj]){H(w)}=[I] (16)
其中机械结构的质量矩阵[MS]和刚度矩阵[KS]都可通过有限元模态分析获得,未知动力学参数只有结合部(1III)单元的刚度矩阵[Kj],这就是我们所需要求的动力学参数。为了方便与机械结构的矩阵组装及参数求解,本发明把结合部的八结点六面体模拟单元的八个结点的编号放在整个有限元模型的最前面即Node 1~Node 8。同时因为这八个结点每个结点只有两个自由度,而机械结构的每一个结点都有三个自由度,为了矩阵的组装方便,先把16×16的矩阵,扩充为24×24的矩阵,其中被约束方向的自由度用0表示,然后再把此24×24的矩阵扩充为与子结构刚度矩阵[KS]同型,方便整个结构矩阵的组装。结合部矩阵扩充为矩阵KJ:
即:
其中,n=24为整个结构自由度总数,KJ为结合面刚度刚度矩阵Kj扩充后的矩阵。
当采用单点激振,多点拾振方法时,由式(16)可知
(-w2[Ms]+jwβ[Kj]+[Ks+Kj]){Hj(w)}=[Ij] (17)
其中,{Hj(w)}、[Ij]分别为{H(w)}、[I]的第j列。
展开、整理得到:
式18中,左侧为未知量,右侧为常数列向量,定义
因结合面刚度矩阵的对称性原理,其中K1为8×8的对称矩阵。取式(19)的第一行进行分析,并转为标准形式:
此方程中共有9个未知参数,包括一个比例阻尼系数β和8个刚度参数。如果能求解出这9个参数的值,就至少需要9个等式方程。所以我们取不同的w值,得到不同的频响函数矩阵,从而得到不同的方程。可以表示为:
其中p为w取值的个数,也是方程的个数。得到了p(p≥9)个方程之后,就可以求解出结合部单元刚度矩阵的第一行的8个参数。然后再取出式(19)的第二行,按照第一行参数相同的求解方法进行求解,只是由于刚度矩阵是对称阵,第二行的未知参数就只有7个,以此类推就可以求出刚度矩阵中的所有未知参数,从而得到结合部1III的刚度矩阵。
建立好方程等式之后对参数的识别过程就是被估参数的优化过程。根据(20)式,定义方程fp(k):
由于阻尼系数β的存在,使得fp(k)为非线性函数,采用非线性最小二乘法(lsqnonlin)对参数进行优化,进而求解最优动力学参数。非线性最小二乘法是以误差的平方和最小为准则来优化非线性模型参数的一种参数估计方法。由于非线性所以不能像线性最小二乘法用偏导数的方法求参数估计值,因此本发明采用的是迭代算法,即从被估计参数的某一初值开始,产生一系列的参数估计值,如果这个估计参数序列使目标函数收敛到极小,则认为寻找到最优解。
基于上述理论分析,本发明在参数求解时选择MATLAB Optimization Toolbox中的非线性最小二乘法(lsqnonlin)来进行参数的优化估计。由于模态试验得到的频率响应函数是复数,所以在参数求解时非线性函数fp(k)取模(abs),目标函数为:
f(k)=[abs(f1(k))abs(f2(k))abs(f3(k))…abs(fp(k))]T (23)
在MATLAB软件中使用最小二乘法求解f(k)时,目标函数表达为:
其中k为未知列向量,f(k)为函数列向量。
利用MATLAB中的lsqnonlin对结合部动力学参数求解的流程图如图5所示,其步骤为:
步骤1:提取滚动导轨可动结合部的频率响应函数的列矩阵;
步骤2:选择属于可动结合部的多个列向量以及提取轨道和滑块的质量矩阵和刚度矩阵和计算出可动结合部刚度矩阵初值导入算式;
步骤3:利用常数矩阵,系数矩阵,目标函数表达式,刚度矩阵的边界条件得到参数;
步骤4:确定参数值;
步骤5:将参数值插入有限元模型;
步骤6:进行有限元模型的模态试验;
步骤7:得到限元模型的分析结果。
基于滚动直线导轨副的模态试验还使用了参数识别验证有效性的方法:
上述对滚动直线导轨副的可动结合部参数识别的方法中,为了求解出结合部1III单元刚度矩阵的未知参数,我们需要得到其整体结构的频率响应函数矩阵中的列向量,因此我们采用单点激励多点响应的(SIMO)模态试验方法来获取。
本发明对滚动直线导轨副进行模态试验的目的概括为两个方面:其一为提取整体结构(包含结合部)的可靠、高精度的频响函数矩阵,为结合部单元刚度矩阵的动力学参数识别提供原始数据。其二为得到反映结合部特性的整体结构的前几阶固有频率和相应的模态振型,为基于参数识别结果的有限元仿真验证提供比较基准。如图6所示,模态试验测试系统由3部分组成,即激振系统、拾振系统、数据采集与处理系统。
实验过程中,数据采集与处理系统使用的是比利时的LMS Test.Lab模态试验测试与分析系统,力锤激振系统使用的是美国PCB公司生产的086C04型脉冲力锤,拾振系统使用的是美国PCB公司生产的356A16型ICP三向加速度传感器。滚动直线导轨副的模态试验为了尽量减少其它因素的影响都是直接在导轨副的滑块上进行测试,但是由于滚动直线导轨副滑块上的面积有限,测点比较少而且不便于试验,难以充分的表征试件结构的固有模态属性,特别是结合部的变形形态。本实验为加大滑块及导轨的面积,在滑块及导轨上连接相应配重,其试件设计如图7所示,试件配重1a通过4个螺栓与滚动直线导轨副滑块2a紧密相连,导轨副结合部3a位于导轨副导轨4a与导轨副滑块2a之间,导轨副导轨4a安装在试件底座5a上面,试件底座5a的大小、材料都与试件配重1a一样,通过螺栓与导轨副导轨4a紧密相连。这样的结构不仅能充分的反映结合部的特征,而且试件结构简单,便于试验。
在试验时为了更好的排除外界因素对滚动直线导轨副试件模态测试的干扰,采取把导轨副试件用钢丝绳悬挂起来(导轨副的两边不与导轨副试件接触),来模拟导轨副试件的自由几何边界条件。一般钢丝绳在滚动直线导轨副试件上悬挂的位置选择在振幅较小的位置,最佳悬挂点应该是某阶振型的节点处。
基于上述验证方法的LMS Test.Lab模态试验测试与分析系统,在进行模态分析试验时,首先需要建立测试模型。测试模型由许多个测试点构成,测试点的多少与分布应该根据实体模型和试验需求决定。一般情况下,比较好的方法是把测试点在试件上按照某种程度地均匀分布,同时可以在关键区域更多的布置测试点。试验模型测点越多,越能充分表征结构的固有模态属性。不过试验测点的过度增多对表征结构的固有模态属性并没有太大的改善,反而增加了测试的困难和任务量,同时测试点不应该靠近节点,这样的测试点得到的信息才会有更高的信噪比。本文在进行测试时,导轨及其底座总共布置了56个测试点,滑块及其配重总共布置了52个测试点,整个滚动直线导轨副试件布置了108个测试点,但是由于滑块及其配重组件的最下面的4个点无法测量,所以总共测试的点数为104个,其LMS模态测试模型如图8所示。
建立好模态试验的测试模型后,下一步是进行试验的参数设置和激励点的选择。激励点选择的原则是以能有效激起各阶模态,若激振点正好落在某阶模态的反节点或者附近,则激励力就能有效地激起该阶的模态。试验采用铝制锤头,灵敏度设置为0.25mV/N;带宽fmax取2048Hz,采样频率即为fs=2fmax取4096Hz,谱线数L取2048,此时频率分辨率Δf=fmax/L。为了能激发出滚动直线导轨副更多的模态,试验采取在锤击点的两个方向进行激振,-Z和-X方向。
完成整个滚动直线导轨副试验104个点的锤击模态测试之后,用LMS软件包自带的数据处理分析系统Modal Analysis进行模态分析,以LGS20H为例,其试验模态分析结果如图9所示。
通过对滚动直线导轨副的锤击模态试验得到了四阶模态,分别为第一阶侧翻运动模态,固有频率为225.989Hz;第二阶偏航运动模态,固有频率为:468.007Hz;第三阶俯仰运动模态,固有频率为:616.490Hz;第四阶上下运动模态,固有频率为:1392.187Hz。模态测试完成后,提取相应的频响函数列进行可动结合部单元刚度矩阵的识别,按照图5所示的参数识别流程图进行可动结合部单元刚度矩阵的识别,其识别结果为:
由于Y方向为滚动直线导轨副的工作运动方向,在分析时视为固定约束,所以在刚度矩阵中其值为0。
表1 16组导轨副的试验和计算结果
/>
其他类型的滚动直线导轨副按照上述处理方法,分别得到了模态试验的结果和参数识别有限元计算的结果,如表1所示。根据上面的试验结果和计算结果可以看出,平均误差在10%左右,计算结果和试验结果比较吻合。但少数试验误差相对较大,如LG25H的第一阶侧翻模态、LG30的第三阶俯仰模态、LG35的第三阶俯仰模态等。分析其原因主要是可能是其一为试件中除了有滚动导轨副的可动结合部之外,导轨与底座、滑块与配重之间两个固定结合部的存在对其影响较大,同时由于导轨底面和滑块顶面的都有凹槽,使得和底座和配置并不能完全的接触,其影响就更大了。其二是由于直线滚动导轨副的高度较小(底座和配重之间的距离),使得在测试可动结合部的点的时候,传感器不能贴到与模型相应的准确位置,只能在测试点的附近位置,从而产生了测试误差,从而对结果产生了影响。其三是由于滑块几何尺寸太小(无法贴传感器),因而无法得到该部分的频响函数,所以本文在测试模型和有限元模型中都忽略了这部分,而实测信号包含了整个滚动直线导轨副试件,从而对识别结果产生了影响。
以上对本发明的具体实施实例进行了详细的阐述,但本发明并不限制与以上描述的具体实施例,其只作为范例。任何等同修改和替代也都在本发明的范畴之中。因此,在不脱离本发明的精神和范围下所作出的均等变换和修改,都应涵盖在本发明的范围内。
本发明提供了适用于导轨结合部的八节点六面体模型,该模型在充分考虑结合部各个节点之间的耦合关系的同时,将每个节点沿着导轨运动方向的自由度释放,所建立的模型与导轨的实际运动及受力情况相吻合。并根据所建立的模型,提出了相应的参数识别方法。通过搭建导轨结合部的实验平台,进行参数识别并验证了模型的有效性。最后,采用所提出的方法,可以建立不同型号的滚动导轨结合部动力学模型参数库,为滚动导轨结合部动力学模型在机床动力学建模中的应用奠定基础。
Claims (1)
1.一种基于滚动直线导轨可动结合部的模态参数识别方法,包括以下步骤:
步骤1、构建模拟直线滚动导轨的计算模型,所述计算模型包括滑块部(1I)、导轨部(1II)和导轨与滑块之间的结合部(1III),所述滑块部(1I)下面为结合部(1III),所述结合部(1III)下面为导轨部(1II),所述结合部(1III)为三维空间八节点六面体单元的有限元,所述结合部(1III)为没有质量属性只有刚度和阻尼属性的单元,沿导轨方向视为约束,其中,每个节点具有2个自由度,一个三维空间八节点六面体单元共计16个自由度;
步骤2、建立节点位移与节点受力之间的关系,构建力学模型,所述力学模型采用基于柔度影响系数的方法完成,所述节点的位移用所述八节点六面体单元上相应的两结点位移之差描述;
步骤3、利用MATLAB中的lsqnonlin对所述力学模型的结合部(1III)的动力学参数求解;
利用MATLAB中的lsqnonlin对结合部动力学参数求解的流程步骤为:
步骤(1):提取滚动导轨可动结合部的频率响应函数的列矩阵;
步骤(2):选择属于可动结合部的多个列向量以及提取轨道和滑块的质量矩阵和刚度矩阵和计算出可动结合部刚度矩阵初值导入算式;
步骤(3):利用常数矩阵,系数矩阵,目标函数表达式,刚度矩阵的边界条件得到参数;
步骤(4):确定参数值;
步骤(5):将参数值插入有限元模型;
步骤(6):进行有限元模型的模态试验;
步骤(7):得到限元模型的分析结果;
关于用于验证基于滚动直线导轨可动结合部的模态参数识别方法有效性的方法,包括以下步骤:
步骤a:采用单点激励多点响应的模态试验方法来获取其整体结构的频率响应函数矩阵中的列向量;
步骤b:提取整体结构,即包含结合部的可靠、高精度的频响函数矩阵,采集结合部单元刚度矩阵的动力学参数的原始数据;
步骤c:得到反映结合部特性的整体结构的前几阶固有频率和相应的模态振型,进行参数识别结果的有限元仿真比较验证;
用于步骤a中的滚动直线导轨副的模态试验所使用的导轨副试件包括试件配重、导轨副滑块、导轨副的结合部,导轨副导轨和试件底座,试件配重置于导轨副滑块上,导轨副滑块置于导轨副导轨之上,导轨副滑块与导轨副导轨之间有一导轨副的结合部,导轨副滑块与导轨副导轨滑动连接,导轨副导轨下面连接试件底座,试件底座的大小、材料都与试件配重一样;
导轨副试件被悬挂放置,导轨副的两边不与导轨副试件接触,以模拟导轨副试件的自由几何边界条件;
导轨副试件被钢丝绳悬挂,最佳悬挂点应该是某阶振型的节点处;
用于步骤a中模态试验方法的测试模型包括多个测试点,测试点的多少与分布应该根据实体模型和试验需求决定;
测试点的布置方法为:
把测试点在试件上按照一定程度地均匀分布,同时在关键区域布置更多的测试点,测试点不靠近节点,同时测试点不应该靠近节点;
用于步骤a中模态试验方法的激励点选择为正好落在某阶模态的反节点或者附近;
步骤a中模态试验的具体方法为:试验采用铝制锤头,灵敏度设置为0.25mV/N;带宽fmax取2048Hz,采样频率即为fs=2fmax取4096Hz,谱线数L取2048,此时频率分辨率Δf=fmax/L;
试验采取在锤击点的两个方向进行激振,即-Z和-X方向;
用力锤锤击试验试件获取激励,通过传感器连接试验试件将激励信息传导进入LMS系统,再将LMS系统接入试验电脑,通过试验电脑处理并导出数据信息;
关于基于滚动直线导轨可动结合部的模态参数识别方法有效性的验证方法的模态测试系统,包括激振系统、拾振系统、数据采集与处理系统,三者电性连接,数据采集与处理系统可以接收并处理激振系统和拾振系统传递的信息,数据采集与处理系统使用的是LMSTest.Lab模态试验测试与分析系统,激振系统为力锤激振系统,力锤激振系统使用的是脉冲力锤,拾振系统使用的是ICP三向加速度传感器;
力锤激振系统包括力锤传感器和第一电荷放大器,拾振系统包括加速度传感器和第二电荷放大器,第一电荷放大器电性连接在力锤传感器和频谱分析仪之间,第二电荷放大器电性连接在加速度传感器和频谱分析仪之间,频谱分析仪与数据采集和数据处理系统电性连接。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010149334.0A CN111368478B (zh) | 2020-03-04 | 2020-03-04 | 一种基于滚动直线导轨可动结合部的模态参数识别方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010149334.0A CN111368478B (zh) | 2020-03-04 | 2020-03-04 | 一种基于滚动直线导轨可动结合部的模态参数识别方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111368478A CN111368478A (zh) | 2020-07-03 |
CN111368478B true CN111368478B (zh) | 2023-12-15 |
Family
ID=71210330
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010149334.0A Active CN111368478B (zh) | 2020-03-04 | 2020-03-04 | 一种基于滚动直线导轨可动结合部的模态参数识别方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111368478B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112364547B (zh) * | 2020-12-03 | 2022-07-08 | 天津大学 | 一种机床整机动力学性能全域快速预估方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2007147634A (ja) * | 2006-12-22 | 2007-06-14 | Rikogaku Shinkokai | 振動解析の方法および装置ならびにコンピュータ読み取り可能な記録媒体 |
CN101804464A (zh) * | 2010-02-24 | 2010-08-18 | 华中科技大学 | 一种机床锥配合固定结合部动力学参数识别方法 |
CN106815407A (zh) * | 2016-12-22 | 2017-06-09 | 四川大学 | 一种数控机床整机动态性能优化方法 |
-
2020
- 2020-03-04 CN CN202010149334.0A patent/CN111368478B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2007147634A (ja) * | 2006-12-22 | 2007-06-14 | Rikogaku Shinkokai | 振動解析の方法および装置ならびにコンピュータ読み取り可能な記録媒体 |
CN101804464A (zh) * | 2010-02-24 | 2010-08-18 | 华中科技大学 | 一种机床锥配合固定结合部动力学参数识别方法 |
CN106815407A (zh) * | 2016-12-22 | 2017-06-09 | 四川大学 | 一种数控机床整机动态性能优化方法 |
Non-Patent Citations (1)
Title |
---|
吴鹏 ; 尹玲 ; 罗卫强 ; 朱睿 ; 梁振锋 ; 陈义 ; .基于动力学分析的滚动导轨可动结合部参数识别方法研究.东莞理工学院学报.2018,(05),第61-65页. * |
Also Published As
Publication number | Publication date |
---|---|
CN111368478A (zh) | 2020-07-03 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Law et al. | Moving force identification: optimal state estimation approach | |
CN107330176B (zh) | 一种基于结构模态估计的应变计与加速度计联合布设方法 | |
Kumar et al. | Spindle dynamics identification for receptance coupling substructure analysis | |
US20130297266A1 (en) | Method For Improving Determination Of Mode Shapes For A Mechanical Structure And Applications Hereof | |
CN109948207A (zh) | 一种航空发动机高压转子装配误差预测方法 | |
CN101458205B (zh) | 一种机床固定结合部动力学参数的识别方法 | |
JP7346462B2 (ja) | ターゲット物体との関係において複数のトランスデューサの空間的構成を決定する方法 | |
JP2015032295A (ja) | 減衰振動解析方法 | |
CN111368478B (zh) | 一种基于滚动直线导轨可动结合部的模态参数识别方法 | |
CN110008520B (zh) | 基于位移响应协方差参数和贝叶斯融合的结构损伤识别方法 | |
CN115270296A (zh) | 一种商用车驾驶室疲劳耐久性分析方法及系统 | |
CN109684663B (zh) | 铁路货车车体焊缝疲劳寿命的评估方法及装置、系统 | |
Westin et al. | Cable-pulley interaction with dynamic wrap angle using the absolute nodal coordinate formulation | |
Jeong et al. | Equivalent stiffness modeling of linear motion guideways for stage systems | |
Nandakumar et al. | Structural parameter identification using damped transfer matrix and state vectors | |
CN111209694B (zh) | 一种桁架结构刚度和轴力的结构识别方法 | |
CN114091300A (zh) | 一种滚珠丝杠进给系统滚动结合部动态特性参数识别方法 | |
CN109916584B (zh) | 基于土-结构-消能减震装置相互作用的子结构试验方法 | |
CN107967387A (zh) | 一种汽车球铰工作力矩的有限元设计方法 | |
Adams et al. | Validation principles of agricultural machine multibody dynamics models | |
Zhu et al. | Sensor placement optimization of vibration test on medium-speed mill | |
JP3847264B2 (ja) | 地震応答解析方法 | |
Hongxia et al. | Research of optimal placement of gearbox sensor based on particle swarm optimization | |
Zhang | Effect of damping coefficients on structural damage identification | |
Lin et al. | Damage assessment of seismically excited buildings through incomplete measurements |
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 |