统 仿 真 学 报 V ol. 22 No. 6
2010年6月 Journal of System Simulation Jun., 2010
第22卷第6期 系
应力波传播模型及其在木材检测中的应用
冯海林, 李光辉, 方益明, 李 剑
（浙江林学院信息工程学院，杭州 311300）
摘 要：木材无损检测技术是充分有效利用木材资源的重要手段之一。依据各向同性和各向异性材
料中的机械波传播理论，研究了应力波在木材中的传播过程。分析了直角坐标系下的应力波传播方程，研究了木材中应力波传播的弹性常数之间的关系。在圆柱极坐标下研究了传播模型，对模型的参数进行了分析，用幂级数的形式给出模型的部分解，分析了木材中横向传播速度。最后对两类木材中的应力波传播进行了计算。
关键词：木材检测；应力波传播；弹性常数；波传播模型
中图分类号：TP391.9 文献标识码：A 文章编号：1004-731X (2010) 06-1490-04
Stress Wave Propagation Modeling and Application in Wood Testing
FENG Hai-lin, LI Guang-hui, FANG Yi-ming, LI Jian
(School of Information Engineering, Zhejiang Forestry University, Hangzhou 311300, China)
Abstract: Wood nondestructive testing technology plays an important role in the utilization of wood resources effectively. According to the theory of mechanical wave transmitting in isotropic and anisotropic material, the stress wave propagation in wood was studied. The differential equations of stress wave propagation were analyzed in rectangular coordinate system. The relationship among the elasticity constants of stress wave propagation in wood was given. The parameters and possible solutions of the model were studied with the power series. The radial transmitting velocity in wood was analyzed, and two woods were used as examples to simulate the stress wave propagation.
Key words: wood testing; stress wave propagation; elastic constant; wave propagation model
引 言木材检测是一个挑战性的问题。传统的人工木材检测难以满足非破坏性快速检测的需要，因此无损检测得到快速发展，并形成一门实用的木材检测技术。应力波无损检测技术是多种木材无损检测技术的一种, 可以通过测量应力波传播时间的变化来确定木材材料的性质。与其他无损检测技术相比, 应力波无损检测技术具有其独特优势。应力波的传播距离远, 抗干扰能力较强, 不受被测木材形状和尺寸的限制，不损坏被测木材，同时应力波无损检测设备小巧, 可方便携带, 因此得到了快速的发展，并在木材检测等行业得到广泛的应用。
基于应力波的木材无损检测有很多的研究，但是使用理论的模型来对应力波在原木材中的传播进行的分析还比较少。杨学春对一维应力波的轴向、径向和弦向给出了一定的分析，在其推导过程中做出了一些假设，如：在原木内部的微小单元运动过程中，截面上应力分布均匀，在截面上的每一个点都向x 方向运动，而且微小单元沿横向的变形是忽略不计的。而原木的尺寸并不能严格满足这些要求。极端情
[1]
况下，一个原木的横向尺寸非常大，由于它与平面纵波传播方向相垂直，横向运动会被阻止。
国外学者利用应力波无损检测技术对木材内部的腐朽和材质也有很多研究，如利用应力波计时器计算木材试件的弹性模量[3]，通过测试应力波传播时间的变化判断木材的内部缺陷[4]。测试应力波传播特征和活立木弹性模量之间的关系[5]。还有一些关于波在各向同性或各向异性材料中的传播模型及应用，如各向同性的材料中的弹性波传播的动力学过程
[6]，应力波的弹性形变对多层材质的影响[7]，非线性lamp 波方法描述固体表面的弹性各向异性特征[8]，圆形管中的各向异性弹性波特性[9,10]，在正交各向异性材料中的三维振动分析[13]，固体中的弹性波传播过程及其反射现象[14]，利用Lord-Shulman 理论在对各项同性含空洞材料中的波传播中的研究[15]等等。
本文依据各向同性和各向异性材料中的机械波传播理论，研究了应力波在木材中传播过程。分析了直角坐标系下的应力波传播方程，给出了木材中应力波传播的弹性常数之间的关系。在圆柱极坐标下研究了传播模型，对模型的参数进行了分析，用幂级数的形式给出模型的部分解，分析了木材中横向传播速度，最后对木材中的应力波传播进行了仿真计算。
收稿日期：2008-12-25 修回日期：2009-03-01
基金项目：国家自然科学基金 (60903144)，浙江省科技厅重点项目(2007C21045)，浙江省自然科学基金项目 (Y1090766,Y107612)，浙江省高校优秀青年教师计划 (2273000011)，校科研基金项目 (2351000848) 作者简介：冯海林(1980-), 男, 安徽怀宁人, 博士，副教授，研究方向为复杂系统建模，智能信号处理，无线传感器网络; 李光辉(1970-), 男, 湖南资兴人, 博士后，教授，研究方向为容错计算，复杂系统建模，无线传感器网络。
1 应力波检测基本原理
基于应力波的木材无损检测的工作原理为用脉冲锤以一定的频率撞击木材，使木材内部产生应力波，当应力波为
2010年6月 冯海林,
等：应力波传播模型及其在木材检测中的应用 Jun., 2010
⎧∂v ∂σ
ρ0=⎪⎪∂t ∂A (3b) ⎨∂∂σv ⎪=ρ0⋅c 2⋅⎪∂A ⎩∂t
接收器收到之后，根据发射和接收的时间间隔，利用计算机和其他工具，分析应力波在木材内部的传播时间和速度的参数的变化情况，从而确定木材的材料性质和内部腐朽情况。
3 木材中的应力波传播与弹性常数关系
3.1 弹性常数关系
引入极坐标(r , θ, z ) ，其中，x =r cos θ，y =r sin θ。如果用C ijkl 来表示和(r , θ, z ) 相关的弹性刚度，那么就有：C ijkl (θ) =Ωip Ωjq Ωkr Ωlt C pqrt (θ) 。
⎛cos θ
⎜
其中：Ωij (θ) =⎜−sin θ
⎜0⎝
sin θcos θ0
0⎞0⎟⎟1⎟⎠
图1 应力波木材无损检测基本原理
由于原木常常被看成一种圆柱形的正交各向异性材料，
2 应力波传播方程
和王礼立[2]一样，在研究的木材中取一个等截面的微小单元体，质点以A 来命名，它在空间中所占的位置称为x ，选取木杆的轴为X 轴，杆的密度，原始截面积等都和坐标的选取无关。
因此，介质的运动可以看作：质点A 在不同的时间t 取不同的空间位置x ，所以有x =x (A,t ) 。现在纵向的应力应变关系设为σx =σx (εx ) ，杆在变形时横截面保持为平面，沿截面只有均匀分布的轴向应力，于是各个运动参量都只是应变ε，质点速度ν，应力σX 和t 的函数。注意：位移u ，等均直接表示X 方向的分量，得到：
∂u ∂u
由：ε=和v =
∂A ∂t ∂v ∂ε
得到： (1) =
∂A ∂t
考察微小单元体A ，由两边总力平衡，得到：
∂v ∂P
ρ0A 0dA =P (A +dA , t ) −P (A , t ) =dA
∂t ∂A 引入：σ=P /A 0
∂v ∂σ
则得到：ρ0 (2) =
∂t ∂A 由应变率无关理论，材料本构关系可以写成：
对于这样的材料，弹性常数是一个四阶对称张量，完善的写法应有四个足标。因此本文只研究在正交各向异性条件下的应力波传播模型，此时C ijkl (θ) 是常量，所以有C ijkl (θ) =C ijkl (0)。木材常为了方便，常用二个缩写下标代替
四个足标[12]，如下表所示：
表1 足标的表示
11
22
33
23 13 12
1 2 3 4 5 6
在圆柱正交各向异性的介质中，用C αβ表示弹性刚度，由于对称性，独立于z 的运动仅仅依赖6个弹性刚度，即：
C 11, C 12, C 22, C 44, C 55, C 66。因此，矩阵简化为：
⎡C 11C 12
⎢C C 22⎢21
⎢C 31C 32⎢
0⎢0
⎢00⎢
0⎣0
C 13
C 23C 33000
000C 4400
0000C 550
0⎤0⎥⎥0⎥ ⎥0⎥0⎥⎥C 66⎦
如果是各向同性的材料，C 12=C 23=C 13=λ，
C 11=C 22=C 33=λ+2μ且C 44=C 55=C 66=μ，其中，λ和μ是拉梅常数。
木材是典型的高度各向异性的，比如对于云杉，给出
σx =σx (εx )
由于σ(ε) 一般是连续可微的
1d σ
引入：C 2=⋅
ρ0d ε
此处我们考虑的应力总是随着应变单调上升函数，所以
d σ2
>0，密度也为正，所以C >0 d ε
C 11 2C 22，C 55 2C 66
在计算弹性刚度之前，先利用木材的正交各向异性给出弹性顺度。
⎡S 11⎢S ⎢21⎢S 31⎢⎢0⎢0⎢⎣0
S 12S 22S 32000
S 13S 23S 33000
000S 4400
0000S 550
⎤⎥⎥⎥ ⎥⎥⎥⎥S 66⎦00000
代入得到：
∂v ∂ε∂σ∂v
，以及 =ρ0=ρ0⋅c 2⋅
∂t ∂A ∂t ∂A
因此问题就转化为：
⎧∂v ∂ε
=⎪⎪∂A ∂t (3a) ⎨
ε∂∂v ⎪=ρ0
⎪∂t ∂A ⎩
弹性刚度可以利用弹性顺度表示，如： C 11=(S 22S 33−(S 23) 2) /S
C 22=(S 11S 33−(S 13) 2) /S (4)
C 33=(S 22S 11−(S 13) 2) /S
2010年6月 系 统 仿 真 学 报 Jun., 2010
其他的参数也可以类似得出，其中：
S =S 11S 22S 33+2S 12S 23S 31−S 11S 232−S 22S 132−S 33S 122
4 应力波传播计算结果
下面以云杉和松木为对象对应力波在木材中的传播速度进行计算。选择木材中的微小单元坐标如图2所示：
L
R
L T
LT
TR RL
R
3.2 木材中的应力波模型
用u 表示位移向量，ρ表示微小单元密度，t 表示时间，C ijkl 表示弹性刚度常数，木材在压力下微小形变的应变张量
[εij ]和位移的关系为：
εij =(
1∂u i ∂u j
+) (5) 2∂x j ∂x i
结合前面的等式，应力波传播方程变为： ρ
∂2u i ∂2u k
−C =0 (6) ijkl ∂t 2∂x l ∂x j
图2 木材中的微小单元坐标
由于检测木材的应力波横向传播是检测腐朽的较好途径。因此我们仅仅考虑垂直于横截面的运动。类似于
云杉的弹性常数，杨氏模量E i ，泊松比μij ，和剪切弹性模量G ij 如表2所示:
表2 木材的弹性常数
木材密度
Bucur [12]，给出问题的一种解形式：
u (r , θ, t ) =Re i {u m (r ) e jm θe −i ωt } (7)
其中，i 和j 是两个非交互的复杂单元，m 是一个整数，ω是振动频率，Re i 是i 的实部，令u m =(u m , v m , w m ) 。首先对弹性刚度的表示进行简化[11]：
c 1=C 11/C 66，c 2=C 22/C 66，c 12=C 12/C 66
E L E R E T
496573
G L T G L R G TR μRT
690 758 39 0.43676 1172 66 0.68
μLR μLT
0.370.42
0.470.51
云杉0.[1**********]松木0.[1**********]3
通过计算，得到云杉的弹性刚度和弹性顺度矩阵：
⎛ 1.00 0.23 0.48 0 0 0⎞
⎜⎟⎜ 0.23 0.55 0.35 0 0 0⎟ ⎜ 0.48 0.35 11.97 0 0 0⎟⎜⎟⎜ 0 0 0 0.69 0 0⎟⎜ 0 0 0 0 0.76 0⎟⎜⎜ 0 0 0 0 0 0.039⎟⎟⎝⎠
由于我们考虑的是横向传播，所以根据前面的等式有：
'' ' '
c 1(r 2u m +r m ) +jm (1+c 12) rv m +(k 2r 2−m 2−c 2) u m −
jm (1+c 2) v m =0
'' ' ' r 2v m +rv m +jm (1+c 12) ru m +(k 2r 2−1−m 2c 2) v m +
(8)
jm (1+c 2) u m =0
(9)
其中，很多研究都是使用Frobenius 级数[13]求解这类问题。如果假设m >0，同样可以用以下幂级数形式来找出解。
∞
∞
u m (r ) =∑a n (kr ) 2n +α, v m (r ) =j ∑b n (kr ) 2n +α
n =0
n =0
其中，系数a n , b n , α待定。下面分析几个待定系数[11,13]。通再令(kr ) 2n +α的系数为0，不妨假设其过把u 和v 代入原式，
因此可以得到A 的两个根。令其为α1和判别式是大于0的，
⎛ 1.120 -0.450 -0.032 0 0 0⎞
⎜⎟⎜ -0.450 2.020 -0.041 0 0 0⎟ ⎜ -0.032 -0.041 0.086 0 0 0⎟⎜⎟⎜ 0 0 0 1.45 0 0⎟⎜ 0 0 0 0 1.32 0⎟⎜⎜ 0 0 0 0 0 25.64⎟⎟⎝⎠
同样得到松木的弹性刚度和弹性顺度矩阵：
⎛ 1.45 0.51 0.34 0 0 0⎞⎜⎟⎜ 0.51 0.75 0.23 0 0 0⎟ ⎜ 0.34 0.23 6.35 0 0 0⎟⎜⎟⎜ 0 0 0 0.68 0 0⎟⎜ 0 0 0 0 1.18 0⎟⎜⎜ 0 0 0 0 0 0.066⎟⎟⎝⎠⎛ 0.91 -0.61 -0.03 0 0 0⎞⎜⎟⎜ -0.61 1.76 -0.03 0 0 0⎟ ⎜ -0.03 -0.03 0.16 0 0 0⎟⎜⎟⎜ 0 0 0 1.47 0 0⎟⎜ 0 0 0 0 0.85 0⎟⎜⎜ 0 0 0 0 0 15.15⎟⎟⎝⎠
α2，选择α1和α2的正值作为我们的解。注意：这里得到的α2的根是依赖于系数m 的。当确定一个a 的值之后，代入
可以得出a 0和b 0，通过a n 和b n 的递推关系可以找到系数。
3.3 木材中的应力波传播速度
假设一个平面谐波在单位矢量方向传播的位移为u ，和波阵面垂直
[12]
，于是方程可以表示为：
(C ijkl n j n k −δik ρV 2) P m =0 (10) 其中，δik 为 Kronecker张量，P m 是在位移方向的单位向量的分量，n k 为所考虑平面的入射方向同x 、y 、z 轴夹角的方向余弦。ρ为密度。v 表示速度。
考虑木材的径向或弦向传播速度，在径向和弦向面上。可以得到如下等式：
22
(C 13+C 55) 2n 12n 3=(C 11n 12+C 55n 3−ρ⋅V 2) 2(C 55n 12+C 33n 3−ρ⋅V 2)
从云杉的数据中可以得到：c 1=25.64，c 2=5.96，
c 12=14.2，γ=0.744。从而可以根据前面的步骤求得α1(0)=0.74，α2(0)=1.00，继而递推得到幂级数的系数。
同时，利用上述的模型和数据，对云杉和松木中的应力波传播速度进行了计算。
(11)
2010年6月 冯海林,
等：应力波传播模型及其在木材检测中的应用 Jun., 2010
Oriented Strandboard [J]. Forest Products Journal (S0015-7473), 2003, 53(3): 33-35. [4]
Lee J, Kim K, Bae M. Patterns of Resistographs for Evaluating Deteriorated Structural Wood Members [J]. Journal of the Korean Wood Science and Technology (S1017-0715), 2003, 31(6): 45-54. [5]
Udaya B Halabe, Gangadhar M Bidigalu, Hota V S GangaRao, Robert J Ross. Nondestructive Evaluation of Green Wood Using Stress Wave and Transverse Vibration Techniques [J]. Materials Evaluation (S0025-5327), 1997, 55(9): 1013-1018. [6]
Bostro¨m A, Johansson G, Olsson P. On the Rational Derivation of a Hierarchy of Dynamic Equations for a Homogeneous, Isotropic, Elastic Plate [J]. Int. J. Solids Struct (S0020-7683), 2001, 38(12): 2487-2501. [7]
Tasdemirci A, Hall I W. The effects of plastic deformation on stress wave propagation in multi-layer materials [J]. International Journal of Impact Engineering (S0734-743X), 2007, 34(11): 1797-1813. [8]
Deng M X, Yang, J. Characterization of elastic anisotropy of a solid plate using nonlinear Lamb wave approach [J]. Journal of Sound and Vibration (S0022-460X), 2007, 308(1-2): 201-211. [9]
Ting T C T. Pressuring, Shearing,Torsion and Extension of a Circular Tube or Bar of Cylindrically Anisotropic Material [J]. Proc.Roy. Soc. A (S0962-8452), 1996, 452: 2397-2421.
[10] Ting T C T. Anisotropic Elasticity [M]. UK: Oxford University Press,
1996.
[11] Martin P A, Berger J R. Waves in Wood: Free Vibrations of a Wooden
Pole [J]. J. Mech. Phys. Solids (S0022-5096), 2001, 49: 1155-1178. [12] V Bucur. Acoustics of Wood [M]. Boca Raton, FL, USA: CRC Press,
1995.
[13] Chou F H, Achenbach J D. Three-dimensional Vibrations of
Orthotropic Cylinders [J]. Proc. ASCE J. Eng. Mech. Div (S0044-7951), 1972, 98(4): 813-822.
[14] Tomar S K, Khurana A. Elastic Waves in an electro-microelastic solid
[J]. International Journal of Solids and Structures (S0020-7683), 2008, 45(1): 276-302.
[15] Singh B. Wave propagation in a generalized thermoelastic material
with voids [J]. Applied Mathematics and Computation (S0096-3003), 2007, 189(1): 698-709.
由于木材发生腐朽或虫蛀等缺陷时，垂直于木材纹理方向的传播速度发生改变，因此，横向（径向或弦向）测量是检测腐朽的较好途径，所以我们计算的是木材的横向传播速度。同时，我们选用德国生产的应力波检测仪Electronic
Hammer 对两种木材的横向传播速度进行了50次单独的检测，其中每次检测结果均由3次检测取平均值得到。两种木材的测量结果和计算出的结果分别有2.05% 和2.94%的平均误差。云杉中的应力波传播速度要比松木小大约140m/s，它的速度随着传播角度的改变而改变。松木的横向传播速度更易受偏振角的影响，速度也更快。
5 结论
木材无损检测技术可以对木材内部的腐朽和材质进行检测，对于有效和充分的利用木材资源有着重要的意义。应力波无损检测技术是多种木材无损检测技术的一种, 可以通过测量应力波传播时间的变化来确定木材材料的性质。与其他无损检测技术相比, 应力波无损检测技术具有其独特优势。本文依据各向同性和各向异性材料中的机械波传播理论，研究了应力波在木材中传播过程。分析了直角坐标系下的应力波传播方程，给出了木材中应力波传播的弹性常数之间的关系。在圆柱极坐标下研究了传播模型，对模型的参数进行了分析，用幂级数的形式给出模型的部分解，分析了木材中横向传播速度，最后对木材中的应力波传播进行了计算。
参考文献:
[1] [2] [3]
杨学春, 王立海. 应力波在原木中传播理论的研究[J]. 林业科学. 2005, 41(5): 132-138.
王礼立. 应力波基础 [M]. 北京: 国防工业出版社, 1985.
Ross R J, Yang V, Illman B L, et al. Relationship between Stress Wave Transmission Time and Bending Strength of Deteriorated
统 仿 真 学 报 V ol. 22 No. 6
2010年6月 Journal of System Simulation Jun., 2010
第22卷第6期 系
应力波传播模型及其在木材检测中的应用
冯海林, 李光辉, 方益明, 李 剑
（浙江林学院信息工程学院，杭州 311300）
摘 要：木材无损检测技术是充分有效利用木材资源的重要手段之一。依据各向同性和各向异性材
料中的机械波传播理论，研究了应力波在木材中的传播过程。分析了直角坐标系下的应力波传播方程，研究了木材中应力波传播的弹性常数之间的关系。在圆柱极坐标下研究了传播模型，对模型的参数进行了分析，用幂级数的形式给出模型的部分解，分析了木材中横向传播速度。最后对两类木材中的应力波传播进行了计算。
关键词：木材检测；应力波传播；弹性常数；波传播模型
中图分类号：TP391.9 文献标识码：A 文章编号：1004-731X (2010) 06-1490-04
Stress Wave Propagation Modeling and Application in Wood Testing
FENG Hai-lin, LI Guang-hui, FANG Yi-ming, LI Jian
(School of Information Engineering, Zhejiang Forestry University, Hangzhou 311300, China)
Abstract: Wood nondestructive testing technology plays an important role in the utilization of wood resources effectively. According to the theory of mechanical wave transmitting in isotropic and anisotropic material, the stress wave propagation in wood was studied. The differential equations of stress wave propagation were analyzed in rectangular coordinate system. The relationship among the elasticity constants of stress wave propagation in wood was given. The parameters and possible solutions of the model were studied with the power series. The radial transmitting velocity in wood was analyzed, and two woods were used as examples to simulate the stress wave propagation.
Key words: wood testing; stress wave propagation; elastic constant; wave propagation model
引 言木材检测是一个挑战性的问题。传统的人工木材检测难以满足非破坏性快速检测的需要，因此无损检测得到快速发展，并形成一门实用的木材检测技术。应力波无损检测技术是多种木材无损检测技术的一种, 可以通过测量应力波传播时间的变化来确定木材材料的性质。与其他无损检测技术相比, 应力波无损检测技术具有其独特优势。应力波的传播距离远, 抗干扰能力较强, 不受被测木材形状和尺寸的限制，不损坏被测木材，同时应力波无损检测设备小巧, 可方便携带, 因此得到了快速的发展，并在木材检测等行业得到广泛的应用。
基于应力波的木材无损检测有很多的研究，但是使用理论的模型来对应力波在原木材中的传播进行的分析还比较少。杨学春对一维应力波的轴向、径向和弦向给出了一定的分析，在其推导过程中做出了一些假设，如：在原木内部的微小单元运动过程中，截面上应力分布均匀，在截面上的每一个点都向x 方向运动，而且微小单元沿横向的变形是忽略不计的。而原木的尺寸并不能严格满足这些要求。极端情
[1]
况下，一个原木的横向尺寸非常大，由于它与平面纵波传播方向相垂直，横向运动会被阻止。
国外学者利用应力波无损检测技术对木材内部的腐朽和材质也有很多研究，如利用应力波计时器计算木材试件的弹性模量[3]，通过测试应力波传播时间的变化判断木材的内部缺陷[4]。测试应力波传播特征和活立木弹性模量之间的关系[5]。还有一些关于波在各向同性或各向异性材料中的传播模型及应用，如各向同性的材料中的弹性波传播的动力学过程
[6]，应力波的弹性形变对多层材质的影响[7]，非线性lamp 波方法描述固体表面的弹性各向异性特征[8]，圆形管中的各向异性弹性波特性[9,10]，在正交各向异性材料中的三维振动分析[13]，固体中的弹性波传播过程及其反射现象[14]，利用Lord-Shulman 理论在对各项同性含空洞材料中的波传播中的研究[15]等等。
本文依据各向同性和各向异性材料中的机械波传播理论，研究了应力波在木材中传播过程。分析了直角坐标系下的应力波传播方程，给出了木材中应力波传播的弹性常数之间的关系。在圆柱极坐标下研究了传播模型，对模型的参数进行了分析，用幂级数的形式给出模型的部分解，分析了木材中横向传播速度，最后对木材中的应力波传播进行了仿真计算。
收稿日期：2008-12-25 修回日期：2009-03-01
基金项目：国家自然科学基金 (60903144)，浙江省科技厅重点项目(2007C21045)，浙江省自然科学基金项目 (Y1090766,Y107612)，浙江省高校优秀青年教师计划 (2273000011)，校科研基金项目 (2351000848) 作者简介：冯海林(1980-), 男, 安徽怀宁人, 博士，副教授，研究方向为复杂系统建模，智能信号处理，无线传感器网络; 李光辉(1970-), 男, 湖南资兴人, 博士后，教授，研究方向为容错计算，复杂系统建模，无线传感器网络。
1 应力波检测基本原理
基于应力波的木材无损检测的工作原理为用脉冲锤以一定的频率撞击木材，使木材内部产生应力波，当应力波为
2010年6月 冯海林,
等：应力波传播模型及其在木材检测中的应用 Jun., 2010
⎧∂v ∂σ
ρ0=⎪⎪∂t ∂A (3b) ⎨∂∂σv ⎪=ρ0⋅c 2⋅⎪∂A ⎩∂t
接收器收到之后，根据发射和接收的时间间隔，利用计算机和其他工具，分析应力波在木材内部的传播时间和速度的参数的变化情况，从而确定木材的材料性质和内部腐朽情况。
3 木材中的应力波传播与弹性常数关系
3.1 弹性常数关系
引入极坐标(r , θ, z ) ，其中，x =r cos θ，y =r sin θ。如果用C ijkl 来表示和(r , θ, z ) 相关的弹性刚度，那么就有：C ijkl (θ) =Ωip Ωjq Ωkr Ωlt C pqrt (θ) 。
⎛cos θ
⎜
其中：Ωij (θ) =⎜−sin θ
⎜0⎝
sin θcos θ0
0⎞0⎟⎟1⎟⎠
图1 应力波木材无损检测基本原理
由于原木常常被看成一种圆柱形的正交各向异性材料，
2 应力波传播方程
和王礼立[2]一样，在研究的木材中取一个等截面的微小单元体，质点以A 来命名，它在空间中所占的位置称为x ，选取木杆的轴为X 轴，杆的密度，原始截面积等都和坐标的选取无关。
因此，介质的运动可以看作：质点A 在不同的时间t 取不同的空间位置x ，所以有x =x (A,t ) 。现在纵向的应力应变关系设为σx =σx (εx ) ，杆在变形时横截面保持为平面，沿截面只有均匀分布的轴向应力，于是各个运动参量都只是应变ε，质点速度ν，应力σX 和t 的函数。注意：位移u ，等均直接表示X 方向的分量，得到：
∂u ∂u
由：ε=和v =
∂A ∂t ∂v ∂ε
得到： (1) =
∂A ∂t
考察微小单元体A ，由两边总力平衡，得到：
∂v ∂P
ρ0A 0dA =P (A +dA , t ) −P (A , t ) =dA
∂t ∂A 引入：σ=P /A 0
∂v ∂σ
则得到：ρ0 (2) =
∂t ∂A 由应变率无关理论，材料本构关系可以写成：
对于这样的材料，弹性常数是一个四阶对称张量，完善的写法应有四个足标。因此本文只研究在正交各向异性条件下的应力波传播模型，此时C ijkl (θ) 是常量，所以有C ijkl (θ) =C ijkl (0)。木材常为了方便，常用二个缩写下标代替
四个足标[12]，如下表所示：
表1 足标的表示
11
22
33
23 13 12
1 2 3 4 5 6
在圆柱正交各向异性的介质中，用C αβ表示弹性刚度，由于对称性，独立于z 的运动仅仅依赖6个弹性刚度，即：
C 11, C 12, C 22, C 44, C 55, C 66。因此，矩阵简化为：
⎡C 11C 12
⎢C C 22⎢21
⎢C 31C 32⎢
0⎢0
⎢00⎢
0⎣0
C 13
C 23C 33000
000C 4400
0000C 550
0⎤0⎥⎥0⎥ ⎥0⎥0⎥⎥C 66⎦
如果是各向同性的材料，C 12=C 23=C 13=λ，
C 11=C 22=C 33=λ+2μ且C 44=C 55=C 66=μ，其中，λ和μ是拉梅常数。
木材是典型的高度各向异性的，比如对于云杉，给出
σx =σx (εx )
由于σ(ε) 一般是连续可微的
1d σ
引入：C 2=⋅
ρ0d ε
此处我们考虑的应力总是随着应变单调上升函数，所以
d σ2
>0，密度也为正，所以C >0 d ε
C 11 2C 22，C 55 2C 66
在计算弹性刚度之前，先利用木材的正交各向异性给出弹性顺度。
⎡S 11⎢S ⎢21⎢S 31⎢⎢0⎢0⎢⎣0
S 12S 22S 32000
S 13S 23S 33000
000S 4400
0000S 550
⎤⎥⎥⎥ ⎥⎥⎥⎥S 66⎦00000
代入得到：
∂v ∂ε∂σ∂v
，以及 =ρ0=ρ0⋅c 2⋅
∂t ∂A ∂t ∂A
因此问题就转化为：
⎧∂v ∂ε
=⎪⎪∂A ∂t (3a) ⎨
ε∂∂v ⎪=ρ0
⎪∂t ∂A ⎩
弹性刚度可以利用弹性顺度表示，如： C 11=(S 22S 33−(S 23) 2) /S
C 22=(S 11S 33−(S 13) 2) /S (4)
C 33=(S 22S 11−(S 13) 2) /S
2010年6月 系 统 仿 真 学 报 Jun., 2010
其他的参数也可以类似得出，其中：
S =S 11S 22S 33+2S 12S 23S 31−S 11S 232−S 22S 132−S 33S 122
4 应力波传播计算结果
下面以云杉和松木为对象对应力波在木材中的传播速度进行计算。选择木材中的微小单元坐标如图2所示：
L
R
L T
LT
TR RL
R
3.2 木材中的应力波模型
用u 表示位移向量，ρ表示微小单元密度，t 表示时间，C ijkl 表示弹性刚度常数，木材在压力下微小形变的应变张量
[εij ]和位移的关系为：
εij =(
1∂u i ∂u j
+) (5) 2∂x j ∂x i
结合前面的等式，应力波传播方程变为： ρ
∂2u i ∂2u k
−C =0 (6) ijkl ∂t 2∂x l ∂x j
图2 木材中的微小单元坐标
由于检测木材的应力波横向传播是检测腐朽的较好途径。因此我们仅仅考虑垂直于横截面的运动。类似于
云杉的弹性常数，杨氏模量E i ，泊松比μij ，和剪切弹性模量G ij 如表2所示:
表2 木材的弹性常数
木材密度
Bucur [12]，给出问题的一种解形式：
u (r , θ, t ) =Re i {u m (r ) e jm θe −i ωt } (7)
其中，i 和j 是两个非交互的复杂单元，m 是一个整数，ω是振动频率，Re i 是i 的实部，令u m =(u m , v m , w m ) 。首先对弹性刚度的表示进行简化[11]：
c 1=C 11/C 66，c 2=C 22/C 66，c 12=C 12/C 66
E L E R E T
496573
G L T G L R G TR μRT
690 758 39 0.43676 1172 66 0.68
μLR μLT
0.370.42
0.470.51
云杉0.[1**********]松木0.[1**********]3
通过计算，得到云杉的弹性刚度和弹性顺度矩阵：
⎛ 1.00 0.23 0.48 0 0 0⎞
⎜⎟⎜ 0.23 0.55 0.35 0 0 0⎟ ⎜ 0.48 0.35 11.97 0 0 0⎟⎜⎟⎜ 0 0 0 0.69 0 0⎟⎜ 0 0 0 0 0.76 0⎟⎜⎜ 0 0 0 0 0 0.039⎟⎟⎝⎠
由于我们考虑的是横向传播，所以根据前面的等式有：
'' ' '
c 1(r 2u m +r m ) +jm (1+c 12) rv m +(k 2r 2−m 2−c 2) u m −
jm (1+c 2) v m =0
'' ' ' r 2v m +rv m +jm (1+c 12) ru m +(k 2r 2−1−m 2c 2) v m +
(8)
jm (1+c 2) u m =0
(9)
其中，很多研究都是使用Frobenius 级数[13]求解这类问题。如果假设m >0，同样可以用以下幂级数形式来找出解。
∞
∞
u m (r ) =∑a n (kr ) 2n +α, v m (r ) =j ∑b n (kr ) 2n +α
n =0
n =0
其中，系数a n , b n , α待定。下面分析几个待定系数[11,13]。通再令(kr ) 2n +α的系数为0，不妨假设其过把u 和v 代入原式，
因此可以得到A 的两个根。令其为α1和判别式是大于0的，
⎛ 1.120 -0.450 -0.032 0 0 0⎞
⎜⎟⎜ -0.450 2.020 -0.041 0 0 0⎟ ⎜ -0.032 -0.041 0.086 0 0 0⎟⎜⎟⎜ 0 0 0 1.45 0 0⎟⎜ 0 0 0 0 1.32 0⎟⎜⎜ 0 0 0 0 0 25.64⎟⎟⎝⎠
同样得到松木的弹性刚度和弹性顺度矩阵：
⎛ 1.45 0.51 0.34 0 0 0⎞⎜⎟⎜ 0.51 0.75 0.23 0 0 0⎟ ⎜ 0.34 0.23 6.35 0 0 0⎟⎜⎟⎜ 0 0 0 0.68 0 0⎟⎜ 0 0 0 0 1.18 0⎟⎜⎜ 0 0 0 0 0 0.066⎟⎟⎝⎠⎛ 0.91 -0.61 -0.03 0 0 0⎞⎜⎟⎜ -0.61 1.76 -0.03 0 0 0⎟ ⎜ -0.03 -0.03 0.16 0 0 0⎟⎜⎟⎜ 0 0 0 1.47 0 0⎟⎜ 0 0 0 0 0.85 0⎟⎜⎜ 0 0 0 0 0 15.15⎟⎟⎝⎠
α2，选择α1和α2的正值作为我们的解。注意：这里得到的α2的根是依赖于系数m 的。当确定一个a 的值之后，代入
可以得出a 0和b 0，通过a n 和b n 的递推关系可以找到系数。
3.3 木材中的应力波传播速度
假设一个平面谐波在单位矢量方向传播的位移为u ，和波阵面垂直
[12]
，于是方程可以表示为：
(C ijkl n j n k −δik ρV 2) P m =0 (10) 其中，δik 为 Kronecker张量，P m 是在位移方向的单位向量的分量，n k 为所考虑平面的入射方向同x 、y 、z 轴夹角的方向余弦。ρ为密度。v 表示速度。
考虑木材的径向或弦向传播速度，在径向和弦向面上。可以得到如下等式：
22
(C 13+C 55) 2n 12n 3=(C 11n 12+C 55n 3−ρ⋅V 2) 2(C 55n 12+C 33n 3−ρ⋅V 2)
从云杉的数据中可以得到：c 1=25.64，c 2=5.96，
c 12=14.2，γ=0.744。从而可以根据前面的步骤求得α1(0)=0.74，α2(0)=1.00，继而递推得到幂级数的系数。
同时，利用上述的模型和数据，对云杉和松木中的应力波传播速度进行了计算。
(11)
2010年6月 冯海林,
等：应力波传播模型及其在木材检测中的应用 Jun., 2010
Oriented Strandboard [J]. Forest Products Journal (S0015-7473), 2003, 53(3): 33-35. [4]
Lee J, Kim K, Bae M. Patterns of Resistographs for Evaluating Deteriorated Structural Wood Members [J]. Journal of the Korean Wood Science and Technology (S1017-0715), 2003, 31(6): 45-54. [5]
Udaya B Halabe, Gangadhar M Bidigalu, Hota V S GangaRao, Robert J Ross. Nondestructive Evaluation of Green Wood Using Stress Wave and Transverse Vibration Techniques [J]. Materials Evaluation (S0025-5327), 1997, 55(9): 1013-1018. [6]
Bostro¨m A, Johansson G, Olsson P. On the Rational Derivation of a Hierarchy of Dynamic Equations for a Homogeneous, Isotropic, Elastic Plate [J]. Int. J. Solids Struct (S0020-7683), 2001, 38(12): 2487-2501. [7]
Tasdemirci A, Hall I W. The effects of plastic deformation on stress wave propagation in multi-layer materials [J]. International Journal of Impact Engineering (S0734-743X), 2007, 34(11): 1797-1813. [8]
Deng M X, Yang, J. Characterization of elastic anisotropy of a solid plate using nonlinear Lamb wave approach [J]. Journal of Sound and Vibration (S0022-460X), 2007, 308(1-2): 201-211. [9]
Ting T C T. Pressuring, Shearing,Torsion and Extension of a Circular Tube or Bar of Cylindrically Anisotropic Material [J]. Proc.Roy. Soc. A (S0962-8452), 1996, 452: 2397-2421.
[10] Ting T C T. Anisotropic Elasticity [M]. UK: Oxford University Press,
1996.
[11] Martin P A, Berger J R. Waves in Wood: Free Vibrations of a Wooden
Pole [J]. J. Mech. Phys. Solids (S0022-5096), 2001, 49: 1155-1178. [12] V Bucur. Acoustics of Wood [M]. Boca Raton, FL, USA: CRC Press,
1995.
[13] Chou F H, Achenbach J D. Three-dimensional Vibrations of
Orthotropic Cylinders [J]. Proc. ASCE J. Eng. Mech. Div (S0044-7951), 1972, 98(4): 813-822.
[14] Tomar S K, Khurana A. Elastic Waves in an electro-microelastic solid
[J]. International Journal of Solids and Structures (S0020-7683), 2008, 45(1): 276-302.
[15] Singh B. Wave propagation in a generalized thermoelastic material
with voids [J]. Applied Mathematics and Computation (S0096-3003), 2007, 189(1): 698-709.
由于木材发生腐朽或虫蛀等缺陷时，垂直于木材纹理方向的传播速度发生改变，因此，横向（径向或弦向）测量是检测腐朽的较好途径，所以我们计算的是木材的横向传播速度。同时，我们选用德国生产的应力波检测仪Electronic
Hammer 对两种木材的横向传播速度进行了50次单独的检测，其中每次检测结果均由3次检测取平均值得到。两种木材的测量结果和计算出的结果分别有2.05% 和2.94%的平均误差。云杉中的应力波传播速度要比松木小大约140m/s，它的速度随着传播角度的改变而改变。松木的横向传播速度更易受偏振角的影响，速度也更快。
5 结论
木材无损检测技术可以对木材内部的腐朽和材质进行检测，对于有效和充分的利用木材资源有着重要的意义。应力波无损检测技术是多种木材无损检测技术的一种, 可以通过测量应力波传播时间的变化来确定木材材料的性质。与其他无损检测技术相比, 应力波无损检测技术具有其独特优势。本文依据各向同性和各向异性材料中的机械波传播理论，研究了应力波在木材中传播过程。分析了直角坐标系下的应力波传播方程，给出了木材中应力波传播的弹性常数之间的关系。在圆柱极坐标下研究了传播模型，对模型的参数进行了分析，用幂级数的形式给出模型的部分解，分析了木材中横向传播速度，最后对木材中的应力波传播进行了计算。
参考文献:
[1] [2] [3]
杨学春, 王立海. 应力波在原木中传播理论的研究[J]. 林业科学. 2005, 41(5): 132-138.
王礼立. 应力波基础 [M]. 北京: 国防工业出版社, 1985.
Ross R J, Yang V, Illman B L, et al. Relationship between Stress Wave Transmission Time and Bending Strength of Deteriorated