基于徑向基點插值法的旋轉Mindlin 板高次剛柔耦合動力學模型1)

杜超凡 * 鄭燕龍 * 章定國 周曉婷 *

* (揚州大學建筑科學與工程學院,江蘇揚州 225000)

? (南京理工大學理學院,南京 210094)

引言

在現代工程領域中,隨著航空航天和民用機械工程等領域的快速發展,對柔性多體系統在高速化及大變形工況下的可靠性要求越來越高.許多構件如太陽能帆板、風力發電機葉片、直升機旋翼及渦輪機葉片等都屬于柔性構件搭載于剛性主體的剛-柔耦合結構.這類結構在經歷大范圍剛體運動的同時還存在著柔性構件的變形,且相互耦合,呈現出的動力學行為十分復雜.由于對這類經歷大范圍運動柔性梁的剛-柔耦合動力學行為的研究比較成熟,因此常將這類板式結構簡化為一維梁進行建模分析[1-3].這種簡化往往只適用于長寬比很大的細長構件,對于長寬比相差不大的板式結構,如太陽能帆板等這種簡化顯然是不合理的,仿真結果會有很大誤差,必須用板的理論來建模分析[4-11].從國內外學者研究中發現,對做大范圍運動矩形板的研究有兩個共同點:一是變形場的離散方法都是假設模態法或有限元法[12-13];二是在建模過程中,大都采用忽略剪切變形的薄板理論,同時引入小變形假設,認為與非線性耦合項相關的一些高階項可以直接去掉,并不影響最終結果.這就導致兩點不足:一是變形場離散方法的單一,在旋轉柔性梁的分析中已經說明了假設模態法的不足[14-15],而有限元法又存在網格畸變和難以構造高階形函數的問題;二是隨著柔性構件尺寸、厚度的增大或構件材料柔性的增加,被省略的剪切變形和非線性耦合項相關的高階項的影響將逐漸增大[16-19],必須考慮其影響建立更加精確的高次耦合動力學模型.陳思佳等[20]考慮了非線性耦合項的高階項,建立了既適用于小變形又適用于大變形問題的旋轉柔性梁的高次剛柔耦合動力學模型,彌補了一次近似耦合模型在處理柔性梁大變形問題上的不足.因此,為了得到更精確的仿真結果,使用新的變形場離散方法和依據中厚板理論建模[21-23]獲得更精確的動力學方程將顯得很有必要.無網格法作為一種新的變形場離散方法,已廣泛應用于板的靜力學和自由振動分析中[24-30],但在柔性多體系統領域鮮有報道[31-32].將其應用到做大范圍轉動的矩形板中,將極大豐富柔性多體系統領域中的變形場離散方法.

本文采用無網格徑向基點插值法(RPIM)描述柔性板的變形場,考慮經典薄板理論中被忽略的剪切變形的影響,保留動能中與非線性耦合變形量相關的高階項,同時考慮了板的平動,建立了較為完整的旋轉矩形板高次剛-柔耦合動力學模型.該模型既能處理薄板問題,又能處理中厚板問題,既能處理小變形問題,又能處理大變形問題.將仿真結果與傳統的假設模態法和有限元法對比,表明徑向基點插值法的正確性,同時說明無網格法用于該領域的可推廣性.

1 旋轉Mindlin 板動力學建模

1.1 系統的動能和勢能

圖1 所示為做大范圍運動的矩形板模型,其中坐標系O-XYZ為慣性坐標系,o-xyz為連體坐標系,且連體坐標系所在平面與薄板未變形的中面重合,其3 個方向的單位矢量分別為e1,e2,e3.板長為L,寬為H,厚度為h,密度為ρ,彈性模型為E,泊松比為μ.

圖1 作大范圍運動的矩形板Fig.1 Rectangular plate with large overall motion

P0為變形前板中面上一點,變形后移至P點,變形位移矢量為u,在連體坐標系下各分量為(u1,u2,u3).基于Mindlin 板理論,板中任意一點變形位移矢量為u,在連體坐標系下為各分量為(u1,u2,u3).其中縱向位移u1,u2可用板中面上(z=0)的變形u01,u02表示為

