0 引言

　　 木材外露的美学特性、优良的力学性能和材料的绿色可再生使其成为贯穿整个历史的最理想建筑材料之一。传统木结构建筑就是在历史长河中形成的一类具有典型特色的建筑形式，此类建筑往往具有极高的历史和人文价值，且很多被列为文物保护建筑或优秀历史建筑。与现代木结构有所区别的是传统木结构主要采用木-木连接节点来传递竖向和水平荷载。因此，准确把握木-木连接节点的受力行为对传统木结构的结构安全至关重要。

　　 斗拱节点是传统木结构中典型的木-木连接节点之一，其受力行为备受关注。目前，国内外针对斗拱节点的受力行为已经积累了一定的试验研究成果。Chen Z Y等 ^[1]和周乾等 ^[2,3]先后开展了斗拱节点竖向加载试验研究，发现斗拱节点的破坏模式主要包括木材的横纹受压屈服、开槽处木材横纹劈裂破坏和枋的顺纹受弯破坏等。袁建力等 ^[4]开展了竖向荷载作用下斗拱节点水平低周反复加载试验研究，发现泥道栱损伤严重，各组件之间的摩擦耗能是斗拱节点的主要耗能方式。程小武等 ^[5]开展了宋式带“昂”斗栱节点抗震性能试验研究和理论分析，发现斗栱节点的水平承载力与底部暗销锚固力和接触面摩擦力之和密切相关。Wu Y J等 ^[6]开展了偏心荷载作用下双斗拱水平低周反复加载试验研究，发现双斗拱相互连接的枋和剪力键的协同工作机制可显著提高结构的整体抗侧刚度和承载力。

　　 斗拱节点复杂的结构组成决定了各组件的损伤演化过程一般难以通过缩尺试验完整跟踪。为进一步掌握各组件的破坏机理，学者们还开展了斗拱节点的有限元模拟分析。陈志勇 ^[7]和谢启芳等 ^[8]先后开展了斗拱节点竖向加载有限元分析，发现有限元模型能合理表征竖向荷载作用下的斗拱节点非线性受力行为。袁建力等 ^[9]采用有限元软件开展了水平荷载作用下斗拱节点单调加载数值模拟分析，发现模拟曲线与低周反复加载试验的骨架曲线接近。薛建阳等 ^[10]开展了斗拱节点滞回性能的有限元分析，发现基于弹塑性本构模型的有限元模型会高估节点的初始刚度、极限承载力和延性比。主要原因是弹塑性模型难以表征反复荷载作用下木材的损伤演化规律。

　　 鉴于此，本文采用三维弹塑性损伤模型跟踪木材的损伤演化过程，建立斗拱节点精细化有限元分析模型。通过已有竖向单调加载试验结果 ^[6]校验有限元模型参数选取的合理性。在此基础上，采用校验后的有限元模型开展节点水平滞回性能模拟分析。基于模拟结果，对不同组件的损伤进行量化，根据应力发展过程分析损伤产生原因。

1 木材三维弹塑性损伤模型

　　 为准确表征木材的损伤演化规律，采用三维弹塑性损伤模型 ^[11]对其进行模拟。该模型主要包括弹塑性模型和损伤模型两部分。其中，弹塑性模型主要用于计算有效应力(无损伤截面材料的应力)、弹性应变和塑性应变。损伤模型主要用于追踪木材的损伤演化过程并计算柯西应力。

1.1 弹塑性模型

　　 有效应力与弹性应变之间的表达式为：

　　 σ¯=E0∶εeσ¯=E0∶εe (1)

　　 式中：σ¯¯¯σ¯为有效应力张量；E₀为弹性模量张量；ε^e为弹性应变张量。

　　 采用Hill屈服准则描述有效应力的屈服面形状 ^[12],其表达式为：

　　 f=σ¯¯¯T:A:σ¯¯¯−−−−−−−−√−σ¯iso=0A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢(H+G)−H−G−H(F+H)−F−G−F(F+G)2N2M2L⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥         (2)f=σ¯Τ:A:σ¯-σ¯iso=0A=[(Η+G)-Η-G-Η(F+Η)-F-G-F(F+G)2Ν2Μ2L]         (2)

