基于APDL的ANSYS嵌入式配筋实现

杜永峰+张玉星+朱前坤

摘要： 以ANSYS为平台，应用APDL实现钢筋混凝土有限元模型中钢筋单元嵌入混凝土单元的功能.应用该方法可以在ANSYS中方便地建立复杂的钢筋混凝土分离式模型.利用NewtonRaphson方法求解等参数单元的坐标插值函数中被约束节点在自然坐标系下的坐标值，得到自然坐标系下节点对应形函数值，并将对应的形函数值作为约束方程中混凝土单元节点位移的分项因数，从而实现嵌入节点与被嵌入单元位移协调.计算结果与Marc的算例结果吻合良好，说明该方法合理.

关键词： 钢筋混凝土； 嵌入式配筋； 分离式模型； NewtonRaphson方法； 约束方程； ANSYS； APDL

中图分类号： TU311； P315.9文献标志码： B

0引言

通用有限元程序ANSYS和Marc等可以比较方便地实现嵌入式配筋.ANSYS在分析钢筋混凝土构件时有2种方法：一种为弥散式配筋，另一种为分离式配筋.弥散式配筋的优点为建立模型简单：将钢筋本构矩阵D乘以体积配筋率，叠加在混凝土应力、应变矩阵上，形成弥散单元的本构矩阵.这种方法需要计算单元的体积配筋率，如果单元划分不规则或配筋复杂，不容易对单元设定相应的配筋率，并且在后处理可视化中不容易观察钢筋的受力情况.分离式配筋在ANSYS中实现有3种方法，分别为实体切分法、节点耦合法和约束方程法.[1]这3种方法的实现需要建立的混凝土单元节点与钢筋节点在空间中的距离较为接近：实体切分需要完全共用节点，节点耦合是将距离最近的钢筋节点与混凝土节点耦合在一起，约束方程将钢筋节点位移用距离最近的单元面上的节点位移通过建立线性约束方程表示.这些方法在ANSYS中的实现都要考虑钢筋的位置，否则建立的模型计算结果偏差较大.本文利用有限元等参数单元的特点，应用NewtonRaphson方法求解被约束节点在自然坐标系下的坐标值，进而求得形函数在该点的值，并将所求值作为约束方程中的节点位移分项系数，以此实现嵌入节点与被嵌入单元的位移协调，实现单元嵌入的功能.

1计算方法

1.1等参数单元形函数

在有限元法中，一般用等参数单元对所要分析的模型进行离散计算，将物理坐标系下单元的坐标值和位移值用相同的形函数插值表示，形函数由自然坐标系下的坐标值表示.对于八节点六面体单元，其形函数基本形式[2]为Ni=18（1±s）（1±t）（1±r）（1）形函数通过物理坐标系下的单元节点坐标和单元节点位移，建立自然坐标系中标准单元与物理坐标系中实际划分单元之间的映射关系，这种映射关系一一对应，插值函数的雅克比矩阵的行列式值在求解域内不能为0，保证求解非线性方程组坐标值的唯一性.等参元的优点是保证整体坐标系中不规则网格边界位移的连续性，并通过映射关系将物理坐标系中的不定积分域转化为自然坐标系中的定域积分.[3]

单元坐标插值函数为f1=ni=1Nixi-x=0

f2=ni=1Niyi-y=0

f3=ni=1Nizi-z=0 （2）单元位移插值函数为u=ni=1Niui

v=ni=1Nivi

w=ni=1Niwi（3）通过求解式（2）中被约束节点在自然坐标系中的s，t和r坐标值，应用式（3）求解形函数在约束节点的值，将所求值作为所要建立约束方程中各节点位移的分项系数，实现嵌入节点与被嵌入单元位移的协调.

1.2NewtonRaphson方法

NewtonRaphson方法是求解非线性方程的方法，该方法通过迭代求解线性化的方程解决非线性方程的求解问题.[4]其基本步骤为：将非线性方程略去高阶项并展开为一阶线性方程，确定求解的初始值，在初始值基础上迭代求解不平衡项，当计算增量在设定许可范围时即可停止迭代；对非线性方程组，需要求解方程组的1阶偏导数矩阵即雅可比矩阵.设定初值为自然坐标下的坐标原点，八节点六面体单元的雅可比矩阵为2次的，同时NewtonRaphson方法具有2阶收敛速度[4]，所以对选定的形函数只需迭代1次即可获得足够精度的解.通过求解结果向量的2范数判断收敛.该方法的迭代格式为 f1sf1tf1r

