夹层组合梁具有可以充分利用组成材料的优点，在保证强度与刚度的基础上，节省材料、减轻重量，在正常使用状态和承载能力极限状态都能显示出其优势. 组合梁的连接件往往不足以保证各层之间实现完美连接（full-interaction），导致“部分作用（partial-interaction）”^[1]的现象，即组合梁受力时各层之间将会产生层间的黏结滑移作用. 层间滑移的存在大大增加了组合梁计算分析的难度，国内外学者对此展开了广泛的研究.

以钢-混凝土组合梁为代表的2层组合梁的研究为工程应用提供了宝贵的解决方法. 在实际工程中，3层组合梁也被广泛应用于工程之中，如钢-混-钢夹心组合梁^[2]、双钢板-混凝土组合防护结构^[3]、夹心约束阻尼梁^[4]、波纹钢腹板混凝土梁^[5]、夹层玻璃柱^[6]等，因此，对3层组合梁的研究具有重要意义.

Newmark^[7]通过实验提出层间滑移与剪力之间的线性假设，基于Euler-Bernoulli梁理论建立钢-混凝土组合梁的控制方程，为部分作用组合梁的研究提供了值得借鉴的研究思路. 目前针对3层组合梁的研究较少，Chui等^[8]给出了在均布荷载和集中力荷载作用下3层简支组合梁的解析解. Sousa等^[9]研究基于Euler-Bernoulli梁和Timoshenko梁理论的多层组合梁的解析解. 在数值分析研究方面，Sousa^[10]基于之前的解析方法提出新的有限元公式. Keo等^[11]利用控制方程导出精确的刚度矩阵并提出精确有限元公式. Lin等^[12]基于Timoshenko梁理论提出考虑独立截面转角和横向剪切变形的3层部分作用组合梁有限元分析方法.

数值分析方法为组合梁变形分析带来了便利，但目前关于3层组合梁有限元数值分析中可以显著提高计算效率的显式刚度矩阵的研究相对缺乏，这为其工程应用带来不便. 另外，基于位移的组合梁有限元模拟中“滑移/曲率自锁”^{[9, 12-14]}是常见的问题，它会引起层间滑移场的振荡从而导致收敛性问题.

本研究在Xu等^[15]工作的基础上，采用Timoshenko梁理论，将高阶插值函数引入内部自由度以避免滑移/曲率锁定问题，并利用MATLAB编译有限元程序. 为了提高计算效率，推导并给出了3层组合梁的8×8刚度矩阵、质量矩阵以及几何刚度矩阵元素的显式表达式. 进行算例验证和应用，通过与文献和ABAQUS模拟结果对比分析，验证本研究提出的新型夹层组合梁有限元方法的准确性、收敛性和稳定性.

夹层组合梁剖面图如图1所示. 图中，E_i、G_i、I_i、A_i和κ_i$ (i = 1,2,3) $分别表示每层梁的弹性模量、剪切模量、惯性矩、截面积和剪切修正系数；L₀和H表示组合梁的长度和高度；$ h'_1 $、$ h'_2 $、$ h'_3 $为每层形心轴与相邻界面的距离；$ {h_1} $、$ {h_2} $为相邻层结构形心轴的距离，$ h $为顶层与底层形心轴距离，$ h = {h_1}+{h_2} $. 在x-z平面内，夹层组合梁整个横截面的形心轴与x轴重合.

推导建立在以下几个基本假设的基础上^[15-16]：1）各层梁均符合线性变形假设；2）夹层结构的层间剪力与对应的层间滑移成正比；3）忽略界面法向分离现象；4）根据Timoshenko梁理论，允许存在横向剪切变形.

如图2所示为层间滑移、截面转角和轴向位移的几何关系，夹层梁的层间滑移与横截面转角$ \psi $的关系^[17]如下：

(1)$ {u_{{\mathrm{s}}1}} = {u_2} - {u_1}+\psi {h_1} \text{，} {u_{{\mathrm{s}}2}} = {u_3} - {u_2}+\psi {h_2}. $

式中：$ {u_1} $、$ {u_2} $和$ {u_3} $为各层的轴向位移，u_s1、u_s2分别为中间层与顶层和底层之间的相对滑移.

如图3所示为长度为dx的微小单元. 图中，M、Q、F分别表示整个横截面的弯矩、剪力和轴力，N_i、M_i和Q_i(i=1,2,3) 为各层梁的轴力、弯矩和剪力，V_i为层间剪力，q、m分别表示作用在该单元上的均布载荷和弯矩.

根据单元在横截面上的平衡条件，可以得到

(2)$ \frac{{{{\mathrm{d}}}M}}{{{{\mathrm{d}}}x}} = Q+m \text{，} $