　　 式中：σ¯isoσ¯iso为各向同性强化等效有效应力； f为屈服函数；F,G,H分别为与材料抗压、抗拉屈服强度和抗剪屈服强度对应的正应力分量的组合系数；L,M,N分别为与材料抗压、抗拉屈服强度和抗剪屈服强度对应的剪应力分量的组合系数 ^[11]。

　　 将描述有效应力σ¯¯¯σ¯与塑性应变ε^p增量之间关系的流动法则设置为关联流动法则，其表达式为：

　　 dεp=dλ∂f∂σ¯¯¯         (3)dεp=dλ∂f∂σ¯         (3)

　　 由弹塑性理论可知，加卸载过程应满足Kuhn-Tucker条件，即：

　　 f≤0,λ˙≥0,λ˙f=0         (4)f≤0,λ˙≥0,λ˙f=0         (4)

　　 式中λ˙λ˙为塑性流动因子率。

　　 由于Voce模型的导数连续性较好，且应用较为广泛，故而选取Voce模型作为木材的硬化模型，其表达式为：

　　 σ¯iso=σ¯0+Q(1−e−bε¯p)         (5)σ¯iso=σ¯0+Q(1-e-bε¯p)         (5)

　　 式中：σ¯0σ¯0为初始屈服应力；Q和b为曲线形状参数，一般通过材性试验确定；ε¯pε¯p为等效塑性应变，是反映塑性应变累积的变量，其数值并不会减小。

1.2 损伤模型

　　 损伤截面材料的应力-应变关系表达式为：

　　 σ=E(d)∶εe         (6)σ=E(d)∶εe         (6)

　　 式中：σ为柯西应力张量；E(d)为损伤材料的刚度张量，其Voigt矩阵形式的表达式为：

　　 E(d)=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢α2C11αβC21αγC31000αβC12β2C22βγC32000αγC13βγC23γ2C33000000βγC44000000αγC55000000αβC66⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥α=1−d1,β=1−d2,γ=1−d3         (7)E(d)=[α2C11αβC12αγC13000αβC21β2C22βγC23000αγC31βγC32γ2C33000000βγC44000000αγC55000000αβC66]α=1-d1,β=1-d2,γ=1-d3         (7)

　　 式中：α,β,γ分别为纵向、径向和切向的损伤变量，取受损截面材料残余刚度与原有完好刚度的比值；d_i为i方向的损伤变量，通常包括受拉损伤变量和受压损伤变量两部分，具体表达式为：

　　 di={dt,i (σ¯¯¯≥0)dc,i (σ¯¯¯<0)   (i=1,2,3)         (8)di={dt,i (σ¯≥0)dc,i (σ¯<0)   (i=1,2,3)         (8)

　　 式中下标t和c分别表示受拉损伤和受压损伤。

　　 为有效识别不同破坏模式，拉应力和剪应力作用下的破坏准则采用Sandhaas ^[13]提出的破坏准则：

　　 Fft,1=σ¯11f ft,1≥1         (σ¯11>0)Fft,2=(σ¯22f ft,2)2+(σ¯12f f12)2+(σ¯23f f23)2≥1 (σ¯22>0)Fft,3=(σ¯33f ft,3)2+(σ¯13f f13)2+(σ¯23f f23)2≥1 (σ¯33>0)         (9)Ft,1f=σ¯11f t,1f≥1         (σ¯11>0)Ft,2f=(σ¯22f t,2f)2+(σ¯12f 12f)2+(σ¯23f 23f)2≥1 (σ¯22>0)Ft,3f=(σ¯33f t,3f)2+(σ¯13f 13f)2+(σ¯23f 23f)2≥1 (σ¯33>0)         (9)

　　 式中：f^f_t为木材的抗拉强度；f^f₁₂,f^f₁₃和f^f₂₃为抗剪强度；σ¯11,⋯,σ¯23σ¯11,⋯,σ¯23为应力张量沿某一特定方向的分量；下标1,2和3分别代表木材的纵向、径向和切向。

　　 在顺纹压应力作用下的破坏准则选取为：

　　 F fc,1=−σ¯11f fc,1≥1   (σ¯11<0)         (10)F c,1f=-σ¯11f c,1f≥1   (σ¯11<0)         (10)

　　 式中f^f_{c, 1}为木材的极限抗压强度。

　　 采用指数型损伤演化模型控制木材受拉损伤演化过程，表达式为：

　　 dt,i=1−1F ft,iexp[(1−F ft,i)Lcσft,iεft,iGft,i] (i=1,2,3)dt,i=1-1F t,ifexp[(1-F t,if)Lcσt,ifεt,ifGft,i] (i=1,2,3) (11)

　　 式中：L_c为单元的特征长度；G_{ft, i}为木材的断裂能；σ^f_{t, i}为损伤起始应变ε^f_{t, i}对应的应力。

　　 木材顺纹受压损伤变量的表达式取为：

　　 dc,1=1−1F fc,1(1−A)−Aexp[B(1−F fc,1)]dc,1=1-1F c,1f(1-A)-Aexp[B(1-F c,1f)] (12)

　　 式中A和B为曲线形状参数，可以通过木材顺纹受压材性试验确定。

1.3 本构模型数值算法的实现

　　 基于应变增量法求解弹塑性损伤模型的数值解。其中，有效应力、弹性应变、塑性应变和一致切线刚度采用Simo等 ^[14]提出的最近投影点算法(CPPA)进行求解。其中，有效应力的更新步骤如下：1)在给定应变增量的前提下求解弹性试算应力。2)根据屈服准则判断试算应力点是否在屈服面以内。3)当试算应力在屈服面以内时，试算应力即为有效应力；当试算应力不在屈服面以内时，通过向试算屈服面做投影求解该荷载步的有效应力。

