一种考虑木材拉压不等强的非线性统一本构模型及试验验证

A Nonlinear Unified Constitutive Model Considering the Unequal Tensile and Compressive Strength of Wood and Its Experimental Verification

摘要: 木材作为一种天然的非均质材料，在荷载作用下同一纹理方向会表现出拉、压性能不对称的现象。为了更好地描述拉压性能不同的木材非线性力学行为，在一种用于描述木材受压的双参数本构模型的基础之上，通过引入Drucker-Prager塑性理论的静水压概念，建立了能够描述拉压性能不同的木材非线性统一本构模型。实验验证结果表明：模型描述的应力–应变曲线与实验结果吻合较好，运用本文中所建立的本构模型能够很好地描述木材在载荷作用下的非线性力学行为，而且考虑了木材同一纹理方向拉伸和压缩性能的差异。

Abstract: As a natural heterogeneous material, wood exhibits asymmetric tensile and compressive properties in the same texture direction under load. In order to better describe the nonlinear mechanical behavior of wood with different tensile and compressive properties, based on a two parameter constitutive model used to describe wood under compression, a nonlinear unified constitutive model of wood that can describe different tensile and compressive properties was established by introducing the concept of hydrostatic pressure from Drucker Prager plasticity theory. The experimental verification results show that the stress-strain curve described by the model is in good agreement with the experimental results. The constitutive model established in this paper can well describe the nonlinear mechanical behavior of wood under load, and consider the differences in tensile and compressive properties of wood in the same texture direction.

1. 引言

木材作为一种建筑材料，被广泛应用于宫殿、庙宇、民居等传统建筑领域[1] [2]。近年来，随着现代工程技术的不断提高，木材在高层木结构建筑[3]-[5]和复杂木结构桥梁[6] [7]等新型木结构领域展现出巨大潜力。随着现代木结构的高度和跨度不断突破，深入理解木材的力学性质对木结构的安全性、经济性和可靠性变得愈加重要。为了能够模拟木结构及其构件在荷载作用下的力学响应，需要建立一套能够准确反映木材应力–应变关系的本构模型。

大量的试验结果表明，木材在偏轴载荷作用下具有显著的物理非线性力学特性。多年来，国内外学者针对木材非线性力学行为的研究取得了许多成就，已建立起木材在不同受力方向和不同受力方式下的本构模型，并在木结构数值分析中进行了初步应用。如：Yoshihara [8]对贝壳杉和连香树沿晶粒不同取向的试样进行了压缩试验，将获得的离轴应力–应变关系回归到Ramberg-Osgood方程中，在通过Hankinson公式建立了模型参数与偏轴角度之间的关系后，成功地根据正交各向异性轴上的应力–应变关系预测了斜方向的应力–应变关系。杨娜等[9]在Sun和Chen [10]的单参数塑性本构模型基础上合理假设木材是横观各向同性材料，提出一种包含双参数的弹塑性本构模型，准确地反映了木材受压的非线性应力–应变关系。此外，Oudjene [11]等基于Hill强度准则和各向异性弹塑性理论，建立了能够反映木材横纹受压致密化效应的本构模型，并将其在商用ABAQUS软件中成功实现。

不应忽略的是，作为一种天然的生物质材料，木材具有复杂的组织构造。这不仅使其在纤维方向(顺纹)和垂直于纤维方向(横纹)的力学性能表现出显著差异。更为复杂的是，同一纹理方向也会表现出类似于“Bauchinger效应”的拉、压性能不对称的现象[12]。尽管上述的非线性本构模型能够合理表征木材不同方向的非线性力学行为，但不能反映木材同一纹理方向拉压不等强的特性。

近年来，研究者们开始认识到这种拉压不同性质对木材力学性能影响的重要性，并致力于建立可反映木材拉压性能不同的统一本构模型。Xu等[13]将Hoffman破坏准则与Hill准则[14]相结合来控制木材的塑性屈服和损伤演化，建立了三维非线性有限元模型。然而，Hoffman准则[15]虽能体现出材料在同一方向上不同的拉压响应，但缺乏系数限制条件，从而导致强度包络线不闭合。陈志勇等[16]和张慎等[17]基于弹塑性理论和连续介质损伤力学先后(分别)建立了木材弹塑性本构模型，均考虑了木材同一纹理方向的拉压力学性能差异，但是却存在方程过于复杂、待测参数多等缺点，不利于实际应用。