(3)$ \frac{{{{\mathrm{d}}}Q}}{{{{\mathrm{d}}}x}} = - q . $

同时，根据x方向的平衡条件可以得到

(4)$ {N_1}+{N_2}+{N_3} = 0\text{，} $

(5)$ {V_1} = - \frac{{{\rm{d}}{N_1}}}{{{\rm{d}}x}}，\qquad{V_2} = {V_1} - \frac{{{\rm{d}}{N_2}}}{{{\rm{d}}x}}.$

整个单元横截面的弯矩可以表示为

(6)$ M = {M_1}+{M_2}+{M_3} - {N_1}{h_1}+{N_3}{h_2}. $

各层间纵向剪力与连接件的抗剪刚度呈线性关系：

(7)$ {V_1} = {k_{{\mathrm{s}}1}}{u_{{\mathrm{s}}1}}\text{，}{V_2} = {k_{{\mathrm{s}}2}}{u_{{\mathrm{s}}2}}. $

根据Timoshenko梁理论，内力与变形关系为

(8)$ {M_1} = - {E_1}{I_1}\frac{{{{\mathrm{d}}} \psi }}{{{{\mathrm{d}}} x}}\text{，}{M_2} = - {E_2}{I_2}\frac{{{{\mathrm{d}}} \psi }}{{{{\mathrm{d}}} x}}\text{，}{M_3} = - {E_3}{I_3}\frac{{{{\mathrm{d}}} \psi }}{{{{\mathrm{d}}} x}} ,$

(9)$ Q = C(\frac{{{{\mathrm{d}}}w}}{{{{\mathrm{d}}}x}} - \psi ). $

式中：$C = {\kappa _1}{G_1}{A_1}+{\kappa _2}{G_2}{A_2}+{\kappa _3}{G_3}{A_3}$，w为梁的挠度.

各层轴向力${N_i}$(i =1，2，3)与相应轴向位移的关系为

(10)$ {N_i} = {E_i}{A_i}\frac{{{{\mathrm{d}}} {u_i}}}{{{{\mathrm{d}}} x}}. $

为了推导夹层梁的有限元单元，首先须建立夹层梁的能量表达式：

(11)$ \pi = {\pi _1}+{\pi _2}+{\pi _3} - {W_q}. $

式中：$ {\pi _1} $为弯曲应变能，$ {\pi _2} $为剪切应变能，$ {\pi _3} $为剪切连接应变能，$ {W_q} $为外力做功. 表达式分别如下：

(12)$ \begin{split} {\pi _1}& = \int_0^{L_0} {\frac{1}{2}E{I_0} {{\left(\frac{{{{\mathrm{d}}} \psi }}{{{{\mathrm{d}}} x}}\right)}^2}{{\mathrm{d}}} x} + \int_0^{L_0} {\frac{1}{2}{E_1}{A_1}{{\left(\frac{{{{\mathrm{d}}} {u_1}}}{{{{\mathrm{d}}} x}}\right)}^2}{{\mathrm{d}}} x} \\ &+ \int_0^{L_0} {\frac{1}{2}{E_2}{A_2} {{\left(\frac{{{{\mathrm{d}}} {u_2}}}{{{{\mathrm{d}}} x}}\right)}^2}{{\mathrm{d}}} x} + \int_0^{L_0} {\frac{1}{2}{E_3}{A_3}{{\left(\frac{{{{\mathrm{d}}} {u_3}}}{{{{\mathrm{d}}} x}}\right)}^2}{{\mathrm{d}}} x} ,\end{split} $

(13)$ {\pi _2} = \int_0^{L_0} {\frac{1}{2}} C{(\frac{{{{\mathrm{d}}}w}}{{{{\mathrm{d}}}x}} - \psi )^2}{{\mathrm{d}}}x, $

(14)$ {\pi _3} = \int_0^{L_0} {\frac{1}{2}{k_{{\mathrm{s}}1}}u_{{\mathrm{s}}1}^2{\mathrm{d}}x} +\int_0^{L_0} {\frac{1}{2}{k_{{\mathrm{s}}2}}u_{{\rm{s}}2}^2{\mathrm{d}}x}, $

(15)$ {W_q} = \int_0^{L_0} {(qw+m\psi ){{\mathrm{d}}}x+(\overline Q w - \overline M \psi - {{\overline N }_1}{u_{\mathrm{s}}})_0^{L_0}}. $

式中：$E{I_0} = {E_1}{I_1}+{E_2}{I_2}+{E_3}{I_3}$，$ \overline M $、$ \overline Q $、$ {\overline N _1} $为施加在边界上的外力.