　　 在迭代计算完成后，通过破坏准则进行损伤判别。若满足破坏准则，则基于式(11)和(12)计算相应的损伤变量，通过式(6)更新柯西应力。采用粘滞正则化技术提升损伤模型的收敛速率。

　　 通过Fortran语言编写用户自定义子程序(UMAT)将弹塑性损伤模型嵌入有限元软件ABAQUS。

2 偏心受力斗拱节点有限元分析

2.1 偏心受力节点模型及试验介绍

　　 基于Wu Y J等 ^[6]开展的缩尺比例为1∶3.4的单榀偏心受力斗拱节点(图1)竖向单调加载试验和水平方向滞回性能试验开展有限元模拟分析。图1中各组件几何尺寸如表1所示。

　　 图1 偏心受力斗拱节点有限元模型   

　　  

　　 偏心受力斗拱节点各组件几何尺寸 表1

组件名称	几何尺寸/mm	组件名称	几何尺寸/mm
普拍枋	245×80	枋	55×80
栌斗	175×175	柱	ϕ165
华拱	55×100	剪力键	10×20
泥道栱	55×114	暗销	30×30
散斗	93×93

 

　　  

　　 斗拱节点采用非洲红花梨木制成。其中，红花梨木的弹性模量、屈服强度和极限强度根据材性试验确定。红花梨木的泊松比和断裂能等参数按照文献[13,15]确定，红花梨木的材性参数详见表2。

　　 试验时将偏心斗拱节点放置于试验台上，对下端的普拍枋进行固定。通过球铰约束对垂直于加载方向的枋进行约束。竖向单调加载试验采用千斤顶加载(加载速率控制为约2mm/min),通过位移计量测柱顶竖向位移。水平滞回性能试验采用量程为100kN的千斤顶在斗拱节点顶部的柱端施加45kN的竖向荷载。根据CUREE加载制度 ^[16]采用水平放置的伺服作动器施加水平位移，直至斗拱节点产生明显破坏或水平承载力降低到峰值荷载的80%以下，加载速率为5mm/min。水平滞回性能试验中主要量测作动器的水平荷载和加载头的水平位移。

　　 非洲红花梨木的材性参数 表2

参数		取值	参数		取值
弹性模量 /MPa	E11	15 207	泊松比	ν12	0.033
	E22	1 977		ν31	0.047
	E33	1 450		ν23	0.373
剪切模量 /MPa	G12	1141	抗拉强度 /MPa	ft, 1	101.99
	G13	912		ft, 2	3.53
	G23	274		ft, 3	2.83
抗压强度 /MPa	fc, 1	50.44	抗剪强度 /MPa	f12	10.32
	fc, 2	9.35		f13	11.71
	fc, 3	8.96		f23	2.00
断裂能 /(N/mm)	Gft, 1	12.0	硬化参数	Q	200.0
	Gft, 2	0.75		b	20.0
	Gft, 3	0.75	受压损伤 参数	A	0.50
