文章快速检索    
  同济大学学报(自然科学版)  2017, Vol. 45 Issue (8): 1108-1113.  DOI: 10.11908/j.issn.0253-374x.2017.08.002
0

引用本文  

陈亮. 组合梁几何非线性与长期效应的同步算法[J]. 同济大学学报(自然科学版), 2017, 45(8): 1108-1113. DOI: 10.11908/j.issn.0253-374x.2017.05.002.
CHEN Liang. Algorithm of Composite Beam Synchronously Considering Geometrical Nonlinearity and Long Term Effects[J]. Journal of Tongji University (Natural Science), 2017, 45(8): 1108-1113. DOI: 10.11908/j.issn.0253-374x.2017.08.002.

基金项目

国家“九七三”重点基础研究发展规划(2013CB036300)

第一作者

陈亮(1975-),男,工学博士,高级工程师,主要研究方向为钢与组合结构桥梁. E-mail: chenliang@smedi.com

文章历史

收稿日期:2016-09-11
组合梁几何非线性与长期效应的同步算法
陈亮1,2    
1. 同济大学 桥梁工程系,上海 200092;
2. 上海市政工程设计研究总院(集团)有限公司,上海 200092
摘要:首先根据建立的随转坐标系下考虑滑移组合梁单元平衡方程,结合初应变法得到组合梁收缩徐变等效节点力有限元列式.然后基于随转列式全量法,利用总体坐标系与随转坐标系下变量间转换矩阵,获得界面滑移平面组合梁单元收缩徐变与几何非线性同步考虑的实用算法.该方法的优点是可实现几何非线性与收缩徐变材料非线性算法上的分离,收缩徐变仅在随转坐标系下考虑,几何非线性通过转换矩阵考虑.最后,结合算例说明考虑滑移的组合梁几何非线性与收缩徐变之间存在着耦合作用.
关键词组合梁    界面滑移    收缩徐变    几何非线性    有限元    
Algorithm of Composite Beam Synchronously Considering Geometrical Nonlinearity and Long Term Effects
CHEN Liang1,2     
1. Department of Bridge Engineering, Tongji University, Shanghai 200092, China;
2. Shanghai Municipal Engineering Design Institute (Group) Co., Ltd., Shanghai 200092, China
Abstract: Firstly, balance equation of composite beam element with interlayer slips under corotational coordinate system was deduced. Combining initial strain method, equivalent nodal force formulation of shrinkage and creep was obtained. Based on the corotational method, using the transformation matrices relating local and global quantities, a practical algorithm for analysis of composite beam with interlayer slips considering geometric nonlinearity, concrete shrinkage and creep synchronously was put forward. The advantage of this approach is that it leads to the separation of geometrical and material non-linearities such as shrinkage and creep. The shrinkage and creep is only present at the level of the local element, whereas the geometrical non-linearity is included in the transformation matrices. Finally, it was illustrated that the geometric nonlinearity is coupled with shrinkage and creep for composite beam with interlayer slip by an example.
Key words: composite beam    interlayer slips    shrinkage and creep    geometric nonlinearity    finite element    

随着科学技术的发展,组合梁、柱在大跨斜拉桥、高层建筑等柔性结构中有越来越多的应用.组合梁、柱钢混界面之间绝对刚性的连接件并不存在,另外混凝土材料不可避免存在收缩徐变现象,因此针对此类柔性组合结构收缩徐变问题的研究具有重要的现实意义.目前,有较多文献研究考虑滑移的组合梁几何非线性有限元列式[1-2],且有较多文献研究考虑滑移的组合梁收缩徐变效应[3-4].根据普通混凝土梁的研究,发现几何非线性与收缩徐变之间存在明显的耦合效应[5].但是,根据笔者掌握的资料,尚未发现针对界面滑移组合梁同时考虑几何非线性、收缩徐变方面的研究报道.本文针对考虑滑移的组合梁,采用基于Newmark模型的有限单元法[6-7],采用随转(CR)列式法[8]考虑几何非线性效应,采用基于初应变法的逐步计算方法[9]考虑收缩徐变效应,获得界面滑移平面组合梁同步考虑几何非线性以及收缩徐变的实用算法,并给出了详细的计算步骤.最后,利用自编程序CBS(组合桥梁分析系统),在文献[4]算例基础上,进行相关分析说明.

1 随转坐标系下单元列式 1.1 单元位移

图 1所示的典型的钢-混组合梁截面,顶板是混凝土桥面板,下托梁是钢梁,两者通过剪力连接件相连.本文计算采用如下基本假定:① 不考虑组合梁的竖向掀起;② 忽略梁体剪切变形影响;③ 剪力连接件沿梁均匀分布.图中随转局部坐标xl轴为顺桥向,yl轴为横桥向,zl轴为竖向.组合单元桥面板和钢主梁截面上任一点的顺桥向、竖向位移uαvα可分别表示为:uα=uα+(zzα)θvα=v;组合梁纵向滑移g可表达为:g=ubua.其中:下标α=a, b,分别表示下层钢主梁及上层桥面板分部(下同),uα表示各分部形心轴处顺桥向的位移,竖向转角θ=-dv/dx,上下分部形心间距h=zbzazbza分别表示上、下分部形心的竖向坐标.

