层合结构是由不同材料属性的薄片通过某种工艺(如粘接)复合而成。现有的层合结构包括纤维增强复合材料^[1-2]、三明治夹层结构^[3-4]和钛合金层合结构^[5]等多种形式，具有参数可设计性强、性价比优等特点，被广泛应用于航空航天领域^[6-7]。在层合结构中，通常不同材料之间连接界面的力学性能最为薄弱。过大的层间应力会导致层合结构损伤破坏，成为其失效的主要形式^[8-9]。如纤维增强复合材料损伤失效有50%以上是分层引起的^[10]。因此，准确预测层合结构的层间应力尤为重要。

层合梁是航空航天领域典型的承力构件，已经发展了多种分析模型，用于计算其横向剪应力，从而获得层间应力。现有的分析模型主要分为3类：等效单层理论^[11-12]、分层理论^[13-14]和锯齿理论^[15-18]。等效单层理论的计算效率较高，然而由于无法满足层间剪应力连续条件，计算误差较大。如Timoshenko梁理论预测的三明治夹层梁的横向剪应力，计算误差高达351%^[18]。理论上，分层理论能够准确预测复合材料层合结构的层间剪应力，但所需的计算量十分庞大。如300层的层合板，每个节点自由度数高达1 803个^[19]。锯齿理论是在整体高阶理论的基础上添加合适的局部锯齿函数构造而成的。通常整体理论不超过三阶，而锯齿函数能够描述面内位移沿厚度方向的锯齿^[20]变化。显然，如果选择合适的锯齿函数，锯齿理论只需较少的未知变量就能够准确计算层合结构的横向剪应力，从而很好地兼顾计算效率和计算精度，成为当前研究的主要方向。

锯齿理论最早由Lekhnitskii^[21]于1935年提出，提出后并未引起重视，直到1986年Murakami锯齿理论^[15]和Di Sciuva锯齿理论^[22]的提出，锯齿理论才逐渐被认可。由于Di Sciuva锯齿理论能够预先满足层间应力连续条件，得到快速发展。基于Di Sciuva锯齿理论，Cho和Oh^[23]提出了一种三阶锯齿模型，并将之用于预测承受力、热和电载荷的复合材料厚板的变形和应力。早期的锯齿理论^{[15-17, 22-23]}仅适用于预测层数较少的层合梁和材料属性差异不大的三明治夹层梁的横向剪应力。2015年，Tessler等^[18]基于经典的修正锯齿理论(RZT)，提出了一种混合修正锯齿理论(RZTM)^[24-25]。贺丹和杨万里^[26]应用该理论分析了层合梁的弯曲和自由振动问题。该理论能够准确计算出材料属性差异较大的夹层梁的横向剪应力。然而对于层数较多的层合梁，预测精度较差。常见的T300/QY8911复合材料的单层厚度为0.12 mm^[27]，在实际的航空结构中，超过25层的纤维增强复合材料(即厚度超过3 mm)十分常见，如机翼蒙皮^[28]。因此，发展能够准确预测较多层数的层合结构横向剪应力的锯齿理论是非常必要的。基于整体高阶锯齿理论^[29]，Wu和Chen^[30]提出了一种混合整体高阶锯齿理论(GHZTM)，能够准确预测大多数的多层厚梁的横向剪应力，然而个别多层梁横向剪应力预测精度稍差(如3层层合梁^[30])。另外，该理论的挠度场要求满足C¹连续，为C¹型锯齿理论，不利于梁单元的构造。早期的锯齿理论大都是C¹型锯齿理论。2011年，Ren等^[31]提出了一种C⁰型锯齿理论的构造方法，并提出了一种C⁰板单元。由于C⁰型锯齿理论具有计算效率高，容易在商业软件实现等优点，而得到快速发展。随后，许多C⁰型锯齿理论^[32-34]和C⁰单元^[35-38]被提出。Wu等^[39]比较了大量的C⁰和C¹型锯齿理论模型，发现相比于C¹型锯齿理论C⁰型锯齿理论不仅有更高的计算效率还具有更高的计算精度。

本文提出了一个新的线性分段锯齿函数，采用Ren等^[31]提出的C⁰型锯齿理论构造方法，建立了一种面向层合梁结构的C⁰型新锯齿理论模型。该模型不仅能够准确预测层数较多的纤维增强复合材料层合梁的横向剪应力，同时能够准确预测芯层与上下面板材料属性差异较大的三明治夹层梁的横向剪应力。

1 锯齿理论描述

锯齿理论的核心在于锯齿函数的构造。评价锯齿理论的优劣，主要看其锯齿函数能否准确地描述面内位移沿厚度方向的C_z^0[20]变化(沿厚度z方向满足C⁰连续条件)。根据锯齿函数的不同，众多的锯齿理论被提出^[15-18]。下面通过目前应用较为广泛的2种锯齿理论(Murakami型锯齿理论(MZZT)^[15]和修正锯齿理论(RZT)^[18])的介绍，论述锯齿理论的基本思想和存在的问题。

锯齿层合梁理论的位移场可统一表示为

