V型波紋夾層板彎曲問題的高階剪切變形理論計算方法

王小明，潘 曼，魏 強，2

（1.中國艦船研究設計中心，武漢430064；2.船舶振動噪聲重點實驗室，武漢430064）

0 引 言

夾層板不僅具有比強度高、比剛度高的優點，還具有吸聲隔振的效果，一直以來就是工程關注的焦點和科學研究的熱點。夾層板中的芯層可設計性非常強，可以是PVC泡沫、玻璃纖維、金屬薄片等，這種特點為夾層板的工程應用提供了更多選擇。其中，金屬薄片又有很多的形式，可以是波紋型（V型）、I 型、Y 型、O 型、Z 型等。長期以來，面板與芯層的連接加工技術制約著金屬夾層板的生產效率；近年來，工程師們將激光焊接技術應用到金屬夾層板的生產制造中，使得金屬夾層板的生產加工更加便捷精準。關于它的工程應用，人們對它充滿期待，前景看好，航空航天、船舶海洋工程都是其應用的方向[1-2]。關于波紋夾層板的彎曲問題，在長期的研究探索中逐漸形成了兩個研究方向：第一個研究方向是把波紋夾層板等效成正交異性板，然后按照正交異性薄板的彎曲求解；第二個研究方向是根據一階剪切變形理論模型，分別計算上下面板的變形能和芯層的變形能，然后采用最小勢能原理，列出波紋夾層板的平衡微分方程[3]。第一個研究方向中等效方式又可以分為兩類：一類是運用變形等效或者能量等效的原理求解芯層的等效彈性模量或者等效剛度[4-9]，夾層板就等效成層合板，上下面板是各向同性板，芯層是正交異性板；另一類是根據芯層的實際形狀，不需要單獨求解芯層的等效參數，直接根據整體變形等效的方法求解整體彎曲剛度和剪切剛度[10-18]，夾層板整體等效成正交異性板。整體等效方式最為著名的是Libove等[11]的等效剛度方法，后來，Fung和Tan等在Libove的等效剛度法基礎上，推導了Z 型、C 型和V 型夾層板的等效剛度[12-15]。Chang 等[19]采用了Libove 的方法計算出等效剛度后，研究了夾層板內力、變形與夾層板的截面參數之間的變化規律；Wang等[7]把波紋夾層板芯層看作正交異性體，采用“同外力，等變形”的方法分別求解了芯層三個方向的等效彈性常數，最后采用正交異性板理論計算了夾層板的位移，與有限元算法比較，差別小于5%。文獻[20]提出波紋夾層板在垂直波紋母線的平面內，其抗剪剛度無限大，垂直于波紋方向平面內的抗剪剛度為有限值；在這一前提下，推導了波紋夾層板的彎曲微分方程，但是并未介紹如何計算垂直波紋方向平面內的抗剪剛度。實際工程應用中，為了達到理想的強度和剛度，必然會將波紋夾層板的整體厚度、上下面板厚度和芯層板厚度都增大，倘若夾層板的整體厚度達到厚板的范疇，則不能忽略芯層的抗剪作用，即不能采用經典薄板理論，也不能采用一階剪切變形理論，因為一階剪切變形理論不滿足上下自由表面條件，這時，應該采用高階剪切變形理論求解。目前，關于波紋夾層板彎曲問題的文獻多數是采用正交異性薄板理論或者一階剪切變形理論算法，采用高階剪切變形理論求解方法的公開文獻很少。對于大厚度的波紋夾層板彎曲問題，倘若不采用高階理論算法，將導致計算誤差偏大。因此，有必要研究大厚度波紋夾層板的高階剪切變形理論求解方法，為工程設計和計算提供指導。本文提出一種“同外力，等能量”的等效方法計算波紋芯層等效彈性參數；求出等效彈性參數之后，把夾層板等效成層合板，應用高階剪切變形理論，列出波紋夾層板的彎曲微分方程，應用傅立葉級數法求解該方程，最終確定變形和應力。

