Nonlinear Numerical Simulation of Finite Elements Based on Fiber Beam Elements With Shear Effects for Structures
-
摘要:
基于Euler-Bernoulli梁理论的经典纤维模型忽略了剪切变形给截面带来的影响，为了得到更加精确的梁单元模型，该文基于考虑剪切效应的纤维梁单元，根据Timoshenko梁理论，推导了该纤维梁单元的刚度矩阵，并结合弹塑性增量理论，同时考虑了几何非线性和材料非线性的双重影响，建立了压弯剪复杂应力状态下结构非线性有限元分析理论。该文最后利用MATLAB编制了相关计算程序，对钢筋混凝土和矩形钢管混凝土的典型压弯剪构件进行有限元数值模拟，得到了构件的荷载-位移非线性全过程曲线。典型算例的验证结果表明：该文建立的非线性有限元分析理论是通用、可行和正确的。
Abstract:The classical fiber beam model based on the Euler-Bernoulli beam theory ignores the influence of shear deformation on the section. To get a more accurate beam element model, based on the fiber beam element with shear effects and the Timoshenko beam theory, the stiffness matrix of the fiber beam element was deduced, and the dual effects of geometric nonlinearity and material nonlinearity were considered at the same time, combined with the elastoplastic incremental theory. Then, the nonlinear finite element analysis theory for the structure under the complex stress state of compression, bending and shear was established. Finally, a program was coded on MATLAB to conduct finite element numerical simulation of the typical compression-bending-shear members of reinforced concrete and rectangular concrete-filled steel tube, and the nonlinear full-process load-displacement curves were obtained. The analysis of the numerical examples show that, the established nonlinear finite element analysis theory is universal, feasible and correct.
-
Key words:
- fiber beam element /
- structural nonlinearity /
- shear effect /
- finite element analysis
-
引 言
基于梁单元的有限元分析模型是近年来人们关注的热点，使用不同梁单元理论的有限元模型对所分析结构的计算精度和计算效率有显著影响。早期的Euler-Bernoulli经典梁理论[1]应用广泛，在一般情况下可以得到比较满意的结果，但由于未考虑剪切变形给截面带来的影响，在梁高度和跨度相差不大的情况下，会有较大的误差。为了考虑剪切带来的影响，Timoshenko[2]于1921年提出了考虑剪切效应的修正梁理论，该理论仍保留了平截面假定，但实际上梁截面的剪应力与剪应变呈抛物线分布，是不均匀的，这与平截面假定相矛盾，于是引入了截面修正系数，将截面上不均匀的剪应力与剪应变等效为均匀分布，这样既能较好地考虑剪切的影响同时又简化了计算模型。Mari和Scordelis[3]提出了基于有限元刚度法的纤维梁模型，该模型基于Euler-Bernoulli梁理论，以单元积分点截面处的力学性能代表整个单元的力学性能，使用线性插值构造积分点的轴向位移场，用Hermite插值函数构造横向位移场。胡郑州等[4]在UL列式下，根据连续性介质力学和虚位移原理，并结合Timoshenko梁理论，推导出了考虑剪切效应的三维纤维梁单元模型。顾观东[5]在Euler-Bernoulli梁单元模型的基础上，修正了塑性阶段的平截面假定，并引入剪切位移，得到了更精确的梁单元模型。郑欣怡[6]将Timoshenko梁理论引入单轴本构的弯剪耦合纤维模型中，这种改进后的纤维梁单元模型考虑了弯剪耦合效应和剪切效应对轴向变形的间接影响。Thai等[7]利用纤维梁单元，对单向轴力和弯矩作用下的桥墩进行了双非线性有限元分析。张传超等[8]基于OPENSEES软件，对纤维梁单元在钢管混凝土柱应用方面的建模方法进行了研究。马银[9]基于纤维梁单元，对型钢混凝土结构进行建模分析，并进一步推导了考虑黏结滑移关系的型钢混凝土梁柱单元模型。蔺鹏臻等[10]基于ABAQUS软件， 采用精细化的纤维梁单元模型， 开发了钢筋和混凝土纤维梁单元材料用户子程序。由此可见， 纤维梁单元在理论研究方面和工程应用方面都得到了广泛的应用[11-13]。
目前，基于Euler-Bernoulli梁理论的纤维梁单元其纤维只能考虑单轴受力状态，无法考虑剪切效应；基于考虑剪切效应的Timoshenko梁理论，多以使用材料多维本构模型来考虑材料非线性，需计算的参数较多，且运算复杂。为此，本文基于考虑剪切效应的三维纤维梁单元有限元模型，推导了该纤维梁单元的截面刚度矩阵、单元弹性刚度矩阵和单元几何刚度矩阵，同时补充了文献[4]在推导插值函数矩阵时，矩阵元素
$ {N_{2,6}} $ 的遗漏，并根据UL列式法和弹塑性增量理论，考虑结构的几何非线性和材料非线性，旨在建立压弯剪复杂应力状态下结构非线性有限元分析理论，能在合理描述压弯剪作用的同时得到准确的计算结果。最后，使用MATLAB编制了相关计算程序，结合钢筋混凝土和矩形钢管混凝土的典型压弯剪构件进行有限元数值模拟，验证了本文所建立理论的有效性、正确性和通用性。1. 考虑剪切效应的纤维模型
1.1 基本假设
本文基于经典纤维模型，并结合Timoshenko梁理论，建立了考虑剪切效应的纤维模型。纤维梁单元模型如图1所示，并做如下假设：
1） 满足平截面假定，但截面不再与梁中轴线垂直，考虑剪切效应导致的截面翘曲；
2） 钢筋与混凝土充分黏结，无相对滑移，变形协调；
3） 剪应力与剪应变均匀分布。
1.2 刚度矩阵
1.2.1 截面刚度矩阵
对于截面内任一纤维单元，其空间几何变换矩阵为
$$ {{\boldsymbol{b}}_i} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0&{{z_i}}&{ - {y_i}} \\ 0&1&0&{ - {z_i}}&0&0 \\ 0&0&1&{{y_i}}&0&0 \end{array}} \right]{\text{.}} $$ (1) 与其材料性质相关的本构关系矩阵为
$$ {{\boldsymbol{D}}_i} = \left[ {\begin{array}{*{20}{c}} {{E_i}}&0&0 \\ 0&{{G_i}}&0 \\ 0&0&{{G_i}} \end{array}} \right], $$ (2) 其中，
$ {E_i} $ 和$ {G_i} $ 分别为第$ i $ 个纤维单元的弹性模量和剪切模量。通过积分运算，得到任一纤维单元的刚度：
$$ {\boldsymbol{k}}_i^e = \int_A {{\boldsymbol{b}}_i^{\rm{T}}{{\boldsymbol{D}}_i}{{\boldsymbol{b}}_i}{\text{d}}A} = {\boldsymbol{b}}_i^{\rm{T}}{{\boldsymbol{D}}_i}{{\boldsymbol{b}}_i}{A_i} ,$$ (3) 其中，
$ {A_i} $ 为纤维单元的截面面积。对截面内的全部纤维单元进行积分，得到截面刚度矩阵：
$$ {{\boldsymbol{k}}_{\rm{s}}} = \sum\limits_{i = 1}^n {{\boldsymbol{k}}_i^e} = \sum\limits_{i = 1}^n {\int_A {{\boldsymbol{b}}_i^{\rm{T}}{{\boldsymbol{D}}_i}{{\boldsymbol{b}}_i}{\text{d}}A} } = \sum\limits_{i = 1}^n {{\boldsymbol{b}}_i^{\rm{T}}{{\boldsymbol{D}}_i}{{\boldsymbol{b}}_i}{A_i}} {\text{.}}$$ (4) 将式(1)、(2)代入式(4)中，展开并进行积分运算，可以得到截面刚度矩阵的显式为
$$ {{\boldsymbol{k}}_{\rm{s}}} = \left[ {\begin{array}{*{20}{c}} {\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}} }&0&0&0&{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{z_i}} }&{ - \displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}} } \\ 0&{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }&0&{ - \displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}} }&0&0 \\ 0&0&{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }&{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}} }&0&0 \\ 0&{ - \displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}} }&{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}} }&{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} ( {y_i^2 + z_i^2} )}&0&0 \\ {\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{z_i}} }&0&0&0&{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2} }&{ - \displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} } \\ { - \displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}} }&0&0&0&{ - \displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }&{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2} } \end{array}} \right]{\text{.}} $$ (5) 1.2.2 弹性刚度矩阵
在得到截面刚度矩阵的基础上，可采用位移形函数将结构的位移与截面的变形联系起来，从而推导三维纤维梁单元的刚度矩阵。每个梁单元共有2个节点，每个节点有6个自由度，单元节点位移向量为
$$ {{{\boldsymbol{u}}^e}} = {( { {{u_1}}\;\;\;{{v_1}}\;\;\;{{w_1}}\;\;\;{{\theta _{1x}}}\;\;\;{{\theta _{1y}}}\;\;\;{{\theta _{1z}}}\;\;\;{{u_2}}\;\;\;{{v_2}}\;\;\;{{w_2}}\;\;\;{{\theta _{2x}}}\;\;\;{{\theta _{2y}}}\;\;\;{{\theta _{2z}}} } )^{\rm{T}}}，$$ (6) 其中，
$\left({u}_{1},{u}_{2}\right),\left({v}_{1},{v}_{2}\right),\left({w}_{1},{w}_{2}\right)$ 分别为节点在x轴、y轴和z轴方向的线位移，$ \left( {{\theta _{1x}},{\theta _{2x}}} \right) $ 为节点的扭转角位移，$({\theta }_{1y},{\theta }_{2y})和\left({\theta }_{1z},{\theta }_{2z}\right)$ 分别为节点在xOz和xOy平面内的转动角位移。单元轴线上任一点x处截面的位移用单元两端节点位移可表示为
$$ {\boldsymbol{u}}\left( x \right) = {\boldsymbol{N}}\left( x \right){{\boldsymbol{u}}^e} ,$$ (7) 式中，
$ {\boldsymbol{N}}\left( x \right) $ 为单元形函数矩阵。轴向位移及扭转角位移均采用相同的线性位移模式，横向弯曲位移均采用三次式的位移模式[14]，此时单元位移插值函数可写为$$ \left\{ \begin{aligned} & u\left( x \right) = {a_0} + {a_1}x, \\& v\left( x \right) = {b_0} + {b_1}x + {b_2}{x^2} + {b_3}{x^3}, \\& w\left( x \right) = {c_0} + {c_1}x + {c_2}{x^2} + {c_3}{x^3}, \\& {\theta _x} = {a_0} + {a_1}x, \\& {\theta _y} = {c_1} + 2{c_2}x + 3{c_3}{x^2} - {\gamma _{xz}}, \\& {\theta _z} = {b_1} + 2{b_2}x + 3{b_3}{x^2} - {\gamma _{xy}}, \end{aligned} \right. $$ (8) 其中，
$ {\gamma _{xy}} = - \dfrac{{6E{I_z}}}{{{\kappa _y}GA}}{b_3} $ ，$ {\gamma _{xz}} = - \dfrac{{6E{I_y}}}{{{\kappa _z}GA}}{c_3} $ ，$ {\kappa _y},{\kappa _z} $ 为截面剪切系数，$ {I_y} $ 和$ {I_z} $ 分别为绕$ y $ 轴和$ z $ 轴的截面惯性矩。单元形函数矩阵为
$$ {\boldsymbol{N}}\left( x \right) = \left[ {\begin{array}{*{20}{c}} {{N_{1,1}}}&0&0&0&0&0&{{N_{1,7}}}&0&0&0&0&0 \\ 0&{{N_{2,2}}}&0&0&0&{{N_{2,6}}}&0&{{N_{2,8}}}&0&0&0&{{N_{2,12}}} \\ 0&0&{{N_{3,3}}}&0&{{N_{3,5}}}&0&0&0&{{N_{3,9}}}&0&{{N_{3,11}}}&0 \\ 0&0&0&{{N_{4,4}}}&0&0&0&0&0&{{N_{4,10}}}&0&0 \\ 0&0&{{N_{5,3}}}&0&{{N_{5,5}}}&0&0&0&{{N_{5,9}}}&0&{{N_{5,11}}}&0 \\ 0&{{N_{6,2}}}&0&0&0&{{N_{6,6}}}&0&{{N_{6,8}}}&0&0&0&{{N_{6,12}}} \end{array}} \right], $$ (9) 其中
$$ {N_{1,1}} = 1 - \frac{x}{L} ,\; {N_{1,7}} = \frac{x}{L} ,\; {N_{2,2}} = \frac{{( {L - x} )( {{L^2} + Lx - 2{x^2} + 2{\varphi _z}} )}}{{L( {{L^2} + 2{\varphi _z}} )}} ,\; {N_{2,6}} = \frac{{x( {L - x} )( {{L^2} - xL + {\varphi _z}} )}}{{L( {{L^2} + 2{\varphi _z}} )}} , $$ $$ \begin{split} & {N_{2,8}} = \frac{{x( { - 2{x^2} + 3Lx + 2{\varphi _z}} )}}{{L( {{L^2} + 2{\varphi _z}} )}} ,\; {N_{2,12}} = - \frac{{x( {L - x} )( {Lx + {\varphi _z}} )}}{{L( {{L^2} + 2{\varphi _z}} )}} ,\; {N_{3,3}} = \frac{{( {L - x} )( {{L^2} + Lx - 2{x^2} + 2{\varphi _y}} )}}{{L( {{L^2} + 2{\varphi _y}} )}} ,\\& {N_{3,5}} = \frac{{x( {L - x} )( {{L^2} - xL + {\varphi _y}} )}}{{L( {{L^2} + 2{\varphi _y}} )}} ,\; {N_{3,9}} = \frac{{x( { - 2{x^2} + 3Lx + 2{\varphi _y}} )}}{{L( {{L^2} + 2{\varphi _y}} )}} ,\; {N_{3,11}} = - \frac{{x( {L - x} )( {Lx + {\varphi _y}} )}}{{L( {{L^2} + 2{\varphi _y}} )}} , \end{split}$$ $$ {N_{4,4}} = {N_{1,1}} ,\; {N_{4,10}} = {N_{1,7}} ,\; {N_{5,3}} = - \frac{{6x( {L - x} )}}{{L( {{L^2} + 2{\varphi _y}} )}} ,\; {N_{5,5}} = \frac{{( {L - x} )( {{L^2} - 3xL + 2{\varphi _y}} )}}{{L( {{L^2} + 2{\varphi _y}} )}} ,\; {N_{5,9}} = \frac{{6x( {L - x} )}}{{L( {{L^2} + 2{\varphi _y}} )}} , $$ $$ {N_{5,11}} = \frac{{x( { - 2{L^2} + 3xL + 2{\varphi _y}} )}}{{L( {{L^2} + 2{\varphi _y}} )}} ,\; {N_{6,2}} = - \frac{{6x\left( {L - x} \right)}}{{L\left( {{L^2} + 2{\varphi _z}} \right)}} ,\; {N_{6,6}} = \frac{{\left( {L - x} \right)( {{L^2} - 3xL + 2{\varphi _z}} )}}{{L\left( {{L^2} + 2{\varphi _z}} \right)}} , $$ $$ {N_{6,8}} = \frac{{6x( {L - x} )}}{{L( {{L^2} + 2{\varphi _z}} )}} ,\; {N_{6,12}} = \frac{{x( { - 2{L^2} + 3xL + 2{\varphi _z}} )}}{{L( {{L^2} + 2{\varphi _z}} )}} ,\; {\varphi _y} = \frac{{6E{I_z}}}{{{\kappa _y}GA}} ,\; {\varphi _z} = \frac{{6E{I_y}}}{{{\kappa _z}GA}} {\text{.}} $$ 对单元形函数矩阵进一步推导，可得到线性应变矩阵
$ {{\boldsymbol{B}}_{\rm{L}}} $ ，再由虚功原理和变分原理得单元弹性刚度矩阵为$$ {\boldsymbol{k}}_{\rm{L}}^e = \int_l {{\boldsymbol{B}}_{\rm{L}}^{\rm{T}}{{\boldsymbol{k}}_{\rm{s}}}{{\boldsymbol{B}}_{\rm{L}}}{\text{d}}x} {\text{.}}$$ (10) 做积分运算后，得到单元弹性刚度矩阵的显式为
$$ {\boldsymbol{k}}_{\rm{L}}^e = {\left[ \begin{gathered} {{\boldsymbol{k}}_{11}}{\text{ }}{{\boldsymbol{k}}_{12}}{\text{ }}{{\boldsymbol{k}}_{13}}{\text{ }}{{\boldsymbol{k}}_{14}} \hfill \\ {{\boldsymbol{k}}_{21}}{\text{ }}{{\boldsymbol{k}}_{22}}{\text{ }}{{\boldsymbol{k}}_{23}}{\text{ }}{{\boldsymbol{k}}_{24}} \hfill \\ {{\boldsymbol{k}}_{31}}{\text{ }}{{\boldsymbol{k}}_{32}}{\text{ }}{{\boldsymbol{k}}_{33}}{\text{ }}{{\boldsymbol{k}}_{34}} \hfill \\ {{\boldsymbol{k}}_{41{\text{ }}}}{\text{ }}{{\boldsymbol{k}}_{42}}{\text{ }}{{\boldsymbol{k}}_{43}}{\text{ }}{{\boldsymbol{k}}_{44}} \end{gathered} \right]_{12 \times 12}} ,$$ (11) 其中，单元弹性刚度矩阵为对称矩阵，即
$$ {{\boldsymbol{k}}_{21}} = {\boldsymbol{k}}_{12}^{\rm{T}},\;{{\boldsymbol{k}}_{31}} = {\boldsymbol{k}}_{13}^{\rm{T}},\;{{\boldsymbol{k}}_{41}} = {\boldsymbol{k}}_{14}^{\rm{T}},\;{{\boldsymbol{k}}_{32}} = {\boldsymbol{k}}_{23}^{\rm{T}},\;{{\boldsymbol{k}}_{42}} = {\boldsymbol{k}}_{24}^{\rm{T}},\;{{\boldsymbol{k}}_{43}} = {\boldsymbol{k}}_{34}^{\rm{T}},\; $$ $$ {{\boldsymbol{k}}_{11}} = \left[ {\begin{array}{*{20}{c}} \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}} }}{L}& 0& 0 \\ 0 & \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}( {\varPhi _Z^2 - 2{\varPhi _Z} + 1} )} {L^2} + 12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2} }}{{\varPhi _Z^2{L^3}}}&\dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^3}}} \\ 0 & \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^3}}}&\dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2} }}{{\varPhi _Y^2{L^3}}} + \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{{\left( {{\varPhi _Z} - 1} \right)}^2}} }}{{\varPhi _Z^2L}} \end{array}} \right] ,$$ $$ {{\boldsymbol{k}}_{12}} = \left[ {\begin{array}{*{20}{c}} 0&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{z_i}} }}{L}& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}} }}{L} \\[-2pt] \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}\left( {1 - {\varPhi _Z}} \right)} }}{{{\varPhi _Z}L}}& - \dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}}&\begin{array}{c}\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }}{2} +\\[-2pt] \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2} + {L^2}\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}\left( {1 - 2{\varPhi _Z}} \right)} }}{{2\varPhi _Z^2{L^2}}} \end{array}\\[-2pt] \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}\left( {1 - {\varPhi _Z}} \right)} }}{{{\varPhi _Z}L}}&\begin{array}{c}\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} \left( {1 - {\varPhi _Z}} \right)}}{{2{\varPhi _Z}}} -\\[-2pt] \dfrac{{12{\varPhi _Z}\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2} + {\varPhi _Y}\left( {1 - {\varPhi _Z}} \right){L^2}\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }}{{2\varPhi _Y^2{\varPhi _Z}{L^2}}}\end{array}&\dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}} \end{array}} \right], $$ $$ {{\boldsymbol{k}}_{13}} = \left[ {\begin{array}{*{20}{c}} - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}} }}{L}&0&0 \\ 0& - \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2 + \displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{L^2}{{\left( {{\varPhi _Z} - 1} \right)}^2}} } }}{{\varPhi _Z^2{L^3}}}& - \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^3}}} \\ 0& - \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^3}}}& - \dfrac{{12{\varPhi _Z}\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2 + \displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{L^2}{\varPhi _Y}\left( {1 - {\varPhi _Y} - {\varPhi _Z} + {\varPhi _Y}{\varPhi _Z}} \right)} } }}{{\varPhi _Y^2{\varPhi _Z}{L^3}}} \end{array}} \right], $$ $$ {{\boldsymbol{k}}_{14}} = \left[ {\begin{array}{*{20}{c}} 0& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{z_i}} }}{L}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}} }}{L} \hfill \\[-1pt] \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}\left( {{\varPhi _Z} - 1} \right)} }}{{{\varPhi _Z}L}}& - \dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}}&\begin{array}{c}\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }}{2} +\\ \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2 + \displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} {L^2}\left( {1 - {\varPhi _Z}} \right)} }}{{2\varPhi _Z^2{L^2}}}\end{array} \hfill \\[-1pt] \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}\left( {{\varPhi _Z} - 1} \right)} }}{{{\varPhi _Z}L}}&\begin{array}{c}\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} \left( {1 - {\varPhi _Z}} \right)}}{{2{\varPhi _Z}}} -\\[-1pt] \dfrac{{12{\varPhi _Z}\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2 + {\varPhi _Y}{L^2}\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}\left( {1 - {\varPhi _Z}} \right)} } }}{{2\varPhi _Y^2{\varPhi _Z}{L^2}}}\end{array}&\dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}} \end{array}} \right] ,$$ $$ {{\boldsymbol{k}}_{22}} = \left[ {\begin{array}{*{20}{c}} \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}( {y_i^2 + z_i^2} )} }}{L}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}( {{\varPhi _Y} - 1} )} }}{{2{\varPhi _Y}}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}( {1 - {\varPhi _Z}} )} }}{{2{\varPhi _Z}}} \hfill \\ \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}( {{\varPhi _Y} - 1} )} }}{{2{\varPhi _Y}}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2( {\varPhi _Y^2 + 3} )} }}{{\varPhi _Y^2L}} + \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} L{{( {{\varPhi _Y} - 1} )}^2}}}{{4\varPhi _Y^2}}& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} ( {{\varPhi _Y}{\varPhi _Z} + 3} )}}{{{\varPhi _Y}{\varPhi _Z}L}} \hfill \\ \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}( {1 - {\varPhi _Z}} )} }}{{2{\varPhi _Z}}}& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} ( {{\varPhi _Y}{\varPhi _Z} + 3} )}}{{{\varPhi _Y}{\varPhi _Z}L}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2( {\varPhi _Z^2 + 3} )} }}{{\varPhi _Z^2L}} + \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} L{{( {{\varPhi _Z} - 1} )}^2}}}{{4\varPhi _Z^2}}\end{array}} \right], $$ $$ $$ $$ {{\boldsymbol{k}}_{23}} = \left[ {\begin{array}{*{20}{c}} 0&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}\left( {{\varPhi _Z} - 1} \right)} }}{{{\varPhi _Z}L}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}\left( {{\varPhi _Y} - 1} \right)} }}{{{\varPhi _Y}L}} \hfill \\ - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{z_i}} }}{L}&\dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }}{2} + \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2 + \displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} {L^2}\left( {1 - {\varPhi _Y}} \right)} }}{{2\varPhi _Y^2{L^2}}} \hfill \\ \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}} }}{L}&\begin{array}{c} - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }}{2} - \\\dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2 + \displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} {L^2}\left( {1 - {\varPhi _Z}} \right)} }}{{2\varPhi _Z^2{L^2}}}\end{array}& - \dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}} \end{array}} \right], $$ $$ {{\boldsymbol{k}}_{24}} = \left[ {\begin{array}{*{20}{c}} - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}( {y_i^2 + z_i^2} )} }}{L}& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}} ( {1 - {\varPhi _Y}} )}}{{2{\varPhi _Y}}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}} ( {1 - {\varPhi _Z}} )}}{{2{\varPhi _Z}}} \hfill \\ \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}} ( {1 - {\varPhi _Y}} )}}{{2{\varPhi _Y}}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2( {3 - \varPhi _Y^2} )} }}{{\varPhi _Y^2L}} + \dfrac{{L\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} {{( {1 - {\varPhi _Y}} )}^2}}}{{4\varPhi _Y^2}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} ( {{\varPhi _Y}{\varPhi _Z} - 3} )}}{{{\varPhi _Y}{\varPhi _Z}L}} \hfill \\ - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}} ( {1 - {\varPhi _Z}} )}}{{2{\varPhi _Z}}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} ( {{\varPhi _Y}{\varPhi _Z} - 3} )}}{{{\varPhi _Y}{\varPhi _Z}L}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2( {3 - \varPhi _Z^2} )} }}{{\varPhi _Z^2L}} + \dfrac{{L\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} {{( {1 - {\varPhi _Z}} )}^2}}}{{4\varPhi _Z^2}} \end{array}} \right], $$ $$ {{\boldsymbol{k}}_{33}} = \left[{\begin{array}{*{20}{c}} \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}} }}{L}&0&0 \hfill \\ 0&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}( {\varPhi _Z^2 - 2{\varPhi _Z} + 1} )} {L^2} + 12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2} }}{{\varPhi _Z^2{L^3}}}&\dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^3}}} \hfill \\ 0&\dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^3}}}&\dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2} }}{{\varPhi _Y^2{L^3}}} + \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{{( {{\varPhi _Y} - 1} )}^2}} }}{{\varPhi _Y^2L}} \end{array}} \right] ,$$ $$ {{\boldsymbol{k}}_{34}} = \left[{\begin{array}{*{20}{c}} 0&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{z_i}} }}{L}& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}} }}{L} \hfill \\ \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}\left( {1 - {\varPhi _Z}} \right)} }}{{{\varPhi _Z}L}}&\dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}}& \begin{array}{c}- \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }}{2} - \\\dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2} + {L^2}\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}\left( {1 - 2{\varPhi _Z}} \right)} }}{{2\varPhi _Z^2{L^2}}} \end{array}\hfill \\ \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}\left( {1 - {\varPhi _Y}} \right)} }}{{{\varPhi _Y}L}}&\begin{array}{c}\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} }}{2} +\\ \dfrac{{12\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2} + {L^2}\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}\left( {1 - 2{\varPhi _Y}} \right)} }}{{2\varPhi _Y^2{L^2}}}\end{array}& - \dfrac{{6\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} }}{{{\varPhi _Y}{\varPhi _Z}{L^2}}} \end{array}} \right], $$ $$ {{\boldsymbol{k}}_{44}} = \left[ {\begin{array}{*{20}{c}} \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}( {y_i^2 + z_i^2} )} }}{L}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}( {1 - {\varPhi _Y}} )} }}{{2{\varPhi _Y}}}& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}( {1 - {\varPhi _Z}} )} }}{{2{\varPhi _Z}}} \hfill \\ \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{y_i}( {1 - {\varPhi _Y}} )} }}{{2{\varPhi _Y}}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}z_i^2( {\varPhi _Y^2 + 3} )} }}{{\varPhi _Y^2L}} + \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} L{{( {{\varPhi _Y} - 1} )}^2}}}{{4\varPhi _Y^2}}&- \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} ( {{\varPhi _Y}{\varPhi _Z} + 3} )}}{{{\varPhi _Y}{\varPhi _Z}L}} \hfill \\ - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}{z_i}( {1 - {\varPhi _Z}} )} }}{{2{\varPhi _Z}}}& - \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}{y_i}{z_i}} ( {{\varPhi _Y}{\varPhi _Z} + 3} )}}{{{\varPhi _Y}{\varPhi _Z}L}}&\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{E_i}{A_i}y_i^2( {\varPhi _Z^2 + 3} )} }}{{\varPhi _Z^2L}} + \dfrac{{\displaystyle\sum\limits_{i = 1}^n {{G_i}{A_i}} L{{( {{\varPhi _Z} - 1} )}^2}}}{{4\varPhi _Z^2}} \end{array}} \right] ,$$ $ {\varPhi _Y} = 1 + \dfrac{{12E{I_y}}}{{{\kappa _z}GA{L^2}}} $ ，$ {\varPhi _Z} = 1 + \dfrac{{12E{I_z}}}{{{\kappa _y}GA{L^2}}} $ ，L为单元长度。1.2.3 几何刚度矩阵
在两节点12个自由度单元线弹性分析力学模型的基础上，本文将建立单元的几何非线性分析力学模型，并推导出单元的几何刚度矩阵。
非线性有限元方程为
$$ \int\nolimits_V {{\boldsymbol{B}}_{\rm{L}}^{\rm{T}}{\boldsymbol{D}}{{\boldsymbol{B}}_{\rm{L}}}{\text{d}}v} \Delta {{\boldsymbol{u}}^e} + \int\nolimits_V {{\boldsymbol{B}}_{{\rm{NL}}}^{\rm{T}}{\boldsymbol{\sigma }}{\text{d}}v} = \Delta {\boldsymbol{P}} - \int\nolimits_V {{\boldsymbol{B}}_{\rm{L}}^{\rm{T}}{\boldsymbol{\sigma }}{\text{d}}v}, $$ (12) 其中，
$ {{\boldsymbol{B}}_{{\rm{NL}}}} $ 为非线性应变矩阵，$ \Delta {\boldsymbol{P}} $ 为等效荷载。在式(12)中，左边第一项积分为弹性刚度矩阵，通过第二项积分，可推导几何刚度矩阵。
$$ \int\nolimits_V {{\boldsymbol{B}}_{\rm{NL}}^{\rm{T}}{\boldsymbol{\sigma }}{\text{d}}v = } \int\nolimits_V {{{\left( {{\boldsymbol{AG}}} \right)}^{\rm{T}}}{\boldsymbol{\sigma }}{\text{d}}v = } \int\nolimits_V {{{\boldsymbol{G}}^{\rm{T}}}{{\boldsymbol{A}}^{\rm{T}}}{\boldsymbol{\sigma }}{\text{d}}v}, $$ (13) 其中，
$ {{\boldsymbol{A}}^{\rm{T}}}{\boldsymbol{\sigma }} = {\boldsymbol{M\theta }} = {\boldsymbol{MG}}\Delta {{\boldsymbol{u}}^e} $ ，$ {\boldsymbol{A}} $ 为非线性应变与位移梯度矢量之间的转换矩阵，$ {\boldsymbol{M}} $ 为应力分量矩阵，$ {\boldsymbol{\theta }} $ 为位移梯度矢量，$ {\boldsymbol{G}} $ 为位移梯度矢量与单元节点位移矢量之间的转换矩阵，故有$$ \int\nolimits_V {{\boldsymbol{B}}_{\rm{NL}}^{\rm{T}}{\boldsymbol{\sigma }}{\text{d}}v = } \int\nolimits_V {{{\boldsymbol{G}}^{\rm{T}}}{\boldsymbol{MG}}\Delta {{\boldsymbol{u}}^e}{\text{d}}v = } \int\nolimits_V {{{\boldsymbol{G}}^{\rm{T}}}{\boldsymbol{MG}}{\text{d}}v \cdot \Delta {{\boldsymbol{u}}^e}} {\text{.}}$$ (14) 因此，单元几何刚度矩阵为
$$ {\boldsymbol{k}}_{\rm{G}}^e = \int\nolimits_V {{{\boldsymbol{G}}^{\rm{T}}}{\boldsymbol{MG}}{\text{d}}v = \int\nolimits_l {{{\boldsymbol{G}}^{\rm{T}}}{\boldsymbol{\hat MG}}{\text{d}}x} }, $$ (15) 式中
${\boldsymbol{G}}和\hat{{\boldsymbol{M}}}$ 的显式为[4]$$ {\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} {{G_{1,1}}}&0&0&0&0&0&{{G_{1,7}}}&0&0&0&0&0 \\ 0&{{G_{2,2}}}&0&0&0&{{G_{2,6}}}&0&{{G_{2,8}}}&0&0&0&{{G_{2,12}}} \\ 0&0&{{G_{3,3}}}&0&{{G_{3,5}}}&0&0&0&{{G_{3,9}}}&0&{{G_{3,11}}}&0 \\ 0&0&0&{{G_{4,4}}}&0&0&0&0&0&{{G_{4,10}}}&0&0 \\ 0&0&{{G_{5,3}}}&0&{{G_{5,5}}}&0&0&0&{{G_{5,9}}}&0&{{G_{5,11}}}&0 \\ 0&{{G_{6,2}}}&0&0&0&{{G_{6,6}}}&0&{{G_{6,8}}}&0&0&0&{{G_{6,12}}} \\ 0&0&0&{{G_{7,4}}}&0&0&0&0&0&{{G_{7,10}}}&0&0 \\ 0&0&{{G_{8,3}}}&0&{{G_{8,5}}}&0&0&0&{{G_{8,9}}}&0&{{G_{8,11}}}&0 \\ 0&{{G_{9,2}}}&0&0&0&{{G_{9,6}}}&0&{{G_{9,8}}}&0&0&0&{{G_{9,12}}} \end{array}} \right], $$ $$ {G_{1,1}} = - \frac{1}{L} ,\; {G_{1,7}} = \frac{1}{L}, $$ $$ {G_{2,2}} = \frac{{\left( {1 - {\varPhi _Z}} \right){L^2} - 6x\left( {L - x} \right)}}{{{\varPhi _Z}{L^3}}} ,\; {G_{2,6}} = \frac{{ - 2{\varPhi _Z}Lx - 6x\left( {L - x} \right) + \left( {1 + {\varPhi _z}} \right){L^2}}}{{2{\varPhi _Z}{L^2}}} ,$$ $$ {G_{2,8}} = - {G_{2,2}} ,\; {G_{2,12}} = \frac{{2{\varPhi _Z}Lx - 6x\left( {L - x} \right) - \left( {{\varPhi _z} - 1} \right){L^2}}}{{2{\varPhi _Z}{L^2}}}, $$ $$ {G_{3,3}} = \frac{{\left( {1 - {\varPhi _Y}} \right){L^2} + 6x\left( {x - L} \right)}}{{{\varPhi _Y}{L^3}}} ,\; {G_{3,5}} = \frac{{2{\varPhi _Y}Lx - 6x\left( {x - L} \right) - \left( {1 + {\varPhi _Y}} \right){L^2}}}{{2{\varPhi _Y}{L^2}}}, $$ $$ {G_{3,9}} = - {G_{3,3}} ,\; {G_{3,11}} = \frac{{ - 2{\varPhi _Y}Lx - 6x\left( {x - L} \right) - \left( {1 - {\varPhi _Y}} \right){L^2}}}{{2{\varPhi _Y}{L^2}}}, $$ $$ {G_{4,4}} = 1 - \frac{x}{L} ,\; {G_{4,10}} = \frac{x}{L} ,$$ $$ {G_{5,3}} = \frac{{6x\left( {L - x} \right)}}{{{\varPhi _Y}{L^3}}} ,\; {G_{5,5}} = \frac{{\left( {3x - {\varPhi _Y}L} \right)\left( {x - L} \right)}}{{{\varPhi _Y}{L^2}}} ,\; {G_{5,9}} = - {G_{5,3}} ,\; {G_{5,11}} = \frac{{\left( {3x - 3L + {\varPhi _Y}L} \right)x}}{{{\varPhi _Y}{L^2}}}, $$ $$ {G_{6,2}} = - \frac{{6x\left( {L - x} \right)}}{{{\varPhi _Z}{L^3}}} ,\; {G_{6,6}} = \frac{{\left( {3x - {\varPhi _Z}L} \right)\left( {x - L} \right)}}{{{\varPhi _Z}{L^2}}} ,\; {G_{6,8}} = - {G_{6,2}} ,\; {G_{6,12}} = \frac{{\left( {3x - 3L + {\varPhi _Z}L} \right)x}}{{{\varPhi _Z}{L^2}}} ,$$ $$ {G_{7,4}} = {G_{1,1}} ,\; {G_{7,10}} = {G_{1,7}}, $$ $$ {G_{8,3}} = \frac{{6\left( {L - 2x} \right)}}{{{\varPhi _Y}{L^3}}} ,\; {G_{8,5}} = \frac{{6x - 3L - {\varPhi _Y}L}}{{{\varPhi _Y}{L^2}}} ,\; {G_{8,9}} = - {G_{8,3}} ,\; {G_{8,11}} = \frac{{6x - 3L + {\varPhi _Y}L}}{{{\varPhi _Y}{L^2}}}, $$ $$ {G_{9,2}} = - \frac{{6\left( {L - 2x} \right)}}{{{\varPhi _Z}{L^3}}} ,\; {G_{9,6}} = \frac{{6x - 3L - {\varPhi _Z}L}}{{{\varPhi _Z}{L^2}}} ,\; {G_{9,8}} = - {G_{9,2}} ,\; {G_{9,12}} = \frac{{6x - 3L + {\varPhi _Z}L}}{{{\varPhi _Z}{L^2}}} ,$$ $$ {\boldsymbol{\hat M}} = \left[ {\begin{array}{*{20}{c}} {{M_{1,1}}}&0&0&0&{{M_{1,5}}}&{{M_{1,6}}}&0&{{M_{1,8}}}&{{M_{1,9}}} \\ 0&{{M_{2,2}}}&0&{{M_{2,4}}}&0&0&{{M_{2,7}}}&0&0 \\ 0&0&{{M_{3,3}}}&{{M_{3,4}}}&0&0&{{M_{3,7}}}&0&0 \\ 0&{{M_{4,2}}}&{{M_{4,3}}}&0&0&0&{{M_{4,7}}}&0&0 \\ {{M_{5,1}}}&0&0&0&0&{{M_{5,6}}}&0&{{M_{5,8}}}&{{M_{5,9}}} \\ {{M_{6,1}}}&0&0&0&{{M_{6,5}}}&0&0&{{M_{6,8}}}&{{M_{6,9}}} \\ 0&{{M_{7,2}}}&{{M_{7,3}}}&{{M_{7,4}}}&0&0&{{M_{7,7}}}&0&0 \\ {{M_{8,1}}}&0&0&0&{{M_{8,5}}}&{{M_{8,6}}}&0&{{M_{8,8}}}&{{M_{8,9}}} \\ {{M_{9,1}}}&0&0&0&{{M_{9,5}}}&{{M_{9,6}}}&0&{{M_{9,8}}}&{{M_{9,9}}} \end{array}} \right], $$ $$ {M_{1,1}} = \sum\limits_{i = 1}^n {{\sigma _{11}}{A_i}} ,\; {M_{1,5}} = \sum\limits_{i = 1}^n {{\tau _{xz}}{A_i}} ,\; {M_{1,6}} = - \sum\limits_{i = 1}^n {{\tau _{xy}}{A_i}} ,\; {M_{1,8}} = \sum\limits_{i = 1}^n {{\sigma _{11}}{z_i}{A_i}} ,\; {M_{1,9}} = - \sum\limits_{i = 1}^n {{\sigma _{11}}{y_i}{A_i}}, $$ $$ {M_{2,2}} = \sum\limits_{i = 1}^n {{\sigma _{11}}{A_i}} ,\; {M_{2,4}} = - \sum\limits_{i = 1}^n {{\tau _{xz}}{A_i}} ,\; {M_{2,7}} = - \sum\limits_{i = 1}^n {{\sigma _{11}}{z_i}{A_i}}, $$ $$ {M_{3,3}} = \sum\limits_{i = 1}^n {{\sigma _{11}}{A_i}} ,\; {M_{3,4}} = \sum\limits_{i = 1}^n {{\tau _{xy}}{A_i}} ,\; {M_{3,7}} = \sum\limits_{i = 1}^n {{\sigma _{11}}{y_i}{A_i}}, $$ $$ {M_{4,2}} = - \sum\limits_{i = 1}^n {{\tau _{xz}}{A_i}} ,\; {M_{4,3}} = \sum\limits_{i = 1}^n {{\tau _{xy}}{A_i}} ,\; {M_{4,7}} = \sum\limits_{i = 1}^n {( {{\tau _{xz}}{z_i} + {\tau _{xy}}{y_i}} ){A_i}}, $$ $$ {M_{5,1}} = \sum\limits_{i = 1}^n {{\tau _{xz}}{A_i}} ,\; {M_{5,6}} = - \sum\limits_{i = 1}^n {{\tau _{yz}}{A_i}} ,\; {M_{5,8}} = \sum\limits_{i = 1}^n {{\tau _{xz}}{z_i}{A_i}} ,\; {M_{5,9}} = -\sum\limits_{i = 1}^n {{\tau _{xz}}{y_i}{A_i}} ,$$ $$ {M_{6,1}} = - \sum\limits_{i = 1}^n {{\tau _{xy}}{A_i}} ,\; {M_{6,5}} = - \sum\limits_{i = 1}^n {{\tau _{yz}}{A_i}} ,\; {M_{6,8}} = - \sum\limits_{i = 1}^n {{\tau _{xy}}{z_i}{A_i}} ,\; {M_{6,9}} = \sum\limits_{i = 1}^n {{\tau _{xy}}{y_i}{A_i}}, $$ $$ {M_{7,2}} = - \sum\limits_{i = 1}^n {{\sigma _{11}}{z_i}{A_i}} ,\; {M_{7,3}} = \sum\limits_{i = 1}^n {{\sigma _{11}}{y_i}{A_i}} ,\; {M_{7,4}} = \sum\limits_{i = 1}^n {( {{\tau _{xz}}{z_i} + {\tau _{xy}}{y_i}} ){A_i}} ,\; {M_{7,7}} = \sum\limits_{i = 1}^n {( {z_i^2 + y_i^2} ){\sigma _{11}}{A_i}} ,$$ $$ {M_{8,1}} = \sum\limits_{i = 1}^n {{\sigma _{11}}{z_i}{A_i}} ,\; {M_{8,5}} = \sum\limits_{i = 1}^n {{\tau _{xz}}{z_i}{A_i}} ,\; {M_{8,6}} = - \sum\limits_{i = 1}^n {{\tau _{xy}}{z_i}{A_i}} ,\; {M_{8,8}} = \sum\limits_{i = 1}^n {{\sigma _{11}}z_i^2{A_i}} ,\; {M_{8,9}} = - \sum\limits_{i = 1}^n {{\sigma _{11}}{z_i}{y_i}{A_i}}, $$ $$ {M_{9,1}} = - \sum\limits_{i = 1}^n {{\sigma _{11}}{y_i}{A_i}} ,\; {M_{9,5}} = - \sum\limits_{i = 1}^n {{\tau _{xz}}{y_i}{A_i}} ,\; {M_{9,6}} = \sum\limits_{i = 1}^n {{\tau _{xy}}{y_i}{A_i}} ,\; {M_{9,8}} = - \sum\limits_{i = 1}^n {{\sigma _{11}}{z_i}{y_i}{A_i}} ,\; {M_{9,9}} = \sum\limits_{i = 1}^n {{\sigma _{11}}y_i^2{A_i}}, $$ 其中，
$ {\sigma _{11}} $ 为沿x轴的轴向应力，$ {\tau _{xy}} $ 为垂直x轴沿y轴的剪应力，$ {\tau _{xz}} $ 为垂直x轴沿z轴的剪应力，$ {\tau _{yz}} $ 为垂直y轴沿z轴的剪应力。通过式(15)的积分运算可得单元几何刚度矩阵，最终，单元刚度矩阵为弹性刚度矩阵和几何刚度矩阵之和，即
$$ {{\boldsymbol{k}}^e} = {\boldsymbol{k}}_{\rm{L}}^e + {\boldsymbol{k}}_{\rm{G}}^e{\text{.}} $$ (16) 2. 复杂应力状态下纤维单元应力的更新
2.1 等效应力与等效应变
本文将采用Mises屈服准则判断复杂应力状态下纤维单元的屈服状态，并据此理论求等效应力与等效应变。式(17)和(18)为等效应力与等效应变的分量表示方法：
$$ \bar \sigma = \frac{{\sqrt 2 }}{2}\sqrt {{{( {{\sigma _x} - {\sigma _y}} )}^2} + {{( {{\sigma _y} - {\sigma _z}} )}^2} + {{( {{\sigma _z} - {\sigma _x}} )}^2} + 6( {\tau _{xy}^2 + \tau _{yz}^2 + \tau _{zx}^2} )}, $$ (17) $$ \bar \varepsilon = \frac{{\sqrt 2 }}{3}\sqrt {{{( {{\varepsilon _x} - {\varepsilon _y}} )}^2} + {{( {{\varepsilon _y} - {\varepsilon _z}} )}^2} + {{( {{\varepsilon _z} - {\varepsilon _x}} )}^2} + 6( {\gamma _{xy}^2 + \gamma _{yz}^2 + \gamma _{zx}^2} )} {\text{.}}$$ (18) 等效应力与等效应变存在以下关系：
$$ \frac{{{\text{d}}\phi \left( {\bar \varepsilon } \right)}}{{{\text{d}}\bar \varepsilon }} = \frac{{{\text{d}}\bar \sigma }}{{{\text{d}}\bar \varepsilon }} ,$$ (19) 其中，
$ \phi $ 为已知的材料本构关系。2.2 基于增量理论应力的更新
由于塑性本构关系中的应力和应变不存在一一对应的关系，一般只能建立应力与应变增量之间的关系，以增量形式来表示。依据此理论对纤维单元应力状态的更新步骤如下：
1） 根据已知的应力分量
${\boldsymbol{\sigma }}$ 和应变分量$ {\boldsymbol{\varepsilon }} $ 计算出等效应力$ \bar \sigma $ 和等效应变$ \bar \varepsilon $ ，计算应变分量增量$ {\text{d}} {\boldsymbol{\varepsilon }} $ 。2） 在Mises屈服准则下，通过比较等效应力与纤维单元的屈服应力判断是否屈服，若未屈服，按第3）步弹性理论更新应力；若屈服，按第4）步塑性理论更新应力。
3） 根据弹性理论，计算弹性Jacobi矩阵
$ {{\boldsymbol{D}}_{\rm{e}}} $ ：$$ {{\boldsymbol{D}}}_{{\rm{e}}}=\dfrac{E\left(1-\mu \right)}{\left(1+\mu \right)\left(1-2\mu \right)}\left[\begin{array}{cccccc}1& & & & & \\ \dfrac{\mu }{1-\mu }& 1& & & {\rm{Sys}}& \\ \dfrac{\mu }{1-\mu }& \dfrac{\mu }{1-\mu }& 1& & & \\ 0& 0& 0& \dfrac{1-2\mu }{2\left(1-\mu \right)}& & \\ 0& 0& 0& 0& \dfrac{1-2\mu }{2\left(1-\mu \right)}& \\ 0& 0& 0& 0& 0& \dfrac{1-2\mu }{2\left(1-\mu \right)}\end{array}\right], $$ (20) 其中，
$ E $ 为纤维单元的弹性模量，$ \mu $ 为Poisson比。然后，计算应力分量增量
$ {\text{d}} {\boldsymbol{\sigma }} $ ：$$ {\text{d}} {\boldsymbol{\sigma }} = {{\boldsymbol{D}}_{\rm{e}}}{\text{d}} {\boldsymbol{\varepsilon }}{\text{.}} $$ (21) 最后，按第5）步更新应力。
4） 根据塑性理论，计算塑性Jacobi矩阵
$ {{\boldsymbol{D}}_{{\rm{ep}}}} $ ：$$ {{\boldsymbol{D}}}_{{\rm{ep}}}=\dfrac{E}{1+\mu }\left[\begin{array}{cccccc}\dfrac{1-\mu }{1-2\mu }-\omega {{\sigma }}_{x}^{{\prime }2}& & & & & \\ \dfrac{\mu }{1-2\mu }-\omega {{\sigma }}_{x}^{\prime }{\sigma }_{y}& \dfrac{1-\mu }{1-2\mu }-\omega {{\sigma }}_{y}^{{\prime }2}& & & {\rm{Sym}}& \\ \dfrac{1-\mu }{1-2\mu }-\omega {{\sigma }}_{x}^{\prime }{\sigma }_{z}& \dfrac{\mu }{1-2\mu }-\omega {{\sigma }}_{y}^{\prime }{{\sigma }}_{z}^{\prime }& \dfrac{1-\mu }{1-2\mu }-\omega {{\sigma }}_{z}^{{\prime }2}& & & \\ -\omega {{\sigma }}_{z}^{\prime }{{\tau }}_{xy}^{\prime }& -\omega {{\sigma }}_{y}^{\prime }{{\tau }}_{xy}^{\prime }& -\omega {{\sigma }}_{z}^{\prime }{{\tau }}_{xy}^{\prime }& \dfrac{1}{2}-\omega {{\tau }}_{xy}^{{\prime }2}& & \\ -\omega {{\sigma }}_{z}^{\prime }{{\tau }}_{yz}^{\prime }& -\omega {{\sigma }}_{y}^{\prime }{{\tau }}_{yz}^{\prime }& -\omega {{\sigma }}_{z}^{\prime }{{\tau }}_{yz}^{\prime }& -\omega {{\tau }}_{xy}^{\prime }{{\tau }}_{yz}^{\prime }& \dfrac{1}{2}-\omega {{\tau }}_{yz}^{{\prime }2}& \\ -\omega {{\sigma }}_{z}^{\prime }{{\tau }}_{zx}^{\prime }& -\omega {{\sigma }}_{y}^{\prime }{{\tau }}_{zx}^{\prime }& -\omega {{\sigma }}_{z}^{\prime }{{\tau }}_{zx}^{\prime }& -\omega {{\tau }}_{xy}^{\prime }{{\tau }}_{zx}^{\prime }& -\omega {{\tau }}_{xy}^{\prime }{{\tau }}_{zx}^{\prime }& \dfrac{1}{2}-\omega {{\tau }}_{zx}^{{\prime }2}\end{array}\right], $$ (22) 其中
$$ \omega = \frac{{9G}}{{2{{\bar \sigma }^2}( {{E_{\rm{p}}} + 3G} )}}{\text{.}} $$ (23) 式(23)中，
$ G $ 为纤维单元的剪切模量；$ \bar \sigma $ 为通过式(17)所求得的等效应力；$ {E_{\rm{p}}} $ 为塑性模量，它与弹性模量$ E $ 和切线模量$ {E_{\rm{t}}} $ 之间的关系为$ {E_{\rm{p}}} ={{E{E_{\rm{t}}}}}/({{E - {E_{\rm{t}}}}}) $ ，切线模量可通过对本构关系的求导得到。平均应力、应力偏量按式(24)计算：
$$ \left\{ \begin{aligned} & {\sigma _{\rm{m}}} = \frac{1}{3}( {{\sigma _x} + {\sigma _y} + {\sigma _z}} ) ,\\& {{\sigma }'_x} = {\sigma _x} - {\sigma _{\rm{m}}} ,\\& {{\sigma }'_y} = {\sigma _y} - {\sigma _{\rm{m}}} ,\\& {{\sigma }'_z} = {\sigma _z} - {\sigma _{\rm{m}}} ,\\& {{\tau }'_{xy}} = {\tau _{xy}} ,\\& {{\tau }'_{yz}} = {\tau _{yz}} ,\\& {{\tau }'_{zx}} = {\tau _{zx}} {\text{.}} \end{aligned} \right. $$ (24) 然后，计算应力分量增量
$ {\text{d}}{\boldsymbol{\sigma }} $ ：$$ {\text{d}} {\boldsymbol{\sigma }} = {{\boldsymbol{D}}_{{\rm{ep}}}}{\text{d}} {\boldsymbol{\varepsilon }} {\text{.}} $$ (25) 最后，按第5）步更新应力。
5） 得到纤维单元的应力增量后，更新应力：
$$ {{\boldsymbol{\sigma '}}} = {\boldsymbol{\sigma }} + {\text{d}} {\boldsymbol{\sigma }} {\text{.}} $$ (26) 更新材料应力状态的流程图如图2所示。
图 2 增量理论应力更新流程图Figure 2. The flow chart for stress updating with the incremental theory3. 算 例 分 析
3.1 算例一
一正方形截面钢筋混凝土柱[15]，如图3所示，柱高1.65 m，截面尺寸为550 mm×550 mm，纵向受力钢筋直径为20 mm，等距分布，保护层厚度为40 mm，考虑结构自重的影响，将重力等效为柱顶的集中力
$ P $ ，同时，在顶部作用一水平集中荷载$ V $ ，考虑几何非线性和材料非线性，对此结构做水平力-侧移的非线性全过程分析。算例分析采用本文所给出的三维纤维梁单元建模计算，计算时将原结构划分为10个单元，单元的每个节点有6个自由度；划分纤维单元时，尽可能保证原钢筋面积与纤维网格的面积相等，优先考虑此条件；将截面划分为31×31的正方形纤维网格，网格按照结构的实际情况分为钢筋纤维和混凝土纤维，如图4所示，其余参数见表1。
表 1 计算参数Table 1. Calculation parametersparameter value compressive strength of concrete $ f_{\rm{c}}^\prime = 32.1\;{\text{MPa}} $ Poisson’s ratio of concrete $ {\nu _{\rm{c}}} = 0.2 $ yield strength of rebar $ {f_{\rm{y}}} = 510\;{\text{MPa}} $ sectional shear coefficient $ {\kappa _y} = {\kappa _z} = 6/5 $ elastic modulus of concrete $ {E_{\rm{c}}} = 2.64 \times {10^4}\;{\text{MPa}} $ Poisson’s ratio of rebar $ {\nu _{\rm{s}}} = 0.3 $ elastic modulus of rebar $ {E_{\rm{s}}} = 2 \times {10^5}\;{\text{MPa}} $ peak compressive strain of the member $ \varepsilon {\text{ = }}0.005\;2 $ shear modulus of concrete $ {G_{\rm{c}}} = 1.1 \times {10^4}\;{\text{MPa}} $ moment of inertia $ I = 7.63 \times {10^{ - 3}}\;{{\text{m}}^4} $ shear modulus of rebar $ {G_{\rm{s}}} = 7.7 \times {10^4}\;{\text{MPa}} $ equivalent weight of the element $ P = 3\;539.25\;{\text{kN}} $ 本文迭代求解的方法采用杜修力等提出的位移控制新方法[16]，控制点为水平力施加点，位移增量控制步步长
$ \Delta \bar u = 1\;{\text{mm}} $ ，总位移增量步数$ n = 25 $ ，迭代收敛准则采用位移收敛准则，允许误差$ \delta {\text{ = }}{10^{ - 3}} $ 。3.1.1 混凝土本构模型
考虑到混凝土受到箍筋约束作用，故在数值模拟中，受拉区混凝土采用经Scott等[17]修正后的Kent-Park[18]单轴混凝土本构模型，如图5所示。需要说明的是，对于受拉区混凝土，经试算表明，可以忽略其受拉贡献，即取受拉混凝土纤维的拉应力为0，这一简化在不影响计算精度的前提下提高了计算效率。混凝土抗压强度及相应的混凝土峰值应变和极限应变均取文献[15]中所给出的数值。混凝土受压本构模型具体表达式如下：
$$ \sigma = \left\{ \begin{aligned} & K{{f}'_{\rm{c}}}\left[ {\frac{{2\varepsilon }}{{{\varepsilon _{\rm{c}}}}} - {{\left( {\frac{\varepsilon }{{{\varepsilon _{\rm{c}}}}}} \right)}^2}} \right],\,\qquad\quad \varepsilon \leqslant {\varepsilon _{\rm{c}}}, \hfill \\& K{{f}_{\rm{c}}'} \left[ {1 - {Z_{\rm{m}}}\left( {\varepsilon - {\varepsilon _{\rm{c}}}} \right)} \right],\qquad {\varepsilon _{\rm{c}}} \lt \varepsilon \leqslant {\varepsilon _{\rm{u}}}, \hfill \\& 0.2Kf_{\rm{c}}^\prime ,\;\quad\qquad\qquad\qquad \varepsilon > {\varepsilon _{\rm{u}}}, \end{aligned} \right. $$ (27) 其中
$$ K = 1 + \frac{{{\rho _{{\rm{sv}}}}{f_{{\rm{yh}}}}}}{{f_{\rm{c}}^\prime }} ,$$ (28) $$ {Z_{\rm{m}}} = \dfrac{{0.5}}{{\dfrac{{3 + 0.29f_{\rm{c}}^\prime }}{{145f_{\rm{c}}^\prime - 1\;000}} + 0.75{\rho _{{\rm{sv}}}}\sqrt {\dfrac{{h'}}{s}} - 0.002K}} ,$$ (29) 式中，
$ K $ 为考虑箍筋约束效应所引起的混凝土强度增大系数，根据本算例参数可计算得$ K{\text{ = }}1.21 $ ；$ {\rho _{{\rm{sv}}}} $ 为箍筋的体积配筋率；$ {f_{{\rm{yh}}}} $ 为箍筋的屈服强度；$ {\varepsilon _{\rm{u}}} $ 为混凝土极限压应变，本文取$ {\varepsilon _{\rm{u}}} = 0.024\;8 $ ；$ {Z_{\rm{m}}} $ 为约束混凝土应力应变曲线下降段的斜率，经式(29)计算，本文取$ {Z_{\rm{m}}} = 13.66 $ 。3.1.2 钢筋本构模型
钢筋的本构模型采用双折线弹塑性本构模型，如图6所示，此本构模型具体表达式如下：
$$ {\sigma _{\rm{s}}} = \left\{ \begin{aligned} & {E_{\rm{s}}}{\varepsilon _{\rm{s}}},\qquad\qquad\qquad{\varepsilon _{\rm{s}}} \leqslant {\varepsilon _{\rm{y}}}, \hfill \\& {f_{\rm{y}}} + k( {{\varepsilon _{\rm{s}}} - {\varepsilon _{\rm{y}}}} ),\;\qquad{\varepsilon _{\rm{y}}} \lt {\varepsilon _{\rm{s}}} \leqslant {\varepsilon _{\rm{u}}}, \hfill \\& 0,\;\qquad\qquad\quad\qquad{\varepsilon _{\rm{s}}} > {\varepsilon _{\rm{u}}} , \end{aligned}\right. $$ (30) 式中，
$ {E_{\rm{s}}} $ 为钢筋的弹性模量；$ {f_{\rm{y}}} $ 为钢筋的屈服强度；$ {\varepsilon _{\rm{y}}} $ 为钢筋的屈服应变；$ k $ 为钢筋硬化段斜率，本文与文献[15]相同，取$ k{\text{ = }}0.01{E_{\rm{s}}} $ 。3.1.3 数值模拟结果
利用弹塑性增量理论结合考虑剪切效应下的三维纤维梁单元模型，编制相应的MATLAB程序计算，将数值模拟结果与文献[15]的数据进行对比分析，如图7和表2所示。
表 2 数值模拟结果Table 2. Numerical simulation resultsdisplacement u/mm ref. [15] FR/kN numerical simulation FN/kN error δ/% 2 167.16 156.02 6.66 4 313.14 294.85 5.84 6 417.97 409.75 1.97 8 490.91 487.48 0.70 10 540.35 537.67 0.50 12 568.83 558.61 1.80 14 577.99 570.97 1.22 16 573.22 555.43 3.10 18 562.80 534.98 4.94 20 552.04 515.57 6.61 22 538.23 497.12 7.64 24 517.62 479.59 7.35 由图7和表2可见：本文的计算结果与文献[15]结果吻合较好，水平荷载的全过程有限元数值模拟结果略低于文献[15]的结果，其相对误差均在8%以内，结构的弹性阶段、弹塑性阶段以及下降段均较为明显。上述结果说明，本文所建立的复杂应力状态下考虑剪切效应的纤维梁单元理论是正确的，依据该理论所编制的非线性有限元分析程序，能有效地解决实际问题。
3.2 算例二
一矩形钢管混凝土压弯剪构件[19]，试件长度
$ L = 1\;200\;{\text{mm}} $ ，截面高$ D = 300\;{\text{mm}} $ ，截面宽$ B = 200\;{\text{mm}} $ ，外钢管厚度$ t = 5.6\;{\text{mm}} $ ，含钢率$ \alpha = 0.1 $ ，所用混凝土强度等级为C60，钢管为Q345钢，试件端部所受轴压$ N = 0.4{N_{\rm{u}}} $ （$ {N_{\rm{u}}} $ 为该构件的轴压极限承载力），跨中强轴方向受水平荷载$ V $ 。在数值模拟有限元分析模型中，试件两端均为铰接，钢材的本构模型采用双折线模型，混凝土的本构模型采用塑性损伤模型[20]。计算简图如图8所示，采用本文所给出的考虑剪切效应的纤维梁单元模型建模，使用位移控制新方法[16]进行有限元非线性方程的迭代求解，计算水平荷载$ V $ 与支座处转角位移$R(R = 2\varDelta /L ，\varDelta$ 为外荷载$ V $ 产生的水平位移）的全过程关系曲线。
图 8 矩形钢管混凝土压弯剪试件计算模型Figure 8. The calculation model for the rectangular CFST column under the compression-bending-shear loading condition图9为本文数值模拟结果与文献[19]的V-R曲线。由图可知，两者在弹性阶段吻合程度良好，在弹塑性阶段和下降段，由于纤维梁单元模型不能模拟钢管和混凝土的相互作用，本文的计算结果稍滞后于文献[19]的结果，但偏差相对较小，从总体上来说，两者结果较为吻合。此外，本文计算出的构件所受最大横向力为479.8 kN，文献[19]的最大横向力约为480.1 kN，说明本文的纤维梁单元模型能较好地模拟构件的横向极限承载力。
综上所述，本文所建立的复杂应力状态下结构非线性分析理论是正确的，证明了依据此理论所编制的非线性程序是可行、有效的。
4. 结 论
1） 算例分析的结果表明，本文基于考虑剪切效应的纤维梁单元，所建立的压弯剪复杂应力状态下结构非线性有限元分析理论是通用、可行和正确的。
2） 利用本文理论进行结构的非线性全过程分析时，可较好地模拟出结构的弹性阶段、弹塑性阶段和下降段，且特征点较为明显。
3） 依据本文理论所编制的非线性有限元分析程序是有效、可靠的，可以较好地解决理论研究和工程应用中所涉及的结构非线性问题。
-
图 2 增量理论应力更新流程图
Figure 2. The flow chart for stress updating with the incremental theory
图 8 矩形钢管混凝土压弯剪试件计算模型
Figure 8. The calculation model for the rectangular CFST column under the compression-bending-shear loading condition
表 1 计算参数
Table 1. Calculation parameters
parameter value compressive strength of concrete $ f_{\rm{c}}^\prime = 32.1\;{\text{MPa}} $ Poisson’s ratio of concrete $ {\nu _{\rm{c}}} = 0.2 $ yield strength of rebar $ {f_{\rm{y}}} = 510\;{\text{MPa}} $ sectional shear coefficient $ {\kappa _y} = {\kappa _z} = 6/5 $ elastic modulus of concrete $ {E_{\rm{c}}} = 2.64 \times {10^4}\;{\text{MPa}} $ Poisson’s ratio of rebar $ {\nu _{\rm{s}}} = 0.3 $ elastic modulus of rebar $ {E_{\rm{s}}} = 2 \times {10^5}\;{\text{MPa}} $ peak compressive strain of the member $ \varepsilon {\text{ = }}0.005\;2 $ shear modulus of concrete $ {G_{\rm{c}}} = 1.1 \times {10^4}\;{\text{MPa}} $ moment of inertia $ I = 7.63 \times {10^{ - 3}}\;{{\text{m}}^4} $ shear modulus of rebar $ {G_{\rm{s}}} = 7.7 \times {10^4}\;{\text{MPa}} $ equivalent weight of the element $ P = 3\;539.25\;{\text{kN}} $ 
下载: 导出CSV
表 2 数值模拟结果
Table 2. Numerical simulation results
displacement u/mm ref. [15] FR/kN numerical simulation FN/kN error δ/% 2 167.16 156.02 6.66 4 313.14 294.85 5.84 6 417.97 409.75 1.97 8 490.91 487.48 0.70 10 540.35 537.67 0.50 12 568.83 558.61 1.80 14 577.99 570.97 1.22 16 573.22 555.43 3.10 18 562.80 534.98 4.94 20 552.04 515.57 6.61 22 538.23 497.12 7.64 24 517.62 479.59 7.35
下载: 导出CSV
-
[1] 马春来. 考虑剪切位移和多因素耦合影响下的空间梁非线性分析模型研究[D]. 硕士学位论文. 杭州: 浙江大学, 2008.MA Chunlai. Nonlinear analysis of space beam based on shear displacement and multi-factors coupling influence[D]. Master Thesis. Hangzhou: Zhejiang University, 2008. (in Chinese) [2] TIMOSHENKO S P. On the correction for shear of the differential equation for transverse vibrations of prismatic bars[J]. Philosophical Magazine, 1921, 41(245): 744-746. [3] MARI A, SCORDELIS A. Nonlinear geometric material and time dependent analysis of three dimensional reinforced and prestressed concrete frames[R]. SESM Report. Department of Civil Engineering, University of California, 1984. [4] 胡郑州, 吴明儿. 考虑剪切效应三维纤维梁单元的几何非线性增量有限元分析[J]. 工程力学, 2014, 31(8): 134-143. (HU Zhengzhou, WU Ming’er. Geometrically nonlinear incremental finite element analysis of 3D fiber beam element considering shear effect[J]. Engineering Mechanics, 2014, 31(8): 134-143.(in Chinese) doi: 10.6052/j.issn.1000-4750.2013.03.0170 [5] 顾观东. 基于变形位移理论的空间梁柱非线性有限元分析模型研究[D]. 硕士学位论文. 杭州: 浙江大学, 2013.GU Guandong. Study on nonlinear spatial beam element model based on theory of deformation and displacement[D]. Master Thesis. Hangzhou: Zhejiang University, 2013.(in Chinese) [6] 郑欣怡. 考虑粘结滑移效应的新型钢筋混凝土纤维梁柱单元模型[D]. 硕士学位论文. 哈尔滨: 哈尔滨工业大学, 2019.ZHENG Xinyi. A new reinforced concrete fiber beam-column element model considering bond-slip effect[D]. Master Thesis. Harbin: Harbin Institute of Technology, 2019.(in Chinese) [7] THAI H T, KIM S E. Second-order inelastic analysis of cable-stayed bridges[J]. Finite Elements in Analysis & Design, 2012, 53(7): 48-55. [8] 张传超, 郑山锁, 李磊, 等. 基于柔度法的纤维梁柱单元及其参数分析[J]. 工业建筑, 2010, 40(12): 90-94. (ZHANG Chuanchao, ZHENG Shansuo, LI Lei, et al. Flexibility-based fiber beam-column element and its parametric analysis[J]. Industrial Construction, 2010, 40(12): 90-94.(in Chinese) [9] 马银. 基于纤维模型的型钢混凝土梁柱单元理论[D]. 硕士学位论文. 西安: 西安建筑科技大学, 2010.MA Yin. Element theory of steel reinforced concrete beam and column based on fiber model[D]. Master Thesis. Xi’an: Xi’an University of Architecture and Technology, 2010. (in Chinese) [10] 蔺鹏臻, 王富平, 柳兴成. 基于纤维梁模型的钢筋混凝土箱梁非线性分析[J]. 铁道建筑, 2019, 59(5): 22-26, 55. (LIN Pengzhen, WANG Fupin, LIU Xingcheng. Nonlinear analysis of RC box girder based on fiber beam model[J]. Railway Construction, 2019, 59(5): 22-26, 55.(in Chinese) doi: 10.3969/j.issn.1003-1995.2019.05.05 [11] KAGERMANOV A, CERESA P. 3D fiber-based frame element with multiaxial stress interaction for RC structures[J]. Advances in Civil Engineering, 2018, 2018: 8596970. [12] PARK K, KIM H, DAE J, et al. Generalized finite element formulation of fiber beam elements for distributed plasticity in multiple regions[J]. Computer-Aided Civil and Infrastructure Engineering, 2019, 34(2): 146-163. doi: 10.1111/mice.12389 [13] GENDY S S F M, AYOUB A. Explicit fiber beam-column elements for impact analysis of structures[J]. Journal of Structural Engineering, 2018, 144(7): 04018068. [14] BAZOUNE A, KHULIEF Y A, STEPHEN N G. Shape functions of three-dimensional Timoshenko beam element[J]. Journal of Sound and Vibration, 2003, 259(2): 473-480. doi: 10.1006/jsvi.2002.5122 [15] AL-AUKAILY A, SCOTT M H. Sensitivity analysis for displacement-controlled finite-element analyses[J]. Journal of Structural Engineering, 2018, 144(3): 04017222. doi: 10.1061/(ASCE)ST.1943-541X.0001983 [16] 杜修力, 曹惠, 金浏. 力-变位关系全过程模拟的有限元位移控制新方法[J]. 工程力学, 2012, 29(1): 1-6. (DU Xiuli, CAO Hui, JIN Liu. New finite element displacement control method for force-displacement relationship simulation in the whole process[J]. Engineering Mechanics, 2012, 29(1): 1-6.(in Chinese) [17] SCOTT B D, PARK R, PRIESTLEY M J N. Stress-strain behavior of concrete confined by overlapping hoops at low and high strain rates[J]. ACI Journal, 1982, 79(1): 13-27. [18] KENT D C, PARK R. Flexural members with confined concrete[J]. Journal of the Structural Division, 1971, 97(7): 1969-1990. [19] 王亚晋, 范重, 侯超, 等. 矩形钢管混凝土构件双向压弯剪受力性能研究[J]. 建筑结构学报, 2019, 40(S1): 257-263. (WANG Yajin, FAN Zhong, HOU Chao, et al. Research on the behavior of rectangular CFST members under bidirectional compression and bend sears[J]. Journal of Building Structures, 2019, 40(S1): 257-263.(in Chinese) [20] 韩林海. 钢管混凝土结构-理论与实践[M]. 北京: 科学出版社, 2004.HAN Linhai. Concrete-Filled Steel Tubular Structure-Theory and Practice[M]. Beijing: Science Press, 2004. (in Chinese) 期刊类型引用(4)
1. 马今伟，段庆林. 无额外自由度广义有限元的近不可压弹-塑性分析. 应用数学和力学. 2024(02): 220-226 . 
本站查看2. 王云超，师振贵. 风力发电石油电动旋挖钻机功率模糊控制. 炼油技术与工程. 2024(07): 41-45 . 百度学术3. 张涛，吴波，张佳惠. 钢筋混凝土结构在动态载荷下的抗震检测. 砖瓦. 2024(09): 101-104 . 百度学术4. 马文金，曾其权，张淑兴. 基于有限元法的核电蓄电池柜抗震分析与评估. 深圳大学学报(理工版). 2023(06): 732-740 . 百度学术其他类型引用(3)
-
下载:
计量
- 文章访问数: 947
- HTML全文浏览量: 344
- PDF下载量: 97
- 被引次数: 7
渝公网安备50010802005915号