CN104166792A - 一种预应力混凝土连续刚构桥温度作用有限元分析方法 - Google Patents

Abstract

本发明提供一种预应力混凝土连续刚构桥温度作用有限元分析方法，用以分析计算不同温度梯度模型下施工及成桥阶段应力和位移分布，为预应力混凝土连续刚构桥的设计及优化提供参考依据，通过基于有限元和结构力学的计算分析，采用正装计算法按照桥梁结构实际施工加载顺序来进行结构变形和受力分析，使得箱梁温度自应力和预应力混凝土超静定结构中的温度次内力及其次应力这一复杂的力学问题得到了简便、高效的解决，具有重要的实用价值。

Description

一种预应力混凝土连续刚构桥温度作用有限元分析方法

技术领域

本发明涉及一种预应力分析方法，具体是一种预应力混凝土连续刚构桥温度作用有限元分析方法，通过求解有限元方程得到每个节点上的广义节点位移，在每个单元上使用广义节点位移插值求得位移、应力、应变等物理量。

背景技术

近年来交通运输业发展迅速,混凝土的强度越来越高,采用悬臂逐段现浇施工方法修建的预应力混凝土连续刚构桥梁的跨径越来越大,桥墩附近箱梁节段浇筑的体积也越来越大,由于混凝土材料导热性能差,在各种温度变化作用下,预应力混凝土桥梁结构内部会产生相当大的应力、变形、甚至出现温度裂缝,混凝土水化热引起的温度效应非常显著。

现有的设计实践中，由于结构分析的复杂性，现有的计算冗长，虽然设计者主观上希望结构设计尽可能优化，但一方面缺乏有效的计算方法，另一方面，缺乏系统的方法指导桥梁结构设计和改进结构设计，使得结构的设计优化主要依靠累计的工程经验，使得现有的设计过程带有较大的主观性，且工作量大，设计周期长，因而，迫切希望对现有的设计计算方法进行改进优化，在降低计算量的前提下，按照桥梁结构实际施工加载顺序来进行结构变形和受力分析，有效计算箱梁温度自应力和预应力混凝土超静定结构中的温度次内力及其次应力。

发明内容

针对现有技术的缺点和不足，本发明旨在提供一种预应力混凝土连续刚构桥温度作用有限元分析方法，采用正装计算法按照桥梁结构实际施工加载顺序来进行结构变形和受力分析，可有效计算不同温度梯度模型下施工及成桥阶段应力和形变分布，有效解决箱梁温度自应力和预应力混凝土超静定结构中的温度次内力及其次应力计算，为预应力混凝土连续刚构桥的设计及优化提供重要的参考依据。

为实现上述目标，本发明采用正装计算法按照桥梁结构实际施工加载顺序来进行结构变形和受力分析，它能较好地模拟桥梁结构的实际施工过程，能得到桥梁结构在各个施工阶段的位移和受力状态，这不仅可用来指导桥梁设计和施工，而且为桥梁施工控制提供依据。同时，正装计算法能较好地考虑一些与桥梁结构形成历程有关的因素，比如结构的非线性问题、预应力损失问题以及混凝土收缩徐变问题等。正装计算法能根据实际的配筋情况和施工方案设计逐步逐阶段地进行计算，能得到各施工阶段的内力和变形值，最终能得到成桥结构的受力状态。这种计算方法的特点是：随着施工阶段的推进，结构形式、边界约束、荷载形式在不断地改变，前期结构会发生收缩、徐变等，其几何位置也在改变，因而，前一阶段结构状态是本次施工阶段结构分析的基础。

本发明为解决其技术问题所采用的技术方案为：一种预应力混凝土连续刚构桥温度作用有限元分析方法，所述预应力混凝土连续刚构桥包括主梁、纵梁、跨梁和桥墩，所述主梁为预应力混凝土连续刚构结构，所述桥墩与主梁刚性连接，所述有限元分析方法具体包括如下步骤，

SS1.定义有限元模型的长度单位、力的单位和温度单位；

SS2.定义纵梁和桥墩的梁单元类型；

SS3.建立有限元计算物理模型，对全桥进行单元划分，设定边界约束；

SS4.定义混凝土和预应力钢材的材料参数，包括弹性模量、容重、线膨胀系数、剪切模量、泊松比、预应力管道偏差系数、预应力钢筋与管道的摩阻系数、预应力钢筋松弛率、锚具变形及钢束回缩值；

SS5.设置边界条件：设定温度梯度模型及加载，加载包括恒载、预应力荷载、徐变和干缩、挂篮荷载、温度荷载；

SS6.建立偏微分控制方程组并对其进行离散化，得到代数方程组并利用步骤SS5设定的边界条件对所述代数方程组进行封闭，用以计算应力和位移；

SS7.对计算域内的所述代数方程组反复进行迭代计算，直到满足所设定的迭代精度为止，得到应力和位移分布。

进一步地，本发明的预应力混凝土连续刚构桥温度作用有限元分析方法，还包括步骤SS8，对计算结果进行后处理，绘制出相关曲线。

进一步地，步骤SS1中，定义有限元模型的长度单位为m、力的单位为牛(N)和温度单位为摄氏度(℃)。

进一步地，步骤SS3中，包括以下三个子步骤：

a.建立有限元计算物理模型：纵梁和桥墩均采用空间梁单元模拟，整体坐标系以X向为桥纵向，Y向为桥横向，Z向为竖向；

b.对全桥进行单元划分：全桥共划分为m+n个单元，其中主梁划分为m个单元，桥墩划分为n个单元，其中，m、n为自然数；

c.设定边界约束：主梁与桥墩的约束关系通过刚性连接模拟；两边跨梁端只有Y-Z平面内的角位移和水平线位移2个自由度，其余4个方向均被约束；不考虑桩土作用，将墩底直接固结。

进一步地，步骤SS4具体为：

a.主梁采用箱梁结构，定义箱梁及桥墩采用的混凝土的弹性模量、剪切模量、泊松比、轴心抗压强度标准值、轴心抗拉强度标准值、线膨胀系数；

b.定义预应力钢材的材料参数，所述预应力钢材包括纵向预应力钢材和竖向预应力钢材，其中，

--所述纵向预应力钢材采用低松弛钢绞线，定义所述低松弛钢绞线的弹性模量、抗拉强度标准值、张拉控制应力、预应力钢束与管道的摩阻系数、预应力管道偏差系数、预应力钢束松弛率、锚具变形及钢束回缩值；优选地，纵向预应力钢材采用Φj15.24mm低松弛钢绞线，技术标准符合ASTM A416-97(270级)标准要求；

--竖向预应力钢材采用预应力粗钢筋，定义预应力粗钢筋的弹性模量、抗拉强度标准值、张拉控制应力、预应力钢筋与管道的摩阻系数、预应力管道偏差系数、预应力钢筋松弛率、锚具变形及钢束回缩值；优选地，竖向预应力钢材采用32mm的精轧粗钢筋。

进一步地，步骤SS5按照如下方式设置温度梯度模型及加载：

a.设置温度梯度模型

选择不同的温度梯度模型，以计算预应力混凝土连续刚构桥在施工过程中，最大悬臂阶段在不同施工阶段时温度应力及应变的变化情况，从以下六种温度梯度模型进行选择，

1)温度梯度模型1：桥面板均匀升温，温差为5℃，箱梁腹板和底板无温差；

2)温度梯度模型2：桥面板表面的最高温度取20℃；

3)温度梯度模型3：按照美国AASHTO规范对温度梯度的规定建立的温度梯度模型；

4)温度梯度模型4：按照英国BS5400规范升温时的温度梯度建立的模型；

5)温度梯度模型5：温度梯度是一条高1200mm的五次抛物线，混凝上表面的温度取32℃，在截面厚度为200mm的底板上采用从0℃到1.5℃的线性温度增长；

6)温度梯度模型6：根据实测结果提出的温度梯度模型，温度梯度曲线为T_y＝T₀e^-αy，T₀为顶板温度，α按最小二乘法进行非线性拟合，y为高度；优选地，根据边跨、中跨1/4截面实测温度梯度值，按最小二乘法进行非线性拟合得到α＝1.57，高度y以米计，顶板温度T₀为20℃；

b.设置的荷载包括恒载、预应力荷载、徐变和干缩、挂篮荷载、温度荷载。

进一步地，步骤SS6中，按照如下方式建立偏微分控制方程组并对其进行离散化：

(1)梁单元刚度矩阵

在局部坐标系下，非线性有限元的梁单元刚度矩阵可以表示为：

[K_T]_B＝[K_E]_B+[K_G]_B      (1)