1 芯層等效彈性模量和泊松比

關于芯層等效彈性參數的計算方法很多，Bartolozzi[8-9]等推導的方法被諸多學者采用，在文獻[8-9]中推導出的芯層為正弦曲線的波紋板的等效模量。Cheon 和Park 等[5-6]采用了“同外力，等能量”的等效方法推導正弦波紋的等效剛度。本文則采用“同外力，等能量”的等效方法推導V 形波紋的等效彈性模量和泊松比。波紋夾層板和坐標系如圖1所示，x軸沿著芯層的母線方向，y軸沿著芯層的波紋方向，z軸垂直于xoy平面，指向下方；xoy平面位于夾層板芯層的中面。z＜0一側的面板叫上面板，z＞0一側的面板叫下面板。上面板厚度為tt，下面板厚度為tb，芯層板厚tc，芯層凈高hc，芯層周期長度為lc，芯層半周期斜面邊長為l。芯層傾斜面與面板夾角為θ。分析時，以x方向長度為單位長度作為研究對象。

圖1 夾層板與坐標系統Fig.1 Sandwich panel and coordinate system

1.1 yz平面等效模量

在夾層板芯層頂點施加一水平力F，頂點將發生微小的水平位移δ（如圖2 所示）。將這個物理模型單獨研究，選取其中的半單元芯層（圖2中的①）進行分析，如圖3 所示。為了保證芯層與純剪切狀態等效，則在芯層頂點加載水平載荷F時，芯層上方頂點（圖3中的B點）不能有垂向位移，下方的頂點（圖3中的A點）不能發生彎曲轉角。這就意味著A點是固定鉸支座，提供水平方向和豎直方向的支反力，B點是活動鉸支座，僅提供豎直方向的支反力，簡化后的物理模型如圖3 所示。為了簡化公式推導，把坐標系原點移到下方芯層頂點的位置，同時轉換視角，以符合坐標習慣。

圖2 芯層xz面剪切變形Fig.2 Core sheets shear deformation in xz plane

研究中認為，芯層沿x方向（垂直于紙面方向），長度為單位長度，根據芯層的平衡條件，容易求出下方頂點的位置處支反力R=Ftan(θ)。

根據圖3的坐標系，容易得出芯層①中心線方程為

圖3 計算力學模型Fig.3 Mechanics model for calculating

因為載荷作用在芯層上方頂點B的位置，則芯層的受力狀態相當于二力桿，即芯層內力只有軸力，沒有剪力和彎矩。

芯層①中的軸力（拉為正，壓為負）表示為

以上表明芯層中僅存在軸力，則芯層半單元中的總應變能為

式中，L為半周期芯層橫截面中心線，如芯層①的中心線，l為芯層①的長度，A為橫截面面積，即A=tc。與芯層單元等效的模型如圖4所示，則等效模型的應變能為

圖4 計算等效均質模型Fig.4 Equivalent homogenous model for calculating

根據外力相等，應變能相等的原理，式（3）與式（4）相等，即可算出。

1.2 x方向等效彈性模量

均勻等效模型的應變能為

1.3 等效模量和

圖5 計算力學模型Fig.5 Mechanics model for calculating

原模型與等效模型具有相同的剪切應變能，即式（9）與式（10）相等，求出。

圖6 計算力學模型Fig.6 Mechanics model for calculating

1.4 y方向等效彈性模量

設想芯層和等效均質模型均受到彎矩M作用后在yz平面內彎曲（如圖7和圖8所示），則兩者的彎曲應變能應當相等。夾層板芯層的應變能為

圖7 計算力學模型Fig.7 Mechanics model for calculating

圖8 計算等效均質模型Fig.8 Equivalent homogenous model for calculating y

等效均質模型的彎曲應變能為

根據變形能相等，即式（13）與式（14）相等，求出

1.5 泊松比和