f2sf2tf2r

f3sf3tf3rnΔs

Δt

Δr=-f1（sn）

f2（tn）

f3（rn），

sn+1=sn+Δs

tn+1=tn+Δt

rn+1=rn+Δr （4） 1.3约束方程

通过约束方程将被约束节点自由度用约束节点自由度线性表示，求解过程将整体刚度矩阵中的被约束节点自由度凝聚，只保留其刚度和节点载荷的影响，约束方程的基本形式为Cons t=ni=1（Coefficient（i）·u（i）） （5）Coefficient（i）为节点位移u（i）的分项系数，分项系数为被约束节点在自然坐标系下形函数的值.[1]

1.4APDL实现步骤

实现过程应用到ANSYS程序内部的函数命令，其中ENEARN是获得节点临近单元号的内部函数命令，*MOPER为求解线性方程组的命令.

基本步骤如下：

（1）首先选择嵌入单元和被嵌入单元，应用CM命令建立集合.

（2）应用*GET命令获取嵌入节点的数量，并建立嵌入节点号数组.

（3）应用*DO循环和内部函数命令ENEARN得到每个嵌入节点临近的嵌入单元号，并存入数组.

（4）获得嵌入单元节点与被嵌入单元节点在物理坐标系下的坐标值，应用*IF命令判断单元形函数所属类型，选择相应的形函数.

（5）将被嵌入单元节点与嵌入单元节点坐标值代入雅克比矩阵和位移插值函数形成的向量，应用*MOPER，*DO和*IF命令语句建立NewtonRaphson求解过程.[5]endprint

（6）设定初值0，求解s，t和r值代入位移插值函数形成各位移的分项系数，应用CE命令建立节点约束方程.

实现流程见图1.

图 1实现流程

Fig.1Implementation process

NR迭代APDL过程，应用*MOPER计算线性方程组，需要将雅可比矩阵和不平衡项FFZERO的第一项写入命令相应位置.[5]可以设定循环求解的次数，超出求解迭代次数认为不收敛，对本问题可以取值为3，观察最终的ITERTIME（I），即每个嵌入节点最终的求解迭代次数，结果显示均为1，由此可以说明NR方法的2阶收敛速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，FFZERO（1）

S（I）=S（I）+RES（1）

T（I）=T（I）+RES（2）

R（I）=R（I）+RES（3）

*IF，RES（1）**2+RES（2）**2+RES（3）**2，LE，1E6，THEN

*EXIT

*ELSE

*ENDIF

ITERTIME（I）=ITERTIME（I）+1

*ENDDO

2算例分析

2.1计算模型和参数

为验证计算方法的合理性，将应用Marc计算的结果与本文方法应用ANSYS计算的结果进行比较，通过设定具有相同的几何模型和物理参数的悬臂柱，对柱顶进行位移加载，比较柱顶的反力位移曲线.

柱截面取为300 mm×400 mm，保护层厚度取为20 mm，混凝土采用我国《混凝土结构设计规范》[6]规定的混凝土轴心受压应力应变关系，混凝土材料强度取C30，钢筋为Q235钢材，箍筋间距为100 mm，钢筋混凝土构件尺寸见图2.有限元模型见图3.

图 2混凝土悬臂柱尺寸， mm

Fig.2Concrete cantilever column size， mm

（a） Marc模型（b） ANSYS模型图 3有限元模型

Fig.3Finite element model

在有限元计算中，Marc需要输入减去弹性应变的等效塑性应力应变关系.[7]对于小应变问题，工程应变与真实应变的计算结果基本一致.本文输入为工程应力应变关系，输入的本构关系见图4.设定随动强化和等向强化材料属性，对单向加载中单元无卸载的情况计算结果没有影响.本文设定为随动强化模型，屈服准则均设定为von Mises屈服准则，钢筋采用理想弹塑性模型.