其中：

[K_T]_B表示局部坐标系下的梁单元切线刚度矩阵；

[K_E]_B表示局部坐标系下的梁单元弹性刚度矩阵；

[K_G]_B表示局部坐标系下的梁单元几何刚度矩阵。

对于梁单元弹性刚度矩阵来说，有：

其中：EI_Y、EI_Z为沿Y轴、Z轴的抗弯刚度；GI_X为沿X轴的抗扭刚度；

应用虚功原理可以推导出非线性有限元平衡方程如下：

[K]{U}＝{P}      (4)

其中，[K]表示结构的整体刚度；

{U}表示全部自由度的位移向量；

{P}表示荷载向量；

(2)温度荷载

混凝土结构内部的温度场是确定温度荷载的关键。

箱梁结构的主梁体内任意一点的温度T是坐标x，y，z和时间t的函数，设定混凝土为均质、各向同性、无内热源，得三维非稳态导热方程：

$\mathrm{\λ}(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{{\∂y}^2}+\frac{{\∂}^{2}T}{{\∂z}^2})=\mathrm{c\γ}\frac{\∂T}{\∂t}---\left(5\right)$

式中：λ——混凝土的导热系数；

c——混凝土的比热；

γ——混凝土的容重；

略去桥长方向温差的影响；

在梁高较小时，略去水平方向很小的热传导作用，用垂直方向的一维热传导状态来分析；

对于梁高较大的箱梁，忽略角隅区附近的热传导状态，用垂直和水平两个方向各自的一维导热状态分别计算，然后再叠加起来；

在箱梁结构计算中简化为一维热传导方程：

$\mathrm{\α}\frac{{\∂}^{2}T}{{\∂x}^2}=\frac{\∂T}{\∂t}---\left(6\right)$

第一类边界条件：混凝土表面温度T是时间的已知函数，

当t＝0时，T(t)＝f(t)      (7)

对于壁板结构，将其近似为一块半无限厚板，并假定气温变化为谐波形式的情况下，根据式(6)可得上述边界条件的弹性力学解：

$T(x,t)=A_0{e}^{-\sqrt{\mathrm{\ω}/2\mathrm{\αx}}}\sin (\mathrm{\ωt}-\sqrt{\mathrm{\ω}/2\mathrm{\αx}})---\left(8\right)$

式中：A₀——壁板表面温度波动峰值；

α——热扩散系数；

ω——圆频率2π/24；

x——计算点至板表面距离；

t——时间；

以某一特定时刻最大温差分布相应的温差荷载作为控制荷载，式(8)表示为温度分布包络线的形式：

$T\left(x\right)=A_0{e}^{-\sqrt{\mathrm{\ω}/2\mathrm{\αx}}}---\left(9\right)$

用函数式表示沿板厚的温度分布，C_x为实验参数；对高架桥箱梁模型采用来计算沿板厚的温度分布；

对混凝土箱形桥墩，采用下式来分析壁厚方向的温度分布：

$T\left(x\right)=T_0{e}^{-C_{X}X}---\left(10\right)$

其中T₀为墩壁内外表面的温差；沿箱梁高、梁宽方向的温差分布按下式计算：

$T\left(y\right)=T_{0y}{e}^{-C_{Y}Y},T\left(x\right)=T_{0X}{e}^{-C_{X}X}---\left(11\right)$

式中：T_0y、T_0X——沿梁高、梁宽方向的温差；

y、X——计算点至受热表面的距离；

C_X、C_y——指数系数，随结构形式、部位、计算时刻而异；

以上式(9)至式(11)，都是采用第一类边界条件、考虑温度分布包络线的函数式，影响第一类边界条件的外部主要因素为太阳辐射强度、气温变化、风速；

混凝土内部的温度分布确定后，根据混凝土的热物理性能，利用线膨胀系数，形成温度荷载，考虑桥面受日照后形成的沿箱梁高度变化的温度梯度T_y＝T₀e^-αy；

(3)温差应力

确定温度梯度模型及温度设计值后，温度应力按结构力学和有限元方法进行计算，计算时假定：

1)桥长方向的温度分布是均匀的；

2)混凝土是弹性匀质材料；

3)梁变形服从平面假定；

4)按单向温差荷载计算温差应力，然后叠加组合多向温差荷载状下的温差应力；

温度应力有由两部分组成：a)梁的温度变形受到纵向纤维之间的相互约束，在截面上产生自平衡的纵向约束应力，即自应力；b)梁的温度上拱变化受到支承条件约束的温度次应力；

①箱梁温度自应力

设温度梯度沿梁高按任意曲线t(y)分布，取单位梁长dy＝1的微分段，当纵向纤维之间不受约束且自由伸缩时，沿梁高的自由应变ε_t(y)与温度梯度一致，即：

ε_t(y)＝α_ct_(y)      (12)

由于纵向纤维之间的相互约束，梁截面应变应符合平面假定，梁截面上的最终应变ε_f(y)为直线分布，即：

ε_f(y)＝ε₀+ψy      (13)

式中ε₀——基轴y＝0处应变；

ψ——截面变形曲率；

y——基轴以下任一点求应变的坐标；

α_c——混凝土线膨胀系数；

自由应变与最终应变之差，系纤维之间的约束产生，其值为：

ε_σ(y)＝ε_t(y)-ε_f(y)＝α_ct_(y)-(ε₀+ψy)      (14)

自应力为：

σ_s(y)＝E_cε_σ(y)＝E_c[α_ct_(y)-(ε₀+ψy)]      (15)

全截面上轴力N和弯矩M

$\begin{array}{c}N=E_c\underset{h}{\∫}{\mathrm{\ϵ}}_{\mathrm{\σ}\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}=E_c\underset{h}{\∫}({\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy})b_{\left(y\right)}\mathrm{dy}\\ =E_c[{\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}-\mathrm{\ψ}\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}]\end{array}---\left(16\right)$

$\begin{array}{c}M=E_c\underset{h}{\∫}{\mathrm{\ϵ}}_{\mathrm{\σ}\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=E_c\underset{h}{\∫}({\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy})b_{\left(y\right)}(y-y_c)\mathrm{dy}\\ =E_c[{\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}-{\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}-\mathrm{\ψ}\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}]\end{array}---\left(17\right)$

式中E_c——混凝土材料弹性模量；

b_(y)——y处的梁宽；

对于任何截面，N＝0，M＝0，即内力总和为零；

公式(16)、(17)分别改写为：

${\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}+\mathrm{\ψ}\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}---\left(18\right)$

${\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}+\mathrm{\ψ}\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}---\left(19\right)$

在公式(18)、(19)内

$\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}=A---\left(20\right)$

$\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}={\mathrm{Ay}}_c---\left(21\right)$

$\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}=\underset{h}{\∫}{b}_{\left(y\right)}{y}^2\mathrm{dy}-\underset{h}{\∫}{b}_{\left(y\right)}{\mathrm{yy}}_c\mathrm{dy}=I_b-\underset{h}{\∫}{b}_{\left(y\right)}{\mathrm{yy}}_c\mathrm{dy}=I_g---\left(22\right)$

$\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}=0$ (对重心轴的静面积矩为零)

式中A——截面面积；

I_b——截面面积对基轴惯性矩；

I_g——截面面积对重心轴惯性矩；

将公式(20)～(22)代入公式(18)、(19)内，得：

${\mathrm{\ϵ}}_{0}A+\mathrm{\ψ}{\mathrm{Ay}}_c={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}---\left(23\right)$

ψI_g＝α_c∫t_(y)b_(y)(y-y_c)dy      (24)

由公式(23)、(24)可得：

${\mathrm{\ϵ}}_0=\frac{{\mathrm{\α}}_c}{A}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c---\left(25\right)$

$\mathrm{\ψ}=\frac{{\mathrm{\α}}_c}{I_g}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}---\left(26\right)$

设在坐标y处，截面内一厚度为i的微小单元面积A_y处温度梯度值为t_y，以t_y为常值代入公式(25)、(26)，积分区段仅在i厚度范围内有值：t_(y)＝t_y,y-y_c＝e_y(单元面积A_y对全面积重心的偏心距)。

${\mathrm{\ϵ}}_0=\frac{{\mathrm{\α}}_c}{A}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c=\frac{{\mathrm{\α}}_c}{A}\underset{i}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c=\frac{{\mathrm{\α}}_c{t}_y{A}_y}{A}-\frac{{\mathrm{\α}}_c{t}_y{A}_y{e}_y{y}_c}{I_g}---\left(27\right)$

