CN112100876A - 一种纸页结构建模及抗张强度模拟的方法 - Google Patents

Abstract

本发明公开了一种纸页结构建模及抗张强度模拟的方法，步骤如下：通过纵向分段、截面离散以及改进的质点‑弹簧模型建立造纸纤维静态几何模型及造纸纤维特性模型；基于对纤维形态因子、纸页结构模型的三维尺寸及克重定义，生成纸页结构模型；使用分段方向包围盒检测纸页结构模型中造纸纤维特性模型之间的接触位置；使用结合弹簧模拟造纸纤维间结合力；最后通过恒速拉伸法，模拟纸页结构模型的抗张强度，得到与实际测量结果及理论经验趋势相符的模拟结果。与有限元方法相比，本发明中的方法具有速度快、复杂度低的特点，为相关研究工作，特别是在新型纤维材料开发工作中，研究新型纤维材料特性与产品性能之间的关系，提供了一种便捷有效的新方法。

Description

一种纸页结构建模及抗张强度模拟的方法

技术领域

本发明涉及制浆造纸及计算机模拟技术领域，具体涉及一种纸页结构建模及抗张强度模拟的方法。

背景技术

在制浆造纸领域，纸页抗张强度作为纸张最重要的性能之一被广泛关注，而纸页抗张强度主要依赖于其原材料，即造纸纤维的特性。为了表征纸页抗张强度特性，当前被广泛认可及使用的是半理论半经验的Page模型。然而，由于造纸纤维特性较多，不同特性之间存在相互影响，并且纸页抗张强度性能往往由多个造纸纤维特性及抄造工艺共同决定，Page模型虽然能较好的通过造纸纤维的部分特性及抄造工艺值预测纸页抗张强度，但难以分析造纸纤维特性与纸页抗张强度间的关系，特别是机理关系分析。因此，近年来有研究报道使用计算机模拟的方法，建立纸页结构模型并对其进行抗张强度模拟。

在已有报道中，为了提高计算机模拟速度，有一类建模模拟方法，通常使用极简化的造纸纤维特性模型模拟纸页结构及抗张强度，例如使用球-棍模型、编码模型等简化方法。这类模拟造纸纤维的模型虽然简单易用且模拟速度较快，但由于模型与真实造纸纤维的形态差异较大，所以较难真实反映造纸纤维特性与纸页抗张强度性能间的关系。另一类建模模拟方法为有限元方法，虽然有限元方法可以较为真实的模拟造纸纤维、纸页结构以及纸页抗张强度性能，但由于有限元方法在模拟计算过程中结合了各种物理场函数及函数插值等方法，导致其模拟速度较慢。

发明内容

本发明的目的在于弥补现有技术对纸页结构建模及抗张强度模拟过程中，真实性较差或模拟速度较慢的问题，基于使用改进质点-弹簧模型对造纸纤维建模的方法，提供一种纸页结构建模及抗张强度模拟的方法。该方法能在加快模拟速度的同时尽可能真实的反映造纸纤维类型及特性对纸页抗张强度的影响，从而为纸页结构建模与性能模拟的研究工作，以及在新型纤维材料开发中，研究新型纤维材料特性与产品物理强度性能之间的关系，提供了一种便捷有效的新方法，可以提供各种条件下，甚至极端条件下的模拟结果，有利于加快研究进程。

本发明的目的可以通过采取如下技术方案达到：

一种纸页结构建模及抗张强度模拟的方法，所述方法包括以下步骤：

S1、使用改进质点-弹簧模型建立造纸纤维特性模型，过程如下：

S11、定义造纸纤维特性参数作为建立造纸纤维特性模型的输入参数，

S12、基于步骤S11中定义的纤维长度L、纤维宽度W、纤维壁厚WT建立造纸纤维静态几何模型，

S13、使用改进质点-弹簧模型建立造纸纤维特性模型，

S14、计算改进质点-弹簧模型参数；

S2、生成纸页结构模型，过程如下：

S21、定义纤维形状因子及计算造纸纤维特性模型中相邻截面法向量偏转角度，

S22、定义纸页结构模型的尺寸及纸页克重，

S23、模拟生成纸页结构模型；

S3、采用分段OBB包围盒检测纸页结构模型中造纸纤维特性模型之间的接触区域，过程如下：

S31、使用二叉树结构分解造纸纤维特性模型，

S32、计算二叉树结构中每个造纸纤维特性模型局部的OBB包围盒，

S33、计算纸页结构模型中造纸纤维特性模型之间的接触区域；

S4、采用结合弹簧模型模拟纤维间结合力；

S5、使用恒速拉伸法模拟纸页结构模型的抗张强度。

本发明相对于现有技术具有如下的优点及效果：

本发明使用了一种基于改进质点-弹簧模型构建的造纸纤维特性模型，建立纸页结构模型并对其进行抗张强度模拟。与以往使用简化造纸纤维特性模型建立纸页结构模型并模拟纸页抗张强度性能的方法相比，本发明方法具有真实度更高的优点；与另一种使用较多的有限元模拟方法相比，本发明方法具有模拟速度较快，且消耗计算机内存较少的优点。

通过使用本方法，可以很好的模拟不同类型造纸纤维及不同造纸纤维特性对纸页抗张强度性能的影响趋势，从而为纸页结构建模与性能模拟的研究工作，以及在新型纤维材料开发中，研究新型纤维材料特性与产品物理强度性能之间的关系，提供了一种便捷有效的新方法，可以提供各种条件下，甚至极端条件下的模拟结果，有利于加快研究进程。

附图说明

图1是本发明实施例中公开的纸页结构建模及抗张强度模拟方法流程图；

图2是造纸纤维静态几何模型示意图；

图3是改进质点-弹簧模型建立造纸纤维特性模型示意图；

图4是造纸纤维特性模型分段受到拉伸力时的形变及受力分析图；

图5是造纸纤维特性模型分段受到扭转力时的形变及受力分析图；

图6是造纸纤维形状因子特性描述示意图；

图7是造纸纤维特性模型相邻截面法向量偏转角度定义示意图；

图8是造纸纤维特性模型实现接触区域检测的分段OBB包围盒结构示意图。

图9是两根造纸纤维特性模型之间的结合弹簧BS受到拉伸力时的形态示意图；

图10是两根造纸纤维特性模型之间的结合弹簧BS受到剪切力时的形态示意图；

图11是基于恒速拉伸法模拟纸页结构模型抗张强度的示意图；

图12是通过恒速拉伸法模拟纸页抗张强度所得到的抗张力-伸长量关系曲线示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

本实施例公开了一种纸页结构建模及抗张强度模拟方法，包括5个步骤，如图1所示，具体描述如下：

S1、使用改进质点-弹簧模型建立造纸纤维特性模型：

为了实现对纸页结构的建模，需要首先建立纸页结构模型的基本组成模型，即造纸纤维特性模型，造纸纤维特性模型的具体建模方法如下：

S11、定义造纸纤维特性参数作为建立造纸纤维特性模型的输入参数：

根据具体模拟需要，被定义的造纸纤维特性参数包括纤维长度L、纤维宽度W、纤维壁厚WT、纤维粗度c、纤维弹性模量E、纤维泊松比v以及纤维剪切模量G。

S12、建立造纸纤维静态几何模型：

