考虑滑移钢⁃混组合梁单元刚度矩阵推导的柔度法

Flexibility method for deriving element stiffness matrix of steel-concrete composite beams with slip

摘要

钢-混组合梁抗剪连接件剪切变形引起的界面滑移会减小截面弯曲刚度、增加挠度，对挠度要求严格的大跨度、轨道交通等结构不利，而有限元设计软件中的两节点梁单元还不能较好地考虑界面滑移。本文采用柔度法推导梁单元刚度矩阵。增加界面滑移自由度，采用梁段平衡微分方程推导平面简支梁四维柔度矩阵，柔度矩阵求逆得到同一节点位移-节点力的简支梁四维刚度矩阵，再进行几何关系转化，得到六维刚度矩阵。推导过程仅需组合轴力的二阶微分方程。均布荷载作用下，考虑滑移与不考虑滑移的杆端等效节点力几乎相同。应用推导的刚度矩阵自编有限元程序，通过了解析法、空间有限元数值模拟两种方法的验证。本文推导的解析型刚度矩阵在其假定范围内为精确解，较少的单元便能获得与解析解相同的结果。

Abstract

The interface slip caused by shear deformation of steel-concrete composite beam shear connectors will reduce the bending stiffness of the section and increase the deflection， which is unfavorable for large-span， rail transit and other structures with strict deflection requirements. However， the two-point beam element in finite element design software cannot effectively consider interface slip. The stiffness matrix of beam elements is derived using the flexibility method in this paper. Increasing the degree of freedom for interface slip， the four-dimensional flexibility matrix of a planar simply supported beam is derived using the equilibrium differential equations of the beam segment. The flexibility matrix is then inverted to obtain the four-dimensional stiffness matrix of the simply supported beam with the same node displacement and node force. The geometric relationship transformation is then performed to obtain the six-dimensional stiffness matrix. The derivation process only requires the second-order differential equation of the composite axial force. Under uniformly distributed loads， the equivalent nodal forces at the end of the rod considering slip and not considering slip are almost the same. The stiffness matrix derived from the application was self-written into a finite element program and was verified through analytical methods and spatial finite element numerical simulations. The derived analytical stiffness matrix is an exact solution within its assumed range， and the same results as the analytical solution can be obtained with fewer elements.

Graphical abstract

关键词

钢⁃混组合梁 / 有限元 / 柔度法 / 解析法 / 单元刚度 / 滑移

Key words

steel-concrete composite beam / finite element method / flexibility method / analytic method / element stiffness / slip

0 引 言

工程结构大多由一维构件组成，即便多维也常将其简化为一维，梁单元运行效率高、使用方便，是有限元分析应用最为广泛的单元类型之一。钢⁃混组合梁有3种剪切变形：①竖向剪切变形^［1］，由Timoshenko梁已能够较好地解决（如商业软件在单元定义中输入等效剪切面积，或勾选考虑剪切变形）；②混凝土板纵向剪力滞效应^［2］，由折减后的截面有效宽度简化考虑，桥梁钢-混组合梁按照现行规范大多不需要折减；③界面纵向剪切变形（界面滑移）^［3］，在有限元设计软件梁单元中还不能较好地考虑，如Midas Civil软件中需要用双单元分别模拟混凝土板、钢梁两个部分截面，界面用弹性连接模拟。考虑界面滑移会减小截面刚度，增加挠度，对挠度要求严格的大跨度、轨道交通等结构影响不能忽略。因此，对单独考虑界面滑移的钢-混组合梁，精度高、应用方便的单元刚度还亟须深入研究。

国内外学者对钢⁃混组合梁考虑滑移已开展大量有益研究。基于微分方程的解析法是主要方法^［4，5］，聂建国等^［6］对组合梁挠度给出了简化的折减刚度计算公式（成果已由相关规范所采纳），童根树^［7］利用桁架模型给出了等效刚度的一种简化公式，Zhou等^［8］利用弹簧串联模型给出了等效刚度的一种简化公式，冀伟等^［9］考虑了长期效应对挠度的影响，肖岩等^［10］推导了单跨超静定梁解析解，Nan等^［11］推导了连续组合梁解析解，Sun等^［12］分析了预应力连续组合梁，Francesco等^{［13，14］}给出了考虑滑移、竖向掀起等耦合的闭合解，刘世忠等^［15］基于能量变分原理推导了波形钢腹板-钢底板组合梁解析解。虽然解析法不具有通用性，条件不同解析解不同，但解析法研究为有限元分析提供了一种验证手段。

有限元分析单元刚度形成方法层出不穷^［16］。两节点刚度法梁单元应用方便，但钢-混组合梁考虑滑移后的两节点刚度法梁单元精度偏低，为了提高精度，孙飞飞等^［17］给出了三节点、四节点力法梁单元，周凌宇等^［18］给出了多自由度钢⁃混组合箱梁空间有限梁段法，Ataei^［19］建立了可拆卸连接件钢⁃混组合梁的有限元模型。柔度法是推导刚度矩阵的一种重要的方法^［20］，Faella等^［21］用虚功原理推导了柔度、刚度矩阵解析解，王文炜等^［22］用试验结果验证了所推导的刚度矩阵。解析型梁单元刚度矩阵精度高，需要的单元数量少。因此，本文单独考虑界面滑移，采用梁段平衡微分方程精确推导简支梁柔度矩阵、等效节点力解析解，进而获得解析型梁单元刚度矩阵，自编Matlab有限元程序，与解析法、Ansys空间有限元数值模拟结果吻合。