$\mathrm{\ψ}=\frac{{\mathrm{\α}}_c}{I_g}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=\frac{{\mathrm{\α}}_c}{I_g}\underset{i}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=\frac{{\mathrm{\α}}_c{t}_y{A}_y{e}_y}{I_g}---\left(28\right)$

自公式(15)可求得任意点应力σ_s(y)：

${\mathrm{\σ}}_{s\left(y\right)}=E_c[{\mathrm{\α}}_c{t}_{\left(y\right)}-({\mathrm{\ϵ}}_0+\mathrm{\ψy})]=E_c{\mathrm{\α}}_c{t}_y-\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y}{A}+\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y{e}_y{y}_c}{I_g}-\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y{e}_{y}y}{I_g}---\left(29\right)$

令：N_ti＝A_yt_yα_cE_c,M_ti＝-N_tie_y＝-A_yt_yα_cE_ce_y

${\mathrm{\σ}}_{s\left(y\right)}=-\frac{N_{\mathrm{ti}}}{A}+\frac{M_{\mathrm{ti}}}{I_g}(y-y_c){+t}_y{\mathrm{\α}}_c{E}_c---\left(30\right)$

公式(30)是由于一个单元面积A_y内的温度作用，在截面任一点产生的应力；对于分为很多块单元面积上不同t_y的作用，应用分段总和法；公式(30)使用于正温差；如为反温差则整个公式前冠以负号；

②超静定结构中的温度次内力及其次应力

在预应力混凝土超静定结构中，前述温度变形ε₀及曲率ψ将受到超静定赘余约束的制约，引起温度次内力，两端固定杆单元的节点荷载向量{F}^e由截面变形曲率及沿梁高y＝0处的变形ε₀直接写出：

${\left\{F\right\}}^e=\left\{\begin{array}{c}{N}_i\\ Q_i\\ M_i\\ N_j\\ Q_j\\ M_j\end{array}\right\}=\left\{\begin{array}{c}\mathrm{EA}({\mathrm{\ϵ}}_0+{\mathrm{\ψy}}_c)\\ 0\\ \mathrm{EI\ψ}\\ -\mathrm{EA}({\mathrm{\ϵ}}_0+{\mathrm{\ψy}}_c)\\ 0\\ -\mathrm{EI\ψ}\end{array}\right\}---\left(31\right)$

杆件单元节点力应以结构坐标系表示，然后分别组集各个杆件单元的结点荷载，从而得到节点外力向量{F}，矩阵位移方程为：

[K]{Δ}+{F}＝0      (32)

式中[K]——结构总刚度矩阵；

{Δ}——单元节点位移向量；

在求得结构各单元因温度变化引起的结点位移后，由单元的杆端力与单元刚度矩阵、单元结点位移的关系{f}^e＝[K]{Δ}^e求得结构的温度次内力N_T、Q_T、M_T及其次应力；在超静定结构中，总的温度应力为：

纵向弯曲应力：

${\mathrm{\σ}}_{t\left(y\right)}=\frac{N_T}{A}+\frac{M_T}{I}y+E_c[{\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy}]---\left(33\right).$

同现有技术相比，本发明具有显著的技术效果：

(1)本发明方法通过基于有限元和结构力学的计算程序，采用正装计算法按照桥梁结构实际施工加载顺序来进行结构变形和受力分析，使得箱梁温度自应力和预应力混凝土超静定结构中的温度次内力及其次应力这一复杂的力学问题得到了简便、高效的解决，具有重要的实用价值。

(2)本发明采用正装计算法按照桥梁结构实际施工加载顺序来进行结构变形和受力分析，它能较好地模拟桥梁结构的实际施工过程，能得到桥梁结构在各个施工阶段的位移和受力状态，这不仅可用来指导桥梁设计和施工，而且为桥梁施工控制提供依据。同时，正装计算法能较好地考虑一些与桥梁结构形成历程有关的因素，比如结构的非线性问题、预应力损失问题以及混凝土收缩徐变问题等。正装计算法能根据实际的配筋情况和施工方案设计逐步逐阶段地进行计算，能得到各施工阶段的内力和变形值，最终能得到成桥结构的受力状态。这种计算方法的特点是：随着施工阶段的推进，结构形式、边界约束、荷载形式在不断地改变，前期结构会发生收缩、徐变等，其几何位置也在改变，因而，前一阶段结构状态是本次施工阶段结构分析的基础。

附图说明

图1为温度梯度模型1的温度梯度曲线；

图2为温度梯度模型2的温度梯度曲线；

图3为温度梯度模型3的温度梯度曲线；

图4为温度梯度模型4的温度梯度曲线；

图5为温度梯度模型5的温度梯度曲线；

图6为温度梯度模型6的温度梯度曲线；

图7为箱梁温度计算模型，其中，图7(a)为箱梁截面示意图；图7(b)为温度梯度沿梁高按任意曲线t(y)分布；图7(c)为箱梁平面变形示意图；图7(d)为箱梁沿梁高自应力应变示意图。

图8为预应力混凝土连续刚构桥正装计算法计算程序流程图；

图9为在最大悬臂施工阶段时六种温度梯度模型下主梁上缘、下缘应力分布图，其中，图9(a)、(b)为温度梯度模型1下主梁上缘、下缘应力分布图，图9(c)、(d)为温度梯度模型2下主梁上缘、下缘应力分布图，图9(e)、(f)为温度梯度模型3下主梁上缘、下缘应力分布图，图9(g)、(h)为温度梯度模型4下主梁上缘、下缘应力分布图，图9(i)、(j)为温度梯度模型5下主梁上缘、下缘应力分布图，图9(k)、(l)为温度梯度模型6下主梁上缘、下缘应力分布图；

图10为最大悬臂施工阶段时六种温度梯度模型下主梁上下缘应力分布对比示意图；

图11为最大悬臂施工阶段时六种温度梯度模型引起的主梁1/4截面的应力沿梁高的分布示意图；

图12为不同温度梯度模型下主梁位移(主梁变形)示意图；

图13为在成桥阶段时六种温度梯度模型下主梁上缘、下缘应力分布图，其中，图13(a)、(b)为温度梯度模型1下主梁上缘、下缘应力分布图，图13(c)、(d)为温度梯度模型2下主梁上缘、下缘应力分布图，图13(e)、(f)为温度梯度模型3下主梁上缘、下缘应力分布图，图13(g)、(h)为温度梯度模型4下主梁上缘、下缘应力分布图，图13(i)、(j)为温度梯度模型5下主梁上缘、下缘应力分布图，图13(k)、(l)为温度梯度模型6下主梁上缘、下缘应力分布图；

图14为在成桥阶段时六种温度梯度模型下主梁上下缘应力分布对比示意图；

图15为在成桥阶段时六种温度梯度模型下主梁各截面最大应力分布对比示意图；

图16为在成桥阶段时六种温度梯度模型下主梁各截面竖向位移分布对比示意图；

图17为系统升、降温20℃引起的主梁竖向位移分布示意图；

图18为三种工况下主梁竖向位移分布对比示意图，其中，工况1为自重、二期恒载、收缩徐变、预应力及施工荷载的共同作用；工况2为工况1+体系升温；工况3为工况1+体系降温；

图19为四种工况下主梁竖向位移分布对比示意图，其中，工况1为合计(自重、二期恒载、收缩徐变、预应力及施工荷载的共同作用)+体系升温20℃；工况2为合计(自重、二期恒载、收缩徐变、预应力及施工荷载的共同作用)+体系降温20℃；工况3为工况1+温度梯度模型6；工况4为工况2+温度梯度模型6；

图20为体系降温20℃时主梁上下缘应力分布对比示意图；

图21为体系降温20℃时主梁上下缘应力分布对比示意图；

图22为有限元计算物理模型；

图23为悬臂浇筑桥梁结构施工阶段定义；

图24为徐变系数计算结果；

图25为混凝土干缩应变计算结果；

图26为挂篮结构简图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下参照附图并举实施例，对本发明进一步详细说明。需要说明的是，以下所述仅为本发明的较佳实施例，并不因此而限定本发明的保护范围。

本发明涉及的计算物理模型系工程实例中的水磨湾特大桥，该桥位于少林寺至洛阳高速公路K21+910m处，主桥为预应力混凝土连续刚构桥，跨经组成为65+110+65＝240m，空心薄壁式桥墩、钻孔灌注桩基础。上部结构为单箱单室断面，顶板宽度为12.75m，底板宽度为6.5m，箱梁根部梁高6.0m，跨中及边跨合拢段梁高为2.3m，箱梁底板下缘按二次抛物线变化。主桥主梁采用三向预应力体系：纵向预应力氛围顶板束、底板束边跨合拢钢束三种，采用ASTM-92标准的15Φ^j 15.24mm的270级钢绞线，OVM15-15锚具，张拉力为2930kN；竖向预应力采用Φ^j 32mm的高强精轧螺纹粗钢筋，YGM锚具，张拉力为452kN；横向预应力采用ASTM-92标准的3Φ^j 15.24毫米的270级钢绞线，BM15-3锚具，张拉力为586kN。