1.2節中的受力狀態正好是x方向受均勻拉力，則x方向的應變為

沿芯層傾斜方向的應變為

沿y方向的應變為

正交異性體彈性模量與泊松比滿足方程

2 夾層板彎曲方程

2.1 位移模式

關于板的高階理論有多種表達形式，本文采用Touratier提出的正弦函數表達式[21-22]。夾層板的坐標系統如圖1所示，則夾層板的位移模式可以表示成

式中，h為整個夾層板厚度，即h=hc+tt+tb。

2.2 本構方程

夾層板的正應變和剪應變可以采用彈性力學幾何方程表達式求解：

上下面板是各向同性板，本構方程直接采用各向同性體的胡克定律，上下面板的正應力和剪應力可以表示為

式（30）-（32）中，E、G和μ分別代表上下面板材料的彈性模量、剪切模量和泊松比。這里并沒有考慮面板的橫向剪應力是注意到上下面板處于遠離中和軸位置，根據式（28）和式（29）計算的剪應變很小，可以忽略不計。

芯層等效為正交異性板，則其應力可以表達為

2.3 最小勢能原理

根據最小勢能原理，結構系統處于穩定狀態，則它的位移必然使得整個系統的動勢（動能與勢能之差）達到最小值。對于夾層板的靜力分析，則系統沒有動能，僅有勢能。勢能包括兩部分，一部分是系統的應變能，另一部分是外力的功。

式中，A表示夾層板所在的xoy平面范圍，U1和W1分別表示應變能面密度和單位面積上外力所做的功，δ表示變分。

系統的應變能包括三部分，分別是上、下面板應變能和芯層應變能。

其中，

對于作用在坐標(x0,y0)的集中力，可以表達成為

式中，δ為狄拉克函數。

集中力所做的功為

對于式（36），根據變分法，令

則式（36）的歐拉方程為

系數A1、B1、C1等與夾層板的截面參數相關。該方程非常復雜，篇幅太長，不便于把系數表達式全部列出，在末尾附錄A中給出了方程推導過程的Mathematica程序。

3 彎曲方程與求解

對于x方向長為a，y方向長為b的四邊簡支的夾層板，邊界條件為

可以設方程（47）-（49）的解為雙傅立葉級數形式

其實，式（51）-（53）已經滿足了式（50）的邊界條件。

集中力P作用點的坐標為（x0，y0），可以把集中力展開成雙傅立葉級數：

將式（51）-（54）代入式（47）-（49），然后比較兩邊的系數，可以分別求出上述式（51）-（53）中雙傅立葉級數的系數wmn、gmn和hmn。

4 夾層板的應力

求出夾層板的變形之后，還需要繼續求解夾層板的應力。

將第3章中求出的三個未知函數?x、?y和w對應代入到式（30）-（32）中則可以求出夾層板上下面板的應力。

芯層的實際應力則不能直接由式（33）-（35）求出，也沒有方法準確求出，這是等效方法的共同點。夾層板彎曲問題，最大正應力一定出現在上下面板上。

5 算例與討論

5.1 薄板模型

夾層板四邊簡支，長邊a=2000 mm，短邊b=1500 mm，上下面板厚度tt=tb=3 mm，芯層板厚度tc=2 mm，芯層凈高度hc=40 mm，芯層周期間距lc=50 mm，上下面板和夾心都是同樣材料，彈性模量E=Ec=2.1×105MPa，泊松比μ=0.3。集中力作用于上面板的中心點，作用點的坐標為（1000,750），集中力的大小為P=2×104N。計算模型即為在圖1的基礎上增加四邊簡支邊界條件和板中心增加集中力作用。

為了形成對比，驗證方法的準確性，分別用2種方法計算這個算例。方法1：本文計算方法，雙傅里葉級數的項數迭代到m=n=15，此時位移和應力（遠離載荷作用點的應力）均收斂；方法2：ANSYS 有限元計算方法，面板和芯層都采用Shell 181單元，網格劃分成19 200個單元。

