曹大志 ， 趙治華 ， 強(qiáng)洪夫 ， 任革學(xué)

（1.清華大學(xué) 工程力學(xué)系，北京 100084；2.第二炮兵工程大學(xué) 201室，西安 710025）

在一些如傳送機(jī)械、太陽能帆板伸展機(jī)構(gòu)、機(jī)械臂等系統(tǒng)中，柔性體的變形對系統(tǒng)的動(dòng)力學(xué)行為產(chǎn)生很大的影響，傳統(tǒng)的多剛體系統(tǒng)動(dòng)力學(xué)理論已經(jīng)不能滿足人們對設(shè)備越來越高的精度要求。此外隨著材料技術(shù)和制造技術(shù)的進(jìn)步，越來越多的柔性部件應(yīng)用到存在大位移和大變形的系統(tǒng)中。于是，柔性體的大變形多體動(dòng)力學(xué)問題不斷出現(xiàn)，傳統(tǒng)的基于小變形理論建立動(dòng)力學(xué)模型的方法已不能得到這類大變形柔性多體動(dòng)力學(xué)問題的精確解。Shabana［1-2］提出的絕對節(jié)點(diǎn)坐標(biāo)法可以很好的解決這類大變形柔性多體動(dòng)力學(xué)問題，并且在多體動(dòng)力學(xué)領(lǐng)域受到廣泛的重視。不同于傳統(tǒng)的浮動(dòng)坐標(biāo)系法，共旋坐標(biāo)法或相對節(jié)點(diǎn)位移法等，基于絕對節(jié)點(diǎn)坐標(biāo)的單元節(jié)點(diǎn)的廣義坐標(biāo)直接采用的是空間坐標(biāo)及其梯度，因此在解決眾多動(dòng)力學(xué)問題中具有很大優(yōu)勢。絕對節(jié)點(diǎn)坐標(biāo)法的出現(xiàn)促進(jìn)了柔性多體動(dòng)力學(xué)與有限元理論的進(jìn)一步整合，基于絕對節(jié)點(diǎn)坐標(biāo)法的梁［3］，板/殼［4］和三維實(shí)體［5］等單元陸續(xù)的提出。在我國對絕對節(jié)點(diǎn)坐標(biāo)法的研究起步較晚。劉錦陽等［6-11］基于絕對節(jié)點(diǎn)坐標(biāo)做了以下工作：考慮了熱載荷與結(jié)構(gòu)變形的耦合，分析了復(fù)合材料以及柔性板的多體系統(tǒng)動(dòng)力學(xué)；提出一套混合坐標(biāo)體系，采用絕對節(jié)點(diǎn)坐標(biāo)法進(jìn)行了驗(yàn)證；此外還考慮了材料的彈塑性，求解了梁的剛?cè)狁詈蟿?dòng)力學(xué)特性。張?jiān)魄宓龋?2-14］基于絕對節(jié)點(diǎn)坐標(biāo)法研究了含分?jǐn)?shù)階阻尼特性元件的多體系統(tǒng)建模和求解問題，并針對約束方程的違約問題提出高效的數(shù)值計(jì)算方法。蔡逢春等［15］基于絕對節(jié)點(diǎn)坐標(biāo)法建立了一種新的一維二節(jié)點(diǎn)輸流管道單元，計(jì)算和分析了不同邊界條件下的輸流管道的非線性動(dòng)力學(xué)行為。

在金屬冶煉或印刷等軋制工業(yè)的加工和處理流程中，帶材會自行偏離生產(chǎn)線的中心，向輥?zhàn)拥囊贿呉苿?dòng)，這稱為帶材跑偏。國外眾多學(xué)者對軋制過程的出現(xiàn)的帶材跑偏的動(dòng)力學(xué)行為進(jìn)行了研究。Shelton等［16］基于歐拉梁假設(shè)推導(dǎo)得到了帶材橫向位移的數(shù)學(xué)模型。Laukkanen［17］提出四節(jié)點(diǎn)的 Mindlin-Reissner彎曲板單元，并對帶紙進(jìn)行了靜力學(xué)和動(dòng)力學(xué)分析。Kulachenko等［18］基于混合坐標(biāo)法提出四節(jié)點(diǎn)20自由度的殼單元，采用縮減積分的辦法克服了剪切鎖死的問題。Yu等［19］對帶鋼的跑偏糾偏裝置進(jìn)行了多體動(dòng)力學(xué)仿真，此外還有學(xué)者對帶材的跑偏行為和糾偏裝置進(jìn)行了實(shí)驗(yàn)研究［20］。