该工程物理模型数据详尽且具有代表性，本发明对该工程物理模型进行数值计算，采用空间有限元分析程序，其主要实现步骤如下：

SS1.定义有限元模型的长度单位为m、力的单位为N和温度单位为℃，其它单位均由以上单位换算得出；

SS2.主梁采用三维变截面梁单元，桥墩采用三维等截面梁单元，底部和顶部局部采用变截面梁单元，预应力钢束直接模拟成预应力荷载加入模型中；

SS3.建立有限元计算物理模型：如图22所示，纵梁和桥墩均采用空间梁单元模拟，整体坐标系以X向为桥纵向，Y向为桥横向，Z向为竖向。对全桥进行单元划分：全桥共划分为78个单元，其中主梁划分为66个单元，主墩划分为12个单元。设定边界约束：悬臂浇筑桥梁结构施工阶段定义与实际施工阶段划分一致，如图23所示；两个边主梁支座模拟为竖向位移受到约束的铰支座，中间两个刚构墩和主梁之间采用弹性连接中的刚接。

4.定义混凝土和钢材的材料参数包括弹性模量、容重、线膨胀系数、剪切模量、泊松比、预应力管道偏差系数、预应力管道偏差系数、预应力钢筋与管道的摩阻系数、预应力钢筋松弛率、锚具变形及钢束回缩值等；

(1)混凝土

箱梁及主墩采用C50混凝土：

☆弹性模量：            E＝3.45×10⁴MPa

☆剪切模量：            G＝1.38×10⁴MPa

☆泊松比：              v＝0.2

☆轴心抗压强度标准值：  f_ck＝32.4MPa

☆轴心抗拉强度标准值：  f_tk＝2.65MPa

☆线膨胀系数：          k＝1.0×10^-5(1/℃)

(2)预应力钢材

纵向预应力采用Φj15.24mm低松弛钢绞线，技术标准符合ASTM A416-97(270级)标准要求。竖向预应力采用32mm的精轧粗钢筋。

①钢绞线

☆弹性模量：             E＝1.95×10⁵MPa

☆抗拉强度标准值：          f_pk＝1860MPa

☆张拉控制应力：                 1395MPa

☆预应力钢束与管道的摩阻系数：   0.25

☆预应力管道偏差系数：        0.0015/m

☆预应力钢束松弛率：          0.035

☆锚具变形及钢束回缩值：      0.006m

②预应力粗钢筋

☆弹性模量：            E＝2.0×10⁵MPa

☆抗拉强度标准值：         f_pk＝750MPa

☆张拉控制应力：                675MPa

☆预应力钢筋与管道的摩阻系数：  0.4

☆预应力管道偏差系数：          0.0015

☆预应力钢筋松弛率：            0.03

☆锚具变形及钢束回缩值：        0.002m

SS5.设置温度梯度模型及加载，加载包括恒载、预应力荷载、徐变和干缩、挂篮荷载、温度荷载等：

(1)温度梯度模型

为了分析不同温度梯度模型作用下梁内温度应力和应变的变化情况，本发明选择了六种不同的温度梯度模型，来计算预应力混凝土连续刚构桥在施工过程中，最大悬臂阶段不同施工阶段时温度应力及应变的变化情况。

1)温度梯度模型1：桥面板均匀升温，温差为5℃，箱梁腹板和底板无温差，即《公路钢筋混凝土及预应力钢筋混凝土桥梁设计规范》(JTJ023-85)中规定的温度梯度模型。温度梯度曲线如图1所示。

2)温度梯度模型2：桥面板表面的最高温度取20℃，该温度梯度模型与我国《公路桥涵通用设计规范》(JTG D60-2004)中对温度梯度的规定一致。具体温度梯度曲线见图2。

3)温度梯度模型3：按照美国AASHTO规范对温度梯度的规定建立的温度梯度模型，该模型与我国《公路桥涵通用设计规范》(JTG D60-2004)中对温度梯度的规定的差别在于，在截面厚度为200mm的底板上采用从0℃到2.5℃的线性温度增长。具体温度梯度曲线见图3。

4)梯度模型4：按照英国BS5400规范升温时的温度梯度建立的模型。具体温度梯度曲线见图4。

5)梯度模型5：温度梯度是一条高1200mm的五次抛物线，混凝上表面的温度取32℃。在截面厚度为200mm的底板上采用从0℃到1.5℃的线性温度增长。梯度模型与新西兰规范相近。具体温度梯度曲线见图5。

6)梯度模型6：根据实测结果提出的温度梯度模型。该模型与我国铁路规范的温度梯度模型相似^[6]，温度梯度曲线为T_y＝T₀e^-αy，根据边跨、中跨1/4截面实测温度梯度值，按最小二乘法进行非线性拟合得到α＝1.57，高度y以米计，顶板温度T₀为20℃。具体温度梯度曲线见图6。

(2)设置荷载，包括恒载、预应力荷载、徐变和干缩、挂篮荷载、温度荷载等：

1)恒载

一期恒载：重量按设计尺寸计算，混凝土容重按26kN/m³计。

二期恒载按50kN/m计。

2)预应力荷载

钢束

截面面积:A_p＝1.387×15＝2080.5mm²

孔道直径:100/103mm

张拉控制力:施加75％抗拉强度的张力：

f_pj＝0.72f_pu＝1395MPa

P_j＝A_p·f_pj＝2930kN

张拉初期的损失：

摩擦损失:P_(X)＝P₀·e^-(μα+kL)

顶板束:μ＝0.20,k＝0.001

底板束:μ＝0.20,k＝0.006

锚固端滑移量:ΔI_c＝6mm

混凝土弹性压缩预应力损失:损失量,ΔP_E＝Δf_P·A_SP

预应力长期损失：

应力松弛

徐变和干缩引起的损失

3)徐变和干缩

水泥:普通水泥

施加持续荷载时混凝土的材龄:t₀＝5日

混凝土暴露在大气中时的材龄:t_s＝3日

相对湿度:RH＝55％

大气或养生温度:T＝20℃

适用标准:道桥设计标准(CEB-FIP)

徐变系数:由程序计算，结果如图24

混凝土干缩应变:由程序计算，结果如图25

4)挂篮荷载

挂篮结构简图如图26所示，挂篮自重如下：

P＝650kN

e＝2.50m

M＝P×e＝1625kN

5)温度荷载

主桥箱梁合拢温度取10～15℃，升温按20℃计，降温按20℃计。

SS6.建立偏微分控制方程组并对其进行离散化：

(1)梁单元刚度矩阵

在局部坐标系下，非线性有限元的梁单元刚度矩阵可以表示为：

[K_T]_B＝[K_E]_B+[K_G]_B      (1)

其中：

[K_T]_B表示局部坐标系下的梁单元切线刚度矩阵；

[K_E]_B表示局部坐标系下的梁单元弹性刚度矩阵；

[K_G]_B表示局部坐标系下的梁单元几何刚度矩阵。

对于梁单元弹性刚度矩阵来说，有：

其中：EI_Y、EI_Z为沿Y轴、Z轴的抗弯刚度；GI_X为沿X轴的抗扭刚度；

应用虚功原理可以推导出非线性有限元平衡方程如下：

[K]{U}＝{P}      (4)

其中，[K]表示结构的整体刚度；

{U}表示全部自由度的位移向量；

{P}表示荷载向量；

(2)温度荷载

混凝土结构内部的温度场是确定温度荷载的关键。

箱梁结构的主梁体内任意一点的温度T是坐标x，y，z和时间t的函数，设定混凝土为均质、各向同性、无内热源，得三维非稳态导热方程：

$\mathrm{\λ}(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{{\∂y}^2}+\frac{{\∂}^{2}T}{{\∂z}^2})=\mathrm{c\γ}\frac{\∂T}{\∂t}---\left(5\right)$

式中：λ——混凝土的导热系数；

c——混凝土的比热；

γ——混凝土的容重；

略去桥长方向温差的影响；

在梁高较小时，略去水平方向很小的热传导作用，用垂直方向的一维热传导状态来分析；