構件的變形和強度是工程關注的焦點，最關注的則是最大變形量和最大應力，以及出現的具體位置。上述算例中，根據構件和載荷的對稱性，最大變形和最大應力可能出現在板的中心點。位移的級數表達式（51）-（53）是收斂的，一般m、n分別取到5 之后則變化非常小。應力表達式（30）～（32）在集中載荷作用點附近是不收斂的，應力為無限大，這一性質與單層板是相類似的。因此不能試圖去比較x=a/2和y=b/2位置處的應力分布。

上述2種計算方法的結果中y=b/2和x=a/2兩個典型位置的位移分布如圖9和圖10所示；4個典型位置的應力分布如圖11～14所示。

由圖9和圖10可以看出，上述2種方法關于位移的計算結果極為接近，誤差普遍很小。圖11和圖13 中關于x向應力的分布，本文方法計算結果與有限元計算結果也吻合一致。圖12 和圖14 中關于y向應力的分布，本文方法是光滑曲線，而有限元計算結果則是周期振蕩的曲線，本文方法正好是有限元方法振蕩波峰的插值光滑連線。經過反復的試算發現，振蕩的周期與網格尺寸沒有關系，僅與夾層板的截面參數有關，并且振蕩周期等于夾層板波紋芯層周期。

圖9 y=b/2位置處的位移Fig.9 Deflection at y=b/2

圖10 x=a/2位置處的位移Fig.10 Deflection at x=a/2

圖11 y=b/3位置處下面板下表面x向應力Fig.11 x direction stress of bottom facesheet lower surface at y=b/3

圖12 x=a/4位置處下面板下表面y向應力Fig.12 y direction stress of bottom facesheet lower surface at x=a/4

圖13 y=b/4位置處上面板上表面x向應力Fig.13 x direction stress of top facesheet upper surface at y=a/4

圖14 x=a/4位置處上面板上表面y向應力Fig.14 y direction stress of top facesheet upper surface at x=a/4

實際上，x向應力隨y坐標變化時，也表現出圖12 和圖14 的振蕩特性，并且與圖12 和圖14 一樣，面板與芯層頂點交匯的位置出現波峰，頂點中心的位置則出現波谷。這種振蕩正是非連續芯層對面板的反作用力所致。

應力沿波紋方向分布表現出振蕩性質，文獻[3]也有相同的結論。目前解析法尚無法構造出夾層板的位移模型，使得其同時滿足芯層板兩側的自由表面條件、上下面板的自由表面條件和芯層頂點與面板間的位移（應力）單值條件，因此解析法計算的結果也就不可能出現與有限元結果完全一致的振蕩周期和振蕩振幅。文獻[3]中的載荷是均布載荷，文中采用了整體應力與局部應力疊加的方式計算上面板應力，即夾層板彎曲應力作為整體應力，兩邊芯層頂點固支板彎曲應力作為局部應力。均布載荷工況下，這種算法計算的上面板應力分布與有限元結果吻合一致。對于本文中集中載荷，這種方法則無法使用，因為集中載荷作用在芯層頂點時，局部應力不存在；倘若集中載荷作用點并非是芯層頂點，則僅一塊局部板存在局部應力，也不能表現出周期振蕩性質。對于下面板應力，上述疊加方法也是無效的，文獻[3]并沒有討論下面板的應力，圖12顯示下面板應力沿波紋方向分布也存在振蕩性質。

5.2 厚板模型

夾層板四邊簡支，長邊a=2000 mm，短邊b=1500 mm，上下面板厚度tt=tb=10 mm，芯層板厚度tc=8 mm，芯層凈高度hc=300 mm，芯層周期間距lc=50 mm，上下面板和夾心都是同樣材料，彈性模量E=Ec=2.1×105MPa，泊松比μ=0.3。集中力作用于上面板的中心點，作用點的坐標為（1000，750），集中力的大小為P=2×106N。

