基于木材弹塑性损伤本构的古建木结构残损梁柱构件损伤非线性分析

张利朋，谢启芳，吴亚杰，刘伊津

(西安建筑科技大学 土木工程学院；结构工程与抗震教育部重点实验室，西安 710055)

古建筑木结构是中国重要的文化遗产。木梁和木柱组成古建筑木结构的承重骨架，并维持其整体稳定。然而，由于数百年来环境与荷载的共同作用，木梁和木柱发生了不同程度的残损，如整体老化、局部腐朽、虫蛀和开裂等，不仅导致其受力性能退化，也使古建筑木结构的整体力学性能受到不同程度的削弱。因此，研究残损梁柱构件的受力行为可为现存残损古建筑木结构的保护提供科学依据[1]。

近年来，学者们围绕残损梁柱构件的调查实测和受力性能分析开展了深入研究，并取得了一定的研究成果。在残损梁柱构件的调查实测方面，周乾[2]通过对故宫古建筑木结构柱子的勘察实测，指出木柱的典型残损问题主要为局部糟朽和开裂，并进行了成因分析；王晓丽[3]利用快速普查法调查了飞云楼的承重木柱和梁枋的残损情况，确定了各构件的残损等级，为古建筑木结构的保护及修缮提供了依据；乔冠峰[4]从承重木柱、承重梁枋、结构整体变形三个方面入手，对古建筑木结构的残损点评估界限进行了探讨，并对残损界限评定的适用条件提出了建议。在残损梁柱构件的受力分析方面，Guan等[5-6]基于弹塑性本构模型，采用单元生死技术研究了跨中开方孔木梁的受力性能和方孔角部的损伤演化；汤永红[7]基于ANSYS模拟了纵向开裂木柱的轴压性能，并将一阶屈曲特征位移向量的2%作为初始几何变形进行了纵向开裂木柱的非线性分析；Zhang等[8]通过人工预制贯通缝的方式模拟了木柱的干缩裂缝，并通过轴心和偏心受压试验研究了失效模式和残余极限承载力的变化规律；朱忠漫[9]基于ABAQUS对纵向开裂木柱的受力性能进行了参数分析，首先对木柱开展线性屈曲模拟，进一步采用弧长法进行非线性屈曲分析，将第一、二阶位移特征向量按照一定比例加权后，作为初始几何缺陷输入有限元模型；Mosallam[10]对人工预制残损木梁进行了抗弯试验，分析了其力学性能退化规律；陈立涛[11]针对实际结构中存在的局部腐朽和虫蛀残损，通过预制局部缺口受弯木梁的受弯试验和数值模拟，研究了残损对木梁受力性能的影响规律；路鹏[12]采用材料性能折减法在ABAQUS中建立了老化木柱的有限元模型，并分析了不同老化程度的影响；Li等[13]通过人工开槽方式开展了柱脚糟朽木构架的试验，并采用有限单元法对其进行了数值模拟。谢启芳等[14]基于ABAQUS模拟了局部带缺口木柱的轴压性能，通过分析缺口的大小和位置，考察了初始缺口缺陷对木柱承载力和非线性的影响。

可见，有关古建筑木结构梁柱构件的研究主要集中在开裂残损方面，而关于木梁、木柱局部糟朽和虫蛀影响下受力性能方面的研究较少，并且该方面的研究大多基于试验方法开展，仅少数研究采用了有限单元法，且其采用的本构模型均为ABAQUS自带的弹塑性模型，导致其结果无法合理反映局部残损部位与完好部位交界处的材料损伤行为。通常情况下，局部残损梁柱构件在局部交界面的破坏涉及损伤的演化、应力的软化，而这是ABAQUS自带弹塑性模型所无法直接反映的，所建有限元模型便无法较好地反映整个梁柱构件的受力非线性和损伤演化过程。笔者通过建立木材的弹塑性损伤本构模型，分析带局部缺陷梁柱构件的损伤演化、发展过程，并基于课题组试验进行验证。研究成果可为残损梁柱构件的受力分析及性能评价提供借鉴。