式中,w1和w2分別為板中面內沿x和y方向的實際伸長量,為二次耦合變形量,傳統的零次近似模型建立在結構動力學小變形假設的基礎上,并沒有考慮二次耦合變形量.φx和 φy分別為橫截面相對于y軸和x軸的轉角.任意P點在慣性坐標系O-XYZ下的速度矢量可表示為

式中,Vo和 ωA分別為連體坐標系相對于慣性坐標系的速度和角速度矢量,ρ0為點P0在連體坐標系中位置矢量,VPA為P點相對連體坐標系的速度矢量.各矢量在連體坐標系的分量形式為

將以上各式代入式(2)中,可得

系統的動能為

考慮剪切變形的影響,任一點處的應變為

對各項同性材料,其應力為

系統勢能為

式中,U1和U2分別為板面內變形能、板橫向彎曲和剪切變形能,表達式分別為

1.2 系統的動力學方程

將問題域離散為若干個節點,將矩形板問題域用N個場節點表示,如圖2 所示.形函數由支持域中所包含的節點形成,具體表達式在文獻[31]中有詳細闡述.

圖2 矩形板無網格法離散Fig.2 Meshless discretization of rectangular plates

與薄板理論建模不同,Mindlin 板理論由于考慮橫向剪切應力對變形的影響,因此變形前垂直于中面的直法線變形后將變為曲線,橫截面相對于x軸和y軸的轉角不再是薄板理論中撓度u3的一階導數,為獨立變量,表示為

式中,n為計算點支持域內的節點數,φ1i(x,y),φ2i(x,y),φ3i(x,y) ,φ4i(x,y) ,φ5i(x,y) 為支持域內對應節點的形函數,通常令 φji(x,y) (j=1,2,···,5) 相同.q1(t)和q2(t) 為相應節點面內縱向變形位移列陣,q3i(t)為相應節點橫向變形及截面轉角組成的列陣.對 應u3和兩個轉角的形函數行陣及q3i(t) 的表達式為

為表述方便,以下表達式中略去自變量x,y,t.

令 Φ4=zφ4,Φ5=zφ5,將式(15)代入式(1),得變量u1,u2及其速度為

式中,H1(x,y)和H2(x,y) 為耦合形函數,表示如下

下標中“,”表示對坐標求偏導數.

式中,a01,a02,a03為基點加速度在連體坐標系下的分量,表達式為

各常數陣為

式(24)即為系統的高次剛-柔耦合動力學方程,該模型既能處理大變形問題又能處理小變形問題,彌補了一次近似模型在處理大變形問題時的不足.其中帶雙下劃線的項是由于考慮了非線性二次耦合變形量的高階項而推導出的.若去掉這些雙下劃線的項,則式(24)退化為一次近似模型;若去掉單下劃線的項,則退化為零次近似模型.

2 數值仿真

2.1 剪切閉鎖

當使用Mindlin 厚板理論分析薄板問題時,往往會發生剪切閉鎖現象.為了使推導的動力學方程能正確仿真薄板問題,應首先對剪切閉鎖現象進行研究.

考慮如圖3 所示的矩形板,其參數為:長a=10.0 m,寬b=10.0 m,彈性模量E=1.0 GPa,泊松比μ=0.3,厚度為h.定義一個撓度系數 ξ,表達式為

圖3 矩形板受力圖Fig.3 Force diagram of rectangular plate

式中,wmax為板的最大撓度值,P為集中力載荷,D為板的彈性剛度,表達式為

表1 為不同節點個數下,使用徑向基點插值法求得的四邊簡支、厚長比h/a=0.01 且中心點受集中力P=100 N 作用的薄板撓度系數.CPT 和FSDT 分別表示經典薄板理論和Mindlin 板理論.從表中可知,當節點個數為17 × 17 左右,可得到與文獻基本相同的解,基于Mindlin 板理論的解比經典薄板理論解更精確.在以下仿真中,離散節點個數都選取17 ×17.

