邊界元法計算復合材料的雙向拉伸強度①

蔡登安，周光明，曹 然，黃 翔

(南京航空航天大學機械結構力學及控制國家重點實驗室，南京 210016)

0 引言

纖維增強復合材料以其輕質、高強、抗疲勞、減振、耐高溫及可設計等一系列優點，在航空航天、能源、交通、機械和建筑等部門日益獲得了廣泛應用。一般而言，構件中復合材料多受雙向或多向應力狀態［1］。復合材料的測試技術已能準確表征纖維增強復合材料的單向力學行為。然而，同樣的試驗方法并不能完全適用于復合材料多向甚至只是雙向響應問題［2］。對于受復雜應力狀態的復合材料，僅僅研究其單向載荷下的強度是遠遠不夠的。因此，不管是對強度準則進行試驗驗證，還是從損傷角度探討復合材料在復雜應力狀態下的破壞，都需要對復合材料試樣進行雙向加載研究。獲得雙向應力狀態的傳統方式是采用圓管形試件，目前更多的研究者把焦點放到了十字型試樣上［3－4］，其優點在于十字型試樣雙向拉伸試驗能直觀地反映板殼復合材料的雙向受力狀態，因而成為當下最受重視的一種雙向拉伸試驗方法。

邊界元法(BEM)是一種將研究的連續體區域的邊界劃分成有限個單元的數值方法，在某些領域已經成為工程設計的重要工具［5－6］。相比有限元法(FEM)，邊界元法針對一些問題的分析有著獨特的優點。Lachat和Watson［7］為邊界元法作為一種有效數值方法的研究作出了貢獻:他們提出了一種類似于有限元法的等參方法，并證明邊界元法可作為解決具有復雜構型問題的一種有效工具。借助邊界元方法，Pineda León等［8］研究了在不同溫度和應力狀態下材料的蠕變分析模型;劉云忠等［9－10］提出了正交各向異性復合材料裂紋問題的邊界元分析方法;商欣萍等［11］用邊界元法從理論上分析非織造土工織物的拉伸性能，與實驗結果有較好的一致性;Araújo等［12］針對碳納米管增強復合材料的微觀結構分析，編寫了邊界元并行計算程序，并據此提出了一種可擴展的復合材料邊界元并行計算代碼;基于線性夾雜模型，Wang等［13］將邊界元方法應用于纖維增強復合材料的3D大尺度熱分析問題的研究中。

針對十字型試樣雙向拉伸實驗研究，本文采用邊界元法，從理論上研究了交織玻璃纖維增強復合材料的雙向拉伸強度。建立復合材料雙向拉伸邊界元分析模型，并與有限元結果及試驗結果進行比較，以期為復合材料雙向拉伸強度研究提供一些理論依據。

1 邊界元模型

1．1 基本方程

與長、寬尺寸相比，交織玻璃纖維增強復合材料十字型試樣的厚度很小，故可將其看作是平面正交各向異性材料。平面正交各向異性問題的基本方程可描述如下:

平衡方程:

幾何方程:

物理方程:

或

式中 ［S］和［Q］分別表示柔度矩陣和剛度矩陣，且［Q］=［S］－1。

平面正交各向異性材料有4個獨立的彈性常數，即 E1、E2、ν12和 G12，分別表示交織玻璃纖維增強復合材料2個方向拉伸時的彈性模量、泊松比和剪切模量。交織玻璃纖維增強復合材料的柔度矩陣［S］可表示為

邊界條件:

應力邊界條件

位移邊界條件

1．2 邊界積分方程與基本解

邊界元法的基礎是用Fredholm［14］積分方程的位勢積分來表示問題的解，它是將支配物理現象的微分方程轉化成邊界上的積分方程，然后再進行離散求解的數值計算方法。

圖1 平面區域V和邊界SFig．1 Planar region V and boundary S

對于平面正交各向異性材料，如圖1所示，由邊界S所包圍的平面區域上的任意點x的位移分量ui(x)，在不考慮體積力時，可由邊界上各點y的位移分量和面力分量通過以下邊界積分方程求得［5，11］:

式中 ui和ti分別為位移和面力;Uij(x，y)和Tij(x，y)分別為平面正交各向異性材料的位移基本解分量和面力基本解分量。

