响应耦合子结构方法（精选5篇）
响应耦合子结构方法 篇1
0引言
机床的切削颤振会降低被加工工件的表面质量并影响机床的使用寿命。目前避免切削颤振最有效的方法是通过加工系统的稳定性叶瓣图选取稳定的切削加工参数[1]。稳定性叶瓣图由加工系统中刀尖点的频响函数计算得到[1,2],刀尖点频响函数可通过实验模态测试确定,但在实际加工中经常更换刀具导致的重复模态测试不仅耗时,而且会引入人 为误差。针 对这个问 题,Schmitz等[3,4]提出了一种基于RCSA的刀尖点频响函数预测方法,该方法将整个机床划分为主轴-刀柄、刀具两个子结构,二者频响函数分别通过模态锤击实验确定和梁模型理论计算得到,采用最小二乘法识别出刀柄-刀具结合面的平移、旋转刚度和阻尼系数后,将两个子结构进行耦合预测出刀尖点的频响函数。Schmitz等[5]对上述方法进行了改进,提出了一种三子结构响应耦合分析法,将机床系统划分为主轴-刀柄(包括刀柄的圆锥部分和凸缘部分)、剩余刀柄、刀具三部分,通过对刀柄模型进行模态锤击实验确定其频响函数,最终通过RCSA法预测出刀尖点频响函数。近年来,很多学者在此基础上对该方法进行了改进,并研究了主轴结构中轴承的刚度、阻尼参数等对刀尖点频响函数预测结果的影响[6,7,8,9,10,11,12]。闫蓉等[13]以均匀分布的弹簧-阻尼单元模拟刀具-刀柄结合部之间的柔性连接,采用实验和仿真相结合的方法辨识刀具-刀柄结合部的刚度和阻尼系数。Mohammad等[14]使用两个线性刚度-阻尼单元计算刀柄和刀具结合面的参数,避免了主轴-刀柄结构的旋转频响函数的计算,并使用遗传算法对线性刚度和阻尼参数进行辨识。王二化等[15]利用传递矩阵法与RCSA耦合算法预测刀尖点频响函数,采用粒子群算法辨识主轴-刀柄、刀柄-刀具结合面的参数。Park等[16]通过反向耦合的方法推导出主轴刀柄结构旋转频响函数的表达式,并结合实验通过非线性方程组计算出主轴的旋转频响函数,然后由RCSA方法预测出刀尖点频响函数。在此基础上,Mancisidor等[17]使用固定边界条件的方法来计算两端自由的刀具子结构端点响应,进而通过RCSA方法预测刀尖点频响函数。
上述研究中,刀尖点频响函数预测方法或者需要辨识主轴-刀柄、刀柄-刀具结合面的参数,或者需要自制刀柄模型来求取主轴-刀柄处的转动频响函数。而结合面参数的辨识计算量大、容易陷入局部最小值;自制刀柄模型费时、费力。针对这些问题,本文提出一种改进的基于RCSA的刀尖点频响函数预测方法,该方法改进机床-主轴刀柄-刀具系统子结构的划分,将机床系统划分为主轴-刀柄-部分刀杆、剩余刀杆和刀齿三部分,改进主轴-刀柄处转动频响函数的计算方法,直接由铣刀的锤击模态实验求出,通过理论计算出刀杆和刀齿部分的频响函数,最后将三个子结构进行刚性耦合确定刀尖点的预测频响函数。
1机床系统子结构划分和响应耦合分析
1.1机床系统子结构划分
机床-主轴-刀柄-刀具系统(以下简称机床系统)主要包括机床-主轴、刀柄、刀具等几个部分,刀尖点的频响函数反映了机床加工系统在刀尖点处的动态特性,它是找出机床稳定切削区域、避免切削颤振的基础。为了快速、方便地预测出刀尖点频响函数,本文对现有机床系统子结构划分方法进行改进,将其划分为主轴-刀柄-部分刀杆、剩余刀杆和刀齿三个子结构。图1以立式加工中心为例,将其划分为子结构Ⅰ、子结构Ⅱ和子结构Ⅲ,各子结构划分及坐标点的定义如图1所示。
图1中,子结构Ⅰ为铣刀的刀齿部分,其两端坐标分别记为1、2a。由于刀齿包含螺旋槽等复杂结构,为便于频响函数的理论计算,可采用等质量法将其等效为一个圆柱梁,等效直径de由下式确定:
式中,M为铣刀的质量;ρ为刀具材料的密度;ds和ls分别为刀杆的直径和长度;lf为刀齿长度。
子结构Ⅱ为去除插入刀柄部分后剩余的刀杆,为一个等直径圆柱梁,其两端坐标分别记为2b、3a,它和等效后的 刀齿可使 用Euler梁或者Timoshenko梁模型计算其频响函数。子结构Ⅲ包括机床-主轴、刀柄以 及插入刀 柄中的部 分刀杆,其端点坐标记为3b。由于坐标点3b非常靠近刀柄端面,相对于整个机床 -主轴 -刀柄结构,插入刀柄中的刀杆部分对其动态特性的影响较小,在更换不同刀具时,可忽略刀具直径的变化(较小范围内)以及刀具材料的变化对子结构Ⅲ动态特性的影响,即认为子结构Ⅲ的动态特性保持不变。另外,由于子结构Ⅲ中3b点处的动态特性已包含了主轴 -刀柄以及刀柄 -刀杆两处结合面的动态特性,因此,本文提出的子结构划分方法无需单独考虑这些结合面的动态特性,避免了这两处结合面参数的辨识。
1.2机床系统的子结构响应耦合分析
子结构响应耦合分析中,一般先通过实验或计算得到各子结构的频响函数,然后根据各子结构之间的兼容条件和平衡条件将它们耦合,最后通过计算确定整个机床系统中刀尖点处的频响函数。图1中划分的子结构Ⅰ、Ⅱ、Ⅲ相邻之间均为刚性耦合,考虑平动和转动两个自由度,施加在各子结构坐标点i(i=1,2a,2b,3a,3b)处的载荷表示为qi = (fi,mi)T,其中,fi为力向量,mi为力矩向量。各子结构 坐标点i的位移响 应记为ui=(xi,θi)T,其中,xi为平动位移响应,θi为转动位移响应。
图1中,对于子结构Ⅰ,其两端1、2a处的频响函数矩阵可记为Rij(i,j = 1,2a),ui = Rijqj,Rij 可统一写成如下形式:
其中,hij、lij、nij、pij分别为在子结构Ⅰ的点j处施加单位力引起的点i处平动位移响应、点j处施加单位力矩引起的点i处平动位移响应、点j处施加单位力引起的点i处转角响应以及点j处施加单位力矩引起的点i处转角响应。i=j时Rij为原点频响函数,i≠j时Rij为跨点频响函数。对子结构Ⅱ,同样可得到其频响函数矩阵Rmn(m,n =2b,3a)。子结构Ⅰ 与Ⅱ刚性 耦合组成 装配体Ⅰ?Ⅱ,如图2所示。
装配体Ⅰ-Ⅱ两端的 原点频响 函数记为RS11、RS3a3a,跨点频响函数记为RS13a、RS3a1。为了求取RS11和RS3a3a,在装配体Ⅰ-Ⅱ的点1处施加载荷Q1,如图2所示,各子结构端点处受到的力分别为q1、q2a、q2b,则各端点处的位移响应可写成如下形式:
根据力平衡条件
以及相容条件
可得
则RS11可表示为
即
同理可求得RS3a1为
即
同理,假设在装配体Ⅰ-Ⅱ的点3a处施加Q3a的载荷,使用与上述相同的分析方法可求出RS3a3a、RS13a:
子结构Ⅰ和Ⅱ 耦合成装配体Ⅰ?Ⅱ后,再与子结构Ⅲ在3a、3b点处进行刚性耦合,组成图3所示的整个机床系统。
根据上述子结构Ⅰ、Ⅱ的刚性耦合 分析原理,同理可得机床系统在刀尖点1处的原点频响函数Gt11为
式(14)中,仅有R3b3b未知,其余均可由梁模型和耦合原理计算得到。R3b3b为子结构Ⅲ在点3b处的频响函数,它包含了机床 -主轴本身的动态特性以及主轴与刀柄、刀柄与刀具两处结合面的动态特性。由上述分析,令
式(14)可表示为
在实际的机床切削稳定性研究中,只需考虑刀尖点的原点位移频响函数ht11,根据频响函数矩阵的对称性有l33=n33,由式(16)得
式(17)中,只有h33、n33、p33 未知,其为R3b3b频响函数矩阵的元素,求出主轴 -刀柄频响函数矩阵R3b3b 后即可通过式(17)预测出不同悬伸长度下刀具刀尖点的位移频响函数ht11。
2子结构Ⅰ、Ⅱ频响函数计算
由上述分析,子结构Ⅰ和子结构Ⅱ均为两端自由的匀质等直径圆柱梁,其两自由端点的原点和跨点频响函数矩阵可以通过Euler或者Timoshenko梁模型计算[18,19,20]。Timoshenko梁考虑了梁的剪切效应,在高频处计算结果比Euler模型更加准确,但其计算过程复杂。在实际切削加工关心的频率范围内,二者计算结果基本没有区别。为计算简便,本文使用Euler梁模型对子结构Ⅰ和Ⅱ的频响函数进行计算。对于两端自由的均质圆柱梁,设其端点分别为j、k,其平动、转动频响函数可通过下式计算:
其中,E为梁材料的弹性模量,I为梁的二阶惯性矩,i为虚部符号,d为梁的直径,η为梁材料的阻尼因子,λ为频率参数,由下式计算:
其中,m为梁的质量,ρ为梁的材料密度,l为梁的长度。式(18)~式(25)中的F1、F3、F5、F6、F7、F8、F10均为λ的函数,它们分别由下式计算:
计算出F1、F3、F5、F6、F7、F8、F10后就可根据式(18)~式(25)计算出子结构Ⅰ、Ⅱ的各自由端点处的频响函数矩阵。
3 子结构 Ⅲ 频响函数的计算
由上文分析,子结构Ⅲ在点3b处的频响函数矩阵R3b3b未知,其平动位移频响函数h3b3b可通过模态锤击实验得到,而其转动频响函数n3b3b、p3b3b难以直接测量,可通过间接测量的方法得到。本文对现有子结构响应耦合分析中主轴 -刀柄结构转动频响函数的求取方法进行了改进,首先根据整体机床系统点3处的频响函数Gt33由反向RCSA推导出R3b3b的表达式,然后通过铣刀的模态实验和有限差分法求出R3b3b中转动频响函数,避免了结合面参数的辨识和实验刀柄模型的制作。
首先求取Gt33,使用与1.2节相似的方法推导出Gt33的表达式:
由式(34)可反向求得R3b3b表达式:
其中,RS3a3a为装配体Ⅰ-Ⅱ在点3a处的频响函数,可通过式(12)计算得到。Gt33可表示为如下形式:
其中,Ht33可通过模态锤击实验获得。转动频响函数Nt33可使用有限差分法计算[21]。有限差分法求取转动频响函数方法中,最简单的为一阶有限差分,该方法只需要考虑两个坐标点处的直接位移频响函数,计算简单,使用的数据量也最少。首先将一把立铣刀安装在机床上,确定子结构划分点3的坐标,然后在距离点3右侧s处位置选择一个坐标点,记为点4,如图3所示。对铣刀刀杆上的点3、4分别进行模态锤击实验,得到点3的原点平动位移频响函数Ht33和点3到点4的跨点平动位移频响函数Ht34,然后根据下式计算出转动频响函数Nt33:
利用已得到的频响函数Ht33和Nt33,通过下式的合成得出Pt33:
由式(37)可知,Nt33由测量得 到的Ht33和Ht34 之间的差值计算得到,因此,在实验中要尽可能选择较大的距离s来保证测量到的Ht33和Ht34之间的幅值存在明显的差异,提高计算的准确性。实际模态实验中,点3、点4均在铣刀的刀杆上,距离s可根据铣刀刀杆的实际长度进行选取。同时在模态实验时对点3、点4进行多次测量取平均值作为Ht33和Ht34的结果,以降低测量中随机噪声对测量结果的影响。
至此,Gt33频响函数中的4个元素均已求出,通过式(35)就可以确定频响函数矩阵R3b3b,最终通过式(15)、式(17)即可预测出任意悬伸长度下刀具的刀尖点位移频响函数ht11,具体步骤如图4所示。
4 实验研究
为验证本文方法的有效性,在VMC850E加工中心上进行实验研究,本文研究中使用的刀柄型号为BT-40 -ER32弹簧夹头刀柄。实验中,每把刀具都使用较大力度进行装夹,可近似认为实验刀柄与刀具之间的装夹力恒定。
4.1求取子结构Ⅲ点3b处频响函数
首先通过一把实验刀具来求取子结构Ⅲ点3b处频响函数R3b3b,将此刀具记为T0,其几何参数如图5所示。T0的材料为高速钢,其弹性模量为210GPa,密度为7800kg/m3,材料阻尼 为0.0015,泊松比为0.3。子结构划分中的各相关坐标点如图5所示。
由式(1)计算出T0刀齿部分的等效直径de=8.6mm,根据式(18)~ 式(25)、式(9)~ 式(13)计算出装配体Ⅰ?Ⅱ两自由端的频响函数矩阵RS11、RS3a3a、RS13a和RS3a1。其中RS11频响函数矩阵中的位移/力频响函数的实部Re、虚部Im如图6所示。
在刀具T0的点3、4处分别进行模态锤击实验,使用Kistler9724A激振力锤对测量点施加瞬时激振力,并使用B&K4525B三向加速度传感器记录此时振动加速度响应。激振力信号和加速度信号使用NI9234数据采集卡进行同步采集。模态锤击实验装置照片如图7所示。
进行多次测量,使用ModalView模态分析软件得到频响函数Ht33和Ht34,然后通过式(18)、式(19)计算出坐标点3处的Nt33和Pt33,三者幅值H如图8所示。
最终由式(35)计算出频 响函数矩 阵R3b3b,其中位移/力频响函 数H3b3b的实部、虚部 如图9所示。
4.2刀尖点频响函数预测与分析
根据本文预测方法对三把不同的刀具进行刀尖点频响函 数的预测,预测实验 刀具分别 记为T1、T2和T3,其参数如表1所示。
在实际的刀具更换时,装夹的刀具长度可能会有所不同,但是其长度变化在较小的范围内,且刀具装夹长度的变化对整个子结构Ⅲ(包括整个机床-主轴结构和刀 柄结构)动态特性 的影响很小,故可忽略刀具装夹长度的变化对预测结果的影响。本文实验验证中各刀具装夹长度与实验刀具T0保持一致,均为20mm。
首先预测较长悬伸长度的刀具T1的刀尖点频响函数,为验证预测的准确性,在T1刀尖点处进行模态锤击实验,获得其实测刀尖点频响函数,二者实部、虚部的对比如图10所示。
从图10中可以看出,预测与实测频响函数的实部、虚部曲线在整体上符合较好,前三阶固有频率处具有很好的一致性。
然后预测较短悬伸长度的刀具T2的刀尖点频响函数,同时为了验证不同刀具直径对预测结果的影响,T2刀具的直径选为8mm。使用上述相同方法得出T2刀尖点的预测、实测频响函数,二者实部、虚部的对比如图11所示。
从图11可以看出,T2刀尖点的预测频响函数与实测频响函数的实部、虚部曲线在整体上符合较好,虽在第一阶固有频率处预测曲线具有一定的误差,但其总体趋势符合较好,且二者曲线中固有频率基本相同。第二阶、第三阶固有频率处预测曲线与实测曲线符合程度较好。说明刀具直径在较小范 围内变化 时,本文预测 方法仍是 准确的。
上述实验刀具T0和预测刀具T1、T2材料均为高速钢,为研究不同材料对预测结果的影响,对硬质合金刀具T3的刀尖点 频响函数 进行预测。使用与上述相同的方法得到T3刀尖点预测频响函数和实测频响函数,二者的实部、虚部对比如图12所示。
从图12可以看出,T3刀尖点的预测频响函数与实测频响函数曲线符合也很好,前三阶固有频率基本相同。由于刀具材料通常为硬质合金和高速钢两种,所以,本文方法适用于不同材料刀具的刀尖点频响函数的预测。
为进一步比较上述预测结果的准确性,从图10~图12的T1、T2、T3刀具的预测、实测刀尖点频响函数曲线中分别求出各自的前3阶固有频率fn,如表2所示。
从表2中可以看出,刀尖点预测和实测固有频率均比较接近,最大误差为6.9%,最小误差仅为0.8%,可认为本文预测的频响函数是准确的。
从上述的实验分析结果可知,通过本文方法预测得到的刀尖点频响函数与实测刀尖点频响函数具有较好的一致性,实测与预测的刀尖点频响函数的前三阶固有频率误差较小,同时,相比于其他刀尖点频响函数预测方法,本文预测结果准确性较高,且无需考虑主轴-刀柄和刀柄-刀具两处结合面的建模及其参数的识别,计算过程更加简单方便。本文刀尖点频响函数预测方法是可行的、准确的,可以应用于实际的切削加工稳定性研究中。
5结论
(1)提出一种改进的基于RCSA方法的铣刀刀尖点频响函数的预测方法。该方法首先将机床系统划分为机床-主轴-刀柄-部分刀杆、剩余刀杆和刀齿三个子结构,通过铣刀的模态锤击实验由反向RSCA和一阶有限差分法得出机床-主轴-刀柄结构的转动频响函数,然后由Euler梁模型计算出两端自由的刀杆和刀齿子结构的频响函数,最后将三者进行耦合预测出刀尖点的频响函数。
(2)以加工中心上三把不同铣刀为研究对象,分别得出其预测和实测的刀尖点的频响函数并进行对比,对比结果表明:刀尖点的预测频响函数与实测频响函数符合程度较高,其预测、实测固有频率之间的误差在6.9% 以内,具有较高的预测精度,本文预测方法是有效的、准确的,适用于具有相似刀柄安装情况的任意悬伸长度的铣刀刀尖点频响函数的预测。
(3)本文将机床-主轴、刀柄和一部分刀杆作为一个整体,直接通过铣刀的模态实验求取机床主轴-刀柄处的转动频响函数,避免了主轴-刀柄和刀柄-刀具两处结合面参数的辨识和实验刀柄模型的制作,计算简单,预测快速、准确,便于实际应用。
响应耦合子结构方法 篇2
车-桥耦合系统动力学建模与响应分析
车辆行驶过桥时,车与桥之间存在相互耦合作用.本文根据车的三种简化模型,分别建立了车与桥结构相互耦合作用的动力学模型,并给出桥结构在两轴车载荷作用下的`动响应计算方法,该方法很容易推广到更一般的多轴车载荷作用的情况.文中通过数值算例计算了基于车的三种简化模型桥梁的响应,将三种计算结果与实验数据比较,证明二自由度的车简化模型为最优,此时桥梁的弯矩响应与实验结果能较好地吻合.
作 者：秦远田 陈国平余岭 张方 Qin Yuantian Chen Guoping Yu Ling Zhang Fang 作者单位：秦远田,陈国平,张方,Qin Yuantian,Chen Guoping,Zhang Fang(南京航空航天大学,210016,南京)余岭,Yu Ling(长江科学院爆破与振动研究所,430010,武汉)
刊 名：应用力学学报 ISTIC PKU英文刊名：CHINESE JOURNAL OF APPLIED MECHANICS 年，卷(期)： 25(1) 分类号：O322 关键词：车辆 桥梁 动力学模型 响应响应耦合子结构方法 篇3
碰撞冲击过程是十分复杂的现象,严格的理论分析涉及撞击体和靶板的几何形状及尺寸、它们的弹塑性性质、冲击波的传播和流动、有限变形和应变、应变率、加工硬化、热性质、摩擦效应,以及断裂的发生和扩展等。因此,长期以来人们主要是依靠大量的试验研究来处理这些问题[1]。
飞机鸟撞在航空工业中是常见的事故,鸟撞过程是一个极其复杂的过程,鸟体会发生剧烈的变形和飞溅,表现出液体的流动性。对其完善的描述应包括连续介质力学的各个方面,而且还要考虑到溅射现象,因此材料参数和模型的选取是鸟撞分析非常重要的环节。在鸟撞问题材料参数和模型选取上,主要考虑以下三方面的内容:
1)蒙皮材料弹塑性本构关系模型,即大应变率下材料的应力-应变关系描述;
2)物态方程,即高压下鸟体的流体压力和密度、内能等的关系描述;
3)材料的失效准则,对于有限元计算主要是考虑单元的失效标准(如蒙皮被穿透)。
本文将利用大型专业分析冲击碰撞有限元程序PAM-CRASH进行弹体冲击靶板的仿真模拟研究。通过SPH单元为鸟体建模,由于SPH方法和有限元求解原理的差异,利用SPH方法求解时不需要定义材料的失效准则,材料是否失效,完全根据SPH节点的求解结果来确定。因此,和有限元方法相比,SPH方法减少了对试验数据的依赖性。另外SPH方法可以克服Euler方法难于跟踪物质变形和不能识别材料界面位形的缺点;解决了传统Lagrange方法在大变形下的网格扭曲(Distortion)问题;并克服ALE方法无法模拟鸟体穿透蒙皮的问题。SPH方法在分析两种强度相差很大的物体相撞上,有独特的优势。
1 算法简介
1.1 显式时间积分算法
PAM-CRASH采用显式时间积分求解瞬态动力学问题,既要把分布空间域进行离散,把连续的微分方程转换成有限阶代数方程,还要在时间域进行离散。算法如下[2]:
如果当前系统状态已知的时刻tn,则各物理参数的近似值应满足以下方程:
Man+Cvn+Kdn=Fundefined (1)
(1)式中,M:结构质量矩阵,C:结构阻尼矩阵,K:结构刚度矩阵,Fundefined:外加载荷列阵,an:加速度在tn时刻的近似值,vn:速度在tn时刻的近似值,dn:位移在tn时刻的近似值。
将(1)式改写为:
Man=Fundefined-Cvn-Kdn (2)
an=M-1Fundefined (3)
其中,Fundefined=Fundefined-Cvn-Kdn,M-1为M的逆矩阵。加速度可以通过对质量矩阵求逆乘以剩余力矢量求得。
在时间推进上采用中心差分法:
undefined
dn+1=dn+vn+1/2Δtn+1/2 (5)
即假设加速度在一个时间步长内是恒定的。
1.2 SPH算法
SPH法全称为“光滑粒子流体动力学”方法,是Lucy于1977年提出的用于天体物理计算的一种无网格化的Lagrange计算方法,其最初模拟天体物理现象,如两星相撞、超新星爆炸以及月球的形成等等。近年来被广泛应用于处理结构动力学问题,主要是用于结构大变形解体、碎裂等分析(如高速碰撞、流固耦合碰撞等)[3]。其特点是在模拟物体大变形时,既可以克服Euler方法难于跟踪结构变形和不能识别材料界面位形的缺点,同时也解决了传统Lagrange方法在大变形下的网格扭曲(Distortion)问题,因此在工程上有很大的应用潜力。
SPH方法的核心是一种插值技术[4,5]。在SPH方法中,任意宏观变量(如密度、压力、温度、内能等)能方便地借助于一组无序点(disordered points)上的值表示成积分插值计算得到。
理论上,任意粒子的连续函数的值或其导数可以利用周围粒子的已知值通过未知函数 δ(x-x′)精确表示:
f(x)=∫Ωf(x′)δ(x-x′)dx′ (6)
(6)式中:f是三维空间位置矢量x的函数,dx′表示体积微分, δ(x-x′)为未知函数。
δ(x-x′)函数在数值计算难以实现,需要其它连续函数来近似,如果将δ(x-x′)函数用w(x-x′)来代替,则上式就不再严格成立。
即,函数f(x) 可以近似为:
f(x)≈∫Ωf(x′)w(x-x′)dx′ (7)
(7)式中,称w(x-x′)为核函数或光滑函数,它应满足下列条件:
1)半正定性,即在紧支子域内满足w(x-x′)≥0;
2)紧支性,即在紧支子域外满足w(x-x′)=0;
3)归一性,即∫Ωf(x′)w(x-x′)dx′=1;
4)w(x-x′)是在d=‖(x-x′)‖的单调递减函数;
5)当支撑域大小h→0时,w(x-x′)→δ(x-x′))。
2 仿真分析模型
2.1 有限元模型[6]
PAM-CRASH前处理工具中,为机翼前缘及鸟体建模,其中,鸟体采用SPH单元建模,铆钉采用专用的铆钉单元建模,共389个;翼肋、蒙皮采用shell单元建模。由于蒙皮在撞击过程中,变形剧烈,故蒙皮的单元尺寸较小,采用2.5 mm×2.5 mm的shell单元,共97 600个。翼肋采用10 mm×10 mm的shell单元,共7 191个。
由于鸟体在鸟撞过程中会发生剧烈的变形、碎裂,故采用12 800个SPH单元进行建模。PAM-CRASH可以通过SPH单元转换器将普通的Lagrange体单元装换成SPH单元。真实鸟的外形比较复杂,本文采用两端半球的圆柱体作为鸟体几何形状,鸟的重量为1.8 kg,直径为113 mm。如图2所示。
蒙皮和翼肋之间的铆钉是由自由的PLINK单元建立模型,也就是说,建立的铆钉与蒙皮和翼肋的网格没有几何上拓扑关系,是独立于蒙皮和翼肋的单元。
2.2 材料模型
PAM-CRASH提供丰富的材料本构模型,对于蒙皮我们可以选用金属弹塑性模型(shell103,105,171等等),也可以选用复合材料模型(shell130.131.132)。这些材料模型结合最广泛地应变硬化理论。本文中,蒙皮选用EWK(ESI-Wilkins-Kamoulakos)材料断裂模型。EWK采用等效塑性应变累积损伤模型,可以表征裂纹尖端附近塑性流状态,并能表征材料的塑性不稳定性。模型的基本定义形式为[2]:
DP=∫w1w2dundefinedP, (8)
(8)式中,undefined为流体静压加权项;参数a表征裂纹尖端附近塑性流状态:undefined;w2=(2-A)为非对称应变加权项,undefined;参量α和β表征材料的塑性不稳定性。
裂纹出现的条件为:DP>DC,距离范围r>rC,DC与rC为破坏门槛值,均由试验获得。
本文中蒙皮采用T2024硬铝,其EWK材料参数为:密度2 800 kg/m3,弹性模量72.4 GPa,泊松比0.33,应力应变硬化规律由曲线给出。
对于SPH鸟体我们选用SPH材料卡(solid28),材料Murnaghan状态方程的压力由下式计算[2]:
P=P0+B((ρ/ρ0)γ-1) (9)
(9)式中,ρ/ρ0现在的密度与原始密度的比值。γ为放大系数,B是体积模量。本文中的鸟体的模型材料参数为B=128 MPa,γ=7.98。
本文采用钛合金铆钉,对于铆钉,需要定义铆钉的断裂准则,钛合金铆钉的轴向极限拉力为2 814 N,极限剪切力为15 900 N。铆钉的失效准则可以定义为[2]:
undefined
即:当左边表达式大于等于1时,铆钉失效,当小于1时,铆钉正常工作。
2.3 接触模型
在经典的有限元方法中,结构之间的相互作用要有网格的连贯性作为保障的。这就使经典有限元方法能很好的处理结构小变形问题,但是很难处理大变形问题,以及网格不连贯,结构之间相离的问题。
PAM-CRASH提供多种接触模型来模拟不同形式的接触(比如,点面接触、面面接触、自接触等等)。接触检查分成两步,第一步进行全局检查,第二步进行接近检测。全局检查用于检查潜在接触对,接近检测检查潜在接触对的接触,并计算接触力及接触力方向。PAM-CRASH接触力的计算使用罚函数法,有三种惩罚方式:线性惩罚、非线性惩罚、运动叠加惩罚。
在鸟撞模拟中,接触模型选用点面主从接触。鸟是从接触,蒙皮是主接触面。
3 仿真分析结果
使用上述模型,分析1.8 kg鸟体以150 m/s的速度撞击蒙皮厚为2 mm的机翼前缘。分析显示鸟体在撞击过程中剧烈变形,碎裂飞散。鸟体在撞击过程中的表现如图4。
在鸟体的猛烈撞击下,蒙皮发生剧烈的塑性变形,并开裂。蒙皮裂纹扩展过程及等效应力云图如图5。
由于鸟体的剧烈剪切作用,蒙皮在和翼肋连接处最先发生剪切失效,并开裂,如图8(a)所示。肋耳片在鸟体的剪切作用下,在和肋腹板连接位置也发生剪切失效,如图8(b)所示。肋腹板在鸟体的剧烈挤压作用下,发生受压屈曲失效,如图8(c)和8(d)所示。
鸟撞过程鸟体和蒙皮的接触力(撞击力)时间历程曲线如图6所示。
在撞击点附近,铆钉会出现失效,随着撞击过程的发展,撞击面的不断扩大,铆钉会不断的出现失效。在3.3×10-3s时,铆钉失效图如图7。
4 结论
(1) 采用SPH方法,能够准确地分析鸟体的剧烈的大变形和飞散现象,避免了传统FEM单元鸟体在撞击过程中因单元失效引起的质量损失,更为准确的模拟了鸟撞的冲击载荷;解决了传统Lagrange方法在大变形下的网格扭曲问题,使计算更加稳定。
(2) 应用应力-应变曲线硬化规律定义的材料参数可由试验直接获得,避免了材料数据在拟合材料本构的材料参数时引起的数据丢失;文中采用EWK材料断裂模型准确的分析了蒙皮在冲击载荷下的断裂失效及裂纹扩展情况 。
(3) 铆钉是机翼结构的主要连接元件,文中采用独立的铆钉连接单元解决了机翼鸟撞模拟建模的难点,更好地模拟铆钉的连接作用及失效机理,使模型更加接近真实的机翼结构。
摘要：鸟体的强度远比铝合金小,鸟体在撞击过程中表现强烈的流体流动性。介绍利用PAM-CRASH软件,采用SPH(Smooth Particle Hydrodynamics)和FEM耦合技术分析鸟撞机翼前缘问题。通过将鸟体离散成光滑粒子(SPH)单元,并且采用铆钉连接单元简化蒙皮和肋板的铆钉连接,准确分析了鸟撞动态特性。SPH方法在鸟撞仿真中,克服了网格稳定性问题,并能有效分析鸟体的飞散现象。
关键词：PAM-CRASH,鸟撞,机翼前缘,有限元,SPH,铆钉
参考文献
[1]钱伟长.穿甲力学.北京:国防工业出版社,1984:22—25
[2]PAM-CRASH2G solver note manual,ESI group,2006
[3]张雄,刘岩.无网格法.北京:清华大学出版社,2004
[4]闫晓军,张玉珠,聂景旭.空间碎片超高速碰撞数值模拟的SPH方法.北京航空航天大学报,2005;31(3):351—354
[5]Haug E,Groenenboom P,Kamoulakos A,et al.Application of SPH techniques in the PAM-SCL code family.In:PUCA’97.ESI Group,1997
响应耦合子结构方法 篇4
结构在外界激励下会产生振动, 不必要的振动会降低结构的性能和使用寿命, 通过结构形状优化能有效减小结构振动。结构形状优化设计建立在最优理论和有限元分析的基础上, 它是根据给定的设计约束条件, 求解满足约束要求且某种结构性能达到最优的结构形状。结构形状优化通过改变结构的形状来改善结构性能和减小结构质量[1], 在工程应用领域中占有非常重要的地位。
近来来, 国内外学者已在结构形状优化方面作了大量研究。文献[2]选择二维结构边界的一组关键节点坐标构建B样条函数, 结构边界上的其他节点和内部节点坐标通过B样条函数插值获取, 以应力为设计约束, 以质量为设计目标, 结合二次规划优化方法和有限元算法对平面结构进行形状优化。文献[3]针对薄壁梁结构有限元模型, 以梁截面关键节点为坐标建立了梁截面形状与截面惯性矩的关系, 优化了梁截面惯性矩参数。文献[4]利用NURBS曲线定义结构形状, 提出了边界曲线的变形能模型, 以结构应力为设计约束, 通过优化降低了结构的质量。文献[5]提出基于结构边界关键节点的形状参数化方法, 结构内部节点通过等参形函数插值获取, 能有效优化二维平面结构应力。文献[6]提出了一种等几何结构形状优化方法, 通过非均匀有理B样条曲线控制结构形状, 将控制点坐标作为设计变量, 实现结构边界光滑。文献[7]通过对二维结构边界逐步寻优, 结合模拟退火方法优化二维结构的应力, 得到了光滑的边界和全局最优解。
在现有的研究工作中, 结构形状优化主要集中在结构静力分析领域, 以结构动力响应为优化目标和约束的结构形状优化并不多见。结构形状描述方法则广泛基于B样条曲线或曲面, 用多项式插值函数描述结构边界, 以控制点坐标为设计变量, 这对复杂结构来说, 会导致设计变量过大或难以对几何形状参数化。针对上述问题, 本文建立了基于有限元网格的参数化形状优化模型, 提出一种基于振动响应场的结构形状优化方法以控制结构振动。基于振动响应场建立结构形状参数化模型, 可减少网格变形失真, 形状变化不需要指定扰动方向, 能实现结构形状自动寻优。研究了外部激励下的结构振动优化问题, 以分析频率带内的最大振动位移为设计目标, 将有限元结构坐标的最大变动值及一阶模态频率作为设计约束, 利用遗传算法进行优化求解。数值算例表明该方法能有效降低结构振动位移幅值, 改善结构动力学特性。
1 结构振动响应的有限元计算
结构有限元的单元质量矩阵和刚度矩阵分别为
式中, ρ为材料密度;N为单元形函数;D为材料参数本构矩阵;B为应变矩阵。
dV和应变矩阵B与结构几何形状相关。
以等厚的线性四边形壳单元为例, 通过坐标映射, 等参单元的质量矩阵和刚度矩阵可简化为
其中, t0为单元的厚度, 雅可比矩阵的行列式|J|与结构几何形状相关, 单元质量矩阵和刚度矩阵采用高斯积分计算。
单元质量矩阵和刚度矩阵装配后, 在外部稳态激励力作用下结构的振动方程为
式中, M、C、K分别为结构的质量矩阵、阻尼矩阵和刚度矩阵;u为位移响应向量;f (t) 为外力向量。
采用比例阻尼假定, 结构阻尼矩阵可表示为
式中, α、β为比例系数。
对式 (5) 两边作傅里叶变换, 得
式中, ω为圆频率;U为位移响应向量u的傅里叶变换;F0为外力向量f (t) 的傅里叶变换。
从式 (6) 可求得结构在外力作用下的频域响应。修改结构的几何形状可改变单元质量矩阵和刚度矩阵, 从而改变整个结构的质量矩阵和刚度矩阵, 进而改变结构在外力作用下的频域响应, 因此可通过优化结构形状来改善结构动力学特性。
2 基于振动响应场的形状参数化模型
结构有限元模型形状由各节点坐标决定, 为了便于描述有限元的形状参数化, 用一个N维列向量X定义节点的位置, 向量包含每个节点的x、y、z方向坐标值, 则第i个节点x、y、z方向的坐标值分别对应列向量X第3i-2、3i-1、3i元素的值, 用一个n维列向量b定义设计变量, 则向量X描述了结构的形状, 结构形状可通过一组向量场的线性组合定义, 结构形状参数化模型可表示为
式中, X0为结构初始形状;bi为设计变量;qi为向量场。
若直接以结构修改区域内的所有节点坐标作为设计变量, 虽然可以建立形状参数化模型, 但由于参数设计变量庞大而无法应用于工程问题, 而且会导致有限元网格严重变形而使计算结果失真, 故必须寻求一种形状映射函数, 以尽可能少的设计变量来描述各节点坐标, 建立有效的结构形状参数化模型。
结构振动响应场是在外部激励下的振动分布, 根据振动响应的类型可分为位移场、速度场和加速度场。考虑结构在外部宽频动载荷激励下产生频域响应U0, U0是频率的函数, 每一频率对应一个振动响应场, 振动位移响应场向量包含了每个节点的x、y、z方向的振动响应位移幅值, 因此振动位移响应场与结构形状参数化模型中的向量场具有相同维数, 可用振动位移响应场作为结构形状参数化模型的向量场。我们采用一组振动位移响应场作为结构形状参数化模型的一组向量场, 以振动响应场的系数作为优化设计变量, 不同的设计变量可以映射出不同的结构形状。通过对初始结构进行动力学分析, 得到结构频域响应U0, 在分析频带内选择若干个与振动峰值对应的位移振动响应场作为形状参数化模型的向量场, 不同的激励特性可得到不同的向量场。由于位移振动响应场是连续变化的, 故采用位移振动响应场可减少形状优化过程中的网格失真变形。
为了便于描述形状改变后几何坐标系下节点新的坐标值, 将振动响应场进行归一化, 归一化后的振动响应场最大振幅为1。给定一组振动响应场的系数, 则对应一种唯一结构形状。每一结构形状下的各节点通过三个方向的坐标描述。以x0和x1表示初始形状和新形状下所有节点在x方向的坐标值、y0和y1表示初始形状和新形状下所有节点在y方向的坐标值、z0和z1表示初始形状和新形状下所有节点在z方向的坐标值, 则结构形状参数化模型的节点坐标可以表示为
式中, bi (i=1, 2, …, n) 为选取的第i个振动响应场的系数, 即将此系数作为设计变量, 取值范围与结构节点坐标波动限值相关;q (x) i、q (y) i、q (z) i分别为选取的第i个振动响应场x、y、z方向的振动位移幅值。
结构节点坐标在x、y、z方向的波动值可以表示为
3 结构形状优化
在所分析的频带内, 振动响应是频率的函数。以单一固定频率处响应为优化目标无法控制整个频带内的振动响应, 甚至可能导致频带内的最大振动峰值增大, 而以整个分析频带内的结构振动最大响应值为优化目标则可控制整个频带内的振动响应。为了能获取整个分析频率带内的所有振动峰值, 取分析频率步长为1Hz, 计算振动响应, 识别出分析频带内的最大响应值, 则目标函数可表示为
式中, fmin、fmax分别为分析频带的下限频率和上限频率;resp为振动响应 (如振动位移、速度或加速度响应) ;F为目标函数, 是设计变量b1、b2、b3、…、bn的函数。
以最小化振动响应为目标, 结构动力学形状优化模型可表示为
式中, biL、biU分别为设计变量的下限和上限。
遗传算法能在较大的设计变量空间内迅速寻优, 有较强的全局优化性能[8], 本文采用遗传算法对结构形状进行全局优化。取设计变量bi为遗传算法的群体, 目标函数F为遗传算法的适应度函数, 选择遗传算法的交叉概率为0.5, 变异概率为0.01。
本文方法的流程如图1所示。
4 数值算例
四边简支矩形平板, 长宽分别为0.8m和0.4m, 板厚为1mm。结构材料性能参数为:弹性模量210GPa, 泊松比0.3, 密度7800kg/m3。由于在形状优化过程中, 结构会变化成空间壳结构, 故结构采用四节点四边形壳单元进行计算, 共有200个单元, 231个节点。板在节点 (0.48, 0.28, 0) (m) 位置处作用一个1N的竖向简谐激振力, 计算结构激励点的位移响应函数, 分析频率为15～300Hz, 步长为1Hz, 平板的初始形状有限元模型如图2所示。
由于平板横向振动主要由板的面外弯曲振动引起, 结构节点在面内x坐标和y坐标的变化对横向振动影响很小, 故不考虑。则节点坐标x、y在迭代过程中保持不变, 只考虑由节点z坐标引起的结构形状变化, 结构形状参数化模型的节点坐标可以简化为
本文以分析频带内加载点最大振动位移响应值作为目标函数, 为了提高结构的动态性能, 优化后结构一阶模态频率应高于结构初始形状一阶模态频率值, 为了保证结构网格单元的质量和结构的可加工性, 节点坐标波动值不宜过大, 取最大波动值Δzmax=10mm, 则设计变量范围为0～10即可满足节点波动值要求, 该算例的具体优化模型如下:
先用有限元计算初始结构的位移响应函数, 激励点的位移频率响应函数如图3所示, 位移频率响应曲线的前5个峰值对应的频率分别为19Hz、30Hz、49Hz、65Hz、76Hz, 选取这5个频率处的结构位移振动响应场作为形状参数化模型的一组向量场。
采用遗传算法求解上述优化问题, 求得设计变量优化解为b*1=2.88, b*2=1.57, b*3=1.38, b*4=0.28, b*5=7.90。采用本文方法优化后的结构有限元模型如图4所示, 为了看清结构形状的变化, 图4是z坐标放大了3倍的效果图。结构形状优化前后各参数结果对比如表1所示, 优化前后整个分析频率范围内的结构振动位移响应如图5所示。优化结果表明结构形状优化后的结构动态位移响应函数曲线更为平坦, 结构一阶模态频率提高了近110Hz, 最大位移响应幅值由3.914mm降到0.009mm, 频率带峰值大大减小, 优化效果明显。
5 结论
结构形状优化设计是一种有效改善结构动力性能的方法, 本文将形状参数化模型、有限元分析和遗传算法相结合, 提出一种基于振动响应场的结构形状优化方法。与B样条曲线、Bezier曲面和多项式插值参数化形状优化方法相比, 本文方法具有如下特点:能减少有限元网格变形失真;可用更少的设计变量映射结构形状;形状的变化不通过几何控制点移动实现, 而是通过向量场的系数变化实现;不需要指定扰动方向;优化设计变量与结构位移响应场的选取个数有关, 与结构的复杂程度无关等。数值算例表明采用一组结构振动位移响应场作为结构形状参数化模型的向量场是可行的;用本文方法进行结构形状优化设计能有效改变结构形状, 降低结构的振动位移幅值, 改善结构的动力学特性。
参考文献
[1] Mackerle J.Topology and Shape Optimization of Structures Using FEM and BEM:a Bibliography[J].Finite Elements in Analysis and Design, 2003, 39 (3) :243-253.
[2] Hyun S, Kim C.An Efficient Shape Optimization Method Based on FEM and B-spline Curves and Shaping a Torque Converter Clutch Disk[J].Finite Elements in Analysis and Design, 2004, 33 (13) :1803-1815.
[3] Vinot P, Cogan S.Shape Optimization of Thin-walled Beam-like Structures[J].Thin-walled Structures, 2001, 39 (7) :611-630.
[4]Wang Xuelin, Zhou Ji.A Physics-based Parameter-ization Method for Shape Optimization[J].Com-puter Methods in Applied Mechanics and Engineer-ing, 1999, 175 (1) :41-51.
[5]Song X, Baldwin J D.A Novel Node-based Struc-tural Shape Optimization Algorithm[J].Computers&Structures, 1999, 70 (5) :569-581.
[6]Wolfgang A W, Moritz A Z, Chirstian C.Isogeo-metric Structural Shape Optimization[J].ComputerMethods in Applied Mechanics and Engineering, 2008, 197 (33) :2976-2988.
[7]Frazil O S.Shape Optimization of 2D Structures UsingSimulated Annealing[J].Computer Methods in Ap-plied Mechanics and Engineering, 2007, 196 (35) :3279-3299.
响应耦合子结构方法 篇5
关键词：减基法,贪婪算法,对数采样法,谐响应,阻尼
0 引言
结构谐响应分析是结构分析很重要的一方面, 它可以确保一个给定的结构能经受住各种不同频率的正弦载荷, 有助于结构设计人员掌握结构的持续动力特性, 并在必要时避免或利用共振。结构谐响应分析需要求解大型多自由度系统方程, 这个循环求解的过程非常耗时。目前有许多降阶方法可用来解决这类问题, 如Guyan静力凝聚法[1,2]、动力子结构法[3]等, 但这些方法在提高效率和保持原模型的物理特性上仍存在一定的局限性, 即使是降阶后, 求解仍比较耗时。
减基法[4,5]是20世纪70年代提出的一种降阶方法。其基本思想是, 当系统由多个参数描述时, 这些参数的不同组合会使系统方程有不同的解, 而系统在新参数下的解可以由事先设计的样本参数组所对应解的线性组合来得到。近年来, 不少学者对该方法进行了研究。国外, Maday等[6]从理论上证明了减基法的一致指数收敛性质;Rovas[7]把减基法应用于不同种类偏微分方程的求解中;国内, 减基法在结构动力学中的应用研究处于起步阶段, 文献[8]引入减基法, 在波数域内快速分析了复合材料层合板的瞬态响应问题。李永红等[9]将减基法用于结构特征值分析与参数优化中, 求解了无阻尼结构的固有频率与振型。黄永辉等[10]通过减基法在实数范围内分析了无阻尼结构的谐响应问题。然而, 在目前可查资料中, 尚未见减基法用于有阻尼结构谐响应分析的报导。
本文在有阻尼的线弹性结构谐响应分析中引入减基法, 在离散化的频域中进行对数预采样[8], 求出样本点的解向量, 得到一系列的复数解向量。而无阻尼的结构谐响应分析得到的样本点的解向量是一系列的实数解向量, 后续工作都是在实数域中进行的。有阻尼的情况要考虑直接用复数解向量或是向量的模来构造减基空间。复数向量包含幅值和相位两个属性, 用模来构造减基空间, 把复数域的工作转换到实数域进行, 会丢失向量的相位属性, 因而考虑直接利用复数解向量, 通过贪婪算法[6]构造减基空间, 把系统的刚度矩阵、质量矩阵、阻尼矩阵以及载荷列向量分别投影到减基空间进行降阶, 得到减缩系统。在减基空间求解原问题的减缩解, 并把减缩解还原到原空间得到原问题的近似解。通过误差和计算耗时的对比, 验证了这种方法在有阻尼的结构谐响应分析中的可靠性和有效性。
1 有阻尼结构谐响应分析的有限元形式
谐响应分析是用来计算结构在简谐载荷激励下响应的方法。在谐响应分析中, 激励外载是已知的, 可以是力也可以是位移、速度或加速度。
在谐响应分析中, 激励载荷在本质上为正弦载荷, 在最简单的情况下, 这种载荷可通过特定频率下的幅值来定义。简谐响应与载荷以相同的频率出现, 由于系统阻尼的影响, 响应在时间上可能移位, 响应移位也称相位移位, 因为载荷峰值与响应峰值不是同时出现的。
在谐响应分析中, 一般有两种不同的数值方法可供选择:直接法和模态法。直接法根据外载荷频率求解耦合的运动方程;模态法利用结构的振型缩减和解耦运动方程, 对各个模态响应进行叠加得到一特定外载荷频率的解。
简谐载荷作用下的有阻尼结构的运动方程为
其中, F0 (t) 为简谐激励载荷, F0 (t) =F (ω) ei ω t;F为载荷幅值;ω为激励载荷的频率, 当载荷幅值不变时, F与ω无关, 可写成F0 (t) =Fei ω t的形式。K、M、C分别为结构的整体刚度矩阵、质量矩阵和阻尼矩阵, d、
d=d0 (ω) ei ω t (2)
式中, d0 (ω) 为复位移向量。
将式 (2) 代入式 (1) 并化简得
(K+i ωC-ω2M) d0 (ω) =F (3)
当考虑阻尼时, 式 (3) 代表复系数方程系统。当激励频率范围较大时, 式 (3) 的求解次数也随着增多, 反复求解工作量大。另一方面, 当系统自由度n很大时, 求式 (3) 非常耗时, 需要考虑计算成本问题。
在实际问题中, 要精确地计算结构受到的阻尼是很困难的。通常从整体上考虑阻尼的影响, 近似地估计阻尼力做功消耗的能量。一般有限元程序中, 假定整体阻尼矩阵C是整体刚度矩阵K与整体质量矩阵M的线性组合, 成为Rayleigh阻尼 (或比例阻尼) , 其表达式为
其中, Ω为结构固有频率;ξ为阻尼比;比例系数α、β可通过测试结构自由振动的衰减率经过换算得到。当阻尼比ξ确定时, α、β也可通过式 (5) 得到。
2 有阻尼结构的减基法谐响应分析
减基法用于结构谐响应分析的目的是缩减刚度矩阵、质量矩阵和阻尼矩阵的阶数, 同时保证减缩后的系统能很好地近似原系统, 以提高求解效率。计算过程分为离线和在线两阶段。
2.1 离线阶段
结构谐响应分析的主要变化参数是载荷频率ω。本文将频域离散成m个频率点, 对频率点进行预采样, 目前预采样方法有对数采样法、等间隔均匀采样法、Chebychev采样方法等。本文用对数采样方法采取N个样本点, 采样公式如下:
式中, λ为常数;μj为采样参数的值;μmax为采样参数的最大值。
N<m时, 初步构造频率的样本空间如下:
SN={ω1, ω2, …, ωN} (7)
在每一个样本点求式 (3) , 由于本文考虑的是有阻尼的系统, 式 (3) 中存在复数项, 所以得到的N个解向量也是复向量, 样本点的复数解向量构成解空间如下:
WN=span (d0 (ω1) , d0 (ω2) , …, d0 (ωN) ) (8)
将d0 (ω) 分离成与变化参数有关部分α (ω) 和无关部分D, 则有
d0 (ω) =Dα (ω) (9)
那么式 (3) 可写成
(K+i ωC-ω2M) Dα=F (10)
式中, α为原系统在减基空间的减缩解。
两边乘以DT得
DT (K+i ωC-ω2M) Dα=DTF (11)
当D为n×N1阶正交矩阵时, 称KN1、CN1、MN1、FN1分别为K、C、M、F关于D的正交投影, 且K、M、C、F的阶数由n×n变成N1×N1, 当KN1=DTKD, CN1=DTCD, MN1=DTMD, FN1=DTF, N1<n时, 矩阵得到降阶。
为了得到正交矩阵D, 构造减基空间, 本文引入贪婪算法, 对解空间WN进行自适应训练。这个过程中, 与无阻尼情况的区别是, 样本点的解向量都是复数向量, 构成的空间也是复空间, 贪婪算法是在复数域进行的。为了保证解向量的正交性, 同时避免减缩矩阵的病态问题, 每将一个解向量加入减基空间, 就要采用施密特正交化方法对减基空间的向量进行正交归一化。贪婪算法的具体步骤如下:
(1) 从预采样点的解空间中选择一个解向量加入减基空间。
(2) 当构造减基空间的向量个数大于1时, 用施密特正交化方法把解向量正交归一化, 构造新的减基空间, 并得到正交矩阵D。
(3) 通过伽辽金映射将与参数ω无关的K、M、C、F矩阵投影到减基空间, 得到缩减后的矩阵和系统控制方程, 在预采样点进行求解, 得到N个解向量。
(4) 利用2范数, 求N个样本点的减基误差和投影误差。减基误差向量er的第j个分量定义为
er (j) =‖d0j-drj‖2/‖d0j‖2 (12)
投影误差向量ep的第j个分量定义为
式中, d0j为直接求解的第j个样本点的解向量;αj为减基空间中的第j个样本点的解向量。
(5) 判断最大减基误差是否小于给定的误差限ε, 若满足, 停止贪婪算法;若不满足, 将最大投影误差对应的预采样点的解向量放入减基空间, 重复步骤 (2) ～步骤 (4) 。
通过贪婪算法构造的减基空间能很好地近似原来的空间, 可表示为
WN1=span (η1, η2, …, ηN1) (14)
式中, ηi (i=1, 2, …, N1) 为构成减基空间的样本点的解向量。
且一般情况下N1<N。
2.2 在线阶段
利用正交矩阵D将K、M、C、F降阶得减缩矩阵KN1、CN1、MN1、FN1。
在有阻尼的情况下, D= (η1, η2, …, ηN1) 为列正交复数矩阵, 降阶后KN1、CN1、MN1、FN1均为复数矩阵, 在减基空间求解m个频率点的响应:
(KN1+i ωCN1-ω2MN1) α=FN1 (15)
由于减基空间中的各个向量线性无关, 故这些向量可以作为任意参数下解空间中N1维子空间的基, 则原系统的解向量可以用减基空间解向量的线性叠加来近似:
drk=Dαkk=1, 2, …, m (16)
其中, d
3 减基法与直接法之间的误差
为了验证减基法的有效性, 需要定义减基法求解与直接求解之间的误差, 利用欧几里德范数定义相对误差:
e (k) =‖d0k-drk‖2/‖d0k‖2 (17)
其中, d0k为直接求解式 (3) 得到的解向量, 根据数值分析中方程组近似解可靠性判别法[11], 可利用系数矩阵条件数以及残余向量来判断减基法得到的近似解是否可靠。令
当
e≤Cond (A) ‖F-Adr‖2/‖F‖2 (18)
时, 近似解是可靠的。其中, Cond (*) 为求矩阵A条件数的函数。
4 算例
算例一 以一L形振动试验夹具为例, 材料常数如下:弹性模量E=70GPa, 泊松比μ=0.33, 密度ρ=2.8×103kg/m3, 夹具板厚度0.01m。该夹具通过4个螺栓连接在振动平台台面上, 其几何结构尺寸及螺栓位置如图1所示, 有限元模型如图2所示, 共有791个节点、2373个自由度, 在747节点y方向施加简谐载荷, 幅值1000N, 频率3500～4500Hz, 系统受瑞利阻尼作用, 阻尼比ξ=0.05。
将频域离散成200个子步, 用对数预采样采集31个频率点, 给定的减基误差限ε1=1×10-5, 经过贪婪算法, 最终选出7个点即可满足精度要求。减基法把原来2373×2373的系数矩阵缩减为7×7的矩阵, 大幅度提高了求解效率。计算用MATLAB程序实现, 在CPU主频为2.01GHz、内存为1GB的电脑上运行, 直接求解与减基法方法求解时间对比见表1。减基法求解所耗时间几乎是直接求解的1/5, 效率明显提高。
贪婪算法过程中, 减基误差、投影误差分别如图3、图4所示, 减基误差从第1代到第2代迅速减小, 随后逐渐收敛于误差限10-5。投影误差变化趋势与减基误差基本相同, 第3代后缓慢收敛, 继续增加基数量对误差影响不大, 以减基误差判断, 当基数量达到7时, 误差小于给定误差限。
相对误差如图5所示, 其中最大相对误差为3.73×10-11, 满足式 (18) , 图6给出了减基法求解与直接求解得到的响应幅值, 找到一频率在4000Hz附近的共振点, 两种所得结果一致, 说明减基法能得到高精度的解, 且所耗时间比直接求解所耗时间少。
算例二 考虑如图7所示的板梁支架结构, 梁底端固支, 弹性模量E=200GPa, 泊松比μ=0.3, 密度ρ=7.8×103kg/m3, 板长0.8m, 宽0.4m, 厚0.01m。梁截面为实心, 宽0.03m, 高0.04m。支架的每段支撑高0.4m。有限元模型如图7所示, 共有405个节点、2430个自由度, 分别在支架的219节点、270节点、347节点加载z方向同频率简谐载荷, 幅值分别为1000N、750N、500N, 频率为0～300Hz, 同样受瑞利阻尼作用, 阻尼比ξ=0.06。
与算例一同理, 在离散为300个子域的频域中采取16个点, 贪婪算法选出10个点即可满足精度要求, 本例减基误差限ε=1×10-3, 两种方法耗时如表2所示。减基误差、投影误差的变化如图8、图9所示, 趋势基本相同, 从第1代到第7代迅速减小, 随后逐渐收敛于误差限10-2。当基数量达到11时, 误差小于给定误差限, 此时, 构造的减基空间满足精度要求, 减缩系统的刚度矩阵、质量矩阵、阻尼矩阵阶次由2430×2430降到11×11。
最后的相对误差如图10所示, 其中最大误差为1.01×10-5, 满足式 (18) , 同时由表2知, 减基法耗时只是直接求解耗时的1/7。图11给出了减基法与直接法求解的响应幅值, 实线为有限元法直接求解的结果, 点线为减基法求解结果, 由图11可见减基法与原有限元法的结果一致, 通过曲线找到共振点频率约为80Hz。
以上两个算例说明在有阻尼的线弹性结构谐响应分析中, 利用减基法求解可以明显提高效率, 并保证解的精度, 最大误差均满足误差限, 说明减基法和原方法的结果是可靠的, 也说明了该方法的有效性。
5 结束语
本文在有限元的基础上, 把减基法用于有阻尼的线弹性结构谐响应分析。分析过程中, 把简谐载荷频率离散, 通过对数预采样和贪婪算法构造减基空间, 把原系统映射到减基空间进行降阶, 在很大程度上缩减了系统刚度矩阵、质量矩阵、阻尼矩阵的阶数, 从而节省了求解资源, 有效提高效率。本文算例表明, 减基法在解决无阻尼的结构谐响应分析的基础上, 同样适用于有阻尼的结构谐响应分析。