王婷，楊先海，李永波

(山東理工大學(xué) 機(jī)械工程學(xué)院，山東 淄博 255049)

近年來，隨著塑料制品的不斷使用，塑料垃圾也越來越多，廢塑料污染問題已成為一個(gè)非常嚴(yán)重的環(huán)境問題。在解決生態(tài)環(huán)境污染的問題中，整治廢棄塑料所帶來的環(huán)境污染是最重要的方面之一，而將廢塑料垃圾從混合生活垃圾中分離出來是首先應(yīng)該被實(shí)施并完成的。廢塑料薄膜等柔性材料在實(shí)際分選過程中本身的振動(dòng)特性會(huì)直接影響其運(yùn)動(dòng)軌跡，從而影響分選純度和效率；因此,需對(duì)廢塑料薄膜進(jìn)行振動(dòng)特性分析，進(jìn)而提高廢塑料薄膜分選的精度和效率。

現(xiàn)有的薄膜振動(dòng)特性研究大多利用線性振動(dòng)動(dòng)力學(xué)方程來進(jìn)行分析[1-2]。Tian等[3]在假定來流為均勻不可壓縮理想勢(shì)流的前提下，研究了空氣中薄膜振動(dòng)產(chǎn)生的附加空氣質(zhì)量以及氣動(dòng)力對(duì)膜振動(dòng)頻率和振幅的影響；Zhang等[4]提出了高功率脈沖激光器產(chǎn)生等同于聲激勵(lì)的理想脈沖激勵(lì)對(duì)薄膜進(jìn)行振動(dòng)激勵(lì)的方法來研究薄膜振動(dòng)特性，并利用實(shí)驗(yàn)驗(yàn)證了該方法的可行性；Si等[5]采用有限元方法對(duì)變壓器油膜的振動(dòng)特性進(jìn)行了數(shù)值研究，并考慮了粘度對(duì)能量損失的影響；邵明月等[6]、武吉梅等[7]基于Von Karman薄板理論推導(dǎo)出軸向運(yùn)動(dòng)薄膜的非線性振動(dòng)方程，得出薄膜穩(wěn)定工作區(qū)間和發(fā)散失穩(wěn)區(qū)間；劉明君[8]以薄膜的褶皺形態(tài)為新的平衡狀態(tài)進(jìn)行動(dòng)力學(xué)分析，建立了不同載荷下褶皺薄膜的動(dòng)力學(xué)模型，并準(zhǔn)確地描述出褶皺薄膜的振動(dòng)特性。上述計(jì)算分析方法都有效地描述出不同薄膜的振動(dòng)特性，但計(jì)算大多較為復(fù)雜，效率較低。

絕對(duì)節(jié)點(diǎn)坐標(biāo)法(absolute nodal coordinate formulation，ANCF)是Shabana教授為解決柔性體的變形問題而提出的[9]，與傳統(tǒng)有限元方法不同，該方法中節(jié)點(diǎn)坐標(biāo)選取的是全局坐標(biāo)系下的位置及其梯度，從而避開了小轉(zhuǎn)角的約束[10]。一些學(xué)者基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法討論了不同維度下的單元模型，進(jìn)而求解動(dòng)力學(xué)問題[11-14]。上述的研究集中于對(duì)非柔性體單元模型的建模仿真和彈性力的研究，將絕對(duì)節(jié)點(diǎn)坐標(biāo)法應(yīng)用于塑料薄膜振動(dòng)特性的研究較少。

由于塑料薄膜的振動(dòng)特性對(duì)分選效率的影響較大，本文在利用振動(dòng)對(duì)廢塑料薄膜進(jìn)行分選的基礎(chǔ)上，基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立薄膜單元的三維柔性模型，分析廢塑料薄膜的振動(dòng)特性，并通過MATLAB、ANSYS等軟件驗(yàn)證該方法的正確性，以期為后續(xù)的振動(dòng)分選提供理論基礎(chǔ)。

1 基于ANCF的薄膜單元模型

ANCF通常用于求解梁結(jié)構(gòu)的變形和位移問題，在傳統(tǒng)算法中，由于薄膜的厚度較小，通常忽略厚度的影響，本文在薄板單元模型的基礎(chǔ)上研究薄膜的振動(dòng)特性，同時(shí)考慮x、y和z的梯度，進(jìn)而增加了計(jì)算精度。三維柔性薄膜單元模型如圖1所示，圖中O-XYZ為全局坐標(biāo)系，與單元模型相連結(jié)，P0和P分別為初始構(gòu)型與當(dāng)前構(gòu)型時(shí)單元模型上的質(zhì)點(diǎn)，r0和r為節(jié)點(diǎn)的位移向量，?rij/?x、?rij/?y和?rij/?z均為節(jié)點(diǎn)坐標(biāo)。