图 4混凝土单轴受压本构关系

Fig.4Concrete uniaxial compression constitutive relation

2.2计算结果分析

分3组模型进行计算：第一组为不嵌入钢筋的素混凝土，本构关系设定为不考虑开裂的弹塑性本构模型；第二组为考虑嵌入钢筋不设定开裂影响的计算模型；第三组为考虑嵌入钢筋和考虑开裂的计算模型.ANSYS混凝土为SOLID65单元，钢筋为LINK8单元，Marc实体选用八节点7号单元，杆单元为二节点9号单元.

2.2.1第一、二组计算对比

第一、二组计算结果为不考虑或考虑钢筋嵌入计算结果.通过将柱顶的节点进行耦合，对保留节点施加水平方向40 mm位移进行计算，结果见图5.

图 5第一组与第二组结果对比

Fig.5Result comparison of group 1 and group 2

对没有嵌入钢筋的素混凝土，2个程序的计算结果相同；对嵌入钢筋的计算结果，ANSYS的峰值点略高于Marc的计算结果，在峰值之前两者的计算结果相同.由此可以说明，通过本文的方法可以实现在ANSYS程序中嵌入配筋的计算.

2.2.2第三组计算对比

为进一步计算考虑混凝土加载过程受拉开裂的结果，设定Marc开裂的极限拉应力为2 MPa，软化模量为弹性模量的1/10，极限压应变为0.003 3，开裂单元的剪力传递因数为0.3.ANSYS通过设定破坏准则设定相同的开裂拉应力，张开裂缝剪力传递因数取0.3，闭合裂缝剪力传递因数取0.95.关闭压碎，同时设定SOLID65单元的单元选项为1和7，不考虑位移形函数附加项和考虑开裂后的拉应力释放，并将收敛准则CNVTOL设定为位移判定，保证计算收敛.[8]第三组结果对比见图6.在没有开裂前的弹性段，2个程序的计算结果相同；在开裂之后，计算的反力ANSYS要小于Marc，计算的结果趋势相近.同时比较开裂与不考虑开裂的刚度和峰值，结果相差较大，开裂峰值反力在20 kN左右，不开裂峰值反力在80 kN左右.

图 6第三组结果对比

Fig.6Result comparison of group 3

3结束语

通过有限元理论，应用等参数单元的坐标插值函数，利用NR方法求解自然坐标系中的节点坐标值，并将求解结果代入位移插值函数，应用线性约束方程将嵌入单元节点的位移用被嵌入单元节点的位移表示，应用ANSYS实现实体单元嵌入杆单元的功能.与Marc的计算结果进行对比，结果表明：在不考虑混凝土开裂情况下，两者计算结果吻合良好，验证本文方法的合理性.2种软件对开裂的处理有一定的差别，考虑开裂后计算结果虽然趋势相近，但计算结果仍有一定差距.参考文献：

[1]王新敏. ANSYS工程结构数值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李义强， 许宏伟. ANSYS结构分析单元与应用[M]. 北京： 人民交通出版社， 2011： 187195.

[3]李人宪. 有限元法基础[M]. 2版. 北京： 国防工业出版社， 2004.

[4]殷有泉. 非线性有限元基础[M]. 北京： 北京大学出版社， 2007： 18.

[5]博弈创作室. APDL参数化有限元分析技术及其应用实例[M]. 北京： 中国水利水电出版社， 2004： 9094.

[6]GB 50010—2010混凝土结构设计规范[S].

[7]陆新征， 叶列平， 缪志伟， 等. 建筑抗震弹塑性分析[M]. 北京： 中国建筑工业出版社， 2009： 169202.

[8]陆新征， 江见鲸. 用ANSYS SOLID65单元分析混凝土组合构件复杂应力[J]. 建筑结构， 2003， 33（6）： 2224.

LU Xinzheng， JIANG Jianjing. ANSYS SOLID65 element analysis in concrete composite components with complex stress state[J]. Building Struct， 2003， 33（6）： 2224.

（编辑武晓英）endprint

（6）设定初值0，求解s，t和r值代入位移插值函数形成各位移的分项系数，应用CE命令建立节点约束方程.

