胡宇達(dá), 秦曉北
(1. 燕山大學(xué) 建筑工程與力學(xué)學(xué)院,河北 秦皇島 066004;2. 燕山大學(xué) 河北省重型裝備與大型結(jié)構(gòu)力學(xué)可靠性重點(diǎn)實(shí)驗(yàn)室,河北 秦皇島 066004)
導(dǎo)電導(dǎo)磁結(jié)構(gòu)件廣泛應(yīng)用于航空航天、超導(dǎo)發(fā)電機(jī)、電磁傳感器、發(fā)動(dòng)機(jī)等眾多高新科技裝置與設(shè)備中,處于電磁場中的導(dǎo)電彈性體因?yàn)槭艿诫?、磁、力等多種因素的作用,而表現(xiàn)出復(fù)雜的耦合動(dòng)力學(xué)行為,或?qū)⒅苯佑绊懴到y(tǒng)運(yùn)動(dòng)的安全性與可靠性。近年來,許多學(xué)者展開了這一領(lǐng)域的研究。Moon 等[1]較早的對薄板的磁彈性屈曲問題進(jìn)行了研究,提出了Moon-Pao模型,為后面研究的學(xué)者提供了很多的理論參考;鄭曉靜等[2-3]對磁場中鐵磁梁式板、軟鐵磁殼的動(dòng)力學(xué)特性進(jìn)行了分析;Hasanyan等[4-5]研究了載流平板的磁彈性屈曲和后屈曲問題,并對處于傾斜磁場中的導(dǎo)電板的非線性振動(dòng)問題進(jìn)行了數(shù)學(xué)建模與研究;Takagi等[6]在考慮磁黏滯阻尼的影響下,采用了新的數(shù)值分析方法,對處于磁場環(huán)境中的薄板的動(dòng)力學(xué)行為進(jìn)行了研究;胡宇達(dá)等[7-10]研究了磁場中軸向運(yùn)動(dòng)導(dǎo)電薄板的非線性動(dòng)力學(xué)行為以及旋轉(zhuǎn)圓板的磁彈性主共振問題;Bayat等[11]對變厚度功能梯度環(huán)形旋轉(zhuǎn)板的磁熱機(jī)械響應(yīng)進(jìn)行了研究;Zenkour[12]研究了功能梯度環(huán)形夾層板的磁熱彈性響應(yīng);高原文等[13]研究了橫向磁場中鐵磁梁式板的混沌運(yùn)動(dòng)問題;王省哲等[14]針對磁場中鐵磁梁式板的磁彈性屈曲問題進(jìn)行了研究;Bhandari等[15]研究了軸向非均勻磁場中旋轉(zhuǎn)盤引起的鐵磁流體的流動(dòng)問題;Rad等[16]針對均布外激勵(lì)作用下變厚度環(huán)形FGM板的磁熱動(dòng)力響應(yīng)進(jìn)行了分析。
在旋轉(zhuǎn)結(jié)構(gòu)動(dòng)力學(xué)問題的研究中,Pei等[17]針對臨界和超臨界速度下的旋轉(zhuǎn)柔性薄板進(jìn)行了研究,分析了溫度分布、熱應(yīng)力等對旋轉(zhuǎn)盤動(dòng)態(tài)穩(wěn)定性及穩(wěn)態(tài)振幅的影響。Sharma等[18]研究了變厚度及變密度彈性旋轉(zhuǎn)環(huán)形薄板的有限差分解。Hussain等[19]研究了加速旋轉(zhuǎn)盤在黏性流體中的數(shù)值解。Jabbari等[20]針對變厚度FGM旋轉(zhuǎn)圓板進(jìn)行研究,分析了厚度、角速度等參數(shù)對旋轉(zhuǎn)盤力學(xué)性能的影響。
本文針對磁場環(huán)境中兩頻激勵(lì)作用下的旋轉(zhuǎn)運(yùn)動(dòng)導(dǎo)電圓板的諧波-組合共振問題進(jìn)行研究,得到系統(tǒng)穩(wěn)態(tài)解穩(wěn)定的判定條件及幅頻響應(yīng)方程,繪制幅頻特性曲線、動(dòng)相平面軌跡圖并分析系統(tǒng)不同物理參量對諧波-組合共振特性的影響。
圖1為旋轉(zhuǎn)運(yùn)動(dòng)導(dǎo)電圓板在外加磁場B0作用下的示意圖。圓板的板厚為h,半徑為R,密度為ρ。
圖1 旋轉(zhuǎn)圓板示意圖Fig.1 Rotating circular plate diagram
設(shè)旋轉(zhuǎn)運(yùn)動(dòng)圓板內(nèi)任意一點(diǎn)的位移為
ur1=ur+zu1
(1)
uθ1=uθ+zv1
(2)
uz=w
(3)
用位移表達(dá)的中面內(nèi)力
(4)
(5)
(6)
則旋轉(zhuǎn)圓板的中面應(yīng)變勢能表示為
(7)
用位移表示的彎矩,扭矩
(8)
(9)
(10)
則旋轉(zhuǎn)圓板的彎曲應(yīng)變勢能表示為
(11)
板上任意一點(diǎn)的速度表達(dá)式
(12)
式中:i,j,k分別為坐標(biāo)軸r,θ,z方向上的單位向量。
則旋轉(zhuǎn)圓板的動(dòng)能為
(13)
旋轉(zhuǎn)運(yùn)動(dòng)導(dǎo)電圓板電流密度表達(dá)式
J=σ0(E+V×B)
(14)
式中:σ0為電導(dǎo)率;E為電場強(qiáng)度矢量;B為磁感應(yīng)強(qiáng)度矢量。
磁場環(huán)境中旋轉(zhuǎn)運(yùn)動(dòng)圓板受到的洛倫茲力矢量表達(dá)式
f(fr,fθ,fz)=J×B0
(15)
(16)
(17)
(18)
則電磁力、電磁力矩在虛位移上所做的虛功為
(19)
假設(shè)圓板在橫向位移w上產(chǎn)生了微小的變化量δw,則Fz在虛位移δw上所做的功為
(20)
應(yīng)用Hamilton原理建立旋轉(zhuǎn)運(yùn)動(dòng)圓板的振動(dòng)微分方程
(21)
將式(7)、式(11)、式(13)、式(19)和式(20)代入式(21)中,得到旋轉(zhuǎn)運(yùn)動(dòng)導(dǎo)電圓板軸對稱非線性磁彈性振動(dòng)微分方程
(22)
受組合動(dòng)載荷Fz=F1cosω1t+F2cosω2t作用下做旋轉(zhuǎn)運(yùn)動(dòng)的圓形薄板,設(shè)滿足周邊夾支邊界條件的位移函數(shù)為
(23)
將式(23)代入式(22),應(yīng)用伽遼金法進(jìn)行積分,推得旋轉(zhuǎn)運(yùn)動(dòng)圓板無量綱化振動(dòng)方程
(24)
對于兩頻激勵(lì),考慮超諧波共振和組合共振同時(shí)發(fā)生的聯(lián)合共振情況,系統(tǒng)的無量綱化頻率為
(25)
式中:σ1,σ2為頻率調(diào)諧參數(shù)。
(26)
采用多尺度法求解方程,取時(shí)間尺度T0=τ,T1=ετ,討論一次近似解,可將式(26)的解表示為
q(τ,ε)=q0(T0,T1)+εq1(T0,T1)
(27)
將式(27)代入式(26),令方程兩端的ε0和ε1項(xiàng)的系數(shù)相等,可得到如下線性偏微分方程分別為
(28)
(29)
式(28)的通解為
q0=A(T1)eiT0+Λ1eiΩ1T0+Λ2eiΩ2T0+cc
(30)
將式(30)的通解代入式(29),根據(jù)消除久期項(xiàng)的條件,得到關(guān)于A的關(guān)系式
(31)
(32)
(33)
分析式(32)和式(33),當(dāng)且僅當(dāng)σ1T1-β和σ2T1-β都為常數(shù)時(shí)才能存在穩(wěn)態(tài)運(yùn)動(dòng),即β′=σ1=σ2=σ,同時(shí),令γ=σT1-β,對于穩(wěn)態(tài)運(yùn)動(dòng)a′=γ′=0,則由式(32)和式(33)可得到磁場環(huán)境下旋轉(zhuǎn)運(yùn)動(dòng)圓板諧波-組合共振幅頻響應(yīng)方程
(34)
研究系統(tǒng)穩(wěn)態(tài)運(yùn)動(dòng)下解的穩(wěn)定性問題,設(shè)
a=a0+a1,γ=γ0+γ1
(35)
式中:a0,γ0為系統(tǒng)穩(wěn)態(tài)運(yùn)動(dòng)下的穩(wěn)態(tài)解;a1,γ1為小的擾動(dòng)量。
將穩(wěn)態(tài)解a,γ代入式(32)和式(33),并對小的擾動(dòng)量a1,γ1進(jìn)行Taylor展開,可得到關(guān)于a1和γ1的一階線性近似式
(36)
(37)
根據(jù)Lyapunov穩(wěn)定性近似理論,穩(wěn)態(tài)運(yùn)動(dòng)下解的穩(wěn)定性依賴于式(36)和式(37)右端系數(shù)矩陣的特征值。通過式(36)和式(37)的Jacobi矩陣并考慮到穩(wěn)態(tài)運(yùn)動(dòng)下a′=γ′=0,可得到特征方程
(38)
再根據(jù)Hurwitz判據(jù),可得系統(tǒng)穩(wěn)態(tài)解穩(wěn)定的充要條件
(39)
針對磁場環(huán)境下受橫向動(dòng)載荷作用的鋁制旋轉(zhuǎn)圓形薄板進(jìn)行數(shù)值分析,得到反映旋轉(zhuǎn)運(yùn)動(dòng)圓板諧波-組合共振特性的圖形。主要參數(shù):密度ρ=2 670 kg/m3,泊松比μ=0.34,彈性模量E=71 GPa,電導(dǎo)率σ0=3.63×107(s/m),圓板半徑R=0.4 m。
圖2為轉(zhuǎn)速Ω=8 000 r/min時(shí)共振振幅a隨頻率調(diào)諧參數(shù)εσ的變化規(guī)律曲線圖。由圖2可知,在給定εσ范圍內(nèi),共振幅頻響應(yīng)曲線向右偏移,呈現(xiàn)硬彈簧特性,隨著εσ的改變,系統(tǒng)出現(xiàn)多值性和跳躍現(xiàn)象,呈現(xiàn)典型的非線性振動(dòng)特性;隨著磁感應(yīng)強(qiáng)度的增加,共振曲線呈現(xiàn)明顯的內(nèi)縮趨勢,共振區(qū)域變窄,振幅減小,可見磁場起到了電磁阻尼的作用;圖2(b)、圖2(c)和圖2(d)中顯示,隨著板厚的減小,外激勵(lì)力幅值F1,F(xiàn)2的增大,在給定的εσ范圍內(nèi),系統(tǒng)共振區(qū)域變寬,激發(fā)多解區(qū)域的臨界點(diǎn)呈滯后右移趨勢。
圖2 幅頻a-εσ響應(yīng)圖Fig.2 Amplitude frequency response
圖3為振幅a隨磁感應(yīng)強(qiáng)度B0z變化曲線圖,取F1=500 N/m2,F(xiàn)2=1 000 N/m2,Ω=8 000 r/min和h=4.5 mm。從圖3可知,在給定的B0z范圍內(nèi),圖中曲線呈現(xiàn)關(guān)于B0z=0左右對稱形式,在給定調(diào)諧參數(shù)εσ范圍內(nèi),曲線呈現(xiàn)從單值到多值并逐漸縮頸最終分離出上部的封閉曲線趨勢;調(diào)諧參數(shù)的改變對曲線變化規(guī)律影響顯著。圖中由虛線分割的多值解區(qū)域中,區(qū)域1、區(qū)域3代表穩(wěn)定解區(qū)域,區(qū)域2代表不穩(wěn)定解區(qū)域。
圖4為不同調(diào)諧參數(shù)下振幅a-力幅F1,F(xiàn)2的三維變化圖,取B0z=0.5 T,h=4.5 mm及Ω=8 000 r/min。從圖4可知,在給出的幾種調(diào)諧參數(shù)下,激勵(lì)力幅值F1=0時(shí),無論F2取何值,共振振幅均為0,而當(dāng)F≠0時(shí),F(xiàn)2=0可能激發(fā)振動(dòng),激發(fā)的共振是超諧波共振。從F1等于某一值的截面看,共振振幅a隨F2的變化曲線可以出現(xiàn)持續(xù)增加、先增加后減小或先減小后增加等多種復(fù)雜情況,調(diào)諧參數(shù)εσ的微小改變對曲線的影響顯著。
圖5為不同磁感應(yīng)強(qiáng)度下振幅a-力幅F1,F(xiàn)2的三維變化圖,取h=4.5 mm,εσ=0.001及Ω=8 000 r/min。從圖5可知,磁感應(yīng)強(qiáng)度取較小值時(shí),振幅a隨激勵(lì)力幅值F1,F(xiàn)2變化曲線較復(fù)雜,振幅有先增加后減小或先減小后增加等多種形式,增大磁感應(yīng)強(qiáng)度到某一值后,振幅呈單一增大趨勢。
圖6為h=4.5 mm時(shí)振幅a隨外激勵(lì)力幅值F1變化規(guī)律曲線圖。從圖6(a)可知,B0z=0.5 T時(shí),小的F1就可以激發(fā)系統(tǒng)的多值解,當(dāng)磁感應(yīng)強(qiáng)度B0z增大到一定值后,系統(tǒng)的解退化為單值,共振振幅隨之減小;從圖6(b)可知,系統(tǒng)的解呈現(xiàn)從一個(gè)到三個(gè)再到一個(gè)的轉(zhuǎn)變,隨著轉(zhuǎn)速的增大,系統(tǒng)的多值解區(qū)域變窄,共振曲線向左偏移;圖6(c)為不同調(diào)諧參數(shù)下的a-εσ變化規(guī)律曲線,從圖6(c)可知,εσ取值的大小決定了單值曲線到多值曲線的變化,隨著εσ的增大,系統(tǒng)的多值解區(qū)域變寬,多值解區(qū)域共振振幅隨之增加;圖6(d)顯示激勵(lì)力幅值對共振曲線影響較大,當(dāng)F2=0時(shí),系統(tǒng)解呈現(xiàn)從單值到多值,后又變?yōu)閱沃档淖兓厔?,此時(shí)系統(tǒng)僅激發(fā)超諧波共振,增大F2,共振曲線繼續(xù)向左偏移,解由多值逐漸退化為單值。
圖7為Ω=8 000 r/min時(shí)振幅a隨激勵(lì)力幅值F2變化規(guī)律曲線圖。圖7(a)顯示隨著F1的增加,共振曲線向左偏移,呈現(xiàn)軟彈簧特性,且隨著磁感應(yīng)強(qiáng)度B0z的增大,多值解區(qū)域減小,共振曲線呈明顯內(nèi)縮趨勢;圖7(b)、圖7(c)顯示隨著板厚、調(diào)諧參數(shù)的增大,共振曲線向右偏移,激發(fā)多值解的臨界點(diǎn)呈滯后右移趨勢;圖7(d)顯示,隨著F1的增大,共振曲線向左偏移,激發(fā)超諧-組合共振多值解所需要的激勵(lì)力F2值減小,F(xiàn)1=0時(shí),a-F2共振曲線不會(huì)被激發(fā)。
圖3 振幅a-磁感應(yīng)強(qiáng)度B0z響應(yīng)圖Fig.3 Amplitude-magnetic density response
圖4 振幅a-力幅F1,F(xiàn)2的三維變化圖Fig.4 Three-dimensional variations of amplitude a-excitation amplitude F1, F2
圖5 振幅a-力幅F1、F2的三維變化圖Fig.5 Three-dimensional variations of amplitude a-excitation amplitude F1,F(xiàn)2
圖6 振幅a-激勵(lì)力幅值F1響應(yīng)圖Fig.6 Amplitude a-excitation amplitude F1
圖8為改變初始條件得到的動(dòng)相平面軌跡圖,箭頭指出了軌跡的運(yùn)動(dòng)方向,取參數(shù)F1=500 N/m2,F(xiàn)2=1 000 N/m2,h=4.5 mm,B0z=0.5 T,Ω=8 000 r/min。圖中顯示當(dāng)調(diào)諧參數(shù)εσ取不同值時(shí),穩(wěn)態(tài)解的個(gè)數(shù)也不同。圖8(a)只有一個(gè)穩(wěn)定焦點(diǎn)S1,其共振振幅為aS=0.157,記為S1(aS=0.157)。圖8(b)中顯示有兩個(gè)穩(wěn)定焦點(diǎn)S1(aS=0.052),S3(aS=0.393)和一個(gè)不穩(wěn)定鞍點(diǎn)S2(aS=0.338)??梢?,圖8(b)中的穩(wěn)態(tài)解與圖8(c)中的結(jié)果是一致的,且圖8(c)中的上支和下支曲線是穩(wěn)定解,中部曲線是不穩(wěn)定解,其解的不同取決于初始條件。
圖9為取板厚h=4.5 mm,εσ=0.08時(shí)的動(dòng)相平面軌跡圖。圖9(a)給出了不同磁感應(yīng)強(qiáng)度下對應(yīng)的兩個(gè)穩(wěn)定焦點(diǎn),其中B0z=2 T對應(yīng)點(diǎn)S1(aS=0.393),B0z=0.5 T對應(yīng)點(diǎn)S2(aS=0.381);圖9(b)給出了不同力幅F2下對應(yīng)的兩個(gè)穩(wěn)定焦點(diǎn),其中F2=800 N/m2和F2=1 000 N/m2分別對應(yīng)點(diǎn)S1(aS=0.452)和S2(aS=0.393)??梢姡瑓⒘康母淖儗ο到y(tǒng)非線性穩(wěn)態(tài)解有顯著影響。
圖7 振幅a-激勵(lì)力幅值F2響應(yīng)圖Fig.7 Amplitude a-excitation amplitude F2
圖8 動(dòng)相平面軌跡Fig.8 Moving phase plane trajectory
圖9 動(dòng)相平面軌跡Fig.9 Moving phase plane trajectory
本文針對磁場環(huán)境中受兩頻激勵(lì)作用旋轉(zhuǎn)運(yùn)動(dòng)圓板的超諧-組合聯(lián)合共振問題進(jìn)行了研究,研究結(jié)果表明:
(1)在給定頻率參數(shù)范圍內(nèi),共振曲線向右偏移,呈硬彈簧特性,并出現(xiàn)多值性和跳躍現(xiàn)象,呈現(xiàn)典型的非線性振動(dòng)特性。
(2)隨著力幅的增加,共振曲線向左偏移,呈現(xiàn)軟彈簧特性,且隨著磁感應(yīng)強(qiáng)度的增大,多值解區(qū)域減小,共振曲線呈明顯內(nèi)縮趨勢。
(3)初始條件對非線性系統(tǒng)的穩(wěn)態(tài)響應(yīng)影響較大,幅頻響應(yīng)曲線上支和下支部分是穩(wěn)定解,中間部分是不穩(wěn)定解。由此可見,通過改變不同參量,可實(shí)現(xiàn)穩(wěn)態(tài)解從多值到單值、穩(wěn)定到不穩(wěn)定的轉(zhuǎn)變,從而達(dá)到控制系統(tǒng)共振現(xiàn)象的目的。
參 考 文 獻(xiàn)
[ 1 ] MOON F C, PAO Y H. Magnetoelastic buckling of a thin plate [J]. ASME Journal of Applied Mechanics, 1968, 35(1): 53-68.
[ 2 ] 鄭曉靜,劉信恩. 鐵磁導(dǎo)電梁式板在橫向均勻磁場中的動(dòng)力特性分析[J]. 固體力學(xué)學(xué)報(bào), 2000,21(3): 243-250.
ZHENG Xiaojing, LIU Xinen. Analysis on dynamic characteristics for ferromagnetic conducting plates in a transverse uniform magnetic field [J]. Acta Mechanica Solida Sinica, 2000, 21(3): 243-250.
[ 3 ] ZHENG X J, WANG X Z. A magnetoelastic theoretical model for soft ferromagnetic shell in magnetic field [J]. International Journal of Solids and Structures, 2003, 40(24): 6897-6912.
[ 4 ] HASANYAN D J, LIBRESCU L, AMBUR D R. Buckling and postbuckling of magnetoelastic flat plates carrying an electric current [J]. International Journal of Solids and Structures, 2006, 43(16): 4971-4996.
[ 5 ] HASANYAN D J, KHACHATURYAN G M, PILIPOSYAN G T. Mathematical modeling and investigation of nonlinear vibration of perfectly conductive plates in an inclined magnetic field [J]. Thin-Walled Structures, 2001, 39(1): 111-123.
[ 6 ] TAKAGI T, TANI J. New numerical analysis method of dynamic behavior of a thin plate under magnetic field considering magnetic viscous damping effect [J]. International Journal of Applied Electromagnetics in Materials, 1993, 4(1): 35-42.
[ 7 ] 胡宇達(dá),胡朋. 軸向運(yùn)動(dòng)導(dǎo)電板磁彈性非線性動(dòng)力學(xué)及分岔特性[J]. 計(jì)算力學(xué)學(xué)報(bào),2014,31(2):180-186.
HU Yuda, HU Peng. Magneto-elastic nonlinear dynamics and bifurcation of axially moving current-conducting plate [J]. Chinese Journal of Computational Mechanics, 2014, 31(2): 180-186.
[ 8 ] HU Y D, LI J. The magneto-elastic subharmonic resonance of current-conducting thin plate in magnetic field [J]. Journal of Sound and Vibration, 2009, 319(3/4/5): 1107-1120.
[ 9 ] HU Y D, WANG T. Nonlinear free vibration of a rotating circular plate under the static load in magnetic field [J]. Nonlinear Dynamics, 2016, 85(3): 1825-1835.
[10] 胡宇達(dá),王彤. 磁場中導(dǎo)電旋轉(zhuǎn)圓板的磁彈性非線性共振[J]. 振動(dòng)與沖擊,2016,35(12):177-181.
HU Yuda, WANG Tong. Nonlinear resonance of a conductive rotating circular plate in magnetic field [J]. Journal of Vibration and Shock, 2016, 35(12): 177-181.
[11] BAYAT M, RAHIMI M, SALEEM M. One-dimensional analysis for magneto-thermo-mechanicl response in a functionally graded annular variable-thickness rotating disk [J]. Applied Mathematical Modelling, 2014, 38(19/20): 4625-4639.
[12] ZENKOUR A. On the magneto-thermo-elastic responses of fg annular sandwich disks [J]. International Journal of Engineering Science, 2014, 75(2): 54-66.
[13] 高原文, 周又和, 鄭曉靜. 橫向磁場激勵(lì)下鐵磁梁式板的混沌運(yùn)動(dòng)分析[J]. 力學(xué)學(xué)報(bào), 2002, 34(1): 101-108.
GAO Yuanwen, ZHOU Youhe, ZHENG Xiaojing. Analysis of chaotic motions of gemetrically nonliear ferromagnetic beam-plates exicited by transverse magnetic fields [J]. Acta Mechanica Sinica, 2002, 34(1): 101-108.
[14] 王省哲, 鄭曉靜. 鐵磁梁式板磁彈性初始后屈曲及缺陷敏感性分析[J]. 力學(xué)學(xué)報(bào), 2006, 38(1): 33-40.
WANG Xingzhe, ZHENG Xiaojing. Analysis on magnetoelastic initial post-buckling and sensitivity to imperfection for ferromagnetic beam-plates [J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(1): 33-40.
[15] BHANDARI A, KUMAR V. Ferrofluid flow due to a rotating disk in the presence of a non-uniform magnetic field [J]. International Journal of Applied Mechanics & Engineering, 2016, 21(1): 273-283.
[16] RAD A B, SHARIYAT M. Thermo-magneto-elasticity analysis of variable thickness annular fgm plates with asymmetric shear and normal loads and non-uniform elastic foundations [J]. Archives of Civil & Mechanical Engineering, 2016, 16(3): 448-466.
[17] PEI Y C, CHATWIN C, HE L, et al. A thermal boundary control method for a flexible thin disk rotating over critical and supercritical speeds [J]. Meccanica, 2017, 52(1/2): 383-401.
[18] SHARMA S, SANEHLATA Y. Finite difference solution of elastic-plastic thin rotating annular disk with exponentially variable thickness and exponentially variable density [J]. Journal of Materials, 2013 (5): 809205.
[19] HUSSAIN S, AHMAD F, SHAFIQUE M, et al. Numerical solution for accelerated rotating disk in a viscous fluid [J]. Applied Mathematics, 2013, 4(6): 899-902.
[20] JABBARI M, GHANNAD M, NEJAD M Z. Effect of thickness profile and fg function on rotating disks under thermal and mechanical loading [J]. Journal of Mechanics, 2016, 32(1): 35-46.