图 1 组合梁几何模型以及位移 Fig.1 Geometric model and displacement of composite beam

本文采用10自由度的单元模式,如图 2所示.单元节点位移矢量为:p=[ub1 ua1 v1 θ1 ub2 ua2 v2 θ2 ub3 ua3]T.引入位移矢量u,则单元内各点位移可由式(1) 表示:

图 2 组合梁单元 Fig.2 Element of composite beam
$\mathit{\boldsymbol{u}}={{\left[ {{{\bar{u}}}_{\rm{b}}}\ \ \ {{{\bar{u}}}_{\rm{a}}}\ \ \ v \right]}^{\rm{T}}}=\mathit{\boldsymbol{Np}}$ (1)

其中形函数为

$\begin{array}{l} \mathit{\boldsymbol{N}} = \left[ \begin{array}{l} {\mathit{\boldsymbol{N}}_{{\rm{ub}}}}\\ {\mathit{\boldsymbol{N}}_{{\rm{ua}}}}\\ {\mathit{\boldsymbol{N}}_{\rm{v}}} \end{array} \right]{\rm{ = }}\\ \;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {{t_1}}&0&0&0&{{t_2}}&0&0&0&{{t_3}}&0\\ 0&{{t_1}}&0&0&0&{{t_2}}&0&0&0&{{t_3}}\\ 0&0&{{s_1}}&{{s_2}}&0&0&{{s_3}}&{{s_4}}&0&0 \end{array}} \right] \end{array}$

式中:NubNua分别为上层分部、下层分部顺桥向位移形函数;Nv为竖向位移形函数.

t1=1-3ξ+2ξ2t2=-ξ+2ξ2t3=4ξ-4ξ2s1=1-3ξ2+2ξ3s2=--2ξ2+ξ3s3=3ξ2-2ξ3s4=-L(ξ3ξ2),而ξ=x/LL为单元长度.

1.2 平衡方程

组合梁单元结点力矢量可表示为:f=[Nb1 Na1 V1 M1 Nb2 Na2 V2 M2 Nb3 Na3]T,根据虚功原理,可得平衡方程如下:

${{\mathit{\boldsymbol{\tilde p}}}^{\rm{T}}}f = \int_L {\left( {\sum\nolimits_\alpha {\int_{A\alpha } {{\sigma _{xa}}{{\tilde \varepsilon }_{xa}}} {\rm{d}}{A_\alpha } + {f_s}\tilde g} } \right){\rm{d}}x} $ (2)

上标“~”表示虚位移或虚应变,ε(z, x)=ε(x)+(zzα)к(x),其中, ε(x)为分部α形心处的轴向应变,к=d2v/dx2, 为组合梁曲率.对线弹性单元,钢、混凝土材料及剪力连接件的本构关系为σ=Eαεfs=ksg,其中, σεEα分别为分部α中任一点的轴向应力、轴向应变及分部材料的弹性模量, fsks分别为界面滑移力及界面连接刚度,下标s表示滑移.

引入广义弹性应变矢量εeT=[εxb εxa к g],引入广义应力矢量rT=[Nb Na Mab fs],其中:${N_a} = \int_{{A_a}} {{\sigma _{xa}}{\rm{d}}{A_a}} ,{M_a} = \int_{{A_a}} {{\sigma _{xa}}(z - {z_a}){\rm{d}}{A_a},{M_{{\rm{ab}}}} = {M_{\rm{a}}} + {M_{\rm{b}}}} $

应变位移关系为

${\varepsilon _{\rm{e}}}{\rm{ = }}D\mathit{\boldsymbol{u}} = \mathit{\boldsymbol{DNp}} = \mathit{\boldsymbol{Bp}}$ (3)

其中,