实现流程见图1.

图 1实现流程

Fig.1Implementation process

NR迭代APDL过程，应用*MOPER计算线性方程组，需要将雅可比矩阵和不平衡项FFZERO的第一项写入命令相应位置.[5]可以设定循环求解的次数，超出求解迭代次数认为不收敛，对本问题可以取值为3，观察最终的ITERTIME（I），即每个嵌入节点最终的求解迭代次数，结果显示均为1，由此可以说明NR方法的2阶收敛速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，FFZERO（1）

S（I）=S（I）+RES（1）

T（I）=T（I）+RES（2）

R（I）=R（I）+RES（3）

*IF，RES（1）**2+RES（2）**2+RES（3）**2，LE，1E6，THEN

*EXIT

*ELSE

*ENDIF

ITERTIME（I）=ITERTIME（I）+1

*ENDDO

2算例分析

2.1计算模型和参数

为验证计算方法的合理性，将应用Marc计算的结果与本文方法应用ANSYS计算的结果进行比较，通过设定具有相同的几何模型和物理参数的悬臂柱，对柱顶进行位移加载，比较柱顶的反力位移曲线.

柱截面取为300 mm×400 mm，保护层厚度取为20 mm，混凝土采用我国《混凝土结构设计规范》[6]规定的混凝土轴心受压应力应变关系，混凝土材料强度取C30，钢筋为Q235钢材，箍筋间距为100 mm，钢筋混凝土构件尺寸见图2.有限元模型见图3.

图 2混凝土悬臂柱尺寸， mm

Fig.2Concrete cantilever column size， mm

（a） Marc模型（b） ANSYS模型图 3有限元模型

Fig.3Finite element model

在有限元计算中，Marc需要输入减去弹性应变的等效塑性应力应变关系.[7]对于小应变问题，工程应变与真实应变的计算结果基本一致.本文输入为工程应力应变关系，输入的本构关系见图4.设定随动强化和等向强化材料属性，对单向加载中单元无卸载的情况计算结果没有影响.本文设定为随动强化模型，屈服准则均设定为von Mises屈服准则，钢筋采用理想弹塑性模型.

图 4混凝土单轴受压本构关系

Fig.4Concrete uniaxial compression constitutive relation

2.2计算结果分析

分3组模型进行计算：第一组为不嵌入钢筋的素混凝土，本构关系设定为不考虑开裂的弹塑性本构模型；第二组为考虑嵌入钢筋不设定开裂影响的计算模型；第三组为考虑嵌入钢筋和考虑开裂的计算模型.ANSYS混凝土为SOLID65单元，钢筋为LINK8单元，Marc实体选用八节点7号单元，杆单元为二节点9号单元.

2.2.1第一、二组计算对比

第一、二组计算结果为不考虑或考虑钢筋嵌入计算结果.通过将柱顶的节点进行耦合，对保留节点施加水平方向40 mm位移进行计算，结果见图5.

图 5第一组与第二组结果对比

Fig.5Result comparison of group 1 and group 2

对没有嵌入钢筋的素混凝土，2个程序的计算结果相同；对嵌入钢筋的计算结果，ANSYS的峰值点略高于Marc的计算结果，在峰值之前两者的计算结果相同.由此可以说明，通过本文的方法可以实现在ANSYS程序中嵌入配筋的计算.

2.2.2第三组计算对比

为进一步计算考虑混凝土加载过程受拉开裂的结果，设定Marc开裂的极限拉应力为2 MPa，软化模量为弹性模量的1/10，极限压应变为0.003 3，开裂单元的剪力传递因数为0.3.ANSYS通过设定破坏准则设定相同的开裂拉应力，张开裂缝剪力传递因数取0.3，闭合裂缝剪力传递因数取0.95.关闭压碎，同时设定SOLID65单元的单元选项为1和7，不考虑位移形函数附加项和考虑开裂后的拉应力释放，并将收敛准则CNVTOL设定为位移判定，保证计算收敛.[8]第三组结果对比见图6.在没有开裂前的弹性段，2个程序的计算结果相同；在开裂之后，计算的反力ANSYS要小于Marc，计算的结果趋势相近.同时比较开裂与不考虑开裂的刚度和峰值，结果相差较大，开裂峰值反力在20 kN左右，不开裂峰值反力在80 kN左右.

