钢筋混凝土结构多尺度建模界面连接方法研究

吴怡凡， 沈建， 金世杰， 夏晋

（1.浙江大学结构工程研究所，杭州 310058；2.浙江大学建筑设计研究院有限公司，杭州 310028）

0 引言

随着有限元技术的不断发展，在大型土木工程结构的模拟分析中，为了满足整体结构的力学响应，并反应关键节点的塑性发展，多尺度建模成为了一种常用策略，即在整体结构的非关键部位使用计算效率较高的宏观单元，而在需要精细分析的关键节点部位使用精度较高的微观单元，并通过一定的连接方式将两种类型的单元耦合起来，以达到既兼顾计算效率又保证精度的目的［1-3］。多点约束法是目前土木工程领域进行多尺度界面连接最常用的一种方法，它的核心思想是通过建立不同尺度界面节点自由度间的协调方程使两界面达到有效的耦合［4］。在过去的研究中，学者们往往采用三种理论来构建多点约束方程，分别为：位移协调、能量协调以及位移协调与能量协调相结合。林旭川等运用平截面假定列出位移协调方程，采用MSC.MARC模拟了钢筋混凝土框架多尺度模型在地震下的受力情况，并通过试验验证［5］。Yue和Wu等［6，7］利用做功相等能量协调原理，分别通过ABAQUS和ANSYS建立了钢筋混凝土框架多尺度模型，并与实验比较，结果表明该方法可以准确模拟结构局部破坏。但在钢筋混凝土结构的数值模拟分析中，由于材料非线性、几何非线性等因素，基于单一理论的多点约束法存在一定的局限性［8-10］。基于位移协调与能量协调相结合的多点约束法可以避开两种单一理论的缺点得到更好的模拟结果，但需要插入数量庞大的约束方程，建模前处理繁琐复杂［11］。

文中提出了一种基于连续分布耦合理论的多点约束法以实现不同尺度单元间的连接，针对钢筋混凝土节点，使用ABAQUS比较分析了连续分布耦合与位移协调、能量协调、位移协调与能量协调相结合4种多尺度界面连接方法的优劣性。

1 多尺度界面连接方式分析

文中以梁-实体单元的连接为例介绍不同理论实现多尺度界面连接的基本原理和方法，梁-实体单元连接如图1所示。

图1 梁-实体单元连接

1.1 基于位移协调的界面连接法

连接界面上实体单元共有n个节点，每个节点有x、y、z3个自由度；梁单元与实体单元通过节点B连接，节点B有xB、yB、zB、θxB、θyB、θzB6个自由度。根据梁、实体单元可能产生的4种位移模式，可知界面节点间产生的变形协调如图2所示。将精细实体单元节点的3个自由度分别用宏观梁单元节点的6个自由度表示，一共形成3n个位移协调方程如式（1）。

图2 梁-实体单元多尺度界面变形协调

式中，bxi、byi为实体单元节点i在X、Y坐标轴上的坐标。

1.2 基于能量协调的界面连接法

在轴力、弯矩、剪力、扭矩的作用下界面节点变形示意如图3所示，根据虚功原理可列出做功相等方程如式（2）。

图3 梁-实体单元多尺度界面受力变形

式中，Fj为作用在宏观单元B节点方向的力；σi1,Fj为精细单元i节点处由Fj引起的轴向应力；τi2,Fj、τi3,Fj为精细单元i节点处由Fj引起的切向应力；u1S、u2S、u3S为精细单元i节点处轴向及切向的位移。

结合材料力学公式，并用形函数来表示节点位移X={N}{X}，Y={N}{Y}，Z={N}{Z}，代入四节点四边形单元形函数对积分方程求解，可得到约束方程如式（3）。

式中，Cxi、Cyi、Czi、Cθxi、Cθyi为与实体单元截面尺寸及连接界面点数个数相关的系数。

1.3 基于位移协调与能量协调相结合的界面连接法

基于位移协调的多点约束法在切向存在过约束的缺陷，而基于能量协调的多点约束法在转角方向存在材料非线性和几何非线性的缺陷。

因此将两种方法在不同自由度上进行方程联立，在切向上使用基于能量协调的多点约束方程，而在轴向和转角方向上使用基于位移协调的多点约束方程，联立得到的方程组如式（4）。

1.4 基于连续分布耦合的界面连接法

连续分布耦合示意如图4所示，连续分布耦合允许约束面上各节点发生相对移动，通过对约束面上各节点的运动进行加权处理，使在此区域上受到的合力和合力距与施加在参考点上的力和力矩相等效。因此，使用连续分布耦合理论建立的多尺度模型在理论上可以避免位移协调理论产生的切向过约束。ABAQUS连续分布耦合有4种加权系数计算方式，包括一致、线性、二次方和三次方，计算公式如式（5）～式（8）。经过反复试算文中采用一致加权系数计算方式对界面节点进行加权处理。