对于梁高较大的箱梁，忽略角隅区附近的热传导状态，用垂直和水平两个方向各自的一维导热状态分别计算，然后再叠加起来；

在箱梁结构计算中简化为一维热传导方程：

$\mathrm{\α}\frac{{\∂}^{2}T}{{\∂x}^2}=\frac{\∂T}{\∂t}---\left(6\right)$

第一类边界条件：混凝土表面温度T是时间的已知函数，

当t＝0时，T(t)＝f(t)      (7)

对于壁板结构，将其近似为一块半无限厚板，并假定气温变化为谐波形式的情况下，根据式(6)可得上述边界条件的弹性力学解：

$T(x,t)=A_0{e}^{-\sqrt{\mathrm{\ω}/2\mathrm{\αx}}}\sin (\mathrm{\ωt}-\sqrt{\mathrm{\ω}/2\mathrm{\αx}})---\left(8\right)$

式中：A₀——壁板表面温度波动峰值；

α——热扩散系数；

ω——圆频率2π/24；

x——计算点至板表面距离；

t——时间；

以某一特定时刻最大温差分布相应的温差荷载作为控制荷载，式(8)表示为温度分布包络线的形式：

$T\left(x\right)=A_0{e}^{-\sqrt{\mathrm{\ω}/2\mathrm{\αx}}}---\left(9\right)$

用函数式表示沿板厚的温度分布，C_x为实验参数；对高架桥箱梁模型采用来计算沿板厚的温度分布；

对混凝土箱形桥墩，采用下式来分析壁厚方向的温度分布：

$T\left(x\right)=T_0{e}^{-C_{X}X}---\left(10\right)$

其中T₀为墩壁内外表面的温差；沿箱梁高、梁宽方向的温差分布按下式计算：

$T\left(y\right)=T_{0y}{e}^{-C_{Y}Y},T\left(x\right)=T_{0X}{e}^{-C_{X}X}---\left(11\right)$

式中：T_0y、T_0X——沿梁高、梁宽方向的温差；

y、X——计算点至受热表面的距离；

C_X、C_y——指数系数，随结构形式、部位、计算时刻而异；

以上式(9)至式(11)，都是采用第一类边界条件、考虑温度分布包络线的函数式，影响第一类边界条件的外部主要因素为太阳辐射强度、气温变化、风速；

混凝土内部的温度分布确定后，根据混凝土的热物理性能，利用线膨胀系数，形成温度荷载，考虑桥面受日照后形成的沿箱梁高度变化的温度梯度T_y＝T₀e^-αy；

(3)温差应力

确定温度梯度模型及温度设计值后，温度应力按结构力学和有限元方法进行计算，计算时假定：

1)桥长方向的温度分布是均匀的；

2)混凝土是弹性匀质材料；

3)梁变形服从平面假定；

4)按单向温差荷载计算温差应力，然后叠加组合多向温差荷载状下的温差应力；

温度应力有由两部分组成：a)梁的温度变形受到纵向纤维之间的相互约束，在截面上产生自平衡的纵向约束应力，即自应力；b)梁的温度上拱变化受到支承条件约束的温度次应力；

①箱梁温度自应力

图7为箱梁温度计算模型，其中，图7(a)为箱梁截面示意图；图7(b)为温度梯度沿梁高按任意曲线t(y)分布；图7(c)为箱梁平面变形示意图；图7(d)为箱梁沿梁高自应力应变示意图。设温度梯度沿梁高按任意曲线t(y)分布，取单位梁长dy＝1的微分段，当纵向纤维之间不受约束且自由伸缩时，沿梁高的自由自由应变ε_t(y)与温度梯度一致，即：

ε_t(y)＝α_ct_(y)      (12)

由于纵向纤维之间的相互约束，梁截面应变应符合平面假定，梁截面上的最终应变ε_f(y)应为直线分布，即：

ε_f(y)＝ε₀+ψy      (13)

式中ε₀——基轴y＝0处应变；

ψ——截面变形曲率；

y——基轴以下任一点求应变的坐标；

α_c——混凝土线膨胀系数。

自由应变与最终应变之差，系纤维之间的约束产生，其值为：

ε_σ(y)＝ε_t(y)-ε_f(y)＝α_ct_(y)-(ε₀+ψy)      (14)

自应力为：

σ_s(y)＝E_cε_σ(y)＝E_c[α_ct_(y)-(ε₀+ψy)]      (15)

全截面上轴力N和弯矩M

$\begin{array}{c}N=E_c\underset{h}{\∫}{\mathrm{\ϵ}}_{\mathrm{\σ}\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}=E_c\underset{h}{\∫}({\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy})b_{\left(y\right)}\mathrm{dy}\\ =E_c[{\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}-\mathrm{\ψ}\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}]\end{array}---\left(16\right)$

$\begin{array}{c}M=E_c\underset{h}{\∫}{\mathrm{\ϵ}}_{\mathrm{\σ}\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=E_c\underset{h}{\∫}({\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy})b_{\left(y\right)}(y-y_c)\mathrm{dy}\\ =E_c[{\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}-{\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}-\mathrm{\ψ}\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}]\end{array}---\left(17\right)$

式中E_c——混凝土材料弹性模量；

b_(y)——y处的梁宽；

对于任何截面，N＝0，M＝0，即内力总和为零；

公式(16)、(17)分别改写为：

${\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}+\mathrm{\ψ}\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}---\left(18\right)$

${\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}+\mathrm{\ψ}\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}---\left(19\right)$

在公式(18)、(19)内

$\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}=A---\left(20\right)$

$\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}={\mathrm{Ay}}_c---\left(21\right)$

$\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}=\underset{h}{\∫}{b}_{\left(y\right)}{y}^2\mathrm{dy}-\underset{h}{\∫}{b}_{\left(y\right)}{\mathrm{yy}}_c\mathrm{dy}=I_b-\underset{h}{\∫}{b}_{\left(y\right)}{\mathrm{yy}}_c\mathrm{dy}=I_g---\left(22\right)$

$\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}=0$   (对重心轴的静面积矩为零)

式中A——截面面积；

I_b——截面面积对基轴惯性矩，参见图7；

I_g——截面面积对重心轴惯性矩，参见图7；

将公式(20)～(22)代入公式(18)、(19)内，得：

${\mathrm{\ϵ}}_{0}A+\mathrm{\ψ}{\mathrm{Ay}}_c={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}---\left(23\right)$

ψI_g＝α_c∫t_(y)b_(y)(y-y_c)dy      (24)

由公式(23)、(24)可得：

${\mathrm{\ϵ}}_0=\frac{{\mathrm{\α}}_c}{A}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c---\left(25\right)$

$\mathrm{\ψ}=\frac{{\mathrm{\α}}_c}{I_g}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}---\left(26\right)$

设在坐标y处，截面内一厚度为i的微小单元面积A_y处温度梯度值为t_y，以t_y为常值代入公式(25)、(26)，积分区段仅在i厚度范围内有值：t_(y)＝t_y,y-y_c＝e_y(单元面积A_y对全面积重心的偏心距)。

${\mathrm{\ϵ}}_0=\frac{{\mathrm{\α}}_c}{A}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c=\frac{{\mathrm{\α}}_c}{A}\underset{i}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c=\frac{{\mathrm{\α}}_c{t}_y{A}_y}{A}-\frac{{\mathrm{\α}}_c{t}_y{A}_y{e}_y{y}_c}{I_g}---\left(27\right)$

$\mathrm{\ψ}=\frac{{\mathrm{\α}}_c}{I_g}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=\frac{{\mathrm{\α}}_c}{I_g}\underset{i}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=\frac{{\mathrm{\α}}_c{t}_y{A}_y{e}_y}{I_g}---\left(28\right)$

自公式(15)可求得任意点应力σ_s(y)：

${\mathrm{\σ}}_{s\left(y\right)}=E_c[{\mathrm{\α}}_c{t}_{\left(y\right)}-({\mathrm{\ϵ}}_0+\mathrm{\ψy})]=E_c{\mathrm{\α}}_c{t}_y-\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y}{A}+\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y{e}_y{y}_c}{I_g}-\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y{e}_{y}y}{I_g}---\left(29\right)$

令：N_ti＝A_yt_yα_cE_c,M_ti＝-N_tie_y＝-A_yt_yα_cE_ce_y

${\mathrm{\σ}}_{s\left(y\right)}=-\frac{N_{\mathrm{ti}}}{A}+\frac{M_{\mathrm{ti}}}{I_g}(y-y_c){+t}_y{\mathrm{\α}}_c{E}_c---\left(30\right)$

