[发明专利]一种适于中国古建筑木结构用材本构关系的数值模拟方法

权利要求书

1.一种适于中国古建筑木结构用材本构关系的数值模拟方法，其特征在于，包括以下步骤：

(1)弹性阶段：中国古建筑木结构用材弹性阶段用以下公式描述：

其中D_ijkl为用材弹性刚度矩阵，i，j，k，l＝1，2，3；σ₁₁，σ₂₂，σ₃₃为用材的主轴正应力，σ₁₂，σ₂₃，σ₁₃为用材的主面剪应力；ε₁₁，ε₂₂，ε₃₃为用材的主轴正应变，ε₁₂，ε₂₃，ε₁₃为用材的主面剪应变；在后续描述中，LR平面内的应力使用σ进行表达，应变使用ε进行表达；RT平面内的应力使用Δ进行表达，应变使用ξ进行表达；TL平面内的应力使用α进行表达，应变使用β进行表达；

用材弹性刚度矩阵D_ijkl中的参数的表达具体如公式(2)所示；

其中E₁，E₂，E₃为用材主轴弹性模量，G₁₂，G₂₃，G₁₃为用材主面剪切模量，△表示中间过渡计算量，ν₁₂，ν₂₃，ν₁₃，ν₂₁，ν₃₂，ν₃₁为用材主面泊松比；

(2)屈服方程：

A.LR平面屈服方程：

将SunChen模型的衍生模型与能够描述用材LR平面内受力响应的力学假设相结合，得到LR平面屈服方程：

其中a₁₁，a₂₂，b₁₁，b₂₂为LR平面内用材塑性特征参数，μ为LR平面内用材弹性特征参数，k₁是与用材应力有关的参数，a₁₁，a₂₂，b₁₁，b₂₂通过中国古建筑木结构修缮替换件用材材料偏轴受拉和受压应力-应变曲线获取，将应力-应变曲线转换为应力-塑性应变曲线，使用如下公式获取有效应力-有效塑性应变曲线；

其中σ_x，分别为用材LR平面内受拉或受压试验中沿加载方向的试验应力和试验塑性应变，和分别为该模型中的有效应力和有效塑性应变，θ为用材偏轴受拉或受压试验中的偏轴角度，偏轴角度为L轴与受力轴间的夹角，h(θ)表示中间过渡计算量，通过试值法确定a₁₁，a₂₂，b₁₁，b₂₂的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R²最大，幂函数主曲线形式见公式(5)；

B.RT平面屈服方程：

将SunChen模型与能够描述用材RT平面内受力响应的力学假设相结合，得到RT平面屈服方程；

其中a₄₄为RT平面内用材塑性特征参数，k₂是与用材应力有关的参数，a₄₄通过中国古建筑木结构修缮替换件用材材料偏轴受压应力-应变曲线获取，将应力-应变曲线转换为应力-塑性应变曲线；使用如下公式获取有效应力-有效塑性应变曲线；

其中Δ_x，为用材RT平面内受压试验中沿加载方向的试验应力和试验塑性应变，和分别为该模型中的有效应力和有效塑性应变，φ为用材偏轴受压试验中的偏轴角度，偏轴角度为R轴与受力轴间的夹角；i(φ)表示中间过渡计算量，通过试值法确定a₄₄的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R²最大，幂函数主曲线形式见公式(8)；

C.TL平面屈服方程：

将SunChen模型的衍生模型与能够描述用材TL平面内受力响应的力学假设相结合，得到TL平面屈服方程；

其中c₁₁，c₂₂，d₁₁，d₂₂为TL平面内用材塑性特征参数，ν为TL平面内用材弹性特征参数，k₃是与用材应力有关的参数，c₁₁，c₂₂，d₁₁，d₂₂通过中国古建筑木结构修缮替换件用材材料偏轴受拉和受压应力-应变曲线获取，将应力-应变关系转换为应力-塑性应变曲线，使用如下公式获取有效应力-有效塑性应变曲线：