$D=\left[ \begin{matrix} \partial & 0 & 0 \\ 0 & \partial & 0 \\ 0 & 0 & {{\partial }^{2}} \\ 1 & -1 & h\partial \\ \end{matrix} \right],\mathit{\boldsymbol{B}}=D\mathit{\boldsymbol{N}}=\left[ \begin{matrix} \mathit{\boldsymbol{N}}_{\rm{ub}}^{'} \\ \mathit{\boldsymbol{N}}_{\rm{ua}}^{'} \\ \mathit{\boldsymbol{N}}_{\rm{v}}^{''} \\ {{\mathit{\boldsymbol{N}}}_{\rm{ub}}}-{{\mathit{\boldsymbol{N}}}_{\rm{ua}}}+h\mathit{\boldsymbol{N}}_{\rm{v}}^{'} \\ \end{matrix} \right]$

则式(2) 可写为如下紧凑形式:

${{\mathit{\boldsymbol{\tilde{p}}}}^{\rm{T}}}\mathit{\boldsymbol{f=}}\int_{L}{{\mathit{\boldsymbol{\tilde{ }\!\!\varepsilon\!\!\rm{ }}}}}_{\rm{e}}^{\rm{T}}\mathit{\boldsymbol{r}}\rm{d}\mathit{x=}{{\mathit{\boldsymbol{\tilde{p}}}}^{\rm{T}}}\int_{L}{{{\mathit{\boldsymbol{B}}}^{\rm{T}}}}\mathit{\boldsymbol{r}}\rm{d}\mathit{x}$ (4)

式(4) 对于任何虚位移都必须成立,于是存在如下关系:

$\mathit{\boldsymbol{f = }}\int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} \mathit{\boldsymbol{r}}{\rm{d}}\mathit{x}$ (5)

考虑到仅桥面板发生收缩、徐变,初始收缩广义应变εsT=[εs 0 0 0],初始徐变广义应变εcT=[εc 0 кc 0],其中εs为收缩引起的桥面板形心处轴向应变,εcкc分别为徐变引起的桥面板形心处轴向应变及弯曲应变.令:

$\mathit{\boldsymbol{D = }}\left[ {\begin{array}{*{20}{c}} {E{A_{\rm{b}}}}&0&0&0\\ 0&{E{A_{\rm{a}}}}&0&0\\ 0&0&{E{I_{{\rm{ab}}}}}&0\\ 0&0&0&{{k_{\rm{s}}}} \end{array}} \right],$
${\mathit{\boldsymbol{D}}_{\rm{b}}}\mathit{\boldsymbol{ = }}\left[ {\begin{array}{*{20}{c}} {E{A_{\rm{b}}}}&0&0&0\\ 0&0&0&0\\ 0&0&{E{I_{\rm{b}}}}&0\\ 0&0&0&{{k_{\rm{s}}}} \end{array}} \right]$

其中,

$E{A_\alpha } = \int_{{A_\alpha }} {{E_a}{\rm{d}}{A_\alpha }} ,EI\alpha = \int_{{A_\alpha }} {{E_a}{{(y - {y_\alpha })}^2}{\rm{d}}{A_\alpha }} $,且有,EIab=EIa+EIb.则总广义应力表达为

$r = \left[ {\begin{array}{*{20}{c}} {E{A_{\rm{b}}}{{\bar \varepsilon }_{x{\rm{b}}}} - E{A_{\rm{b}}}{{\bar \varepsilon }_{\rm{s}}} - E{A_{\rm{b}}}{{\bar \varepsilon }_{\rm{c}}}}\\ {E{A_{\rm{a}}}{{\bar \varepsilon }_{x{\rm{a}}}}}\\ {E{I_{{\rm{ab}}}}\kappa - E{I_\mathit{b}}{\kappa _{\rm{c}}}}\\ {{k_{\rm{s}}}g} \end{array}} \right] = \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\varepsilon }}_{\rm{e}}} - {\mathit{\boldsymbol{D}}_{\rm{b}}}{\mathit{\boldsymbol{\varepsilon }}_{\rm{s}}} - {\mathit{\boldsymbol{D}}_{\rm{b}}}{\mathit{\boldsymbol{\varepsilon }}_{\rm{c}}}$ (6)

将式(6) 代入式(5),得:

$\begin{array}{l} \mathit{\boldsymbol{f = }}\int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} \mathit{\boldsymbol{DB}}{\rm{d}}\mathit{x}\mathit{\boldsymbol{p}} - \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_{\rm{b}}}{\mathit{\boldsymbol{\varepsilon }}_{\rm{s}}}{\rm{d}}\mathit{x} - \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_{\rm{b}}}{\mathit{\boldsymbol{\varepsilon }}_{\rm{c}}}{\rm{d}}\mathit{x} = \\ \;\;\;\;\;\;{\mathit{\boldsymbol{f}}_{\rm{e}}} + {\mathit{\boldsymbol{f}}_{\rm{s}}} + {\mathit{\boldsymbol{f}}_{\rm{c}}} \end{array}$ (7)

对于式(7),显然$\mathit{\boldsymbol{f}}_{\rm{e}} = \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{DB}}{\rm{d}}x\mathit{\boldsymbol{p}}} = \mathit{\boldsymbol{Kp}}$为随转坐标系下节点位移引起的杆端力,$\mathit{\boldsymbol{f}}_{\rm{s}} =- \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{D}}_{\rm{b}}\mathit{\boldsymbol{\varepsilon}}_{\rm{s}}{\rm{d}}x} $,为随转坐标系下收缩初应变引起的杆端力,$\mathit{\boldsymbol{f}}_{\rm{c}} =- \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{D}}_{\rm{b}}\mathit{\boldsymbol{\varepsilon}}_{\rm{c}}{\rm{d}}x} $,为随转坐标系下徐变初应变引起的杆端力.

1.3 收缩徐变逐步计算法

在混凝土分部常应力rb作用下,以徐变系数来表示混凝土徐变应变的表达式为