鉴于此，本文从木材偏轴拉伸和压缩的宏观试验出发，旨在建立描述拉压性能不同的木材非线性行为的唯象理论，以期进一步拓展非线性本构理论在工程实际中的应用。具体思路为：首先，在杨娜的双参数塑性模型基础上，引入Drucker-Prager [18]塑性理论的静水压概念，并据此深入考虑了木材塑性变形对载荷方向的依赖性，建立能够描述木材拉压性能不同的非线性本构方程。其次，对木材的离轴试样进行了压缩和拉伸试验，并基于试验结果，获取模型的各向异性系数和加载方向效应系数。最后，运用本文所建立的本构理论，预测木材在偏轴拉伸和拉压缩荷作用下的非线性应力–应变响应，并与试验结果进行对比，证验本文所提出的理论模型的正确性和有效性。

2. 本构模型

增量弹塑性本构模型可表示为：

$\text{d}{\epsilon }_{ij}=S_{ijkl}^{ep}\text{d}{\sigma }_{kl}$ (1)

式中：$S_{ijkl}^{ep}$ 为弹塑性柔度张量。

材料加载过程中，总的应变增量$\text{d}{\epsilon }_{ij}$ 可以看作是由弹性增量$\text{d}{\epsilon }_{ij}^e$ 和塑性增量$\text{d}{\epsilon }_{ij}^p$ 两部分组成，如下式所示

$\text{d}{\epsilon }_{ij}=\text{d}{\epsilon }_{ij}^e+\text{d}{\epsilon }_{ij}^p$ (2)

又可表示为

$\text{d}{\epsilon }_{ij}=\left(S_{ijkl}^e+S_{ijkl}^p\right)\text{d}{\sigma }_{kl}$ (3)

式中：$S_{ijkl}^e$ 与$S_{ijkl}^p$ 分别为弹性柔度和塑性柔度。

根据广义胡克定律，弹性应变增量$\text{d}{\epsilon }_{ij}^e$ 与应力增量$\text{d}{\sigma }_{kl}$ 的关系为：

$\text{d}{\epsilon }_{ij}^e=S_{ijkl}^e\text{d}{\sigma }_{kl}$ (4)

平面应力状态下弹性本构方程的增量形式表示为

$\left\{\begin{array}{c}\text{d}{\epsilon }_{11}^e\\ \text{d}{\epsilon }_{22}^e\\ \text{d}{\epsilon }_{12}^e\end{array}\right\}=\left[\begin{array}{ccc}1/E_1& -{\nu }_{12}/E_1& 0\\ & 1/E_2& 0\\ sym& & 1/G_{12}\end{array}\right]\left\{\begin{array}{c}\text{d}{\sigma }_{11}\\ \text{d}{\sigma }_{22}\\ \text{d}{\sigma }_{12}\end{array}\right\}$ (5)

式中：1，2，3方向分别对应木材的顺纹方向(L)，径向(R)和弦向(T)。

杨娜等在一种用于纤维复合材料的单参数本构模型的基础之上，通过考虑木材在顺纹方向产生塑性应变并假设木材是横观各向同性材料，提出了一个适用于描述木材受压性能的双参数弹塑性本构模型。他们假设平面应力状态下的等效应力为：

${\tilde{\sigma }}^{2D}=\sqrt{\frac{3}{2}\left({\sigma }_{11}^2+a_{22}{\sigma }_{22}^2-{\sigma }_{11}{\sigma }_{22}+2a_{66}{\sigma }_{12}^2\right)}$ (6)

式中${\sigma }_{ij}\left(i,j=1,2,3\right)$ 指材料主方向的各应力分量，$a_{22},a_{66}$ 为各向异性系数。

值得说明的是，由于式(6)中不包含静水压力项，因此不能反映材料拉压不等强的性质。本研究中，采用与Drucker-Prager [18]屈服条件类比的方法，在式(6)的基础上引入静水压力项反映木材拉压性能差异。平面应力状态下的等效应力为