$ \left\{ {\begin{array}{*{20}{l}} {{u^k}(x,z) = {u_0}(x) + \sum\limits_{i = 1}^j {{z^i}} {u_i}(x) + {u_z}(x,z)}\\ {w = {w_0}(x)} \end{array}} \right. $	（1）

式中：z为厚度坐标；k表示第k层；u₀和w₀为中面位移；u₁为中面法线方向绕y轴的转角, u_i为泰勒级数展开的高阶项；j表示整体部分的阶次，如j=1表示整体一阶，j=3表示整体三阶；u_z(x)为锯齿函数。

MZZT锯齿理论^[15]的锯齿函数可写为

$ {u_z}(x,z) = {( - 1)^k}{\zeta _k}\psi (x) $	（2）

式中：ζ_k为每层z方向的局部坐标，ζ_k=a_kz－b_k, ζ_k∈(－1, 1)，a_k=2/(z_k+1－z_k)，b_k=(z_k+1+z_k)/(z_k+1－z_k)；ψ(x)为锯齿转角，度量锯齿函数对沿厚度分布的轴向位移的影响程度。对于任意x=x_a，MZZT锯齿函数沿厚度的变化如图 1所示^[15]。从图 1和式(2)可知，MZZT锯齿函数在层和层交界面处的值，只能从－ψ(x_a)和ψ(x_a)中选取。此外，该锯齿函数不依赖于每层材料的性质，无法描述每层材料属性的变化，导致该理论无法准确预测横向各向异性的层合梁的变形和应力^[40]。尤其对于面板软夹层硬的夹层结构，其预测精度较差^[41]。

RZT锯齿理论^[18]的锯齿函数可写为

$ {u_z}(x,z) = \left( {N_1^k\left( {{\zeta _k}} \right){\phi _k} + N_2^k\left( {{\zeta _k}} \right){\phi _{k + 1}}} \right)\psi (x) $	（3）

式中：N₁^k和N₂^k是线性插值函数：N₁^k(ζ_k)=(1－ζ_k)/2, N₂^k(ζ_k)=(1+ζ_k)/2；ϕ_k为界面位移。对于任意x=x_a，RZT锯齿函数沿厚度的变化如图 2^[18]所示。相比于MZZT锯齿函数，该锯齿函数在层和层交界面处可取任意值。由于ϕ_k依赖于每层的材料属性，因此RZT锯齿函数能够描述每层材料的变化，能够准确预测横向各向异性的层合梁的变形和应力。图 2中给出了层合梁不同横截面处(x_a和x_b处)的锯齿函数沿厚度方向的变化情况。从图 2中可以看出，RZT锯齿函数在不同横截面处具有相同的变化规律，区别只是幅值不同。因此，RZT锯齿理论不适用于计算不同横截面处面内位移沿厚度方向变化规律相差较大的情形。Wu和Chen^[30]指出随着层合梁层数的增加，RZT锯齿理论预测精度会逐渐降低。因此需要发展一种既能够描述每层材料变化，又能够描述不同横截面处变化规律不同的新锯齿函数。构造一种既能够准确预测材料属性差异较大夹层梁的变形和应力，又能够准确预测层数较多的层合梁的变形和应力的新锯齿理论模型。

2 新锯齿理论模型 2.1 新锯齿函数

新锯齿函数为

$ {u_z}(x,z) = N_1^k\left( {{\zeta _k}} \right){{\bar u}_k}(x) + N_2^k\left( {{\zeta _k}} \right){{\bar u}_{k + 1}}(x) $	（4）

式中：u_k(x)为层与层交界面上的修正位移，与每层的材料性质相关。图 3给出了层合梁不同横截面处(x_a和x_b处)的新锯齿函数沿厚度方向的变化情况。对比图 3和图 2可知，新锯齿函数能够描述不同横截面处不同的变化规律。综上，本文提出的新锯齿函数既能描述每层材料变化，又能描述不同横截面处变化规律的不同。

2.2 位移场

位移场的整体部分仍取三阶，因此，构造的初始位移场为

$ \left\{ \begin{array}{l} {u^k}(x,z) = {u_0}(x) + \sum\limits_{i = 1}^3 {{z^i}} {u_i}(x) + \\ \;\;\;\;\;N_1^k\left( {{\zeta _k}} \right){{\bar u}_k}(x) + N_2^k\left( {{\zeta _k}} \right){{\bar u}_{k + 1}}(x)\\ w = {w_0}(x) \end{array} \right. $	（5）

式(5)中共含有n+6(n为层数)个未知变量，并且未知量个数与层数相关。下面推导的目的是消去部分未知量，得到最终位移场。

首先，令u₁(x)=u_n+1(x)=0，可消去2个未知量，还剩余n+4个未知量。

其次，使用横向剪应力层间连续条件

$ \tau _{xz}^{k + 1}\left( {x,{z_k}} \right) = \tau _{xz}^k\left( {x,{z_k}} \right) $	（6）

和横向剪应力自由表面条件

$ \tau _{xz}^1\left( {x,{z_1}} \right) = \tau _{xz}^n\left( {x,{z_{n + 1}}} \right) = 0 $	（7）

可得n+1个方程可以消去n+1个未知量。消去u₂、w_{0, x}、u₂~u_k。n+1个方程可写成如下形式：

$ \mathit{\boldsymbol{Ku}} = \mathit{\boldsymbol{A}}{u_1} + \mathit{\boldsymbol{B}}{u_3} $	（8）

式中：u=[u₂   w_{0, x}   u₂   …   u_n]^T; K、A和B是系数矩阵。u可表示为

$ \boldsymbol{u}=\boldsymbol{K}^{-1} \boldsymbol{A} u_{1}+\boldsymbol{K}^{-1} \boldsymbol{B} u_{3} $	（9）

式(9)可重写为

$ \boldsymbol{u}=\boldsymbol{C}_{1} u_{1}+\boldsymbol{C}_{3} u_{3} $	（10）

式中：C₁=K^－1A，C₃=K^－1B。

将式(10)代入式(5)中，高阶锯齿理论最终位移场可写为

$ \left\{ {\begin{array}{*{20}{l}} {{u^k}(x,z) = {u_0}(x) + \mathit{\varPhi }_1^k(z){u_1}(x) + \mathit{\varPhi }_2^k(z){u_3}(x)}\\ {w = {w_0}(x)} \end{array}} \right. $	（11）

当k=1时：

$ \mathit{\varPhi }_1^1(z) = z + {\mathit{\boldsymbol{C}}_1}\left( 1 \right){z^2} + {\mathit{\boldsymbol{C}}_1}\left( 3 \right)N_2^1 $

$ \mathit{\varPhi }_2^1(z) = {\mathit{\boldsymbol{C}}_3}(1){z^2} + {z^3} + {\mathit{\boldsymbol{C}}_3}(3)N_2^1 $

当k=2~n-1时：

$ \begin{array}{l} \mathit{\varPhi }_1^k(z) = z + {\mathit{\boldsymbol{C}}_1}(1){z^2} + {\mathit{\boldsymbol{C}}_1}(k + 1)N_1^k + \\ \;\;\;\;\;\;\;\;{\mathit{\boldsymbol{C}}_1}(k + 2)N_2^k \end{array} $

$ \begin{array}{l} \mathit{\varPhi }_2^k(z) = {\mathit{\boldsymbol{C}}_3}(1){z^2} + {z^3} + {\mathit{\boldsymbol{C}}_3}(k + 1)N_1^k + \\ \;\;\;\;\;\;\;\;{\mathit{\boldsymbol{C}}_{\rm{3}}}(k + 2)N_2^k \end{array} $

当k=n时：

$ \mathit{\varPhi }_1^n(z) = z + {\mathit{\boldsymbol{C}}_1}(1){z^2} + {\mathit{\boldsymbol{C}}_1}(n + 1)N_1^n $

$ \mathit{\varPhi }_2^n(z) = {\mathit{\boldsymbol{C}}_3}(1){z^2} + {z^3} + {\mathit{\boldsymbol{C}}_3}(n + 1)N_1^n $

该理论仅含有4个与层数无关的未知量。不含有横向位移的一阶导数，构造有限单元时，仅需满足C⁰连续。

2.3 几何方程和本构方程

对于线弹性问题，层合梁的几何方程为

$ \varepsilon _x^k = \mathit{\boldsymbol{B}}_x^k{\mathit{\boldsymbol{\delta }}_x},\gamma _{xz}^k = \mathit{\boldsymbol{B}}_{xz}^k{\mathit{\boldsymbol{\delta }}_{xz}} $	（12）

式中：ε_x^k为第k层沿x方向的正应变；γ_xz^k为第k层在垂直于x轴的平面上沿z方向的剪应变；B_x^k=[1   Φ₁^k  Φ₂^k]，δ_x=[u_{0, x}  u_{1, x}  u_{3, x}]^T，B_xz^k=[Φ_{1, z}^k   Φ_{2, z}^k  1]，δ_xz=[u₁  u₃  w_{0, x}]^T。

第k层层合梁的本构方程为

$ \left[ {\begin{array}{*{20}{l}} {\sigma _x^k}\\ {\tau _{xz}^k} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{E^k}}&0\\ 0&{{G^k}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\varepsilon _x^k}\\ {\gamma _{xz}^k} \end{array}} \right] $	（13）

式中：σ_x^k和τ_xz^k分别为对应的正应力和剪应力；E^k和G^k分别为第k层的杨氏模量和剪切模量。

2.4 横向剪应力的预处理

为了获得更为精确的横向剪应力，Tessler^[24]提出了一种通过局部平衡方程获得横向剪应力的方法。忽略体力，第k层局部平衡方程可表示为

$ \sigma _{x,x}^k + \tau _{xz,z}^k = 0 $	（14）

将式(14)沿z向积分，并强迫横向剪应力满足层间连续条件。横向剪应力可重写为

$ \bar \tau _{xz}^k(x,z) = - {T_{\rm{b}}}(x) - \int_{{z_1}}^z {\sigma _{x,x}^k} {\rm{d}}z $	（15）

式中：T_b为下表面的横向剪应力。

将式(13)代入式(15)中，可得

$ \bar \tau _{xz}^k(x,z) = - {T_{\rm{b}}}(x) - {\mathit{\boldsymbol{C}}_z}{\mathit{\boldsymbol{U}}_x} $	（16）

式中：U_x=[u_{0, xx}   u_{1, xx}  u_{3, xx}]^T；C_z=$\left[\int_{z_{1}}^{z} E^{k} \mathrm{d} z \quad \int_{z_{1}}^{z} \varPhi_{1}^{k} E^{k} \mathrm{d} z \quad \int_{z_{1}}^{z} \varPhi_{2}^{k} E^{k} \mathrm{d} z\right]$。

当z=z_n+1时，有

$ {T_{\rm{t}}}(x) = - {T_{\rm{b}}}(x) - \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{U}}_x} $	（17）