公式(30)是由于一个单元面积A_y内的温度作用，在截面任一点产生的应力；对于分为很多块单元面积上不同t_y的作用，应用分段总和法；公式(30)使用于正温差；如为反温差则整个公式前冠以负号；

②超静定结构中的温度次内力及其次应力

在预应力混凝土超静定结构中，前述温度变形ε₀及曲率ψ将受到超静定赘余约束的制约，引起温度次内力，两端固定杆单元的节点荷载向量{F}^e由截面变形曲率及沿梁高y＝0处的变形ε₀直接写出：

${\left\{F\right\}}^e=\left\{\begin{array}{c}{N}_i\\ Q_i\\ M_i\\ N_j\\ Q_j\\ M_j\end{array}\right\}=\left\{\begin{array}{c}\mathrm{EA}({\mathrm{\ϵ}}_0+{\mathrm{\ψy}}_c)\\ 0\\ \mathrm{EI\ψ}\\ -\mathrm{EA}({\mathrm{\ϵ}}_0+{\mathrm{\ψy}}_c)\\ 0\\ -\mathrm{EI\ψ}\end{array}\right\}---\left(31\right)$

杆件单元节点力应以结构坐标系表示，然后分别组集各个杆件单元的结点荷载，从而得到节点外力向量{F}，矩阵位移方程为：

[K]{Δ}+{F}＝0      (32)

式中[K]——结构总刚度矩阵；

{Δ}——单元节点位移向量；

在求得结构各单元因温度变化引起的结点位移后，由单元的杆端力与单元刚度矩阵、单元结点位移的关系{f}^e＝[K]{Δ}^e求得结构的温度次内力N_T、Q_T、M_T及其次应力；在超静定结构中，总的温度应力为：

纵向弯曲应力：

${\mathrm{\σ}}_{t\left(y\right)}=\frac{N_T}{A}+\frac{M_T}{I}y+E_c[{\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy}]---\left(33\right).$

SS7.对计算域内的所述代数方程组反复进行迭代计算，直到满足所设定的迭代精度为止，得到温度应力和位移分布。图8为预应力混凝土连续刚构桥正装计算法计算程序流程图。

SS8.对计算结果进行后处理，绘制出相关曲线。

(1)最大悬臂施工阶段计算结果

在最大悬臂施工阶段时六种模型下的主梁上下缘应力图如图9，六种模型下主梁上下缘应力对比图如图10，应力单位为MPa。

由图9、10可见，在梯度升温的情况下，主梁上下缘均出现压应力，且上缘应力比下缘应力大，其原因主要是由于从梁底到梁顶温度逐步升高所引起的。由图10可见，模型2与模型3引起的主梁上缘应力基本一致，而引起起的下缘应力却大不相同，原因是模型3与模型2相比，区别仅在于模型3考虑了底板的升温，而模型2未考虑。模型3、4、5引起的下缘应力连续性较差，其原因是这三个模型均考虑了底板的升温。

最大悬臂阶段，六种温度梯度模型引起的主梁1/4截面的应力沿梁高的分布如图11所示。

表1 跨中合拢口长度缩短实测值与计算结果对比(单位mm)

实测值	模型1	模型2	模型3	模型4	模型5	模型6
							10.25	3.29	5.22	5.12	3.64	13.84	13.78

由图12及表1可见，不同温度梯度模型对主梁变形也有较大影响。尤其对最大悬臂端，模型5与模型1引起的端部竖向位移相差37mm多，而实测端部竖向位移也有26mm。由梯度升温引起的主梁水平伸长量，模型5最大，为13.84mm；模型1最小，为3.29mm。而实测值为10.25mm，与实测值最为接近的是模型6，为13.78mm。如果设计时对梯度温差考虑不足的话，很可能会造成梁端部下挠大于设计计算值，从而造成全桥合龙及线形控制的困难。因此，在连续刚构桥的施工过程中，针对梯度温差的影响同步进行施工控制是十分必要的。

(2)成桥阶段计算结果

在成桥阶段时六种模型下的主梁上下缘应力图如图13，六种模型下主梁上下缘应力对比图如图14，应力单位为MPa。

由图13、14可见，在成桥阶段，六种模型下主梁上缘均为压应力，且最大压应力出现的位置均为跨中或边跨根部。主梁下缘应力出现两种情况，一种是几乎所有截面均受拉，如模型1、模型2、模型6；另一种是部分截面受压部分截面受拉，如模型3、模型4、模型5。其原因是模型3、4、5均考虑了底板升温，从而减小了下缘的拉应力，使得部分截面的下缘出现了压应力，这对于混凝土抗裂设计是偏于不安全的。模型5与模型6主梁下缘出现的拉应力较其他模型的大，其原因是这两个模型的顶底板之间的温差较大，且升温面积较大。由图15可见，梯度升温引起的主梁最大拉应力出现在模型6的中跨跨中截面，为4.54MPa。这种较大的主拉应力通常是引起混凝土箱梁开裂的主要原因，因此，为了有效地防止混凝土箱梁裂缝的出现，减小温度影响是十分必要的。

表2 成桥阶段梯度升温引起的墩顶弯矩(单位：kN·m)

模型	1	2	3	4	5	6
							弯矩	6254.1	10000.5	9326.07	6508.81	26055.9	25302.2

六种模型下主梁竖向挠度与引起墩顶弯矩规律一致，由大到小依次为模型5、6、2、3、4、1。由表2知，主梁梯度温差引起的主墩墩顶弯矩最大达到了26055.9kN·m，与最小的相差了将近20000kN·m，可见对于连续刚构桥，主梁的梯度温差对主墩的内力影响是很大的。

以上计算结果显示，按照我国公路桥梁规范(JTJ023-85)建立的温度梯度模型1所引起的连续刚构桥的应力及位移值与实测值以及其他模型结果相比均偏小，这就预示着按照我国公路桥梁规范(JTJ023-85)设计施工的混凝土连续刚构桥，由于设计时对日照温差的影响考虑不足可能造成运营时出现实际应力过大而引起混凝土箱梁开裂或跨中下挠过大的问题。现实运营中的某些连续刚构桥确实出现了箱梁开裂和跨中下挠过大的问题，所以在分析这类问题时可将温度的影响作为一个方面的因素来考虑。

系统升降温20℃引起的主梁位移如图17、18所示。在图17中，工况1为自重、二期恒载、收缩徐变、预应力及施工荷载的共同作用；工况2为工况1+体系升温；工况3为工况1+体系降温。

由图17、18可见，体系降温引起主梁竖向位移为：中跨跨中下挠，边跨靠近边墩1/4处附近上拱。当系统降温20℃时，中跨跨中下挠最大值可达21.66mm，边跨最大上拱1.56mm，当系统降温时主梁变形的总趋势为下挠；体系升温与体系降温引起的主梁位移趋势正好相反。

在图19中工况1：合计(自重、二期恒载、收缩徐变、预应力及施工荷载的共同作用)+体系升温20℃；工况2：合计+体系降温20℃；工况3：工况1+梯度升温(模型6)；工况4：工况2+梯度升温(模型6)。由图19可见，不论体系升温还是降温，当有梯度升温与其共同作用时，都将引起主梁竖向挠度增大，且当体系降温与梯度升温共同作用时竖向挠度最大。

体系升温与体系降温引起的主梁上下缘应力大小相等，方向相反。体系升温使得主梁中跨部分下缘受压，上缘受拉；边跨部分下缘受拉，上缘受压。系统降温时与之相反。由图21可见，当体系降温20℃时，可在中跨跨中下缘引起将近1.0MPa的拉应力，当其与梯度升温共同作用于连续刚构时，在中跨跨中下缘引起的拉应力将是十分惊人的，也是引起箱梁下缘开裂的主导因素，因此，在设计及施工控制中对这种荷载组合一定要引起足够的重视。

通过上述实施例，完全有效地实现了本发明的目的。该领域的技术人员可以理解本发明包括但不限于附图和以上具体实施方式中描述的内容。虽然本发明已就目前认为最为实用且优选的实施例进行说明，但应知道，本发明并不限于所公开的实施例，任何不偏离本发明的功能和结构原理的修改都将包括在权利要求书的范围中。

Claims (7)

1.一种预应力混凝土连续刚构桥温度作用有限元分析方法，所述预应力混凝土连续刚构桥包括主梁、纵梁、跨梁和桥墩，所述主梁为预应力混凝土连续刚构结构，所述桥墩与主梁刚性连接，其特征在于，所述有限元分析方法具体包括如下步骤，

SS1.定义有限元模型的长度单位、力的单位和温度单位；

SS2.定义纵梁和桥墩的梁单元类型；