${\mathit{\boldsymbol{\varepsilon }}_{\rm{c}}}{\rm{ = }}{\mathit{\boldsymbol{r}}_{\rm{b}}}/{\mathit{\boldsymbol{D}}_{\rm{b}}}\varphi \left( {t,\tau } \right)$ (8)

式中:εcT为时刻τ始作用于混凝土的常应力rb至时刻t所产生的徐变应变;φ(t, τ)为加载龄期为τ; 计算龄期为t的混凝土徐变系数.

由于大跨径桥梁都是采用分节段多阶段施工完成,一般采用逐步计算法来进行徐变分析,将整个求解时间轴分为t0t1,…,tntn+1等时刻,在各时刻作用有Δrb0,Δrb1,…,Δrbn,Δrbn+1等应力增量.一般而言,一个施工阶段对应于一个计算时步,对运营期等时间较长的阶段,可划分为几个时间步,以提高计算精度.下面以在tn时刻,在Δtn+1,即tntn+1时段内为例,介绍基于初应变法的单元徐变等效节点力逐步计算方法.

该时段内,在结构的应力和材料常数不随时间而变化的假定下,由式(8) 得徐变增量为

$\Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{c}}^{n + 1} = \sum\limits_{i = 0}^n {\Delta _{\rm{b}}^i} /{\mathit{\boldsymbol{D}}_{\rm{b}}}\varphi \left( {{t^{n + 1}},{\mathit{t}^i}} \right)$ (9)

式中:Δrbi为时间轴上ti时刻混凝土桥面板的应力增量;Δφ(tn+1, ti)为加载龄期为ti,在Δtn+1,即tntn+1时段内的徐变系数.

将Δεcn+1作为初应变,计算出其等效结点力增量,并将其作用于Δtn+1时段末,即tn+1时刻,按有限元的基本方法,单元徐变等效结点力增量可表示为

$\begin{array}{l} \Delta f_{b,c}^{n + 1} = - \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{c}}^{n + 1}{\rm{d}}\mathit{x}\\ \;\;\;\;\;\;\;\; = - \sum\nolimits_{i = 0}^n {\int_L {\left[ {{\mathit{\boldsymbol{B}}^{\rm{T}}}\Delta \mathit{\boldsymbol{r}}_{\rm{b}}^i{\rm{d}}x\Delta \varphi \left( {{t^{n + 1}},{t^i}} \right)} \right]} } \end{array}$ (10)

式(10) 中Δrbi可表达为

$\Delta \mathit{\boldsymbol{r}}_{\rm{b}}^i = {\mathit{\boldsymbol{D}}_{\rm{b}}}\left( {\Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{e}}^i - \Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{s}}^i - \Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{c}}^i} \right)$ (11)

其中Δεeiti时刻结点位移增量引起的应变增量,Δεsi、Δεci分别为ti时刻混凝土收缩、徐变初应变增量.

$\begin{array}{l} \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} \Delta \mathit{\boldsymbol{r}}_{\rm{b}}^i{\rm{d}}x{\rm{ = }}\int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{e}}^i{\rm{d}}x - \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{s}}^i{\rm{d}}x - \\ \;\;\;\;\int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{c}}^i{\rm{d}}x{\rm{ = }}\Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{e}}}^i + \Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{s}}}^i + \Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{c}}}^i{\rm{ = }}\Delta \mathit{\boldsymbol{f}}_{\rm{b}}^i \end{array}$ (12)

式中:$\Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{e}}}^i = \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon}} _{\rm{e}}^i{\rm{d}}x} $ti时刻结点位移增量引起的对应于混凝土桥面板分部的杆端力增量,$\Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{s}}}^i = \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon}} _{\rm{s}}^i{\rm{d}}x} $$\Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{c}}}^i = \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon}} _{\rm{c}}^i{\rm{d}}x} $ti时刻由收缩、徐变初应变Δεsi、Δεci引起的固端力增量.

将式(12) 代入式(10) 得:

$\Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{c}}}^{n + 1} = \sum\limits_{i = 0}^n {\Delta \mathit{\boldsymbol{f}}_{\rm{b}}^i} \Delta \varphi \left( {{t^{n + 1}},{t^i}} \right)$ (13)

假定Δεsn+1为Δtn+1,即tntn+1时步内收缩应变增量,令Δεsn+1=[Δεsn+1 0 0 0]T,则对应的等效结点力增量为