1 基本假定

（1）按照Newmark假定^［23］，钢⁃混组合梁连接界面上的抗剪连接件产生剪切变形，平截面假定不成立，但在混凝土板、钢梁两个部分截面上仍满足平截面假定，混凝土板与钢梁截面曲率相同。

（2）不考虑混凝土板竖向掀起，混凝土板、钢梁两个部分截面，无相互分离和嵌入。

（3）混凝土、钢、抗剪连接件为弹性材料，不考虑混凝土开裂。

（4）忽略混凝土板纵向剪切变形（剪力滞效应）。

（5）忽略钢梁腹板竖向剪切变形，梁体变形为小变形，转角为曲率沿纵向的一次积分，挠度为曲率的二次积分。

2 从柔度矩阵推导刚度矩阵

以平面梁单元为研究对象，选择两节点单元，有两个节点自由度（挠度w、转角θ）。考虑钢⁃混组合梁界面滑移后，需增加滑移自由度，可用混凝土板形心位置轴向位移、钢梁形心轴向位移两个自由度表示，也可用界面滑移（ξ）一个自由度。选择后者，则梁单元为两节点六自由度，单元刚度矩阵为6×6。节点力-节点位移设定如图１所示。有限元分析中采用的正负号规定为：挠度（剪力）向下为正，转角（弯矩）顺时针方向为正，界面滑移（钢梁顶部相对混凝土板底部位移差）向右为正，反之为负。

图1中，\overline{\xi }、\overline{w}、\overline{\theta }、分别为界面滑移、挠度、转角，{\overline{N}}_{\mathrm{i}}、\overline{V}、\overline{M}分别为组合轴力、剪力、弯矩，下标“1”、“2”分别表示节点1、2。组合轴力{\overline{N}}_{\mathrm{i}}为节点内力，{\overline{N}}_{\mathrm{i}}为分别作用于混凝土板形心、钢梁形心的一对力偶，而其对应的位移（界面滑移）位于混凝土板、钢梁连接界面位置。

单元刚度矩阵推导方法很多，采用刚度法，精度取决于位移形函数，如挠度三次多项式函数（Hermite插值形函数）满足平衡方程，为Euler⁃Bernoulli梁、Timoshenko梁的精确解。钢⁃混组合梁考虑滑移自由度后，刚度法的两节点单元精度偏差，一般用增加节点数量（即高次单元）方式提高精度。而用柔度法推导单元刚度矩阵，基本结构为静定结构，弯矩分布为线性，所导出的柔度矩阵在其假定范围内是精确解。

首先，对图1位移法的无约束梁单元施加约束，转化为柔度法基本体系（简支梁），节点力⁃节点位移见图2。其中，x为简支梁纵向长度坐标。

约束后，表示同一节点力-节点位移关系的简支梁刚度矩阵（ k_r）、简支梁柔度矩阵（ f_r）互为逆矩阵。节点位移表示为：

\left\{\begin{array}{l}{\xi }_{\mathrm{1}}\\ {\theta }_{\mathrm{1}}\\ {\xi }_{\mathrm{2}}\\ {\theta }_{\mathrm{2}}\end{array}\right\}={\mathit{f}}_{\mathrm{r}}\left\{\begin{array}{l}{N}_{\mathrm{i}\mathrm{1}}\\ M_{\mathrm{1}}\\ N_{\mathrm{i}\mathrm{2}}\\ M_{\mathrm{2}}\end{array}\right\}（1）

式中：ξ₁、θ₁分别为简支梁节点1的界面滑移、转角；ξ₂、θ₂分别为简支梁节点2的界面滑移、转角；N_i1、M₁分别为节点1的组合轴力、弯矩；N_i2、M₂分别为节点2的组合轴力、弯矩。

柔度矩阵求逆，得到简支梁的4×4刚度矩阵解析解。

比较简支梁单元与位移法无约束单元可知，简支梁少了刚体平移、刚度转角，如图3所示。因此进行位移的几何关系换算，如转角关系为：

{\theta }_{\mathrm{1}}=\frac{{\overline{w}}_{\mathrm{1}}-{\overline{w}}_{\mathrm{2}}}{L}+{\overline{\theta }}_{\mathrm{1}}（2）

式中：L为单元的长度。

简支梁单元与位移法无约束单元的节点位移换算关系为：

\left\{\begin{array}{l}{\xi }_{\mathrm{1}}\\ {\theta }_{\mathrm{1}}\\ {\xi }_{\mathrm{2}}\\ {\theta }_{\mathrm{2}}\end{array}\right\}=\underset{\mathit{S}}{\underbrace{\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{L}& \mathrm{1}& \mathrm{0}& -\frac{\mathrm{1}}{L}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{L}& \mathrm{0}& \mathrm{0}& -\frac{\mathrm{1}}{L}& \mathrm{1}\end{array}\right]}}\left\{\begin{array}{l}{\overline{\xi }}_{\mathrm{1}}\\ {\overline{w}}_{\mathrm{1}}\\ {\overline{\theta }}_{\mathrm{1}}\\ {\overline{\xi }}_{\mathrm{2}}\\ {\overline{w}}_{\mathrm{2}}\\ {\overline{\theta }}_{\mathrm{2}}\end{array}\right\}（3）

