预应力混凝土梁桥由于具有良好的工作性能在我国得到了广泛地使用，为更好地了解预应力梁桥在正常使用阶段及破坏阶段的工作性能，很多学者基于各种有限元模型对其进行了非线性受力过程分析。目前已有文献中的有限元模型主要有分层梁单元、梁段单元、退化板壳单元及体单元等。对于预应力混凝土T梁桥而言，上述有限元模型各有优缺点，难分伯仲；但对于预应力混凝土箱梁桥而言，由于截面为典型的薄壁结构，综合比较而言，壳单元模型更为适用^[1]。目前已有的采用壳单元模型进行预应力混凝土梁桥非线性分析的研究文献^{[2, 3, 4, 5, 6, 7, 8]}，对于混凝土及普通钢筋的模拟，一般都是采用分层实体退化壳单元来模拟，各模型的差别主要在于预应力钢筋特别是曲线预应力钢筋的模拟上。文献^[2]处理预应力钢筋的方法与普通钢筋一样，将其视为分层实体退化壳单元的一层。文献^[3]则按预应力钢筋的几何形状是平直还是弯曲分别考虑，对于前者，处理方法与文献^[2]相同，也是将预应力钢筋视为分层壳单元中的一层；对于后者，介绍较笼统，只称用组合壳单元模型模拟弯曲预应力钢筋与周围混凝土，用虚功原理导出弯曲预应力钢筋对组合壳单元刚度矩阵的贡献，但并未详细介绍弯曲预应力钢筋本身的非线性刚度矩阵怎么导出。文献^[4]则基于全拉格朗日列式导出的三维杆单元来模拟所有预应力钢筋，混凝土梁及与预应力筋结合的模拟则与文献^[3]完全相同。文献^{[5, 6, 7, 8]}的单元模型与文献^{[2, 3, 4]}完全相同，只是在预应力梁具体性能的研究方面有一定区别，这些研究成果均有助于了解预应力梁在非线性阶段的工作性能。在非线性分析中，基于全拉格朗日列式或增量拉格朗日列式进行非线性有限元推导时，材料本构关系与单元的运动描述是紧密耦合的，不仅需要专门的应力更新算法(即使对于线弹性材料也不例外)，还需要不断更新单元形函数，导致有限元列式复杂(因此上述大部分文献对预应力钢筋都是采用理想弹塑性等较简单的材料特性)，计算效率也不高^{[9, 10, 11]}。鉴于此，本文基于随转坐标法和场一致性原则，导出能考虑几何及材料双非线性的预应力钢筋空间杆单元，该单元具有几何与材料两种非线性不耦合的特点。采用已有的实体退化壳单元模拟混凝土，根据预应力钢筋空间杆单元与混凝土实体退化壳单元在单元内的位移协调条件和虚功原理将两者组合成一个混合壳单元，基于该混合壳单元模型编制了预应力钢筋混凝土梁非线性分析程序，对已有文献介绍的算例进行分析，并用其他文献的计算结果及试验结果进行验证。

本文的实体退化壳单元如图 1所示，为四边形八节点等参壳单元^{[12, 13]}。

图 1 四边形等参壳体单元 Fig. 1 quadrilateral isoparametric shell

图选项

在结构坐标系下壳单元每个节点上有5个自由度，分别为3个位移自由度(u_k^中,v_k^中,w_k^中)，及两个节点处垂直于中面法线的旋转自由度(β_1k,β_2k)，单元位移场u，v，w可表示为：

点击浏览公式

为计入几何非线性，引入Von karman假设，基于全拉格朗日列式来描述几何非线性，可得到单元的应变位移矩阵B由线性部分B₀和非线性部分B_L组成，即：

点击浏览公式

为考虑材料非线性，实体退化壳单元引入文献^[12]的分层模型，采用由混凝土层和等效厚度的均匀钢层组成的分层模型来计算单元切线刚度矩阵，其中模拟普通钢筋的钢层只具有轴向强度和刚度，并假设混凝土层和钢层之间无相对滑移。