由于计算机较难解决连续问题(包括连续长度、连续时间等)，所以需要对真实造纸纤维形态进行离散处理并建模。考虑到真实造纸纤维形态属于空心束状结构，所以使用纵向分段，截面离散的方式建立造纸纤维静态几何模型，如图2所示。为了简化模型，将造纸纤维静态几何模型沿纵向均匀分段，即定义所有分段长度一致，为l_seq；造纸纤维特性模型的截面定义为同心圆环，但截面离散点仅用单层表示圆形轮廓，圆形轮廓直径为Φ＝W-2·WT。定义造纸纤维静态几何模型纵向共有N_seq个分段，则有N_seq+1个截面。记录第i(0≤i≤N_seq)个截面的中心点位置为P_(i)，截面法向量为

通过上述定义可得，造纸纤维静态几何模型中第i个分段由序号为i-1和i的截面定义，P_(i)相对于P_(i-1)的位置由

和分段长度l_seq计算得到，如式(1)所示：

定义造纸纤维特性模型中的每一个截面被离散为N_dis个点，可定义三维空间中截面离散点的坐标为：p_(i)(j)＝(x_(i)(j),y_(i)(j),z_(i)(j))，0≤i≤N_seq；1≤j≤N_dis)，表示第i个横截面上的第j个离散点，x_(i)(j),y_(i)(j)和z_(i)(j)分别表示离散点p_(i)(j)在空间坐标系中对应x、y和z轴的坐标值。

本实施例中，序号为0的截面上的离散点可取圆周上的任意位置为起点，等距离定义截面其它离散点位置坐标；后续截面上的点可由前一个截面上的离散点根据截面中心点坐标位置进行平移，再基于截面法向量的方向进行旋转计算得到。

S13、使用改进质点-弹簧模型建立造纸纤维特性模型：

对于使用传统质点-弹簧模型建立造纸纤维特性模型，主要包含质点模型和弹簧模型的建立，其方法是将造纸纤维特性模型的轮廓和质量简化为质点模型，内部力使用弹簧模型进行模拟。其中，质点模型的坐标与步骤S12中建立的造纸纤维静态几何模型中的截面离散点坐标一致，即造纸纤维特性模型中第i个截面上的第j个质点模型坐标mp_(i)(j)＝p_(i)(j)。弹簧模型分为纵向结构弹簧模型(ST1)、截面结构弹簧模型(ST2)和剪切弹簧模型(SS)。其中，ST1用于连接造纸纤维特性模型前后截面中同一序号的质点模型，即连接mp_(i)(j)与mp_(i+1)(j)和mp_(i)(j)与mp_(i-1)(j)，其中mp_(i+1)(j)和mp_(i-1)(j)分别代表造纸纤维特性模型中第i+1个截面和第i-1个截面上的第j个质点模型；ST2用于连接造纸纤维特性模型的同一截面中相邻质点模型，即连接mp_(i)(j)与mp_(i)(j+1)和mp_(i)(j)与mp_(i)(j-1)，其中mp_(i)(j+1)和mp_(i)(j-1)分别代表造纸纤维特性模型中第i个截面上的第j+1个和第j-1个质点模型；SS用于连接造纸纤维特性模型中，相邻截面上相邻位置的质点模型，即连接mp_(i)(j)与mp_(i+1)(j+1)和mp_(i)(j)与mp_(i+1)(j-1)，其中，mp_(i+1)(j+1)和mp_(i+1)(j-1)分别代表造纸纤维特性模型中第i+1个截面上的第j+1个和第j-1个质点模型。

由于使用上述传统质点-弹簧模型模拟造纸纤维极易导致截面崩溃，从而难以获得模拟结果，所以本实施例通过在传统质点-弹簧模型的基础上，添加中心弹簧模型CS以保证造纸纤维特性模型截面强度。改进的质点-弹簧模型结构如图3所示，CS用于连接造纸纤维特性模型中，同一截面上位置相对的一组质点模型，即连接mp_(i)(j)与

其中，

代表造纸纤维特性模型中第i个截面上的第

个质点模型。

S14、计算改进质点-弹簧模型参数：

首先计算改进质点-弹簧模型中，质点模型的质量。为了简化计算过程，假设造纸纤维特性模型的质量均匀分布，则可以通过式(2)计算质点模型的质量：

上式中，c和L分别表示用于建立造纸纤维特性模型的纤维粗度和纤维长度参数，其乘积表示造纸纤维特性模型的总质量；m_(i)(j)表示造纸纤维特性模型中第i个截面上第j个质点模型的质量；m_(0)(j)和

分别表示造纸纤维特性模型中，序号为0和N_seq截面上的质点模型质量，由于这些质点模型只属于造纸纤维特性模型中的一个分段，而造纸纤维特性模型中间截面的质点模型同时属于前后两个分段，所以在序号为0和N_seq截面上的质点模型，其质量为中间截面上质点模型质量的一半。

弹簧模型的弹性系数通过假设造纸纤维特性模型分段受到拉伸力或扭转力时的受力分析进行求解。如图4所示，是造纸纤维特性模型分段受到拉伸力作用发生形变时，形变量及对应质点-弹簧模型中任一质点模型的受力分析结果。定义造纸纤维特性模型分段的一个截面受到大小为F_s的拉伸力，另一个截面固定。根据步骤S11中对纤维弹性模量E和纤维泊松比v的设定，以及基于材料力学中对弹性模量及泊松比的定义，可以计算得到造纸纤维特性模型分段在拉伸力F_s作用下的纵向形变量和截面形变量，从而得到形变后各弹簧模型的长度及偏转角度，分别包括在受到拉伸力后，连接目标质点模型的纵向结构弹簧长度l′_ST1、截面结构弹簧长度l′_ST2、剪切弹簧长度l′_SS和中心弹簧长度l′_CS，以及剪切弹簧在纵向受力分析中弹簧力与水平方向的夹角α和在截面受力分析中弹簧力与垂直方向的夹角β；受力截面中每个质点模型受到的拉伸力f_s可由截面受到拉伸力F_s对截面总质点模型数求平均得到。最终，列出截面质点模型的受力平衡方程如式(3)所示：

其中：

上式中，f_ST1、f_ST2、f_SS和f_CS分别代表纵向结构弹簧、截面结构弹簧、剪切弹簧以及中心弹簧的弹簧力；k_ST1、k_ST2、k_SS和k_CS分别代表纵向结构弹簧、截面结构弹簧、剪切弹簧以及中心弹簧的弹性系数；l_ST1、l_ST2、l_SS和l_CS分别代表纵向结构弹簧、截面结构弹簧、剪切弹簧以及中心弹簧的原始长度，等于该弹簧模型所连接的两个质点模型在受到拉伸力前的距离。

如图5所示，是造纸纤维特性模型分段受到扭转力F_T时的质点模型位移及受力分析示意图。由材料力学理论可知，造纸纤维截面所受到扭转力会转换为质点模型受到切向力。定义每个质点模型受到切向力大小为f_T，则基于步骤S11中对剪切模量G的设定，以及基于材料力学中对剪切模量的定义，可以计算得到每个质点模型在切向力f_T作用下的偏转角度γ，从而得到形变后各弹簧的长度及偏转角度，包括连接目标质点模型的纵向结构弹簧受到切向力后的长度l″_ST1，左、右两侧剪切弹簧受到切向力后的长度l″_SS和l″′_SS，以及连接目标质点模型的纵向结构弹簧及左、右两侧剪切弹簧受到切向力后的弹簧力在纵向受力分析时与水平方向的夹角θ″_ST1、θ″_SS和θ″′_SS；。