式中： S 为几何关系矩阵。

最后，得到6×6单元刚度矩阵（ k_e）：

{\mathit{k}}_{\mathrm{e}}={\mathit{S}}^{\mathrm{T}}{\mathit{k}}_{\mathrm{r}}\mathit{S}（4）

综上所述，任意梁单元的节点力⁃节点位移关系可由一个相应的简支梁表达出来，简支梁柔度矩阵求逆，得到简支梁刚度矩阵，再进行位移几何关系换算，得到位移法梁单元的刚度矩阵。

3 柔度法

3.1　梁段平衡微分方程

钢-混组合梁在刚性连接情况下（无界面滑移），弯矩（M）作用可以等效为混凝土板弯矩（M_c）、钢梁弯矩（M_s）、组合弯矩（N_id）三部分之和：

M_{\mathrm{c}}+M_{\mathrm{s}}+N_{\mathrm{i}}d=M（5）

其中： \left\{\begin{array}{l}{M}_{\mathrm{c}}=\frac{E_{\mathrm{c}}{I}_{\mathrm{c}}}{E{I}_{\mathrm{i}}}M,M_{\mathrm{s}}=\frac{E_{\mathrm{s}}{I}_{\mathrm{s}}}{E{I}_{\mathrm{i}}}M\\ E{I}_{\mathrm{i}}=EI+EA{d}^{\mathrm{2}},EI=E_{\mathrm{c}}{I}_{\mathrm{c}}+E_{\mathrm{s}}{I}_{\mathrm{s}}\\ N_{\mathrm{i}}d=\frac{EA{d}^{\mathrm{2}}}{E{I}_{\mathrm{i}}}M,\frac{\mathrm{1}}{EA}=\frac{\mathrm{1}}{E_{\mathrm{c}}{A}_{\mathrm{c}}}+\frac{\mathrm{1}}{E_{\mathrm{s}}{A}_{\mathrm{s}}}\end{array}\right.

式中：N_i为组合轴力；d为混凝土板形心至钢梁形心距离；E_c、E_s分别为混凝土、钢弹性模量；I_c、I_s分别为混凝土板、钢梁惯性矩；A_c、A_s分别为混凝土板、钢梁横截面面积；EI为部分截面叠加刚度；EA为部分截面轴向串联刚度；EI_i为组合截面刚度。

从式（5）可以看出，刚性连接情况下，组合截面内力计算只需要进行截面分析，部分截面内力直接按照刚度进行分配，组合后截面刚度（EI_i）实现了“1+1>2”。

考虑界面滑移，部分截面内力无法由截面分析（式（5））直接求解，需解微分方程，并结合梁的纵向边界条件给出定解常数，而纵向边界条件与荷载类型与分布、支座、简支还是连续梁等许多因素有关，不同梁具有不同的解析解，解析法不具有通用性。

按照Newmark假定条件下的梁段分析（见图4），将离散分布抗剪连接件等效为连续分布，可推导出组合轴力二阶微分方程：

{N_{\mathrm{i}}}^{\mathrm{\text{'}}\mathrm{\text{'}}}-{\alpha }^{\mathrm{2}}{N}_{\mathrm{i}}=-\frac{k_{\mathrm{s}}d}{EI}M（6）

式中：{\alpha }^{\mathrm{2}}=k_{\mathrm{s}}\left(\frac{d^{\mathrm{2}}}{EI}+\frac{\mathrm{1}}{EA}\right)，k_{\mathrm{s}}=\frac{n_{\mathrm{s}}{K}_{\mathrm{s}}}{a}；α为反映组合梁材料、几何特性的参数；n_s为连接件排数；K_s为单个连接件抗剪刚度；a为连接件间距。

分析可知，如果采用直接刚度法推导刚度矩阵，难度在于确定无约束单元的边界条件，而柔度法在简支梁条件下则容易给出。

按照Newmark假定条件，截面曲率κ为：

\kappa =\frac{M_{\mathrm{c}}+M_{\mathrm{s}}}{EI}（7）

将式（5）代入，有：

\kappa =\frac{M}{EI}-\frac{N_{\mathrm{i}}d}{EI}（8）

有了曲率，按照小变形假定，转角与挠度为：

\left\{\begin{array}{l}{\theta }^{\text{'}}=-\kappa \\ w^{\text{'}}=\theta \end{array}\right.（9）

需要注意的是，梁段分析遵循的正负号规定与有限元不同，在计算过程中应注意两套规定的转化关系。对梁段分析，本文按照轴力受拉为正，剪力以绕脱离体顺时针方向为正，弯矩以梁体下侧受正，曲率以逆时针方向为正，截面高度坐标以下为正。

3.2　柔度矩阵解析解

3.2.1　柔度矩阵第一列

按照柔度矩阵定义，第一列元素为基本体系（见图2）在N_{\mathrm{i}\mathrm{1}}=\mathrm{1}、其他节点力为0作用下节点的位移（节点界面滑移、节点转角），此时M=\mathrm{0}，式（6）改写为：

{N_{\mathrm{i}}}^{\mathrm{\text{'}}\mathrm{\text{'}}}-{\alpha }^{\mathrm{2}}{N}_{\mathrm{i}}=\mathrm{0}（10）

由边界条件：

x=\mathrm{0}，N_{\mathrm{i}}=\mathrm{1}；x=L，N_{\mathrm{i}}=\mathrm{0}

得到组合轴力解析解：

N_{\mathrm{i}}=\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}+\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}（11）