系數Cij(x)的值與x點處邊界形狀有關，邊界S光滑時，有

對于平面正交各向異性材料，位移和面力基本解可表示為［15］

其中，Kα、Ai、αi、Mi(i=1，2)為與彈性常數有關的參數;ri和θi分別為復平面中的矢徑和幅角，且有

式中 l1、l2分別為y點與x點橫、縱坐標之差;n1、n2分別為邊界外法線的方向余弦。

1．3 積分方程離散化

通過邊界積分方程，可求出邊界上所有未知的位移和表面力分量。然而，通常情況下很難用解析方法求解積分方程，只能采用數值解法。將區域的邊界劃分成N(N=N1+N2)個簡單單元。其中，N1個邊界單元屬于表面力已知，N2個邊界單元屬于位移已知，以N個邊界單元上的積分之和表示整個邊界上的積分。各邊界單元上的位移和表面力，可通過插值函數與邊界單元上有限個節點的位移和表面力以多項式近似。本文采用常單元將邊界積分方程離散，單元節點取為邊界單元中點，邊界單元的位移和表面力等于單元節點的位移和表面力。

設 yi(i=1，2，…，N)為邊界 Sj上的任意一點，對任意節點xi，邊界積分方程可離散為

式中 下標i，j分別為第i、第j個單元。

設:

則對任意節點，邊界積分方程可離散為

對于N個節點，根據各節點的離散化方程，可得到包含2N個方程的一次方程組，用矩陣形式可表示為

式中 ［H］、［G］均為2N×2N階的位移和面力系數矩陣;{U}、{T}分別為邊界單元節點的位移分量和表面力分量。

2 計算方法

根據以上理論分析，利用FORTRAN語言編寫計算程序，來分析交織玻璃纖維增強復合材料雙向拉伸強度，程序的計算流程見圖2。計算前，要輸入試樣4個獨立的彈性常數(縱向和橫向的彈性模量、泊松比及剪切模量);同時，要輸入試樣的強度參數、幾何尺寸、各單元的坐標以及邊界條件。

圖2 程序計算流程圖Fig．2 Flow-process diagram of program

交織玻璃纖維增強復合材料十字型試樣的幾何尺寸(單位:mm)，如圖3所示。十字型試樣中心36 mm×36 mm的正方形區域為中心試驗區，該區域滿足應力分布均勻，且應力水平較高，以保證初始破壞發生在中心試驗區;同時，試驗區剪應力遠小于正應力，以至可忽略不計。受雙向拉伸時，十字型試樣在實際拉伸試驗過程中，試樣的AB端和CD端由夾頭夾持在拉力作用下移動，EF端和GH端由夾頭固定。與試驗不同的是邊界元理論計算時，為避免附加彎矩的影響，需釋放EF端在1方向的位移約束和GH端在2方向的位移約束，如圖4所示。若分別給AB端和CD端一個2方向和1方向上的伸長量，通過邊界元法計算程序，便可得出十字型試樣各自由邊上邊界單元的位移、各內點的位移和應力，以及達到該伸長量所施加的雙向拉伸載荷。

對十字型試樣進行單元劃分，由于試件完全對稱，為減少計算量，只需對1/4試樣劃分單元，如圖5所示。

圖3 十字型試樣幾何形狀及尺寸Fig．3 Geometry and dimensions of cruciform specimen

圖4 十字型試樣的載荷與約束Fig．4 Loads and constraints of cruciform specimen

圖5 十字型試樣邊界單元劃分Fig．5 Boundary element of cruciform specimen

試樣的OI邊單元約束1方向的位移，OJ邊單元約束2方向的位移;于IB與JC端作用原十字型試樣拉伸載荷的一半。圖中標記的各點為指定的內點。通過計算各指定內點的應力，可了解十字型試樣在雙向拉伸過程中的內力分布情況。

要計算交織玻璃纖維增強復合材料雙向拉伸強度，需確定十字型試樣斷裂的條件，即試件在雙向載荷下，拉伸到何種程度達到斷裂。本文采用的強度準則為Tsai-Wu張量準則。平面應力狀態下，Tsai-Wu張量準則可表述為