$\Delta \mathit{\boldsymbol{f}}_{{\rm{b}},{\rm{s}}}^{n + 1}{\rm{ = }} - \int_L {{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_{\rm{b}}}\Delta \mathit{\boldsymbol{\varepsilon }}_{\rm{s}}^{n + 1}{\rm{d}}x$ (14)
1.4 静力凝聚

为便于单元到系统的组装,使用静力凝聚法[10],可消除中间节点k的自由度.用下标i表示内部中间结点自由度,下标o表示外部结点自由度,下标l表示随转坐标下静力凝聚后相关量,则内部结点位移矢量pi=[ub3 ua3]T,静力凝聚后随转坐标系下单元局部位移矢量pl=po=[ub1 ua1 v1 θ1 ub2 ua2 v2 θ2]T.将单元刚度矩阵K按外部结点自由度和内部结点自由度分块,即$\mathit{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_{{\rm{oo}}}}}& {{\mathit{\boldsymbol{K}}_{{\rm{oi}}}}}\\ {{\mathit{\boldsymbol{K}}_{{\rm{io}}}}}& {{\mathit{\boldsymbol{K}}_{{\rm{ii}}}}} \end{array}} \right]$,则可得对应于交界结点位移的刚度矩阵Kl=KooKoiKii-1Kio.当刚度矩阵Kl中对应于钢梁及连接件分部的刚度项为零时,可得对应于混凝土桥面板分部的刚度矩阵Kl, b.则由结点位移增量可计算得到对应于混凝土桥面板分部的单元内力增量Δfl, b, ei=Kl, bΔPli其中, Δpli为随转坐标系下单元局部位移矢量增量.也可得到对应于单元整体的内力增量Δfl, ei=KlΔPli.

徐变固端力增量同样需要静力凝聚,可以证明,根据静力凝聚后的单元内力,可自然得到静力凝聚后的徐变等效结点力.即有:

$\Delta \mathit{\boldsymbol{f}}_{{\rm{1,b}},{\rm{c}}}^{n + 1} = \sum\limits_{i = 0}^n {\Delta \mathit{\boldsymbol{f}}_{{\rm{1,b}}}^i} \Delta \varphi \left( {{t^{n + 1}},{t^i}} \right)$ (15)
$\Delta \mathit{\boldsymbol{f}}_{{\rm{1,b}},{\rm{e}}}^i + \Delta \mathit{\boldsymbol{f}}_{{\rm{1,b}},{\rm{s}}}^i + \Delta \mathit{\boldsymbol{f}}_{{\rm{1,b}},{\rm{c}}}^i = \Delta \mathit{\boldsymbol{f}}_{{\rm{1,b}}}^i$ (16)

另由于$\int_L {{{t'}_3}{\rm{d}}x = 0} $故收缩应变增量引起的等效结点力增量在中间结点处为零,静力凝聚后收缩等效结点力增量可由式(17) 直接得到:

$\Delta \mathit{\boldsymbol{f}}_{{\rm{1,b}},s}^{n + 1} = {\left[ {E{A_{\rm{b}}}\;\;\;0\;\;\;0\;\;\;0\;\;\; - E{A_{\rm{b}}}\;\;\;0\;\;\;0\;\;\;0} \right]^{\rm{T}}}\Delta \mathit{\varepsilon }_{\rm{s}}^{n + 1}$ (17)
2 总体坐标系下单元列式 2.1 梁体运动描述

图 3图 4表示总体坐标系(x, z)中一个梁单元在初始时刻及t时刻的位置和形状,(xlzl)是t时刻单元CR坐标系.t时刻总体坐标系下位移可表示为:${\mathit{\boldsymbol{p}}_{\rm{g}}} = {[{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} }_{{\rm{a}}1}}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} }_1}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \theta } }_1}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over g} }_1}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} }_{{\rm{a}}2}}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over v} }_2}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \theta } }_2}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over g} }_2}]^{\rm{T}}}$.则t时刻的刚体转角ɸ可由式(18)、(19) 得到:

$\sin \phi = {c_0}s - {s_0}c$ (18)
$\cos \phi = {c_0}c + {s_0}s$ (19)

其中,

${c_0} = \cos {\beta _0} = \frac{1}{{{l_0}}}\left( {{x_{{\rm{a}}2}} - {x_{{\rm{a}}1}}} \right)$
${s_0} = \sin {\beta _0} = \frac{1}{{{l_0}}}\left( {{z_{{\rm{a}}2}} - {z_{{\rm{a2}}}}} \right)$
$c = \cos \beta = \frac{1}{{{l_n}}}\left( {{x_{{\rm{a}}2}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} }_{{\rm{a}}2}} - {x_{{\rm{a}}1}} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} }_{{\rm{a1}}}}} \right)$
$s = \sin \beta = \frac{1}{{{l_n}}}\left( {{z_{{\rm{a1}}}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over v} }_1} - {z_{{\rm{a2}}}} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over v} }_2}} \right)$
${l_0} = {\left[ {{{\left( {{x_{{\rm{a}}2}} - {x_{{\rm{a}}1}}} \right)}^2} + {{\left( {{z_{{\rm{a1}}}} + {z_{{\rm{a}}2}}} \right)}^2}} \right]^{1/2}}$
${l_n} = {\left[ {{{\left( {{x_{{\rm{a}}2}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} }_{{\rm{a}}2}} - {x_{{\rm{a}}1}} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} }_{{\rm{a1}}}}} \right)}^2} + {{\left( {{z_{{\rm{a1}}}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over v} }_1} - {z_{{\rm{a2}}}} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over v} }_2}} \right)}^2}} \right]^{1/2}}$