式中：T_t为上表面的横向剪应力；

$ \mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {\int_{{z_1}}^{{z_{n + 1}}} {{E^k}} {\rm{d}}z}&{\int_{{z_1}}^{{z_{n + 1}}} {\mathit{\varPhi }_1^k} {E^k}{\rm{d}}z}&{\int_{{z_1}}^{{z_{n + 1}}} {\mathit{\varPhi }_2^k} {E^k}{\rm{d}}z} \end{array}} \right]。$

由横向剪应力自由表面条件T_t=T_b=0，消去面内位移的二阶导可得

$ {u_{0,xx}} = {\overline {\mathit{\boldsymbol{CU}}} _x} $	（18）

式中：$\overline{\boldsymbol{U}}_{x}=\left[\begin{array}{ll}{u_{1, x x}} & {u_{3, x x}}\end{array}\right]^{\mathrm{T}} ;\overline{\boldsymbol{C}}=\left[-\frac{\boldsymbol{C}(2)}{\boldsymbol{C}(1)}-\frac{\boldsymbol{C}(3)}{\boldsymbol{C}(1)}\right]$。

将式(18)代入式(16)中，可得横向剪应力表达式为

$ \bar \tau _{xz}^k(x,z) = {\mathit{\boldsymbol{\overline C}} _z}{\mathit{\boldsymbol{\overline U}} _x} $	（19）