将组合梁沿轴向离散为若干个无穷小单元，如图4所示为其中第e个单元. 记单元节点位移为$ {{{\boldsymbol{u}}}_e} = {[{w_i},{\psi _i},{u_{{\mathrm{s}}i1}},{u_{{\mathrm{s}}i2}},{w_j},{\psi _j},{u_{{\mathrm{s}}j1}},{u_{{\mathrm{s}}j2}}]^{\text{T}}} $. 单元起点和终点在x方向的坐标分别为$ {x_i} $和$ {x_j} $，因此在局部单元内任意点的坐标可以表示为

(16)$ \alpha = \frac{1}{l}({x_j} - x) \text{，} \beta = \frac{1}{l}(x - {x_j}). $

式中：l为该单元的长度.

(17)$ l = {x_j} - {x_i}. $

所假设位移的插值函数表达式如下：

(18)$ \left.\begin{aligned} &w(x) = {w_i}\alpha +{w_j}\beta +{\xi _1}\alpha \beta +{\xi _2}({\alpha ^2}\beta - \alpha {\beta ^2}), \\& \psi (x) = {\psi _i}\alpha +{\psi _j}\beta +{\eta _1}\alpha \beta, \\ &{u_{{\rm{s}}1}}(x) = {u_{{\rm{s}}i1}}\alpha +{u_{{\rm{s}}j1}}\beta +{\eta _2}\alpha \beta, \\& {u_{{\rm{s}}2}}(x) = {u_{{\rm{s}}i2}}\alpha +{u_{{\rm{s}}j2}}\beta +{\eta _3}\alpha \beta. \end{aligned}\right\} $

式中：$ {\xi }_{1}、{\xi }_{2}、\text{ }{\eta }_{1}、{\eta }_{2}、\text{ }{\eta }_{3} $为新引入的未知参数，这些参数只影响单元的局部变形，而非外部节点或其他单元的变形，可称为“内部自由度”^[18]，记$ {{{\boldsymbol{u}}}^i} = {[{\xi _1},{\eta _1},{\xi _2},{\eta _2},{\eta _3}]^{\text{T}}} $. 内部自由度与节点自由度的关系可以通过单元的最小势能原理来确定.

将式（18）代入式（11），可以得到新的能量表达式：

(19)$ {\pi _e}' = \frac{1}{2}{({{{\boldsymbol{u}}}_e}')^{\text{T}}}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}_0^e}&{{{\boldsymbol{L}}^{\text{T}}}} \\ {\boldsymbol{L}}&{\boldsymbol{H}} \end{array}} \right]{{{\boldsymbol{u}}}_e}'. $

式中：$ {{{\boldsymbol{u}}}_e}' = {\left[ {{{\boldsymbol{u}}}_e^{\text{T}},{\xi _1},{\eta _1},{\xi _2},{\eta _2},{\eta _3}} \right]^{\text{T}}} $；$ {\boldsymbol{K}}_0^e $、$ {\boldsymbol{H}} $和$ {\boldsymbol{L}} $为含内部自由度单元刚度矩阵的子块，具体表达式见附录A.

为了消除内部的自由度，利用内部自由度应使单元的势能最小化的关系：