另外，按照梁段平衡关系，剪力p为：

p=-N_{\mathrm{i}}^{\text{'}}（12）

界面抗剪本构关系为：

\xi =\frac{p}{k_{\mathrm{s}}}（13）

得到界面滑移解析解：

\xi =\frac{\alpha }{k_{\mathrm{s}}}\left(\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}-\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}\right)（14）

则简支梁的节点界面滑移为：

\left\{\begin{array}{l}{\xi }_{\mathrm{1}}=\frac{\alpha }{k_{\mathrm{s}}}\frac{{\mathrm{e}}^{\alpha L}+{\mathrm{e}}^{-\alpha L}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\\ {\xi }_{\mathrm{2}}=\frac{\alpha }{k_{\mathrm{s}}}\frac{\mathrm{2}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\end{array}\right.（15）

由式（8）（11），得到曲率解析解：

\kappa =-\frac{d}{EI}\left(\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}+\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}\right)（16）

再结合式（9）及边界条件：

x=\mathrm{0}，w=\mathrm{1}；x=L，w=\mathrm{0}

得到转角解析解：

\theta =\frac{d}{EI{\alpha }^{\mathrm{2}}L}+\frac{d}{EI\alpha }\left(\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}-\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}\right)（17）

则简支梁的节点转角为：

\left\{\begin{array}{l}{\theta }_{\mathrm{1}}=\frac{L}{EI}\frac{d}{\alpha L}\left(\frac{\mathrm{1}}{\alpha L}-\frac{{\mathrm{e}}^{\alpha L}+{\mathrm{e}}^{-\alpha L}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\right)\\ {\theta }_{\mathrm{2}}=\frac{L}{EI}\frac{d}{\alpha L}\left(\frac{\mathrm{1}}{\alpha L}-\frac{\mathrm{2}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\right)\end{array}\right.（18）

3.2.2　柔度矩阵第二列

同理，第二列元素为M_{\mathrm{1}}=\mathrm{1}条件下节点的位移，此时M=\mathrm{1}-\frac{x}{L}，式（6）改写为：

{N_{\mathrm{i}}}^{\mathrm{\text{'}}\mathrm{\text{'}}}-{\alpha }^{\mathrm{2}}{N}_{\mathrm{i}}=-\frac{k_{\mathrm{s}}d}{EI}+\frac{x}{L}\frac{k_{\mathrm{s}}d}{EI}（19）

由边界条件：

x=\mathrm{0}，N_{\mathrm{i}}=\mathrm{1}；x=L，N_{\mathrm{i}}=\mathrm{0}

得到组合轴力解析解：

N_{\mathrm{i}}=\frac{k_{\mathrm{s}}d}{{\alpha }^{\mathrm{2}}EI}-\frac{x}{L}\frac{k_{\mathrm{s}}d}{{\alpha }^{\mathrm{2}}EI}-

\frac{k_{\mathrm{s}}d}{{\alpha }^{\mathrm{2}}EI}\left(\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}-\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}\right)（20）

界面滑移解析解：

\xi =\frac{d}{\alpha EI}\left(\frac{\mathrm{1}}{\alpha L}+\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}-\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}\right)（21）

节点界面滑移：

\left\{\begin{array}{l}{\xi }_{\mathrm{1}}=\frac{L}{EI}\frac{d}{\alpha L}\left(\frac{\mathrm{1}}{\alpha L}-\frac{{\mathrm{e}}^{\alpha L}+{\mathrm{e}}^{-\alpha L}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\right)\\ {\xi }_{\mathrm{2}}=\frac{L}{EI}\frac{d}{\alpha L}\left(\frac{\mathrm{1}}{\alpha L}-\frac{\mathrm{2}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\right)\end{array}\right.（22）

式（22）与式（18）结果相同，与柔度矩阵对称特性一致，说明2个计算公式推导正确。

曲率解析解为：

\kappa =\frac{\mathrm{1}}{EI}\left(\mathrm{1}-\frac{x}{L}\right)-\frac{k_{\mathrm{s}}{d}^{\mathrm{2}}}{{\alpha }^{\mathrm{2}}E{I}^{\mathrm{2}}}\left(\mathrm{1}-\frac{x}{L}\right)+

\frac{k_{\mathrm{s}}{d}^{\mathrm{2}}}{{\alpha }^{\mathrm{2}}E{I}^{\mathrm{2}}}\left(\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}+\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}\right)（23）

再结合式（9）及边界条件：

x=\mathrm{0}，w=\mathrm{1}；x=L，w=\mathrm{0}

得到转角解析解：

\theta =\frac{L}{\mathrm{3}EI}\left(\mathrm{1}-\frac{\mathrm{3}x}{L}+\frac{\mathrm{3}{x}^{\mathrm{2}}}{\mathrm{2}{L}^{\mathrm{2}}}\right)-\frac{k_{\mathrm{s}}{d}^{\mathrm{2}}}{{\alpha }^{\mathrm{2}}E{I}^{\mathrm{2}}}\left(\frac{x^{\mathrm{2}}}{\mathrm{2}L}-x+\frac{L}{\mathrm{3}}\right)-\frac{k_{\mathrm{s}}{d}^{\mathrm{2}}}{{\alpha }^{\mathrm{4}}E{I}^{\mathrm{2}}L}-\frac{k_{\mathrm{s}}{d}^{\mathrm{2}}}{{\alpha }^{\mathrm{3}}E{I}^{\mathrm{2}}}\left(\frac{{\mathrm{e}}^{\alpha x}}{\mathrm{1}-{\mathrm{e}}^{\mathrm{2}\alpha L}}-\frac{{\mathrm{e}}^{-\alpha x}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\alpha L}}\right)（24）