十字型試樣中心試驗區的剪應力較小，忽略不計，因而準則方程可簡化為

式中 Xt、Xc、Yt、Yc分別為材料縱向的拉壓強度和橫向的拉壓強度，可由單向拉伸(壓縮)試驗測得;強度參數F12由試驗獲得雙向強度后計算得到。

邊界元程序計算過程中，若中心試驗區的縱向和橫向平均應力使準則方程左邊大于等于1，則判定十字型試樣斷裂。

3 實驗驗證

3．1 基本力學性能測試

在雙向拉伸試驗之前，對交織玻璃纖維增強復合材料進行基本力學性能測試，從而獲得交織玻璃纖維增強復合材料的工程彈性常數和縱、橫向拉壓強度，所得材料特性見表1。本文研究的復合材料在縱、橫向具有相同的力學性能，增強相為交織玻璃纖維織物，基體為WSR618型環氧樹脂。

表1 材料特性Table 1 Material properties

3．2 試驗結果與分析

雙向拉伸試驗在SDS100電液伺服動靜試驗機上進行。主要進行4個載荷比(f=縱向載荷∶橫向載荷=1∶1，2∶1，3∶1，4∶1)的雙向拉伸試驗，試驗采用位移控制方式，按f×0．5mm/min的加載速率加載。

同時，對十字型試樣的1/4模型，建立了有限元分析模型，獲得其有限元分析結果。表2給出了邊界元模型與有限元模型的比較。

表2 邊界元與有限元分析模型的比較Table 2 Comparison for the models of BEM and FEM

由表2可見，針對同一個分析模型，邊界元法在單元數目的劃分上要遠少于有限元法。因此，其模型的自由度數目也存在較大懸殊，邊界元法求解方程數目明顯較少，從而節約占機內存及計算時間。由此可見，邊界元法在計算效率上比有限元法更具優勢。

交織玻璃纖維增強復合材料的雙向拉伸縱向應力-應變曲線如圖6所示。從圖6中可看出，不同載荷比下，材料的雙向拉伸應力、應變保持著線性關系。邊界元法、有限元法計算結果與試驗結果趨勢一致，且較為接近。說明2種數值方法的模擬結果都是比較可信的。因此，可利用邊界元法，從理論上分析復合材料雙向拉伸力學性能。

圖6 材料雙向拉伸應力-應變曲線(縱向)Fig．6 Stress-strain curve of the material under biaxial tensile(longitudinal direction)

交織玻璃纖維增強復合材料不同載荷比的拉伸強度結果比較見表3。分析表3中數據可知，十字型試樣邊界元模型與有限元模型計算所得的雙向拉伸強度結果較接近;與試驗值比較，不同載荷比下2種計算方法的誤差范圍分別為5．9% ～6．8%和6．4% ～7．3%，且均為正誤差。可認為2種數值方法在計算精度上相當，2種方法的計算值均比試驗值偏大。誤差來源可能有2種:一是實際試樣存在一定的缺陷和細微損傷，造成拉伸強度的降低;二是理論分析模型與實際情況不完全吻合。

同時，由表3中結果可見，交織玻璃纖維增強復合材料的雙向拉伸行為存在明顯的雙向弱化現象:在載荷比為1∶1的加載條件下，材料的雙向拉伸強度(縱向強度)僅為單向拉伸強度的60．5%;載荷比為2∶1時，雙向拉伸強度為單向拉伸強度的71．5%;載荷比為3∶1和4∶1時，雙向拉伸強度分別為單向拉伸強度的77．5%和81．2%。由此可見，縱、橫向載荷比越小，該材料的雙向弱化效應越明顯，雙向加載為等雙拉時，雙向拉伸強度最小，材料抵抗雙向拉伸破壞的能力最弱。

表3 拉伸強度結果及比較Table 3 Results and comparison of tensile strength

十字型試樣雙向拉伸的破壞形式如圖7所示，破壞屬于脆性斷裂。