图 6第三组结果对比

Fig.6Result comparison of group 3

3结束语

通过有限元理论，应用等参数单元的坐标插值函数，利用NR方法求解自然坐标系中的节点坐标值，并将求解结果代入位移插值函数，应用线性约束方程将嵌入单元节点的位移用被嵌入单元节点的位移表示，应用ANSYS实现实体单元嵌入杆单元的功能.与Marc的计算结果进行对比，结果表明：在不考虑混凝土开裂情况下，两者计算结果吻合良好，验证本文方法的合理性.2种软件对开裂的处理有一定的差别，考虑开裂后计算结果虽然趋势相近，但计算结果仍有一定差距.参考文献：

[1]王新敏. ANSYS工程结构数值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李义强， 许宏伟. ANSYS结构分析单元与应用[M]. 北京： 人民交通出版社， 2011： 187195.

[3]李人宪. 有限元法基础[M]. 2版. 北京： 国防工业出版社， 2004.

[4]殷有泉. 非线性有限元基础[M]. 北京： 北京大学出版社， 2007： 18.

[5]博弈创作室. APDL参数化有限元分析技术及其应用实例[M]. 北京： 中国水利水电出版社， 2004： 9094.

[6]GB 50010—2010混凝土结构设计规范[S].

[7]陆新征， 叶列平， 缪志伟， 等. 建筑抗震弹塑性分析[M]. 北京： 中国建筑工业出版社， 2009： 169202.

[8]陆新征， 江见鲸. 用ANSYS SOLID65单元分析混凝土组合构件复杂应力[J]. 建筑结构， 2003， 33（6）： 2224.

LU Xinzheng， JIANG Jianjing. ANSYS SOLID65 element analysis in concrete composite components with complex stress state[J]. Building Struct， 2003， 33（6）： 2224.

（编辑武晓英）endprint

（6）设定初值0，求解s，t和r值代入位移插值函数形成各位移的分项系数，应用CE命令建立节点约束方程.

实现流程见图1.

图 1实现流程

Fig.1Implementation process

NR迭代APDL过程，应用*MOPER计算线性方程组，需要将雅可比矩阵和不平衡项FFZERO的第一项写入命令相应位置.[5]可以设定循环求解的次数，超出求解迭代次数认为不收敛，对本问题可以取值为3，观察最终的ITERTIME（I），即每个嵌入节点最终的求解迭代次数，结果显示均为1，由此可以说明NR方法的2阶收敛速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，FFZERO（1）

S（I）=S（I）+RES（1）

T（I）=T（I）+RES（2）

R（I）=R（I）+RES（3）

*IF，RES（1）**2+RES（2）**2+RES（3）**2，LE，1E6，THEN

*EXIT

*ELSE

*ENDIF

ITERTIME（I）=ITERTIME（I）+1

*ENDDO

2算例分析

2.1计算模型和参数

为验证计算方法的合理性，将应用Marc计算的结果与本文方法应用ANSYS计算的结果进行比较，通过设定具有相同的几何模型和物理参数的悬臂柱，对柱顶进行位移加载，比较柱顶的反力位移曲线.

柱截面取为300 mm×400 mm，保护层厚度取为20 mm，混凝土采用我国《混凝土结构设计规范》[6]规定的混凝土轴心受压应力应变关系，混凝土材料强度取C30，钢筋为Q235钢材，箍筋间距为100 mm，钢筋混凝土构件尺寸见图2.有限元模型见图3.

图 2混凝土悬臂柱尺寸， mm

Fig.2Concrete cantilever column size， mm

（a） Marc模型（b） ANSYS模型图 3有限元模型

Fig.3Finite element model

在有限元计算中，Marc需要输入减去弹性应变的等效塑性应力应变关系.[7]对于小应变问题，工程应变与真实应变的计算结果基本一致.本文输入为工程应力应变关系，输入的本构关系见图4.设定随动强化和等向强化材料属性，对单向加载中单元无卸载的情况计算结果没有影响.本文设定为随动强化模型，屈服准则均设定为von Mises屈服准则，钢筋采用理想弹塑性模型.

