早期复合材料的应用主要局限于薄壁非承力构件如夹层结构的表皮(壁厚3 mm左右),目前已逐步应用于承力构件且已从次承力构件转向主承力构件,例如一些大型飞机的复合材料厚壁结构(壁厚50 mm左右)^[1].复合材料厚壁管是一种典型的主承力构件,在航空和航天器结构中有着广泛的应用^[2],同时在民用方面也十分普遍,如跳高运动员撑竿^[3]和自行车支架^[4]等.作为主承力构件复合材料厚壁管在设计中更需要准确计算,但由于截面厚且曲率大,经典层合板壳理论^[5]及其他近似方法^[6]和实验研究^[7-8]难以满足精度要求.实际工程中由于层数众多且单层极薄,精度较高有限元等数值模拟中难以在厚度方向上进行精细建模^[9],因此需要采用弹性理论精确求解; 而大曲率复合材料层导致的高度各向异性材料结合较为复杂的弯曲荷载工况的理论分析又具有很大的挑战^[10-11].因此,进一步研究具有重要的经济价值和深远的社会意义.

复合材料结构的三维弹性力学理论求解主要包括Lekhnitskii^[12]的柔度法和Stroh^[13]的刚度法.Ting^[14]对两种方法作了较为详细的总结,Tarn等^[15]、Jolicoeur等^[16]、Zhang等^[17]和刘沛^[18]进一步发展了Lekhnitskii^[12]的双应力函数柔度法.对全由0°和90°特殊缠绕层以及各向同性材料层组合的特殊复合材料管,由于不存在面外应力而只需要单个应力函数进行求解^[12,19-20].但当有其他普通缠绕层与之相邻时,即便是特殊缠绕层和各向同性层也需要采用双应力函数进行求解^[21].对复合材料管的普通缠绕层可以采用统一联系参数法^[17]求解双应力函数,但当普通缠绕层趋向于0°或90°特殊缠绕层时需要分析统一参数的极限,由于该极限分析与具体复合材料的弹性系数导数有关而十分复杂,特别是对于各向同性材料层该极限分析中由于不存在缠绕角度而难以进行导数计算.

本文中采用独立双应力函数法直接求解特殊层,代替了统一参数法中当普通层趋向于特殊缠绕层的极限分析,避免了复杂的弹性系数求导,特别有利于与缠绕角无关而难以求导的各向同性材料层情况.

图1 纯弯曲载荷作用下的复合材料管
Fig.1 Composite tube subjected to pure bending loading

复合材料管是由多层纤维以不同角度缠绕而成的高度各向异性的复杂结构,纯弯曲荷载M_x和M_y作用下的复合材料管如图1所示,其中,b_n和b_n+1以及φ_n(n=1,2,…,N)为第n层的内外半径及其纤维缠绕角.柱坐标系下复合材料管内各层均具有以下形式的本构关系^[12]

{ε_r

ε_θ

ε_z

γ_θz

γ_rz

γ_rθ}=[C₁₁ C₁₂ C₁₃ C₁₄ 0 0

C₁₂ C₂₂ C₂₃ C₂₄ 0 0

C₁₃ C₂₃ C₃₃ C₃₄ 0 0

C₁₄ C₂₄ C₃₄ C₄₄ 0 0

0 0 0 0 C₅₅ C₅₆

0 0 0 0 C₅₆ C₆₆]{σ_r

σ_θ

σ_z

τ_θz

τ_rz

τ_rθ},(1)

式中,C_ij为圆柱坐标系下的弹性系数.复合材料管中的主要应力分量为沿着轴线的弯曲应力分量σ_z,主要受复合材料管整体弯曲变形影响,而次要应力分量包括平行于横截面的面内应力分量σ_r、σ_θ和τ_rθ以及垂直于横截面的面外应力分量τ_θz和τ_rz,仅与横截面的局部变形包括形状改变和翘曲有关(应力单元如图1所示).Lekhnitskii提出如下缩减弹性系数^[12]

β_ij=C_ij-(C_i3C_3j)/(C₃₃),(i,j=1,2,4,5,6).(2)

当缠绕角为0°或90°时材料主坐标系与圆柱坐标系重合,为圆柱型正交各向异性材料.在实际工程应用中,复合材料管可以根据需要采用数量不同的纤维材料和铺层角,从而形成性能各异的复合材料管,由于90°特殊缠绕角有利于复合材料管整体稳定,因而在工程中应用广泛.当缠绕角为0°或90°时

C₁₄=C₂₄=C₃₄=C₅₆=0.(3)

把式(3)带入式(2)中可得该特殊层的缩减弹性系数如下:

β₁₄=β₂₄=β₅₆=0.(4)

Lekhnitskii^[12]首先提出了包含两个应力函数的柔度法,其中,对于复合材料管纯弯曲问题,将面内应力分量表示成应力函数F(r,θ)的如下形式:

σ_r=1/r(F)/(r)+1/(r²)(²F)/(θ²),σ_θ=(²F)/(r²),

τ_rθ=-(²)/(rθ)(F/r).(5)

同时,将面外应力分量表示成应力函数Ψ(r,θ)的如下形式:

τ_rz=1/r(Ψ)/(θ),τ_θz=-(Ψ)/(r).(6)

复合材料管纯弯曲的三维弹性力学问题是一个广义平面应变问题,其弯曲应力可以表示成如下形式^[12]

σ_z=1/(C₃₃)(κ_xr sin θ-κ_yrcos θ-C₁₃σ_r-

C₂₃σ_θ-C₃₄τ_θz).(7)

由式(7)可见,弯曲应力与曲率κ_x和κ_y以及面内应力式(5)和面外应力式(6)有关,但体现整体弯曲变形的轴线曲率是主要的,而与横截面形状改变和翘曲的局部变形有关的面内应力和面外应力是次要的.在纯弯曲复合材料管分析中可以进一步采用变量分离方法将双应力函数假设为如下形式

{F(r,θ)=(κ_xsin θ-κ_ycos θ)f(r),

Ψ(r,θ)=(κ_xsin θ-κ_ycos θ)ψ(r).(8)

将式(8)带入相容方程^[12]可得到以下两个常微分方程

{ ((β₂₂)/r(r(d⁴)/(dr⁴)+2(d³)/(dr³))-(β₁₁+2β₁₂+β₆₆)/(r⁴)

(r²(d²)/(dr²)-r(d)/(dr)+1))f(r)+((β₁₄)/(r²)((d)/(dr)+

r(d²)/(dr²))-(β₂₄)/r(r(d³)/(dr³)+2(d²)/(dr²))+(β₅₆)/(r²)(d)/(dr))

ψ(r)=2/r(C₁₃-C₂₃)/(C₃₃),

((β₁₄)/(r³)(r²(d²)/(dr²)-r(d)/(dr)-1)+(β₂₄)/r(-r(d³)/(dr³)+

(d²)/(dr²))+(β₅₆)/(r³)(-r(d)/(dr)-1))f(r)+

((β₄₄)/r(r(d²)/(dr²)+(d)/(dr))-(β₅₅)/(r²))ψ(r)=(2C₃₄)/(C₃₃).(9)

对复合材料组合管纯弯曲的求解可以将联立的相容方程组(9)转化为更加高阶的微分方程^[12],但在柱坐标系下比较困难,为此,Jolicoeur等^[16]通过引进双应力函数之间的联系参数建立起它们之间的桥梁,但存在无穷大参数,Zhang等^[17]和刘沛^[18]对联系参数法^[16]进行了初步改进,从而提出了统一联系参数法:

{ f(r)=∑²_i=1(K_i)/(m_i)r^m_i+1+∑⁴_i=3(K^*_i)/(m_i)g^*_i r^m_i+1+

1/2μ₁r³,

ψ(r)=∑²_i=1K_ig_ir^m_i+∑⁴_i=3K^*_i r^m_i+μ₂r².(10)

式中:K_i和K^*_i为组合系数; μ₁和μ₂为与式(9)两个特解有关的参数; m_i(i=1,2,3,4)为与式(2)中缩减弹性系数有关的特征参数; g₁,g₂,g^*₃,g^*₄ 为特征参数m_i(i=1,2,3,4)所对应的统一联系参数; K₁和K₂以及K^*₃和K^*₄为组合系数.

Lekhnitskii^[12]在求解复合材料管纯弯曲问题时只考虑90°特殊缠绕单层的情况,由于没有面外剪切应力,事实上属于单应力函数问题.Zhang等^[20]考虑了全部由特殊缠绕层组合的情况如[0°/90°],但也属于单应力函数问题.Xia等^[19]分析了夹层复合材料管纯弯曲问题,由于把交叉缠绕面层等价为90°缠绕层,所以求解的是一个包含等价的90°特殊缠绕层和各向同性芯层的夹层组合管,由于没有面外剪切应力,从而只需要单应力函数即可.然而,并非所有特殊缠绕层和各向同性层都不存在面外切向应力而仅需要单应力函数,事实上当有其他缠绕层与之相邻时,不仅有面内应力而且有面外应力,所以,即便是特殊缠绕层和各向同性层,也需要引进双应力函数并采用统一参数法进行求解.但是,在统一参数法中当普通缠绕层趋向于特殊缠绕层时,式(10)中的联系参数g^*₃ 与g^*₄ 为0/0不定型,可以利用洛必塔法则,通过材料参数关于缠绕角的导数求解它们的极限,但比较复杂; 特别是对于与普通缠绕层相邻的各向同性层,由于仍然存在面外剪切应力而必须采用双应力函数求解,而该各向同性层又不存在任何缠绕角度,所以难以采用关于缠绕角导数来求解统一联系参数的极限.为此,需要进一步对统一联系参数法进行研究.

对于0°和90°特殊缠绕层以及各向同性材料,把式(4)代入式(9)可得

{ ((β₂₂)/r(r(d⁴)/(dr⁴)+2(d³)/(dr³))-(β₁₁+2β₁₂+β₆₆)/(r⁴)

(r²(d²)/(dr²)-r(d)/(dr)+1))f(r)=2/r(C₁₃-C₂₃)/(C₃₃),

((β₄₄)/r(r(d²)/(dr²)+(d)/(dr))-(β₅₅)/(r²))ψ(r)=0.(11)

由式(11)可见双应力函数f(r)和ψ(r)互相独立,所以,虽然与普通缠绕层比邻的特殊缠绕层,特别是各向同性层,由于仍然存在面外剪切应力而必须采用双应力函数求解,但可以独立求解这两个应力函数而不需引进它们之间的任何联系参数.假设f(r)和ψ(r)具有以下形式的解

f(r)=K^Ir^mI+1,ψ(r)=K^IIr^mII,(12)

代入式(11)可得两个独立特征方程

{(m^I)²(β₂₂(m^I)²-(β₂₂+β₁₁+2β₁₂+

β₆₆))=0,

β₄₄(m^II)²-β₅₅=0.(13)

由式(13)求得零特征根之外的如下特征根:

{m^I₁=((β₁₁+2β₁₂+β₂₂+β₆₆)/(β₂₂))^1/2,

m^I₂=-((β₁₁+2β₁₂+β₂₂+β₆₆)/(β₂₂))^1/2;

{m^II₁=((β₅₅)/(β₄₄))^1/2,

m^II₂=-((β₅₅)/(β₄₄))^1/2.(14)

容易验证微分方程(11)中的第1式的特解可表示如下:

f₁(r)=1/2μ₁r²,(15)

式中

μ₁=1/(-β₁₁-2β₁₂+3β₂₂-β₆₆)(C₁₃-C₂₃)/(C₃₃).(16)

于是微分方程(11)的通解为:

f(r)=∑²_i=1(K^I_i)/(m^I_i )r^mI_i+1+1/2μ₁r³,

ψ(r)=∑²_i=1K^II_i r^mII_i.(17)

显然,在统一联系参数法中当特殊缠绕层特别是各向同性层时,直接在式(10)中令

g₁=g₂=g^*₃=g^*₄=μ₂=0.(18)

即可得到与式(17)完全等价的结果,换言之,式(17)就是统一联系参数法中当普通缠绕层趋向于特殊缠绕层的极限,可以在统一联系参数法中代替该极限,从而避免了与弹性系数复杂导数有关的极限分析过程,这对于与缠绕角无关从而难以求导的各向同性层尤其重要,推广了统一参数法的应用范围,是统一参数法的一个重要改进.

利用式(17)可以由式(5)求得应力分量,再利用所求应力分量由式(1)得到应变,并进一步利用几何方程得到关于位移的微分方程,最后对位移微分方程进行积分可以求得位移表达式.在此基础上,采用完全类似于常规联系参数法^[16]或者统一联系参数法^[17-18]的求解步骤,可以由复合材料管表面自由边界条件和层间应力和位移连续性得到每个普通缠绕层的K_1,n,K_2,n和K^*_3,n,K^*_4,n 或者特殊层的K^I_1,n,K^I_2,n 和K^II_1,n,K^II_2,n,其中n=1,2,…,N.最终,利用弯曲应力式(7)和杆端载荷条件求得复合材料管轴线的曲率如下:

κ_x=(M_x)/((EI)),κ_y=(M_y)/((EI)).(19)

其中,(EI)为复合材料管等效抗弯刚度,可以表示成

(EI)=∑^N_n=1(EI)_n.(20)

表1 复合材料(Carbon/Epoxy)的材料参数
Tab.1 Composite material(Carbon/Epoxy)properties

式中对普通缠绕层的等效抗弯刚度^[17-18]表示如下

(EI)_n=(π)/(C_33.n){∑²_i=1K_i.n[C_13,n+C_23,n(m_i,n+1)-

C_34,nm_i,ng_i,n](b^m_i,n+2_n-b^m_i,n+2_n+1)/(m_i,n+2)+

∑⁴_j=3K^*_j,n[C_13,ng^*_j,n+C_23,n(m_j,n+1)g^*_j,n-

C_34,nm_j,n](b^m_j,n+2_n-b^m_j,n+2_n+1)/(m_j,n+2)+[μ_1,n(C_13,n+

3C_23,n)-2μ_2,nC_34,n-1](b⁴_n-b⁴_n+1)/4}.(21)

而对于特殊缠绕层的等效抗弯刚度表示如下:

(EI)_n=(π)/(C_33,n){∑²_i=1K^I_i,n[C_13,n+C_23,n(m^I_i,n+1)]

(b^{mI_i,n+2_n-bmI_i,n+2_n+1)/(mI_i,n+2)+[μ_1,n(C_13,n+3C_23,n)-1]}

(b⁴_n-b⁴_n+1)/4}.(22)

特别对于复合材料组合管中的各向同性层,由于

C₁₁=C₂₂=C₃₃=1/E,C₁₂=C₁₃=C₂₃=-ν/E,

C₄₄=C₅₅=C₆₆=(2(1+ν))/E,(23)

式中E和ν分别是杨氏弹性模量和泊松比,把式(23)代入式(2)可得

β₁₁=β₂₂=(1-ν²)/E,β₁₂=β₂₁=-(ν(1+ν))/E,

β₄₄=β₅₅=β₆₆=(2(1+ν))/E.(24)

把式(24)代入式(14)和(16)可得

{m^I₁=2,

m^I₂=-2, {m^II₁=1,

m^II₂=-1,μ₁=0.(25)

因此,各向同性层等效抗弯刚度可表示如下

(EI)_n=(4K^I_1,n ν_n+1)EI_n,(26)

式中

I_n=(π(b⁴_n+1-b⁴_n))/4,(27)

由式(26)可见,由于通常K^I_1,n≠0,所以(EI)_n≠EI_n,说明与材料力学中理论所得结果略有不同.

考虑纯弯曲荷载作用下的8根简单复合材料组合管,左端固支而右端作用线性分布荷载,等效弯矩为M_x=8 kN·m和M_y=0.复合材料组合管的内部为钢管而外部缠绕复合材料,钢管的弹性模量和泊松比分别为E=210 GPa和ν=0.25,采用表1所示复合材料(Carbon/Epoxy),长度为600 mm,铺层顺序及其几何参数见图2,其中,复合材料组合管[Steel/90°]只包含90°特殊缠绕层和各向同性材料层,而[Steel/25°]和[Steel/±25°]在特殊的各向同性材料层之外还包含有一般的25°缠绕层和±25°交叉缠绕层.采用本文中改进的统一联系参数法进行计算,并与精细有限元数值模拟结果进行对比,其中精细有限元采用三维六面体单元,以复合材料组合管[Steel/90°]为例,有限元模型如图3所示,包括单元个数24 000,节点个数28 633,本文的理论解析求解和有限元数值求解的应力结果对比如图4所示.

图2 复合材料组合管横截面几何参数
Fig.2 Geometrical dimension of cross section for composite tubes

图3 复合材料组合管有限元模型
Fig.3 FEM model for composite tubes

图4 复合材料组合管的典型应力规律
Fig.4 Typical stresses for composite tubes

通过复合材料组合管应力分量随半径的变化曲线可见,本文中方法与精细有限元数值模拟结果吻合较好,验证了本文中方法是准确和可靠的.从计算结果可见:复合材料组合管[Steel/90°]由于不包含其他普通缠绕层,面外剪切应力为零τ_θz=τ_rz=0,而其他复合材料组合管由于同时包含了普通层,使得与之相邻的各向同性层(钢管)也存在面外剪切应力τ_θz≠0和τ_rz≠0.但复合材料组合管[Steel/±25°]由于交叉缠绕层在一定程度上可以近似等价为90°特殊缠绕层,该钢管面外剪切应力虽不为零但比较小.

在上述算例中需要指出:1)因为特殊层外紧邻有普通缠绕层,所以必须采用双应力函数求解,以往的单应力函数法^[12,19-20]难以求解; 2)由于包含了各向同性材料层,因此,以往的统一参数法^[17-18]也难以有效处理; 3)本文中虽然对层数较少且单层相对厚度较大的简单复合材料管采用有限元进行精细分析以验证本文中理论方法的准确性和可靠性,但一般实际工程中的复合材料管由于层数众多且单层相对层厚极小所以很难在厚度方向上采用有限元进行精细建模,因此通常需要采用经典层合板理论或有限元分析中的层合单元进行近似; 然而,一旦通过简单复合材料管验证了本文中理论方法的准确性和可靠性,则可以采用这种理论方法去准确分析层数众多且单层相对层厚极小的实际工程复合材料管了,而这种实际复合材料管是有限元等数值方法难以精细模拟的.

本文中采用独立双应力函数求解复合材料管中0°和90°特殊缠绕层以及各向同性层,由于可以代替统一参数的极限分析,避免了弹性系数关于缠绕角的复杂导数,特别有利于不存在缠绕角而难以计算导数的各向同性材料层情况,从而很好地改进了传统的统一参数法.