圖1 三維柔性薄膜單元模型

三維柔性薄膜單元中每個(gè)節(jié)點(diǎn)含有3個(gè)位置坐標(biāo)和9個(gè)位移斜率坐標(biāo)，單元的自由度為48，因此單元上任意一點(diǎn)的位置坐標(biāo)可利用節(jié)點(diǎn)坐標(biāo)和形函數(shù)來表示，即

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

(1)

式中：S為單元的形函數(shù)；x和y為任意的局部坐標(biāo)；qe(t)為單元節(jié)點(diǎn)坐標(biāo)。

單元的節(jié)點(diǎn)坐標(biāo)可定義為

(2)

其中任一節(jié)點(diǎn)坐標(biāo)可表示為

(3)

形函數(shù)可表示為

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

(4)

各個(gè)分量如下：

1.1 單元質(zhì)量矩陣

薄膜單元的質(zhì)量矩陣可以利用薄膜的動(dòng)能來獲得，薄膜單元的動(dòng)能可以表示為

(6)

根據(jù)節(jié)點(diǎn)速度矢量，可以求得絕對(duì)速度矢量的表達(dá)式

(7)

代入薄膜單元?jiǎng)幽鼙磉_(dá)式中可得

(8)

從而得到薄膜單元質(zhì)量矩陣

(9)

1.2 彈性力矩陣

由于薄膜為柔性材料，利用傳統(tǒng)的線性求解方法已無(wú)法滿足單元模型對(duì)于材料屬性的精確描述，因此在非線性連續(xù)介質(zhì)力學(xué)的基礎(chǔ)上，利用Neo-Hookean本構(gòu)模型來求解三維柔性薄膜單元的彈性力矩陣，進(jìn)而求解出單元的剛度矩陣。將廢塑料薄膜材料變形視為各向同性，基于Neo-Hookean本構(gòu)模型的薄膜單元應(yīng)變能密度函數(shù)可以表示為

(10)

式中：變形張量不變量I1可以定義為右柯西-格林變形張量C的跡，即

I1=tr(C)；

(11)

J為位置坐標(biāo)求偏導(dǎo)的矩陣，可以表示為

J=det(J)=|J|。

(12)

將應(yīng)變能密度函數(shù)對(duì)節(jié)點(diǎn)坐標(biāo)進(jìn)行微分，可以得到三維柔性薄膜單元的彈性力矩陣

(13)

根據(jù)右柯西-格林變形張量定義得

因此C的跡為

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

(15)

式(15)對(duì)節(jié)點(diǎn)坐標(biāo)求微分可得

(16)

薄膜單元的變形梯度為

(17)

將式(17)代入J的矩陣行列式可得

J=S1xqeS2yqeS3zqe+S3xqeS1yqeS2zqe+

S2xqeS3yqeS1zqe-S1xqeS3yqeS2zqe-

S2xqeS1yqeS3zqe-S3xqeS2yqeS1zqe。

(18)

由于柔性薄膜材料的泊松比與不可壓縮材料的泊松比極為接近，因此在以往的計(jì)算中都將其視為不可壓縮材料。但在實(shí)際操作條件下，材料的不可壓縮條件是無(wú)法滿足的，上述假設(shè)會(huì)影響計(jì)算結(jié)果的精確性，因此需在推導(dǎo)公式時(shí)增加一個(gè)補(bǔ)償方程。薄膜單元的應(yīng)變能密度函數(shù)可以表示為

(19)

式中

(20)

取系數(shù)k=1×109N/m2。

求得薄膜單元的彈性力矩陣為

(21)

薄膜單元的彈性力矩陣還可表示為

(22)

式中K為薄膜單元的剛度矩陣，可表示為

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

(23)

式中

1.3 單元廣義外力

當(dāng)有外力F作用于薄膜單元上時(shí)，該力的虛功可表示為

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

(25)

則薄膜單元的廣義外力Qj為

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

(26)

薄膜結(jié)構(gòu)中所有單元的廣義外力為

Q=∑n1×n2HTQj。

(27)

2 廢塑料薄膜動(dòng)力學(xué)方程及求解