式中：

$ \begin{array}{l} {\mathit{\boldsymbol{\overline C}} _z} = \\ \;\;\;\;\;\;\;\left[ { - \frac{{\mathit{\boldsymbol{C}}(2){\mathit{\boldsymbol{C}}_z}(1)}}{{\mathit{\boldsymbol{C}}(1)}} - {\mathit{\boldsymbol{C}}_z}(2) - \frac{{\mathit{\boldsymbol{C}}(3){\mathit{\boldsymbol{C}}_z}(1)}}{{\mathit{\boldsymbol{C}}(1)}} - {\mathit{\boldsymbol{C}}_z}(3)} \right] \end{array} $

2.5 Reissner混合变分原理和平衡方程

Reissner混合变分原理假定位移和横向剪应力和横向法应力是相互独立。Reissner混合变分原理可以看作Hellinger-Reissner变分原理^[42]的特例，可由Hellinger-Reissner变分原理退化得到。Hellinger-Reissner变分原理的总势能函数可表示为

$ \mathit{\Pi }\left( {{u_i},{\varepsilon _{ij}}} \right) = \int_V {\left( {W\left( {{\varepsilon _{ij}}} \right) - {f_i}{u_i}} \right)} {\rm{d}}V - \int_{{S_\sigma }} {{{\bar t}_i}} {u_i}{\rm{d}}S $	（20）

式中：u_i和ε_ij分别为位移分量和应变张量；W(ε_ij)为应变能密度；f_i为体力分量；S_σ为已知边界力(t_i)的边界；V为体积；S为表面积。

应用拉格朗日乘子法修正式(20)可得

$ \begin{array}{l} \mathit{\bar \Pi }\left( {{u_i},{\varepsilon _{ij}},{\sigma _{ij}}} \right) = \int_V {\left( {W\left( {{\varepsilon _{ij}}} \right) - {f_i}{u_i}} \right)} {\rm{d}}V + \\ \;\;\;\;\;\;\;\;\int_V {{\sigma _{ij}}} \left[ {{\varepsilon _{ij}} - \frac{1}{2}\left( {{u_{i,j}} + {u_{j,i}}} \right)} \right]{\rm{d}}V - \\ \;\;\;\;\;\;\;\;\int_{{S_\sigma }} {{{\bar t}_i}} {u_i}{\rm{d}}S - \int_{{S_u}} {{t_i}} \left( {{u_i} - {{\bar u}_i}} \right){\rm{d}}S \end{array} $

（21）

式中：σ_ij为应力张量；S_u为已知边界位移(u_i)的边界。

根据Reissner的假定，应变和应力满足本构关系

$ {\varepsilon _{ij}} = \frac{{\partial {W_{\rm{c}}}\left( {{\sigma _{ij}}} \right)}}{{\partial {\sigma _{ij}}}} $	（22）

式中：W_c(σ_ij)为应变余能密度。

将式(22)代入式(21)中，可得Reissner混合变分原理的总势能函数

$ \begin{array}{l} {\mathit{\Pi }_R}\left( {{u_i},{\sigma _{ij}}} \right) = \\ \;\;\;\;\;\;\;\;\int_V {\left[ {\frac{1}{2}\left( {{u_{i,j}} + {u_{j,i}}} \right){\sigma _{ij}} - {W_{\rm{c}}}\left( {{\sigma _{ij}}} \right) - {f_i}{u_i}} \right]{\rm{d}}V} - \\ \;\;\;\;\;\;\;\;\int_{{S_\sigma }} {{{\hat t}_i}} {u_i}{\rm{d}}S - \int_{{S_u}} {{t_i}} \left( {{u_i} - {{\hat u}_i}} \right){\rm{d}}S \end{array} $

（23）

根据式(22)，横向切应变和横向正应变可以写成式(24)的形式^[43]

$ \varepsilon _{\alpha z}^{\rm{a}} = \frac{{\partial {W_{\rm{c}}}\left( {{\sigma _{ij}}} \right)}}{{\partial {\tau _{az}}}} = \frac{1}{2}\gamma _{\alpha z}^{\rm{a}},\varepsilon _{zz}^{\rm{a}} = \frac{{\partial {W_{\rm{c}}}\left( {{\sigma _{ij}}} \right)}}{{\partial {\sigma _{zz}}}} $	（24）

式中：α=x, y；上标a表示由应变余能密度得到的横向应变。

将式(24)代入式(23)，并对式(23)进行变分可得

$ \begin{array}{l} \delta {\mathit{\Pi }_R} = \int_V {\left[ {\left( {{\sigma _{\alpha \beta }}{\rm{ \mathit{ δ} }}{\varepsilon _{\alpha \beta }} + \tau _{\alpha z}^{\rm{a}}{\rm{ \mathit{ δ} }}{\gamma _{\alpha z}} + \sigma _{zz}^{\rm{a}}{\rm{ \mathit{ δ} }}{\varepsilon _{zz}}} \right) + } \right.} \\ \;\;\;\;\;\;\;\;\left. {{\rm{ \mathit{ δ} }}\tau _{\alpha z}^{\rm{a}}\left( {{\gamma _{\alpha z}} - \gamma _{\alpha z}^{\rm{a}}} \right) + {\rm{ \mathit{ δ} }}\sigma _{zz}^{\rm{a}}\left( {{\varepsilon _{zz}} - \varepsilon _{zz}^{\rm{a}}} \right)} \right]{\rm{d}}V - \\ \;\;\;\;\;\;\;\;\int_{{S_\sigma }} {{{\hat t}_i}} {u_i}{\rm{d}}S - \int_{{S_u}} {{t_i}} \left( {{u_i} - {{\hat u}_i}} \right){\rm{d}}S \end{array} $