大型復(fù)雜結(jié)構(gòu)本身受諸多因素的限制，其實(shí)驗(yàn)往往很難進(jìn)行，因此有必要借助先進(jìn)的仿真手段來得到較精確的系統(tǒng)模型。現(xiàn)代柔性多體動(dòng)力學(xué)理論特別是絕對節(jié)點(diǎn)坐標(biāo)法的發(fā)展為復(fù)雜多體系統(tǒng)的建模和仿真提供了有利的工具。本文借助現(xiàn)代柔性多體動(dòng)力學(xué)理論，以某型號瓦楞機(jī)的瓦楞成型系統(tǒng)為對象，建立了系統(tǒng)的柔性多體動(dòng)力學(xué)仿真模型，并對傳送帶的跑偏等動(dòng)力學(xué)行為進(jìn)行了分析。

1 柔性多體動(dòng)力學(xué)理論

本節(jié)簡要介紹基于絕對節(jié)點(diǎn)坐標(biāo)殼單元的有限元格式，赫茲碰撞理論和柔性多體系統(tǒng)方程的求解。

1.1 絕對節(jié)點(diǎn)坐標(biāo)的殼單元描述

Dufva等［14］基于絕對節(jié)點(diǎn)坐標(biāo)體系提出了一種可精確描述任意剛體位移和大變形的四節(jié)點(diǎn)36自由度的殼單元。如圖1所示。其中a、b和t分別為板殼單元的長度、寬度和高度。

圖1 基于絕對節(jié)點(diǎn)坐標(biāo)的四節(jié)點(diǎn)殼單元Fig.1 The ANC-based four nodes shell element

此單元含有四個(gè)節(jié)點(diǎn)，單元的廣義坐標(biāo)表示為：

其中ei為第i個(gè)節(jié)點(diǎn)的廣義坐標(biāo)，包含有空間位置向量ri及其在中性面上兩個(gè)方向的一階導(dǎo)數(shù)和，故每個(gè)節(jié)點(diǎn)共有9個(gè)自由度：

單元中任一點(diǎn)的絕對坐標(biāo)r可表示為單元節(jié)點(diǎn)坐標(biāo)e的有限元插值形式：

其中：ξ為實(shí)體單元母單元的自然坐標(biāo)，N為單元的形函數(shù)［4］。單元的動(dòng)能可表示為：

其中：ε和κ分別為殼單元面內(nèi)的膜應(yīng)變和橫向曲率，具體定義見文獻(xiàn)［4］。Eε，Eκ分別為殼單元的彈性矩陣，并且能夠描述材料的正交各向異性。

基于虛功原理，無約束的單元的動(dòng)力學(xué)方程可表達(dá)為：

式中：Fe為單元彈性力列陣，Qe為廣義外力列陣。

1.2 碰撞檢測和赫茲碰撞理論

接觸碰撞分析是機(jī)械動(dòng)力學(xué)的重要環(huán)節(jié)，很多情況下，機(jī)械中各零件之間力的傳遞是通過兩零件的接觸碰撞來實(shí)現(xiàn)的。在瓦楞機(jī)工作過程中，柔性的傳送帶和剛性的各支撐輥之間存在大面積的接觸碰撞，因此需要采用合理、高效的接觸算法對其進(jìn)行建模。

本文采用在單元內(nèi)布置若干檢測點(diǎn)的方法，將柔性傳送帶與剛體支撐輥側(cè)面的面-面之間的碰撞檢測簡化為點(diǎn)-面之間的檢測，如圖2所示。

圖2 殼單元與剛體的碰撞Fig.2 The shell-to-rigid contact

由赫茲理論可得碰撞點(diǎn)上的摩擦碰撞力f為：

其中n為碰撞點(diǎn)的法線方向，如圖2所示。fn、fτ和τ分別為法向碰撞力、摩擦力和切向單位向量，其具體計(jì)算公式見文獻(xiàn)［21］。

利用虛功原理，作用在檢測點(diǎn)p處的碰撞摩擦力fp可以轉(zhuǎn)化為單元廣義節(jié)點(diǎn)力Qp：