表1 四邊簡支、中心點受集中力作用的薄板撓度系數 (h/a=0.01)Table 1 Deflection coefficient of thin plates subjected to concentrated force at the center point with SSSS boundary condition (h/a =0.01)

表2 為不同邊界條件的厚長比h/a=0.1,中心點受集中力P=100 N 作用的厚板撓度系數,分別使用經典薄板理論(CPT)和Mindlin 板理論(FSDT).為表述方便,做如下規定,板四邊按左上右下順序排列,簡支邊界用S 表示,固支邊界用C 表示,如SCSC 表示左右邊界簡支,上下邊界固支,SCCS 表示左下簡支,右上固支.從表中可知,基于Mindlin 板理論的結果與文獻[33]提供的解析解結果基本一致,且總比基于經典薄板理論的結果大.同時說明經典薄板理論只能處理薄板問題,而Mindlin 板理論即可以處理厚板問題,又可以處理薄板問題,并且精度更高,因此具有更廣泛的應用價值.

表2 不同邊界條件和受力情況下厚板撓度系數(h/a=0.01)Table 2 Deflection coefficient with different boundary conditions and forces (h/a=0.1)

當板厚度更薄即厚長比更小時,基于Mindlin 板理論將出現剪切閉鎖現象.在有限元法中,通常通過構造高階單元形函數來消除剪切閉鎖現象.根據這一思想,在無網格法中同樣可以通過構造高階形函數來消除剪切閉鎖現象,而可方便的構造高階形函數正是無網格法的優勢,通過添加高階多項式基函數即可完成.表3 表示受集中力作用,四邊簡支的矩形板在不同厚長比下的撓度系數 ξ .從表中可明顯看出,當厚長比小于0.001 時,添加低階多項式的徑向基點插值法出現剪切閉鎖現象,求得的撓度系數比解析解小很多.當添加的多項式為15 和18 項時,可消除剪切閉鎖現象.厚長比為0.000 1 的板是極端情況,現實中并不存在,因此對于基于Mindlin 板理論的徑向基點插值法,添加15 項多項式基函數即可消除剪切閉鎖現象.需要注意的是,此時需增大支持域尺寸,否則將導致力矩矩陣不可求逆,支持域的無量綱尺寸通常取為 αs=4.1 .當厚長比較大即不發生剪切閉鎖現象時,添加多項式基函數對解的精度影響并不大,因此對于中厚板,為提高計算效率,可選擇不添加多項式基函數.圖4 所示為受集中力作用四邊簡支板在不同厚長比下,數值解與文獻提供的解析解的比值.從圖中可更直觀地看出,當厚長比小于0.001 時,添加低階多項式基函數發生明顯的剪切閉鎖現象,當多項式基函數為15 和18 項時,可消除剪切閉鎖現象.當厚長比較大時,多項式基函數對結果影響并不明顯.

圖4 多項式基函數對剪切閉鎖的影響Fig.4 Influence of polynomial basis function on shear locking

表3 受集中力作用四邊簡支板不同厚長比下的撓度系數Table 3 Deflection coefficient with different aspect ratios

2.2 大范圍運動已知的動力學仿真

考慮如圖5 所示的作定軸轉動的中心剛體-矩形懸臂板,與中心剛體固連.中心剛體以角速度 ω 繞y軸旋轉,板的材料參數為:長度a=1.828 8 m,寬度b=1.219 2 m,厚度h=2.54 mm,彈性模量E=70 GPa,密度ρ=2000 kg/m3,泊松比取0.3,中心剛體半徑R取0.ω1=ω3=0, ω2=ω , ω˙1=ω˙3=0 ,ω˙2=,給定的角速度規律為

圖5 作定軸轉動矩形薄板Fig.5 Configuration of rotating cantilever plate

式中,T=30 s .