$\begin{array}{c}\overline{\sigma }=\sqrt{\frac{3}{2}\left({\sigma }_{11}^2+a_{22}{\sigma }_{22}^2-{\sigma }_{11}{\sigma }_{22}+2a_{66}{\sigma }_{12}^2\right)}+b\left({\sigma }_{11}+{\sigma }_{22}\right)\\ ={\sigma }_{dev}+b\left({\sigma }_{11}+{\sigma }_{22}\right)\end{array}$ (7)

式中：b为加载方向效应系数。

塑性应变增量表示为

$\text{d}{\epsilon }_{ij}^p=\frac{\partial \overline{\sigma }}{\partial {\sigma }_{ij}}\text{d}{\overline{\epsilon }}^p$(8)

式中：$\text{d}{\overline{\epsilon }}^p$为等效塑性应变增量

由式(7)，(8)可得

$\left\{\begin{array}{c}\text{d}{\epsilon }_{11}^p\\ \text{d}{\epsilon }_{22}^p\\ \text{d}{\epsilon }_{12}^p\end{array}\right\}=\left\{\begin{array}{c}\frac{3\left(2{\sigma }_{11}-{\sigma }_{22}\right)}{4{\sigma }_{dev}}+b\\ \frac{3\left(2a_{22}{\sigma }_{22}-{\sigma }_{11}\right)}{4{\sigma }_{dev}}+b\\ \frac{3a_{66}{\sigma }_{12}}{{\sigma }_{dev}}\end{array}\right\}\text{d}{\overline{\epsilon }}^p$(9)

Sun和Chen [10]将等效应力–等效塑性应变关系近似为幂律，允许材料加载的初始阶段便开始产生塑性应变，但由于一开始施加的荷载很小，材料的塑性应变处在一个很低的水平。因此，假设等效应力与等效塑性应变的关系为幂函数

${\overline{\epsilon }}^p=A{\overline{\sigma }}^n$ (10)

其中，A、n为表征材料塑性的相关参数，可通过试验取得。则硬化系数可为

$H_p=\frac{\text{d}\overline{\sigma }}{\text{d}{\overline{\epsilon }}^p}=\frac{1}{An{\overline{\sigma }}^{n-1}}$(11)

代入式(9)得

$\left\{\begin{array}{c}\text{d}{\epsilon }_{11}^p\\ \text{d}{\epsilon }_{22}^p\\ \text{d}{\epsilon }_{12}^p\end{array}\right\}=\left\{\begin{array}{c}\frac{3\left(2{\sigma }_{11}-{\sigma }_{22}\right)}{4{\sigma }_{dev}}+b\\ \frac{3\left(2a_{22}{\sigma }_{22}-{\sigma }_{11}\right)}{4{\sigma }_{dev}}+b\\ \frac{3a_{66}{\sigma }_{12}}{{\sigma }_{dev}}\end{array}\right\}\frac{\text{d}\overline{\sigma }}{H_P}$(12)

将式(7)代入式(12)得塑性柔度为

$S_{ijkl}^p=\frac{1}{H_P}\frac{\partial \overline{\sigma }}{\partial {\sigma }_{ij}}\frac{\partial \overline{\sigma }}{\partial {\sigma }_{kl}}$ (13)

结合式(1)、(5)、(13)可知，考虑拉压不同强的弹塑性应力–应变关系的完整表达式为

$\left\{\begin{array}{c}\text{d}{\epsilon }_{11}\\ \text{d}{\epsilon }_{22}\\ \text{d}{\epsilon }_{12}\end{array}\right\}={\left[\begin{array}{ccc}{S}_{11}& S_{12}& S_{16}\\ & S_{22}& S_{26}\\ sym& & S_{66}\end{array}\right]}^{ep}\left\{\begin{array}{c}\text{d}{\sigma }_{11}\\ \text{d}{\sigma }_{22}\\ \text{d}{\sigma }_{12}\end{array}\right\}$ (14)

其中：