1 木材弹塑性损伤本构模型

1.1 理论模型

σ=Cd(d):(ε-εp)

(1)

式中：Cd(d)为木材的弹性损伤刚度张量；ε、εp和σ分别为木材的总应变张量、塑性应变张量和应力张量；d为木材的损伤变量集，如式(2)所示。

d={d1t,d1c,d2t,d2c,d3t,d3c,d12s,d13s,,drs}

(2)

式中：dit和dic(i=1,2,3)为木材在各主轴方向的受拉和受压损伤变量；d12s、d13s和d23s分别为各剪切面内的剪切损伤变量。

由式(1)可知，弹塑性损伤本构模型的关键在于确定塑性应变和损伤变量的演化规律。

1.1.2 塑性部分 木材的塑性屈服函数为

(3)

木材强度准则选用Hoffman准则[15]。

(4)

式中：Pα和pα可展开为式(5)和式(6)。

(5)

pα=[α11α22α330 0 0]T

(6)

式中：αij为与木材单轴受拉和受压强度有关的参数，其表达形式和确定方法参见文献[16]。

(7)

采用相关性流动法则，则木材的塑性应变演化方程可表达为

(8)

1.1.3 损伤部分 木材弹塑性损伤本构模型的损伤部分需给出损伤演化起始条件和损伤变量演化方程。

参考复合材料的损伤建模理论，木材损伤准则可以考虑失效机制并表达为[17]

fI(φI,rI)=φI-rI≤0

(I=1t,1c,2t,2c,2v,3t,3c,3v)

(9)

式中：φI为与木材损伤机制有关的函数，采用Sandhaas[18]提出的木材损伤准则，如式(10)所示。

(10)

式(9)中，rI为损伤阈值函数，其表达式为

(11)

本文模型采用的损伤演化方程为[15-16]

J={1t;1c;2t;2c;3t;3c;12s;13s;rs}

(12)

gJ=GJ/Lc

(13)

式中：L为单元特征长度，由ABAQUS主程序计算后赋值给变量CELENT，具体计算方法参见ABAQUS帮助文档[19]；GL为沿单元特征长度的应变能释放，N/mm。

(14)

对于增量过程，式(14)可显化表达为

(15)

1.2 本构模型的数值算法

建立的弹塑性损伤本构模型由隐式的后退欧拉算法进行求解，该方法通过在每个增量步结束时更新内变量以避免求解结果从屈服面漂移。

有效应力空间中的积分算法可以表达为

εn+1=εn+Δεn+1

(16)

1.2.1 弹性预测 根据式(16)，可得

(17)

(18)

增量步开始时，总应变增量全部被视作弹性应变，据此计算试探应力，如果试探应力在屈服面内，则说明试算过程为弹性过程，有效应力就等于试探应力，否则应进行塑性修正。

1.2.2 塑性修正 经过初始弹性预测，可得塑性应变增量

(19)

将式(19)代入式(16)，可得

(20)

通过求解式(20)方程组，可得

(21)

式中:

(22)

将式(20)代入式(19)，可得

(23)

因此，可求得塑性乘数增量

(24)

进一步可得有效应力张量、塑性应变张量和等效塑性应变的表达式。

(25)

1.2.3 损伤修正 基于塑性修正过程求得的有效应力状态，可计算损伤内变量，进而对有效应力进行折减，即可求得名义应力。由于该过程不涉及迭代计算，计算过程简单，不再展开。

1.3 本构模型验证

通过将所建立的木材本构模型基于ABAQUS的材料二次开发接口UMAT进行程序实现，可借助其计算模块实现二次开发。文献[16]给出了详细的本构算法和程序开发过程，且已基于多种加载工况下的材性试验(顺纹及横纹方向的单轴受拉和受压、顺纹及横纹重复受压等[20-21])进行了模型验证，模拟结果与试验结果吻合较好。基于三点受弯木梁在构件层面对所建本构模型和所开发的本构程序进行模型验证。