另外，当造纸纤维特性模型截面受到扭转力作用时，虽然可能伴随纵向长度形变及截面形变，但由于形变量较小，可忽略不计；截面受力分析中，质点模型所受到的纵向结构弹簧力f″_ST1和剪切弹簧力f″_SS与f″′_SS，在截面上的投影力方向与受力正交分解后水平方向的夹角较小(本实施例中，上述三个力与水平方向的夹角约为5°～15°)，其余弦函数值约为1，为简化计算，这些夹角忽略不计。因此，得到受力平衡方程如式(4)所示：

f_T＝f″_SScos(θ″_SS)+f″′_SScos(θ″′_SS)+f″_ST1cos(θ″_ST1)

其中：

考虑到纵向结构弹簧ST1和截面结构弹簧ST2均属于结构弹簧，具有相同性质，即相同弹性模量。根据材料力学中弹性模量的定义和胡克定律，可以得到两种结构弹簧间弹性系数的关系如式(5)所示：

上式中，E_f为两种结构弹簧的弹性模量；F和S分别表示纵向结构弹簧及截面结构弹簧所受到的拉伸力和受力面积，Δl_ST1和Δl_ST2表示纵向结构弹簧和截面结构弹簧受到拉伸力F后的形变量。为了便于计算，定义两种结构弹簧的受力大小F及受力面积S一致。

在式(3)-(5)中，除了k_ST1、k_ST2、k_SS和k_CS分别表示纵向结构弹簧、截面结构弹簧、剪切弹簧和中心弹簧的弹性系数，为待求未知参数外，其余参数均可根据设定的受力大小及材料力学中关于弹性模量、泊松比及剪切模量的定义计算得到。因此，将式(3)-(5)联立，可以求解k_ST1、k_ST2、k_SS和k_CS四个未知参数的值，即质点-弹簧模型中弹簧模型的弹性系数值。

S2、生成纸页结构模型：

S21、定义纤维形状因子及计算造纸纤维特性模型中相邻截面法向量偏转角度：

生成纸页结构模型时，需要考虑真实造纸纤维的弯曲形态，该弯曲形态由纤维形状因子定义。如图6所示，假设纤维弯曲后形成一个圆弧，其对应弧线段和弦线段分别为AOB和AO′B，其中，A和B为弧线段及弦线段的两个端点，O和O′分别为弧线段和弦线段的中点，O同时为弧线段的最高点，则纤维形状因子C_f的定义如式(6)所示：

通过对纤维形状因子的定义，可以基于简单几何方法计算得到∠OAB的角度

进一步与造纸纤维特性模型对应，如图7所示，在一个分段数为N_seq的造纸纤维特性模型中，序号为0和

的截面法向量之间的偏转角度等于图6中的

定义造纸纤维特性模型中所有相邻截面法向量偏转角度大小一致，则可以通过式(7)计算得到两个相邻截面法向量的偏转角度θ_N。

S22、定义纸页结构模型的尺寸及纸页克重：

除了定义造纸纤维特性模型弯曲形态，在生成纸页结构模型前，还需要定义纸页结构模型的尺寸，包括：长l_p、宽w_p和厚度t_p，以及纸页克重BW。上述l_p、w_p、t_p和BW参数可以根据实际模拟需求及计算机运算能力进行定义。如运算能力强的计算机可以定义较大尺寸、较大克重的纸页结构模型参数，否则，定义较小尺寸、较小克重的纸页结构模型参数。

S23、模拟生成纸页结构：

在纸页结构的生成过程中，主要基于如下步骤：

S231、采用随机位置、随机朝向的方法逐一生成并放置造纸纤维特性模型，并在每放置一根造纸纤维特性模型后，通过删除超出纸页结构模型的尺寸长l_p、宽w_p和厚度t_p范围的造纸纤维特性模型纵向分段，以满足纸页结构模型的尺寸要求；

S232、在步骤S231不断添加造纸纤维特性模型的过程中，造纸纤维特性模型的纸页克重逐步提升。通过计算当前纸页结构模型的克重，判断是否已达到定义值BW的要求，即与BW的定义值相比，偏差不超过5％。如果未达到要求，则通过步骤S231继续添加造纸纤维特性模型；否则，停止添加造纸纤维特性模型。

S233、在通过步骤S231和S232构建了符合纸页结构模型尺寸和纸页克重要求的成纸结构模型后，进一步计算当前纸页结构模型中所有造纸纤维特性模型的平均长度，以判断是否满足步骤S11中对纤维长度L的定义要求，即与L的定义值相比，偏差不超过5％。如果纸页结构模型中造纸纤维特性模型的平均长度不满足定义值L的要求，可以通过逐个随机添加较长或较短的造纸纤维特性模型的方法，提高或降低造纸纤维特性模型的平均长度，满足定义值L的要求。其中，造纸纤维特性模型长度的取值范围为定义值L的0.5-1.5倍。如果在上述添加造纸纤维特性模型的过程中，导致纸页结构模型的纸页克重超出了定义值BW的要求，则需要逐个随机删除步骤S231添加的造纸纤维特性模型，直到满足纸页克重BW的定义值要求。通过反复调节，使最终生成的纸页结构模型同时满足定义值BW和L的要求。

S3、采用分段OBB包围盒检测纸页结构模型中，造纸纤维特性模型之间的接触区域：

在纸页结构中，造纸纤维间结合力大小对纸页抗张强度影响较大，所以需要对纸页结构模型中的造纸纤维特性模型间结合力进行模拟。为了实现对造纸纤维特性模型间结合力的模拟，需要首先采用碰撞检测算法检测造纸纤维特性模型之间的接触位置。基于造纸纤维特性模型形态特点，本实施例使用分段方向包围盒(Oriented Bounding Box,OBB)方法检测造纸纤维特性模型之间的接触区域，具体步骤如下：

S31、使用二叉树结构分解造纸纤维特性模型：

如图8所示，为了提高OBB包围盒碰撞检测算法的计算准确度，需要首先使用二叉树结构分解造纸纤维特性模型，即构建分段OBB包围盒结构，其具体方法如下：

S311、将完整造纸纤维特性模型作为二叉树结构根节点；

S312、将造纸纤维特性模型沿纵向分段平均分为两个部分，并将这两个部分作为二叉树结构根节点的左、右子节点，插入到二叉树结构中；

S313、针对二叉树结构中每个新添加的子节点，将其所代表的造纸纤维特性模型局部再次沿纵向分段平分为两部分，添加为对应子节点的左、右子节点，重复该过程直到当前子节点只表示造纸纤维特性模型的一个分段，即无法再沿纵向平分为两个部分为止；

S314、对于当前子节点只表示造纸纤维特性模型一个纵向分段的情况，需要基于造纸纤维特性模型截面离散点进一步平分成两个部分，添加到对应子节点下作为新的左、右子节点，重复该过程直到当前子节点所代表的造纸纤维特性模型截面也无法再继续平分为止，即每个截面只包含2个离散点为止。