其中N是單元的形函數(shù)。

1.3 基于絕對節(jié)點(diǎn)坐標(biāo)的柔性多體系統(tǒng)動(dòng)力學(xué)方程及求解

基于絕對節(jié)點(diǎn)坐標(biāo)柔性體的總動(dòng)能和總勢能可由單元?jiǎng)幽苁剑?）和勢能式（5）累加而成：

其中n和ne分別代表節(jié)點(diǎn)個(gè)數(shù)和絕對節(jié)點(diǎn)坐標(biāo)單元的個(gè)數(shù)。利用第一類Lagrange方程，可得到含約束的基于絕對節(jié)點(diǎn)坐標(biāo)的多體系統(tǒng)動(dòng)力學(xué)方程為：

式中：q為柔性的廣義坐標(biāo)；λ為拉格郎日乘子列陣；C（q，t）為約束方程；Cq為約束方程對q的偏導(dǎo)數(shù)矩陣；M為柔性體的質(zhì)量陣；Q為廣義外力列陣；F為柔性體的彈性力列陣，是柔性體總勢能對q的偏導(dǎo)數(shù)：

式（10）是一組典型的微分代數(shù)方程組（DAEs），可由隱式的向后差分法［12］（Backward Differentiation Formulas，BDF）求解。

2 瓦楞成型系統(tǒng)的多體動(dòng)力學(xué)建模

瓦楞機(jī)是生產(chǎn)瓦楞紙板的機(jī)構(gòu)，它是由傳送帶，張力輥，上下瓦楞輥等部件組成（見圖3）。瓦楞機(jī)的基本工作原理：經(jīng)預(yù)熱和預(yù)濕后的瓦楞原紙（紙2）進(jìn)入上下瓦楞輥的通道受擠壓熨燙成型為瓦楞芯紙。由涂膠輥上膠后隨瓦楞輥繼續(xù)運(yùn)轉(zhuǎn)，與同時(shí)經(jīng)過預(yù)熱及溫控后到達(dá)的面紙1在繃緊的傳送帶的壓力作用下進(jìn)行熱壓復(fù)合干燥，由此高速運(yùn)轉(zhuǎn)，生產(chǎn)出二層瓦楞紙板。

圖3 瓦楞機(jī)Fig.3 The corrugating machine

如果兩個(gè)旋轉(zhuǎn)張力輥由于生產(chǎn)或是安裝誤差導(dǎo)致不平行，或是輥?zhàn)愚D(zhuǎn)動(dòng)的方向與傳送帶運(yùn)行的方向不垂直，就會造成原紙跟隨傳送帶的跑偏、瓦楞芯紙與面紙未對齊、紙板的各種翹曲等各種影響生產(chǎn)的異常現(xiàn)象。因此有必要掌握瓦楞機(jī)的動(dòng)力學(xué)行為。

2.1 瓦楞成型系統(tǒng)的尺寸及仿真模型

瓦楞成型系統(tǒng)是瓦楞機(jī)的核心，主要是由傳送帶與三個(gè)支撐輥（兩個(gè)張力輥和上瓦楞輥）組成，模型如圖4所示。其中將傳送帶考慮成柔性體，將支撐輥考慮成剛體，在本文中將瓦楞輥考慮成表面光滑的圓柱體，不考慮瓦楞輥的鋸齒。模型參數(shù)見表1，其中碰撞相關(guān)參數(shù)由赫茲碰撞理論通過實(shí)驗(yàn)而得到。

圖4 瓦楞成型系統(tǒng)Fig.4 Modeling image of corrugating machine

表1 模型參數(shù)Tab.1 Model＇s parameters

采用1.1節(jié)中所描述的基于絕對節(jié)點(diǎn)坐標(biāo)法的36自由度的殼單元對傳送帶進(jìn)行有限元離散，如圖5所示。值得注意的是基于絕對節(jié)點(diǎn)坐標(biāo)單元的位移模式相對傳統(tǒng)的有限元含有更高階的完備多項(xiàng)式，因此網(wǎng)格劃分的可以相對稀疏。經(jīng)收斂性測試實(shí)驗(yàn)，將傳送帶劃分為40×1的網(wǎng)格密度可保證數(shù)值解的收斂。