SS3.建立有限元计算物理模型，对全桥进行单元划分，设定边界约束；

SS4.定义混凝土和预应力钢材的材料参数，包括弹性模量、容重、线膨胀系数、剪切模量、泊松比、预应力管道偏差系数、预应力钢筋与管道的摩阻系数、预应力钢筋松弛率、锚具变形及钢束回缩值；

SS5.设置边界条件：设定温度梯度模型及加载，加载包括恒载、预应力荷载、徐变和干缩、挂篮荷载、温度荷载；

SS6.建立偏微分控制方程组并对其进行离散化，得到代数方程组并利用步骤SS5设定的边界条件对所述代数方程组进行封闭，用以计算应力和位移；

SS7.对计算域内的所述代数方程组反复进行迭代计算，直到满足所设定的迭代精度为止，得到应力和位移分布。

2.根据权利要求1所述的有限元分析方法，其特征在于，还包括步骤SS8，对计算结果进行后处理，绘制出相关曲线。

3.根据权利要求1所述的有限元分析方法，其特征在于，步骤SS1中，定义有限元模型的长度单位为m、力的单位为牛(N)和温度单位为摄氏度(℃)。

4.根据权利要求1所述的有限元分析方法，其特征在于，步骤SS3中，包括以下三个子步骤：

a.建立有限元计算物理模型：纵梁和桥墩均采用空间梁单元模拟，整体坐标系以X向为桥纵向，Y向为桥横向，Z向为竖向；

b.对全桥进行单元划分：全桥共划分为m+n个单元，其中主梁划分为m个单元，桥墩划分为n个单元，其中，m、n为自然数；

c.设定边界约束：主梁与桥墩的约束关系通过刚性连接模拟；两边跨梁端只有Y-Z平面内的角位移和水平线位移2个自由度，其余4个方向均被约束；不考虑桩土作用，将墩底直接固结。

5.根据权利要求1所述的有限元分析方法，其特征在于，步骤SS4具体为：

a.主梁采用箱梁结构，定义箱梁及桥墩采用的混凝土的弹性模量、剪切模量、泊松比、轴心抗压强度标准值、轴心抗拉强度标准值、线膨胀系数；

b.定义预应力钢材的材料参数，所述预应力钢材包括纵向预应力钢材和竖向预应力钢材，其中，

--所述纵向预应力钢材采用低松弛钢绞线，定义所述低松弛钢绞线的弹性模量、抗拉强度标准值、张拉控制应力、预应力钢束与管道的摩阻系数、预应力管道偏差系数、预应力钢束松弛率、锚具变形及钢束回缩值；优选地，纵向预应力钢材采用Φj15.24mm低松弛钢绞线，技术标准符合ASTM A416-97(270级)标准要求；

--竖向预应力钢材采用预应力粗钢筋，定义预应力粗钢筋的弹性模量、抗拉强度标准值、张拉控制应力、预应力钢筋与管道的摩阻系数、预应力管道偏差系数、预应力钢筋松弛率、锚具变形及钢束回缩值；优选地，竖向预应力钢材采用32mm的精轧粗钢筋。

6.根据权利要求1所述的有限元分析方法，其特征在于，步骤SS5按照如下方式设置温度梯度模型及加载：

a.设置温度梯度模型

选择不同的温度梯度模型，以计算预应力混凝土连续刚构桥在施工过程中，最大悬臂阶段在不同施工阶段时温度应力及应变的变化情况，从以下六种温度梯度模型进行选择，

1)温度梯度模型1：桥面板均匀升温，温差为5℃，箱梁腹板和底板无温差；

2)温度梯度模型2：桥面板表面的最高温度取20℃；

3)温度梯度模型3：按照美国AASHTO规范对温度梯度的规定建立的温度梯度模型；

4)温度梯度模型4：按照英国BS5400规范升温时的温度梯度建立的模型；

5)温度梯度模型5：温度梯度是一条高1200mm的五次抛物线，混凝上表面的温度取32℃，在截面厚度为200mm的底板上采用从0℃到1.5℃的线性温度增长；

6)温度梯度模型6：根据实测结果提出的温度梯度模型，温度梯度曲线为T_y＝T₀e^-αy，T₀为顶板温度，α按最小二乘法进行非线性拟合，y为高度；优选地，根据边跨、中跨1/4截面实测温度梯度值，按最小二乘法进行非线性拟合得到α＝1.57，高度y以米计，顶板温度T₀为20℃；

b.设置的荷载包括恒载、预应力荷载、徐变和干缩、挂篮荷载、温度荷载。

7.根据权利要求1所述的有限元分析方法，其特征在于，步骤SS6中，按照如下方式建立偏微分控制方程组并对其进行离散化：

(1)梁单元刚度矩阵

在局部坐标系下，非线性有限元的梁单元刚度矩阵可以表示为：

[K_T]_B＝[K_E]_B+[K_G]_B      (1)

其中：

[K_T]_B表示局部坐标系下的梁单元切线刚度矩阵；

[K_E]_B表示局部坐标系下的梁单元弹性刚度矩阵；

[K_G]_B表示局部坐标系下的梁单元几何刚度矩阵。

对于梁单元弹性刚度矩阵来说，有：

其中：EI_Y、EI_Z为沿Y轴、Z轴的抗弯刚度；GI_X为沿X轴的抗扭刚度；

应用虚功原理可以推导出非线性有限元平衡方程如下：

[K]{U}＝{P}      (4)

其中，[K]表示结构的整体刚度；

{U}表示全部自由度的位移向量；

{P}表示荷载向量；

(2)温度荷载

混凝土结构内部的温度场是确定温度荷载的关键。

箱梁结构的主梁体内任意一点的温度T是坐标x，y，z和时间t的函数，设定混凝土为均质、各向同性、无内热源，得三维非稳态导热方程：

$\mathrm{\λ}(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{{\∂y}^2}+\frac{{\∂}^{2}T}{{\∂z}^2})=\mathrm{c\γ}\frac{\∂T}{\∂t}---\left(5\right)$

式中：λ——混凝土的导热系数；

c——混凝土的比热；

γ——混凝土的容重；

略去桥长方向温差的影响；

在梁高较小时，略去水平方向很小的热传导作用，用垂直方向的一维热传导状态来分析；

对于梁高较大的箱梁，忽略角隅区附近的热传导状态，用垂直和水平两个方向各自的一维导热状态分别计算，然后再叠加起来；

在箱梁结构计算中简化为一维热传导方程：

$\mathrm{\α}\frac{{\∂}^{2}T}{{\∂x}^2}=\frac{\∂T}{\∂t}---\left(6\right)$

第一类边界条件：混凝土表面温度T是时间的已知函数，

当t＝0时，T(t)＝f(t)      (7)

对于壁板结构，将其近似为一块半无限厚板，并假定气温变化为谐波形式的情况下，根据式(6)可得上述边界条件的弹性力学解：

$T(x,t)=A_0{e}^{-\sqrt{\mathrm{\ω}/2\mathrm{\αx}}}\sin (\mathrm{\ωt}-\sqrt{\mathrm{\ω}/2\mathrm{\αx}})---\left(8\right)$

式中：A₀——壁板表面温度波动峰值；

α——热扩散系数；

ω——圆频率2π/24；

x——计算点至板表面距离；

t——时间；

以某一特定时刻最大温差分布相应的温差荷载作为控制荷载，式(8)表示为温度分布包络线的形式：

$T\left(x\right)=A_0{e}^{-\sqrt{\mathrm{\ω}/2\mathrm{\αx}}}---\left(9\right)$

用函数式表示沿板厚的温度分布，C_x为实验参数；对高架桥箱梁模型采用来计算沿板厚的温度分布；

对混凝土箱形桥墩，采用下式来分析壁厚方向的温度分布：

$T\left(x\right)=T_0{e}^{-C_{X}X}---\left(10\right)$

其中T₀为墩壁内外表面的温差；沿箱梁高、梁宽方向的温差分布按下式计算：

$T\left(y\right)=T_{0y}{e}^{-C_{Y}Y},T\left(x\right)=T_{0X}{e}^{-C_{X}X}---\left(11\right)$

式中：T_0y、T_0X——沿梁高、梁宽方向的温差；

y、X——计算点至受热表面的距离；

C_X、C_y——指数系数，随结构形式、部位、计算时刻而异；

以上式(9)至式(11)，都是采用第一类边界条件、考虑温度分布包络线的函数式，影响第一类边界条件的外部主要因素为太阳辐射强度、气温变化、风速；

