﻿ 组合梁几何非线性与长期效应的同步算法
 文章快速检索
 同济大学学报(自然科学版)  2017, Vol. 45 Issue (8): 1108-1113.  DOI: 10.11908/j.issn.0253-374x.2017.08.002 0

### 引用本文

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.

### 文章历史

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 随转坐标系下单元列式 1.1 单元位移

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

 图 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}$

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 平衡方程

 ${{\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)

 ${\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]$

 ${{\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)

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

 $\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)

 $\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)

1.3 收缩徐变逐步计算法

 ${\mathit{\boldsymbol{\varepsilon }}_{\rm{c}}}{\rm{ = }}{\mathit{\boldsymbol{r}}_{\rm{b}}}/{\mathit{\boldsymbol{D}}_{\rm{b}}}\varphi \left( {t,\tau } \right)$ (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)

 $\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)

 $\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)

 $\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{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)

 $\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 静力凝聚

 $\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)

 $\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 梁体运动描述

 $\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}}$

 $\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)

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

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

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

 ${\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)

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

 ${\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)

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

(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，其中对应于混凝土桥面板分部的内力历史需要记录.

4 算例分析

 图 5 简支组合梁 Fig.5 Simply supported composite beam

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.