基于绝对节点坐标的废塑料薄膜振动特性研究

2022-04-27王婷杨先海李永波

山东理工大学学报（自然科学版） 2022年4期

王婷，杨先海，李永波

(山东理工大学机械工程学院，山东淄博 255049)

近年来，随着塑料制品的不断使用，塑料垃圾也越来越多，废塑料污染问题已成为一个非常严重的环境问题。在解决生态环境污染的问题中，整治废弃塑料所带来的环境污染是最重要的方面之一，而将废塑料垃圾从混合生活垃圾中分离出来是首先应该被实施并完成的。废塑料薄膜等柔性材料在实际分选过程中本身的振动特性会直接影响其运动轨迹，从而影响分选纯度和效率；因此,需对废塑料薄膜进行振动特性分析，进而提高废塑料薄膜分选的精度和效率。

现有的薄膜振动特性研究大多利用线性振动动力学方程来进行分析[1-2]。Tian等[3]在假定来流为均匀不可压缩理想势流的前提下，研究了空气中薄膜振动产生的附加空气质量以及气动力对膜振动频率和振幅的影响；Zhang等[4]提出了高功率脉冲激光器产生等同于声激励的理想脉冲激励对薄膜进行振动激励的方法来研究薄膜振动特性，并利用实验验证了该方法的可行性；Si等[5]采用有限元方法对变压器油膜的振动特性进行了数值研究，并考虑了粘度对能量损失的影响；邵明月等[6]、武吉梅等[7]基于Von Karman薄板理论推导出轴向运动薄膜的非线性振动方程，得出薄膜稳定工作区间和发散失稳区间；刘明君[8]以薄膜的褶皱形态为新的平衡状态进行动力学分析，建立了不同载荷下褶皱薄膜的动力学模型，并准确地描述出褶皱薄膜的振动特性。上述计算分析方法都有效地描述出不同薄膜的振动特性，但计算大多较为复杂，效率较低。

绝对节点坐标法(absolute nodal coordinate formulation，ANCF)是Shabana教授为解决柔性体的变形问题而提出的[9]，与传统有限元方法不同，该方法中节点坐标选取的是全局坐标系下的位置及其梯度，从而避开了小转角的约束[10]。一些学者基于绝对节点坐标法讨论了不同维度下的单元模型，进而求解动力学问题[11-14]。上述的研究集中于对非柔性体单元模型的建模仿真和弹性力的研究，将绝对节点坐标法应用于塑料薄膜振动特性的研究较少。

由于塑料薄膜的振动特性对分选效率的影响较大，本文在利用振动对废塑料薄膜进行分选的基础上，基于绝对节点坐标法建立薄膜单元的三维柔性模型，分析废塑料薄膜的振动特性，并通过MATLAB、ANSYS等软件验证该方法的正确性，以期为后续的振动分选提供理论基础。

1 基于ANCF的薄膜单元模型

ANCF通常用于求解梁结构的变形和位移问题，在传统算法中，由于薄膜的厚度较小，通常忽略厚度的影响，本文在薄板单元模型的基础上研究薄膜的振动特性，同时考虑x、y和z的梯度，进而增加了计算精度。三维柔性薄膜单元模型如图1所示，图中O-XYZ为全局坐标系，与单元模型相连结，P0和P分别为初始构型与当前构型时单元模型上的质点，r0和r为节点的位移向量，∂rij/∂x、∂rij/∂y和∂rij/∂z均为节点坐标。

图1 三维柔性薄膜单元模型

三维柔性薄膜单元中每个节点含有3个位置坐标和9个位移斜率坐标，单元的自由度为48，因此单元上任意一点的位置坐标可利用节点坐标和形函数来表示，即

r(x,y,z)=S(x,y,z)qe(t)，

(1)

式中：S为单元的形函数；x和y为任意的局部坐标；qe(t)为单元节点坐标。

单元的节点坐标可定义为

(2)

其中任一节点坐标可表示为

(3)

形函数可表示为

S(x,y,z)=[S1IS2I…S16I]，

(4)

各个分量如下：

1.1 单元质量矩阵

薄膜单元的质量矩阵可以利用薄膜的动能来获得，薄膜单元的动能可以表示为

(6)