极限抗压强度/MPa	ffc,1c,1f	69.20		B	1.15

 

　　 注：下标1,2和3分别代表纵向、径向和切向。

　　  

　　 图3 偏心受力斗拱节点竖向荷载-位移曲线 

　　  

　　 图4 偏心受力斗拱节点滞回曲线

　　  

2.2 斗拱节点有限元模型的建立

　　 根据Wu Y J等 ^[6]提供的几何尺寸，采用有限元软件ABAQUS建立偏心受力斗拱节点精细化有限元分析模型。

　　 为提高模型的计算收敛性，上部柱脚和底部普拍枋(图1中浅色组件)采用弹塑性模型进行模拟，其余组件(图1中深色组件)均采用损伤模型进行模拟。损伤模型的部分参数参考文献[11]进行设定。损伤模型和弹塑性模型的材性参数如表2所示。

　　 为考虑组件之间几何间隙的影响，在上部木柱和与其垂直的木枋两侧各设置了3mm间隙。此外，在与普拍枋连接的暗销两侧各设置了1mm间隙。有限元模型中单元的特征尺寸通常取为20mm, 木枋沿长度方向的单元尺寸取为100mm。

　　 为提高计算效率，有限元模型中仅对普拍枋与栌斗之间的暗销进行了精细化建模。暗销的几何尺寸为30mm×30mm×30mm, 各散斗底部与木枋之间的相互作用通过Tie约束命令进行考虑。其余组件之间的相互作用采用接触进行模拟，接触面之间的摩擦系数参考已有的同类型木-木节点研究成果确定 ^[11],取为0.20。

　　 对普拍枋底部施加固定约束，然后对垂直于加载方向的木枋端部施加滑动约束。竖向加载试验模拟中竖向荷载通过参考点以竖向位移的形式进行施加。滞回性能试验模拟中荷载分析步设置为两步：第一步以压强的形式对柱顶施加45kN恒定竖向荷载，第二步参考CUREE加载制度 ^[16]通过参考点对柱顶施加水平位移(图2)。为保证有限元模型的计算收敛性，仅对CUREE加载制度的主循环进行模拟。控制位移Δ取为40mm。

　　 图2 滞回性能试验加载 方式示意图  

　　  

2.3 计算结果及验证

　　 首先对建好的有限元模型进行竖向单调荷载作用分析，以考察模型的计算精度和参数设置的合理性。然后采用验证后的模型开展节点水平低周反复荷载作用分析，并与已有试验结果进行对比。

　　 偏心受力斗拱节点竖向荷载-位移曲线的试验和数值模拟结果对比如图3所示。由图3可知，有限元模型能合理反映斗拱节点在竖向荷载作用下的非线性受力行为，数值模拟曲线与试验曲线十分接近。

　　 根据竖向荷载-位移曲线可进一步计算节点的初始刚度k₀和竖向承载力P₀,如表3所示。由表3可知，有限元模型对节点初始刚度和竖向承载力具有较好的预测精度(模拟误差在11.4%以内)。这说明有限元模型参数的选取是合理的，可用于斗拱节点滞回性能的模拟分析。

　　 偏心受力斗拱节点滞回曲线的试验和数值模拟结果对比如图4所示。

　　 斗拱节点竖向单调加载试验和数值模拟结果 表3

斗拱节点	k0/(kN/mm)	δk0/%	P0/kN	δM0/%
试验结果	13.77	—	210.90	—
数值模拟结果	13.71	-0.4%	234.90	11.4%

 

　　 注：δ_k0,δ_M0分别表示k₀,P₀的模拟误差。

　　  

　　 图5 栌斗和与加载方向垂直的木枋的转动变形 ^[6]   

　　  

　　 图6 关键组件顺纹方向的损伤云图  

　　  

　　 图7 关键组件横纹方向的损伤云图(水平位移10mm)  

　　  

　　 图8 关键组件横纹方向的损伤云图(水平位移20mm) 

　　  

　　 由图4可知，由于充分考虑了木材顺纹方向和横纹方向的损伤演化过程，有限元模型能合理表征滞回曲线的捏拢效应、强度软化行为和刚度退化行为。