在进行了上述几何与材料非线性分析后，由虚功原理可得到混凝土及普通钢筋分层壳元包含几何与材料双非线性的切线刚度矩阵K^cs(上标‘cs’代表混凝土层及普通钢筋层，以下均同)为：

点击浏览公式

显然，混凝土及普通钢筋分层壳元等效节点力F^cs的计算公式为^[13]：

点击浏览公式

在计算混凝土及普通钢筋分层壳元的切线刚度矩阵K^cs及等效节点力F^cs时所需的钢筋物理特性、混凝土的本构关系及开裂与拉伸刚化作用如何考虑等均可见文献^[12]。

壳单元中带有一根预应力钢筋如图 2所示，假设杆元端点与周围混凝土完全黏结，也就是说杆元端点位移和周围混凝土的位移是完全相同。为导出结构坐标系下预应力钢筋对混合壳单元切线刚度矩阵的贡献，分两步进行推导：第一步基于随转坐标法和场一致性原则，导出结构坐标系下考虑初应变的预应力钢筋空间杆元几何与材料双非线性切线刚度矩阵；第二步根据预应力钢筋空间杆元节点在壳单元中的坐标，由壳单元内的位移协调条件和虚功原理，导出结构坐标系下预应力钢筋对混合壳单元切线刚度矩阵的贡献。

图 2 带预应力钢筋的四边形等参壳体单元 Fig. 2 quadrilateral isoparametric shell element

图选项

图 3为初始时刻与计算t时刻的空间杆元ij。设结构坐标系为xyz，设单元的随转坐标系为XYZ，该坐标系是随单元变形而转动的，它始终以节点i为原点，以节点i到j的连线方向为X轴，由节点i与无穷远处点形成Y轴，再通过右手法则形成Z轴(图中未示意出随转坐标系XYZ)。

图 3 各时刻的空间杆单元 Fig. 3 Spatial truss element at different time

图选项

设初始时刻杆单元节点在结构坐标系里的坐标为(⁰x_i,⁰y_i,⁰z_i)和(⁰x_j,⁰y_j,⁰z_j)，在荷载作用下计算t时刻结构坐标系中的的坐标为(^tx_i,^ty_i,^tz_i)和(^tx_j,^ty_j,^tz_j)；节点的位移向量为d=[u_i v_i w_i u_j v_j w_j]^T，有：

点击浏览公式

设初始时刻及计算t时刻杆单元长度分别为⁰l,^tl,显然有：

点击浏览公式

从初始时刻到计算t时刻杆单元的轴向伸长量为：

点击浏览公式

计算t时刻杆单元与结构坐标系的x，y，z轴的方向余弦值分别为：

点击浏览公式

对式(7)微分，可得到δΔl用结构坐标系下的位移微分δd=[δu_iδv_iδw_iδu_jδv_jδw_j]^T表示为：

点击浏览公式

计算t时刻杆单元的应变值为ε=Δl/⁰l+ε₀(ε₀为预应力张拉引起的初应变)，通过预应力钢筋的应力-应变关系由ε可求得杆单元的应力σ及切线模量E^T，可知单元i节点和j节点的轴力N_i和N_j为：

点击浏览公式

设随转坐标系下的杆端力向量为f=[N_i N_j]^T，对式(10)微分，并写成矩阵形式有：

点击浏览公式

设预应力筋杆元在结构坐标系下杆端力向量F=[H_iV_iW_iH_jV_jW_j]^T，由场一致性原则有^{[14, 15]}：

点击浏览公式

微分式(12)可得：

点击浏览公式

为得到δt^Tf用δd来表述的计算式，可先微分式(8)得到：

点击浏览公式

对式(12)中t的微分，并结合式(14)，最终可将δt^Tf用δd表述为：

点击浏览公式

为得到t^Tδf用δd表达的形式，联立式(11)、(9)有：

点击浏览公式

联立式(15)、(16)，可得到结构坐标系下预应力钢筋空间杆元考虑几何与材料双重非线性的切线刚度矩阵K^T为：