傳送帶與三個(gè)支撐輥的碰撞采用1.2節(jié)所述的點(diǎn)-面檢測的辦法進(jìn)行建模。由于基于絕對節(jié)點(diǎn)坐標(biāo)的單元的網(wǎng)格密度相對稀疏，因此在每個(gè)板殼單元內(nèi)部布置7×7=49個(gè)檢測點(diǎn)，考慮到每個(gè)檢測點(diǎn)需要和三個(gè)支撐輥的圓柱表面進(jìn)行點(diǎn)-面碰撞檢測，因此在模型中共建立49×40×3=5 880個(gè)點(diǎn)-面碰撞檢測對。

圖5 傳送帶的有限元模型Fig.5 FEM model of resin belt

2.2 瓦楞機(jī)工況的邊界條件

邊界條件主要包含有運(yùn)動(dòng)約束和載荷，其中運(yùn)動(dòng)約束如圖4所示，載荷邊界條件如圖6所示。

圖6 載荷邊界條件Fig.6 The load boundary condition

3 仿真結(jié)果分析

采用自主研發(fā)的多體動(dòng)力學(xué)程序?qū)ν呃愠尚拖到y(tǒng)的仿真模型進(jìn)行動(dòng)力學(xué)分析，仿真時(shí)間取為15 s。采用基于浮動(dòng)坐標(biāo)系（floating frame of reference formulation，F(xiàn)FRF）的小變形的殼單元建立相同的動(dòng)力學(xué)模型，并進(jìn)行對比分析。值得注意的是，由于基于小變形的殼不能較準(zhǔn)確的求解大變形的問題，因此本文采用基于FFRF的小變形單元的網(wǎng)格更加精細(xì)。

從瓦楞機(jī)的工況可知，由于編號為2的張力輥兩端所受力大小不同，因此它與編號為1的張力輥不再平行，從而導(dǎo)致了傳送帶受力不均，發(fā)生跑遍?；贏NCF和FFRF兩種不同建模方法下，得到編號為2的張力輥沿X方向的位移和傳送帶中心（見圖5）沿Z方向的偏心位移，如圖7和圖8所示。從兩圖的對比可知，計(jì)算結(jié)果基本一致，從而也驗(yàn)證了基于絕對節(jié)點(diǎn)坐標(biāo)建模方法的準(zhǔn)確性。此外，采用基于絕對節(jié)點(diǎn)坐標(biāo)的建模方法，需要的單元個(gè)數(shù)少，相對浮動(dòng)坐標(biāo)系下的小變形單元建模的方法，本文中基于絕對節(jié)點(diǎn)坐標(biāo)系的建模方法具有更高的計(jì)算效率和計(jì)算精度。

在瓦楞機(jī)的工作過程中，瓦楞芯紙與面紙的粘合是通過上瓦楞輥與傳送帶的接觸碰撞實(shí)現(xiàn)的。適合的碰撞力可以在不壓碎瓦楞芯紙的情況下，讓面紙與瓦楞更充分的粘合。利用傳送帶上的檢測點(diǎn)與圓柱側(cè)面的解析檢測方法，能準(zhǔn)確描述帶與支撐輥之間的面與面的摩擦碰撞，從而較準(zhǔn)確的模擬瓦楞機(jī)在壓制過程中的工作狀態(tài)。計(jì)算得到的摩擦碰撞力分布如圖9所示。

在瓦楞機(jī)的工作過程中，研究傳送帶的最大Mises應(yīng)力或最大主應(yīng)力等，可以分析傳送帶的強(qiáng)度是否滿足要求，并且傳送帶的應(yīng)力分布規(guī)律也是進(jìn)行其疲勞分析的基礎(chǔ)。圖10為帶中心Mises應(yīng)力的時(shí)間歷程曲線，圖11為傳送帶在穩(wěn)定工作過程的應(yīng)力分布圖。從圖11可知由于編號為2的張力輥兩端受力不均，因此導(dǎo)致帶的應(yīng)力分布同樣不均。在受力較大一端的應(yīng)力明顯比受力小的一端大。帶的內(nèi)部應(yīng)變分布的不均也直接導(dǎo)致帶子朝張力小的方向移動(dòng)。

圖9 碰撞力顯示Fig.9 The display of contact force