图 4混凝土单轴受压本构关系

Fig.4Concrete uniaxial compression constitutive relation

2.2计算结果分析

分3组模型进行计算：第一组为不嵌入钢筋的素混凝土，本构关系设定为不考虑开裂的弹塑性本构模型；第二组为考虑嵌入钢筋不设定开裂影响的计算模型；第三组为考虑嵌入钢筋和考虑开裂的计算模型.ANSYS混凝土为SOLID65单元，钢筋为LINK8单元，Marc实体选用八节点7号单元，杆单元为二节点9号单元.

2.2.1第一、二组计算对比

第一、二组计算结果为不考虑或考虑钢筋嵌入计算结果.通过将柱顶的节点进行耦合，对保留节点施加水平方向40 mm位移进行计算，结果见图5.

图 5第一组与第二组结果对比

Fig.5Result comparison of group 1 and group 2

对没有嵌入钢筋的素混凝土，2个程序的计算结果相同；对嵌入钢筋的计算结果，ANSYS的峰值点略高于Marc的计算结果，在峰值之前两者的计算结果相同.由此可以说明，通过本文的方法可以实现在ANSYS程序中嵌入配筋的计算.

2.2.2第三组计算对比

为进一步计算考虑混凝土加载过程受拉开裂的结果，设定Marc开裂的极限拉应力为2 MPa，软化模量为弹性模量的1/10，极限压应变为0.003 3，开裂单元的剪力传递因数为0.3.ANSYS通过设定破坏准则设定相同的开裂拉应力，张开裂缝剪力传递因数取0.3，闭合裂缝剪力传递因数取0.95.关闭压碎，同时设定SOLID65单元的单元选项为1和7，不考虑位移形函数附加项和考虑开裂后的拉应力释放，并将收敛准则CNVTOL设定为位移判定，保证计算收敛.[8]第三组结果对比见图6.在没有开裂前的弹性段，2个程序的计算结果相同；在开裂之后，计算的反力ANSYS要小于Marc，计算的结果趋势相近.同时比较开裂与不考虑开裂的刚度和峰值，结果相差较大，开裂峰值反力在20 kN左右，不开裂峰值反力在80 kN左右.

图 6第三组结果对比

Fig.6Result comparison of group 3

3结束语

通过有限元理论，应用等参数单元的坐标插值函数，利用NR方法求解自然坐标系中的节点坐标值，并将求解结果代入位移插值函数，应用线性约束方程将嵌入单元节点的位移用被嵌入单元节点的位移表示，应用ANSYS实现实体单元嵌入杆单元的功能.与Marc的计算结果进行对比，结果表明：在不考虑混凝土开裂情况下，两者计算结果吻合良好，验证本文方法的合理性.2种软件对开裂的处理有一定的差别，考虑开裂后计算结果虽然趋势相近，但计算结果仍有一定差距.参考文献：

[1]王新敏. ANSYS工程结构数值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李义强， 许宏伟. ANSYS结构分析单元与应用[M]. 北京： 人民交通出版社， 2011： 187195.

[3]李人宪. 有限元法基础[M]. 2版. 北京： 国防工业出版社， 2004.

[4]殷有泉. 非线性有限元基础[M]. 北京： 北京大学出版社， 2007： 18.

[5]博弈创作室. APDL参数化有限元分析技术及其应用实例[M]. 北京： 中国水利水电出版社， 2004： 9094.

[6]GB 50010—2010混凝土结构设计规范[S].

[7]陆新征， 叶列平， 缪志伟， 等. 建筑抗震弹塑性分析[M]. 北京： 中国建筑工业出版社， 2009： 169202.

[8]陆新征， 江见鲸. 用ANSYS SOLID65单元分析混凝土组合构件复杂应力[J]. 建筑结构， 2003， 33（6）： 2224.

LU Xinzheng， JIANG Jianjing. ANSYS SOLID65 element analysis in concrete composite components with complex stress state[J]. Building Struct， 2003， 33（6）： 2224.

（编辑武晓英）endprint