S32、计算二叉树结构中每个造纸纤维特性模型局部的OBB包围盒：

针对二叉树结构中单个节点所代表的造纸纤维特性模型局部，计算OBB包围盒，其方法如下：

S321、依次获取二叉树结构中每个节点所代表的造纸纤维特性模型局部，统计其表面轮廓点坐标，即当前造纸纤维特性模型局部中所包含的截面离散点坐标，定义当前获取的所有表面轮廓点坐标按x、y和z轴分别表示为X＝(x₁,x₂,...,x_n)，Y＝(y₁,y₂,...,y_n)，Z＝(z₁,z₂,...,z_n)，其中n表示当前获取的所有表面轮廓点坐标个数，X，Y和Z分别表示所有表面轮廓点坐标在x，y和z轴上的坐标值集合，x_q、y_q、z_q(1≤q≤n)分别表示第q个表面轮廓点在x、y和z轴的坐标值。

S322、计算步骤S321中所获取的造纸纤维特性模型局部表面轮廓点坐标相对于x、y和z轴的协方差矩阵，协方差矩阵定义如式(8)所示:

上式中，C表示协方差矩阵，cov表示协方差函数，cov(X,X)，cov(X,Y)，cov(X,Z)，cov(Y,X)，cov(Y,Y)，cov(Y,Z)，cov(Z,X)，cov(Z,Y)以及cov(Z,Z)分别表示括号内第一个参数所代表的表面轮廓点坐标在x、y或z轴上的坐标值集合，相对与括号内第二个参数所代表的x、y或z轴上坐标值集合之间的协方差；以cov(X,Y)为例，协方差函数计算方法如式(9)所示：

上式中，

和

分别代表X和Y中所有坐标值的平均值。

S323、求步骤S322中所得到协方差矩阵C的特征向量。根据特征向量的定义，对于协方差矩阵C，若存在数λ和非零矩阵A，使得CA＝λA，则称λ为矩阵C的特征值，A为矩阵C的特征向量，列出表达式如式(10)所示：

上式中，E为单位矩阵。

分析式(10)可知，该行列式是一个关于特征值λ的一元三次方程，通常有3个解，如果最后得到的结果有2个甚至3个同值解，则说明步骤S321中所选择的造纸纤维特性模型局部表面轮廓点坐标全部共面，甚至共线。根据本实施例步骤S12中对造纸纤维静态几何模型的定义可知，步骤S321中所选择的造纸纤维特性模型表面轮廓点坐标不可能全部共面或共线，所以特征值λ一定能得到3个互不相同的解，记为λ₁、λ₂和λ₃。

将特征值λ₁、λ₂和λ₃分别带入特征向量定义式：CA＝λA，可以求到三个不同的特征向量A₁、A₂及A₃。

S323、计算造纸纤维特性模型的OBB包围盒参数。根据步骤C3中对协方差矩阵C求解得到的特征向量A₁、A₂及A₃，以及协方差矩阵特征向量在空间几何中的几何意义可知，特征向量A₁、A₂及A₃可以更好的描述步骤S321中所获取的所有轮廓点坐标的特征信息，并对这些轮廓坐标点特征进行正交分解。因此，可以使用特征向量A₁、A₂及A₃作为OBB包围盒的三维坐标轴，计算一个长方体包围盒模型，使该长方体包围盒模型的边与A₁、A₂或A₃平行，并且该长方体包围盒模型能刚好将步骤S321中获取的所有轮廓点坐标包围。将该长方体包围盒模型作为步骤S321中所选取的造纸纤维特性模型局部OBB包围盒，并记录该长方体包围盒模型在A₁、A₂及A₃方向上的长度参数以表征该OBB包围盒的形状特点。

S33、计算纸页结构模型中造纸纤维特性模型之间的接触区域：

通过步骤S32得到的二叉树结构中不同造纸纤维特性模型局部OBB包围盒，即分段OBB包围盒结构后，采用如下方法计算纸页结构模型中造纸纤维特性模型之间的相互接触区域：

S331、逐一遍历选取纸页结构模型中两根不同的造纸纤维特性模型；

S332、获取步骤S331中所选择的两根造纸纤维特性模型对应的分段OBB包围盒结构，并比较两个分段OBB包围盒结构的根节点OBB包围盒的位置关系，即判断根节点所表示的两个长方体包围盒模型之间是否包含重叠区域。如果不包含，则说明两根造纸纤维特性模型未发生接触，返回步骤S331，重新选择两根造纸纤维特性模型进行判断；如果包含重叠则继续进行后续步骤；

S333、依次比较两根造纸纤维特性模型的分段OBB包围盒结构中，根节点的左、右子节点所包含的OBB包围盒是否存在重叠关系。选择存在重叠关系的一对OBB包围盒所在节点，依次比较其左、右子节点所包含的OBB包围盒是否仍然存在重叠关系，重复该步骤，直到当前发生重叠关系的OBB包围盒所属二叉树结构节点已没有左、右子节点，则记录当前节点所代表的造纸纤维特性模型局部为造纸纤维特性模型之间的接触区域。若通过执行本步骤，已无法找到具有重叠关系的一对OBB包围盒，则返回步骤S331，重新选取两根造纸纤维特性模型进行步骤S332和步骤S333。

S4、采用结合弹簧模型模拟纤维间结合力：

在通过步骤S3计算得到造纸纤维特性模型的接触区域后，需要将相互接触的造纸纤维特性模型，在接触区域使用结合弹簧模型BS连接，BS的连接端点分别为两根造纸纤维特性模型接触区域包含的质点模型。考虑到在后续模拟过程中，不同造纸纤维特性模型可能在受到拉伸力或剪切力作用下分离，需要针对不同受力情况对造纸纤维特性模型在分离过程中纤维间结合力大小进行模拟。如图9和图10所示，分别为结合弹簧BS受到拉伸力和剪切力时的形态示意图，其中，斜线阴影区域表示造纸纤维特性模型之间的接触区域，面积分别为S_A和S_B，法向量分别为

和

用于区分BS受到拉伸力或剪切力的判断标准为：BS伸长方向与所对应造纸纤维特性模型接触区域的法向量

和

之间的平均夹角大小，若该平均夹角小于45°，则认为BS受到拉伸力，否则BS受到剪切力。

定义BS受到拉伸力和剪切力时的最大受力强度分别为T_N和T_S；BS的最大受力伸长长度为D_bond，则在不同受力情况下，结合弹簧模型的弹性系数如式(11)所示：

上式中，k_N和k_S分别代表结合弹簧模型BS在受到拉伸力或剪切力时的弹性系数。

当结合弹簧长度超过D_bond时，标记该结合弹簧模型断裂，后续模拟将不再计算该结合弹簧模型的形变量及对应弹簧力。

S5、使用恒速拉伸法模拟纸页结构模型的抗张强度