圖12（a）～（e）顯示的是傳送帶的變形歷程，在t=0 s時(shí)，傳送帶處于未變形的構(gòu)型，是一個(gè)圓柱殼。在t=0～1 s時(shí)間段，由于受張力輥的拉伸以及上瓦楞輥的擠壓開始繃緊，然后在逐漸增大的張力的作用下發(fā)生彈性變形，隨后受瓦楞輥摩擦接觸力驅(qū)動(dòng)傳送帶開始旋轉(zhuǎn)，如圖12（d）和（e）中所示。圖12（f）為傳送帶中心點(diǎn)在三維空間中的運(yùn)動(dòng)軌跡圖，可明顯看出傳送帶經(jīng)歷了大位移、大轉(zhuǎn)動(dòng)的運(yùn)動(dòng)過程。

圖12 傳送帶變形歷程Fig.12 Snapshots of resin belt animation

對機(jī)械的動(dòng)力學(xué)行為分析最終的目的是對系統(tǒng)進(jìn)行動(dòng)力學(xué)控制。由于現(xiàn)代造紙以及軋制過程的控制非常復(fù)雜，設(shè)計(jì)到速度、張力、壓力等大量物理參數(shù)，以及結(jié)構(gòu)的變形及部件之間的摩擦碰撞，因此是一個(gè)非線性、強(qiáng)耦合的復(fù)雜過程?；诮^對節(jié)點(diǎn)坐標(biāo)的建模方法，可以研究不同邊界條件下如施加的張力，帶材的驅(qū)動(dòng)速度等對帶材的動(dòng)力學(xué)行為的影響。此外還可方便地獲得如偏移量，碰撞力，應(yīng)力分布等動(dòng)力學(xué)響應(yīng)參數(shù)，為下一步的動(dòng)力學(xué)控制打下基礎(chǔ)，為控制參數(shù)的選擇和優(yōu)化提供依據(jù)。

4 結(jié)論

瓦楞機(jī)軋制紙板的過程是個(gè)非線性、強(qiáng)耦合的復(fù)雜過程，本文采用基于絕對節(jié)點(diǎn)坐標(biāo)的板殼單元離散傳送帶，既精確地描述了傳送帶的大轉(zhuǎn)動(dòng)和大位移的剛體運(yùn)動(dòng)，又描述了傳送帶的非線性變形。采用點(diǎn)-面檢測的方法和赫茲碰撞理論描述了支撐輥與傳送帶之間的大面積接觸碰撞、強(qiáng)耦合變形的動(dòng)力學(xué)行為。在自主研發(fā)的多體動(dòng)力學(xué)程序中計(jì)算得到柔性和剛性部件的動(dòng)位移，碰撞力大小以及應(yīng)力分布等動(dòng)響應(yīng)，為進(jìn)一步的動(dòng)力學(xué)控制和系統(tǒng)的設(shè)計(jì)和改進(jìn)提供了理論的基礎(chǔ)。綜上，基于絕對節(jié)點(diǎn)坐標(biāo)的建模方法在處理柔性多體大變形動(dòng)力學(xué)問題具有很大的潛力。

［1］ Shabana A A.Flexible multibody dynamics：review of past and recent developments［J］.Multibody System Dynamics，1997，1（2）：189-222.

［2］Shabana A A.Definition of the slopes and the finite element absolute nodal coordinate formulation［J］.Multibody System Dynamics，1997，1（3）：339-348.

［3］Shabana A A，Yakoub R Y.Three dimensional absolute nodal coordinate formulation for beam elements：theory［J］.Journal of Mechanical Design，2001，123（4）：606-613.

［4］Dufva K，Shabana A A.Analysis of thin plate structures using the absolute nodal coordinate formulation［J］.Proceedings of the Institution of Mechanical Engineers Part K-Journal of Multi-Body Dynamics，2005，219（4）：345-355.

［5］ Gerstmayr J，Schoberl J.A 3D finite element method for flexible multibody systems［J］.Multibody System Dynamics，2006，15（4）：309-324.

［6］Liu J Y，Lu H.Thermal effect on the deformation of a flexible beam with large kinematical driven overall motions［J］.European Journal of Mechanics-A/Solids，2006，26（1）：137-151.

［7］ Liu J Y，Lu H.Rigid-flexible coupling dynamics of threedimensional hub-beamssystem［J］. MultibodySystem Dynamics，2007，18（4）：487-510.