点击浏览公式

如图 2所示，要导出预应力钢筋空间杆单元对混合壳单元切线刚度矩阵的贡献，须先确定预应力钢筋空间杆元两端节点9和10的(ξ,η,ζ)坐标。

显然，节点9和10在壳单元局部坐标系下的坐标(X₉,Y₉,Z₉)和(X₁₀,Y₁₀,Z₁₀)、节点9和10距离壳单元底部的距离h₉和h₁₀以及9和10节点处壳单元的厚度t₉和t₁₀是作为已知条件的(并且假定这些值不随壳单元发生位移而改变)。假定预应力空间杆元两端节点9和10分别位于壳单元两端ξ=1和ξ=－1的端面上，故有：

点击浏览公式

点击浏览公式

由于壳单元局部坐标系下任一点坐标(X,Y,Z)都是η的二次函数，在形函数已知的情况下很容易求得6个η₉和6个η₁₀，即：

(1)解X₉=X(ξ₉,η₉,ζ₉)，得η₉¹和η₉²；

(2)解Y₉=Y(ξ₉,η₉,ζ₉)，得η₉³和η₉⁴；

(3)解Z₉=Z(ξ₉,η₉,ζ₉)，得η₉⁵和η₉⁶；

(4)解X₁₀=X(ξ₁₀,η₁₀,ζ₁₀)，得η₁₀¹和η₁₀²；

(5)解Y₁₀=Y(ξ₁₀,η₁₀,ζ₁₀)，得η₁₀³和η₁₀⁴；

(6)解Z₁₀=Z(ξ₁₀,η₁₀,ζ₁₀)，得η₁₀⁵和η₁₀⁶。

由于(X,Y,Z)和(ξ,η,ζ)是一对一变换，所以必有至少3个位于-1≤η≤+1区域内的相同的η₉和η₁₀，这也就是要求的η₉和η₁₀。

至此得到预应力空间杆元两端节点位移向量d^p=[u₉v₉w₉u₁₀v₁₀w₁₀]^T(上标‘p’代表预应力钢筋，以下均同)与壳单元节点位移向量 d^cs=[d₁^csd₂^csd₃^csd₄^csd₅^csd₆^csd₇^csd₈^cs]^T(其中d_j^cs=u_jv_jw_jβ_1jβ_2j,j=1,2,…,8)的关系为：

点击浏览公式

对式(20)微分有：

点击浏览公式

与前述相同，设预应力钢筋空间杆元在结构坐标系下节点力向量为 F^p=[H₈^p V₈^p W₈^p H₉^p V₉^p W₉^p]^T(F^p实际就是1.2.2节的F)，在混合壳单元中预应力钢筋空间杆元在结构坐标系下节点力向量F_cs^p=[^pF₁^{x p}F₁^{y p}F₁^{z p}M₁^{β₁ p}M₁^β₂ … ^pF₈^{x p}F₈^{y p}F₈^{z p}M₈^{β₁ p}M₈^β₂]^T，由虚功原理有：

点击浏览公式

点击浏览公式

点击浏览公式

对式(23)微分且联立式(24)和(21)，可得：

点击浏览公式

因此，由实体等参退化壳元与预应力钢筋空间杆元形成的混合壳单元在结构坐标系下考虑几何与材料双非线性的切线刚度矩阵K为：

点击浏览公式

混凝土及普通钢筋的材料特性见文献^[3]。

预应力钢筋则采用图 4所示的三折线应力-应变关系曲线^[16],图中各参数的定义及具体取值详见文献^[15]。σ_p，ε_p分别为预应力钢筋的应力和应变；f_pu，ε_pu分别为极限抗拉强度和极限拉应变；f_0.2，ε_0.2分别为条件屈服强度和条件屈服应变；f_e，ε_e分别为弹性应力极限和弹性应变极限。

图 4 预应力钢筋的应力-应变曲线 Fig. 4 Stress-strain curve of prestressed bar

图选项