(20)$ \frac{{\partial {\pi _e}'}}{{\partial {\xi _1}}} = \frac{{\partial {\pi _e}'}}{{\partial {\xi _2}}} = \frac{{\partial {\pi _e}'}}{{\partial {\eta _1}}} = \frac{{\partial {\pi _e}'}}{{\partial {\eta _2}}} = \frac{{\partial {\pi _e}'}}{{\partial {\eta _3}}} = 0. $

将式（19）重新表示为

(21)$ \begin{split} {\boldsymbol{L}}&{[{w_i},{\psi _i},{u_{{\mathrm{s}}1i}},{u_{{\mathrm{s}}2i}},{w_j},{\psi _j},{u_{{\mathrm{s}}1j}},{u_{{\mathrm{s}}2j}}]^{\text{T}}}+\\&\qquad{\boldsymbol{H}}{[{\xi _1},{\xi _2},{\eta _1},{\eta _2},{\eta _3}]^{\text{T}}} = {\bf{0}}.\end{split} $

进而可以得到内部自由度与单元节点上位移的关系：

(22)$ {{{\boldsymbol{u}}}^i} = \dfrac{1}{{1 - \mu }}\left[ {\begin{array}{*{20}{c}} 0&{\dfrac{{l(1 - \mu)}}{2}}&0&0&0&{\dfrac{{ - l(1 - \mu)}}{2}}&0&0 \\ {\dfrac{{Cl}}{{4{\lambda _1}}}}&{\dfrac{{C{l^2}}}{{8{\lambda _1}}}}&{\dfrac{{{k_{\rm{s}}}_1{l^2}}}{{8{\lambda _3}}}}&{\dfrac{{{k_{\rm{s}}}_2{l^2}}}{{8{\lambda _2}}}}&{ - \dfrac{{Cl}}{{4{\lambda _1}}}}&{\dfrac{{C{l^2}}}{{8{\lambda _1}}}}&{\dfrac{{{k_{\rm{s}}}_1{l^2}}}{{8{\lambda _3}}}}&{\dfrac{{{k_{\rm{s}}}_2{l^2}}}{{8{\lambda _2}}}} \\ {\dfrac{{ - C{l^2}}}{{24{\lambda _1}}}}&{\dfrac{{ - C{l^3}}}{{48{\lambda _1}}}}&{\dfrac{{ - {k_{\rm{s}}}_1{l^3}}}{{48{\lambda _3}}}}&{\dfrac{{ - {k_{\rm{s}}}_2{l^3}}}{{48{\lambda _2}}}}&{\dfrac{{C{l^2}}}{{24{\lambda _1}}}}&{\dfrac{{ - C{l^3}}}{{48{\lambda _1}}}}&{\dfrac{{ - {k_{\rm{s}}}_1{l^3}}}{{48{\lambda _3}}}}&{\dfrac{{ - {k_{\rm{s}}}_2{l^3}}}{{48{\lambda _2}}}} \\ {\dfrac{{Cl}}{{4{\lambda _3}}}}&{\dfrac{{C{l^2}}}{{8{\lambda _3}}}}&{\dfrac{{{k_{\rm{s}}}_1{l^2}}}{{8{\lambda _4}}}}&{\dfrac{{ - {k_{\rm{s}}}_2{l^2}}}{{8{\lambda _6}}}}&{\dfrac{{ - Cl}}{{4{\lambda _3}}}}&{\dfrac{{C{l^2}}}{{8{\lambda _3}}}}&{\dfrac{{{k_{\rm{s}}}_1{l^2}}}{{8{\lambda _4}}}}&{\dfrac{{ - {k_{\rm{s}}}_2{l^2}}}{{8{\lambda _6}}}} \\ {\dfrac{{Cl}}{{4{\lambda _2}}}}&{\dfrac{{C{l^2}}}{{8{\lambda _2}}}}&{ - \dfrac{{{k_{\rm{s}}}_1{l^2}}}{{8{\lambda _6}}}}&{\dfrac{{{k_{\rm{s}}}_2{l^2}}}{{8{\lambda _5}}}}&{ - \dfrac{{Cl}}{{4{\lambda _2}}}}&{\dfrac{{C{l^2}}}{{8{\lambda _2}}}}&{\dfrac{{ - {k_{\rm{s}}}_1{l^2}}}{{8{\lambda _6}}}}&{\dfrac{{{k_{\rm{s}}}_2{l^2}}}{{8{\lambda _5}}}} \end{array}} \right]{{{\boldsymbol{u}}}_e}. $

式中：$\mu= {\left[5 C \gamma_5\left(1+12 \gamma_0\right)+6{\overline{E A h_2}}^2 k_{\mathrm{s} 1}\left(1+10 \gamma_2\right)\right.} \left.+6{\overline{E A h_1}}^2 k_{\mathrm{s} 2}\left(1+10 \gamma_3\right)\right]l_e^2 / 120 \gamma_4, $

$ \gamma_4=\overline{E A_3} \;\overline{E A h_1} \;\overline{E A h_2} \text {, } $

$ \begin{split} \gamma_5=&{\overline{E A_3}}^2-k_{\mathrm{s} 1} k_{\mathrm{s} 2} l_e^4\left(1+10 \gamma_2\right)\left(1+10 \gamma_3\right) / 100,\\\overline{E A h_1}=&\overline{E A_1} h_1+\overline{E A_3} h, \overline{E A h_2}=\overline{E A_2} h_2+\overline{E A_3} h,\end{split}$

$ \lambda_1=\frac{\gamma_4}{{\overline{E A_3}}^2-k_{\mathrm{s} 1} k_{\mathrm{s} 2} l_e^4\left(\gamma_2+1 / 10\right)\left(\gamma_3+1 / 10\right)}, $

$ \lambda_2=\frac{\gamma_4}{\overline{E A_3} \times \overline{E A h_1}-k_{\mathrm{s} 1} \overline{E A h_2} l_e^4\left(\gamma_2+1 / 10\right)}, $

$ \lambda_3=\frac{\gamma_4}{\overline{E A_3} \times \overline{E A h_2}-k_{\mathrm{s} 2} \overline{E A h_1} l_e^4\left(\gamma_3+1 / 10\right)}, $

$ \lambda_4=\frac{\gamma_4}{{\overline{E A h_2}}^2-C k_{\mathrm{s} 2} l_e^4\left(\gamma_3+1 / 10\right)\left(\gamma_0+1 / 12\right)}, $

$ \lambda_5=\frac{\gamma_4}{{\overline{E A h_1}}^2-C k_{\mathrm{s} 1} l_e^4\left(\gamma_2+1 / 10\right)\left(\gamma_0+1 / 12\right)}, $

$ \lambda_6=\frac{\gamma_4}{\overline{E A h_1} \times \overline{E A h_2}-C \overline{E A_3} l_e^2\left(\gamma_0+1 / 12\right)} . $

利用式(22)，位移表达式(18)可以重新表示为

$ \begin{aligned} w(x)=&\left[\alpha -m\alpha \beta (\alpha -\beta )\right]{w }_{i}+\left[\alpha \beta -m\alpha \beta (\alpha -\beta )\right]\frac{l}{2}{\psi }_{i}-\\ &lq\alpha \beta (\alpha -\beta ){u}_{{\rm{s}}1i}-lr\alpha \beta (\alpha -\beta ){u}_{{\rm{s}}2i}+\\ &\left[\beta +m\alpha \beta (\alpha -\beta )\right]{w }_{j}-\left[\alpha \beta +m\alpha \beta (\alpha -\beta )\right]\frac{l}{2}{\psi }_{j}-\\ &lq\alpha \beta (\alpha -\beta ){u}_{{\rm{s}}1j}-lr\alpha \beta (\alpha -\beta ){u}_{{\rm{s}}2j}，\end{aligned} $

$ \begin{aligned} \psi (x)&=\frac{6m}{l}\alpha \beta {w }_{i}+\left(\alpha +3m\alpha \beta \right){\psi }_{i}+6q\alpha \beta {u}_{{\rm{s}}1i}+6r\alpha \beta {u}_{{\rm{s}}2i}+\\ & \left( - \frac{6m}{l}\alpha \beta \right) {w }_{j} + \left(\beta + 3m\alpha \beta \right) {\psi }_{j} + 6q\alpha \beta {u}_{{\rm{s}}1j} + 6r\alpha \beta {u}_{{\rm{s}}2j} ，\end{aligned} $

$ \begin{split} {u}_{{\rm{s}}1}(x)&=\frac{6}{l}p\alpha \beta {w }_{i}+3p\alpha \beta {\psi }_{i}+\left(\alpha +s\alpha \beta \right){u}_{{\rm{s}}1i}-v\alpha \beta {u}_{{\rm{s}}2j}-\\ &v\alpha \beta {u}_{{\rm{s}}2i}-\frac{6}{l}p\alpha \beta {\omega }_{j}+3p\alpha \beta {\psi }_{j}+\left(\beta +s\alpha \beta \right){u}_{{\rm{s}}1j}，\end{split} $

(23)$ \begin{split} {u}_{{\rm{s}}2}(x)& = \frac{6}{l}n\alpha \beta {w}_{i} + 3n\alpha \beta {\psi }_{i} - \frac{{k}_{{\rm{s}}1}}{{k}_{{\rm{s}}2}}v\alpha \beta {u}_{{\rm{s}}1i} + \left(\alpha +t\alpha \beta \right){u}_{{\rm{s}}2i}-\\ & \frac{6}{l}n\alpha \beta {w }_{j} + 3n\alpha \beta {\psi }_{j} - \frac{{k}_{{\rm{s}}1}}{{k}_{{\rm{s}}2}}v\alpha \beta {u}_{{\rm{s}}1j} + \left( \beta + t\alpha \beta \right) {u}_{{\rm{s}}2j}.\end{split} $

式中：$ m = \dfrac{1}{{1 - \mu }}\dfrac{{C{l_e}^2}}{{24{\lambda _1}}} $，$ n = \dfrac{1}{{1 - \mu }}\dfrac{{C{l_e}^2}}{{24{\lambda _2}}} $，$ p = \dfrac{1}{{1 - \mu }}\dfrac{{C{l_e}^2}}{{24{\lambda _3}}}, $$ q = \dfrac{1}{{1 - \mu }}\dfrac{{{k_{{\rm{s}}1}}{l_e}^2}}{{48{\lambda _3}}}, $$ r = \dfrac{1}{{1 - \mu }}\dfrac{{{k_{{\rm{s}}2}}{l_e}^2}}{{48{\lambda _2}}}, $$ s = \dfrac{1}{{1 - \mu }}\dfrac{{{k_{{\rm{s}}1}}{l_e}^2}}{{8{\lambda _4}}} $，$ t = \dfrac{1}{{1 - \mu }}\dfrac{{{k_{{\rm{s}}2}}{l_e}^2}}{{8{\lambda _5}}}, $$ v = \dfrac{1}{{1 - \mu }}\dfrac{{{k_{{\rm{s}}2}}{l_e}^2}}{{8{\lambda _6}}} $.

将式(23)回代到能量方程（式(11)）并进行变分可以得到8×8的对称单元刚度矩阵K^e，具体表达式如附录B所示.

广义载荷对应于外力做功表达式（式(15)）中的线性项，即

(24)$ {({{{\boldsymbol{F}}}^e})^{\text{T}}}{{{\boldsymbol{u}}}_e} = \int_{{x_i}}^{{x_j}} {(qw+m\psi ){\mathrm{d}}x}. $

将新的插值函数（式(23)）代入式(24)得到广义载荷向量$ {{{\boldsymbol{F}}}^e} $.

组合梁的整体刚度矩阵和全局载荷向量表达式为

(25)$ {\boldsymbol{K}} = \sum\limits_{e = 1}^N {{{\boldsymbol{K}}^e}} \text{，} {{\boldsymbol{F}}} = \sum\limits_{e = 1}^N {{{{\boldsymbol{F}}}^e}}. $

式中：N为单位数量.

可以得到组合梁的总势能为

(26)$ \pi = \frac{1}{2}{{{\boldsymbol{u}}}^{\text{T}}}{\boldsymbol{K}}{{\boldsymbol{u}}} - {{{\boldsymbol{F}}}^{\text{T}}}{{\boldsymbol{u}}} $

采用最小势能原理可以得到

(27)$ {\boldsymbol{K}}{{\boldsymbol{u}}} = {{\boldsymbol{F}}} .$

式中：$ {{\boldsymbol{F}}} = {\left[ {{Q_1},{M_1},{N_1},\cdots,{Q_{N+1}},{M_{N+1}},{N_{N+1}}} \right]^{\text{T}}} ，\;{{\boldsymbol{u}}} = $$ {\left[ {{w_1},{\psi _1},{u_{{\rm{s}}1(1)}},{u_{{\rm{s}}2(1)}},\cdots,{w_{N+1}},{\psi _{N+1}},{u_{{\rm{s}}1(N+1)}},{u_{{\rm{s}}2(N+1)}}} \right]^{\text{T}}} $.

组合梁的动力行为可以通过哈密顿原理获得^[19]：

(28)$ \delta \int_{{t_1}}^{{t_2}} {(T - {\pi _1}+{\pi _2}+{\pi _3}+{W_q}){\mathrm{d}}t} = 0, $

(29)$ T = \frac{1}{2}\left\{ {\int_0^{L_0} {\rho {A_0}\dot w_{}^2{{\mathrm{d}}} x+\int_0^{L_0} {\rho {I_0}\dot \psi _{}^2{{\mathrm{d}}} x} } } \right\}. $

式中：T为动能，$ \dot w = \dfrac{{{\rm{d}}w}}{{{\rm{d}}t}} $，$ \dot \psi = \dfrac{{{\rm{d}}\psi }}{{dt}} $，t表示时间，$ \rho {A_0} = {\rho _1}{A_1}+{\rho _2}{A_2}+{\rho _3}{A_3} $，$ \rho {I_0} = {\rho _1}{I_1}+{\rho _2}{I_2}+{\rho _3}{I_3} $.

对于自由振动频率$ \omega $，位移分量可以表示为

(30)$ w(x,t) = \hat w(x){{\mathrm{e}}^{{\mathrm{i}}\omega t}} \text{，} \psi (x,t) = \hat \psi (x){{\mathrm{e}}^{{\mathrm{i}}\omega t}}\text{，} {u_{\rm{s}}}(x,t) = {\hat u_{\rm{s}}}(x){{\mathrm{e}}^{{\mathrm{i}}\omega t}}. $

式中：符号“^”表示相应的模态形状函数. 基于哈密顿原理，可以得到自振频率的变分表达式：

(31)$ \begin{split} \omega^2=&\mathrm{st} \sum_{e=1}^N \int_{x_i^e}^{x_j^e}\left[E I_0\left(\frac{\mathrm{d} \hat{\psi}}{\mathrm{d} x}\right)^2+\sum_{i=1}^3 E A_i\left(\frac{\mathrm{d} \hat{u}_i}{\mathrm{~d} x}\right)^2\right.+ \\&\left.C\left(\frac{\mathrm{d} \hat{w}}{\mathrm{~d} x}-\hat{\psi}\right)^2+\sum_{i=1}^2 k_{\mathrm{s}, i} \hat{u}_{\mathrm{s}, i}^2\right] \mathrm{d} x \big/ \\& \sum_{\rho=1}^N \int_{x_i^e}^{x_j^e}\left(\rho I_0 \hat{\psi}^2+\rho A_0 \hat{\omega}^2\right) \mathrm{d} x.\end{split} $

式中：符号“st”表示驻值.

与自振频率分析类似，屈曲荷载$ {P_{{\mathrm{cr}}}} $的变分方程可以写为

(32)$ \begin{split} P_{\mathrm{cr}}= & \mathrm{st} \sum_{e=1}^N \int_{x_i^e}^{x_j^e}\left[E I_0\left(\frac{\mathrm{d} \hat{\psi}}{\mathrm{d} x}\right)^2+\sum_{i=1}^3 E A_i\left(\frac{\mathrm{d} \hat{u}_i}{\mathrm{d} x}\right)^2+\right. \\& \left.C\left(\frac{\mathrm{d} \hat{w}}{\mathrm{d} x}-\hat{\psi}\right)^2+\sum_{i=1}^2 k_{\mathrm{s}, i} \hat{u}_{\mathrm{s}, i}^2\right] \mathrm{d} x\big/ \\&\sum_{e=1}^N \int_{x_i^e}^{x_j^e}(\mathrm{~d} \hat{w} / \mathrm{d} x)^2 \mathrm{d} x .\end{split} $

将式(23)代入式(31)，沿单元长度进行积分得到不含内部自由度的第e个单元的能量方程，根据最小势能原理可以得到

(33)$ ({{\boldsymbol{K}}^e} - {\omega ^2}{{\boldsymbol{M}}^e}){{\hat {\boldsymbol{u}}}^e} = {\bf{0}}. $

式中：$ {{\boldsymbol{M}}^e} $为质量矩阵，其元素表达式见附录C.

类似于自振频率的分析，通过对式(32)进行变分可以得到单元的几何刚度矩阵：

(34)$ ({{\boldsymbol{K}}^e} - {P_{{\mathrm{cr}}}}{\boldsymbol{K}}_{{P_A}}^e){{{\boldsymbol{u}}}^e} = {\bf{0}}. $

式中：$ {\boldsymbol{K}}_{{P_A}}^e $为几何刚度矩阵，其结果见附录D. 通过求解方程（34）可以得到各阶屈曲失稳临界荷载.

首先针对文献[9]中如图5所示的组合梁算例进行分析. 该梁跨径L₀=4 m，横截面宽度为0.1 m，梁高为0.24 m. 顶层和底层为钢，弹性模量为200 GPa，层高为0.02 m，泊松比为0.2；中间层为混凝土，弹性模量为34.5 GPa，层高为0.2 m，泊松比为0.3. 截面抗剪修正系数均为5/6. 剪切连接件抗剪刚度为k_s1=40 MPa，k_s2=k_s1/8=5 MPa.

简支和连续2种不同边界条件下计算所得的层间最大滑移、最大挠度及连续梁跨间最大弯矩结果分别如表1、2所示. 表中，δ₁为相对误差，δ₁=[单元数为100时本研究结果−文献[9]结果)/文献[9]结果]×100%. 当单元数量为100时，计算所得相对误差均不超过0.1%，显示了本研究计算结果的准确性. 在使用100个组合结构单元时2个算例的计算时间分别仅为0.184 7、0.309 3 s（Intel Core i5-11400CPU @ 2.60 GHz），且对于独立转角假设的含三截面转角模型^[12]，计算时间分别为0.661 5、0.795 8 s，本研究方法分别节省了258%、157%的计算时间，表明本研究方法的计算效率高. 对于简支梁算例，计算结果相对误差随着单元数的增加而变化的情况如图6所示. 图中，N_e为单元数量；δ₂为相对误差，δ₂=[(不同单元数下计算结果−单元数为500时计算结果)/单元数为500时结果]×100%. 计算结果显示，采用本研究方法建立的有限元模型收敛速度快，在4个单元时最大相对误差为0.03%；当单元数大于22时，相对误差低于0.00005%.

如图7(a)所示为简支边界时跨中最大挠度随层间抗剪刚度增大而变化的结果. 有限元模型中单元数为100，k_s1由0.1 MPa增至100 GPa，k_s2保持为k_s1/8. 结果显示，当抗剪刚度较小时（k_s1≤0.1 MPa），组合梁可以视为无抗剪刚度，跨中挠度超过14 mm；当抗剪刚度较大时（k_s1≥100 GPa），组合梁可视为完美抗剪连接，跨中挠度为2.84 mm. 如图7(b)所示为变形缩减系数$ \varsigma $（不同抗剪强度下变形计算结果与k_s1=0.1 MPa时变形计算结果的比值）随抗剪刚度增大的变化趋势. 计算结果显示存在有效抗剪刚度范围[1 MPa, 100 GPa]，可以使得夹层组合梁的变形快速降低.

当简支梁模型跨高比L₀/H=5、10和20时，采用Euler-Bernoulli梁（EB）与Timoshenko梁（TB）理论所得跨中挠度结果的相对误差随着抗剪刚度的变化如图8所示. 图中，δ₃为相对误差，δ₃=[(w_max,EB−w_max,TB)/w_max,TB]×100%，其中，w_{max, EB}、w_{max, TB}分别为EB、TB梁理论所得跨中挠度结果的最大值. 结果显示EB梁计算结果的相对误差随着抗剪刚度的增大而增大. 若以相对误差3%为界限，在低抗剪刚度（k_s1=0.1 MPa）情况下，当L₀/H=5时，EB梁计算结果不满足要求；当L₀/H=10时，抗剪刚度超过6.3 GPa，EB梁计算结果误差将超限；对于L₀/H=20的细长梁，EB梁计算结果相对误差不超过1%. 随着跨高比的减小，Euler-Bernoulli组合梁将不再适用，Timoshenko梁理论组合梁适用性更强.

如图9所示为简支梁与连续梁模型4种不同抗剪刚度时，采用Euler-Bernoulli梁（EB）与Timoshenko梁（TB）理论所得跨中挠度结果的相对误差随着跨高比的变化. 结果显示简支梁与连续梁在跨高比较小时，EB梁与TB梁计算结果的相对误差较大. 若以3%幅值为界限，当k_s1=1 MPa时，简支梁和连续梁的跨高比分别须大于5.4和8.8才可满足要求；当k_s1=10 GPa时，简支梁和连续梁的跨高比则分别须大于11.7和18.8才可满足要求.

简支边界时自振频率与屈曲荷载的计算结果如表3所示，其中ABAQUS模拟时采用平面应力有限元模型，δ₄为相对误差，δ₄=[(本研究结果−ABAQUS结果)/ ABAQUS结果]×100%. 结果显示本研究方法结果可与ABAQUS数值模拟结果吻合. 第1阶自振频率与屈曲荷载的相对误差均小于1%；前4阶模态时相对误差均小于3%. 如图10所示为模态1~5的自振频率与屈曲荷载的相对误差δ₅随单元数增加的变化，相对误差δ₅=[计算结果−单元数为500时结果)/单元数为500时结果]×100%. 对于第1阶自振频率与屈曲荷载，4个单元时相对误差即小于0.06%. 随着自振和屈曲失稳模态阶数的提高，所需单元数量增加. 对于第5阶自振频率与屈曲荷载，为了保证相对误差小于0.5%所需的单元数量分别为16和18.