参考国标GB/T 22898-2008中测量纸和纸板抗张强度的恒速拉伸法，模拟预测纸页结构模型的抗张强度，具体方法如图11所示：固定纸页结构模型中最左侧一段距离l_left内的质点模型，给纸页结构模型中最右侧一段距离l_right内的质点模型添加恒速向右的速度v_tensile，动态模拟纸页结构模型被恒速拉伸的形变过程，计算与纸页结构模型中最右侧一段距离l_right内的质点模型相连的所有弹簧模型，包括：纵向结构弹簧模型ST1、截面结构弹簧模型ST2、剪切弹簧模型SS、中心弹簧模型CS，以及结合弹簧模型BS产生的总弹簧力，作为纸页结构模型形变过程的抗张力σ。

通过上述方法描述对纸页结构模型的恒速拉伸过程进行动态模拟，其模拟方法为：将动态模拟过程的连续时间离散为多个时间点，每两个相邻的时间点之间的间隔为Δt，并假设在两个相邻时间点之间，纸页结构模型中所有质点模型的受力不变。通过分析每个质点模型在离散时刻t的受力情况，可以运用牛顿运动定律计算每个质点模型在t+Δt时刻的位移量及新位置，为了提高计算精度，使用基于牛顿运动定律的龙格库塔法迭代计算每个质点模型在不同离散时间点的位置，直到纸页结构模型完全分离为两个部分。

记录纸页结构模型在形变过程中的抗张力σ，以及纸页结构模型中最右侧一段距离l_right内质点模型的移动距离，即伸长量δ，可以绘制出纸页结构模型的抗张力-伸长量关系曲线，如图12所示。该抗张力-伸长量关系曲线的最大值即表示纸页结构模型的最大抗张力，通过式(12)可换算为纸页结构模型的抗张强度。

上式中，T_p为纸页结构模型的抗张强度；max(σ)、δ和w_p分别为纸页结构模型的最大抗张力、最大抗张力所对应的伸长量以及纸页结构模型的宽度。

实施例二

使用不同类型的造纸纤维特性参数建立造纸纤维特性模型，并进行纸页结构建模及抗张强度预测，其步骤如下：

T1、选取五种性能差异较大的造纸纤维并分别获取造纸纤维特性参数，包括桉木纤维、针叶木纤维、棉纤维、蔗渣纤维和竹纤维。分别测量五种造纸纤维的特性参数。被测量的造纸纤维特性参数包括：纤维重均长度L_WA、纤维宽度W、纤维壁厚WT、纤维粗度c以及纤维弹性模量E。其中，纤维重均长度L_WA是指按重量加权后的造纸纤维平均长度，可以与纤维平均长度L进行替换，但是L_WA比L能更好的反映造纸纤维原料的平均长度特性。上述造纸纤维特性参数的测量方法如下：L_WA、W、WT和c可以使用纤维分析仪直接测量获得；E可以基于材料力学中的横向负载梁弯曲实验测量并计算获得，步骤如下：(1)利用纤维自身重量使其悬挂在被固定于测量台面的细金属丝上，结合材料力学挠曲方程计算造纸纤维柔软度Flex；(2)通过电镜扫描获得纤维截面图形并进行二值化处理，基于该图像及截面惯性矩定义，计算纤维截面惯性矩I；(3)由材料力学中对柔软度和刚度的定义可知，柔软度Flex为刚度的倒数，而刚度等于纤维弹性模量E与截面惯性矩I的乘积。因此可以使用前两步实验计算得到的Flex和I反推纤维弹性模量E。

最终得到用于纸页结构模型建模使用的不同类型造纸纤维特性参数，如表1所示：

表1.用于纸页结构模型建模的不同类型造纸纤维特性参数

除了上述造纸纤维特性参数外，还需要获得造纸纤维的泊松比v及剪切模量G特性参数。由于造纸纤维比较细小，常规实验较难测量泊松比v，所以参考现有纸页结构模拟的文献，取v＝0.2；剪切模量G可通过弹性模量E和泊松比v计算得到，公式如式(10)所示：

基于表1及上述对造纸纤维泊松比v和剪切模量G特性参数的定义，以及步骤S12-步骤S14中所述的建立造纸纤维特性模型的方法，建立不同类型造纸纤维特性模型。

为了验证模拟结果，本实施例使用不同类型造纸纤维原料及MESSMER 200纸页手抄机进行抄纸实验，再将抄造好的纸张使用L&W CE062抗张强度仪测量纸页抗张强度。由于实验参数与模拟参数较难达到统一，后续步骤主要对使用不同造纸纤维原料生成纸页结构的抗张强度大小趋势进行比较。通过实验实测得到的纸页结构抗张强度大小趋势为：桉木纤维>竹纤维>针叶木纤维>蔗渣纤维>棉纤维。

T2、参考文献报道，定义纸页结构中的纤维因子C_f＝0.95；结合弹簧BS受到拉伸力和剪切力时的最大受力强度分别为T_N＝36.36MPa和T_S＝4.00MPa，最大拉伸长度为D_bond＝1.1μm；纸页结构模型的长、宽、厚度分别为l_p＝4mm、w_p＝2mm及t_p＝0.057mm，纸页克重BW＝27g/m²。基于上述参数定义，使用步骤T1中所建立的不同造纸纤维特性模型，采用步骤S2-S4中所描述的方法，建立由不同类型造纸纤维特性模型所生成的纸页结构模型。

T3、按照步骤S5中所描述的模拟纸页结构模型抗张强度的方法，需要首先定义参数l_left、l_right和v_tensile的值。本实施例中，l_left和l_right的大小被定义为纸页结构模型长度l_p的10％，即l_left＝l_right＝0.4mm；v_tensile参考国标中测量纸和纸板抗张强度的恒速拉伸法(GB/T 22898-2008)，被定义为16.67mm/ms。基于步骤S5所介绍的方法对使用不同类型造纸纤维特性模型生成的纸页结构模型进行抗张强度模拟，得到不同类型造纸纤维特性模型对应纸页结构模型抗张强度的大小关系趋势为：桉木纤维(1.02N/m)>竹纤维(0.9kN/m)>针叶木纤维(0.81kN/m)>蔗渣纤维(0.80kN/m)>棉纤维(0.45kN/m)，这一趋势与步骤T1中所获得的实测纸页抗张强度的大小关系趋势一致，说明本发明中所描述的纸页结构建模及抗张强度模拟方法可以很好的区分不同类型造纸纤维特性参数的差异，以及对纸页抗张强度的影响。

另外，与相关报道中使用有限元方法模拟纸页结构及抗张强度的结果相比，本实施例方法的模拟速度提高了约7.2倍，内存消耗减少约2.6倍。

实施例三

以桉木纤维重均长度特性为例，模拟造纸纤维特性参数对纸页抗张强度的影响趋势，其步骤如下：

U1、基于表1中测量得到的桉木纤维特性参数，假设其重均长度分别为0.5mm、0.75mm、1.00mm、1.25mm和1.50mm，其余特性参数不变；另外，考虑到表1中实际测得的桉木纤维重均长度为0.64mm，也将该数据用于后续建模模拟。使用步骤S1中的方法生成不同重均长度的桉木纤维特性模型。

U2、基于步骤T2中定义的所有参数值，使用步骤U1中生成的不同重均长度的桉木纤维特性模型，依据步骤S2-S4中描述的方法，建立由不同重均长度的桉木纤维特性模型所生成的纸页结构模型。

U3、基于步骤T3中定义的所有参数值，以及步骤S5中描述的方法，模拟不同重均长度的桉木纤维特性模型所生成纸页结构模型的抗张强度，模拟结果如表2所示：