$\{\begin{array}{l}{S}_{11}=\frac{1}{E_1}+\frac{1}{H_P}{\left(\frac{3\left(2{\sigma }_{11}-{\sigma }_{22}\right)}{4{\tilde{\sigma }}^{2D}}+b\right)}^2\\ S_{12}=\frac{-{\nu }_{12}}{E_1}+\frac{1}{H_P}\left(\frac{3\left(2{\sigma }_{11}-{\sigma }_{22}\right)}{4{\sigma }_{dev}}+b\right)\left(\frac{3\left(2a_{22}{\sigma }_{22}-{\sigma }_{11}\right)}{4{\sigma }_{dev}}+b\right)\\ S_{16}=\frac{3a_{66}{\sigma }_{12}}{H_P{\sigma }_{dev}}\left(\frac{3\left(2a_{22}{\sigma }_{22}-{\sigma }_{11}\right)}{4{\sigma }_{dev}}+b\right)\\ S_{22}=\frac{1}{E_2}+\frac{1}{H_P}{\left(\frac{3\left(2a_{22}{\sigma }_{22}-{\sigma }_{11}\right)}{4{\sigma }_{dev}}+b\right)}^2\\ S_{26}=\frac{3a_{66}{\sigma }_{12}}{H_P{\sigma }_{dev}}\left(\frac{3\left(2a_{22}{\sigma }_{22}-{\sigma }_{11}\right)}{4{\sigma }_{dev}}+b\right)\\ S_{66}=\frac{1}{G_{12}}+\frac{1}{H_P}{\left(\frac{3a_{66}{\sigma }_{12}}{{\sigma }_{dev}}\right)}^2\end{array}$ (15)

3. 木材偏轴实验

3.1. 试验材料

选取白蜡作为本次试验的对象，其密度为0.59 g/cm³。为了确保试件表面纹理清晰、规整，所有受拉及受压试件均取自同一根原木的径切板。依据《木结构试验方法标准》[19]，棱柱体受压试件的高度是40 mm，横截面面积20 mm × 20 mm。受拉试件与Yoshihara和Ohta [20]所采用的试件相类似，试件的总长度是150 mm，中部的直线部分是40 mm，厚度为4 mm。为了防止试件两端被夹具夹碎，试件的端部用环氧树脂粘贴了白蜡薄片。根据荷载作用方向与木材纤维方向夹角的不同，分别制作了偏轴角度为0˚、15˚、30˚、45˚、60˚、75˚、90˚的受拉试件(图1)和受压试件(图2)。受拉及受压试件纤维角的变化始终发生在LR面内，荷载作用方向始终与弦向(T)垂直。在正式加载之前，需进行含水率平衡处理，使其含水率达到12% ± 1%。

Figure 1. Tension specimens with different off-axis angles in ash LR principal plane

图1. 白蜡LR主平面内不同偏轴角度的受拉试件

Figure 2. Compression specimens with different off-axis angles in ash LR principal plane

图2. 白蜡LR主平面内不同偏轴角度的受压试件

3.2. 试验方法

试验采用位移控制方式在100 kN微机控制电液伺服压力试验机上进行，受压加载速率为1 mm/min，受拉试件以0.5 mm/min匀速加载。应变测量通过在试件前后两面均粘贴应变片获取，试件受压时采用BE120-3AA型号(电阻值120 Ω栅格尺寸2.8 mm × 2 mm)应变片，受拉应变片型号为BE120-20AA (电阻值120 Ω栅格尺寸20 mm × 4.0 mm)。不同试件应变片的粘贴方式如图3所示。

(a) 受压试件 (b) 受拉试件

Figure 3. Compressive/tensile specimen dimensions and strain gauge placement

图3. 受压/拉试件尺寸和应变片贴法

假设作用在试件上的荷载为，荷载作用方向为x，则试件的平均应力为

${\sigma }_x=\frac{F}{A}$ (16)

式中：A为试件的横截面面积。考虑到加工带来的误差，A使用加工后的实测值。

3.3. 模型用于偏轴实验

Sun等[10]指出材料等效应力与等效塑性应变之间的关系是固定不变的，可通过偏轴试样的拉伸/压缩试验获得。考虑在偏轴载荷作用下建立等效应力–等效塑性应变关系，假设x轴为单向加载方向，1轴对应木材的顺纹方向，2轴对应木材的径向，$\theta$ 指荷载方向与顺纹方向的夹角。

根据应力坐标转换公式，材料主轴方向(1、2方向)应力分量与载荷应力${\sigma }_x$ 的关系可表示为

$\{\begin{array}{l}{\sigma }_{11}={\sigma }_x{\cos }^2\theta \\ {\sigma }_{22}={\sigma }_x{\sin }^2\theta \\ {\tau }_{12}={\sigma }_x\cos \theta \sin \theta \end{array}$ (17)

将上式代入式(7)中，等效应力可得

$\overline{\sigma }=h\left(\theta \right){\sigma }_x$ (18)

其中$h\left(\theta \right)$ 为角度函数，可表示为

$h\left(\theta \right)=\{\begin{array}{l}\sqrt{\frac{3}{2}\left[{\cos }^4\theta +a_{22}{\sin }^4\theta +\left(2a_{66}-1\right){\cos }^2\theta {\sin }^2\theta \right]}+b\left({\sigma }_x\ge 0\right)\\ -\sqrt{\frac{3}{2}\left[{\cos }^4\theta +a_{22}{\sin }^4\theta +\left(2a_{66}-1\right){\cos }^2\theta {\sin }^2\theta \right]}+b\left({\sigma }_x<0\right)\end{array}$ (19)

同理，可以得到x方向塑性应变增量$\text{d}{\epsilon }_x^p$ 与材料主方向塑形应变增量之间的关系为

$\text{d}{\epsilon }_x^p=d_{11}^p{\cos }^2\theta +d_{22}^p{\sin }^2\theta -\text{d}{\gamma }_{12}^p\sin \theta \cos \theta$ (20)

将式(9)、式(16)代入式(19)，可得沿加载方向的塑性应变增量：

$\text{d}{\epsilon }_x^p=h\left(\theta \right)\text{d}{\overline{\epsilon }}^p$(21)

上式积分可得

${\overline{\epsilon }}^p={\epsilon }_x^p/h\left(\theta \right)$ (22)

由上述可见，可利用试验先确定出${\sigma }_x-{\epsilon }_x^p$ 关系，再根据式(18)和式(22)将其转换为等效应力和等效塑性应变($\overline{\sigma }-\overline{\epsilon }$ )关系。应该注意的是，材料等效应力$\overline{\sigma }$ 与等效塑性应${\overline{\epsilon }}^p$ 之间的关系对于不同纤维角$\theta$ 的试件是固定不变的，因此，如果通过对式(7)中的a₂₂、a₆₆和b进行合理的取值能够使得不同偏轴拉伸/压缩试件的$\overline{\sigma }-\overline{\epsilon }$ 曲线分别聚合成一条主曲线，那么该弹塑性本构模型就是有效的。

4. 木材偏轴实验

4.1. 木材基本参数识别

截取由试验得到的偏轴拉伸/试件应力–应变曲线达到峰值荷载之前的部分，汇总如图4所示。

(a) 拉伸 (b) 压缩

Figure 4. Stress-strain curves of tensile/compressive tests for specimens with different off-axis angles

图4. 不同偏轴角度试样的拉伸/压缩应力–应变曲线

不同试件的弹性模量为应力–应变曲线初始段进行线性拟合得到的斜率值，如式(23)所示：

$E_x=\Delta {\sigma }_x/\Delta {\epsilon }_x$ (23)

式中：${\epsilon }_x$ 为试件加载方向的应变，考虑到可能的偏压/偏拉效应，最后的应变数据采用试件前后两面测试结果的均值。

为了得到真实可靠的弹性模量，拟合范围的选取至关重要。拟合范围选取的方法是：从起始的10个应力–应变数据点开始，逐个增加数据点，取应力和应变之间的相关系数最大的一组数据点进行线性拟合。不同偏轴拉伸/压缩试件弹性模量的计算结果列于表1。

Table 1. Test results of elastic modulus for tensile/compressive specimens with different off-axis angles

表1. 不同偏轴角度拉伸/压缩试件的弹性模量测试结果

纤维角θ (˚)	弹性模量MPa
	拉伸	压缩
0	15,418	13,220
15	7690	6570
30	5250	4650
45	3990	3040
60	2700	2580
75	2540	2400
90	2370	2260

4.2. 模型待定参数识别

(a) 拉伸 (b) 压缩

Figure 5. Aggregation and regression results of equivalent stress-equivalent plastic strain curves for off-axis tension/compression

图5. 偏轴拉伸/压缩等效应力–等效塑性应变曲线聚合及回归效果

获得不同偏轴拉伸/压缩试件弹性模量后，将其分别代入式(24)，可得不同偏轴试件沿x方向的塑形应变${\epsilon }_x^p$ 为：

${\epsilon }_x^p={\epsilon }_x-\frac{{\sigma }_x}{E_x}$ (24)

确定出不同偏轴试件$\overline{\sigma }-\overline{\epsilon }$ 关系后，再根据式(18)和式(22)将其转换为等效应力和等效塑性应变($\overline{\sigma }-\overline{\epsilon }$ )关系，并通过选定合适的a₂₂、a₆₆和b采用直观判断的方法使得不同偏轴拉伸/压缩试件的$\overline{\sigma }-\overline{\epsilon }$ 曲线分别聚合成一条主曲线。为了更加科学地评判曲线的聚合效果，对于a₂₂、a₆₆和b每一个(组)取值得到的$\overline{\sigma }-\overline{\epsilon }$ 曲线分别利用式(10)进行拟合，当非线性拟合决定系数最接近于1时，此时a₂₂、a₆₆和b的取值就是本构模型中未知参数的最终取值。应用上述最优参数值表征出的$\overline{\sigma }-\overline{\epsilon }$ 主曲线如图5所示。基于试验所取得的各参数值如表2所示。

Table 2. Model parameters

表2. 模型参数

荷载	a22	a66	b	A	n
拉伸	14	2.4	0.4	3.34 × 10^−^2	2.33
压缩	14	2.4	0.4	1.49 × 10^−^4	3.85

4.3. 模型对比结果及讨论

利用拟合得到的关系式${\overline{\epsilon }}^p=A{\overline{\sigma }}^n$ 返回模型预测不同偏轴拉伸/压缩试件应力–应变关系时，由前述式(2)可知，加载方向上总应变为弹性应变与塑性应变之和，即：

${\epsilon }_x=\frac{\overline{\sigma }}{h\left(\theta \right)E_x}+A{\overline{\sigma }}^{n}h\left(\theta \right)$ (25)

加载方向上应力大小为：

${\sigma }_x=\frac{\overline{\sigma }}{h\left(\theta \right)}$ (26)

本文分别计算了偏轴角度为15˚、30˚、45˚、60˚、75˚及90˚试件在拉伸、压缩下的应力–应变曲线，并与试验结果进行了对比。图6(a)、图6(b)和图7(a)、图7(b)分别表示在偏轴拉伸和压缩下的结果。

(a) 15˚，30˚，45˚试件 (b) 60˚，75˚，90˚试件

Figure 6. Comparison between measured tensile stress-strain relationship and constitutive model predictions

图6. 受拉应力–应变关系的实测结果和本构模型的预测结果对比

由图6可见，在单调拉伸载荷下，模型能够准确描述木材的非线性行为。观察图7可以发现，模型的预测结果与试验数据在单调压缩载荷下整体上吻合良好，同时存在个别试件在应变大于4000的范围内出现模型的预测结果略低于试验结果，其中以$\theta ={75}^{\circ }$ 和$\theta ={90}^{\circ }$ 最为突出，其偏差约为5%。

(a) 15˚，30˚，45˚试件 (b) 60˚，75˚，90˚试件

Figure 7. Comparison between measured compressive stress-strain relationship and constitutive model predictions

图7. 受压应力–应变关系的实测结果和本构模型的预测结果对比

5. 结论

通过引入Drucker-Prager塑性理论的静水压概念，将杨娜等的双参数塑性模型扩展到包括木材塑性流动对载荷方向的依赖性，提出了一种考虑木材同一纹理拉压不等强的本构模型。对白蜡偏轴试样进行了单轴拉伸/压缩试验，基于实验结果，给出了一个含三参数的非线性本构表达式。将模型预测结果与离轴拉伸/压缩试验结果进行了比较，结果表明两者具有良好的相关性，所建立的塑性本构模型能够很好地描述木材在载荷作用下的模型仅适用于非线性力学行为，而且考虑了木材同一纹理方向拉伸和压缩性能的差异。