图4 连续分布耦合

式中，ωi为界面节点i的加权系数；ri为界面节点i到参考点的径向距离；ω0为最远点的径向距离。

2 不同多尺度界面连接法算例分析

分析模型如图5所示，构件高2m，截面为边长0.4m的正方形，多尺度界面距离底部1m。图5（a）为精细化模型，图5（b）～图5（e）分别为基于位移协调、能量协调、位移协调与能量协调相结合、连续分布耦合理论建立的多尺度模型。分别命名为位移法模型、能量法模型、结合法模型和DC法模型。

图5 钢筋混凝土节点梁-实体单元多尺度连接模型

对于实体单元部分钢筋混凝土采用分离式建模处理，其建模原理是将混凝土和钢筋分别建模，混凝土采用C3D8R单元，钢筋采用T3D2单元，并利用嵌固作用将钢筋嵌入到混凝土中，实现钢筋与混凝土的良好结合。对于梁单元部分的钢筋混凝土采用整体式模型，其建模原理是在同一位置建立两个共节点B31梁单元，梁单元1为混凝土，梁单元2用箱形截面等效模拟钢筋，其建模方式如图6所示。混凝土采用塑性损伤本构模型，钢筋采用无屈服点的弹塑性加工硬化两折线模型。模型网格尺寸为50mm。

图6 箱形截面等效配筋

2.1 轴力工况分析

在分析模型的底部设置固定约束，在顶部设置加载点通过施加位移控制加载，竖向位移从0逐渐增长到10mm。各模型在塑性阶段混凝土单元的应力云图如图7所示，各模型在破坏阶段混凝土受压损伤云图如图8所示。通过观察加载过程中应力云图的变化可以发现，在轴力工况下，随着竖向位移的增大，各模型的单元应力状态逐渐由弹性进入塑性，其中构件的两端和中部最先进入塑性，混凝土在达到极限抗压强度后发生破坏，应力下降。在加载过程中，构件的中部逐渐形成损伤累积区并最终发生破坏。对比发现，除了位移法模型，其余各多尺度模型均能较好地反应构件在弹塑性阶段的应力和损伤分布，印证了基于位移协调理论构建的多尺度模型在仅受轴力作用时，多尺度界面存在侧向过约束。

图7 轴力工况下塑性阶段各模型混凝土压应力云图

图8 轴力工况下各模型混凝土受压损伤累积云图

提取各模型在加载过程中加载点的位移和对应的轴力，作出顶点位移-轴力曲线如图9所示。5组模型加载结束时顶点轴力分别为1510、1550、1610、1500、1530kN。相比精细化模型，各多尺度模型相对误差分别为2.65%、5.88%、0.2%、1.12%。

图9 顶点位移-轴力曲线

按照象限编号法则，由内到外，在各模型多尺度界面上选取16个特征单元，如图10所示。当顶点竖向位移加载到2.15mm时构件处于塑性阶段，分别提取各模型16个特征单元的Min.principle应力，统计精细化模型与多尺度模型在特征单元处的应力值并计算相对误差。观察发现，平均相对误差从大到小分别为位移法模型10.1%、能量法模型1.8%、结合法模型1.6%、DC法模型4.41%。位移法模型的13、14号特征单元相对误差最大，达到32.1%。能量法模型的1、2、3、4号特征单元相对误差最大，达到2.5%。结合法模型的1、2、3、4号特征单元相对误差最大，达到2.1%。DC法模型的1、2、3、4号特征单元相对误差最大，达到5.1%。在轴力工况下，位移法模型的相对误差率高于其他多尺度模型，印证了位移协调法在仅受轴力作用下存在侧向过约束局限性的分析。DC法模型的平均相对误差为4.4%，控制在5%之内，印证了使用连续分布耦合约束可以有效改善侧向过约束降低误差率。

图10 多尺度界面特征单元编号

2.2 弯矩工况分析

在分析模型的底部设置固定约束，在顶部设置加载点通过施加绕X轴的弯矩荷载控制加载，弯矩从0逐渐增长到40kN·m。

各模型在塑性阶段混凝土单元的应力云图如图11所示，各模型在破坏阶段混凝土受拉损伤云图如图12所示。通过观察加载过程中应力云图的变化可以发现，混凝土最大拉、压应力分别沿受拉、受压侧边缘分布。随着不断加载，构件在受拉侧，沿柱身不断形成损伤累积区并最终发生破坏。对比各模型可以发现，能量法模型相比精细化模型应力分布和损伤累积差别都较大。位移法模型和DC发模型与精细化模型在应力分布上吻合度比较高，但在多尺度界面附近的损伤累积有所差别。结合法模型与精细化模型吻合度极高，能很好地反应构件的应力分布及损伤累积。

图11 弯矩工况下塑性阶段各模型混凝土轴向应力云图