随转坐标系下位移可由式(20) 计算得到:

$\begin{array}{l} {u_{{\rm{a}}1}}{\rm{ = }}0\\ {v_1} = 0\\ {v_2} = 0\\ {u_{{\rm{a}}2}}{\rm{ = }}{l_n} - {l_0}\\ {\theta _1} = {\theta _1} - \phi \\ {\theta _2} = {\theta _2} - \phi \\ {u_{{\rm{b}}1}}{\rm{ = }}{g_1} + {h_1}\\ {u_{{\rm{b}}2}}{\rm{ = }}{g_2} + {u_{{\rm{a}}2}} + h{\theta _2} \end{array}$ (20)

其中${g_i} = {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over g} }_i}\cos {\theta _i},i = 1,2$.

图 3 随转运动:滑移 Fig.3 Corotational kinematics: slips
图 4 随转运动:平移及转角 Fig.4 Corotational kinematics: displacements and rotations
2.2 单元列式

一旦通过式(20) 求得随转坐标系下位移,由于相对于局部随转坐标系的位移较小,则随转坐标系下单元内力及刚度矩阵可由随转坐标系下几何线性单元定义求得.在某一坐标系下单元内力矢量fx,切线刚度矩阵Kx,位移矢量px,则存在如下关系:

$\delta {\mathit{\boldsymbol{f}}_x} = {\mathit{\boldsymbol{K}}_x}\delta {\mathit{\boldsymbol{p}}_s}$ (21)

两位移矢量pxpy存在如下关系:

$\delta {\mathit{\boldsymbol{p}}_x} = {\mathit{\boldsymbol{B}}_{xy}}\delta {\mathit{\boldsymbol{p}}_y}$ (22)

其中,Bxy可由式(23) 求得:

${\mathit{\boldsymbol{B}}_{xy}} = \frac{{\partial {\mathit{\boldsymbol{p}}_x}}}{{\partial {\mathit{\boldsymbol{p}}_y}}}$ (23)

由两个不同坐标下虚功相等,可得:

$\delta \mathit{\boldsymbol{p}}_y^{\rm{T}}{\mathit{\boldsymbol{f}}_y} = \delta \mathit{\boldsymbol{p}}_x^{\rm{T}}{\mathit{\boldsymbol{f}}_x}$ (24)

结合式(22) 可得:

${\mathit{\boldsymbol{f}}_y} = \mathit{\boldsymbol{B}}_{xy}^{\rm{T}}{\mathit{\boldsymbol{f}}_x}$ (25)

结合式(25) 和式(21),与位移矢量py对应的切线刚度矩阵Ky可表示为

${\mathit{\boldsymbol{K}}_y} = \mathit{\boldsymbol{B}}_{xy}^{\rm{T}}{\mathit{\boldsymbol{K}}_x}{\mathit{\boldsymbol{B}}_{xy}} + {\mathit{\boldsymbol{H}}_{xy}}$ (26)
${\mathit{\boldsymbol{H}}_{xy}} = {\left. {\frac{{\partial \left( {\mathit{\boldsymbol{B}}_{xy}^{\rm{T}}{\mathit{\boldsymbol{f}}_x}} \right)}}{{\partial {\mathit{\boldsymbol{p}}_y}}}} \right|_{{f_x}}}$ (27)

现引入中间位移矢量pm=[ua2 θ1 θ2 g1 g2]T${\mathit{\boldsymbol{p}}_{\rm{n}}} = {[{u_{{\rm{a}}2}}\;\;{\theta _1}\;\;{\theta _2}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over g} }_1}\;\;{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over g} }_2}]^{\rm{T}}}$.可通过对应于4个位移矢量plpmpnpg相关量之间的变换,即flKlfmKmfnKnfgKg,求得总体坐标系下单元刚度矩阵Kg及内力fg,具体计算过程详见文献[2].

3 同时考虑几何非线性与收缩徐变的计算步骤

从第1~2节推导过程,可知同时考虑几何非线性和收缩徐变效应时的基本计算原理为:将不平衡力(由总外荷载减去收缩徐变等效节点力总量得到)作用下经过几何非线性分析得到的变形作为该时刻总体坐标系下总的弹性变形,转换至随转坐标下的变形,与上一时步随转坐标系下弹性变形差即为该时步随转坐标系下弹性变形增量,徐变变形增量与之成线性关系,将收缩徐变变形作为初应变,计算得到的等效节点力增量总是作用在单元随转坐标系中,因此在各阶段的累加计算过程中无须考虑坐标系的转换,可直接累加形成收缩徐变等效节点力总量,将收缩徐变等效节点力总量由结构当前构形下的随转坐标转换到总体坐标系下就可参与不平衡力的计算.

具体计算中,首先输入基本参数,确定时间轴(一般采用对数时间步划分[9]),然后进行时步循环,时步循环结束则徐变计算结束.其中在Δtn+1,即tntn+1时步内(计算tn+1时刻的平衡)包括5个主要步骤,按顺序依次为:

(1) 累计本时刻前(含本时刻)所有外荷载.

(2) 通过式(15)、(17),可求得本时步随转坐标系下收缩、徐变等效结点力Δfl, b, cn+1、Δfl, b, sn+1,进行累计$\sum\limits_{i = 1}^{n + 1} {(\Delta \mathit{\boldsymbol{f}}_{\rm{l},\rm{b},\rm{c}}^i + \Delta \mathit{\boldsymbol{f}}_{\rm{l},\rm{b},\rm{s}}^i)} $后即为不平衡力计算所需的收缩徐变等效节点力总量.

(3) 利用非线性求解器求解本时步结束时,即tn+1时刻总体坐标系下结点位移.中间过程在计算总体坐标系下切线刚度矩阵时涉及到的fl,须计入本时步前所有时步的收缩徐变等效节点力总量$\sum\limits_{i = 1}^{n + 1} {(\Delta \mathit{\boldsymbol{f}}_{\rm{l},\rm{b},\rm{c}}^i + \Delta \mathit{\boldsymbol{f}}_{\rm{l},\rm{b},\rm{s}}^i)} $不平衡力由总外荷载减去总体坐标系下收缩徐变等效节点力总量得到;可采用力收敛准则,范数类型为2范数,当||不平衡力||2/||外荷载||2 < 容许误差,则非线性迭代结束.

(4) 根据式(20),将本时段末总位移矢量pgn+1转换到随转坐标系下得到pln+1;与保存的上一时段随转坐标系下的总位移pln相减得到Δpln+1.显然,Δfl, b, en+1=Kl, bΔPln+1、Δfl, en+1=KlΔPln+1分别为本时段节点位移增量引起的对应于混凝土桥面板分部及单元整体的杆端力增量.

(5) 由结点位移增量引起的内力增量Δfl, b, en+1、Δfl, en+1分别与本时步收缩徐变等效结点力增量Δfl, b, sn+1fl, b, cn+1相加,可得本时步对应于混凝土桥面板分部及单元整体的杆端内力增量Δfl, b, en+1、Δfln+1,其中对应于混凝土桥面板分部的内力历史需要记录.

梁式结构几何非线性含P-Δ(轴力位移)及大位移效应.P-Δ效应的实质是考虑单元内力对其刚度的影响,大位移的实质是让单元在变形后的位置上发挥其作用,以满足平衡方程必须建立在结构变形后位置上这一重要条件.本文基于CR列式及总体坐标与局部随转坐标之间的转换关系,推导了考虑滑移的组合梁单元的几何非线性列式.其中转换矩阵Bij为位移的函数,Hij为位移及内力的函数.在计算BijHij时,代入零内力则不考虑P-Δ效应,代入零位移则转换矩阵Bij为与位移无关的常矩阵.在求解过程中,如代入零位移及零内力,可实现线性计算;如代入零位移及当前内力,可实现仅考虑P-Δ效应的二阶计算;如代入当前位移、内力,则可实现考虑P-Δ效应及大位移效应的全几何非线性计算.即收缩徐变的线性计算及非线性计算可通过统一算法实现,仅需在程序中进行相应的开关设置.

4 算例分析

文献[4]采用增量法对图 5所示钢混简支组合梁进行收缩徐变分析,考虑了界面滑移效应的影响,但未涉及几何非线性问题.其中简支组合梁钢梁上下翼缘宽0.335 m,厚10 mm,腹板高0.58 m,厚10 mm.混凝土弹性模量Ec28=3.362×104MPa,钢材弹性模量Es=2×105MPa,界面连接刚度为0.15 kN·mm-2.徐变、收缩系数按CEB-90模型[4]计算,其中参数年平均相对湿度RH=80%,依水泥种类而定的系数βsc=4,s=0.25,理论厚度h0=150 mm,混凝土28 d轴心抗压强度标准值fck=30 MPa.假定收缩开始时间为28 d,加载龄期为28 d,收缩徐变计算时间为25 550 d(约70年).文献组合梁作用均布荷载为P=25 kN·m-1,组合梁轴向集中力N=0,本文首先按该荷载进行计算,以与文献结果进行对比.然后在P=25 kN·m-1N=3 000 kN作用下,分别进行线性及非线性计算,以考察非线性影响.表 1给出梁端滑移s,跨中桥面板轴力Nc、桥面板弯矩Mc、钢梁轴力Ns、钢梁弯矩Ms、整体弯矩M、挠度f的短期与长期效应值.简支梁划分为12个单元,时间轴上按对数划分为6步.

图 5 简支组合梁 Fig.5 Simply supported composite beam
下载CSV 表 1 简支组合梁短期、长期计算结果 Tab.1 Calculation result of short term and long term effects of simply supported composite beam