表2.不同重均长度的桉木纤维生成纸页结构模型的抗张强度

表2中，纸页结构模型的抗张强度随着桉木纤维重均长度特性的增加而减小，与纸浆造纸行业现有经验认知基本一致；另外，考虑到造纸纤维的重均长度特性主要与抄造工艺中的打浆工艺相关，打浆度越高，造纸纤维的重均长度越短，所以表2中的模拟结果还可以解释为纸页结构模型的抗张强度随着打浆度的增大而增大，这也纸浆造纸行业现有经验认知基本一致，说明本发明描述的纸页结构建模及抗张强度模拟方法可以较好的模拟不同造纸纤维特性及纸页抄造工艺对纸页抗张强度的影响趋势。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (6)

1.一种纸页结构建模及抗张强度模拟的方法，其特征在于，所述方法包括以下步骤：

S1、使用改进质点-弹簧模型建立造纸纤维特性模型，过程如下：

S11、定义造纸纤维特性参数作为建立造纸纤维特性模型的输入参数，所述造纸纤维特性参数包括纤维长度L、纤维宽度W、纤维壁厚WT、纤维粗度c、纤维弹性模量E、纤维泊松比v以及纤维剪切模量G；

S12、基于步骤S11中定义的纤维长度L、纤维宽度W、纤维壁厚WT建立造纸纤维静态几何模型，造纸纤维静态几何模型基于真实造纸纤维形态沿纵向均匀分段，即定义所有分段长度为l_seq；造纸纤维静态几何模型的截面定义为同心圆环，但截面离散点仅用单层表示圆形轮廓，圆形轮廓直径为Φ＝W-2·WT；定义造纸纤维静态几何模型纵向共有N_seq个分段，则有N_seq+1个截面，记录第i，0≤i≤N_seq个截面的中心点位置为P_(i)，截面法向量为

通过上述定义可得，造纸纤维静态几何模型中第i个分段由序号为i-1和i的截面定义，P_(i)相对于P_(i-1)的位置由

和分段长度l_seq计算得到，如式(1)所示：

定义造纸纤维静态几何模型中的每一个截面被离散为N_dis个点，可定义三维空间中截面离散点的坐标为：p_(i)(j)＝(x_(i)(j),y_(i)(j),z_(i)(j))，0≤i≤N_seq，1≤j≤N_dis，表示第i个横截面上的第j个离散点，x_(i)(j),y_(i)(j)和z_(i)(j)分别表示离散点p_(i)(j)在空间坐标系中对应x、y和z轴的坐标值；其中，序号为0的截面上的离散点可取圆周上的任意位置为起点，等距离定义截面其它离散点位置坐标；后续截面上的点可由前一个截面上的离散点根据截面中心点坐标位置进行平移，再基于截面法向量的方向进行旋转计算得到；

S13、使用改进质点-弹簧模型建立造纸纤维特性模型：

所述改进质点-弹簧模型包括质点模型和弹簧模型，其中，所述质点模型用于模拟造纸纤维特性模型的轮廓和质量，其坐标与步骤S12中建立的造纸纤维静态几何模型中的截面离散点坐标一致，即造纸纤维特性模型中第i个截面上的第j个质点模型坐标mp_(i)(j)＝p_(i)(j)，弹簧模型分为纵向结构弹簧模型ST1、截面结构弹簧模型ST2、剪切弹簧模型SS和中心弹簧模型CS，其中，ST1用于连接造纸纤维特性模型前后截面中同一序号的质点模型，即连接mp_(i)(j)与mp_(i+1)(j)和mp_(i)(j)与mp_(i-1)(j)，其中mp_(i+1)(j)和mp_(i-1)(j)分别代表造纸纤维特性模型中第i+1个截面和第i-1个截面上的第j个质点模型；ST2用于连接造纸纤维特性模型的同一截面中相邻质点模型，即连接mp_(i)(j)与mp_(i)(j+1)和mp_(i)(j)与mp_(i)(j-1)，其中mp_(i)(j+1)和mp_(i)(j-1)分别代表造纸纤维特性模型中第i个截面上的第j+1个和第j-1个质点模型；SS用于连接造纸纤维特性模型中，相邻截面上相邻位置的质点模型，即连接mp_(i)(j)与mp_(i+1)(j+1)和mp_(i)(j)与mp_(i+1)(j-1)，其中，mp_(i+1)(j+1)和mp_(i+1)(j-1)分别代表造纸纤维特性模型中第i+1个截面上的第j+1个和第j-1个质点模型；CS用于连接造纸纤维特性模型中，同一截面上位置相对的一组质点模型，即连接mp_(i)(j)与

其中，

代表造纸纤维特性模型中第i个截面上的第

个质点模型；

S14、计算改进质点-弹簧模型参数：

首先计算改进质点-弹簧模型中质点模型的质量，假设造纸纤维特性模型的质量均匀分布，则通过式(2)计算质点模型的质量：

上式中，c和L分别表示用于建立造纸纤维特性模型的纤维粗度和纤维长度，其乘积表示造纸纤维特性模型的总质量；m_(i)(j)表示造纸纤维特性模型中第i个截面上第j个质点模型的质量；m_(0)(j)和

分别表示造纸纤维特性模型中，序号为0和N_seq截面上的质点模型的质量；

弹簧模型的弹性系数通过假设造纸纤维特性模型分段受到拉伸力或扭转力时的受力分析进行求解，定义造纸纤维特性模型分段的一个截面受到大小为F_s的拉伸力，另一个截面固定，根据步骤S11中对弹性模量E和纤维泊松比v的定义，以及基于材料力学中对弹性模量及泊松比的定义，计算得到造纸纤维特性模型分段在拉伸力F_s作用下的纵向形变量和截面形变量，从而得到形变后各弹簧模型的长度及偏转角度，分别包括在受到拉伸力后，连接目标质点模型的纵向结构弹簧长度l′_ST1、截面结构弹簧长度l′_ST2、剪切弹簧长度l′_SS和中心弹簧长度l′_CS，以及剪切弹簧在纵向受力分析中弹簧力与水平方向的夹角α和在截面受力分析中弹簧力与垂直方向的夹角β；受力截面中每个质点模型受到的拉伸力f_s可由截面受到拉伸力F_s对截面总质点模型数求平均得到，最终，列出截面质点模型的受力平衡方程如式(3)所示：

上式中，f_ST1、f_ST2、f_SS和f_CS分别代表纵向结构弹簧、截面结构弹簧、剪切弹簧以及中心弹簧的弹簧力；k_ST1、k_ST2、k_SS和k_CS分别代表纵向结构弹簧、截面结构弹簧、剪切弹簧以及中心弹簧的弹性系数；l_ST1、l_ST2、l_SS和l_CS分别代表纵向结构弹簧、截面结构弹簧、剪切弹簧以及中心弹簧的原始长度，等于该弹簧模型所连接的两个质点模型在受到拉伸力前的距离；