（25）

式中：β=x, y；τ_αz^a和σ_zz^a分别为假定的横向剪应力和横向正应力。

对于层合梁，忽略正应力，Reissner混合变分原理可写为

$ \begin{array}{l} \int_\mathit{\Omega } {\left( {\sigma _x^k{\rm{ \mathit{ δ} }}\varepsilon _x^k + \bar \tau _{xz}^k{\rm{ \mathit{ δ} }}\gamma _{xz}^k} \right){\rm{d}}\mathit{\Omega }} + \int_\mathit{\Omega } {{\rm{ \mathit{ δ} }}\bar \tau _{xz}^k\left( {\gamma _{xz}^k - \bar \gamma _{xz}^k} \right){\rm{d}}\mathit{\Omega }} - \\ \;\;\;\;\;\;\;{\rm{ \mathit{ δ} }}{W_{\rm{e}}} = 0 \end{array} $	（26）

式中：Ω为面积；γ_xz^k=τ_xz^k/G^k；W_e为外力做的功，

$ \begin{array}{l} {\rm{ \mathit{ δ} }}{W_{\rm{e}}} = \int_0^L {\left[ {q{\rm{ \mathit{ δ} }}w - \left( {{T_{x0}}{\rm{ \mathit{ δ} }}u_x^k(0,z) + {T_{z0}}{\rm{ \mathit{ δ} }}w(0)} \right) + } \right.} \\ \;\;\;\;\;\;\;\left. {\left( {{T_{xL}}{\rm{ \mathit{ δ} }}u_x^k(L,z) + {T_{zL}}{\rm{ \mathit{ δ} }}w(L)} \right)} \right]{\rm{d}}x \end{array} $	（27）

其中：L为梁的长度；q为横向载荷；T_x0和T_z0分别为梁左端的轴向载荷和剪切载荷；T_xL和T_zL分别为梁右端的轴向载荷和剪切载荷。

式(26)可以写成

$ \int_\mathit{\Omega } {\left( {\sigma _x^k{\rm{ \mathit{ δ} }}\varepsilon _x^k + \bar \tau _{xz}^k{\rm{ \mathit{ δ} }}\gamma _{xz}^k} \right){\rm{d}}\mathit{\Omega }} - {\rm{ \mathit{ δ} }}{W_{\rm{e}}} = 0 $	（28）

$ \int_\mathit{\Omega } {{\rm{ \mathit{ δ} }}\bar \tau _{xz}^k\left( {\gamma _{xz}^k - \bar \gamma _{xz}^k} \right){\rm{d}}\mathit{\Omega }} = 0 $	（29）

将式(12)和式(19)代入式(29)中可得

$ {{\mathit{\boldsymbol{\bar U}}}_x} = {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{E}}}{\mathit{\boldsymbol{\delta }}_{xz}} $	（30）

式中：$\boldsymbol{C}_{E}=\left(\int_{z_{1}}^{z_{n+1}}\left(\frac{\overline{\boldsymbol{C}}_{z}^{\mathrm{T}} \overline{\boldsymbol{C}}_{z}}{G^{k}}\right) \mathrm{d} z\right)^{-1} \int_{z_{1}}^{z_{n+1}} \overline{\boldsymbol{C}}_{z}^{\mathrm{T}} \boldsymbol{B}_{x z}^{k} \mathrm{d} z$。

将式(30)代入式(19)，可得最终的横向剪应力的表达式

$ \bar \tau _{xz}^k\left( {x,z} \right) = {{\mathit{\boldsymbol{\bar C}}}_z}{\mathit{\boldsymbol{C}}_E}{\mathit{\boldsymbol{\delta }}_{xz}} = \mathit{\boldsymbol{G}}_{xz}^k{\mathit{\boldsymbol{\delta }}_{xz}} $	（31）

将式(12)、式(13)和式(31)代入式(28)中，分部积分可得平衡方程