图12 弯矩工况下各模型混凝土受拉损伤累积云图

提取各模型在加载过程中加载点绕X轴的弯矩和对应绕X轴的转角，作出顶点转角-弯矩曲线如图13所示。5组模型加载结束时顶点转角分别为0.00718、0.00923、0.01768、0.00676、0.00493。相比精细化模型，各多尺度模型相对误差分别为28.55%、91.55%、5.85%、31.94%。

图13 顶点转角-弯矩曲线

当顶点弯矩加载到36.8kN·m时构件处于塑性阶段，分别提取各模型16个特征单元的S33应力，统计精细化模型与多尺度模型在特征单元处的应力值并计算相对误差，观察发现，平均相对误差从大到小分别为能量法模型78.5%、位移法模型27.7%、DC法模型6.6%、结合法模型0.7%。能量法模型的3、4号特征单元相对误差最大，达到239.2%。位移法模型的1、2、3、4号特征单元相对误差最大，达到50%。DC法模型的1、2、3、4号特征单元相对误差最大，达到7.1%。结合法模型的9、10、11、12号特征单元相对误差最大，达到2.3%。在弯矩工况下，能量法模型的相对误差率远大于其他多尺度模型，印证了能量协调法在弯矩作用下局限性的分析，既当构件进入塑性阶段后，弯矩产生的正应力呈非线性分布，基于能量协调的多点约束方程不再适用。

2.3 剪力工况分析

在分析模型的底部设置固定约束，在顶部设置加载点通过施加沿X轴的剪力荷载控制加载，剪力从0逐渐增长到20kN。

各模型在塑性阶段混凝土单元的应力云图如图14所示，各模型在破坏阶段混凝土受拉损伤云图如图15所示。通过观察加载过程中应力云图的变化可以发现，构件在端部剪力的作用下，底部混凝土首先进入塑性，并在拉压两侧产生应力集中。随着不断加载，构件底部发生破坏，塑性变形快速增大。对比各模型可以发现，位移法模型与结合法模型的应力分布、损伤累积计算结果均较贴合精细化模型所得结果。能量法模型在多尺度界面出现了明显应力突变，损伤积累也与精细化模型有较大差别。

图14 剪力工况下塑性阶段各模型混凝土轴向应力云图

提取各模型在加载过程中加载点沿X轴的剪力和对应沿X轴的位移，作出顶点位移-剪力曲线如图16所示。5组模型加载结束时顶点位移分别为1.92、1.89、1.89、1.87、1.89mm。相比精细化模型，各多尺度模型相对误差分别为1.22%、1.22%、2.25%。

图16 顶点位移-剪力曲线

当顶点剪力加载到20kN时构件处于塑性阶段，分别提取各模型16个特征单元的S33应力，统计精细化模型与多尺度模型在特征单元处的应力值并计算相对误差见表1，观察发现，在剪力荷载工况下，平均相对误差从大到小分别为能量法模型142.2%、位移法模型23.7%、DC法模型15.2%、结合法模型3.6%。能量法模型的1号特征单元相对误差最大，达到341.6%。位移法模型的1、2、3、4号特征单元相对误差最大，达到33.3%。DC法模型的1、2、3、4号特征单元相对误差最大，达到45.1%。结合法模型的6、7、8号特征单元相对误差最大，达到5.4%。在剪力工况下，能量法模型的相对误差率远大于其他多尺度模型，是因为在端部剪力的作用下构件底部产生弯矩，混凝土进入塑性后弯矩产生的正应力呈非线性分布，基于能量协调的多点约束方程不再适用。

表1 钢筋混凝土节点平均应力相对误差统计

3 结语

文中针对钢筋混凝土节点，基于位移协调、能量协调、位移协调与能量协调相结合以及基于连续分布耦合理论实现多尺度界面的连接并进行算例分析，所得具体研究成果如下：

（1） 基于位移协调实现多尺度模型界面连接时，除了在仅受轴力这种工况外，可以较为合理地反应构件的应力分布和损伤累积。

（2） 基于能量协调实现多尺度模型界面连接时，除了在轴力工况下得到的结果较为精确外，在剪力工况、弯矩工况下，由于材料进入塑性，弯矩产生的正应力分布公式不再成立，约束方程不再适用，塑性区域单元迭代运算逐渐出现误差累积。

（3） 基于位移协调和能量协调相结合实现多尺度模型界面连接时，在各工况下均可合理地反应构件应力分布和损伤累积，得到的结果也最接近精细化模型，但需要输入大量约束方程，建模前处理较为复杂。

（4） 基于连续分布耦合理论实现多尺度模型界面连接时，在各工况下均可合理地反应构件的应力分布和损伤累积，相比位移协调理论不仅可以有效改善侧向过约束，在轴力工况下得到较为精确的结果，还能避免输入大量的约束方程等繁琐的建模前处理。