混凝土内部的温度分布确定后，根据混凝土的热物理性能，利用线膨胀系数，形成温度荷载，考虑桥面受日照后形成的沿箱梁高度变化的温度梯度T_y＝T₀e^-αy；

(3)温差应力

确定温度梯度模型及温度设计值后，温度应力按结构力学和有限元方法进行计算，计算时假定：

1)桥长方向的温度分布是均匀的；

2)混凝土是弹性匀质材料；

3)梁变形服从平面假定；

4)按单向温差荷载计算温差应力，然后叠加组合多向温差荷载状下的温差应力；

温度应力有由两部分组成：a)梁的温度变形受到纵向纤维之间的相互约束，在截面上产生自平衡的纵向约束应力，即自应力；b)梁的温度上拱变化受到支承条件约束的温度次应力；

①箱梁温度自应力

设温度梯度沿梁高按任意曲线t(y)分布，取单位梁长dy＝1的微分段，当纵向纤维之间不受约束且自由伸缩时，沿梁高的自由应变ε_t(y)与温度梯度一致，即：

ε_t(y)＝α_ct_(y)      (12)

由于纵向纤维之间的相互约束，梁截面应变应符合平面假定，梁截面上的最终应变ε_f(y)为直线分布，即：

ε_f(y)＝ε₀+ψy      (13)

式中ε₀——基轴y＝0处应变；

ψ——截面变形曲率；

y——基轴以下任一点求应变的坐标；

α_c——混凝土线膨胀系数；

自由应变与最终应变之差，系纤维之间的约束产生，其值为：

ε_σ(y)＝ε_t(y)-ε_f(y)＝α_ct_(y)-(ε₀+ψy)      (14)

自应力为：

σ_s(y)＝E_cε_σ(y)＝E_c[α_ct_(y)-(ε₀+ψy)]      (15)

全截面上轴力N和弯矩M

$\begin{array}{c}N=E_c\underset{h}{\∫}{\mathrm{\ϵ}}_{\mathrm{\σ}\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}=E_c\underset{h}{\∫}({\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy})b_{\left(y\right)}\mathrm{dy}\\ =E_c[{\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}-\mathrm{\ψ}\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}]\end{array}---\left(16\right)$

$\begin{array}{c}M=E_c\underset{h}{\∫}{\mathrm{\ϵ}}_{\mathrm{\σ}\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=E_c\underset{h}{\∫}({\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy})b_{\left(y\right)}(y-y_c)\mathrm{dy}\\ =E_c[{\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}-{\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}-\mathrm{\ψ}\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}]\end{array}---\left(17\right)$

式中E_c——混凝土材料弹性模量；

b_(y)——y处的梁宽；

对于任何截面，N＝0，M＝0，即内力总和为零；

公式(16)、(17)分别改写为：

${\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}+\mathrm{\ψ}\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}---\left(18\right)$

${\mathrm{\ϵ}}_0\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}+\mathrm{\ψ}\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}---\left(19\right)$

在公式(18)、(19)内

$\underset{h}{\∫}{b}_{\left(y\right)}\mathrm{dy}=A---\left(20\right)$

$\underset{h}{\∫}{\mathrm{yb}}_{\left(y\right)}\mathrm{dy}={\mathrm{Ay}}_c---\left(21\right)$

$\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}=\underset{h}{\∫}{b}_{\left(y\right)}{y}^2\mathrm{dy}-\underset{h}{\∫}{b}_{\left(y\right)}{\mathrm{yy}}_c\mathrm{dy}=I_b-\underset{h}{\∫}{b}_{\left(y\right)}{\mathrm{yy}}_c\mathrm{dy}=I_g---\left(22\right)$

$\underset{h}{\∫}{b}_{\left(y\right)}(y-y_c)\mathrm{ydy}=0$ (对重心轴的静面积矩为零)

式中A——截面面积；

I_b——截面面积对基轴惯性矩；

I_g——截面面积对重心轴惯性矩；

将公式(20)～(22)代入公式(18)、(19)内，得：

${\mathrm{\ϵ}}_{0}A+\mathrm{\ψ}{\mathrm{Ay}}_c={\mathrm{\α}}_c\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}---\left(23\right)$

ψI_g＝α_c∫t_(y)b_(y)(y-y_c)dy      (24)

由公式(23)、(24)可得：

${\mathrm{\ϵ}}_0=\frac{{\mathrm{\α}}_c}{A}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c---\left(25\right)$

$\mathrm{\ψ}=\frac{{\mathrm{\α}}_c}{I_g}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}---\left(26\right)$

设在坐标y处，截面内一厚度为i的微小单元面积A_y处温度梯度值为t_y，以t_y为常值代入公式(25)、(26)，积分区段仅在i厚度范围内有值：t_(y)＝t_y,y-y_c＝e_y(单元面积A_y对全面积重心的偏心距)。

${\mathrm{\ϵ}}_0=\frac{{\mathrm{\α}}_c}{A}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c=\frac{{\mathrm{\α}}_c}{A}\underset{i}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}\mathrm{dy}-{\mathrm{\ψy}}_c=\frac{{\mathrm{\α}}_c{t}_y{A}_y}{A}-\frac{{\mathrm{\α}}_c{t}_y{A}_y{e}_y{y}_c}{I_g}---\left(27\right)$

$\mathrm{\ψ}=\frac{{\mathrm{\α}}_c}{I_g}\underset{h}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=\frac{{\mathrm{\α}}_c}{I_g}\underset{i}{\∫}{t}_{\left(y\right)}{b}_{\left(y\right)}(y-y_c)\mathrm{dy}=\frac{{\mathrm{\α}}_c{t}_y{A}_y{e}_y}{I_g}---\left(28\right)$

自公式(15)可求得任意点应力σ_s(y)：

${\mathrm{\σ}}_{s\left(y\right)}=E_c[{\mathrm{\α}}_c{t}_{\left(y\right)}-({\mathrm{\ϵ}}_0+\mathrm{\ψy})]=E_c{\mathrm{\α}}_c{t}_y-\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y}{A}+\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y{e}_y{y}_c}{I_g}-\frac{E_c{\mathrm{\α}}_c{t}_y{A}_y{e}_{y}y}{I_g}---\left(29\right)$

令：N_ti＝A_yt_yα_cE_c,M_ti＝-N_tie_y＝-A_yt_yα_cE_ce_y

${\mathrm{\σ}}_{s\left(y\right)}=-\frac{N_{\mathrm{ti}}}{A}+\frac{M_{\mathrm{ti}}}{I_g}(y-y_c){+t}_y{\mathrm{\α}}_c{E}_c---\left(30\right)$

公式(30)是由于一个单元面积A_y内的温度作用，在截面任一点产生的应力；对于分为很多块单元面积上不同t_y的作用，应用分段总和法；公式(30)使用于正温差；如为反温差则整个公式前冠以负号；

②超静定结构中的温度次内力及其次应力

在预应力混凝土超静定结构中，前述温度变形ε₀及曲率ψ将受到超静定赘余约束的制约，引起温度次内力，两端固定杆单元的节点荷载向量{F}^e由截面变形曲率及沿梁高y＝0处的变形ε₀直接写出：

${\left\{F\right\}}^e=\left\{\begin{array}{c}{N}_i\\ Q_i\\ M_i\\ N_j\\ Q_j\\ M_j\end{array}\right\}=\left\{\begin{array}{c}\mathrm{EA}({\mathrm{\ϵ}}_0+{\mathrm{\ψy}}_c)\\ 0\\ \mathrm{EI\ψ}\\ -\mathrm{EA}({\mathrm{\ϵ}}_0+{\mathrm{\ψy}}_c)\\ 0\\ -\mathrm{EI\ψ}\end{array}\right\}---\left(31\right)$

杆件单元节点力应以结构坐标系表示，然后分别组集各个杆件单元的结点荷载，从而得到节点外力向量{F}，矩阵位移方程为：

[K]{Δ}+{F}＝0      (32)

式中[K]——结构总刚度矩阵；

{Δ}——单元节点位移向量；

在求得结构各单元因温度变化引起的结点位移后，由单元的杆端力与单元刚度矩阵、单元结点位移的关系{f}^e＝[K]{Δ}^e求得结构的温度次内力N_T、Q_T、M_T及其次应力；在超静定结构中，总的温度应力为：

纵向弯曲应力：

${\mathrm{\σ}}_{t\left(y\right)}=\frac{N_T}{A}+\frac{M_T}{I}y+E_c[{\mathrm{\α}}_c{t}_{\left(y\right)}-{\mathrm{\ϵ}}_0-\mathrm{\ψy}]---\left(33\right).$