［8］ Liu J Y，Hong J Z.Nonlinear formulation for flexible multibody system with large deformation［J］.Acta Mechanica Sinica，2007，23（1）：111-119.

［9］石 望，劉錦陽.大變形彈塑性梁的剛-柔耦合動(dòng)力學(xué)特性［J］.上海交通大學(xué)學(xué)報(bào)，2011，45（10）：1469-1474.

SHI Wang， LIU Jin-yang. Investigation on rigid-flexible coupling dynamics for elasto-plastic beam with large deformation［J］.Journal of Shanghai Jiaotong University，2011，45（10）：1469-1474.

［10］ Liu J Y，Hong J Z，Cui L.An exact nonlinear hybridcoordinate formulation for flexible multibody systems［J］.Acta Mechanica Sinica，2007，23（6）：699-706.

［11］Liu Z Y，Hong J Z，Liu J Y.Finite element formulation for dynamics of planar flexible multi-beam system［J］.Multibody System Dynamics，2009，22（1）：1-26.

［12］ Zhang Y Q，Tian Q，Chen L P，et al.Simulation of a viscoelastic flexible multibody system using absolute nodal coordinate and fractional derivative methods［J］.Multibody System Dynamics，2009，21（3）：281-303.

［13］田 強(qiáng)，張?jiān)魄?，陳立平，?含分?jǐn)?shù)阻尼特性元件的多體系統(tǒng)動(dòng)力學(xué)研究［J］.力學(xué)學(xué)報(bào)，2009，41（6）：920-928.

TIAN Qiang，ZHANG Yun-qing，CHEN Li-ping，et al.Dynamics research on the multibody system with fractional derivative damper［J］.Chinese Journal of Theoretical and Applied Mechanics，2009，41（6）：920-928.

［14］ Tian Q，Chen L P，Zhang Y Q，et al.An efficient hybrid method for multibody dynamics simulation based on absolute nodal coordinate formulation［J］.Journal of Computational and Nonlinear Dynamics，2009，4（2）：021009-021014.

［15］蔡逢春，臧峰剛，葉獻(xiàn)輝，等.基于絕對節(jié)點(diǎn)坐標(biāo)法的輸流管道非線性動(dòng)力學(xué)分析［J］.振動(dòng)與沖擊，2011，30（6）：143-146.

CAI Feng-chun，ZANG Feng-gang，YE Xian-hui，et al.Analysis of nonlinear dynamic behavior of pipe conveying fluid based on absolute nodal coordinate formulation［J］.Journal of Vibration and Shock，2011，30（6）：143-146.

［16］ Shelton J J，Reid K N.Lateral dynamics of a real moving web［J］.ASME J.Dyn.Syst.Meas.Control，1971，93（3）：180-186.

［17］ Laukkanen J.Fem analysis of a travelling web［J］.Computers＆ Structures，2002，80（24）：1827-1842.

［18］ Kulachenko A，Gradin P，Koivurova H.Modelling the dynamical behaviour of a paper web.part I［J］.Computers＆ Structures，2007，85（3-4）：131-147.

［19］ Yu L，Zhao Z，Ren G.Multibody dynamic model of web guiding system with moving web［J］.Journal of Dynamic Systems， Measurement， and Control， 2010， 132（5）：051004-051010.

［20］ Ryu J K，Lee S G，Rhim S S，et al.Simulation and experimental methods for media transport system：PartⅡ，effect of normal force on slippage of paper［J］.Journal of Mechanical Science and Technology，2005，19（1）：403-410.

［21］虞 磊，趙治華，任啟鴻，等.基于絕對節(jié)點(diǎn)坐標(biāo)的柔性體碰撞仿真［J］.清華大學(xué)學(xué)報(bào)（自然科學(xué)版），2010，50（7）：1135-1140.

YU Lei，ZHAO Zhi-hua，REN Qi-hong， et al. Contact simulations offlexiblebodiesbased on absolute nodal coordinates［J］.Journal of Tsinghua University（Science and Technology），2010，50（7）：1135-1140.

［22］ Hairer E，Wanner G.Solving ordinary differential equationsⅡ：stiff and differential-algebraic problems ［M］.Heidelberg：Springer-Verlag，1996.