其中α_x，为用材TL平面内受拉或受压试验中沿加载方向的试验应力和试验塑性应变，和分别为该模型中的有效应力和有效塑性应变，γ为用材偏轴受拉或受压试验中的偏轴角度，偏轴角度为T轴与受力轴间的夹角，j(γ)表示中间过渡计算量，通过试值法确定c₁₁，c₂₂，d₁₁，d₂₂的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R²最大，幂函数主曲线形式见公式(11)；

(3)用材受拉：中国古建筑木结构用材LR平面、RT平面和TL平面受拉时均出现脆性断裂，整体应力-应变曲线呈现线性，LR平面屈服方程、RT平面屈服方程和TL平面屈服方程描述的屈服面上的点为应力-应变曲线非线性阶段起始点或线性阶段终点，因此用材受拉状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零，采用如下方程进行描述：

假设第n步迭代得到的应力张量为n+1步中给出一个新的应变增量dε_kl，计算得到试探应力：

将试探应力代入屈服方程，屈服方程根据变形出现的平面拟定，分别为LR平面屈服方程、RT平面屈服方程和TL平面屈服方程；若则说明用材仍处于弹性阶段，其中n'＝1，2，3，n'＝1时对应LR平面屈服方程中的k₁，n'＝2时对应RT平面屈服方程中的k₂，n'＝3时对应TL平面屈服方程中的k₃，令作为新的应力张量，按照公式(12)进行循环；若则说明用材已经进入屈服阶段：用材处于受拉状态时发生断裂，保留该应变增量，应变增量步后对应的应力为零；

(4)用材受剪：中国古建筑木结构用材LR平面、RT平面和TL平面受剪时均出现脆性断裂，整体应力-应变曲线呈现线性，LR平面屈服方程、RT平面屈服方程和TL平面屈服方程描述的屈服面上的点为应力-应变曲线非线性阶段起始点或线性阶段终点，因此用材受剪状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零，同样采用公式(12)进行描述，同样将试探应力代入屈服方程，屈服方程根据变形出现的平面拟定，分别为LR平面屈服方程、RT平面屈服方程和TL平面屈服方程，若则说明用材仍处于弹性阶段，令作为新的应力张量，按照公式(12)进行循环，若则说明用材已经进入屈服阶段：用材处于受剪状态时发生断裂，保留该应变增量，应变增量步后对应的应力为零；

(5)用材LR平面受压：用材LR平面受压，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定；用材LR平面受压时的应力状态达到公式(3)、公式(4)和公式(5)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：

假设第n步迭代得到的应力张量为n+1步中给出一个新的应变增量dε_ij，根据公式(12)计算得到试探应力若在LR平面中，k_n'＝k₁，则说明用材已经进入屈服阶段，该试探应力不是真实用材应力，需要将超出屈服面的部分减去：将总应变增量中的塑性应变增量引起的虚假应力从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，具体方程如下：

总应变增量dε_ij由弹性应变增量和塑性应变增量之和进行表达，

若n+1步的真实应力为：

其中为塑性应变增量，公式(14)能够使用材应力始终处于屈服面上，实现理想塑性流动状态，其中塑性应变增量根据相关联的流动法则求解：

其中△λ为一致性系数，根据公式(15)，LR平面上的塑性应变表示为：

其中为用材沿L轴方向的塑性正应变增量，为用材位于LR平面内的塑性剪应变增量，为用材沿R轴方向的塑性正应变增量，△λ_LR求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：

其中为第n步试探应力代入公式(3)的结果，而为第n+1步试探应力代入公式(3)的结果；

故此

其中为与的差值；

将上式按照LR平面屈服方程进行一阶泰勒展开，得到：

其中△σ₁₁为沿着用材L轴方向的正应力增量，△σ₂₂为沿着用材R轴方向的正应力增量，△σ₁₂为位于用材LR平面内的剪应力增量；

其中M、X、Y、Z为运算过程中的过渡量；根据径向返回算法，公式(19)中的应力增量表述为：