根据节点速度矢量，可以求得绝对速度矢量的表达式

(7)

代入薄膜单元动能表达式中可得

(8)

从而得到薄膜单元质量矩阵

(9)

1.2 弹性力矩阵

由于薄膜为柔性材料，利用传统的线性求解方法已无法满足单元模型对于材料属性的精确描述，因此在非线性连续介质力学的基础上，利用Neo-Hookean本构模型来求解三维柔性薄膜单元的弹性力矩阵，进而求解出单元的刚度矩阵。将废塑料薄膜材料变形视为各向同性，基于Neo-Hookean本构模型的薄膜单元应变能密度函数可以表示为

(10)

式中：变形张量不变量I1可以定义为右柯西-格林变形张量C的迹，即

I1=tr(C)；

(11)

J为位置坐标求偏导的矩阵，可以表示为

J=det(J)=|J|。

(12)

将应变能密度函数对节点坐标进行微分，可以得到三维柔性薄膜单元的弹性力矩阵

(13)

根据右柯西-格林变形张量定义得

因此C的迹为

I1=tr(C)=qeTSaqe+qeTSbqe+qeTScqe。

(15)

式(15)对节点坐标求微分可得

(16)

薄膜单元的变形梯度为

(17)

将式(17)代入J的矩阵行列式可得

J=S1xqeS2yqeS3zqe+S3xqeS1yqeS2zqe+

S2xqeS3yqeS1zqe-S1xqeS3yqeS2zqe-

S2xqeS1yqeS3zqe-S3xqeS2yqeS1zqe。

(18)

由于柔性薄膜材料的泊松比与不可压缩材料的泊松比极为接近，因此在以往的计算中都将其视为不可压缩材料。但在实际操作条件下，材料的不可压缩条件是无法满足的，上述假设会影响计算结果的精确性，因此需在推导公式时增加一个补偿方程。薄膜单元的应变能密度函数可以表示为

(19)

式中

(20)

取系数k=1×109N/m2。

求得薄膜单元的弹性力矩阵为

(21)

薄膜单元的弹性力矩阵还可表示为

(22)

式中K为薄膜单元的刚度矩阵，可表示为

K=(λ+2μ)K1+λK2+4μK3，

(23)

式中

1.3 单元广义外力

当有外力F作用于薄膜单元上时，该力的虚功可表示为

δWj=FTδr=FTSδqe=QjTδqe，

(25)

则薄膜单元的广义外力Qj为

Qj=S(x,y,z)TF，

(26)

薄膜结构中所有单元的广义外力为

Q=∑n1×n2HTQj。

(27)

2 废塑料薄膜动力学方程及求解

由于薄膜单元被离散为多个单元，假设系统中包含n个单元，第i个单元的单元节点坐标为qei，利用布尔矩阵Bi将单元坐标映射到系统坐标，即

qei=Biqe，

(28)

将式(28)代入废塑料薄膜系统的质量、刚度及外力矩阵中，得到

(29)

利用牛顿方程建立废塑料薄膜系统的动力学方程，即

(30)

利用Newmark-β结合牛顿迭代方式来求解薄膜系统的动力学方程，求解流程如图2所示。

图2 动力学方程求解流程

3 数值仿真与对比分析

为了利用MATLAB软件对绝对节点坐标法下薄膜的动力学模型进行自由振动分析，首先应将方程线性化，即

(31)

式中Kτ为系统的切向刚度矩阵，在系统处于静平衡状态时Kτ为系统广义坐标的函数。系统静平衡时的方程为

QK+Kuq=Q，

(32)

可求得系统的切向刚度矩阵Kτ为

(33)

当系统为自由振动时，切向刚度为

(34)

则固有频率以及对应的振型可表示为

(Kτ-ω2M)φ=0，

(35)

式中：ω为固有频率；φ为固有振型。

根据上述研究，对基于绝对节点坐标法的薄膜动力学模型进行自由振动模态分析。选择薄膜尺寸为60 mm×60 mm×0.03 mm，薄膜材料为中密度的聚乙烯材料，弹性模量E=1.72×108Pa，密度ρ=1 390 kg/m3，泊松比v=0.439。利用MATLAB软件编程计算得到的薄膜自由振动第7阶至第15阶振型如图3所示。