三点受弯木梁模型采用的材料和试验数据取自Khennane等[22]的试验，木材为云杉，木梁横截面尺寸为60 mm×60 mm，长1 600 mm，如图1所示。

图1 三点受弯木梁尺寸及试验加载装置[22]Fig.1 Size and test setup of three-point bending timber

为了防止木梁局部压溃，试验过程中，在施加集中力处放置了60 mm × 60 mm钢垫块，采用位移控制加载，直到试件破坏。试验结束时，试件底部跨中发生了脆性受拉断裂。

1.3.1 几何模型、网格划分、接触设置 三点受弯木梁有限元模型采用与试验相同的尺寸，梁顶面的加载垫块和梁底面两端支座处采用截面尺寸为60 mm×60 mm，厚20 mm的垫块进行模拟。选取C3D8R单元类型，网格长度设置为10 mm。确保垫块与梁的网格节点重合以提高收敛性。如果在钢垫块与梁顶部接触面之间采用绑定约束(Tie)，则会在钢垫块边缘与梁顶接触处引起应力集中，且损伤将首先在此产生，并进一步引起收敛性问题，而采用考虑摩擦系数的方法则可以避免这一问题[23]。因此，垫块与梁之间的接触面相互作用采用法向硬接触和切向摩擦系数法。考虑到垫块与梁顶面之间发生相对位移的可能性很小，取摩擦系数为0.45。

表1 本构模型参数确定Table 1 Determination of wood properties

1.3.3 边界条件及加载方式 有限元模型的边界条件按试验情况设为简支，左端约束梁的水平位移和竖向位移，右端仅约束竖向位移。采用位移加载模式加载至计算过程终止。

1.3.4 模拟结果分析 由三点受弯木梁有限元模拟结果的失效模式与试验结果对比(顺纹受拉损伤变量SDV20)，如图2所示。由图2可见，顺纹受拉损伤分布与试验结果较为一致，均为梁跨中截面(图3中的A点)发生。此时，SDV20的最大值为0.81，计算过程中断。通过所建立的木材弹塑性损伤本构模型可以较好地预测三点受弯木梁的脆性断裂失效。

图2 三点受弯木梁有限元模拟与试验结果对比Fig.2 Comparison of finite element simulation and experimental results of three-point bending wood

图3 三点受弯木梁底部损伤单元的顺纹受拉应力应变曲线与应力云图分析Fig.3 Analysis of stress-strain curve and stress nephogram of damage element along grain at the bottom of three-point

1.4 黏性调整对收敛性的影响

以简单受拉试件为例，考察黏性规则化对脆性断裂损伤引起收敛困难的改善，主要讨论两个问题：一是有无黏性调整对收敛性的影响；二是黏性系数取值范围的问题。所选试件模型如图4所示，构件尺寸为210 mm × 10 mm × 10 mm，所采用的主要材料参数选自文献[16]，顺纹抗拉强度为24 MPa，顺纹弹性模量为11 000 MPa，黏性规则化系数为0.000 1。具体操作方法为：沿构件长度方向选定部分单元为初始损伤单元，并通过人工降低损伤单元材料性能(强度设置为23 MPa)的方式形成断裂带，实现初始损伤位置的预设，从而排除软化段对网格敏感性的影响而单独研究黏性调整对收敛性的影响。

图4 黏性调整对单向受拉木构件收敛性的影响Fig.4 Influence of viscosity adjustment on convergence of unidirectional tension wood

由图4可知，未进行黏性调整时，模型仅能预测木材顺纹受拉峰值强度，而无法反映其断裂过程；黏性调整后，则可以较好地预测断裂路径。可见，黏性调整有利于在使用弹塑性损伤本构模型程序计算构件损伤破坏过程时提高其收敛性。

图5 黏性规则化系数的优化取值分析Fig.5 Optimization value analysis of viscosity

2 局部残损梁柱构件的数值模型

2.1 局部残损木梁的数值模型

2.1.1 残损木梁几何模型 建立的残损木梁有限元模型的几何模型选自陈立涛进行的残损木梁试验[11]，其几何尺寸如图6所示。