节点转角为：

\left\{\begin{array}{l}{\theta }_{\mathrm{1}}=\frac{L}{\mathrm{3}E{I}_{\mathrm{i}}}+\frac{L}{E{I}_{\mathrm{i}}}\frac{\beta }{\alpha L}\left(\frac{{\mathrm{e}}^{\alpha L}+{\mathrm{e}}^{-\alpha L}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}-\frac{\mathrm{1}}{\alpha L}\right)\\ {\theta }_{\mathrm{2}}=-\frac{L}{\mathrm{6}E{I}_{\mathrm{i}}}+\frac{L}{E{I}_{\mathrm{i}}}\frac{\beta }{\alpha L}\left(\frac{\mathrm{2}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}-\frac{\mathrm{1}}{\alpha L}\right)\end{array}\right.（25）

式中：刚度系数\beta =\frac{EA{d}^{\mathrm{2}}}{EI}。

3.2.3　柔度矩阵

由式（15）（18）（22）（25）及柔度矩阵对称特性，可以写出简支梁的4×4柔度矩阵（见式（26））。采用梁段平衡微分方程推导的柔度矩阵与虚功原理推导的结果相同。

柔度矩阵求逆，得到简支梁的刚度矩阵。由此得到的刚度矩阵在其假定范围内为精确解。其显示表达式冗长、复杂，因此直接用矩阵求逆函数，在计算程序中刚度矩阵采用隐式表达。

\begin{array}{c}\left\{\begin{array}{l}{N}_{\mathrm{i}\mathrm{1}}\\ M_{\mathrm{1}}\\ N_{\mathrm{i}\mathrm{2}}\\ M_{\mathrm{2}}\end{array}\right\}=\\ \underset{\mathrm{柔度}\mathrm{矩阵}{\mathit{f}}_{\mathrm{r}}}{\underbrace{\left[\begin{array}{cccc}\frac{\alpha }{k_{\mathrm{s}}}\frac{{\mathrm{e}}^{\alpha L}+{\mathrm{e}}^{-\alpha L}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}& f_{\mathrm{12}}=f_{\mathrm{21}}& f_{\mathrm{13}}=f_{\mathrm{31}}& f_{\mathrm{14}}=f_{\mathrm{41}}\\ T_{\mathrm{1}}& T_{\mathrm{2}}& f_{\mathrm{23}}=f_{\mathrm{41}}& f_{\mathrm{24}}=f_{\mathrm{42}}\\ \frac{\alpha }{k_{\mathrm{s}}}\frac{\mathrm{2}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}& f_{\mathrm{32}}=f_{\mathrm{41}}& f_{\mathrm{33}}=f_{\mathrm{11}}& f_{\mathrm{34}}=f_{\mathrm{21}}\\ T_{\mathrm{3}}& T_{\mathrm{4}}& f_{\mathrm{43}}=f_{\mathrm{21}}& f_{\mathrm{44}}=f_{\mathrm{22}}\end{array}\right]}}\left\{\begin{array}{l}{\xi }_{\mathrm{1}}\\ {\theta }_{\mathrm{1}}\\ {\xi }_{\mathrm{2}}\\ {\theta }_{\mathrm{2}}\end{array}\right\}\end{array}（26）

式中：T₁=\frac{L}{EI}\frac{d}{\alpha L}\left(\frac{\mathrm{1}}{\alpha L}-\frac{{\mathrm{e}}^{\alpha L}+{\mathrm{e}}^{-\alpha L}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\right)；T₂=\frac{L}{\mathrm{3}E{I}_{\mathrm{i}}}+\frac{L}{E{I}_{\mathrm{i}}}\frac{\beta }{\alpha L}\left(\frac{{\mathrm{e}}^{\alpha L}+{\mathrm{e}}^{-\alpha L}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}-\frac{\mathrm{1}}{\alpha L}\right)；T₃=\frac{L}{EI}\frac{d}{\alpha L}\times \left(\frac{\mathrm{1}}{\alpha L}-\frac{\mathrm{2}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}\right)；T₄=\frac{L}{E{I}_{\mathrm{i}}}\frac{\beta }{\alpha L}\left(\frac{\mathrm{2}}{{\mathrm{e}}^{\alpha L}-{\mathrm{e}}^{-\alpha L}}-\frac{\mathrm{1}}{\alpha L}\right)-\frac{L}{\mathrm{6}E{I}_{\mathrm{i}}}；f为柔度，下标为“行号+列号”。

3.3　均布荷载作用下等效节点力

按照等效节点力的定义，考虑简支梁单元支座反力对位移法无约束梁单元的修正，由简支梁等效节点力（ Q_r）可以得到等效节点力（ Q ）：

\mathit{Q}={\mathit{S}}^{\mathrm{T}}{\mathit{Q}}_{\mathrm{r}}+{\left\{\begin{array}{cccccc}\mathrm{0}& \frac{qL}{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \frac{qL}{\mathrm{2}}& \mathrm{0}\end{array}\right\}}^{\mathrm{T}}（27）

式中：{\mathit{Q}}_{\mathrm{r}}={\mathit{f}}_{\mathrm{r}}^{-\mathrm{1}}{\mathit{\delta }}_{\mathrm{r}}，其中， δ_r为均布荷载作用下的节点位移。

同理， δ_r由梁段平衡微分方程推导有：

{\mathit{\delta }}_{\mathrm{r}}=\left\{\begin{array}{l}-\frac{qd}{EI{\alpha }^{\mathrm{3}}}\left(\frac{\mathrm{1}-{\mathrm{e}}^{\alpha L}}{\mathrm{1}+{\mathrm{e}}^{\alpha L}}+\frac{\alpha L}{\mathrm{2}}\right)\\ \frac{q{L}^{\mathrm{3}}}{\mathrm{24}E{I}_{{}^{\mathrm{i}}}}+\frac{q{L}^{\mathrm{3}}}{E{I}_{{}^{\mathrm{i}}}}\frac{\beta }{{\alpha }^{\mathrm{3}}{L}^{\mathrm{3}}}\left(\frac{\mathrm{1}-{\mathrm{e}}^{\alpha L}}{\mathrm{1}+{\mathrm{e}}^{\alpha L}}+\frac{\alpha L}{\mathrm{2}}\right)\\ \frac{qd}{EI{\alpha }^{\mathrm{3}}}\left(\frac{\mathrm{1}-{\mathrm{e}}^{\alpha L}}{\mathrm{1}+{\mathrm{e}}^{\alpha L}}+\frac{\alpha L}{\mathrm{2}}\right)\\ -\frac{q{L}^{\mathrm{3}}}{\mathrm{24}E{I}_{{}^{\mathrm{i}}}}-\frac{q{L}^{\mathrm{3}}}{E{I}_{{}^{\mathrm{i}}}}\frac{\beta }{{\alpha }^{\mathrm{3}}{L}^{\mathrm{3}}}\left(\frac{\mathrm{1}-{\mathrm{e}}^{\alpha L}}{\mathrm{1}+{\mathrm{e}}^{\alpha L}}+\frac{\alpha L}{\mathrm{2}}\right)\end{array}\right\}（28）

经过分析与验证，均布荷载作用下，考虑滑移与不考虑滑移的等效节点力几乎相同，即：

\mathit{Q}={\left\{\begin{array}{cccccc}\mathrm{0}& \frac{qL}{\mathrm{2}}& -\frac{q{L}^{\mathrm{2}}}{\mathrm{12}}& \mathrm{0}& \frac{qL}{\mathrm{2}}& \frac{q{L}^{\mathrm{2}}}{\mathrm{12}}\end{array}\right\}}^{\mathrm{T}}（29）

4 算　例

选择均布荷载作用下简支组合梁为研究对象，长度10 m，均布荷载50 kN/m，混凝土弹性模量3×10⁴ N/mm²，钢弹性模量2.1×10⁵ N/mm²。钢-混组合截面如图5所示。抗剪连接件（栓钉）单排布置，间距100 mm，栓钉抗剪刚度K_{\mathrm{s}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}\mathrm{N}/\mathrm{m}\mathrm{m}。无量纲参数\alpha L=\mathrm{8.4}，\beta =\mathrm{1.723}。

如果不考虑界面滑移（刚性连接），跨中挠度为\frac{\mathrm{5}q{L}^{\mathrm{4}}}{\mathrm{384}E{I}_{\mathrm{i}}}=\mathrm{32.25}\mathrm{m}\mathrm{m}。

考虑界面滑移，使得钢⁃混组合梁的组合轴力减小，截面刚度减小，挠度增大。

采用3种方法计算考虑界面滑移的挠度，可相互验证。

4.1　方法一：解析法

按照梁段平衡微分方程，推导均布荷载作用下简支梁跨中挠度w_m，计算公式为：

w_{\mathrm{m}}=

\frac{\mathrm{5}q{L}^{\mathrm{4}}}{\mathrm{384}E{I}_{\mathrm{i}}}\left\{\mathrm{1}+\frac{\mathrm{384}\beta }{\mathrm{5}{\alpha }^{\mathrm{4}}{L}^{\mathrm{4}}}\left[\frac{{\alpha }^{\mathrm{2}}{L}^{\mathrm{2}}}{\mathrm{8}}-\mathrm{1}+\mathrm{c}\mathrm{s}\mathrm{c}\mathrm{h}\left(\frac{\alpha L}{\mathrm{2}}\right)\right]\right\}（30）

从式（30）可以看出，考虑滑移的挠度是在不考虑滑移挠度基础上乘以一个扩大系数，该系数与无量纲参数β、αL有关。

代入已知条件，w_{\mathrm{m}}=\mathrm{38.22}\mathrm{m}\mathrm{m}。考虑界面滑移的扩大系数为1.185。

4.2　方法二：数值模拟

用Ansys建立了一个精细的空间有限元模型，混凝土板采用Solid45实体单元（尺寸x×y×z=100 mm×37.5 mm×220 mm，单元数量为4 000个），钢梁翼板采用Solid45实体单元（单元尺寸100 mm×15 mm×100 mm，单元数量为400个），钢梁腹板采用Plane42平面单元（单元尺寸100 mm×92.5 mm，单元数量为800个）。栓钉采用Cimbine39非线性弹簧单元，忽略了混凝土板与钢梁界面上的粘结力、摩擦力，假定仅有栓钉来抵抗界面滑移。非线性分析采用自动时间步。

跨中挠度计算结果为38.32 mm，位移云图如图6所示。空间有限元数值模拟结果略大于解析法，原因在于解析法未考虑钢梁腹板的竖向剪切变形，而钢梁腹板Plane42平面单元考虑了剪切变形。由此说明结果合理，所建立的空间有限元数值模型正确。

4.3　方法三：自编有限元程序

应用本文两节点梁单元，按照弹性有限元程序既定模块化编制流程，经过调试、测试（包括两套正负号规定的转换等问题），获得了正确的Matlab有限元程序。

考虑滑移钢⁃混组合梁跨中挠度、单元数量等结果，如表1所示。从表1可以看出，随着单元数量增加，结果趋于收敛，当单元数量为4，与解析解偏差几乎为0，说明推导的梁单元解析型刚度矩阵精度较高，较少单元便能得到与解析解相同的结果，而空间有限元数值模拟则需要3种单元类型及庞大的单元数量。

3种方法计算结果吻合，说明本文两节点梁单元通过了解析法、空间有限元数值模拟的验证，单元刚度矩阵结果正确。

另外，需要说明的是，一般情况下解析型刚度矩阵仅需一个单元变可获得与解析解完全相同的结果。本算例求跨中挠度，则必须在跨中设置节点，故在此情况下最小单元数量为2个。表1在单元数量2时出现了偏差0.7%，原因是等效节点力用了近似解式（29）；若用式（27），则偏差为0。

5 结束语

按照梁段平衡微分方程，推导了简支梁节点位移-节点力关系柔度矩阵解析解，柔度矩阵求逆获得刚度矩阵。对简支梁4×4刚度矩阵进行几何关系换算，获得两节点梁单元6×6刚度矩阵。通过计算与分析，得出以下结论：①通过解析法、空间有限元数值模拟验证，由柔度法推导的解析型单元刚度矩阵在其假定范围内为精确解，较少的单元便能获得与解析解相同的结果；②均布荷载作用下，考虑滑移与不考虑滑移的杆端等效节点力几乎相同；③从梁段平衡微分方程推导过程可以看出，仅需使用组合轴力的二阶微分方程，不需要使用更高阶的微分方程。后续研究将在滑移的基础上考虑混凝土板竖向掀起等因素的影响。

参考文献

[1] Zhao X, Shao X, Cao J, et al.Experimental studies on shear behavior of steel-UHPC composite beam with hot rolled shape steel[J]. Engineering Structures,2023, 274(1): No.115160.

[2] He S, Yang G, Zhou W, et al.Evaluation of shear lag effect in HSS-UHPC composite beams with perfobond strip connectors: experimental and numerical studies[J].Journal of Constructional Steel Research, 2022, 194(7): No.107312.

[3] 高策, 解会兵, 韩冰, 等. 疲劳荷载作用下钢-混组合梁界面滑移对其性能影响研究[J]. 桥梁建设, 2024, 54(4): 28-36.

[4] Gao Ce, Hui-bing Jie, Han Bing, et al. Study on the influence of interface slip on the performance of steel-concrete composite beams under fatigue load[J]. Bridge Construction, 2024, 54(4): 28-36.

[5] Ge H, Dai G, Liu W W,et al. Theoretical solution for bond-slip behavior of composite structures consisting of H-shape beam and concrete based on experiment, numerical simulation, and theoretical derivation[J]. Engineering Structures, 2024, 302(1):No.117456.

[6] Peng S, Zhu Z, Wei Y. An analytic solution for bending of multilayered structures with interlayer-slip[J]. International Journal of Mechanical Sciences, 2024, 282(15): No.109642.

[7] 聂建国, 沈聚敏, 余志武. 考虑滑移效应的钢-混凝土组合梁变形计算的折减刚度法[J]. 土木工程学报, 1995, 28(6): 11-17.

[8] Nie Jian-guo, Shen Ju-min, Yu Zhi-wu. Reduced stiffness method for deformation calculation of steel-concrete composite beams considering slip effect[J]. Journal of Civil Engineering,1995,28(6): 11-17.

[9] 童根树. 钢-混凝土组合梁的刚度[J]. 钢结构, 2022, 37(10): 49-52.

[10] Tong Gen-shu. Stiffness of steel-concrete composite beams [J]. Steel Structures,2022,37(10):49-52.

[11] Zhou Dong-hua, Matthias Pahn, Wolfgang Kurz. Contribution to the calculation of composite beams with elastic composite action[J]. Stahlbau, 2014, 83(4): 278-288.

[12] 冀伟, 白倩, 罗奎. 长期效应下不同边界钢-混凝土组合梁挠度分析[J]. 哈尔滨工程大学学报, 2023, 44(3): 434-442.

[13] Ji Wei, Bai Qian, Luo Kui. Analysis of deflection of steel-concrete composite beams with different boundaries under long-term effects[J]. Journal of Harbin Engineering University, 2023, 44(3): 434-442.

[14] 肖岩, 彭罗文, Kunnath S. 组合梁考虑滑移效应的理论分析[J]. 湖南大学学报:自然科学版, 2017, 44(1):77-86.

[15] Xiao Yan, Peng Luo-wen, Kunnath S. Theoretical analysis of composite beams considering slip effects[J]. Journal of Hunan University(Natural Science Edition), 2017, 44(1): 77-86.

[16] Nan H, Wang P, Zhang Q, et al. Study on the mechanical properties of continuous composite beams under coupled slip and creep[J]. Materials, 2023, 16(13): No.4741.

[17] Sun Jun, Yue Zhen-hong, He Yan,et al. Slip analysis of prestressed steel-concrete continuous composite beam[J]. Journal of King Saud University(Engineering Sciences), 2024, 36(8): 535-542.

[18] Francesco S A, Leo K, Gordan J. Closed-form solutions for two-layer Timoshenko beams with interlayer slip, uplift and rotation compliance[J]. Meccanica: Journal of the Italian Association of Theoretical and Applied Mechanics, 2023, 58(8): 893-918.

[19] Oeng T, Keo P, Guezouli S,et al. Large displacement analysis of two-layer beam-columns taking into account slip and uplift[J]. Engineering computations: International Journal for Computer-aided Engineering and Software, 2023, 40(1): 265-295.

[20] 刘世忠, 毛亚娜, 马驰, 等. 单箱双室波形钢腹板-钢底板组合梁滑移效应研究[J]. 桥梁建设, 2023, 53(4): 47-54.

[21] Liu Shi-zhong, Mao Ya-na, Ma Chi, et al. Study on slip effect of single box double chamber corrugated steel web plate steel bottom plate composite beam[J]. Bridge Construction, 2023, 53(4): 47-54.

[22] Fortuna B, Turk G, Schnabl S. A new locking-free finite element for N-layer composite beams with interlayer slips and finger joints[J]. Finite Elements in Analysis and Design, 2023, 220(8): No.103936.

[23] 孙飞飞, 李国强. 考虑剪力滞后效应的钢-混凝土组合梁力法单元[J]. 力学季刊, 2004, 25(1): 78-88.

[24] Sun Fei-fei, Li Guo-qiang. Force method element for steel-concrete composite beams considering shear lag effect[J]. Mechanics Quarterly, 2004, 25(1): 78-88.

[25] 周凌宇, 贺桂超. 考虑界面滑移的多自由度钢-混组合箱梁空间有限梁段法[J]. 铁道科学与工程学报, 2012, 9(4): 1-7.

[26] Zhou Ling-yu, He Gui-chao. Space finite beam section method for multi degree of freedom steel-concrete composite box girder considering interface slip[J]. Journal of Railway Science and Engineering, 2012, 9(4): 1-7.

[27] Ataei A. Modelling of demountable steel-concrete composite connections: validation of finite element model and parametric study[J]. Journal of Constructional Steel Research, 2022, 198(Sup.C): No.107585.

[28] Hosseini S M, Dehcheshmeh M M, Amiri G G.Experimental and numerical structural damage detection using a combined modal strain energy and flexibility method[J]. Structural Engineering and Mechanics, 2023, 87(6): 555-574.

[29] Faella C, Martinelli E, Nigro E. Steel and concrete composite beams with flexible shear connection:"exact" analytical expression of the stiffness matrix and applications[J].Computers & Structures, 2002, 80(11): 1001-1009.

[30] 王文炜, 翁昌年, 万水, 等. 考虑滑移的钢-混凝土组合梁有限单元法[J]. 东南大学学报: 自然科学版, 2007, 37(2): 206-212.

[31] Wang Wen-wei, Weng Chang-nian, Wan Shui, et al. Finite element method for steel-concrete composite beams considering slip[J]. Journal of Southeast University(Natural Science Edition), 2007, 37(2):206-212.

[32] Newmark N M, Siess C P, Viest I M.Tests and analysis of composite beams with incomplete interaction[J]. Society for Experimental Stress Analysis, 1951, 9(1): 75-92.

基金资助

国家自然科学基金项目(51868034)