基于柔性多體動力學的瓦楞成型系統(tǒng)建模與仿真研究

曹大志 ， 趙治華 ， 強洪夫 ， 任革學

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

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

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

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

1 柔性多體動力學理論

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

1.1 絕對節(jié)點坐標的殼單元描述

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

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

此單元含有四個節(jié)點，單元的廣義坐標表示為：

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

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

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

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

基于虛功原理，無約束的單元的動力學方程可表達為：

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

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

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

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

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

由赫茲理論可得碰撞點上的摩擦碰撞力f為：

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

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

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

1.3 基于絕對節(jié)點坐標的柔性多體系統(tǒng)動力學方程及求解

基于絕對節(jié)點坐標柔性體的總動能和總勢能可由單元動能式（4）和勢能式（5）累加而成：

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

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

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

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

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

圖3 瓦楞機Fig.3 The corrugating machine

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

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

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

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

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

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

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

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

2.2 瓦楞機工況的邊界條件

邊界條件主要包含有運動約束和載荷，其中運動約束如圖4所示，載荷邊界條件如圖6所示。

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

3 仿真結(jié)果分析

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

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

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

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

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

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

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

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

4 結(jié)論

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

［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］石 望，劉錦陽.大變形彈塑性梁的剛-柔耦合動力學特性［J］.上海交通大學學報，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］田 強，張云清，陳立平，等.含分數(shù)阻尼特性元件的多體系統(tǒng)動力學研究［J］.力學學報，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］蔡逢春，臧峰剛，葉獻輝，等.基于絕對節(jié)點坐標法的輸流管道非線性動力學分析［J］.振動與沖擊，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é)點坐標的柔性體碰撞仿真［J］.清華大學學報（自然科學版），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.