图6 残损木梁几何模型Fig.6 Geometric model of damaged wooden

图7 底部缺口残损木梁的破坏模式Fig.7 Failure mode of timber beam with bottom

2.1.2 残损木梁有限元模型 根据上述带缺口木梁的几何条件建立有限元模型，其边界条件、加载位置及支座布置、材料模型、分析步设置等与试验基本等效。有限元模型的边界条件按试验情况设为简支，左端约束梁的水平位移和竖向位移，右端仅约束竖向位移。采用位移加载模式加载至计算过程终止。区别之处仅在于网格划分模块，残损木梁有限元模型在缺口附近的网格尺寸划分较小。

2.1.3 残损木梁的破坏模式 底部缺口残损木梁有限元模型的顺纹剪切损伤(SDV24)发展和演化情况如图8所示。由图8可见，使用所开发的木材弹塑性损伤本构模型可较好地预测其损伤的发展和演化过程，但由于使用隐式材料程序(UMAT)，该模型不具有显式有限元程序的网格删除功能，因此，裂缝的开展情况只能通过损伤进行分析而无法直观地在有限元模型中表达。

图8 底部缺口残损木梁有限元模型的剪切损伤发展和演化Fig.8 Shear damage development and evolution of fintte element model of damaged wood beam with bottom

图9 残损木梁荷载位移曲线Fig.9 Load-displacement curves of damaged wooden

由图9可见，建立的木材弹塑性损伤本构模型可以较好地预测带有局部缺口的四点受弯木梁的承载能力，但由于模型假定在材料各个方向达到屈服前均为线性，导致构件模拟的刚度偏大。

2.2 局部残损木柱的数值模型

建立的残损木柱有限元模型所采用的几何模型选自王玄[24]进行的残损木柱试验。边界条件、加载位置及支座布置、材料模型、分析步设置等与试验基本等效。采用位移加载模式加载至计算过程终止。带缺口木柱有限元模型的顺纹受压损伤(SDV21)发展和演化情况如图10所示。由图10可见，使用建立的木材弹塑性损伤本构模型可较好地预测其损伤的发展和演化过程，但由于所使用的是隐式材料程序，不具有网格删除功能，因此，局部压屈的发展情况只能通过损伤进行分析而无法直观地在有限元模型中进行表达。

图10 局部缺口木柱破坏模式的有限元与试验结果对比Fig.10 Comparison of finite element and experimental results of failure modes of wood columns with local

残损木柱竖向荷载位移曲线的有限元分析结果与试验结果的对比情况如图11所示。由图11可见，建立的弹塑性损伤本构模型可以较好地预测其承载能力，且在局部损伤出现之前与试验结果吻合较好。但在木柱损伤逐渐趋近于缺口附近时，有限元结果出现较为明显的软化段，总体带有脆性失效特征，因此无法反映出试验结果的延性失效过程。

图11 残损木柱竖向荷载位移曲线Fig.11 Vertical load-displacement curve of

3 结论

1)在连续介质力学和损伤力学框架内建立了木材的弹塑性损伤本构模型，可用于模拟三点受弯木梁的受弯损伤模式及分布规律。

2)所开发的弹塑性损伤本构模型可用于反映局部残损梁柱构件在残损界面与完好界面处的损伤行为，可较好地反映残损梁柱构件受力过程中的非线性行为和损伤演化。

3)从文中梁柱构件的3个模拟实例来看，建立的木材弹塑性损伤本构模型在应用于残损梁柱构件的数值模拟时，无法完全使分析过程收敛到最终，这是由于建立的木材弹塑性损伤本构模型并未完全解决塑性与损伤的耦合问题：文中采用的强度准则和屈服准则都是有效应力的函数，二者在本质上相同，即损伤准则某种程度上也是强度准则，强度准则某种程度上也是损伤准则，都是对线性区域的限定，但由于其表达式不同，所限定的区域并不相同，从而出现了塑性与损伤的“打架”问题。该问题需在后续研究中进一步解决。