其中为沿着用材L轴方向的试探正应力增量，为沿着用材R轴方向的试探正应力增量，为位于用材LR平面内的试探剪应力增量；

将公式(20)和公式(21)代入公式(19)，△λ_LR被求解：

在总应变达到应力软化阈值时，应力软化阈值为弹性总应变阈值的1.2倍，出现应力软化现象；

(6)用材RT平面受压：用材RT平面受压，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定，用材RT平面受压时的应力状态达到公式(6)、公式(7)和公式(8)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：

假设第n步迭代得到的应力张量为n+1步中给出一个新的应变增量dξ_ij，根据公式(12)计算得到试探应力若则说明用材已经进入屈服阶段，在RT平面中，k_n'＝k₂，该试探应力不是真实用材应力，需要将超出屈服面的部分减去：将总应变增量中的塑性应变引起的虚假应力增量从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，总应变增量由弹性应变增量和塑性应变增量之和进行表达，见公式(13)，真实应力计算见公式(14)，相关联塑性流动法则见公式(15)，RT平面内的塑性应变表示为：

其中用材沿R轴方向的塑性正应变增量，为用材位于RT平面内的塑性剪应变增量，为用材沿T轴方向的塑性正应变增量，△λ_RT的求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：

其中为第n步试探应力代入公式(6)的结果，而为第n+1步试探应力代入公式(6)的结果；故此：

其中为与的差值，将上式按照LR平面屈服方程进行一阶泰勒展开，得到：

其中△△₂₂为沿着用材R轴方向的正应力增量，△△₃₃为沿着用材T轴方向的正应力增量，△△₂₃为位于用材RT平面内的剪应力增量；

其中U、V、W为运算过程中的过渡量，公式(29)中的应力增量表述为：

其中为沿着用材R轴方向的试探正应力增量，为沿着用材T轴方向的试探正应力增量，为位于用材RT平面内的试探剪应力增量；

将公式(30)和公式(31)代入公式(29)，△λ_RT被求解：

在总应变达到应力硬化阈值时，应力硬化阈值为弹性总应变阈值的2倍，用材出现应力硬化现象，木材领域统称为二次强化现象；

(7)用材TL平面受压：用材TL平面受压，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定；用材RT平面受压时的应力状态达到公式(9)、公式(10)和公式(11)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：

假设第n步迭代得到的应力张量为n+1步中给出一个新的应变增量dβ_ij，根据公式(12)计算得到试探应力若在TL平面中，k_n'＝k₃，则说明用材已经进入屈服阶段，该试探应力不是真实用材应力，需要将超出屈服面的部分减去：将总应变增量中的塑性应变引起的虚假应力增量从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，总应变增量由弹性应变增量和塑性应变增量之和进行表达，见公式(13)，真实应力计算见公式(14)，相关联塑性流动法则见公式(15)，TL平面内的塑性应变表示为：

其中为用材沿T轴方向的塑性正应变增量，为用材位于TL平面内的塑性剪应变增量，为用材沿L轴方向的塑性正应变增量，△λ_TL求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：

其中为第n步试探应力代入公式(9)的结果，而为第n+1步试探应力代入公式(9)的结果；

故此：

其中为与的差值；

将上式按照TL平面屈服方程进行一阶泰勒展开，得到：

其中△α₁₁为沿着用材L轴方向的正应力增量，△α₃₃为沿着用材T轴方向的正应力增量，△α₁₃为位于用材TL平面内的剪应力增量；

其中：

其中L、G、H、I为运算过程中的过渡量，公式(40)中的应力增量表述为：

其中为沿着用材T轴方向的试探正应力增量，为沿着用材L轴方向的试探正应力增量，为位于用材TL平面内的试探剪应力增量；

将公式(41)和公式(42)代入公式(40)，△λ_TL被求解

在总应变达到应力硬化阈值时，应力硬化阈值应变为弹性总应变的2倍，应力出现硬化，木材领域统称为二次强化现象。