本发明涉及一种木结构用齿板连接节点受力性能的数值模拟方法，属于工程结构仿真技术领域。
背景技术：
齿板于上世纪50年代初期发明，是通过采用专用模具在由镀锌薄钢板构成的面板上冲压出齿的方法制造的。齿板有多种类型，其形状、尺寸也各异。如图1所示，一个典型的齿板1由面板11、在面板11上开的槽孔12以及槽孔两边竖立的齿13组成。规定平行于开槽方向为齿板主轴方向，垂直于开槽方向为齿板次轴方向。齿板主要用于制作轻型木桁架，将加工后的木构件按图纸摆放好，采用平板压机或滚轮压机将齿板压入木构件连接处形成齿板连接节点，如图2所示，其中2为木构件，1为齿板，齿板的齿13压入木构件2内，将木构件的顺纹方向称为木材纹理方向。多个齿板连接节点和木构件形成轻型木桁架，也称为齿板桁架。
在北美，轻型木桁架工业非常发达，超过60％的木结构住宅的屋顶采用轻型木桁架(canadawood,internationalbuildingseriesno.2.woodtrussesstrength,economy,versatility)。木结构房屋的安全性与轻型木桁架的受力性能密切相关，而轻型木桁架的受力性能一般是由齿板连接节点的受力性能决定。齿板连接节点的工作状态比较复杂，其受力性能随着荷载与木材纹理及荷载与齿板主轴方向的角度不同而不同。实际上齿板连接是一种广义的钉连接，其受力性能受齿板材料、尺寸及被连接木材的材性等多个因素的影响。板齿可看成是一种扁钉，其受力性能受木材纹理方向影响，且面板上开槽造成了面板在主轴方向和次轴方向的性能有差异。在荷载作用下，根据齿板连接节点的构造及受力方式不同，齿板的破坏包括齿拔出、齿板抗拉破坏及齿板剪切破坏等多种形式。
由于轻型木桁架的受力性能一般是由齿板连接节点的受力性能决定的，故对轻型木桁架的模拟的重点即是对齿板连接节点受力性能的模拟。下面列举2种现有的齿板连接节点受力性能模拟方法。
我国《轻型木桁架设计规范》推荐了一种模拟轻型木桁架受力性能的方法：用杆单元(element)模拟轻型木桁架杆件，尽量保证有限元模型中的杆单元与实际桁架的木构件轴线一致，各个杆件之间通过铰连接。图3给出了一个轻型木桁架的局部示意图，图中实线部分表示的是组成轻型木桁架的木构件2，由上弦杆、下弦杆及其连接处设置的一个楔形木垫块21组成，各个杆件之间通过齿板连接联系在一起。按我国《轻型木桁架设计规范》方法进行模拟时，图中的虚线即为模拟木构件的杆单元22，由于上、下弦杆轴线的交点在下弦杆外，且为了考虑垫块的影响，模拟时需在端节点(node)处设置一根竖向虚拟杆单元l1-l2，且在上弦杆端部设置一根虚拟杆单元l1-l3，这些虚拟杆件与模拟木构件的单元22及各个杆单元连接处的铰接节点23共同组成轻型木桁架有限元模型。这种方法的不足主要体现在两方面：1)虚拟杆单元的引入导致有限元模型与桁架实际构造不一致，虚拟杆件的截面、弹性模量等参数确定方法较为复杂，通用性差；2)采用铰接的简化方式不能模拟齿板连接节点的半刚性受力性能，准确度较低。
现介绍一种现有的通过专用有限元程序模拟齿板连接节点受力性能的方法，该方法分别采用平面单元模拟木构件和齿板面板、每个齿连接均用一个虚拟工字梁模拟。图4给出了使用该方法的建立的一个齿板连接节点的有限元模型示意图，其中3为模拟木构件2的平面单元，4为模拟齿板面板11的平面单元，5为模拟单个齿13与木构件2连接的单元。图5给出了该方法中单个齿连接模拟方法的局部放大图，模拟单个齿连接的单元5采用虚拟工字梁，该虚拟工字梁一端与木构件单元3的节点连接，另一端与齿面板单元4的节点连接。该方法通过引入工字梁以及调整工字梁截面的性质，可以模拟齿板连接节点的半刚性特性，也可以较为灵活的模拟不同形状齿板连接节点的受力性能，然而该方法也有如下不足之处：
1)不能模拟齿板连接节点受力性能受木材纹理方向及齿板主轴方向双重影响的特点
工字梁为双轴对称截面，随荷载作用方向与截面对称轴方向夹角的不同，其受力性能不同，具有2个互相垂直的主方向。工字梁截面的放置方向可作为一个参数，表征木纹方向或齿板主轴方向二者之一。然而实际齿板连接节点的受力性能受木纹方向或齿板主轴方向的双重影响，需要两个参数来描述齿板连接节点的构造特点，上述方法只具有一个参数，故该方法不能描述齿板连接节点受力性能受木材纹理方向及齿板主轴方向双重影响的特点，具有一定的局限性。
2)拟合的参数与齿板连接节点试验结果无显式关系，拓展性差
由于模型是平面模型，仅考虑虚拟工字梁在齿板平面内的平动自由度。虚拟工字梁的截面尺寸、弹性模量、流动和硬化参数等需要通过试验进行校准。由于这些参数与齿板连接节点的试验结果并无显式关系，连接构造改变时，已经拟合的参数不再适用，需要重新进行对应试验，并再次拟合，因而该方法的拓展性较差。
3)模型的建立、求解、后处理等全部采用自编程序，实施难度较大，操作较为复杂，可读性及通用性较差
齿板连接节点有限元模拟与一般有限元模型模拟的主要区别在于齿连接受力性能的模拟，木构件等常见构件已有较为成熟的模拟方法，可以采用通用有限元程序进行模拟，而上述方法对木构件、齿板面板和连接全部采用了自编程序，其工作量较大，且降低了自编程序的可读性。此外自编程序与通用有限元软件的数据难以互通，通用性较差。
4)求解计算量大
由于采用一个实体的工字型梁模拟一个齿的连接受力性能，且需要通过工字型梁的弹塑性变形来模拟齿板连接节点的半刚性性能，计算过程中每个增量步均需要计算若干个(数量等于模型中齿的数量)工字梁的弹塑性受力性能，其计算量非常大。
有限单元法是一种数值计算方法，其基本思想是把计算域划分为有限个互不重叠的单元(element)，每个单元由若干节点(node)组成，这些节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。
有限单元法在工程结构领域应用非常广泛，并形成了多种通用有限元软件。通用有限元软件一般由前处理模块、分析模块和后处理模块组成。如图6所示，前处理模块用于建立有限元模型，并生成描述有限元模型的数据文件。分析模块读取模型数据文件，并进行分析计算，将计算结果输出到结果数据文件。后处理模块读取计算结果数据文件，并对计算结果进行进一步地描述和解释。
模型数据文件包括有限元模型及求解所需的全部数据，通常包括：1)节点和单元等几何数据；2)材料定义；3)零件组合方式；4)初始应力、温度；5)边界条件；6)分析类型；7)荷载；8)输出要求；9)子程序信息等。一般有两种生成模型数据文件的方式，一种是用户在有限元软件的图形界面中进行点选、输入、绘图等操作，有限元软件将用户输入转化成为模型的数据文件；另一种是按照有限元软件规定的数据格式直接编写模型数据文件。
有限元软件的核心是其分析模块，其实质为采用数值方法求解关于基本未知量的代数方程组或常微分方程组的过程，目前有中心差分法，牛顿-拉普森法等多种方法，以下以中心差分法为例，结合图6简要介绍通过有限元的求解流程得到有限元模型受力性能的方法如下：
步骤1：读取前处理器建立的有限元模型数据文件，并初始化各变量，包括设置加载时间t＝0。
步骤2：当加载时间为t，时间增量为dt时，通过建立动力方程求解节点的加速度；该动力方程是描述节点质量、作用于节点的内力、作用于节点的外力以及节点加速度之间关系的方程，可表示为：
采用中心差分法，将节点加速度
其中，
步骤3：根据步骤2求解的节点位移建立并求解几何方程，得到单元的应变增量；该几何方程是描述单元应变和节点位移之间关系的方程，可表示为dε＝l(u(t+dt),u(t))，其中l为与单元类型有关的函数，dε为从加载时间t到加载时间t+dt的单元应变增量。
步骤4：根据步骤3求解的应变增量建立并求解物理方程，得到单元的应力；物理方程是描述材料应力和应变之间关系的方程，可表示为σ(t+dt)＝f(σ(t),dε)，其中f为材料的性质有关的函数，也叫本构关系函数，物理方程也可称为本构关系模型；σ(t+dt)和σ(t)分别为加载时间t+dt和加载时间t的单元应力。
步骤5：进行节点内力汇总：一个节点一般与多个单元连接，将与节点连接的各个单元分配到节点的内力进行汇总，得到加载时间为t+dt的节点内力i(t+dt)。
步骤6：计算过程中，在规定的时间点将节点位移、单元应力、应变等写入计算结果数据文件中。若加载时间t小于分析步时间t，则将加载时间t+dt设置为t，返回到步骤1进行计算；否则计算结束。
计算结束后，可以采用通用有限元软件的后处理器读取结果数据文件，进行应力云图绘制，数据导出等后处理操作。
通用有限元软件虽然功能强大，但是在遇到特殊的工程材料、结构时，仅靠通用有限元软件自带的常规材料模型、常规单元种类等并不能满足模拟需求。故通用有限元软件一般都提供了二次开发接口，以便使用者在原有功能基础上进行相关功能的定制和扩展。
技术实现要素：
本发明为了解决现有技术存在的问题，提供了一种木结构用齿板连接节点受力性能的数值模拟方法。本发明提供的方法能在模型中考虑木材的纹理方向、齿板的主轴方向，并在考虑荷载作用方向与二者的夹角的情况下确定齿板连接节点的受力性能。本模型通过建立每个板齿处的连接单元来模拟齿板连接节点的受力性能，可以考虑各种齿板形状、方向以及齿密度的齿板连接节点受力性能，克服了现有模拟方法准确性低、建模方法复杂、参数较难确定等不足。
本发明是通过以下技术方案实现的：
一种木结构用齿板连接节点受力性能的数值模拟方法，其特征在于，该方法具体包括以下步骤：
1)建立有限元模型，生成有限元模型数据文件
1.1)确定有限元模型的单元
建立齿板连接节点的有限元模型，该模型包括三种类型单元：模拟木构件的单元、模拟齿板面板的单元以及模拟齿板面板与木构件间连接的单个齿的自定义单元，木构件和齿板面板在每个齿处均通过一个自定义单元连接；其中，模拟木构件的单元和模拟齿板面板的单元均采用平面应力单元模拟；模拟单个齿的自定义单元为2个节点的自定义单元，一个节点位于齿板面板单元上，另一个节点位于木构件单元上；
确定单元类型后，根据齿板连接节点的实际构造，按有限元模型数据文件编写规则，输入木构件、齿板面板以及自定义单元的各个节点的编号及对应的几何坐标；对于木构件和齿板面板单元，通过按顺序依次指定各单元内各个节点编号定义一个平面应力或单元；对于自定义单元，通过指定各自定义单元2个节点的节点编号定义一个2节点单元；
1.2)确定有限元模型输入参数和输出变量并生成有限元模型数据文件；
1.2.1)确定模拟木构件间的接触性能、加载时间的输入参数和输出变量：
模拟木构件间的接触性能，在有限元模型中两木构件单元的接触边定义接触对，在模型数据文件中输入接触对之间的摩擦系数；
确定输入加载步时长t、施加荷载的类型及大小、各类单元的阻尼，以及结果保存的次数n；输出变量的种类，包括：节点位移、节点速度、节点加速度、单元应力、单元应变；
1.2.2)确定木构件单元和齿板面板单元的材性参数，包括木材的顺纹方向弹性模量和横纹方向弹性模量、齿板主轴方向的弹性模量、次轴方向的弹性模量；
1.2.3)定义自定义单元的材性参数，包括用于描述4种标准工况下齿连接荷载-位移曲线的参数；齿板主轴方向与全局坐标系x轴的夹角βs；木构件的纹理方向与全局坐标系x轴的夹角βw；单个齿的自定义单元的质量m；其中，
用于描述4种标准工况下齿连接荷载-位移曲线的参数按下述方法确定并输入：
所述4种标准工况包括：aa-荷载平行于木材纹理及齿板主轴，ae-荷载平行于木材纹理但垂直于齿板主轴，ea-荷载垂直于木材纹理但平行于齿板主轴，ee-荷载垂直于木材纹理及齿板主轴；通过板齿极限承载力试验分别获得4种标准工况下齿板连接节点的荷载-位移曲线，并将所得曲线中的荷载除以试验时试件的有效齿连接个数，得到考虑群体效应的单个齿连接的荷载-位移曲线；选用公式(1)所示的函数分别拟合4种标准工况的荷载-位移曲线：
式中，f为单个齿连接承担的荷载，δ为齿板连接节点的整体位移，f0、k1、k2分别为3个待拟合的齿连接荷载-位移曲线参数；
通过公式(1)得到4种标准工况的齿连接荷载-位移曲线参数分别为：aa组(f0_aa、k1_aa、k2_aa)，ae组(f0_ae、k1_ae、k2_ae)，ea组(f0_ea、k1_ea、k2_ea)，ee组(f0_ee、k1_ee、k2_ee)；每组3个共计12个参数；每个自定义单元质量m按一个齿连接的实际质量作为输入参数；
将上述有限元模型的输入参数、材性参数及输出变量生成有限元模型数据文件；
2)获取有限元模型数据文件的数据，求解有限元模型并输出结果到结果数据文件中，该输出结果作为木结构用齿板连接节点受力性能的数值模拟；
2.1)读取根据步骤1)确定的有限元模型数据文件，计算有限元模型中各类单元的质量及稳定增量步长，初始化各变量，其中设置加载时间t＝0：
所述木构件和齿板面板单元为常规的平面应力单元，木构件和齿板面板单元的质量通过单元的几何尺寸和材料密度计算得到，木构件和齿板面板单元的稳定增量步长由单元的特征尺寸与材料膨胀波速的比值确定；
计算该自定义单元稳定增量步长：自定义单元的初始刚度取k1aa，自定义单元的稳定增量步长按
取有限元模型中全部单元的稳定增量步长的最小值作为有限元模型的稳定增量步长dt；进行各变量初始化；
2.2)在加载时间为t，增量步为dt时对有限元模型进行计算，完成木结构用齿板连接节点受力性能的数值模拟；具体包括：
2.2.1)建立并求解动力方程，确定加载时间为t时的有限元模型所有节点的加速度，通过中心差分法计算得到t+dt/2时刻所有节点的速度和t+dt时刻所有节点的位移；
2.2.2)根据步骤2.2.1)得到t+dt/2时刻各个节点的速度和t+dt时刻各个节点的位移建立并求解几何方程，得到加载时间t到t+dt时所有单元的应变增量dε；
2.2.3)根据步骤2.2.2)得到的应变增量对木构件和齿板面板单元建立并求解物理方程，计算加载时间为t+dt时的单元应力；根据计算得到木构件和齿板面板单元的应力计算分配到自定义单元2个节点的节点内力；
2.2.4)对各个单元的节点的内力进行汇总，即将与节点连接的各个单元分配到该节点的内力相加，得到t+dt时刻的节点内力；
2.2.5)若加载时间t为t/n的整数倍，则将此时得到的有限元模型的节点位移、节点速度、节点加速度、单元应力、单元应变输入到结果数据文件中，执行步骤2.2.6)；否则，直接执行步骤2.2.6)；
2.2.6)如果t<t，将加载时间t设置为t＝t+dt，返回步骤2.2.1)步进行下一增量步计算；如果t＝t，保存计算结果，该输出结果作为木结构用齿板连接节点受力性能的数值模拟。
与目前大部分模拟齿板连接节点的方法相比，本发明具有如下几个优点：
(1)齿板连接节点有限元模型所需输入参数较易确定，选用我国《轻型木桁架技术规范》规定的齿板连接节点标准试验(与加拿大tpic规范规定相同)结果、齿板用钢材的性能试验结果即可确定齿板连接节点模拟所需的所有参数，且试验工况与模拟工况不需要一致。
(2)本发明可以预测不同齿板形状、木材纹理方向、齿板主轴方向时的齿板连接节点的受力性能。
(3)本发明可以自动判断齿板荷载方向与木材纹理方向和齿板主轴方向的关系，并确定该方向齿板连接节点的荷载-位移曲线参数，无需对齿板的荷载作用方向进行预判。
(4)本发明将构件受力性能及构件接触等的计算交由通用有限元程序处理，通过二次开发模拟齿板连接节点的受力性能。具有更好的通用性，适用性更强。
综上所述，本发明可更为准确地模拟不同齿板形状、各种复杂构造齿板连接节点的受力性能，且与通用有限元结合，具有操作简便，适用性强等优势。
附图说明
图1为典型齿板结构示意图。
图2为典型齿板连接节点示意图。
图3为轻型木桁架虚拟杆件设置示意图。
图4为现有齿板连接节点有限元模型示意图。
图5为图4的局部放大示意图。
图6为通用有限元求解流程框图。
图7为本发明的方法流程框图。
图8齿板连接节点有限元模型示意图。
图9板齿承载力试验标准工况。
具体实施方式
结合附图对本发明做进一步详细说明：本实施例在本发明技术方案前提下实施，给出了详细的实施方式，但本发明的保护范围不限于下述实施例。
本发明提出的一种木结构用齿板连接节点受力性能的数值模拟方法，具体流程如图7所示，包括如下步骤：
1)建立有限元模型，生成有限元模型数据文件
1.1)确定有限元模型的单元
建立齿板连接节点的有限元模型，该模型包括三种类型单元：模拟木构件的单元、模拟齿板面板的单元以及模拟齿板面板与木构件间连接的单个齿的自定义单元，木构件和齿板面板在每个齿处均通过一个自定义单元连接；其中，模拟木构件的单元和模拟齿板面板的单元均采用平面应力单元模拟；模拟单个齿的自定义单元为2个节点的自定义单元，一个节点位于齿板面板单元上，另一个节点位于木构件单元上；
确定单元类型后，根据齿板连接节点的实际构造，按有限元模型数据文件编写规则，输入木构件、齿板面板以及自定义单元的各个节点的编号及对应的几何坐标；对于木构件和齿板面板单元，通过按顺序依次指定各单元内各个节点编号定义一个平面应力或单元；对于自定义单元，通过指定各自定义单元2个节点的节点编号定义一个2节点单元；
1.2)确定有限元模型输入参数和输出变量并生成有限元模型数据文件；
1.2.1)确定模拟木构件间的接触性能、加载时间的输入参数和输出变量：
模拟木构件间的接触性能，在有限元模型中两木构件单元的接触边定义接触对，在模型数据文件中输入接触对之间的摩擦系数；
确定输入加载步时长t、施加荷载的类型及大小、各类单元的阻尼，以及结果保存的次数n；输出变量的种类，包括：节点位移、节点速度、节点加速度、单元应力、单元应变；
1.2.2)确定木构件单元和齿板面板单元的材性参数，包括木材的顺纹方向弹性模量和横纹方向弹性模量、齿板主轴方向的弹性模量、次轴方向的弹性模量；
1.2.3)定义自定义单元的材性参数，包括用于描述4种标准工况下齿连接荷载-位移曲线的参数；齿板主轴方向与全局坐标系x轴的夹角βs；木构件的纹理方向与全局坐标系x轴的夹角βw；单个齿的自定义单元的质量m；其中，
用于描述4种标准工况下齿连接荷载-位移曲线的参数按下述方法确定并输入：
所述4种标准工况包括：aa-荷载平行于木材纹理及齿板主轴，ae-荷载平行于木材纹理但垂直于齿板主轴，ea-荷载垂直于木材纹理但平行于齿板主轴，ee-荷载垂直于木材纹理及齿板主轴；通过板齿极限承载力试验分别获得4种标准工况下齿板连接节点的荷载-位移曲线，并将所得曲线中的荷载除以试验时试件的有效齿连接个数，得到考虑群体效应的单个齿连接的荷载-位移曲线；选用公式(1)所示的函数分别拟合4种标准工况的荷载-位移曲线：
式中，f为单个齿连接承担的荷载，δ为齿板连接节点的整体位移，f0、k1、k2分别为3个待拟合的齿连接荷载-位移曲线参数；
通过公式(1)得到4种标准工况的齿连接荷载-位移曲线参数分别为：aa组(f0_aa、k1_aa、k2_aa)，ae组(f0_ae、k1_ae、k2_ae)，ea组(f0_ea、k1_ea、k2_ea)，ee组(f0_ee、k1_ee、k2_ee)；每组3个共计12个参数；每个自定义单元质量m按一个齿连接的实际质量作为输入参数；
将上述有限元模型的输入参数、材性参数及输出变量生成有限元模型数据文件；
2)获取有限元模型数据文件的数据，求解有限元模型并输出结果到结果数据文件中，该输出结果作为木结构用齿板连接节点受力性能的数值模拟；
2.1)读取根据步骤1)确定的有限元模型数据文件，计算有限元模型中各类单元的质量及稳定增量步长，初始化各变量，其中设置加载时间t＝0：
所述木构件和齿板面板单元为常规的平面应力单元，木构件和齿板面板单元的质量通过单元的几何尺寸和材料密度计算得到，木构件和齿板面板单元的稳定增量步长由单元的特征尺寸与材料膨胀波速的比值确定；
计算该自定义单元稳定增量步长：自定义单元的初始刚度取k1aa，自定义单元的稳定增量步长按
取有限元模型中全部单元的稳定增量步长的最小值作为有限元模型的稳定增量步长dt；进行各变量初始化；
2.2)在加载时间为t，增量步为dt时对有限元模型进行计算，完成木结构用齿板连接节点受力性能的数值模拟；具体包括：
2.2.1)建立并求解动力方程，确定加载时间为t时的有限元模型所有节点的加速度，通过中心差分法计算得到t+dt/2时刻所有节点的速度和t+dt时刻所有节点的位移；
2.2.2)根据步骤2.2.1)得到t+dt/2时刻各个节点的速度和t+dt时刻各个节点的位移建立并求解几何方程，得到加载时间t到t+dt时所有单元的应变增量dε；
2.2.3)根据步骤2.2.2)得到的应变增量对木构件和齿板面板单元建立并求解物理方程，计算加载时间为t+dt时的单元应力；根据计算得到木构件和齿板面板单元的应力计算分配到自定义单元2个节点的节点内力；
2.2.4)对各个单元的节点的内力进行汇总，即将与节点连接的各个单元分配到该节点的内力相加，得到t+dt时刻的节点内力；
2.2.5)若加载时间t为t/n的整数倍，则将此时得到的有限元模型的节点位移、节点速度、节点加速度、单元应力、单元应变输入到结果数据文件中，执行步骤2.2.6)；否则，直接执行步骤2.2.6)；
2.2.6)如果t<t，将加载时间t设置为t＝t+dt，返回步骤2.2.1)步进行下一增量步计算；如果t＝t，保存计算结果，该输出结果作为木结构用齿板连接节点受力性能的数值模拟。
实施例1
本实施例为通用有限元软件abaqus及其提供的二次开发接口的一种木结构用齿板连接节点受力性能的数值模拟方法的实现方式。具体流程如图7所示，具体包括如下步骤：
1)建立有限元模型，生成有限元模型数据文件
结合实际工程中轻型木桁架或齿板连接节点的构造，建立轻型木桁架或齿板连接节点的有限元模型，生成模型数据文件；本实施例以图2所示齿板连接节点为对象建立其有限元模型，该连接节点采用的齿板为如图1所示的典型齿板。具体包括以下步骤：
1.1)确定有限元模型的单元：
本发明采用通用有限元软件(如abaqus、ansys、adina等)结合其二次开发功能建立如图2所示齿板连接节点的有限元模型(本实施例采用的通用有限元软件为abaqus)，如图8所示，该模型主要由三种类型单元组成：模拟木构件的单元(图中两个较大矩形中的小方块表示的每个单元)、模拟齿板面板的单元(图中在两个矩形相邻处的矩形中的小方块表示的每个单元)以及模拟齿板面板与木构件间连接的单个齿的自定义单元(图中用“×”号表示的每个自定义单元)，模拟木构件的有限元模型单元的相邻边分别用gh和ef表示；木构件和齿板面板在每个齿处均通过一个自定义单元连接。其中，模拟木构件的单元和模拟齿板面板的单元均采用平面应力单元模拟(本实施例采用abaqus自带的4节点缩减积分平面应力单元cps4r单元模拟)；模拟单个齿的自定义单元为2个节点的自定义单元，一个节点位于齿板面板单元上，另一个节点位于木构件单元上，该自定义单元采用通用有限元软件相应的二次开发编程语言(本实施例选用abaqus对应二次开发语言fortran语言)编写。
确定单元类型后，根据齿板连接节点的实际构造，按通用有限元软件的有限元模型数据文件编写规则，输入木构件、齿板面板和自定义单元的各个节点的编号及对应的几何坐标，对于木构件和齿板面板单元，按逆时针的方向依次指定四个节点编号即可定义一个cps4r单元，对于每个自定义单元，依次指定该单元的2个节点的节点编号即可；按上述方法输入所有木构件、齿板面板和自定义单元的几何信息。
1.2)确定有限元模型输入参数和输出变量并生成有限元模型数据文件
1.2.1)确定模拟木构件间的接触性能、加载时间等输入参数和输出变量：
按通用有限元软件的方法模拟木构件间的接触性能，本实施例在有限元模型两木构件单元的gh边和ef边定义接触对(contactpair)，在模型数据文件中输入接触对之间的摩擦系数；
按通用有限元软件的要求确定输入加载步时长t(t的取值应保证模型所有单元的惯性力合力小于作用于整个模型的外荷载的5％，由试算确定)、施加荷载的类型及大小、各类单元的阻尼等，以及结果保存的次数n(n>2,本实施例采用等时间间隔输出法，即当加载时间t为t/n的整数倍时输出计算结果到结果数据文件)；输出变量的种类(节点位移、节点速度、节点加速度、单元应力、单元应变)；
1.2.2)确定木构件单元和齿板面板单元的材性参数，包括木材的顺纹方向弹性模量和横纹方向弹性模量、齿板主轴方向的弹性模量、次轴方向的弹性模量：
木构件和齿板面板单元均采用通用有限元软件自带单元cps4r模拟，按通用有限元软件的要求输入该类单元所需要的对应的材性参数；本实施例木材采用正交各向异性线弹性本构关系，其中木材的顺纹弹性模量和横纹弹性模量按照相关标准如《木材抗弯弹性模量测定方法》(gb/t1936.2)和《木材横纹抗压弹性模量测试方法》(gb/t1943)对木构件采用的木材进行测试获得。本实施例齿板钢材采用正交各向异性线弹性本构关系，输入的钢材弹性模量按齿板主轴方向和次轴方向输入，按下述方法确定：按照《金属材料弹性模量和泊松比试验方法》(gb/t22315)测试制造齿板的钢材的弹性模量es。齿板主轴方向的弹性模量取值为
1.2.3)定义自定义单元的材性参数：
本实施例的单个齿的自定义单元为通过对通用有限元软件二次开发编写的单元，用于有限元模型输入的单个齿的自定义单元材性参数包括：用于描述4种标准工况下齿连接荷载-位移曲线的参数；齿板主轴方向与全局坐标系x轴的夹角βs；木构件的纹理方向与全局坐标系x轴的夹角βw；单个齿的自定义单元的质量m；其中，
用于描述4种标准工况下齿连接荷载-位移曲线的参数按下述方法确定并输入：
按照我国《轻型木桁架技术规范》，板齿极限承载力试验的标准工况包括4种，图9中1代表齿板，2代表木构件，图中表示了试验时荷载的施加方向、木构件纹理方向和齿板主轴方向。如图9所示：图(a)中aa-荷载平行于木材纹理及齿板主轴；(b)中ae-荷载平行于木材纹理但垂直于齿板主轴；(c)中ea-荷载垂直于木材纹理但平行于齿板主轴；(d)中ee-荷载垂直于木材纹理及齿板主轴。通过板齿极限承载力试验分别获得4种标准工况下齿板连接节点的荷载-位移曲线，并将这些曲线中的荷载除以试验时试件的有效齿连接个数，得到考虑群体效应的单个齿连接的荷载-位移曲线。选用合适的函数分别拟合四种标准工况的荷载-位移曲线，得到描述齿连接荷载-位移曲线的参数。本实施例采用如式(1)所示函数进行拟合：
式中，f为单个齿连接承担的荷载，δ为齿板连接节点的整体位移，f0、k1、k2分别为3个待拟合的齿连接荷载-位移曲线参数。分别对上述四种标准工况的荷载-位移曲线进行拟合，得到4种标准工况的齿连接荷载-位移曲线参数分别为：aa组(f0_aa、k1_aa、k2_aa)；ae组(f0_ae、k1_ae、k2_ae)；ea组(f0_ea、k1_ea、k2_ea)；ee组(f0_ee、k1_ee、k2_ee)，每组3个共计12个参数；
每个自定义单元质量按一个齿连接的实际质量作为输入参数；
将上述有限元模型的输入参数、材性参数及输出变量生成有限元模型数据文件；
2)获取有限元模型数据文件的数据，求解有限元模型并输出结果到结果数据文件中
本实施例采用通用有限元软件求解器(本实施例为abaqus的显式求解器explicit)求解步骤1)确定的有限元模型，并输出计算结果到结果数据文件中，具体求解步骤如下：
2.1)读取根据步骤1)确定的有限元模型数据文件，计算有限元模型中各类单元的质量及稳定增量步长，初始化各变量，其中设置加载时间t＝0：
对于本实施例中木构件和齿板面板单元采用的通用有限元软件自带单元，木构件和齿板面板单元的质量通过单元的几何尺寸和材料密度计算得到，木构件和齿板面板单元的稳定增量步长由单元的特征尺寸与材料膨胀波速的比值确定；常规单元的单元质量和稳定增量步长计算方法与单元类型有关，可直接通过通用有限元软件求解器计算木构件和齿板面板单元的质量和稳定增量步长；
对于本实施例中模拟单个齿连的自定义单元，该自定义单元质量取输入值m，估算该自定义单元稳定增量步长时，自定义单元的初始刚度取k1aa，自定义单元的稳定增量步长按
根据通用有限元软件求解器获得有限元模型中全部自带单元和自定义单元的稳定增量步长后，取有限元模型中全部单元的稳定增量步长的最小值作为有限元模型的稳定增量步长dt；进行各变量初始化；
2.2)在加载时间为t，增量步为dt时对有限元模型进行计算，完成木结构用齿板连接节点受力性能的数值模拟；具体包括：
2.2.1)通过通用有限元软件求解器建立并求解动力方程，确定加载时间为t时的有限元模型所有节点的加速度，通过中心差分法计算得到t+dt/2时刻各个节点的速度和t+dt时刻各个节点的位移；
2.2.2)通过通用有限元软件求解器根据步骤2.2.1)求解的节点位移建立并求解几何方程，得到加载时间t到t+dt时所有单元的应变增量dε；
2.2.3)通过通用有限元软件求解器根据步骤2.2.2)求解的应变增量对木构件和齿板面板单元建立并求解物理方程，计算加载时间为t+dt时木构件和齿板面板单元的应力；根据计算得到的单元应力计算分配到自定义单元2个节点的节点内力；具体流程如下：
2.2.3.1)读取有限元模型数据文件中输入的自定义单元材性参数、加载时间为t时的内部变量值以及从通用有限元软件求解器传输的所有自定义单元几何信息等3类参数，包括：
从有限元模型数据文件中读取自定义单元的材性参数，包括：4种标准工况的齿连接荷载-位移曲线参数，本实施例为f0_aa、k1_aa、k2_aa、f0_ae、k1_ae、k2_ae、f0_ea、k1_ea、k2_ea、f0_ee、k1_ee、k2_ee等12个参数；齿板主轴方向与全局坐标系x轴的夹角βs；木材纹理方向与全局坐标系x轴的夹角βw；自定义单元质量m等；
从自定义单元计算的内部变量中读取加载时间为t时的各个内部变量的值，包括：自定义单元的单元内力
从通用有限元软件求解器获取所有自定义单元的几何信息，包括：自定义单元在全局坐标系下的节点坐标x(j)、y(j)，节点位移增量dx(j)、dy(j)等信息，其中j取1或2，分别代表自定义单元的第1或第2个节点；
2.2.3.2)计算齿连接荷载作用方向
按式(2)-(5)计算自定义单元2个节点在全局坐标系下x轴方向、y轴方向的坐标差分量dx、dy以及位移增量差分量ddx、ddy：
dx＝x(2)-x(1)(2)
ddx＝dx(2)-dx(1)(3)
dy＝y(2)-y(1)(4)
ddy＝dy(2)-dy(1)(5)
荷载作用方向取为与自定义单元位移方向相同，荷载作用方向与x轴夹角按
2.2.3.3)确定齿连接荷载-位移曲线参数：
将计算得到的荷载作用方向β(t+dt)与加载时间为t时的荷载作用方向β(t)进行比较，如没有变化，则采用原荷载-位移曲线参数
其中：
其中，
2.2.3.4)建立齿连接荷载-位移曲线，并计算自定义单元内力：
按照上述齿连接荷载-位移曲线参数
其中，
将局部坐标系x’轴设置为与荷载作用方向平行，计算局部坐标系下x’轴方向的单元变形δ(t)及单元变形增量dδ，如式(10)所示：
加载时间为t+dt时，局部坐标系x’轴方向的单元总变形为δ(t+dt)＝δ(t)+dδ；
如果δ(t)和dδ的符号相同，则处于加载状态，则将δ(t+dt)和dδ代入荷载-位移曲线函数(9)得到加载时间为t+dt时的单元内力
然后按
将局部坐标系下的单元内力
由自定义单元内力得到自定义单元中节点1和节点2的节点内力分别为
将加载时间为t+dt时的自定义单元的单元内力
2.2.4)由于有限元模型中一个节点往往同时属于多个单元，如本实施例模拟木构件的每个cps4r单元均包含4个节点(即图8中两个较大矩形中每个小方块的四个角点)，构件内部的一个节点同时属于与之连接的4个单元。对各个单元进行计算后，对每个节点的内力进行汇总，即将与节点连接的各个单元分配到该节点的内力相加，得到t+dt时刻的节点内力。
2.2.5)若加载时间t为t/n的整数倍，则将此时得到的有限元模型的节点位移、节点速度、节点加速度、单元应力、单元应变输入到结果数据文件中，执行步骤2.2.6)；否则，直接执行步骤2.2.6)；
2.2.6)如果t<t，将加载时间t设置为t＝t+dt，返回步骤2.2.1)步进行下一增量步计算；如果t＝t，保存计算结果(节点位移、节点速度、节点加速度、单元应力、单元应变)，完成木结构用齿板连接节点受力性能的数值模拟。
实施例2
与实施例1计算相同，只是木材采用理想弹塑性模型。
以上所述，仅为本发明较佳的具体实施方式，这些具体实施方式都是基于本发明整体构思下的不同实现方式，而且本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。