由圖7可見，不同載荷比下，交織玻璃纖維增強復合材料十字型試樣的破壞形式不同∶載荷比為1∶1時，試樣沿試驗區對角線斷裂;載荷比為2∶1和3∶1時，試樣裂紋產生及斷裂方向與試驗區對角線約成10°～15°夾角;載荷比為4∶1時，試樣出現明顯的分層破壞，且分層方向偏向于載荷大的加載方向。

4 結論

(1)建立了交織玻璃纖維增強復合材料雙向拉伸的邊界元理論分析模型，理論分析結果與試驗結果具有較好的一致性。因此，可利用邊界元法，從理論上分析該類復合材料的雙向拉伸性能。

圖7 十字型試樣破壞形貌Fig．7 Fracture morphology of cruciform specimen

(2)在計算交織玻璃纖維增強復合材料雙向拉伸強度時，邊界元法在計算效率上優于有限元法，計算精度上兩者相當;同時，2種方法的計算值均比試驗值偏大。

(3)交織玻璃纖維增強復合材料的雙向拉伸強度，與其單向拉伸強度相比，具有明顯的雙向弱化效應。縱、橫向載荷比越小，材料的雙向弱化現象越明顯，等雙拉時，材料的雙向拉伸強度(縱向強度)僅為單向拉伸強度的60．5%。

(4)不同載荷比下，交織玻璃纖維增強復合材料雙向拉伸的破壞形式有所不同，主要表現為試驗區纖維斷裂方向的不同，載荷比為4∶1時出現了分層破壞。

［1］ Lin W P，Hu H T．Parametric study on failure of fiber-reinforced composite laminates under biaxial tensile load［J］．Journal of Composite Materials，2002，36(12):1481-1503．

［2］ Smits A，Van Hemelrijck D，Philippidis T P，et al．Design of a cruciform specimen for biaxial testing of fibre reinforced composite laminates［J］．Composites Science and Technology，2006，66(7-8):964-975．

［3］ Makris A，Vandenbergh T，Ramault C，et al．Shape optimisation of a biaxially loaded cruciform specimen［J］．Polymer Testing，2010，29(2):216-223．

［4］ Serna Moreno M C，López Cela J J．Failure envelope under biaxial tensile loading for chopped glass-reinforced polyester composites［J］．Composites Science and Technology，2011，72(1):91-96．

［5］ Liu Y J．A new fast multipole boundary element method for solving large-scale two-dimensionalelastostatic problems［J］．International Journal Numerical Methods in Engineering，2006，65(6):863-881．

［6］ 陳龍，朱興一．求解具有弱界面的復合材料有效性能的近似模型與邊界元法［J］．復合材料學報，2013，30(2):137-143．

［7］ Lachat J C，Watson J O．Effective numerical treatment of boundary integral equations:A formulation for three‐dimensional elastostatics［J］．International Journal for Numerical Methods in Engineering，1976，10(5):991-1005．

［8］ Pineda León E，Samayoa Ochoa D，Rodríguez-Castellanos A，et al．Creeping analysis with variable temperature applying the boundary element method［J］．Engineering Analysis with Boundary Elements，2012，36(12):1715-1720．

［9］ 劉云忠，雷海峰．正交各向異性復合材料中裂紋問題的邊界元法分析［J］．固體火箭技術，1995，18(2):52-61．

［10］ 劉云忠．邊界元法計算厚板表面裂紋問題的應力強度因子［J］．固體火箭技術，1996，19(3):64-70．

［11］ 商欣萍，儲才元．用邊界元法分析試樣寬度對非織造土工織物拉伸性能的影響［J］．東華大學學報(自然科學版)，2004，30(1):1-5．

［12］ Araújo F C，d＇Azevedo E F，Gray L J．Boundary-element parallel-computing algorithm for the microstructural analysis of general composites［J］．Computers and Structures，2010，88(11):773-784．

［13］ Wang H T，Yao Z H．Large-scale thermal analysis of fiber composites using a line-inclusion model by the fast boundary element method［J］．Engineering Analysis with Boundary Elements，2013，37(2):319-326．

［14］ Fredholm I．Sur une classe d＇équations fonctionnelles［J］．Acta Mathematica，1903，27(1):365-390．

［15］ Rizzo F J，Shippy D J．A method for stress determination in plane anisotropic elastic bodies［J］．Journal of Composite Materials，1970，4(1):36-61．