本节所分析的连续边界夹层组合梁如图11所示，横截面宽度为0.1 m，截面高度共为0.8 m，3层的厚度分别为0.2、0.4、0.2 m. 3层材料均相同，弹性模量E=8 GPa，泊松比μ=0.3，截面抗剪修正系数κ=5/6，上下层间抗剪刚度相同，k_s1=k_s2.

当抗剪刚度k_s=1 MPa、20 MPa、50 MPa与10 GPa时，跨中挠度计算结果及其相对误差随跨高比的变化如图12所示. 图中，δ₆为相对误差，是相对于采用独立截面转角假设的含三截面转角的有限元模型^[12]的结果，δ₆=[(本研究计算结果−文献[12]计算结果)/文献[12]计算结果]×100%. 当跨高比大于8时，相对误差小于1.3%. 结果显示本研究所提单截面转角模型与三截面转角模型的计算结果基本吻合，对于常规跨高比的夹层组合梁，本研究方法具有足够的计算精度.

采用Timoshenko梁理论，引入含内部变量的高阶插值函数，推导3层组合梁的刚度矩阵、质量矩阵以及几何刚度矩阵元素的显式表达式. 利用MATLAB编译有限元程序，验证本研究方法的准确性和收敛性.

（1）通过算例分析表明，所提方法可以用于考虑界面滑移的夹层组合梁计算分析. 引入含内部变量高阶插值函数，避免了滑移锁定现象. 所提方法收敛速度快、所需单元少，简支梁静力算例中的单元数量取4时便有较高精度，相对误差仅为0.03%.

（2）当夹层梁跨高比相同时，Euler-Bernoulli梁的挠度计算结果相对于Timoshenko梁计算结果的误差随着抗剪刚度的增大而增大. 随着夹层梁跨高比的减小，采用基于Timoshenko梁理论的组合梁单元的必要性和优势逐渐凸显.

（3）相对于独立转角，假设的截面含多转角有限元模型，所推导的单元自由度更少，计算速度更快. 对于常规跨高比的各层材料性能没有较大突变的夹层梁，所提方法具有足够计算精度.

（4）所推导的3层组合梁的刚度矩阵、质量矩阵以及几何刚度矩阵元素的显式表达式具有计算效率高的优点，并可推广用于ABAQUS或ANSYS商业有限元软件程序中.

（5）本研究未考虑结合界面抗剪刚度的非线性力学行为，计算分析过程也未考虑材料非线性及几何非线性，相关问题值得进一步研究.