圖6 為Ω=20 rad/s 時,本文方法與參考文獻[6]中矩形薄板外側角點M的橫向變形位移歷程的對比.從圖中可以看出,本文RPIM 的仿真結果與參考文獻的仿真結果高度吻合,驗證了本文動力學模型的正確性.圖7 所示分別為 Ω=0.2 Hz和Ω=0.75 Hz時板外側角點M的橫向變形.從兩圖中可以看出,當轉速較低時,零次近似模型、一次近似模型和高次模型的差別并不大,但零次近似模型的最大值比一次近似模型和高次模型的稍大.當轉速較高時,零次近似模型與另兩種模型出現很大區別.零次近似模型結果發散,仿真失敗,而一次近似模型和高次模型結果收斂,且仿真結果高度吻合.這與旋轉柔性梁中的分析結論相同,說明零次近似模型忽略耦合變形量存在建模理論的缺陷,其剛度陣隨著轉速的增大而減小,產生動力柔化效應,這也是造成在轉速較低時其變形比一次近似模型大的原因.而一次近似模型和高次模型考慮了耦合變形量,其剛度陣隨著轉速的增大而增大,產生動力剛化效應,符合實際情況,同時也說明本文所建高次模型的正確性.進一步研究發現,對于大范圍運動已知的情況,隨著轉速的增加,高次模型并沒有比一次近似模型優越,相反使計算效率降低很多.因此對于該類問題的研究優先選用一次近似模型.

圖6 外側角點橫向變形比較Fig.6 Comparison of lateral deformation of outer corner point M

圖7 矩形薄板外側角點的橫向變形Fig.7 Lateral deformation of outer corner point M

圖8 為 Ω=0.2 Hz 時使用不同離散方法獲得的一次近似模型板外側角點M的橫向變形,其中假設模態法和有限元法基于經典薄板理論,徑向基點插值法基于Mindlin 板理論.有限元法分別使用8 ×4 個矩形單元和2 × 8 × 4 個三角形單元,徑向基點插值法離散為16 × 8 個節點.從圖中可看出,各離散方法的仿真結果基本一致.圖8(b)為勻速轉動階段的局部放大圖,同樣可看出各離散方法的響應頻率及振幅基本一致.圖9 為 Ω=0.2 Hz 時使用不同離散方法獲得的一次近似模型板外側角點M的橫向變形速度.如圖所示,各離散方法的仿真結果基本一致,從局部放大圖中可知,大范圍運動轉速恒定時的速度響應頻率和振幅也基本一致.

圖8 矩形薄板外側角點橫向變形(Ω=0.2 Hz)Fig.8 Lateral deformation of outer corner point M(Ω=0.2 Hz)

圖9 矩形薄板外側角點橫向變形速度(Ω=0.2 Hz)Fig.9 Lateral deformation rate of outer corner point M(Ω=0.2 Hz)

圖10 為 Ω=0.75 Hz 時使用不同離散方法獲得的一次近似模型板外側角點M的橫向變形.從圖中可看出,各方法的仿真結果基本一致.圖10(b)為勻速轉動階段的局部放大圖,各方法的響應振幅和頻率并沒有 Ω=0.2 Hz 時吻合的好,但也相差不大.其中假設模態法、三角形單元有限元法的仿真結果最吻合,徑向基點插值法的響應振幅最大,而矩形單元有限元法響應振幅最小.從縱坐標的數量級來看,這種差別是極微小的,都能滿足實際的工程需要.圖11為 Ω=0.75 Hz 時使用不同離散方法獲得的一次近似模型板外側角點M的橫向變形速度.同樣的,各方法的仿真結果基本一致.圖11(b)為勻速轉動階段的局部放大圖,如圖所示,與橫向變形的局部放大圖中的變化相同,我們可得出同樣的結論.

圖10 矩形薄板外側角點橫向變形(Ω=0.75 Hz)Fig.10 Lateral deformation of outer corner point M(Ω=0.75 Hz)

圖11 矩形薄板外側角點橫向變形速度(Ω=0.75 Hz)Fig.11 Lateral deformation rate of outer corner point M(Ω=0.75 Hz)

