基于ANSYS的汽車(chē)傳動(dòng)軸有限元分析與優(yōu)化設(shè)計(jì)【含有限元】【說(shuō)明書(shū)+CAD+UG】,含有限元,說(shuō)明書(shū)+CAD+UG,基于ANSYS的汽車(chē)傳動(dòng)軸有限元分析與優(yōu)化設(shè)計(jì)【含有限元】【說(shuō)明書(shū)+CAD+UG】,基于,ansys,汽車(chē),傳動(dòng)軸,有限元分析,優(yōu)化,設(shè)計(jì),有限元,說(shuō)明書(shū),仿單,cad 應(yīng)力為基礎(chǔ)的有限元方法應(yīng)用于靈活的曲柄滑塊機(jī)構(gòu) （多倫多大學(xué)：Y.L. Kuo .L. Cleghorn加拿大） 摘要:本文在歐拉一伯努利梁基礎(chǔ)上提出了一種新的適用于以應(yīng)力為基礎(chǔ)的有限元方法的程序。先選擇一個(gè)近似彎曲應(yīng)力的分布，然后通過(guò)一體化確定近似橫位移。該方法適用于解決靈活滑塊曲柄機(jī)構(gòu)問(wèn)題，制定的依據(jù)是歐拉-拉格朗日方程，而拉格朗日包括與動(dòng)能，應(yīng)變能有關(guān)的組件，并通過(guò)彈性橫向撓度構(gòu)成的軸向負(fù)荷的鏈接來(lái)工作。梁元模型以翻轉(zhuǎn)運(yùn)動(dòng)為基礎(chǔ)，結(jié)果表明以應(yīng)力和位移為基礎(chǔ)的有限元方法。 關(guān)鍵詞：應(yīng)力為基礎(chǔ)的有限元方法，曲柄滑塊機(jī)構(gòu)，拉格-朗日方程 1.前言 以位移為基礎(chǔ)的有限元方法通過(guò)實(shí)行假定位移補(bǔ)充能量。這種方法可能由內(nèi)部因素產(chǎn)生不連續(xù)應(yīng)力場(chǎng)，同時(shí)由于采用了低階元素，邊界條件與壓力不能得到滿足。因此，另一種被成為以應(yīng)力為基礎(chǔ)采用假定應(yīng)力的有限元方法得到了應(yīng)用和發(fā)展。Veubeke和Zienkiewicz[1-2]首先對(duì)應(yīng)力有限元素進(jìn)行了研究。之后，這種方法被廣泛用于解決應(yīng)用程序中的問(wèn)題[3-5]。此外，還有各種書(shū)籍提供更加詳細(xì)的方法[6，7]。 這一高速運(yùn)作機(jī)制采用振動(dòng)，聲輻射，協(xié)同聯(lián)結(jié)，和撓度彈性鏈接的準(zhǔn)確定位。因此，有必要分析靈活的彈塑性動(dòng)力學(xué)這一類(lèi)的問(wèn)題，而不是分析剛體動(dòng)力學(xué)。 靈活的機(jī)制是一個(gè)由無(wú)限多個(gè)自由度組成的連續(xù)動(dòng)力學(xué)系統(tǒng)，其運(yùn)動(dòng)方程是由非線性偏微分方程建立的模型，但得不到分析解決方案。Cleghorn et al[8-10] 闡述了橫向振動(dòng)上的軸向荷載對(duì)靈活四桿機(jī)構(gòu)的影響。并且通過(guò)能有效預(yù)測(cè)橫向振動(dòng)和彎曲應(yīng)力的五次多項(xiàng)式建立了一個(gè)翻轉(zhuǎn)梁?jiǎn)卧? 本文提出了一種新的方法來(lái)執(zhí)行建立在歐拉一伯努利基礎(chǔ)上的以應(yīng)力為基礎(chǔ)的有限元方法。改進(jìn)后的方法首先選定了假定應(yīng)力函數(shù)。然后通過(guò)整合假定應(yīng)力函數(shù)得到橫向位移函數(shù)。當(dāng)然，這種方法能解決沒(méi)有強(qiáng)制制約因素的應(yīng)力集中問(wèn)題。我們可以通過(guò)這種方法解決靈活曲柄滑塊機(jī)構(gòu)體系中存在的問(wèn)題。目的是通過(guò)這種方法提高準(zhǔn)確性，該系統(tǒng)存在的問(wèn)題也可以通過(guò)取代基有限元方法來(lái)解決。結(jié)果可以證明偏差比較。 2.以應(yīng)力為基礎(chǔ)的歐拉一伯努利梁 歐拉一伯努利梁的彎曲應(yīng)力與橫向位移的二階導(dǎo)數(shù)相關(guān)，也就是曲率，可以近似的看做是形函數(shù)和交點(diǎn)變量： 這里[(i)N(c)]是連續(xù)載體的形函數(shù)；{(i)?e} 是列向量的交點(diǎn)函數(shù)，y是關(guān)于中性線的橫向定位，E是楊氏模量，（i）v是橫向位移，x軸向定位函數(shù)。 由方程（1）可以推導(dǎo)出橫向位移轉(zhuǎn)換方程： 橫向位移： 這里 (i)C1和(i)C2是兩個(gè)一體化常數(shù)，可以通過(guò)滿足兼容性來(lái)確定。 將方程（2）和（3）代入（1），可以得到有限元位移和回轉(zhuǎn)曲率，如下所示： 這里下標(biāo)（C），（R）和（D）分別代表曲率，自轉(zhuǎn)和位移。運(yùn)用變分原理，可以得到這些方程[11-13]。 表1 分別比較以位移和應(yīng)力為基礎(chǔ)的有限元方法的歐拉-伯努利梁元素 以位移為基礎(chǔ)的有限元方法 以應(yīng)力為基礎(chǔ)的有限元方法 近似橫向位移自由度 立方米 立方米 近似彎曲應(yīng)力 線性 線性 交點(diǎn)變量 兩端位移和回轉(zhuǎn) 兩端曲率 邊界應(yīng)力滿足條件 位移，回轉(zhuǎn) 位移，回轉(zhuǎn)，彎曲應(yīng)力 自由度數(shù)量 四 二 3.以位移和應(yīng)力為基礎(chǔ)的有限元方法的比較 主要區(qū)別在于以位移為基礎(chǔ)的有限元方法的應(yīng)力場(chǎng)存在不連續(xù)的內(nèi)部因素，同時(shí)具有低階形函數(shù)。主要是因?yàn)椴贿B續(xù)量的產(chǎn)生以及間離散分布。再者，它可能由于使用過(guò)多交點(diǎn)變量而產(chǎn)生剛度矩陣。 以應(yīng)力為基礎(chǔ)的方法與以位移為基礎(chǔ)的方法比較具有很多優(yōu)點(diǎn)。首先，以應(yīng)力為基礎(chǔ)的方法產(chǎn)生的交點(diǎn)變量較少（如表1）。第二，使用以應(yīng)力為基礎(chǔ)的方法時(shí)，彎曲應(yīng)力的邊界條件可以得到滿足。最后，應(yīng)力由體系方程直接計(jì)算得到。 4.方程推導(dǎo) 曲柄滑塊機(jī)構(gòu)如圖1所示，由做剛體運(yùn)動(dòng)的曲柄來(lái)運(yùn)作，該方程由有限元公式推導(dǎo)而得。有限元方程的推導(dǎo)過(guò)程如下：（1）建立剛體運(yùn)動(dòng)學(xué)曲柄滑塊機(jī)構(gòu)；（2）構(gòu)建基于剛體運(yùn)動(dòng)學(xué)機(jī)構(gòu)的翻轉(zhuǎn)梁?jiǎn)卧唬?）確定一套變量用來(lái)描述靈活曲柄滑塊機(jī)構(gòu)的運(yùn)動(dòng)；（4）裝配所有梁?jiǎn)卧?。最后，就可以得到有限元方程，同時(shí)該靈活曲柄滑塊機(jī)構(gòu)的時(shí)間響應(yīng)可以通過(guò)時(shí)間一體化確定。 圖1 靈活曲柄滑塊機(jī)構(gòu) A．翻轉(zhuǎn)梁的元方程 考慮靈活的梁?jiǎn)卧艿絼傮w翻轉(zhuǎn)和回轉(zhuǎn)運(yùn)動(dòng)。疊加在剛體運(yùn)動(dòng)軌跡時(shí)，縱向和橫向方向上允許一些撓度變量。通過(guò)拉格-朗日方程可以得到任意靈活翻轉(zhuǎn)的組件的微分方程。由于彈性變形認(rèn)為是很小的，而且自由度是有限的，這個(gè)方程是線性的并且很容易畫(huà)出來(lái)。推導(dǎo)公式的元素也被很明確的列出來(lái)[8-10]，并且做了簡(jiǎn)要的介紹。 鑒于在軸向有很強(qiáng)的剛度，因此很有必要在縱向方向上合理考慮為剛性梁。所以，縱向方向如一下所示： （5） 這里u1是交點(diǎn)變量，是關(guān)于x軸方向的常數(shù)，如圖2所示。橫向可以表示為： 翻轉(zhuǎn)梁?jiǎn)卧先我恻c(diǎn)的速度可以表示如下： 這里（（i）Vax（i）Vay）是梁?jiǎn)卧贠點(diǎn)的絕對(duì)速度，如圖2所示； 是梁?jiǎn)卧慕撬俣?；（（i）u（i）v）分別是梁?jiǎn)卧先我恻c(diǎn)縱向和橫向的位移，x是梁?jiǎn)卧v向的定位，如圖2所示。 圖2 旋轉(zhuǎn)梁 如果我們把 當(dāng)作組件材料的單位體積；A是組件的橫截面積，L是組件的長(zhǎng)度，組件的動(dòng)能可以表示如下： 均勻剛性組件的軸向彎曲應(yīng)變能量與楊氏模量E有關(guān)，得到二階矩陣I，如下所示： 由縱向拉伸負(fù)荷工作，（i）P,組件的橫向撓度表示如下： 運(yùn)功機(jī)制的縱向負(fù)荷不是一成不變的，與位置和時(shí)間有關(guān)。在忽略縱向彈性形變的前提下，縱向負(fù)荷可能來(lái)自于剛性慣性力，可以表示如下： 這里PR是元件右側(cè)的外部縱向負(fù)載， 是x軸方向上O點(diǎn)的絕對(duì)加速度。如圖2所示。 拉格-朗日形式表示如下： 將公式（5-100）代入（12），并且運(yùn)用歐拉-拉格朗日方程，旋轉(zhuǎn)梁的運(yùn)動(dòng)方程可以表示為一下形式： 這里[Me]、[Ce]和[Ke]分別是元件的質(zhì)量、等效阻尼和等效剛度矩陣；{Fe}是元件的載荷向量。當(dāng)建立質(zhì)量耦合矩陣時(shí)，應(yīng)主要考慮滑塊機(jī)構(gòu)。 B.曲柄滑塊機(jī)構(gòu)方程 提出解決曲柄滑塊機(jī)構(gòu)問(wèn)題的方法，變量是曲率的節(jié)點(diǎn)。裝配所有元件時(shí)，考慮機(jī)構(gòu)的邊界條件是很有必要的。因?yàn)樵搫?dòng)力適用于基礎(chǔ)曲柄結(jié)構(gòu)，在O點(diǎn)存在彎矩，如圖1所示，在O點(diǎn)也存在曲率。如圖1所示的A點(diǎn)和B點(diǎn)，我們假定它們是很小的點(diǎn)。然而，實(shí)際上，彎矩和曲率在這兩個(gè)點(diǎn)上都為零。 因?yàn)楣剑?3）是變量的矩陣表示方式{?} ，這個(gè)公式可以通過(guò)總結(jié)所有的方程來(lái)得到，可以表示如下： 這里[M]、[C]、[K]分別是質(zhì)量、阻尼和剛度矩陣，{F}是負(fù)載向量。 5.穩(wěn)定狀態(tài)基礎(chǔ)上的數(shù)值模擬 曲柄的轉(zhuǎn)速是150rad/s (1432rpm),該靈活曲柄滑塊機(jī)構(gòu)的各項(xiàng)數(shù)值表示如下：R2=0.15(m)，R3=0.30(m), =0.225(kg/m), EI=12.72(N-m2), mB=0.03375(kg)。 這里R2 和R3分別是曲柄和耦合器的長(zhǎng)度，mB是滑塊的質(zhì)量。 通過(guò)曲柄和耦合器的一個(gè)運(yùn)動(dòng)周期，可以看出穩(wěn)態(tài)橫向位移和中點(diǎn)彎曲應(yīng)力的變化情況，以及分析本課題的結(jié)果。可以通過(guò)增加物理阻尼矩陣提高穩(wěn)定性，被稱(chēng)作瑞利阻尼： 這里α和β是兩個(gè)常數(shù)，可以從[15]中對(duì)應(yīng)于兩個(gè)不同頻率的振動(dòng)的阻尼比得到。本文中α和β的值取決于自然頻率。 通過(guò)在運(yùn)動(dòng)方程中增加物理阻尼，也可以通過(guò)Newmark時(shí)間步驟觀測(cè)超過(guò)20個(gè)周期的運(yùn)動(dòng)，從而得到分析結(jié)果。當(dāng)采用數(shù)值時(shí)間積分是出示條件從零開(kāi)始。 誤差可以表示為： 這里QFEQRef 和分別表示以有限元方法和參考方法為基礎(chǔ)的兩個(gè)值，總的來(lái)說(shuō)，可以建立時(shí)間方程，而且很容易被接受，比如能量、位移、彎矩等等。t1 和t2指的是時(shí)間積分的間隔，通常指的是穩(wěn)態(tài)條件下的以個(gè)周期。因?yàn)闆](méi)有一個(gè)合適的準(zhǔn)確的方法，在本文中可以通過(guò)一個(gè)五次多項(xiàng)式表示20個(gè)元件鏈接為基礎(chǔ)的位移有限元方法得到參考值。 Fig. 3. Time responses of the total energy, dimensionless midpoint deflection of the coupler, and the midpoint strain of the coupler at the steady state condition 圖3 總能量的時(shí)間響應(yīng)，耦合器的量綱中點(diǎn)撓度，耦合器在穩(wěn)態(tài)條件下的中點(diǎn)應(yīng)變。 6.數(shù)值模擬 在這一節(jié)中，我們討論剛性曲柄機(jī)構(gòu)。耦合器是唯一的一個(gè)靈活的連桿。在第六節(jié)中以以梁?jiǎn)卧獮榛A(chǔ)，該梁?jiǎn)卧梢宰鰟傂暂S運(yùn)動(dòng)，但是存在橫向撓度。 在第三節(jié)中討論以有限元為基礎(chǔ)的方法時(shí)，很有必要考慮模型的邊界條件和形函數(shù)的相近程度，我們粗略的建立了耦合器應(yīng)變線性分布方程，而且在彎矩不為零的條件下考慮耦合器的邊界條件。 在下面這個(gè)例子中，我們認(rèn)為耦合器是由兩個(gè)、三個(gè)、四個(gè)或者五個(gè)元件組成的，同時(shí)它的曲率分布可以表示為線性方程： 于是，時(shí)間響應(yīng)和總能量誤差，耦合器的中點(diǎn)撓度和應(yīng)變都可以通過(guò)以應(yīng)力為基礎(chǔ)的有限元方法得到。同時(shí)，也評(píng)估了第一自然頻率。 曲柄的轉(zhuǎn)速為150rad/s (1432rpm) ，該靈活的曲柄滑塊機(jī)構(gòu)中各個(gè)部件的值可以表示如下[16]： R2=0.15(m)，R3=0.30(m), =0.225(kg/m), EI=12.72(N-m2), mB=0.03375(kg)。 這里R2 和R3分別是曲柄和耦合器的長(zhǎng)度，mB是滑塊的質(zhì)量。 為了通過(guò)以位移為基礎(chǔ)的有限元方法比較誤差，我們同樣要用它建立一個(gè)機(jī)構(gòu)，結(jié)果可以參考文獻(xiàn)[17]。 表2 兩種有限元方法的第一自然頻率誤差 元件數(shù)目 第一自然頻率 以位移為基礎(chǔ)的有限元方法 以應(yīng)力為基礎(chǔ)的有限元方法 1 1.10E-1(DOF=2) 2 3.91E-3(DOF=4) 7.21E-3(DOF=1) 3 8.10E-4(DOF=6) 1.12E-3(DOF=2) 4 2.60E-4(DOF=8) 3.05E-4(DOF=3) 5 1.07E-4(DOF=10) 1.19E-4(DOF=4) DOF：自由度數(shù)目 表3 兩種有限元方法的總能量誤差 元件數(shù)目 第一自然頻率 以位移為基礎(chǔ)的有限元方法 以應(yīng)力為基礎(chǔ)的有限元方法 1 1.94E-2(DOF=2) 2 3.21E-3(DOF=4) 8.85E-4(DOF=1) 3 1.92E-3(DOF=6) 3.77E-4(DOF=2) 4 1.20E-3(DOF=8) 2.82E-4(DOF=3) 5 9.00E-4(DOF=10) 2.17E-4(DOF=4) DOF：自由度數(shù)目 圖3顯示了總能量的時(shí)間響應(yīng)，耦合器的量綱中點(diǎn)撓度，耦合器在穩(wěn)態(tài)條件下的中點(diǎn)應(yīng)變。表2-5分別比較了以位移為基礎(chǔ)和以應(yīng)力為基礎(chǔ)的有限元方法的第一自然頻率誤差、總能量、耦合器的中點(diǎn)撓度量綱、以及耦合器的中點(diǎn)應(yīng)變。誤差可以由公式（16）得到。結(jié)果表明，當(dāng)兩種方法中的元件數(shù)目相同時(shí)，以應(yīng)力為基礎(chǔ)的方法誤差較以位移為基礎(chǔ)的誤差大。但是，當(dāng)自由度的數(shù)目相同時(shí)，以應(yīng)力為基礎(chǔ)的有限元方法的誤差比以位移為基礎(chǔ)的有限元方法的誤差小很多。同時(shí)，我們注意到當(dāng)元件相同，除去第一自然頻率誤差時(shí)，以應(yīng)力為基礎(chǔ)的有限元方法的誤差也比以位移為基礎(chǔ)的有限元方法的小很多。這說(shuō)明以應(yīng)力為基礎(chǔ)的有限元方法可以提供大量精確的解決動(dòng)態(tài)彈塑性問(wèn)題的方法。 表4 兩種有限元方法的耦合器中點(diǎn)撓度誤差 元件數(shù)目 第一自然頻率 以位移為基礎(chǔ)的有限元方法 以應(yīng)力為基礎(chǔ)的有限元方法 1 3.60E-1(DOF=2) 2 5.27E-2(DOF=4) 1.29E-2(DOF=1) 3 3.26E-2(DOF=6) 8.41E-3(DOF=2) 4 2.14E-2(DOF=8) 6.12E-3(DOF=3) 5 1.57E-2(DOF=10) 4.90E-3(DOF=4) DOF：自由度數(shù)目 表5 兩種有限元方法的耦合器中點(diǎn)應(yīng)變誤差 元件數(shù)目 第一自然頻率 以位移為基礎(chǔ)的有限元方法 以應(yīng)力為基礎(chǔ)的有限元方法 1 4.25E-1(DOF=2) 2 1.65E-1(DOF=4) 2.38E-1(DOF=1) 3 5.35E-2(DOF=6) 2.89E-2(DOF=2) 4 4.38E-2(DOF=8) 2.28E-2(DOF=3) 5 2.27E-2(DOF=10) 1.80E-2(DOF=4) DOF：自由度數(shù)目 7.結(jié)論 本文提出了一種新的以應(yīng)力為基礎(chǔ)的有限元方法來(lái)解決歐拉-拉格朗日梁?jiǎn)栴}。該方法尤其適用于解決動(dòng)態(tài)彈塑性問(wèn)題。并且提出了梁的近似曲率。然后我們可以通過(guò)整合近似曲率得到橫向撓度和應(yīng)力分布。在整合過(guò)程中，有必要使梁?jiǎn)卧倪吔鐥l件得到滿足，從而可以得到整合常數(shù)。本文中，我們提出了在高速運(yùn)作下解決靈活曲柄滑塊機(jī)構(gòu)問(wèn)題。結(jié)果表明，在同樣的自由度下，以應(yīng)力為基礎(chǔ)的有限元方法的誤差小于常規(guī)方法的誤差，常規(guī)方法也就是以位移為基礎(chǔ)的有限元方法。同樣，在元件數(shù)目相同的條件下，以應(yīng)力為基礎(chǔ)的有限元方法可以提供更多準(zhǔn)確的解決方法。 參考文獻(xiàn) [1]B. Fraeijs de Veubeke, “Displacement an equilibrium models in the finite element method”, stress Analysis, edited by O.C. Zienkiewicz, Wiley, New York, 1965. [2]B. Fraeijs de Veubekd and O.C. Zienkiewicz, “Strai-energy bounds in finite-element analysis by slab analogy” J. Strain Analysis, Vol. 2, pp. 265-271, 1967. [3]Z. Wieckowski, S.K. Youn, and B.S. Moon, “Stressedbased finite element analysis of plane plasticity problems”, Int. J Numer. Meth. Engng., Vol. 44, pp. 1505-1525, 1999. [4]H. Chanda and K.K. Tamma, “Developments encompassig stress based finite element formulations for materially nonlinear static dynamic problems”, Comp. Struct.,Vol. 59, No. 3, pp. 583-592, 1996. [5]M. Kaminski, “Stochastic second-order perturbation approach to the stress-based finite element method”, Int. J. Solids and Struct., Vol. 38, No. 21, pp. 3831-3852, 2001. [6]O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, McGraw-Hill, London, 2000. [7]R.H. Gallagher, Finite Element Fundamentals, Prentice-Hall, Englewood Cliffs, 1975. [8]W.L. Cleghorn, 1980, Analysis and design of high-speed flexible mechanism, Ph. D. Thesis, University of Toronto. [9]W.L. Cleghorn, R. G. Fenton, and B. Tabarrok, 1981, “Finite element analysis of high-speed flexible mechanisms”, Mechanism and Machine Theory, 16(4), 407-424. [10]W.L. Cleghorn, R.G. Fenton, and B. Tabarrok, 1984, “Steady-state vibrational response of high-speed fexible mechanisms”, Mechanism and Machine Theory, 19(4/5)417-423. [11]Y.L. Kuo, W.L. Cleghorn and K. Behdinan, “Stress-bsed Finite Element Method for Euler-Bernoulli Beams”,Transactions of the Canadian Society for Mechanical Engineering, Vol. 30(1), pp. 1-6, 2006. [12]Y.L. Kuo, W.L. Cleghorn, and K. Behdinan “Applicatons of Stress-based Finite Element Method on Euler-Bernoulli Beams ”, Proceedings of the 20th Canadian Congress of Applied Mechanics, Montreal, Quebec, Canada, May 30-Jun 2, 2005. [13]Y.L. Kuo, Applications of the h-, p-, and r-refinements of the Finite Element Method on Elasto-dynamic Problems, Ph.D. Thesis, University of Toronto, 2005. [14]L. Meirovitch, 1967, Analytical Methods in Vibrations, Macmillan, New York, 436-463. [15]K.J. Bathe, 1996, Finite Element Procedures, Prentice Hall, Englewood Cliffs, NJ, USA. [16]A.L. Schwab and J.P. Meijaard, 2002, “Small vibratons superimposed on prescribed rigid body motion”, Mulibody System Dynamics, 8, 29-49. [17]Y.L. Kuo and W.L. Cleghorn, “The h-p-r-refinement Finite Element Analysis of a High-speed Flexible Slider Crank Mechanism”, Journal of Sound and Vibration, in press. 英文原稿 12thIFToMM World Congress,Besancon,June 18-21,2007 Application of Stress-based Finite Element Methodto a Flexible Slider Crank Mechanism Y.L.Kuo? W.L.Cleghorn University of Toronto University of Toronto Toronto,Canada Toronto,Canada Abstract—This paper presents a new procedure to apply the stress-based finite element method on Euler-Bernoulli beams.An approximated bending stress distribution is selected,and then the approximated transverse displacement is determined by integration.The proposed approach is applied to solve a flexible slider crank mechanism.The formulation is based on the Euler-Lagrange equation,for which the Lagrangian includes the components related to the kinetic energy,the strain energy,and the work done by axial loads in a link that undergoes elastic transverse deflection.A beam element is modeled based on a translating and rotating motion.The results demonstrate the error comparison obtained from the stress-and displacement-based finite element methods. Keywords:stress-based finite element method;slider crank mechanism;Euler-Lagrange equation. 1.Introduction The displacement-based finite element method employs complementary energy by imposing assumed displacements.This method may yield the discontinuities of stress fields on the inter-element boundary while employing low-order elements,and the boundary conditions associated with stress could not be satisfied.Hence,an alternative approach was developed and called the stress-based finite element method,which utilizes assumed stress functions.Veubeke and Zienkiewicz[1,2]were the first researchers introducing the stress-based finite element method.After that,the method was applied to a wide range of problems and its applications[3-5]In addition,there are various books providing details about the method[6,7]. The operation of high-speed mechanisms introduces vibration,acoustic radiation,wearing of joints,and inaccurate positioning due to deflections of elastic links.Thus,it is necessary to perform an analysis of flexible elasto-dynamics of this class of problems rather than the analysis of rigid body dynamics.Flexible mechanisms are continuous dynamic systems with an infinite number of degrees of freedom,and their governing equations of motion are modeled bynonlinear partial differential equations,but their analytical______________________ ?Email:ylkuo@mie.utoronto.ca solutions are impossible to obtain.Cleghorn et al.[8-10]included the effect of axial loads on transverse vibrations of a flexible four-bar mechanism.Also,they constructed a translating and rotating beam element with a quintic polynomial,which can effectively predict the transverse vibration and the bending stress. This paper presents a new approach for the implementation of the stress-based finite element method on the Euler-Bernoulli beams.The developed approach first selects an assumed stress function.Then,the approximated transverse displacement function is obtained by integrating the assumed stress function.Thus,this approach can satisfy the stress boundary conditions without imposing a constraint.We apply this approach to solve a flexible slider crank mechanism.In order to show the accuracy enhancement by this approach,the mechanism is also solved by the displace-based finite element method.The results demonstrate the error comparison. II.Stress-based Method for Euler-Bernoulli Beams The bending stress of Euler-Bernoulli beams is associated with the second derivative of the transverse displacement,namely curvature,which can be approximated as the product of shape functions and nodal variables: Where is a row vector of shape functions for the ith element; is a column vector of nodal curvatures,y is the lateral position with respect to the neutral line of the beam,E is the Young’s modulus,and is the transverse displacement,which is a function of axial position x. Integrating Eq.(1)leads to the expressions of the rotation and the transverse displacement as Rotation： Transverse displacement: Where and are two integration constants for the ith element,which can be determined by satisfying the compatibility. Substituting Eqs.(2)and(3)into(1),the finite element displacement,rotation and curvature can be expressed as: where the subscripts(C),(R)and(D)refer to curvature,rotation and displacement,respectively.By applying the variational principle,the element and global equations can be obtained[11-13]. Table 1:Comparison of the displacement-and the stress-based finite element methods for an Euler-Bernoulli beam element III.Comparisons of the Displacement-and Stress-based Finite Element Methods The major disadvantage of the displacement-based finite element method is that the stress fields at the inter-element nodes are discontinuous while employing low-degree shape functions.This discontinuity yields one of the major concerns behind the discretization errors.In addition,it might use excessive nodal variables while formulating stiffness matrices. The stress-based method has several advantages over the displacement-based finite method.First of all,the stress-based method produces fewer nodal variables (Table 1).Secondly,when employing the stress-basedfinite method,the boundary conditions of bending stress can be satisfied,and the stress is continuous at theinter-element nodes.Finally,the stress is calculated directly from the solution of the global system equations.However,the only disadvantage of the stress-based finite method is that the integration constants are different for each element. IV.Generation of Governing Equation The slider crank mechanism shown in Fig.1 is operated with a prescribed rigid body motion of the crank,and the governing equations are derived using a finite element formulation.The derivation procedure of the finite element equations involves:(1)deriving the kinematics of a rigid body slider crank mechanism;(2) constructing a translating and rotating beam element based on the rigid body motion of the mechanism;(3)defining a set of global variables to describe the motion of a flexible slider crank mechanism;(4)assembling all beam elements.Finally,the global finite element equations can be obtained,and the time response of a flexible slider crank mechanism can be obtained by time integration. A.Element equation of a translating and rotating beam Consider a flexible beam element subjected to prescribed rigid body translations and rotations.Superimposed on the rigid body trajectory,a finite number of deflection variables in the longitudinal and transverse directions is allowed.The Euler-Lagrange equation is used to derive the governing differential equations for an arbitrarily translating and rotating flexible member.Since elastic deflections are considered small,and there is a finite number of degrees of freedom,the governing equations are linear and are conveniently written in matrix form.The derivation of the element equations has been precisely presented in [8-10],and this section provides a brief summary. In view of high axial stiffness of a beam,it is reasonable to consider the beam as being rigid in its longitudinal direction.Hence,the longitudinal deflection is given as where u1 is a nodal variable,which is constant with respect to the x direction shown in Fig.2.The transverse deflection can be represented as The velocity of an arbitrary point on the beam element with a translating and rotating motion is given as where is the absolute velocity of point O of the beam element shown in Fig.2;θ?is the angular velocity of the beam element; are the longitudinal and transverse displacements of an arbitrary point on the beam element,respectively;x is a longitudinal position on the beam element shown in Fig. 2. If we letρbe the mass per unit volume of element material;A,the element cross-sectional area,and L the element length,then the kinetic energy of an element is expressed as The flexural strain energy of uniform axially rigid element with the Young’s modulus,E,and second moment of area,I,is given as The work done by a tensile longitudinal load,(i)P,in an element that undergoes an elastic transverse deflection is given by[14] Longitudinal loads in a moving mechanism element are not constant,and depend both on the position in the element and on time.With the longitudinal elastic motions neglected,the longitudinal loads may be derived from the rigid body inertia forces,and can be expressed as where PR is an external longitudinal load acting at theright hand end of an element,andox (i )ais the absolute eacceleration of the point O in the x direction shown in Fig.2. The Lagrangian takes the form Substituting Eqs.(5-10)into(12),and employing the Euler-Lagrange equations,the governing equations of motion for a rotating and translating elastic beam can be expressed in the following matrix form: where[Me],[Ce]and[Ke]are mass,equivalent damping,and equivalent stiffness matrices of a element,respectively;{Fe}is a load vector of an element.When formulating the mass matrix of the coupler,the mass of the slider should be taken into account. B.Global equations of slider crank mechanism For the proposed approach to solve a flexible slider crank mechanism,the global variables are the curvatures on the nodes.For assembling all elements,it is necessary to consider the boundary conditions applied to the mechanism.Since a prescribed motion applied to the base of the crank,there is a bending moment at point O shown in Fig.1,i.e.,the curvature at point O exists.For points A and B shown in Fig.1,we presume that both points refer to pin joints.Thus,the bendingmoments and the curvatures at both points are zeros. Since Eq.(13)is a matrix-form expression in terms of the vector of global variables{φ},the global equations can be obtained by directly summing up all of element equations,which can be expressed as where[M],[C],[K]are global mass,damping and stiffness matrices,respectively;{F}is a global load vector. V.Numerical simulation based on steady state The rotating speed of the crank is operating at 150rad/s(1432 rpm),and the system parameters of a flexible slider crank are as follows: R2=0.15(m),R3=0.30(m),ρA=0.225(kg/m),EI=12.72(N-m2),mB=0.03375(kg) where R2 and R3 are the lengths of the crank and coupler,respectively;mB is the mass of the slider. The analytical results of this paper are presented by plotting steady state transverse displacements and bending strains of midpoints on crank and coupler throughout a cycle of motion.The steady state can be obtained by adding a physical damping matrix,namely Rayleigh damping whereαandβare two constants,which can be determined from two given damping ratio that correspond to two unequal frequencies of vibration[15]. In this paper,the values ofαandβare determined based on the first two natural frequencies. By adding physical damping to the equations of motion,the analytical solution is obtained by performing the constant time-step Newmark method over twenty cycles of motion.The initial conditions are set to zeros when performing numerical time integration. The error indicator is defined as where QFE and QRef are two quantities based on a finite element solution and a reference solution,respectively.Generally,they are functions of time,and they can be arbitrarily selected,such as energy,displacement,bending strain,etc.t1 and t2 refer to the interval of timeintegration,which are usually one cycle after steady-state condition has been reached.Since an exact solution is not available,a reference solution is obtained by the displace