由材料力学理论可知，造纸纤维特性模型截面所受到扭转力会转换为质点模型受到切向力，定义每个质点模型受到切向力大小为f_T，则基于步骤S11中对剪切模量G的定义，以及基于材料力学中对剪切模量的定义，计算得到每个质点模型在切向力f_T作用下的偏转角度γ，从而得到形变后各弹簧的长度及偏转角度，包括连接目标质点模型的纵向结构弹簧受到切向力后的长度l″_ST1，左、右两侧剪切弹簧受到切向力后的长度l″_SS和l″′_SS，以及连接目标质点模型的纵向结构弹簧和左、右两侧剪切弹簧受到切向力后的弹簧力在纵向受力分析时与水平方向的夹角θ″_ST1、θ″_SS和θ″′_SS；

在针对造纸纤维特性模型截面受到扭转力进行受力分析时，忽略伴随纵向长度形变的截面形变，忽略截面受力分析时，弹簧力在截面投影方向与受力正交分解后水平方向的较小夹角，得到受力平衡方程如式(4)所示：

f_T＝f″_SScos(θ″_SS)+f″′_SScos(θ″′_SS)+f″_ST1cos(θ″_ST1)

其中：

考虑到纵向结构弹簧ST1和截面结构弹簧ST2均属于结构弹簧，具有相同性质，即相同弹性模量，根据材料力学中弹性模量的定义和胡克定律，得到两种结构弹簧间弹性系数的关系如式(5)所示：

上式中，E_f为两种结构弹簧的弹性模量，F和S分别表示纵向结构弹簧及截面结构弹簧所受到的拉伸力和受力面积，Δl_ST1和Δl_ST2表示纵向结构弹簧和截面结构弹簧受到拉伸力F后的形变量，定义两种结构弹簧的受力大小F及受力面积S一致；

将式(3)-(5)联立，求解纵向结构弹簧弹性系数k_ST1、截面结构弹簧弹性系数k_ST2、剪切弹簧弹性系数k_SS和中心弹簧弹性系数k_CS；

S2、生成纸页结构模型，过程如下：

S21、定义纤维形状因子及计算造纸纤维特性模型中相邻截面法向量偏转角度：

假设真实造纸纤维弯曲后形成一个圆弧，其对应弧线段和弦线段分别为AOB和AO′B，其中，A和B为弧线段及弦线段的两个端点，O和O′分别为弧线段和弦线段的中点，O同时为弧线段的最高点，则纤维形状因子C_f的定义如式(6)所示：

通过纤维形状因子的定义，基于简单几何方法计算得到∠OAB的角度

进一步与造纸纤维特性模型对应，在一个分段数为N_seq的造纸纤维特性模型中，序号为0和

的截面法向量之间的偏转角度等于

定义造纸纤维特性模型中所有相邻截面法向量偏转角度大小一致，则通过式(7)计算得到两个相邻截面法向量的偏转角度θ_N；

S22、定义纸页结构模型的尺寸及纸页克重：

定义纸页结构模型参数，包括纸页尺寸：长l_p、宽w_p和厚度t_p，以及纸页克重BW；

S23、模拟生成纸页结构模型：

采用随机位置、随机朝向的方法逐一生成并放置造纸纤维特性模型，并根据定义的纸页结构模型参数l_p、w_p、t_p和BW，以及步骤S11中定义的纤维长度L，调整纸页结构模型中造纸纤维特性模型的长度及数量；

S3、采用分段OBB包围盒检测纸页结构模型中造纸纤维特性模型之间的接触区域，过程如下：

S31、使用二叉树结构分解造纸纤维特性模型：

将造纸纤维特性模型首先沿纵向分段不断平分为两个部分，再按截面离散点不断平分为两个部分，从而得到多个造纸纤维特性模型局部，将完整造纸纤维特性模型作为二叉树结构的根节点；被平分前的造纸纤维特性模型局部作为被平分后造纸纤维特性模型局部的父节点，从而建立二叉树结构；

S32、计算二叉树结构中每个造纸纤维特性模型局部的OBB包围盒：

针对步骤S31中基于二叉树结构分解后的造纸纤维特性模型，对二叉树结构中每个节点所代表的造纸纤维特性模型局部建立OBB包围盒，即得到分段OBB包围盒结构；

S33、计算纸页结构模型中造纸纤维特性模型之间的接触区域：

通过步骤S32得到二叉树结构中不同造纸纤维特性模型局部OBB包围盒后，采用沿二叉树结构从上到下的顺序，依次比较两根造纸纤维特性模型局部OBB包围盒之间是否发生重叠现象，从而判断两根造纸纤维特性模型是否发生接触并计算发生接触的区域；

S4、采用结合弹簧模型模拟纤维间结合力：

在通过步骤S3计算得到造纸纤维特性模型的接触区域后，将相互接触的造纸纤维特性模型，在接触区域使用结合弹簧模型BS连接，BS的连接端点分别为两根造纸纤维特性模型接触区域包含的质点模型，定义两个相互接触的造纸纤维特性模型，其接触区域面积大小分别为S_A和S_B，接触区域表面法向量分别为

和

则通过判断结合弹簧模型BS的伸长方向与接触区域表面法向量

和

方向之间的平均夹角大小，判断结合弹簧模型BS当前受到的力为拉伸力或剪切力，若平均夹角小于45°，则认为BS受到拉伸力，否则BS受到剪切力；

定义BS受到拉伸力和剪切力时的最大受力强度分别为T_N和T_S，BS的最大受力伸长长度为D_bond，则在不同受力情况下，结合弹簧模型的弹性系数如式(8)所示：

上式中，k_N和k_S分别代表结合弹簧模型BS在受到拉伸力和剪切力时的弹性系数；

当结合弹簧模型长度超过D_bond时，标记该结合弹簧模型断裂，后续模拟将不再计算该结合弹簧模型的形变量及对应弹簧力；

S5、使用恒速拉伸法模拟纸页结构模型的抗张强度：

固定纸页结构模型中最左侧一段距离l_left内的质点模型，给纸页结构模型中最右侧一段距离l_right内的质点模型添加恒速向右的速度v_tensile，动态模拟纸页结构模型被恒速拉伸的形变过程，计算与纸页结构模型中最右侧一段距离l_right内的质点模型相连的所有弹簧模型，包括：纵向结构弹簧模型ST1、截面结构弹簧模型ST2、剪切弹簧模型SS、中心弹簧模型CS，以及结合弹簧模型BS产生的总弹簧力，作为纸页结构模型形变过程的抗张力σ；

记录纸页结构模型在形变过程中的抗张力σ，以及纸页结构模型中最右侧一段距离l_right内质点模型的移动距离，即伸长量δ，绘制出纸页结构模型的抗张力-伸长量关系曲线，再通过式(9)换算为纸页结构模型的抗张强度；

上式中，T_p为纸页结构模型的抗张强度；max(σ)、δ和w_p分别为纸页结构模型的最大抗张力、最大抗张力所对应的伸长量以及纸页结构模型的宽度。

2.根据权利要求1所述的一种纸页结构建模及抗张强度模拟的方法，其特征在于，所述的步骤S23中，调整纸页结构模型中造纸纤维特性模型长度及数量的过程如下：

S231、采用随机位置、随机朝向的方法逐一生成并放置造纸纤维特性模型，并在每放置一根造纸纤维特性模型后，通过删除超出纸页结构模型的尺寸长l_p、宽w_p和厚度t_p范围的造纸纤维特性模型纵向分段，以满足纸页结构模型的尺寸要求；