以施加第n级荷载增量为例，介绍本文算法中预应力钢筋空间杆元的非线性分析过程，混凝土及普通钢筋分层壳元的非线性分析过程可见文献^{[12, 13]}。

(1)根据第n-1级荷载增量末壳单元节点在结构坐标系下的总位移向量d^cs，由式(20)求出预应力钢筋空间杆元节点在结构坐标系下的位移向量d^p，结合式(7)求出预应力钢筋应变，由其本构关系求出应力和切线模量，利用式(17)求出结构坐标系下预应力钢筋空间杆元考虑几何与材料双重非线性的切线刚度矩阵K^T，由式(12)求出预应力钢筋空间杆元在结构坐标系下节点力向量F(F^p)；

(2)由式(25)和式(23)分别求出结构坐标系下预应力钢筋空间杆元对混合壳单元切线刚度矩阵与节点力向量的贡献K^p及F_cs^p；

(3)按文献^{[12, 13]}介绍的方法求出结构坐标系下混凝土及普通钢筋分层壳元切线刚度矩阵K^cs与节点力向量F^cs；

(4)由K^p与K^cs之和形成混合壳单元切线刚度矩阵K，由F_cs^p与F^cs之和形成总等效节点力向量F^int；

(5)对所有单元重复(1)至(4)的步骤，得结构切线刚度矩阵∑K和等效节点力向量∑F^int；

(6)计算不平衡力矩阵ΔR=P－∑F^int，其中P为截止到计算t时刻施加的总外荷载所转换成的等效节点力；

(7)解方程∑K·Δd=ΔR，得到节点位移增量Δd，将其叠加到总位移向量d中；

(8)收敛条件判断，如收敛，则转到第n+1级荷载增量计算，如不收敛，则返回(1)，进行下一次迭代计算。

交通部第二公路勘察设计院曾对图 5所示的预应力混凝土简支T梁进行了2片梁的试验，文献^[3]对试验梁进行过非线性有限元分析。该试验T梁总长为24.96 m，计算跨径为24.30 m，梁高为1.45 m，其他几何尺寸、材料参数及加载过程描述详见文献^[3]。

图 5 预应力钢筋混凝土T梁构造图(单位：m) Fig. 5 Structure of prestressed RC T-beam(unit:m)

图选项

采用本文算法对该T梁进行非线性有限元分析，为使计算结果具有可比性，梁的单元网格划分与文献^[3]完全相同。本文及文献^[3]计算得到的T梁极限承载力分别为1 156 kN和1 130 kN，跨中和1/4跨处梁底的荷载-挠度曲线分别见图 6和图 7，可看出两者与试验结果均基本吻合，说明本文算法是正确的。相对于文献^[3]而言，本文推导预应力钢筋切线刚度矩阵的方法力学概念更明确，过程更简单。

图 6 跨中荷载-挠度曲线 Fig. 6 Mid-span load-displacement curves

图选项

图 7 1/4跨荷载-挠度曲线 Fig. 7 One fourth span load-displacement curves

图选项

不同于已有的非线性组合壳单元模型在推导预应力钢筋的切线刚度矩阵时都是基于全拉格朗日列式或增量拉格朗日列式，存在几何与材料两种非线性因素耦合从而导致有限元列式复杂的缺点，本文在推导预应力钢筋的切线刚度矩阵时基于随转坐标法和场一致性原则，利用随转坐标系下刚体位移被消除后得到的位移为纯变形特点，在随转坐标系下考虑材料非线性，几何非线性则通过随转坐标系与结构坐标系之间的转换予以考虑。这种处理不但能实现几何非线性与材料非线性的解耦使有限元列式简单，而且能有效扩大非线性材料本构关系的使用范围。再结合已有的实体退化壳单元，导出适合于预应力混凝土梁双非线性分析的混合壳单元模型(本文方法也可看成已有线性分析方法^[17]向非线性分析领域的延伸)，基于该模型编制了非线性计算程序，算例结果表明该模型及算法是完全正确的，能用于预应力钢筋混凝土梁的非线性分析。