$ \left\{ {\begin{array}{*{20}{l}} {\delta {u_0}:}&{{N_{x,x}} = 0}\\ {\delta {u_1}:}&{{M_{\mathit{\varPhi }1,x}} - {V_{\mathit{\varPhi }1}} = 0}\\ {\delta {u_3}:}&{{M_{\mathit{\varPhi }2,x}} - {V_{\mathit{\varPhi }2}} = 0}\\ {\delta w:}&{{V_{\mathit{\varPhi }1,x}} + q = 0} \end{array}} \right. $	（32）

和在x=0和x=L处的边界条件：

$ \left\{ {\begin{array}{*{20}{l}} {\delta {u_0}:}&{{u_0} = {{\bar u}_0}\quad 或\quad {N_x}(\bar \alpha ) = {{\bar N}_{x\bar \alpha }}}\\ {\delta {u_1}:}&{{u_1} = {{\bar u}_1}\quad 或\quad {M_{\mathit{\varPhi }1}}(\bar \alpha ) = {{\bar M}_{\mathit{\varPhi }1\bar \alpha }}}\\ {\delta {u_3}:}&{{u_3} = {{\bar u}_3}\quad 或\quad {M_{\mathit{\varPhi }2}}(\bar \alpha ) = {{\bar M}_{\mathit{\varPhi }2\bar \alpha }}}\\ {\delta w:}&{{w_0} = {{\bar w}_0}\quad 或\quad {V_{\mathit{\varPhi }1}}(\bar \alpha ) = {{\bar V}_{\mathit{\varPhi }1\bar \alpha }}} \end{array}} \right. $

（33）

式中：α=0, L。

$ \left( {{N_x},{M_{\mathit{\varPhi }1}},{M_{\mathit{\varPhi }2}}} \right) = \int_{{z_1}}^{{z_{n + 1}}} {\left( {\sigma _x^k,\mathit{\varPhi }_1^k\sigma _x^k,\mathit{\varPhi }_2^k\sigma _x^k} \right){\rm{d}}z} $	（34）

$ \left( {{V_{\mathit{\varPhi }1}},{V_{\mathit{\varPhi }2}}} \right) = \int_{{z_1}}^{{z_{n + 1}}} {\left( {\mathit{\varPhi }_{1,z}^k\bar \tau _x^k,\mathit{\varPhi }_{2,z}^k\bar \tau _x^k} \right){\rm{d}}z} $	（35）

(N_xα, M_Φ1α, M_Φ2α, V_Φ1α)=∫_z1^z_n+1(T_xα, Φ₁^kT_xα, Φ₂^kT_xα, T_zα)dz，式(34)和式(35)可以写成

$ \left[ {\begin{array}{*{20}{c}} {{N_x}}\\ {{M_{\mathit{\varPhi }1}}}\\ {{M_{\mathit{\varPhi }2}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{A_{11}}}&{{B_{12}}}&{{B_{13}}}\\ {{B_{12}}}&{{D_{11}}}&{{D_{12}}}\\ {{B_{13}}}&{{D_{12}}}&{{D_{22}}} \end{array}} \right]{\mathit{\boldsymbol{\delta }}_x} $	（36）

$ \left[ {\begin{array}{*{20}{c}} {{V_{\mathit{\varPhi }1}}}\\ {{V_{\mathit{\varPhi }2}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{H_{11}}}&{{H_{12}}}&{{H_{13}}}\\ {{H_{21}}}&{{H_{22}}}&{{H_{23}}} \end{array}} \right]{\mathit{\boldsymbol{\delta }}_{xz}} $	（37）

其中：(A₁₁, B₁₂, D₁₁)=$\int_{z_{1}}^{z_{n+1}} G^{k}\left(1, \varPhi_{1}^{k}, \varPhi_{2}^{k}\right) \mathrm{d} z$；(B₁₃, D₁₂, D₂₂)=$\int_{z_{1}}^{z_{n+1}} G^{k} \varPhi_{2}^{k}\left(1, \quad \varPhi_{1}^{k}, \quad \varPhi_{2}^{k}\right) \mathrm{d} z$；H_ij=$\int_{z_{1}}^{z_{n+1}} \varPhi_{i, z}^{k} G_{x z}^{k}(j) \mathrm{d} z$, i=1~2, j=1~3。

3 算例

为了测试所提出的新锯齿理论(NZT)的计算精度和稳健性，以层数较多的纤维增强复合材料层合梁和面板，芯层材料属性差异较大的三明治夹层梁以及出现分层的非对称夹层梁为研究对象，考虑静力弯曲特性，给出了上表面承受正弦载荷的简支梁和自由端承受集中载荷的悬臂梁的解析解。计算结果与现有多种理论的计算结果进行了对比。具体包括Pagano^[44]提出的精确三维弹性解(Exact)，Tahani^[45]给出的分层理论解(BLWT)，Tessler等^[46]给出的高精度有限元的解(2D-FEM)和混合修正锯齿理论解(RZTM)^[24]，Wu和Chen^[30]提出的混合整体高阶锯齿理论解(GHZTM)，Ren等^[31]提出的C⁰型锯齿理论解(ZZTC-C⁰)，Han等^[47]提出的C⁰型高阶锯齿理论解(EHOZT)，Oñate等^[48]给出的二维有限元解(PS)和基于修正锯齿理论提出的两节点梁单元的解，基于经典锯齿理论(ZZT)^[49]和混合变分原理提出的混合锯齿理论解(ZZTM)、Reddy梁理论解(Reddy)、一阶剪切变形理论解(FSDT)以及基于FSDT和混合变分原理提出的混合一阶剪切变形理论解(FSDTM)。

具体的材料参数如下：

1) 层合梁材料力学特性^[50]

E₁=25 GPa, E₂=E₃=1 GPa, G₂₃=0.5 GPa, G₁₂=G₁₃=0.2 GPa, ν₁₂=ν₁₃=ν₂₃=0.25。

2) 夹层梁材料力学特性^[46]

下表层：E=73 GPa, G=29.2 GPa；夹芯层：E=0.073 GPa, G=0.029 2 GPa；上表层：E=21.9 GPa, G=8.76 GPa。

3.1 层数较多的简支梁

该算例是承受正弦载荷的简支梁，如图 4所示。层合梁每层的厚度和材料性质相同(材料1)。

简支梁的边界条件如下：

$ x = 0,L,w = {N_x} = {M_{\mathit{\varPhi }1}} = {M_{\mathit{\varPhi }2}} = 0 $	（38）

满足全部边界条件的位移函数可设为

$ \begin{array}{l} {w_0}(x) = {w_{00}}\sin ({\rm{ \mathsf{ π} }}x/L),\\ \left( {{u_0}(x),{u_1}(x),{u_3}(x)} \right) = \\ \;\;\;\;\;\;\left( {{u_{00}},{u_{10}},{u_{30}}} \right)\cos ({\rm{ \mathsf{ π} }}x/L) \end{array} $	（39）

位移和应力的无量纲化为

$ {{\bar u}^k} = \frac{{{E_1}}}{{{q_0}h}}{u^k},\bar w = \frac{{100{E_1}{h^3}}}{{{q_0}{L^4}}}w,\bar \sigma _x^k = \frac{{\sigma _x^k}}{{{q_0}}},\bar \tau _{xz}^k = \frac{{\bar \tau _{xz}^k}}{{{q_0}}} $

式中：h为梁的厚度。

为了验证NZT预测多层梁的变形和应力的能力，分别计算了3层对称铺设层合梁(铺设角：0°/90°/0°)，25层对称铺设层合梁(铺设角：0°/(90°/0°)₁₂)，4层反对称铺设层合梁(铺设角：(0°/90°)₂)，8层反对称铺设层合梁(铺设角：(90°/0°)₄)，50层反对称铺设层合梁(铺设角：(0°/90°)₂₅)和100层反对称铺设层合梁(铺设角：(0°/90°)₅₀)。图 5给出了3层层合梁的轴向位移、轴向应力和横向剪应力沿厚度方向的变化图，并将提出的NZT的计算结果与三维弹性解(Exact)、分层理论解(BLWT)、RZTM、GHZTM、ZZTC-C⁰、Reddy和FSDT的结果进行对比。由图 5(a)可知，提出的NZT能够准确地描述面内位移沿厚度方向的C_z⁰变化。图 6比较了不同梁模型得到的25层对称铺设层合梁的轴向应力和横向剪应力。由图 5(b)、图 5(c)和图 6可以看出，NZT能够准确预测层数不大于25层的对称铺设层合梁的应力，尤其是横向剪应力，和三维弹性解(Exact)、分层理论解(BLWT)吻合较好，预测精度高于RZTM、GHZTM和ZZTC-C⁰。而Reddy梁理论和一阶剪切变形理论(FSDT)预测厚梁的横向剪应力误差较大。图 7和图 8分别给出了4层和8层反对称铺设层合梁的横向剪应力。从图中可知，NZT的预测精度高于RZTM和GHZTM。然而，NZT依旧具有一定的误差，最大误差约为6%。而从图 9和图 10可以看出，随着反对称层合梁的层数增大到50层、100层，NZT的结果反而与三维弹性解(Exact)吻合较好，计算精度逐步提高。

为了定量地对比本文提出的NZT的计算精度。4种层合梁的最大挠度和梁左端最大横向剪应力分别在表 1和表 2中给出，不同梁理论相对三维弹性解(Exact)的相对误差在表 1和表 2中的括号中给出。表 1和表 2中位移和应力的相对误差为：$\frac{|\mathrm{NZT}-\mathrm{Exact}|}{\mathrm{Exact}} \times 100 \%$。

从表 1和表 2可以看出, 本文提出的NZT的计算结果和三维弹性解(Exact)吻合较好，对4种层合梁，挠度和横向剪应力误差都不超过1%。而一阶剪切变形理论(FSDT)预测挠度误差较大，最大达到27.7%，Reddy梁理论预测横向剪应力的误差较大，最大误差超过50%。RZTM的计算精度会随着层数的增加而降低，当层数由3层增加到25层，RZTM的横向剪应力的误差从1.7%增大到9.5%。当层数超过50层时，GHZTM预测横向剪应力的精度降低。当层数较少时，ZZTC-C⁰计算横向剪应力误差较大，超过20%。然而，随着层数的增加，ZZTC-C⁰预测横向剪应力的精度逐渐提高。总结表 1和表 2分析可得如下结论：相比于FSDT、Reddy和RZTM，提出的NZT的变形和应力的计算精度都有大幅度提升。当层数少时，提出的NZT预测横向剪应力的精度高于ZZTC-C⁰，当层数较多时，提出的NZT预测横向剪应力的精度高于GHZTM，NZT和EHOZT具有相当的计算精度，表现出较好的稳健性。

综上，本文提出的NZT能够准确地预测对称铺设多层厚梁和层数较多的反对称铺设多层厚梁的挠度和应力。对于层数较少的反对称铺设层合厚梁横向剪应力，NZT具有一定的误差，但误差随着层数的增加，逐渐减小。

3.2 面板和芯层材料属性差异较大的三明治悬臂夹层梁

准确地预测面板和芯层材料属性差异较大的夹层梁的变形和应力十分困难，Tessler提出的RZTM很好地解决了该难题。本文提出的新锯齿理论模型不仅能够准确地预测层数较多的复合材料层合梁的变形和应力，对该类结构也具有较好的预测精度。这里以3层软核夹层悬臂梁为例，每层厚度为(0.1h/0.8h/0.1h)，夹层梁每层具有不同的材料性质(材料2)。悬臂梁在自由端承受横向载荷F，如图 11所示。

边界条件为

$ \begin{array}{l} x = 0,{u_0} = {u_1} = {u_3} = {w_0} = 0\\ x = L,{N_x} = {M_{\mathit{\varPhi }1}} = {M_{\mathit{\varPhi }2}} = 0,{V_{\mathit{\varPhi }1}} = F \end{array} $	（40）

满足全部边界条件的位移函数^[32]可设为

$ \begin{array}{l} {u_0}(x) = \left( { - {C_8} + {C_3}{C_7} - \frac{{{C_2}{C_7}}}{{{R^2}D_{11}^*}}} \right)\left( {{a_1}\cosh (Rx) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{a_2}\sinh (Rx)} \right) - \frac{{{C_2}{C_7}{a_3}}}{{2D_{11}^*}}{x^2} + {a_6}x + {a_7} \end{array} $

$ \begin{array}{l} {u_1}(x) = \left( { - {C_3} + \frac{{{C_2}}}{{{R^2}D_{11}^*}}} \right)\left( {{a_1}\cosh (Rx) + } \right.\\ \;\;\;\;\;\;\;\;\left. {{a_2}\sinh (Rx)} \right) + \frac{{{C_2}{a_3}}}{{2D_{11}^*}}{x^2} + {a_4}x + {a_5} \end{array} $

$ {u_3}(x) = {a_1}\cosh (Rx) + {a_2}\sinh (Rx) + {a_3} $

$ \begin{array}{l} {w_0}(x) = \\ \left[ {\frac{{{C_3}}}{R} - \frac{{{C_2}}}{{{R^3}D_{11}^*}} + \frac{1}{R}\left( {\frac{{{C_2}{C_5}}}{{D_{11}^*}} - {C_4}} \right) + R\left( {{C_6} - {C_3}{C_5}} \right)} \right] \cdot \\ \;\;\;\;\;\;\;\left( {{a_1}\sinh (Rx) + {a_2}\cosh (Rx)} \right) - \frac{{{C_2}{a_3}}}{{6D_{11}^*}}{x^3}\\ \;\;\;\;\;\;\; - \frac{{{a_4}}}{2}{x^2} + \left[ {\left( {\frac{{{C_2}{C_5}}}{{D_{11}^*}} - {C_4}} \right){a_3} - {a_5}} \right]x + {a_8} \end{array} $

（41）

式中：

$ D_{{\rm{11}}}^ * = \frac{{{D_{11}}{A_{11}} - {{\left( {{B_{12}}} \right)}^2}}}{{{A_{11}}}},D_{{\rm{12}}}^ * = \frac{{{D_{12}}{A_{11}} - {B_{12}}{B_{13}}}}{{{A_{11}}}} $

$ D_{22}^* = \frac{{{D_{22}}{A_{11}} - {{\left( {{B_{13}}} \right)}^2}}}{{{A_{11}}}},{C_4} = \frac{{{H_{12}}}}{{{H_{11}}}},{C_5} = \frac{{D_{11}^ * }}{{{H_{11}}}} $

$ {C_6} = \frac{{D_{12}^*}}{{{H_{11}}}},{C_7} = \frac{{{B_{12}}}}{{{A_{11}}}},{C_8} = \frac{{{B_{13}}}}{{{A_{11}}}} $

$ {C_1} = \frac{{\left[ {{{\left( {D_{12}^*} \right)}^2} - D_{11}^*D_{22}^*} \right]{H_{11}}}}{{D_{12}^*{H_{11}} - D_{11}^*{H_{12}}}} $

$ {C_2} = \frac{{\left[ {{H_{11}}{H_{22}} - {{\left( {{H_{12}}} \right)}^2}} \right]D_{11}^*}}{{D_{12}^*{H_{11}} - D_{11}^*{H_{12}}}} $

$ {C_3} = \frac{{D_{22}^*{H_{11}} - D_{12}^*{H_{12}}}}{{D_{12}^*{H_{11}} - D_{11}^*{H_{12}}}} $

$ R = \sqrt {\left| {{C_2}/{C_1}} \right|} $

其中:C₂/C₁ < 0，如果C₂/C₁>0，则式(41)中的双曲正弦函数和双曲余弦函数分别改为正弦函数和余弦函数即可。

位移和应力的无量纲化为

$ \tilde w = \frac{{{D_{11}}}}{{10F{L^3}}}w,\tilde \sigma _x^k = \frac{{{h^2}}}{{FL}}\sigma _x^k,\tilde \tau _{xz}^k = \frac{{{h^2}}}{{FL}}\bar \tau _{xz}^k $

图 12和图 13分别给出了应力沿厚度方向分布图和中面挠度沿轴向分布图。从图 12和图 13可以看出，一阶剪切变形理论(FSDT)严重低估了最大挠度和最大轴向应力，挠度最大误差可达到82%，轴向应力最大误差可达到76.5%。受FSDT直法线假定的限制，即使使用混合变分原理进行修正，混合一阶剪切变形理论(FSDTM)同样具有较大的误差。图 12(b)中的2D-FEM解是使用Abaqus软件计算得到，梁模型使用20万(1 000(长)×200(高))个4节点，平面应力，完整积分，CPS4单元进行离散。结果显示本文提出的NZT与2D-FEM解吻合较好，具有较高的精度。

3.3 考虑分层的非对称悬臂层合梁

Oñate等^[48]曾指出“复合材料层合梁的分层预测是所有梁模型面临的挑战”。本文借用Oñate给出的分层的层合梁算例来验证提出的新锯齿理论预测层合梁分层的能力。该算例为一端固支一端承受集中载荷的非对称悬臂厚梁(L/h=5)。为了模拟分层，在层合梁的第2层和第4层之间引入了一层非常薄的界面层(见图 14)，并赋予界面层一个非常小的剪切模量(见表 3)。在表 4对比了自由端的挠度，与二维有限元(PS)结果相比，提出的新锯齿理论的相对误差仅为0.7%，而划分100个单元的高精度梁单元(LRZ)的计算误差高达11.4%。图 15和图 16分别给出了沿厚度分布的轴向位移和横向剪应力的对比图。当分层发生时，提出的新锯齿理论模型能够准确预测出界面处的轴向位移发生跳跃，横向剪应力消失的现象，具有更高的计算精度(相比于LRZ后处理结果：LRZ-100-Nz)。

4 结论

本文通过构造一个新的线性分段锯齿函数，提出了一种能够准确预测层合梁结构横向剪应力的新锯齿理论模型(NZT)。为了验证该理论模型精度，计算了层数较多的简支梁、材料属性差异较大的三明治悬臂梁以及发生分层的非对称悬臂梁的位移和应力。可得到如下结论：

1) 新锯齿函数具备描述每层材料变化和不同横截面处变化规律不同的能力，使得NZT能够准确预测横向各向异性的多层梁的变形和应力。

2) NZT位移场中仅含有4个与层数无关的未知量，不含有横向位移的一阶导数，构造梁单元时，仅需使用C⁰插值形函数。

3) NZT能够预先满足横向剪应力层间连续，无需后处理就能够准确预测层合梁的层间应力。

4) 算例结果显示：NZT能够准确预测层数较多的层合梁和材料属性差异较大的三明治夹层梁的变形和应力，能够准确预测层合梁的分层。