文章信息
- 郝晓峰, 俞昌铭, 蒋佳荔, 吕建雄, 徐康
- Hao Xiaofeng, Yu Changming, Jiang Jiali, Xu Kang
- 木材快速生长期内早材与晚材密度分布的数学描述初探
- A Preliminary Study on Modeling of Earlywood and Latewood Density Distribution during the Fast Growth Period
- 林业科学, 2013, 49(10): 118-126
- Scientia Silvae Sinicae, 2013, 49(10): 118-126.
- DOI: 10.11707/j.1001-7488.20131019
-
文章历史
- 收稿日期:2012-11-22
- 修回日期:2013-03-01
-
作者相关文章
2. 北京科技大学机械工程学院 北京 100083;
3. 中国林业科学研究院林业新技术研究所 北京 100091;
4. 中南林业科技大学材料科学与工程学院 长沙 410004
2. School of Mechanical Engineering,University of Science and Technology Beijing Beijing 100083;
3. Research Institute of New Technology,CAF Beijing 100091;
4. College of Material Science and Engineering,Central South University of Forestry and Technology Changsha 410004
木材密度是木材重要的物理性质之一。一般而言,木材密度是指木材的宏观密度,表示为$\bar \rho = m/V = \int_V {\rho {\rm{d}}V/V} $,是木材质量 m 与其体积V的比值;ρ为木材的微观密度,是指木材内部在某一狭小范围内的密度。研究木材生长轮内早、晚材密度的变化规律,可以揭示不同生长条件下木材材质变异规律,为生态木材学、树木年代学等学科提供研究方法,而且也为速生林的人工抚育措施以及最佳轮伐周期的确定提供理论依据,为木材加工提供重要的质量参数(李坚,2002)。
木材宏观密度的研究主要体现在以下几个领域: 在林产品加工利用方面,木材密度影响着纸浆出浆率(Blair et al.,1976)、结构用材强度(Zobel et al.,1989 ; 吕建雄等,2009); 在树木育种改良方面,木材密度具有可遗传的特性(Sprague et al.,1983); 在人工林木材性质方面,针对幼龄材密度的研究也是近些年研究的重点(Zobel et al.,1989 ; Lü et al.,2003 ; 吕建雄等,2005); 在树木年代学方面,木材密度对研究气候变化有辅助作用(Kim,1995 ;Yasue et al.,1995 ; Anta et al.,2006 ; 徐金梅等,2011)。
关于木材微观密度的研究,Ferrand(1982)研究了生长轮的形状与密度分布,建立了回归分析模型,但该模型中没有划分早晚材。Schweingruber(1988)提出了2种确定过渡密度值的方法来定量描述生长轮内早材和晚材的体积比率: 第1种方法是设定一个固定的密度值,在此密度值之上的为晚材,之下的为早材; 第2种方法是根据生长轮内最大密度值与最小密度值的均值为界线来划分早材和晚材。Pernestål等(1995)认为低密度早材与高密度晚材之间的过渡区域可用泊松分布描述,建立了木材生长轮内的密度分布模型,并通过试验确定了模型中的6个参数。
本文基于生长环境影响树木生长速率、生长速率进而影响木材的早晚材密度这一假定,将树木的生长环境、生长速率与木材微观密度相结合,建立生长轮内早晚材密度沿径向分布的数学模型。早晚材的密度变化影响着木材干燥过程的传热传质行为。在以往的研究中,研究者们在建立木材干燥传热传质模型时将木材假定均质材料,即宏观密度等于微观密度(Perre et al.,1993 ; Pang,1998 ; Yu et al.,2007)。事实上,木材微观密度的不均匀性直接影响着木材导热系数、水分扩散率等热物性参数。因此,在构建木材干燥传热传质物理模型及相应数学模型时,为提高模型的精确度,应考虑木材微观密度的不均匀性这一特征(俞昌铭,2011)。本文旨在木材宏观密度与微观密度之间建立一种内在的数学联系,为优化木材干燥过程中传热传质模型提供理论依据。
1 树木生长环境的物理与数学模型树木生长是一种生理现象,由内因与外因共同作用引起,内因是指树木自身的条件,包括树种、种源、家系、枝下高、冠幅等遗传因素; 外因是指林分密度、降雨量、温度、光照、土壤养分等外在因素(Rennols et al.,1985 ; Kim,1995 ; Yasue et al.,1995 ; 胡晓龙,2003 ; 张建国等,2004)。对于同一片人工林,认为影响树木生长的内因相同,则影响树木生长的因素为温度、光照、降雨量、土壤养分等。
图 1是湖南省常德市5年间气温与光照时间(此数据来源于中国气象科学数据共享服务网)的变化曲线。由图可见,温度与光照二者的变化规律基本相同; 温度与光照相比,其分布更规律,易于数学描述。从树木生长的生理角度看,一般情况下,当环境气温低于临界温度Tc(如4.4 ℃)且持续1周以上时,树木停止生长,进入休眠状态(尹思慈,1996)。
水分的供给也是树木生长的基本要素之一。 Fritts(1976)和Kim(1995)认为土壤中的水分是影响树木生长的重要因素之一。土壤中的水分主要来源于大气降水,它直接影响土壤的含水量及肥沃程度,因此降雨量对树木的生长起重要作用。与温度相比,降雨量在1年内的分布具有较大随机性,不能用连续函数进行描述,因此,采用阶梯函数表示。图 2是湖南省常德市5年间的降雨量分布情况(此数据来源于中国气象科学数据共享服务网),从图 2可看出,降雨量主要集中在春、夏两季。
本研究中,为了便于数学描述,对树木生长的环境进行简化,引入以下假定进行物理模型的构建。
1.1 物理模型假定1,将环境温度、光照、降雨量、土壤养分等外部因素对树木生长的影响假定为只有温度与降雨量2个因素的影响。由图 1可看出,温度与光照随时间的变化趋势相似,假定温度的因素涵盖了光照的因素; 同样,假定降雨量的因素涵盖了土壤中水分的因素; 忽略了土壤养分等诸多环境因素对树木生长的影响,完全忽略各种不测气候因素(如暴雨、风灾等)的影响。
假定2,寒温带1年内环境温度的变化按正弦波起伏波动,波动的周期为1年。全年最高温度为正弦波波峰,最低温度为正弦波波谷,振幅为最高(最低)温度与年平均温度的差值。严格地说,环境温度除了按年周期波动外,按天也是周期波动的。这里,忽略了一天之内的温度变化。
假定3,1年内降雨量的分布用分段线性函数表示。
假定4,当环境温度低于某一临界值 Tc 时,树木的形成层细胞停止活动,处于休眠状态,树木停止生长,此时不计降雨量对树木生长的影响。将当年春天超过 Tc 以上之日至夏天达到全年最高温度之日定为早材生长期; 将该日起至冬天气温降到Tc 以下之日定为晚材生长期,其余为休眠期。通常,将休眠期合并在晚材生长期内。
1.2 数学模型 1.2.1 温度的数学模型由假定2可知,以正弦波动的周期函数表示环境温度T随时间t的变化,可写成:
$ T(t)={{T}_{m}}+\vartriangle T\sin [\frac{2\pi }{P}(t+{{t}_{0}})]。 $ | (1) |
表 1是湖南省常德市花岩溪林场在1980 —1989年间月平均气温的数据(此数据来源于中国气象科学数据共享服务网),以此作为建立温度数学模型的基础数据。
由表 1可见,若用月作为时间单位,以此来表示气温年度变化曲线,由于数据量太少致使曲线不光滑。为绘制光滑曲线需增加数据量,对表 1中数据进行处理,采用三次样条插值法对每年按月平均气温进行插值,每月插入5点。相应地,式(1)中的周期P变为72,初相位t0变为50 。由此,式(1)改写为:
$ T(t)={{T}_{m}}+\vartriangle T\sin [\frac{2\pi }{36}(t+50)]。 $ | (2) |
用插值后数据与式(2)所示的数学形式进行拟合,结果如表 2 。
将表 2 中的拟合值代入式(2),得
$ T(t)=16.6+11.36\sin [\frac{2\pi }{36}(t+50)]。 $ | (3) |
将式(3)用图 3表示。由图可知,1年中温度的变化规律为: 年最低温度5.24 ℃,出现在1月20日; 年最高温度27.96 ℃,出现在7月20日; 年平均温度16.6 ℃。该地区的全年温度高于临界值 Tc,树木全年都在生长。由假定4可知,早材生长期为从1月1日—7月20日共计200天; 晚材生长期为7月21日—12月30日共计160天,早材期处在春季与初夏,晚材期处在后夏、秋季与冬季。
根据式(3)计算出月平均温度,将它与湖南省常德市花岩溪林场1990 —1999年间实测的月平均温度进行比较,结果绘于图 4中,将二者进行线性回归分析,决定系数达0.994,说明可用该温度模型准确预测此地区的环境温度变化情况。
降雨量的分段阶梯函数其数学形式为:
$ W(t)=\left\{ \begin{matrix} {{W}_{1}}\ \ \ t={{t}_{1}},\ \begin{matrix} \vdots \ {{W}_{n}}\ \ \ t={{t}_{n}} \ \end{matrix} \ \end{matrix} \right.。 $ | (4) |
表 3是湖南常德花岩溪林场1975 —2010年间月平均降雨量的数据(此数据来源于中国气象科学数据共享服务网),其中,权重系数是指月降雨量占全年降雨量的比重,日均降雨量通过月均降雨量除以30得到。将日均降雨量数据作为构建降雨量数学模型的基础数据。
当权重系数小于0.05时,认为降雨量对树木生长的影响甚微,对应的月份降雨量可以忽略; 当权重系数大于0.1时,认为降雨量对树木生长的影响一致; 当权重系数介于0.05与0.1之间时,降雨量与树木生长具有正相关关系。图 2中的分段阶梯函数在自变量时间区域内是不连续的,为便于计算,以表 3中数据为依据,将图 2中的分段阶梯函数改写成分段线性函数式(5),即降雨量数学模型:
$ W(t)=\left\{ \begin{matrix} 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [0+360n\ \ 30+360n); \\ \begin{matrix} \begin{matrix} \begin{matrix} 5.92(t-360n-30)/60\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [30+360n\ \ 90+360n);\ \ \ \ \ \ \ \ \\ 5.92\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [90+360n\ \ 210+360n); \\ \end{matrix} \\ 5.92-1.89(t-360n-210)/30\ \ \ \ \ \ \ [210+360n\ \ 270+360n); \\ \end{matrix} \\ 2.14-(t-360n-270)/30\ \ \ \ \ \ \ \ \ \ \ \ [270+360n\ \ 330+360n); \\ \end{matrix} \\ 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [330+360n\ \ 360+360n) \\ \end{matrix} \right. $ | (5) |
上述针对外部环境因素(温度与降雨量)建立的数学模型可作为量化分析树木生长的外部条件,为建立树木生长速率模型和早晚材密度分布模型需要引入以下假定。
2.1 物理模型假定1,树干以髓心为轴沿弦向对称生长,如图 6所示。树干的横切面是直径为D的圆,圆的半径为R。生长轮是以髄心为圆心、以不同r为半径的同心圆。树龄为N年的树干直径D可表示为:
$ D=2R=2({{r}_{0}}+\sum\limits_{i=1}^{N}{\vartriangle {{r}_{i}}}) $ | (6) |
假设原木为圆柱体,密度沿径向变化,在弦向与长度方向上均匀分布。这一假定意味着忽略了一些影响生长轮形成的因素,如树木光照的不对称性、株间距对树木生长的影响等。
假定2,在快速生长期内,树干直径随年份的增加呈线性增大的趋势(胡云云等,2009),这里忽略了树木老年期生长速率逐年下降的影响。树干直径D在其生长周期t内按线性增长,树干直径增长率为常数c0,可表达为:
$ \frac{dD}{dt}={{c}_{0}}。 $ | (7) |
用差分形式将式(7)改写为:
$ \begin{align} & \frac{{{D}_{N}}-{{D}_{N-1}}}{\vartriangle t}=\frac{2({{R}_{N}}-{{R}_{N-1}})}{\vartriangle t}=\frac{2\vartriangle {{r}_{N}}}{\vartriangle t}={{c}_{0}} \ &(N=2,3\ldots)。 \ \end{align} $ | (8) |
假定3,每年气候的变化规律是相同的,即年复一年的规律性周期变化; 且树木在一年内新生成一层木材,使生长轮与年份一一对应。这里忽略了一年内会形成多个生长轮的情况。
假定4,环境温度高,降雨量充足,树木生长快;温度低,降雨量少,树木生长慢。即温度与降雨量对树木生长有交互作用,这与 Yasue等(1995)、Kim(1995)和刘刚等(2009)研究结果一致。
假定5,树木生长速率不仅与环境温度有关,还与温度变化率有关。春季和秋季温度大致相同,但温度随时间的变化趋势相反。此外,当环境温度T低于某一临界值 Tc 时,木材形成层细胞停止活动,处于休眠状态,树木停止生长。
假定6,在一个生长轮内,木材的微观密度ρ(r)随树木生长速率的增加而降低,即早材期的树木生长快、密度小,晚材期的树木生长慢、密度大。
2.2 数学模型 2.2.1 树木生长速率的数学模型根据上述假定2,若以年为时间单位,树干直径随时间的增长率为常数; 而在一年内,以天或月为时间单位,树干直径随时间的生长速率不再是常数。根据上述假定4和5,树木生长速率(dr/dt)与环境温度(T)、降雨量(W)及温度变化率(dT/dt)线性相关,可表示为:
$ \frac{\text{d}r}{\text{d}t}={{c}_{1}}T(t)+{{c}_{2}}W(t)+{{c}_{3}}\frac{\text{d}T(t)}{\text{d}t}。 $ | (9) |
$ t=0,\ r={{r}_{0}}。 $ | (10) |
由式(9)与式(10)构成的常微分方程初值问题可用数值方法求解,为此,将式(9)写成差分形式如下:
$ \frac{{{r}^{n+1}}-{{r}^{n}}}{\vartriangle t}={{c}_{1}}{{T}^{n}}+{{c}_{2}}{{W}^{n}}+{{c}_{3}}\frac{{{T}^{n+1}}-{{T}^{n}}}{\vartriangle t}。 $ | (11) |
式(11)经整理可得:
$ \begin{align} & {{r}^{n+1}}={{r}^{n}}+({{c}_{1}}{{T}^{n}}+{{c}_{2}}{{W}^{n}})\vartriangle t+{{c}_{3}}({{T}^{n+1}}-{{T}^{n}}),\ & n=1,2\ldots \ \end{align} $ | (12) |
根据式(10)可得:
$ {{r}^{1}}={{r}_{0}}。 $ | (13) |
上述式(9)中的 c1,c2与 c3是3个可调整的系数,它们与树木生长环境(温度、降雨量)、树木种类等因素有关。
2.2.2 早晚材密度分布的数学模型根据上述假定6,木材的微观密度ρ(r)与树木的生长速率(dr/dt)呈反比,表示为:
$ \rho(r)={{c}_{4}}+{{c}_{5}}\frac{\text{d}r}{\text{d}t}。 $ | (14) |
式(14)的差分形式为:
$ {{\rho }^{n+1}}={{c}_{4}}+{{c}_{5}}\frac{\vartriangle t}{{{r}^{n+1}}-{{r}^{n}}}。 $ | (15) |
由式(15)可见,木材的微观密度ρ(r)是随树干半径变化的单值函数。c4和c5是2个可调整的系数,c4为木材髓心处的微观密度,c5与树木生长速率有关。
3 算例分析从湖南省常德市花岩溪林场采伐杉木(Cunninghamia lanceolata)5株,在胸径处沿树干纵向截取长10 cm 的圆盘,分别从各株髓心至树皮截取无瑕疵试件6个,尺寸为长度 L = 5 mm,宽度R= 165 mm,高度T= 10 mm,如图 7所示,试件在常规干燥条件下干燥,烘至绝干后密封。采用软 X 射线密度分析仪(JOYCE 3CS)测试木材密度沿径向的分布。
图 8是根据湖南省常德市花岩溪林场的气候条件由树木生长速率模型计算得到的杉木快速生长期内的生长曲线。由图 8可知,树干半径(实线)随时间呈波动式上升趋势。若以年为单位时间,将树干半径与时间进行线性回归分析,得回归曲线(虚线),决定系数为0.999 。说明若以年为单位时间,树干半径生长是线性的,这与胡云云等(2009)得出的针叶材和阔叶材的胸径与树龄呈正相关关系的结论相同。此外,树木的生长速率(dr/dt)随时间t呈正弦规律的周期变化特征(双点画线),在一年内,早材期的树木生长速率逐渐上升,晚材期的树木生长速率逐渐下降。
图 9是由早晚材密度分布模型计算得到的杉木快速生长期内的密度分布曲线。从图 9可以看出,木材密度随着树木生长速率的增加而降低,曲线的最低点为早材的最低密度,曲线的最高点为晚材的最大密度。
分别截取图 8与图 9中的一部分,即一个生长轮内树木生长速率、树干半径与木材密度的变化,分别表示在图 10 a 与图 10 b 中。由图 10 a 中树干半径随时间的变化曲线r(t)(实线)可知,树木在全年生长过程中没有休眠期 ; 从1月20起,由于环境温度上升、降雨量增加,树木生长速率迅速加快,表现为图中的虚线陡然上升 ; 同时,树干半径也迅速增加,此阶段为早材生长期。树木生长速率在7月6日达到最大值,此时树干半径生长了5 mm,木材密度为全年密度最小值(153 kg·m-3)(图 10 b); 此后,树木生长速率逐渐降低,树干半径的增加减缓,木材密度逐渐增大,直至密度达到最大值(594 kg·m-3)。从图 10 b 中可知,由于早材和晚材的生长速率不同,它们在同一生长轮内的宽度也不同。
图 11是密度分布实测值与模型计算值的比较。图中的4个样本从30个样本中随机选取,密度实测值取自从髓心至第11个生长轮之间,此阶段可视为杉木的快速生长期。从比较的结果看,实测值与模型计算值在生长轮宽度与早晚材密度分布这2方面基本一致。此外,早晚材密度分布模型还可以准确反映木材的最大密度、最小密度、平均密度及生长轮宽度(表 4)。
本文假定了树木生长环境(温度与降雨量)影响树木生长速率,树木生长速率影响木材早晚材密度。针对这种假定给出了定量描述,利用正弦波函数、分段线性函数和一阶常微分方程构建了树木生长速率模型和早晚材密度分布模型。并以湖南省常德市花岩溪林场人工林杉木为例,一方面,利用软 X 射线密度分析仪实测了杉木密度的径向分布; 另一方面,结合该地区的气象数据利用早晚材密度分布模型预测了杉木早晚材密度沿径向的分布规律。将密度的试验测量值与模型计算值进行比较,对早晚材密度分布模型进行验证。得出以下结论:
1)利用正弦波函数拟合了树木生长环境中环境温度的周期性变化,函数拟合的决定系数为0.907,正弦波函数可以准确反映环境温度的波动; 利用分段线性函数定量描述树木生长环境中降雨量的周期性变化。
2)树木生长速率模型可用于定量分析温度与降雨量对树木生长速率的影响。
3)早晚材密度分布模型准确描述了木材密度沿径向的分布规律,为木材宏观密度与微观密度的内在联系给出了定量的表述。该模型对建立诸如木材干燥、人造板热压等木材加工过程中木材宏观密度变化与内部微观密度分布的变化之间的联系提供理论依据。
[1] | 胡晓龙.2003.长白落叶松林分断面积生长模型的研究.林业科学研究, 16(4): 449-452.(1) |
[2] | 胡云云,亢新刚,赵俊卉.2009.长白山地区天然林林木年龄与胸径的变动关系.东北林业大学学报, 37(11): 38-42.(2) |
[3] | 刘刚,陆元昌,李晓慧,等.2009.六盘山地区气候因子对树木年轮生长的影响.东北林业大学学报,37(4): 1-4.(1) |
[4] | 李坚.2002.木材科学.2版.北京: 高等教育出版社, 184-185.(1) |
[5] | 吕建雄,林志远,赵有科,等.2005.杉木和I-72杨人工林木材干缩性质的研究.林业科学,41(5): 127-131.(1) |
[6] | 吕建雄,王金林,黄安民.2009.中国杨树木材加工利用研究进展.第二届中国林业学术大会-S11,木材及生物质资源高效增值利用与木材安全论文集,广西南宁, 233-241.(1) |
[7] | 徐金梅,吕建雄,鲍甫成,等. 2011.祁连山青海云杉木材密度对气候变化的响应.北京林业大学学报,33(5): 115-121.(1) |
[8] | 尹思慈.1996.木材学.北京: 中国林业出版社, 19.(1) |
[9] | 俞昌铭.2011.多孔材料传热传质及其数值分析.北京: 清华大学出版社, 13.(1) |
[10] | 张建国,段爱国.2004.理论生长方程与直径结构模型的研究.北京: 科学出版社, 2-7.(1) |
[11] | Anta M B,Dorado F C,Aranda U D,et al. 2006.Development of a basal area growth system for maritime pine in northwestern Spain using the generalized algebraic difference approach.Canadian Journal of Forest Research,36(6): 1461-1474.(1) |
[12] | Blair R B,Zobel B,Hitchings R C,et al.1976.Pulp yield and physical properties of young loblolly pine with high density juvenile wood.Applied Polymer Symposium,28: 435-444.(1) |
[13] | Ferrand J C.1982.Considerations on wood density,Part 2.Calculation of density and its heterogeneity.Holzforschung,36(3): 153-158.(1) |
[14] | Fritts,H C.1976.Tree rings and climate.Academic Press,New York,567.(1) |
[15] | Kim E S.1995.The growth of Korean fir trees(Abies koreana Wilson)affected by the changes of soil moisture content on Mt.Halla,Cheju island,Korea.International Workshop on Asian and Pacific Dendrochronology,Tsukuba Japan,58-63.(4) |
[16] | Lü J X,Urakami H,Hirkawa Y,et al.2003.Pilot Study to Examine the Radial Variation in Annual Ring Width,Density and Shrinkage.Chinese Forestry Science and Technology, 2(4): 21-26.(1) |
[17] | Pang S.1998.Relative importance of vapor diffusion and convective flow in modeling of softwood drying.Drying Technology, 16(1/2): 271-281.(1) |
[18] | Pernestål K,Jonsson B,Larsson B.1995.A simple model for density of annual rings.Wood Science and Technology, 29(6): 441-449.(1) |
[19] | Perre P,Moser M,Martin M.1993.Advances in transport phenomena during convective drying with superheated steam or moist air.International Journal of Heat and Mass Transfer,38(11): 2725-2746.(1) |
[20] | Rennols K,Geary D N,Rollinson T J D.1985.Characterizing diameter distributions by the use of the Weibull distribution.Forestry,58(1): 57-66.(1) |
[21] | Schweingruber F H.1988.Tree ring: Basics and applications of dendrochronology.Kluwer Academic Publishers,London,20-83.(1) |
[22] | Sprague J R,Talbert J T,Jett J B,et al.1983.Notes: Utility of the pilodyn in selection for mature wood specific gravity in loblolly pine.Forest Science, 29(4): 696-701.(1) |
[23] | Yasue K,Funada R,Fukazawa K.1995.Effect of climate factors on the maximum density of Picea glehnii.International Workshop on Asian and Pacific Dendrochronology,Tsukuba Japan,58-63.(3) |
[24] | Yu C M,Dai C P,Wang J H.2007.Heat and mass transfer in wood composite panels during hot pressing: Part 3.Predicted variations and interactions of the pressing variables.Holzforschung,61(1): 74-82.(1) |
[25] | Zobel B J,Buijtenen J P.1989.Wood variation: its causes and control.Springer-Verlag,New York,40-63.(2) |