S232、在步骤S231不断添加造纸纤维特性模型的过程中，通过计算当前纸页结构模型的纸页克重BW，判断是否已达到步骤S22中定义值的要求，即与BW的定义值相比，偏差不超过5％；如果未达到要求，则通过步骤S231继续添加造纸纤维特性模型；否则，停止添加造纸纤维特性模型；

S233、通过计算当前纸页结构模型中所有造纸纤维特性模型的平均长度，以判断是否满足步骤S11中对纤维长度L的定义要求，即与L的定义值相比，偏差不超过5％，若纸页结构模型中造纸纤维特性模型的平均长度不满足定义值L的要求，通过逐个随机添加较长或较短的造纸纤维特性模型的方法，提高或降低造纸纤维特性模型的平均长度，以满足纤维长度定义值L的要求，其中，造纸纤维特性模型长度的随机取值范围为定义值L的0.5-1.5倍，在上述添加造纸纤维特性模型的过程中，如果纸页结构模型的纸页克重超出了定义值BW的要求，则需要逐个随机删除步骤S231添加的造纸纤维特性模型，直到满足纸页克重BW的定义值要求，通过反复调节，使最终生成的纸页结构模型同时满足定义值BW和L的要求。

3.根据权利要求2所述的一种纸页结构建模及抗张强度模拟的方法，其特征在于，所述的步骤S31中，建立二叉树结构的过程如下：

S311、将完整的造纸纤维特性模型作为二叉树结构根节点；

S312、将造纸纤维特性模型沿纵向分段平均分为两个部分，并将这两个部分作为二叉树结构根节点的左、右子节点，插入到二叉树结构中；

S313、针对二叉树结构中每个新添加的子节点，将其所代表的造纸纤维特性模型局部再次沿纵向分段平分为两部分，添加为对应子节点的左、右子节点，重复该过程直到当前子节点只表示造纸纤维特性模型的一个分段，即无法再沿纵向分段平分为两个部分为止；

S314、对于当前子节点只表示造纸纤维特性模型一个纵向分段的情况，需要基于造纸纤维特性模型截面离散点进一步平分成两个部分，添加到对应子节点下作为新的左、右子节点，重复该过程直到当前子节点所代表的造纸纤维特性模型截面也无法再继续平分，即每个截面只包含2个离散点为止。

4.根据权利要求3所述的一种纸页结构建模及抗张强度模拟的方法，其特征在于，所述的步骤S32中，对二叉树结构中每个节点所代表的造纸纤维特性模型局部建立OBB包围盒的步骤如下：

S321、依次获取二叉树结构中每个节点所代表的造纸纤维特性模型局部，统计其表面轮廓点坐标，即当前造纸纤维特性模型局部中所包含的截面离散点坐标，定义当前获取的所有表面轮廓点坐标按x、y和z轴分别表示为X＝(x₁,x₂,...,x_n)，Y＝(y₁,y₂,...,y_n)，Z＝(z₁,z₂,...,z_n)，其中n表示当前获取的所有表面轮廓点坐标个数，X，Y和Z分别表示所有表面轮廓点坐标在x，y和z轴上的坐标值集合，x_q、y_q、z_q分别表示第q个表面轮廓点在x、y和z轴的坐标值，1≤q≤n；

S322、计算步骤S321中所获取的所有表面轮廓点坐标相对于x、y和z轴的协方差矩阵，协方差矩阵定义如式(10)所示：

上式中，C表示协方差矩阵，cov表示协方差函数，cov(X,X)，cov(X,Y)，cov(X,Z)，cov(Y,X)，cov(Y,Y)，cov(Y,Z)，cov(Z,X)，cov(Z,Y)以及cov(Z,Z)分别表示括号内第一个参数所代表的表面轮廓点坐标在x、y或z轴上的坐标值集合，相对与括号内第二个参数所代表的x、y或z轴上坐标值集合之间的协方差；以cov(X,Y)为例，协方差函数计算方法如式(11)所示：

上式中，

和

分别代表集合X和Y中所有坐标值的平均值；

S323、求步骤S322中所得到协方差矩阵C的特征向量，过程如下：

定义协方差矩阵C的特征值和特征向量分别为λ和A，根据特征值和特征向量的定义，通过式(12)求解特征值λ的3个互不相同的解，记为λ₁、λ₂和λ₃；

上式中，E为单位矩阵；

将λ₁、λ₂和λ₃分别带入特征向量定义式：CA＝λA，求到三个不同的特征向量，记为A₁、A₂及A₃；

S324、使用步骤S323中计算得到的特征向量A₁、A₂及A₃，作为当前二叉树节点的OBB包围盒三维坐标轴，计算一个长方体包围盒模型，使该长方体包围盒模型的边与A₁、A₂或A₃平行，并且该长方体包围盒模型能刚好将步骤S323中获取的所有表面轮廓点坐标包围；将该长方体包围盒模型作为步骤S323中所选取的造纸纤维特性模型局部的OBB包围盒，并记录该长方体包围盒模型在A₁、A₂及A₃方向上的长度参数以表征该OBB包围盒的形状特点。

5.根据权利要求4所述的一种纸页结构建模及抗张强度模拟的方法，其特征在于，所述的步骤S33中，判断两根造纸纤维特性模型是否发生接触并计算发生接触区域的步骤如下：

S331、逐一遍历选取纸页结构模型中两个不同的造纸纤维特性模型；

S332、获取步骤S331中所选择的两根造纸纤维特性模型对应的分段OBB包围盒结构，并比较两个分段OBB包围盒结构的根节点OBB包围盒的位置关系，即判断根节点所表示的两个长方体包围盒模型之间是否包含重叠区域；如果不包含，则说明两根造纸纤维特性模型未发生接触，返回步骤S331，重新选择两根造纸纤维特性模型进行判断；如果包含重叠则继续进行后续步骤；

S333、依次比较两根造纸纤维特性模型的分段OBB包围盒结构中，根节点的左、右子节点所包含的OBB包围盒是否存在重叠关系，选择存在重叠关系的一对OBB包围盒所在节点，依次比较其左、右子节点所包含的OBB包围盒是否仍然存在重叠关系，重复该步骤，直到当前发生重叠关系的OBB包围盒所属二叉树结构节点已没有左、右子节点，则记录当前节点所代表的造纸纤维特性模型局部为造纸纤维特性模型之间的接触区域；若通过执行本步骤，已无法找到具有重叠关系的一对OBB包围盒，则返回步骤S331，重新选取两根造纸纤维特性模型进行步骤S332和步骤S333。

6.根据权利要求5所述的一种纸页结构建模及抗张强度模拟的方法，其特征在于，所述的动态模拟纸页结构模型被恒速拉伸的形变过程，具体如下：

将动态模拟过程的连续时间离散为多个时间点，每两个相邻的时间点之间的间隔为Δt，并假设在两个相邻时间点之间，纸页结构模型中所有质点模型的受力不变，通过分析每个质点模型在离散时刻t的受力情况，并运用牛顿运动定律计算每个质点模型在t+Δt时刻的位移量及新位置，为了提高计算精度，使用基于牛顿运动定律的龙格库塔法迭代计算每个质点模型在不同离散时间点的位置，直到纸页结构模型完全分离为两个部分。

Images