圖12 為 Ω=10 Hz 時使用不同離散方法獲得的一次近似模型板外側角點M的橫向變形.從圖中可看出,仿真結果在10 s 之前基本一致.隨著仿真時間的增加,有限元法結果發散,計算失敗.而徑向基點插值法結果收斂,體現了無網格法計算上的優勢.

圖12 矩形薄板外側角點橫向變形(Ω=10 Hz)Fig.12 Lateral deformation of outer corner point M(Ω=10 Hz)

圖13 為在三維空間轉動的矩形板,矩形板固結于一繞Y軸作定軸轉動的中心剛體上.慣性坐標系O-XYZ的原點O和矩形板連體坐標系o的距離r=0.2 m.連體坐標系中單位方向矢量e1和e2分別和未變形矩形板長度和寬度平行,e3垂直于矩形板中面.(e1,e2,e3)可由慣性坐標系O-XYZ經過三次繞軸旋轉獲得,3 個歐拉角和系統所有參數與文獻[6]相同.圖14 為Ω=10 rad/s 時,矩形板外側角點P橫向變形同文獻[6]對比圖,從圖中可以看出,本文的仿真結果與參考文獻的仿真結果高度吻合,驗證了本文動力學模型處理三維算例的正確性.

圖13 三維空間轉動板Fig.13 Rotating plate in 3-D

圖14 矩形薄板外側角點P橫向變形(Ω=10 rad/s)Fig.14 Lateral deformation of outer corner point P(Ω=10 rad/s)

2.3 大范圍運動未知的動力學仿真

板的材料參數為:長度a=1.0 m,寬度b=0.5 m,厚度h=2.5 mm,彈性模量E=70 GPa,密度ρ=3000 kg/m3,泊松比取0.3,中心剛體半徑R取0.假定施加在中心剛體繞y軸方向的外力偶矩為

其中,T=2 s.

圖15 所示分別為τ0=5 N·m 和τ0=70 N·m 時一次近似模型和高次模型板外側角點M的橫向變形.從兩圖中可以看出,當τ0=5 N·m 時,一次近似模型和高次模型的仿真結果高度吻合,此時柔性板外側角點最大變形位移約為0.03 m,屬于小變形,說明兩種模型均適于小變形情形.當驅動力偶矩τ0=70 N·m,一次近似模型結果發散,仿真失敗,而高次模型仿真結果仍然收斂.此時柔性板外側角點最大變形位移約為0.4 m,對于長為1 m 的板而言屬于大變形.由此可知,對于大變形情形,一次近似模型并不適用,而高次模型仍然適用.同時也說明二次耦合變形量的高階項在大變形情形下尤為重要,不能忽略.

圖15 矩形薄板外側角點橫向變形Fig.15 Lateral deformation of outer corner point M

3 結論

本文采用無網格徑向基點插值法描述柔性板的變形場,基于一階剪切變形理論,保留動能中有關非線性耦合項的高階量,通過構造高階形函數避免了徑向基點插值法出現剪切閉鎖的現象,建立了能處理不同厚長比的作大范圍運動矩形板的高次剛-柔耦合動力學模型.通過靜力學算例驗證避免剪切閉鎖方法的正確性.分別運用假設模態法、有限元法和徑向基點插值法及兩種板的建模理論對作定軸轉動的中心剛體與懸臂板的動力剛化問題進行動力學仿真,其仿真結果與已有文獻一致,說明本文徑向基點插值法建立的模型是正確的,也說明了在同等條件下,無網格法具有計算上的優勢.比較了傳統的零次近似模型、一次近似模型和高次模型之間的差異,發現傳統的零次近似模型只能適用于轉速很低的情況,適用范圍狹窄,而一次近似模型和高次模型既適用于低轉速又適用于高轉速的情況,證明了考慮耦合變形量的必要性.研究發現,高次模型由于考慮了非線性耦合變形項的高階量,相對于一次近似模型,計算效率降低很多,因此仿真時優先選用一次近似模型;一次近似模型只適用于小變形情況,而高次模型既適用于小變形情況,又適用于大變形情況.