由于薄膜單元被離散為多個(gè)單元，假設(shè)系統(tǒng)中包含n個(gè)單元，第i個(gè)單元的單元節(jié)點(diǎn)坐標(biāo)為qei，利用布爾矩陣Bi將單元坐標(biāo)映射到系統(tǒng)坐標(biāo)，即

qei=Biqe，

(28)

將式(28)代入廢塑料薄膜系統(tǒng)的質(zhì)量、剛度及外力矩陣中，得到

(29)

利用牛頓方程建立廢塑料薄膜系統(tǒng)的動(dòng)力學(xué)方程，即

(30)

利用Newmark-β結(jié)合牛頓迭代方式來求解薄膜系統(tǒng)的動(dòng)力學(xué)方程，求解流程如圖2所示。

圖2 動(dòng)力學(xué)方程求解流程

3 數(shù)值仿真與對(duì)比分析

為了利用MATLAB軟件對(duì)絕對(duì)節(jié)點(diǎn)坐標(biāo)法下薄膜的動(dòng)力學(xué)模型進(jìn)行自由振動(dòng)分析，首先應(yīng)將方程線性化，即

(31)

式中Kτ為系統(tǒng)的切向剛度矩陣，在系統(tǒng)處于靜平衡狀態(tài)時(shí)Kτ為系統(tǒng)廣義坐標(biāo)的函數(shù)。系統(tǒng)靜平衡時(shí)的方程為

QK+Kuq=Q，

(32)

可求得系統(tǒng)的切向剛度矩陣Kτ為

(33)

當(dāng)系統(tǒng)為自由振動(dòng)時(shí)，切向剛度為

(34)

則固有頻率以及對(duì)應(yīng)的振型可表示為

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

(35)

式中：ω為固有頻率；φ為固有振型。

根據(jù)上述研究，對(duì)基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法的薄膜動(dòng)力學(xué)模型進(jìn)行自由振動(dòng)模態(tài)分析。選擇薄膜尺寸為60 mm×60 mm×0.03 mm，薄膜材料為中密度的聚乙烯材料，彈性模量E=1.72×108Pa，密度ρ=1 390 kg/m3，泊松比v=0.439。利用MATLAB軟件編程計(jì)算得到的薄膜自由振動(dòng)第7階至第15階振型如圖3所示。

圖3 MATLAB求解振型圖

選用同樣的薄膜材料，利用有限元軟件ANSYS對(duì)薄膜進(jìn)行自由振動(dòng)模態(tài)分析，得到薄膜自由振動(dòng)的第7階至第15階振型如圖4所示，固有頻率如圖5所示。由于薄膜為自由振動(dòng)，所以前6階為剛體運(yùn)動(dòng)，固有頻率近似為0;隨著階數(shù)的提高，薄膜單元振動(dòng)形式變得更為復(fù)雜，固有頻率會(huì)逐漸增大。

圖4 薄膜自由振動(dòng)振型圖

圖5 固有頻率圖

將MATLAB軟件仿真結(jié)果與ANSYS軟件計(jì)算結(jié)果進(jìn)行對(duì)比分析，結(jié)果見表1。

從表1中數(shù)據(jù)可以看出，基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立的薄膜單元模型的自由振動(dòng)頻率與ANSYS傳統(tǒng)的理論值基本一致，誤差較小，平均誤差約為3.21%。綜上所述，基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立的單元模型可以有效地分析廢塑料薄膜的固有特性，證明利用絕對(duì)節(jié)點(diǎn)坐標(biāo)法進(jìn)行薄膜振動(dòng)特性分析是可行的。

表1 薄膜自由振動(dòng)頻率對(duì)比

4 結(jié)論

通過建立薄膜單元模型并進(jìn)行仿真分析計(jì)算，研究了廢塑料薄膜的振動(dòng)規(guī)律，結(jié)論如下：

1)利用絕對(duì)節(jié)點(diǎn)坐標(biāo)法，提出了三維柔性薄膜單元模型，通過求解矩陣，得到了廢塑料薄膜的動(dòng)力學(xué)方程表達(dá)式。

2)利用MATLAB、ANSYS等軟件對(duì)該方法進(jìn)行的驗(yàn)證結(jié)果表明，基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法求解的薄膜振動(dòng)頻率與傳統(tǒng)理論值相比平均誤差約為3.21%，具有可行性。

3)利用絕對(duì)節(jié)點(diǎn)坐標(biāo)法求得的廢塑料薄膜振動(dòng)規(guī)律，可應(yīng)用于廢塑料薄膜分選設(shè)備的設(shè)計(jì)優(yōu)化中，對(duì)提高廢塑料薄膜的分選精度具有參考意義。