　　 需要说明的是，加载前期有限元模型的卸载路径和耗能能力与试验结果有所区别。斗拱各部件之间的水平荷载主要通过接触面的摩擦力传递，故而接触面之间的间隙和木材初始裂缝对卸载后反向加载路径具有重要影响。为保证有限元模型的计算收敛性，模型仅考虑了上部木柱和与其垂直的木枋之间的间隙以及暗销两侧的间隙，并未考虑其他组件之间间隙的影响以及试件初始裂缝的影响。后续研究可进一步考察不同组件之间的间隙和初始裂缝的影响，以便进一步提高有限元模型的预测精度。

　　 根据滞回曲线可进一步计算节点的初始刚度k和水平承载力P,以便于更为合理地评价有限元模型的计算精度。根据Wu Y J等 ^[6]的研究成果，节点的初始刚度基于第一个加载主循环确定，水平承载力根据滞回曲线的峰值荷载确定。表4给出了滞回性能试验均值和数值模拟结果。由表4可知，有限元模型对节点初始刚度和水平承载力的预测结果与试验结果吻合较好，模拟误差在12.5%以内。

　　 斗拱节点滞回性能试验和数值模拟结果 表4

斗拱节点	k/(kN/mm)	δk/%	P/kN	δM/%
试验结果	0.88	—	7.57	—
数值模拟结果	0.77	-12.5%	8.47	11.9%

 

　　 注：δ_k和δ_M分别表示k,P的模拟误差，k和P均为正反向加载的平均值。

　　  

　　 为进一步校验有限元模型的计算精度，图5给出了栌斗和与加载方向垂直的木枋的转动变形。由图5可知，有限元模型能合理反映栌斗和木枋在加载过程中的转动变形。

2.4 反复荷载作用下关键组件的损伤演化过程分析

　　 为研究各组件的损伤状况，图6和图7～10分别给出了关键组件在反复荷载作用下顺纹方向和横纹方向损伤云图。由图6可知，各组件顺纹方向的损伤并不显著，这与试验中观察到的现象基本一致。由图7～10可知，横纹方向的损伤主要集中于栌斗、泥道栱和上部与木柱接触的木枋。当水平位移加载到10mm时，栌斗和泥道栱开始出现损伤；当水平位移加载到20mm时，上部与木柱接触的木枋开始出现损伤；随着水平位移的不断增加，三个组件横纹方向的损伤不断扩展，并加速了节点承载力的丧失。从图7～10中还可以发现，栌斗横纹方向的损伤位置与试验中栌斗的开裂位置基本一致。综上，有限元模型能合理反映斗拱节点在反复荷载作用下的损伤演化过程。

　　 图9 关键组件横纹方向的损伤云图(水平位移40mm)  

　　  

　　 图10 关键组件横纹方向的损伤云图(水平位移60mm)  

　　  

　　 图11 4种轴压下斗拱节点滞回曲线  

　　  

3 参数分析

　　 为研究斗拱节点中木柱所受轴向压力对节点滞回性能的影响，采用验证后的有限元模型开展参数分析。保持节点几何尺寸、材料力学性能参数和边界条件不变，仅变换木柱所受轴向压力，从而揭示其对节点滞回性能的影响规律。

　　 分别考虑4种轴向压力(20,30,40,45kN)下斗拱节点相应的滞回曲线，如图11所示。由图11可知，斗拱节点的承载力和耗能能力随轴向压力的逐渐增加而增加。当轴向压力超过40kN后，节点的承载力增幅降低。主要原因是轴向压力增加后，组件之间的摩擦力会提高，从而提升节点的承载力。过大的轴向压力会显著增加栌斗承受的弯矩，使得该处木材横纹方向产生损伤，限制了节点承载力的提升。

4 结论

　　 (1)本文建立的有限元模型能合理表征偏心受力节点在低周反复荷载作用下的捏拢效应、强度软化行为和刚度退化行为，且该模型对节点初始刚度和水平承载力的预测精度较高，模拟误差在12.5%以内。

　　 (2)有限元模拟结果表明，有限元模型能合理反映斗拱节点关键部件在反复荷载作用下的转动变形和木材的损伤演化过程。

　　 (3)参数分析结果表明，斗拱节点水平方向的极限承载力和耗能能力随轴向压力的逐渐增加而有所提高。当上部木柱承受的轴向压力超过40kN后，节点的承载力增幅降低。