表 1可见,对应N=0工况,本文计算结果与文献[4]基本相同,证明了本文考虑滑移的组合梁收缩徐变计算的正确性.对应N=3 000 kN工况,在线性计算情况下,对本文外部静定结构,长期、短期跨中总弯矩相同,长期、短期跨中挠度之比为1.90.在短期荷载作用下,几何非线性计算与线性计算得到的跨中挠度比值为1.16,跨中总弯矩比值为1.17.但对计入徐变后的长期效应,几何非线性和线性计算结果的比值变化较大,跨中挠度比值为1.19,跨中总体弯矩比值为1.32.同时还可看出几何非线性长期效应值并非几何线性短期效应值与几何非线性影响系数及徐变影响系数的简单乘积,表明几何非线性与徐变之间存在明显的耦合作用.另外,如按变形后受力平衡条件来计算跨中总弯矩,短期、长期值应分别为525.3、595.8 kN·m,与表 1列出本文的计算值525.3、595.6 kN·m最大相差0.3‰,说明本文计算结果具有很高的精度.

5 结论

(1) 基于界面滑移组合梁Newmark模型、收缩徐变初应变法、几何非线性随转列式全量法,得到界面滑移组合梁同步考虑几何非线性、收缩徐变的算法,为组合梁斜拉桥等大跨柔性结构精细分析打下了理论基础.

(2) 当界面连接刚度输入极大值时,可认为钢混界面间无滑移位移,此时对应于刚性连接的组合梁.

(3) 本文采用的收缩徐变逐步计算方法,具有易与有限元相结合的优点,通过算例得到了很好的验证.

(4) 本文采用的CR列式全量法,没有时步间误差累积的缺点,且求解精度与荷载步个数无关,大大提高了计算精度和效率.

(5) 本文算法中收缩徐变仅在随转坐标系下考虑,几何非线性通过转换矩阵考虑,具有几何非线性与收缩徐变算法上分离的优点.

(6) 本文收缩徐变的线性、非线性计算可通过统一算法实现,仅需程序执行中设置相应开关.

(7) 算例表明,对于考虑界面滑移的组合梁等柔性结构进行几何非线性和收缩徐变效应的耦合分析是很有必要的,收缩徐变线性计算结果可能会有较大偏差.

(8) 建议开展几何非线性界面滑移组合梁收缩徐变作用的试验研究,以进一步证明理论分析的正确性.

参考文献
[1] RANZI G, DALL'ASTA A, RAGNI L, et al. A geometric nonlinear model for composite beams with partial interaction[J]. Engineering Structures, 2010, 32(5): 1384
[2] BATTINI J, NGUYEN Q and HJIAJ M. Non-linear finite element analysis of composite beams with interlayer slips[J]. Computers and Structures, 2009, 87(13-14): 904
[3] DEZI B L, LANNI C and TARANTINO A M. Simplified creep analysis of composite beams with flexible connectors[J]. Journal of Structural Engineering, 1993, 119(5): 1484
[4] JURKIEWIEZ B, BUZON S and SIEFFERT J G. Incremental viscoelastic analysis of composite beams with partial interaction[J]. Computers and Structures, 2005, 83(21-22): 1780 DOI:10.1016/j.compstruc.2005.02.021
[5] 邓继华, 邵旭东, 彭建新. 几何非线性平面梁考虑收缩徐变的算法研究[J]. 湖南大学学报(自然科学版), 2014, 41(9): 14
DENG Jihua, SHAO Xudong and PENG Jianxin. Algorithm study of the geometrica1 nonlinearity plane beam considering creep and shrinkag[J]. Journal of Hunan University(Natura1 Sciences), 2014, 41(9): 14
[6] RANZI G, GARA F, LEONI G, et al. Analysis of composite beams with partial shear interaction using available modelling techniques:a comparative study[J]. Computers and Structures, 2006, 84(13-14): 930
[7] ASTA A D and ZONA A. Slip locking in finite elements for composite beams with deformable shear connection[J]. Finite Elements in Analysis and Design, 2004, 40(13-14): 1907
[8] 梁鹏. 超大跨度斜拉桥几何非线性及随机模拟分析[D]. 上海: 同济大学, 2004.
LIANG Peng. Geometrical nonlinearity and random simulation of super long span cable-stayed bridges[D]. Shanghai: Tongji University, 2004. http://d.wanfangdata.com.cn/Thesis/Y675230
[9] 颜东煌, 田仲初, 李学文, 等. 混凝土桥梁收缩徐变计算的有限元方法与应用[J]. 中国公路学报, 2004, 17(2): 55
YAN Donghuang, TIAN Zhongchu, LI Xuewen, et al. Finite element method and application for the shrinkage and creep of concrete bridges[J]. China Journal of Highway and Transport, 2004, 17(2): 55
[10] 匡文起, 张玉良, 辛克贵.结构矩阵分析和程序设计[M].第1版.北京:高等教育出版社, 1991.
KUANG Wenqi, ZHANG Yuliang, XIN Kegui. Structure matrix analysis and program design[M]. 1st ed. Beijing:Higher Education Press, 1991.