同樣采用5.1 節中的2 種方法計算夾層板的位移，不過在方法1 中，為了保證應力的收斂性，雙傅里葉級數的迭代項數更多，達到m=n=25。圖15～20 分別列出了與5.1 節中薄板對應位置處的位移和應力分布。

圖15 厚板y=b/2位置處的位移Fig.15 Deflection of thick plate at y=b/2

圖16 厚板x=a/2位置處的位移Fig.16 Deflection of thick plate at x=a/2

圖17 厚板y=b/3位置處下面板下表面x向應力Fig.17 Bottom facesheet lower surface x direction stress at y=b/3 for thick plate

圖18 厚板x=a/4位置處下面板下表面y向應力Fig.18 Bottom facesheet lower surface y direction stress at x=a/4 for thick plate

圖19 厚板y=b/4位置處上面板上表面x向應力Fig.19 Top facesheet upper surface x direction stress at y=b/4 for thick plate

圖20 厚板x=a/4位置處上面板上表面y向應力Fig.20 Top facesheet upper surface y direction stress at x=a/4 for thick plate

這一算例中，（hc+tt+tb）/b＞1/5，屬于厚板的范疇。目前，商業有限元計算軟件中，關于厚板的計算都采用一階剪切變形理論，因此，方法2 是采用一階剪切變形理論。在圖15～20 中，2 種方法的計算結果已經不再保持薄板范疇的吻合狀態。位移分布和x向應力分布在遠離載荷作用點時，2種計算方法差別較小，越靠近中心，2種計算方法差別越大。y向應力則僅在少數位置點2種計算方法結果保持一致，多數位置差別較大。兩種計算方法關于峰值有一個共同的特點：集中載荷作用下的厚夾層板，在平行于矩形板邊緣的橫剖面上，本文方法的位移峰值和應力峰值都比有限元結果偏大。

方法2沿波紋方向的應力分布依然表現出振蕩性，其振蕩的波幅僅在中心位置附近較大，遠離中心則振幅很小。眾所周知，厚板單層板彎曲問題需要采用高階剪切理論或三維彈性理論計算，所以計算厚板波紋板彎曲問題也應當采用高階剪切變形理論。以本文計算方法結果為基準，在上述算例中，有限元法位移峰值的差別為-16.3%，x向應力峰值差別為-26.6%，y向應力峰值差別為-55.4%。

6 結 論

本文通過將夾層板中間芯層等效成正交異性板，根據應變能等效原理求解芯層各項等效彈性模量和泊松比。運用高階剪切變形理論和最小勢能原理，推導夾層板的彎曲變形微分方程，然后采用雙傅里葉級數的方法求解該微分方程。

通過算例驗證，在薄板范疇內，本文方法計算的位移與有限元法結果吻合良好；本文方法計算的面板應力在芯層頂點位置與有限元一致，其他位置則比有限元結果偏大；由于波紋芯層與上下面板的相互作用，上下面板彎曲應力沿波紋方向分布表現出振蕩性，越靠近外載荷作用點，振蕩振幅越大，本文方法計算的應力分布恰好是有限元法結果中振蕩波峰點的光滑連線。這一點對工程設計而言，計算結果是偏于安全的；這也說明薄板范圍內的波紋夾層板，高階剪切變形理論計算結果與一階剪切變形理論結果是基本一致的。

通過算例驗證，在厚板范疇內，集中載荷作用下的夾層板位移分布和x向應力分布在遠離載荷作用點時，本文方法和有限元方法計算結果差別較小，越靠近中心，兩種計算方法差別越大；y向應力則僅在少數位置點時本文方法與有限元方法結果保持一致，多數位置差別較大。在平行于矩形板邊緣的橫剖面上，對于位移和應力的分布峰值而言，本文方法結果普遍大于有限元結果。因此，在工程設計中，推薦采用高階理論算法校核波紋夾層板的位移和應力。