CN112541285A - 一种适于中国古建筑木结构用材本构关系的数值模拟方法 - Google Patents

一种适于中国古建筑木结构用材本构关系的数值模拟方法 Download PDF

Info

Publication number
CN112541285A
CN112541285A CN202011256618.6A CN202011256618A CN112541285A CN 112541285 A CN112541285 A CN 112541285A CN 202011256618 A CN202011256618 A CN 202011256618A CN 112541285 A CN112541285 A CN 112541285A
Authority
CN
China
Prior art keywords
stress
plane
strain
yield
increment
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011256618.6A
Other languages
English (en)
Other versions
CN112541285B (zh
Inventor
杨娜
李天昊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Jiaotong University
Original Assignee
Beijing Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Jiaotong University filed Critical Beijing Jiaotong University
Priority to CN202011256618.6A priority Critical patent/CN112541285B/zh
Publication of CN112541285A publication Critical patent/CN112541285A/zh
Application granted granted Critical
Publication of CN112541285B publication Critical patent/CN112541285B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明为一种适于中国古建筑木结构用材本构关系的数值模拟方法，主要分为(1)弹性阶段：将用材视为正交各项异性材料；(2)屈服方程：以宏观唯象的方式建立能够描述中国古建筑木结构用材累计损伤的屈服方程；(3)用材受拉状态：中国古建筑木结构用材在LR、RT、TL平面受拉时均出现脆性断裂，屈服时刻应变增量末端对应应力为零；(4)用材受剪状态：中国古建筑木结构用材在LR、RT、TL平面受剪时均出现脆性断裂，屈服时刻应变增量末端对应应力为零；(5)用材LR平面受压状态：引入径向返回算法、损伤因子和弹性应变能进行描述；(6)用材RT、TL平面受压状态：引入径向返回算法、背应力增量和等效塑性应变进行描述。

Description

一种适于中国古建筑木结构用材本构关系的数值模拟方法
技术领域
本发明涉及中国古建筑木结构用材力学性能领域，尤其涉及中国古建筑木结构复杂三维受力状态下本构关系的数值模拟。
背景技术
中国古建筑木结构属于东方独立系统，它区别于世界所有建筑体系，其具有独特的结构形式、布置规模。中国古建筑木结构是极特殊、极长寿、极体面的，它也是中华民族精神的体现。遗憾的是，在漫长的生命周期中中国古建筑木结构受到环境因素和人为因素的影响，而产生结构层次、构件层次、材料层次的损伤，因此对于中国古建筑木结构的保护迫在眉睫。近年来，越来越多的工程人员和研究人员投入到中国古建筑木结构保护领域，勘察、检测、监测、修缮、机理研究等诸多方面起头并进，共同服务于中国古建筑木结构的保护。
由于中国古建筑木结构极其珍贵，对于该对象的研究多基于数值模拟方法。目前大型有限元分析软件ABAQUS被大多数研究人员使用，汇总相关研究发现，在使用大型有限元分析软件ABAQUS进行数值模拟时，该软件未提供能够适应中国古建筑木结构用材的模型。中国古建筑木结构用材为木材料，目前已有研究人员使用Fortran语言编写适合木材料的本构关系，嵌入大型有限元分析软件ABAQUS的UMAT或VUMAT材料子程序中实现对现代木结构的数值分析。但上述木材料本构关系所采用的屈服方程均来自于理想正交各项异性材料，例如Hill屈服准则、Hoffman屈服准则、Tsai-Wu屈服准则、Norris屈服准则、Hashin屈服准则和Yamad-Sun屈服准则。上述屈服准则中的应力分量系数具有严格的物理意义和比例关系，即上述屈服准则可以很好的描述理想状态下木材料的力学响应。众所周知，木材料宏观力学响应往往与微观构造相关，新木材料可以假定为理性正交各向异性材料。中国古建筑木结构用材，属于老旧木材料，其漫长的生命周期使其在复杂的工作环境中逐渐出现累计损伤，例如木材料的腐朽、虫蛀、裂缝、孔洞等。中国古建筑木结构用材微观构造与新木材完全不同，故此宏观力学响应也存在差别，上述子程序的屈服方程无法适用于中国古建筑木结构用材。工程人员和研究人员需要一种能够描述中国古建筑木结构用材在复杂受力状态下的屈服方程，该屈服方程应能够描述中国古建筑木结构用材累计损伤的特征。已有学者引入严格的损伤力学方程，但是复杂的工作环境带来的随机性使得损伤力学方程中的系数难以确定，引入严格的损伤力学方程会极大提高数值模拟的计算成本。
针对上述问题，本发明公开的方法使用宏观唯象的屈服方程，屈服方程中的参数由中国古建筑木结构修缮替换件用材材性试验结果确定，使用曲线拟合的方式描述用材的累计损伤。参数没有严格的力学定义，应力分量间的比例关系与用材真实响应一致。结合径向返回算法、损伤因子、弹性应变能、背应力增量和等效塑性应变，获取能够描述中国古建筑木结构用材复杂三维受力状态下力学响应的本构关系。本构关系中的关键参数已展示，可通过Fortran语言编写，嵌入大型有限元分析软件ABAQUS的UMAT材料子程序中，应用于实际数值模拟中。
发明内容
针对现有技术中存在的缺陷，本发明的目的在于提供一种适于中国古建筑木结构用材本构关系的数值模拟方法。
为达到以上目的，本发明采取的技术方案是：
本发明公开一种适于中国古建筑木结构用材本构关系的数值模拟方法，包括以下步骤：
(1)弹性阶段：中国古建筑木结构用材弹性阶段用以下公式描述：
Figure BDA0002773316740000031
其中Dijkl为用材弹性刚度矩阵，i，j，k，l＝1，2，3，因为用材被公认为正交各项异性材料，并存在3个方向的材料主轴，因此弹性刚度矩阵中存在零值；σ11，σ22，σ33为用材的主轴正应力，σ12，σ23，σ13为用材的主面剪应力；ε11，ε22，ε33为用材的主轴正应变，ε12，ε23，ε13为用材的主面剪应变；另外，值得注意的是，在后续描述中，LR平面内的应力使用σ进行表达，应变使用ε进行表达；RT平面内的应力使用Δ进行表达，应变使用ξ进行表达；TL平面内的应力使用α进行表达，应变使用β进行表达；
用材弹性刚度矩阵Dijkl中的参数的表达具体如公式(2)所示；
Figure BDA0002773316740000032
其中E1，E2，E3为用材主轴弹性模量，G12，G23，G13为用材主面剪切模量，△表示中间过渡计算量，ν12，ν23，ν13，ν21，ν32，ν31为用材主面泊松比；
(2)屈服方程：
A.LR平面屈服方程：
将Sun&Chen模型的衍生模型与能够描述用材LR平面内受力响应的力学假设相结合，得到LR平面屈服方程：
Figure BDA0002773316740000041
其中a11，a22，b11，b22为LR平面内用材塑性特征参数，μ为LR平面内用材弹性特征参数，k1是与用材应力有关的参数，a11，a22，b11，b22通过中国古建筑木结构修缮替换件用材材料偏轴受拉和受压应力-应变曲线获取，将应力-应变曲线转换为应力-塑性应变曲线，使用如下公式获取有效应力-有效塑性应变曲线；
Figure BDA0002773316740000042
其中σx，
Figure BDA0002773316740000043
分别为用材LR平面内受拉或受压试验中沿加载方向的试验应力和试验塑性应变，
Figure BDA0002773316740000044
和
Figure BDA0002773316740000045
分别为该模型中的有效应力和有效塑性应变，θ为用材偏轴受拉或受压试验中的偏轴角度，偏轴角度为L轴与受力轴间的夹角，h(θ)表示中间过渡计算量，通过试值法确定a11，a22，b11，b22的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R2最大，幂函数主曲线形式见公式(5)；
Figure BDA0002773316740000051
B.RT平面屈服方程：
将Sun&Chen模型与能够描述用材RT平面内受力响应的力学假设相结合，得到RT平面屈服方程，值得注意的是，为区分符号，RT平面内的应力使用Δ进行表达，应变使用ξ进行表达；
Figure BDA0002773316740000052
其中a44为RT平面内用材塑性特征参数，k2是与用材应力有关的参数，a44通过中国古建筑木结构修缮替换件用材材料偏轴受压应力-应变曲线获取，将应力-应变曲线转换为应力-塑性应变曲线；使用如下公式获取有效应力-有效塑性应变曲线；
Figure BDA0002773316740000053
其中Δx，
Figure BDA0002773316740000054
为用材RT平面内受压试验中沿加载方向的试验应力和试验塑性应变，
Figure BDA0002773316740000055
和
Figure BDA0002773316740000056
分别为该模型中的有效应力和有效塑性应变，φ为用材偏轴受压试验中的偏轴角度，偏轴角度为R轴与受力轴间的夹角；i(φ)表示中间过渡计算量，通过试值法确定a44的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R2最大，幂函数主曲线形式见公式(8)；
Figure BDA0002773316740000061
C.TL平面屈服方程：
将Sun&Chen模型的衍生模型与能够描述用材TL平面内受力响应的力学假设相结合，得到TL平面屈服方程；另外，值得注意的是，为区分符号，TL平面内的应力使用α进行表达，应变使用β进行表达；
Figure BDA0002773316740000062
其中c11，c22，d11，d22为TL平面内用材塑性特征参数，ν为TL平面内用材弹性特征参数，k3是与用材应力有关的参数，c11，c22，d11，d22通过中国古建筑木结构修缮替换件用材材料偏轴受拉和受压应力-应变曲线获取，将应力-应变关系转换为应力-塑性应变曲线，使用如下公式获取有效应力-有效塑性应变曲线：
Figure BDA0002773316740000063
其中αx，
Figure BDA0002773316740000064
为用材TL平面内受拉或受压试验中沿加载方向的试验应力和试验塑性应变，
Figure BDA0002773316740000065
和
Figure BDA0002773316740000066
分别为该模型中的有效应力和有效塑性应变，γ为用材偏轴受拉或受压试验中的偏轴角度，偏轴角度为T轴与受力轴间的夹角，j(γ)表示中间过渡计算量，通过试值法确定c11，c22，d11，d22的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R2最大，幂函数主曲线形式见公式(11)；
Figure BDA0002773316740000071
(3)用材受拉：中国古建筑木结构用材LR平面、RT平面和TL平面受拉时均出现脆性断裂，整体应力-应变曲线呈现线性，LR平面屈服方程、RT平面屈服方程和TL平面屈服方程描述的屈服面上的点为应力-应变曲线非线性阶段起始点或线性阶段终点，因此用材受拉状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零，采用如下方程进行描述：
假设第n步迭代得到的应力张量为
Figure BDA0002773316740000072
n+1步中给出一个新的应变增量dεkl，计算得到试探应力：
Figure BDA0002773316740000073
将试探应力代入屈服方程，屈服方程根据变形出现的平面拟定，分别为LR平面屈服方程、RT平面屈服方程和TL平面屈服方程；若
Figure BDA0002773316740000074
则说明用材仍处于弹性阶段，其中n'＝1，2，3，n'＝1时对应LR平面屈服方程中的k1，n'＝2时对应RT平面屈服方程中的k2，n'＝3时对应TL平面屈服方程中的k3，令
Figure BDA0002773316740000075
作为新的应力张量，按照公式(12)进行循环；若
Figure BDA0002773316740000076
则说明用材已经进入屈服阶段，即用材处于受拉状态时发生断裂，保留该应变增量，应变增量步后对应的应力为零；
(4)用材受剪：中国古建筑木结构用材LR平面、RT平面和TL平面受剪时均出现脆性断裂，整体应力-应变曲线呈现线性，LR平面屈服方程、RT平面屈服方程和TL平面屈服方程描述的屈服面上的点为应力-应变曲线非线性阶段起始点或线性阶段终点，因此用材受剪状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零，本发明同样采用公式(12)进行描述，同样将试探应力代入屈服方程，屈服方程根据变形出现的平面拟定，分别为LR平面屈服方程、RT平面屈服方程和TL平面屈服方程，若
Figure BDA0002773316740000081
则说明用材仍处于弹性阶段，令
Figure BDA0002773316740000082
作为新的应力张量，按照公式(12)进行循环，若
Figure BDA0002773316740000083
则说明用材已经进入屈服阶段，即用材处于受剪状态时发生断裂，保留该应变增量，应变增量步后对应的应力为零；
(5)用材LR平面受压：用材LR平面受压(即顺纹受压)，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定；用材LR平面受压时的应力状态达到公式(3)、公式(4)和公式(5)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：假设第n步迭代得到的应力张量为
Figure BDA0002773316740000084
n+1步中给出一个新的应变增量dεij，根据公式(12)计算得到试探应力
Figure BDA0002773316740000085
若
Figure BDA0002773316740000086
在LR平面中，kn'＝k1，则说明用材已经进入屈服阶段，该试探应力不是真实用材应力，需要将超出屈服面的部分减去，即将总应变增量中的塑性应变增量引起的虚假应力从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，具体方程如下：
总应变增量dεij由弹性应变增量
Figure BDA0002773316740000091
和塑性应变增量
Figure BDA0002773316740000092
之和进行表达，
Figure BDA0002773316740000093
若
Figure BDA0002773316740000094
n+1步的真实应力为：
Figure BDA0002773316740000095
其中
Figure BDA0002773316740000096
为塑性应变增量，公式(14)能够使用材应力始终处于屈服面上，实现理想塑性流动状态，其中塑性应变增量
Figure BDA0002773316740000097
根据相关联的流动法则求解：
Figure BDA0002773316740000098
其中△λ为一致性系数，根据公式(15)，LR平面上的塑性应变表示为：
Figure BDA0002773316740000099
其中
Figure BDA00027733167400000910
为用材沿L轴方向的塑性正应变增量，
Figure BDA00027733167400000911
为用材位于LR平面内的塑性剪应变增量，
Figure BDA00027733167400000912
为用材沿R轴方向的塑性正应变增量，△λLR求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：
Figure BDA00027733167400000913
其中
Figure BDA00027733167400000914
为第n步试探应力代入公式(3)的结果，而
Figure BDA00027733167400000915
为第n+1步试探应力代入公式(3)的结果；
故此
Figure BDA0002773316740000101
其中
Figure BDA0002773316740000102
为
Figure BDA0002773316740000103
与
Figure BDA0002773316740000104
的差值；
将上式按照LR平面屈服方程进行一阶泰勒展开，得到：
Figure BDA0002773316740000105
其中△σ11为沿着用材L轴方向的正应力增量，△σ22为沿着用材R轴方向的正应力增量，△σ12为位于用材LR平面内的剪应力增量；
Figure BDA0002773316740000106
其中M、X、Y、Z仅为运算过程中的过渡量，以缩减后续公式的长度；根据径向返回算法，公式(19)中的应力增量表述为：
Figure BDA0002773316740000107
其中
Figure BDA0002773316740000111
为沿着用材L轴方向的试探正应力增量，
Figure BDA0002773316740000112
为沿着用材R轴方向的试探正应力增量，
Figure BDA0002773316740000113
为位于用材LR平面内的试探剪应力增量；
将公式(20)和公式(21)代入公式(19)，△λLR被求解：
Figure BDA0002773316740000114
在总应变达到应力软化阈值时，应力软化阈值
Figure BDA0002773316740000115
为弹性总应变阈值
Figure BDA0002773316740000116
的1.2倍(即
Figure BDA0002773316740000117
)，出现应力软化现象；本发明公开的方法引入损伤因子描述上述应力软化现象，使用应力逐步退化模型来模拟用材应力软化：
Figure BDA0002773316740000118
其中，σij为无损应力，即理想塑性流动末端应力，dLR为损伤因子，当用材没有损伤时，dLR＝0，当木材完全损伤时，dLR＝1，
Figure BDA0002773316740000119
为损伤后的应力；上述损伤因子一般与当前状态应变、应力和能量函数有关，本发明公开的方法中采用一种与总应变相关的损伤累计理论模型，该模型认为材料损伤与总应变和无损弹性模量相关，具体表示如下：
Figure BDA00027733167400001110
其中τ为总应变对应的应变能，εij、εkl为LR平面对应的应变；用材LR平面内应力软化对应的损伤因子如下：
dLR＝1-exp[-(τLR-τLR,0)/τLR,0] (25)
其中，τLR为LR平面内总应变对应的应变能，而τLR,0为应力软化阈值对应的应变能；
(6)用材RT平面受压：首先需要说明的是，为区分符号，RT平面内的应力使用Δ进行表达，应变使用ξ进行表达，上下角标的物理意义保持不变。用材RT平面受压(即横纹受压)，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定，用材RT平面受压时的应力状态达到公式(6)、公式(7)和公式(8)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：假设第n步迭代得到的应力张量为
Figure BDA0002773316740000121
n+1步中给出一个新的应变增量dξij，根据公式(12)计算得到试探应力
Figure BDA0002773316740000122
若
Figure BDA0002773316740000123
则说明用材已经进入屈服阶段，在RT平面中，kn'＝k2，该试探应力不是真实用材应力，需要将超出屈服面的部分减去，即将总应变增量中的塑性应变引起的虚假应力增量从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，总应变增量可由弹性应变增量和塑性应变增量之和进行表达，见公式(13)，真实应力计算见公式(14)，相关联塑性流动法则见公式(15)，RT平面内的塑性应变表示为：
Figure BDA0002773316740000124
其中
Figure BDA0002773316740000125
用材沿R轴方向的塑性正应变增量，
Figure BDA0002773316740000126
为用材位于RT平面内的塑性剪应变增量，
Figure BDA0002773316740000127
为用材沿T轴方向的塑性正应变增量，△λRT的求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：
Figure BDA0002773316740000128
其中
Figure BDA0002773316740000129
为第n步试探应力代入公式(6)的结果，而
Figure BDA00027733167400001210
为第n+1步试探应力代入公式(6)的结果；故此：
Figure BDA0002773316740000131
其中
Figure BDA0002773316740000132
为
Figure BDA0002773316740000133
与
Figure BDA0002773316740000134
的差值，将上式按照LR平面屈服方程进行一阶泰勒展开，得到：
Figure BDA0002773316740000135
其中△△22为沿着用材R轴方向的正应力增量，△△33为沿着用材T轴方向的正应力增量，△△23为位于用材RT平面内的剪应力增量；
Figure BDA0002773316740000136
其中U、V、W仅为运算过程中的过渡量，以缩减后续公式的长度；根据径向返回算法，公式(29)中的应力增量表述为：
Figure BDA0002773316740000137
其中
Figure BDA0002773316740000138
为沿着用材R轴方向的试探正应力增量，
Figure BDA0002773316740000139
为沿着用材T轴方向的试探正应力增量，
Figure BDA00027733167400001310
为位于用材RT平面内的试探剪应力增量；
将公式(30)和公式(31)代入公式(29)，△λRT被求解：
Figure BDA0002773316740000141
在总应变达到应力硬化阈值时，应力硬化阈值
Figure BDA0002773316740000142
为弹性总应变阈值
Figure BDA0002773316740000143
的2倍(即
Figure BDA0002773316740000144
)，用材出现应力硬化现象，木材领域统称为二次强化现象，这种现象是因为用材由大量纤维沿纵向排列，形成了纤维束状材料。在横向受压时，细胞壁被压扁，而出现用材密实现象，导致用材出现塑性流动后并不发生应力软化破坏，而因为压实而出现二次应力硬化。
本发明公开的方法采用设置最初屈服面和最终屈服面的方法来描述上述问题，最初屈服面为用材进入塑性流动阶段的起点，对应的最初的屈服面方程，见公式(6)；引入NRT表示用材RT平面受压时，最终屈服强度和初始屈服强度之比，故此最终屈服面公式如下：
Figure BDA0002773316740000145
其中，NRT-2为RT平面受压时，R轴的最终屈服强度和初始屈服强度之比，数值通过试验曲线确定，NRT-3为RT平面受力时，T轴的最终屈服强度和初始屈服强度之比，数值通过试验曲线确定，k2是与用材应力有关的参数，通过公式(6)确定，屈服面转移通过背应力增量进行描述，背应力增长与应变的增量相关，总背应力增量是最终屈服应力和初始屈服应力之差，RT平面背应力构造方式如下：
Figure BDA0002773316740000146
其中，cRT为RT平面内背应力增长的快慢程度，由RT平面主轴受压试验获取，GRT-2和GRT-3为RT平面二次硬化屈服面转移约束方程，其中ω22为R轴方向的背应力，ω33为T轴方向的背应力，△ω22为R轴方向的背应力增量，△ω33为T轴方向的背应力增量，防止增加的屈服面超过最终屈服面；
Figure BDA0002773316740000151
通过公式(35)可知，最初GRT-2和GRT-3均为1，对背应力增长无限制作用，当背应力总量达到最大值时，即最终屈服强度与初始屈服强度之差，背应力停止增长，公式(35)均为0，而
Figure BDA0002773316740000152
和
Figure BDA0002773316740000153
表征背应力的增长方向，其中
Figure BDA0002773316740000154
和
Figure BDA0002773316740000155
为用材最终屈服强度，用试验获取，
Figure BDA0002773316740000156
为等效塑性应变增量，由如下公式获取；
Figure BDA0002773316740000157
(7)用材TL平面受压：首先需要说明的是，为区分符号，TL平面内的应力使用α进行表达，应变使用β进行表达，上下角标的物理意义保持不变。用材TL平面受压(即横纹受压)，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定；用材RT平面受压时的应力状态达到公式(9)、公式(10)和公式(11)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：假设第n步迭代得到的应力张量为
Figure BDA0002773316740000158
n+1步中给出一个新的应变增量dβij，根据公式(12)计算得到试探应力
Figure BDA0002773316740000161
若
Figure BDA0002773316740000162
在TL平面中，kn'＝k3，则说明用材已经进入屈服阶段，该试探应力不是真实用材应力，需要将超出屈服面的部分减去，即将总应变增量中的塑性应变引起的虚假应力增量从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态。总应变增量可由弹性应变增量和塑性应变增量之和进行表达，见公式(13)，真实应力计算见公式(14)，相关联塑性流动法则见公式(15)，TL平面内的塑性应变表示为：
Figure BDA0002773316740000163
其中
Figure BDA0002773316740000164
为用材沿T轴方向的塑性正应变增量，
Figure BDA0002773316740000165
为用材位于TL平面内的塑性剪应变增量，
Figure BDA0002773316740000166
为用材沿L轴方向的塑性正应变增量，△λTL求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：
Figure BDA0002773316740000167
其中
Figure BDA0002773316740000168
为第n部试探应力代入公式(9)的结果，而
Figure BDA0002773316740000169
为第n+1步试探应力代入公式(9)的结果；
故此：
Figure BDA00027733167400001610
其中
Figure BDA00027733167400001611
为
Figure BDA00027733167400001612
与
Figure BDA00027733167400001613
的差值；
将上式按照TL平面屈服方程进行一阶泰勒展开，得到：
Figure BDA00027733167400001614
其中△α11为沿着用材L轴方向的正应力增量，△α33为沿着用材T轴方向的正应力增量，△α13为位于用材TL平面内的剪应力增量；
其中：
Figure BDA0002773316740000171
其中L、G、H、I仅为运算过程中的过渡量，以缩减后续公式的长度；根据径向返回算法，公式(40)中的应力增量表述为：
Figure BDA0002773316740000172
其中
Figure BDA0002773316740000173
为沿着用材T轴方向的试探正应力增量，
Figure BDA0002773316740000174
为沿着用材L轴方向的试探正应力增量，
Figure BDA0002773316740000175
为位于用材TL平面内的试探剪应力增量；
将公式(41)和公式(42)代入公式(40)，△λTL被求解
Figure BDA0002773316740000176
在总应变达到应力硬化阈值时，应力硬化阈值应变
Figure BDA0002773316740000181
为弹性总应变
Figure BDA0002773316740000182
的2倍(即
Figure BDA0002773316740000183
)，应力出现硬化，木材领域统称为二次强化现象。本发明公开的方法采用设置最初屈服面和最终屈服面的方法来描述上述问题，最初屈服面为用材进入理想塑性流动阶段的起点，对应的最初屈服面方程，见公式(9)；引入NTL表示用材TL平面受压时，最终屈服强度和初始屈服强度之比，故此最终屈服面公式如下
Figure BDA0002773316740000184
其中，NTL-3为TL平面受力时，T轴的最终屈服强度和初始强度之比，数值通过试验曲线确定；NTL-1为TL平面受力时，L轴的最终屈服强度和初始强度之比，数值通过试验曲线确定；k3是与用材应力有关的参数，通过公式(9)确定。屈服面转移通过背应力增量进行描述，背应力增长与应变的增量相关，总背应力增量即是最终屈服应力和初始屈服应力之差。TL平面背应力构造方式如下：
Figure BDA0002773316740000185
其中，cTL为TL平面内背应力增长的快慢程度，由TL平面主轴受压试验获取；GTL-3和GTL-1为TL平面二次硬化屈服面转移约束方程，其中ρ11为L轴方向的背应力，ρ33为T轴方向的背应力，△ρ11为L轴方向的背应力增量，△ρ33为T轴方向的背应力增量，防止增加的屈服面超过最终屈服面；
Figure BDA0002773316740000191
通过公式(46)可知，最初GTL-3和GTL-1均为1，对背应力增长无限制作用。当背应力总量达到最大值时，即最终屈服强度与初始屈服强度之差，背应力停止增长，公式(46)均为0。而
Figure BDA0002773316740000192
和
Figure BDA0002773316740000193
表征背应力的增长方向，其中
Figure BDA0002773316740000194
和
Figure BDA0002773316740000195
为用材最终屈服强度，用试验获取。
Figure BDA0002773316740000196
为等效塑性应变增量，具体表示如下所示：
Figure BDA0002773316740000197
附图说明
本发明有如下附图：
图1-1为中国古建筑木结构基本形式；
图1-2为中国古建筑木结构用材主轴示意图；
图2为中国古建筑木结构用材受拉应力-应变示意图(LR/RT/TL)；
图3为中国古建筑木结构用材受剪应力-应变示意图(LR/RT/TL)；
图4为中国古建筑木结构用材受压应力-应变示意图(LR)；
图5为中国古建筑木结构用材受压应力-应变示意图(RT/TL)。
具体实施方式
以下结合附图1-5对本发明作进一步详细说明。
本方法能够描述中国古建筑木结构用材在三维受力状态下的力学响应，中国古建筑木结构基本形式见图1-1。将上述本构关系分为(1)弹性阶段；(2)屈服方程；(3)用材受拉；(4)用材受剪；(5)用材LR平面受压；(6)用材RT平面和TL平面受压。本发明展示了上述方法中的关键物理量，便于嵌入有限元分析软件ABAQUS中，以实现数值模拟功能。
弹性阶段中，中国古建筑木结构用材可视为正交各向异性材料，中国古建筑木结构用材的三个主轴方向见图1-2，用材受拉、受压、受剪、受弯和受扭等应力-应变关系可由6×6的正交各向异性材料刚度矩阵进行描述。
首先将用材简化为LR平面、RT平面和TL平面状态进行描述，上述平面状态下的屈服方程源于单向纤维复合材料领域的Sun&Chen模型和其衍生模型，模型增加了能够考虑用材受力响应的力学假设。结合中国古建筑木结构修缮替换件材料性能试验结果，使用宏观唯象的方法描述包含当前累计损伤的用材的屈服状态，避免了使用复杂的损伤力学方程。
由于上述唯象屈服方程关注的点为用材应力-应变曲线非线性阶段的起点。用材LR平面、RT平面和TL平面状态受拉时，均发生脆性断裂，应力-应变曲线呈现线性，故用材受拉时，应力水平达到屈服面后，用材断裂，故保留应变增量，而应变增量步对应的应力为零，用材在LR平面、RT平面和TL平面内受拉时的应力-应变曲线示意图见图2。
由于上述唯象屈服方程关注的点为用材应力-应变曲线非线性阶段起点。用材LR平面、RT平面和TL平面状态受剪时，均发生脆性断裂，应力-应变曲线呈现线性，故用材受剪时，应力水平达到屈服面后，用材断裂，故保留应变增量，而应变增量步对应应力为零。用材在LR平面、RT平面和TL平面内受剪时的应力-应变曲线示意图见图3。
用材处于LR平面内受压时，即用材处于顺纹受压状态，用材弹性段使用弹性正交各向异性材料矩阵描述，屈服方程同样采用宏观唯象屈服方程。屈服后，基于径向返回算法实现理想塑性流动阶段。总应变达到顺纹弹性总应变的1.2倍时，引入与用材弹性应变能相关的损伤因子描述用材的应力软化现象。用材在LR平面内受压时的应力-应变曲线示意图见图4。
用材处于RT平面和TL平面状态受压时，即用材处于横纹受压状态，用材弹性段使用弹性正交各向异性材料矩阵描述，屈服方程同样采用宏观唯象屈服方程。屈服后，基于径向返回算法实现理想塑性流动阶段。总应变达到横纹弹性总应变的2倍时，引入与用材等效塑性应变相关的背应力增量，使初始屈服面转移至最终屈服面，以描述应力硬化现象。用材在RT平面和TL平面内受压时的应力-应变曲线示意图见图5。
上述过程可利用Fortran语言编写，嵌入有限元分析软件ABAQUS的UMAT材料子程序中，所需参数的求解方式在本发明公开方法中已展示，便于实际数值模拟使用。
本发明公开的适于中国古建筑木结构用材本构关系的数值模拟方法实施步骤如下：
(1)弹性阶段使用如下公式描述
Figure BDA0002773316740000221
(2)屈服方程：屈服方程分为LR平面、RT平面、TL平面，通过中国古建筑修缮替换件用材本构关系确定参数。
Figure BDA0002773316740000222
(3)用材受拉：用材受拉状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零。假设第n步迭代得到的应力张量为
Figure BDA0002773316740000223
n+1步中给出一个新的应变增量dεkl，计算得到试探应力
Figure BDA0002773316740000224
Figure BDA0002773316740000225
时用材处于受拉状态时发生断裂，保留该应变增量，应变增量步后对应的应力为零。
(4)用材受剪：用材受剪状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零。本发明同样可以采用公式(12)进行描述，同样将试探应力代入上述屈服方程，
Figure BDA0002773316740000231
则说明用材受剪发生断裂，保留该应变增量，应变增量步后对应的应力为零。
(5)用材LR平面受压：弹性阶段和屈服方程在步骤(1)、(2)中给出。屈服后进入理想塑性流动阶段，塑流阶段使用相关联流动方程，一致性系数如下，△λLR可被求解
Figure BDA0002773316740000232
在总应变达到应力软化阈值时，使用应力逐步退化模型
Figure BDA0002773316740000233
其中
Figure BDA0002773316740000234
对于用材LR平面内应力软化对应的损伤因子
dLR＝1-exp[-(τLR-τLR,0)/τLR,0] (25)
(6)用材RT平面受压：弹性阶段和屈服方程在步骤(1)、(2)中给出。屈服后进入理想塑性流动阶段，塑流阶段使用相关联流动方程，一致性系数如下，△λRT可被求解
Figure BDA0002773316740000235
在总应变达到应力硬化阈值时，应力出现硬化，使用设置最初屈服面和最终屈服面的方法来描述上述问题，引入NRT表示用材RT平面受压时，最终屈服强度和初始屈服强度之比，故此最终屈服面公式如下
Figure BDA0002773316740000241
屈服面转移通过背应力增量进行描述，RT平面背应力构造方式如下：
Figure BDA0002773316740000242
其中
Figure BDA0002773316740000243
Figure BDA0002773316740000244
为等效塑性应变增量
Figure BDA0002773316740000245
(7)用材TL平面受压：弹性阶段和屈服方程在步骤(1)、(2)中给出。屈服后进入理想塑性流动阶段，塑流阶段使用相关联流动方程，一致性系数如下，△λTL可被求解
Figure BDA0002773316740000246
在总应变达到应力硬化阈值时，应力出现硬化，设置最初屈服面和最终屈服面的方法来描述上述问题，引入NTL表示用材TL平面受压时，最终屈服强度和初始屈服强度之比，故此最终屈服面公式如下
Figure BDA0002773316740000247
屈服面转移通过背应力增量进行描述，TL平面背应力构造方式如下：
Figure BDA0002773316740000251
其中
Figure BDA0002773316740000252
Figure BDA0002773316740000253
为等效塑性应变增量：
Figure BDA0002773316740000254
本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (4)

1.一种适于中国古建筑木结构用材本构关系的数值模拟方法，其特征在于，包括以下步骤：
(1)弹性阶段：中国古建筑木结构用材弹性阶段用以下公式描述：
Figure FDA0002773316730000011
其中Dijkl为用材弹性刚度矩阵，i，j，k，l＝1，2，3；σ11，σ22，σ33为用材的主轴正应力，σ12，σ23，σ13为用材的主面剪应力；ε11，ε22，ε33为用材的主轴正应变，ε12，ε23，ε13为用材的主面剪应变；在后续描述中，LR平面内的应力使用σ进行表达，应变使用ε进行表达；RT平面内的应力使用Δ进行表达，应变使用ξ进行表达；TL平面内的应力使用α进行表达，应变使用β进行表达；
用材弹性刚度矩阵Dijkl中的参数的表达具体如公式(2)所示；
Figure FDA0002773316730000012
其中E1，E2，E3为用材主轴弹性模量，G12，G23，G13为用材主面剪切模量，△表示中间过渡计算量，ν12，ν23，ν13，ν21，ν32，ν31为用材主面泊松比；
(2)屈服方程：
A.LR平面屈服方程：
将Sun&Chen模型的衍生模型与能够描述用材LR平面内受力响应的力学假设相结合，得到LR平面屈服方程：
Figure FDA0002773316730000021
其中a11，a22，b11，b22为LR平面内用材塑性特征参数，μ为LR平面内用材弹性特征参数，k1是与用材应力有关的参数，a11，a22，b11，b22通过中国古建筑木结构修缮替换件用材材料偏轴受拉和受压应力-应变曲线获取，将应力-应变曲线转换为应力-塑性应变曲线，使用如下公式获取有效应力-有效塑性应变曲线；
Figure FDA0002773316730000022
其中σx，
Figure FDA0002773316730000023
分别为用材LR平面内受拉或受压试验中沿加载方向的试验应力和试验塑性应变，
Figure FDA0002773316730000024
和
Figure FDA0002773316730000025
分别为该模型中的有效应力和有效塑性应变，θ为用材偏轴受拉或受压试验中的偏轴角度，偏轴角度为L轴与受力轴间的夹角，h(θ)表示中间过渡计算量，通过试值法确定a11，a22，b11，b22的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R2最大，幂函数主曲线形式见公式(5)；
Figure FDA0002773316730000031
B.RT平面屈服方程：
将Sun&Chen模型与能够描述用材RT平面内受力响应的力学假设相结合，得到RT平面屈服方程；
Figure FDA0002773316730000032
其中a44为RT平面内用材塑性特征参数，k2是与用材应力有关的参数，a44通过中国古建筑木结构修缮替换件用材材料偏轴受压应力-应变曲线获取，将应力-应变曲线转换为应力-塑性应变曲线；使用如下公式获取有效应力-有效塑性应变曲线；
Figure FDA0002773316730000033
其中Δx，
Figure FDA0002773316730000034
为用材RT平面内受压试验中沿加载方向的试验应力和试验塑性应变，
Figure FDA0002773316730000035
和
Figure FDA0002773316730000036
分别为该模型中的有效应力和有效塑性应变，φ为用材偏轴受压试验中的偏轴角度，偏轴角度为R轴与受力轴间的夹角；i(φ)表示中间过渡计算量，通过试值法确定a44的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R2最大，幂函数主曲线形式见公式(8)；
Figure FDA0002773316730000041
C.TL平面屈服方程：
将Sun&Chen模型的衍生模型与能够描述用材TL平面内受力响应的力学假设相结合，得到TL平面屈服方程；
Figure FDA0002773316730000042
其中c11，c22，d11，d22为TL平面内用材塑性特征参数，ν为TL平面内用材弹性特征参数，k3是与用材应力有关的参数，c11，c22，d11，d22通过中国古建筑木结构修缮替换件用材材料偏轴受拉和受压应力-应变曲线获取，将应力-应变关系转换为应力-塑性应变曲线，使用如下公式获取有效应力-有效塑性应变曲线：
Figure FDA0002773316730000043
其中αx，
Figure FDA0002773316730000044
为用材TL平面内受拉或受压试验中沿加载方向的试验应力和试验塑性应变，
Figure FDA0002773316730000045
和
Figure FDA0002773316730000046
分别为该模型中的有效应力和有效塑性应变，γ为用材偏轴受拉或受压试验中的偏轴角度，偏轴角度为T轴与受力轴间的夹角，j(γ)表示中间过渡计算量，通过试值法确定c11，c22，d11，d22的值，将上述有效应力-有效塑性应变曲线拟合成为一条幂函数主曲线求得最优解，拟合过程中的决定系数R2最大，幂函数主曲线形式见公式(11)；
Figure FDA0002773316730000051
(3)用材受拉：中国古建筑木结构用材LR平面、RT平面和TL平面受拉时均出现脆性断裂，整体应力-应变曲线呈现线性，LR平面屈服方程、RT平面屈服方程和TL平面屈服方程描述的屈服面上的点为应力-应变曲线非线性阶段起始点或线性阶段终点，因此用材受拉状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零，采用如下方程进行描述：
假设第n步迭代得到的应力张量为
Figure FDA0002773316730000052
n+1步中给出一个新的应变增量dεkl，计算得到试探应力：
Figure FDA0002773316730000053
将试探应力代入屈服方程，屈服方程根据变形出现的平面拟定，分别为LR平面屈服方程、RT平面屈服方程和TL平面屈服方程；若
Figure FDA0002773316730000054
则说明用材仍处于弹性阶段，其中n'＝1，2，3，n'＝1时对应LR平面屈服方程中的k1，n'＝2时对应RT平面屈服方程中的k2，n'＝3时对应TL平面屈服方程中的k3，令
Figure FDA0002773316730000055
作为新的应力张量，按照公式(12)进行循环；若
Figure FDA0002773316730000056
则说明用材已经进入屈服阶段：用材处于受拉状态时发生断裂，保留该应变增量，应变增量步后对应的应力为零；
(4)用材受剪：中国古建筑木结构用材LR平面、RT平面和TL平面受剪时均出现脆性断裂，整体应力-应变曲线呈现线性，LR平面屈服方程、RT平面屈服方程和TL平面屈服方程描述的屈服面上的点为应力-应变曲线非线性阶段起始点或线性阶段终点，因此用材受剪状态下的应力水平达到屈服面时，用材出现断裂，断裂时的应变增量对应的应力变为零，同样采用公式(12)进行描述，同样将试探应力代入屈服方程，屈服方程根据变形出现的平面拟定，分别为LR平面屈服方程、RT平面屈服方程和TL平面屈服方程，若
Figure FDA0002773316730000061
则说明用材仍处于弹性阶段，令
Figure FDA0002773316730000062
作为新的应力张量，按照公式(12)进行循环，若
Figure FDA0002773316730000063
则说明用材已经进入屈服阶段：用材处于受剪状态时发生断裂，保留该应变增量，应变增量步后对应的应力为零；
(5)用材LR平面受压：用材LR平面受压，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定；用材LR平面受压时的应力状态达到公式(3)、公式(4)和公式(5)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：
假设第n步迭代得到的应力张量为
Figure FDA0002773316730000064
n+1步中给出一个新的应变增量dεij，根据公式(12)计算得到试探应力
Figure FDA0002773316730000065
若
Figure FDA0002773316730000066
在LR平面中，kn'＝k1，则说明用材已经进入屈服阶段，该试探应力不是真实用材应力，需要将超出屈服面的部分减去：将总应变增量中的塑性应变增量引起的虚假应力从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，具体方程如下：
总应变增量dεij由弹性应变增量
Figure FDA0002773316730000071
和塑性应变增量
Figure FDA0002773316730000072
之和进行表达，
Figure FDA0002773316730000073
若
Figure FDA0002773316730000074
n+1步的真实应力为：
Figure FDA0002773316730000075
其中
Figure FDA0002773316730000076
为塑性应变增量，公式(14)能够使用材应力始终处于屈服面上，实现理想塑性流动状态，其中塑性应变增量
Figure FDA0002773316730000077
根据相关联的流动法则求解：
Figure FDA0002773316730000078
其中△λ为一致性系数，根据公式(15)，LR平面上的塑性应变表示为：
Figure FDA0002773316730000079
其中
Figure FDA00027733167300000710
为用材沿L轴方向的塑性正应变增量，
Figure FDA00027733167300000711
为用材位于LR平面内的塑性剪应变增量，
Figure FDA00027733167300000712
为用材沿R轴方向的塑性正应变增量，△λLR求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：
Figure FDA00027733167300000713
其中
Figure FDA00027733167300000714
为第n步试探应力代入公式(3)的结果，而
Figure FDA00027733167300000715
为第n+1步试探应力代入公式(3)的结果；
故此
Figure FDA00027733167300000716
其中
Figure FDA0002773316730000081
为
Figure FDA0002773316730000082
与
Figure FDA0002773316730000083
的差值；
将上式按照LR平面屈服方程进行一阶泰勒展开，得到：
Figure FDA0002773316730000084
其中△σ11为沿着用材L轴方向的正应力增量，△σ22为沿着用材R轴方向的正应力增量，△σ12为位于用材LR平面内的剪应力增量；
Figure FDA0002773316730000085
其中M、X、Y、Z为运算过程中的过渡量；根据径向返回算法，公式(19)中的应力增量表述为：
Figure FDA0002773316730000086
其中
Figure FDA0002773316730000087
为沿着用材L轴方向的试探正应力增量，
Figure FDA0002773316730000088
为沿着用材R轴方向的试探正应力增量，
Figure FDA0002773316730000089
为位于用材LR平面内的试探剪应力增量；
将公式(20)和公式(21)代入公式(19)，△λLR被求解：
Figure FDA0002773316730000091
在总应变达到应力软化阈值时，应力软化阈值
Figure FDA0002773316730000092
为弹性总应变阈值
Figure FDA0002773316730000093
的1.2倍，出现应力软化现象；
(6)用材RT平面受压：用材RT平面受压，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定，用材RT平面受压时的应力状态达到公式(6)、公式(7)和公式(8)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：
假设第n步迭代得到的应力张量为
Figure FDA0002773316730000094
n+1步中给出一个新的应变增量dξij，根据公式(12)计算得到试探应力
Figure FDA0002773316730000095
若
Figure FDA0002773316730000096
则说明用材已经进入屈服阶段，在RT平面中，kn'＝k2，该试探应力不是真实用材应力，需要将超出屈服面的部分减去：将总应变增量中的塑性应变引起的虚假应力增量从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，总应变增量由弹性应变增量和塑性应变增量之和进行表达，见公式(13)，真实应力计算见公式(14)，相关联塑性流动法则见公式(15)，RT平面内的塑性应变表示为：
Figure FDA0002773316730000097
其中
Figure FDA0002773316730000098
用材沿R轴方向的塑性正应变增量，
Figure FDA0002773316730000099
为用材位于RT平面内的塑性剪应变增量，
Figure FDA00027733167300000910
为用材沿T轴方向的塑性正应变增量，△λRT的求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：
Figure FDA00027733167300000911
其中
Figure FDA0002773316730000101
为第n步试探应力代入公式(6)的结果，而
Figure FDA0002773316730000102
为第n+1步试探应力代入公式(6)的结果；故此：
Figure FDA0002773316730000103
其中
Figure FDA0002773316730000104
为
Figure FDA0002773316730000105
与
Figure FDA0002773316730000106
的差值，将上式按照LR平面屈服方程进行一阶泰勒展开，得到：
Figure FDA0002773316730000107
其中△△22为沿着用材R轴方向的正应力增量，△△33为沿着用材T轴方向的正应力增量，△△23为位于用材RT平面内的剪应力增量；
Figure FDA0002773316730000108
其中U、V、W为运算过程中的过渡量，公式(29)中的应力增量表述为：
Figure FDA0002773316730000109
其中
Figure FDA00027733167300001010
为沿着用材R轴方向的试探正应力增量，
Figure FDA00027733167300001011
为沿着用材T轴方向的试探正应力增量，
Figure FDA00027733167300001012
为位于用材RT平面内的试探剪应力增量；
将公式(30)和公式(31)代入公式(29)，△λRT被求解：
Figure FDA0002773316730000111
在总应变达到应力硬化阈值时，应力硬化阈值
Figure FDA0002773316730000112
为弹性总应变阈值
Figure FDA0002773316730000113
的2倍，用材出现应力硬化现象，木材领域统称为二次强化现象；
(7)用材TL平面受压：用材TL平面受压，简化为3个阶段，弹性阶段由公式(1)和公式(2)确定；用材RT平面受压时的应力状态达到公式(9)、公式(10)和公式(11)确定的屈服面后，用材进入理想塑性流动阶段，过程如下：
假设第n步迭代得到的应力张量为
Figure FDA0002773316730000114
n+1步中给出一个新的应变增量dβij，根据公式(12)计算得到试探应力
Figure FDA0002773316730000115
若
Figure FDA0002773316730000116
在TL平面中，kn'＝k3，则说明用材已经进入屈服阶段，该试探应力不是真实用材应力，需要将超出屈服面的部分减去：将总应变增量中的塑性应变引起的虚假应力增量从试探应力中减去，从而使应力值返回屈服面已实现理想塑性流动状态，总应变增量由弹性应变增量和塑性应变增量之和进行表达，见公式(13)，真实应力计算见公式(14)，相关联塑性流动法则见公式(15)，TL平面内的塑性应变表示为：
Figure FDA0002773316730000117
其中
Figure FDA0002773316730000118
为用材沿T轴方向的塑性正应变增量，
Figure FDA0002773316730000119
为用材位于TL平面内的塑性剪应变增量，
Figure FDA00027733167300001110
为用材沿L轴方向的塑性正应变增量，△λTL求解是上述过程的关键，假设用材在第n步屈服，那么用材在第n+1步同样屈服，应力状态均在屈服面上，即：
Figure FDA0002773316730000121
其中
Figure FDA0002773316730000122
为第n步试探应力代入公式(9)的结果，而
Figure FDA0002773316730000123
为第n+1步试探应力代入公式(9)的结果；
故此：
Figure FDA0002773316730000124
其中
Figure FDA0002773316730000125
为
Figure FDA0002773316730000126
与
Figure FDA0002773316730000127
的差值；
将上式按照TL平面屈服方程进行一阶泰勒展开，得到：
Figure FDA0002773316730000128
其中△α11为沿着用材L轴方向的正应力增量，△α33为沿着用材T轴方向的正应力增量，△α13为位于用材TL平面内的剪应力增量；
其中：
Figure FDA0002773316730000129
其中L、G、H、I为运算过程中的过渡量，公式(40)中的应力增量表述为：
Figure FDA0002773316730000131
其中
Figure FDA0002773316730000132
为沿着用材T轴方向的试探正应力增量，
Figure FDA0002773316730000133
为沿着用材L轴方向的试探正应力增量，
Figure FDA0002773316730000134
为位于用材TL平面内的试探剪应力增量；
将公式(41)和公式(42)代入公式(40)，△λTL被求解
Figure FDA0002773316730000135
在总应变达到应力硬化阈值时，应力硬化阈值应变
Figure FDA0002773316730000136
为弹性总应变
Figure FDA0002773316730000137
的2倍，应力出现硬化，木材领域统称为二次强化现象。
2.如权利要求1所述的适于中国古建筑木结构用材本构关系的数值模拟方法，其特征在于，步骤(5)中，引入损伤因子描述应力软化现象，采用应力逐步退化模型模拟用材应力软化：
Figure FDA0002773316730000138
其中，σij为无损应力，dLR为损伤因子，当用材没有损伤时，dLR＝0，当木材完全损伤时，dLR＝1，
Figure FDA0002773316730000139
为损伤后的应力；损伤因子与当前状态应变、应力和能量函数有关，此处采用一种与总应变相关的损伤累计理论模型，该模型认为材料损伤与总应变和无损弹性模量相关，具体表示如下：
Figure FDA00027733167300001310
其中τ为总应变对应的应变能，εij、εkl为LR平面对应的应变；
用材LR平面内应力软化对应的损伤因子如下：
dLR＝1-exp[-(τLR-τLR,0)/τLR,0] (25)
其中，τLR为LR平面内总应变对应的应变能，而τLR,0为应力软化阈值对应的应变能。
3.如权利要求1所述的适于中国古建筑木结构用材本构关系的数值模拟方法，其特征在于，步骤(6)中采用设置最初屈服面和最终屈服面的方法来描述二次强化现象，最初屈服面为用材进入塑性流动阶段的起点，对应的最初的屈服面方程，见公式(6)；引入NRT表示用材RT平面受压时，最终屈服强度和初始屈服强度之比，故此最终屈服面公式如下：
Figure FDA0002773316730000141
其中，NRT-2为RT平面受压时，R轴的最终屈服强度和初始屈服强度之比，数值通过试验曲线确定，NRT-3为RT平面受力时，T轴的最终屈服强度和初始屈服强度之比，数值通过试验曲线确定，k2是与用材应力有关的参数，通过公式(6)确定，屈服面转移通过背应力增量进行描述，背应力增长与应变的增量相关，总背应力增量是最终屈服应力和初始屈服应力之差，RT平面背应力构造方式如下：
Figure FDA0002773316730000142
其中，cRT为RT平面内背应力增长的快慢程度，由RT平面主轴受压试验获取，GRT-2和GRT-3为RT平面二次硬化屈服面转移约束方程，其中ω22为R轴方向的背应力，ω33为T轴方向的背应力，△ω22为R轴方向的背应力增量，△ω33为T轴方向的背应力增量，防止增加的屈服面超过最终屈服面；
Figure FDA0002773316730000151
通过公式(35)可知，最初GRT-2和GRT-3均为1，对背应力增长无限制作用，当背应力总量达到最大值时，背应力停止增长，公式(35)均为0，而
Figure FDA0002773316730000152
和
Figure FDA0002773316730000153
表征背应力的增长方向，其中
Figure FDA0002773316730000154
和
Figure FDA0002773316730000155
为用材最终屈服强度，用试验获取，
Figure FDA0002773316730000156
为等效塑性应变增量，由如下公式获取；
Figure FDA0002773316730000157
4.如权利要求1所述的适于中国古建筑木结构用材本构关系的数值模拟方法，其特征在于，步骤(7)中采用设置最初屈服面和最终屈服面的方法来描述二次强化现象，最初屈服面为用材进入理想塑性流动阶段的起点，对应的最初屈服面方程，见公式(9)；引入NTL表示用材TL平面受压时，最终屈服强度和初始屈服强度之比，故此最终屈服面公式如下：
Figure FDA0002773316730000158
其中，NTL-3为TL平面受力时，T轴的最终屈服强度和初始强度之比，数值通过试验曲线确定；NTL-1为TL平面受力时，L轴的最终屈服强度和初始强度之比，数值通过试验曲线确定；k3是与用材应力有关的参数，通过公式(9)确定，屈服面转移通过背应力增量进行描述，背应力增长与应变的增量相关，总背应力增量是最终屈服应力和初始屈服应力之差，TL平面背应力构造方式如下：
Figure FDA0002773316730000161
其中，cTL为TL平面内背应力增长的快慢程度，由TL平面主轴受压试验获取；GTL-3和GTL-1为TL平面二次硬化屈服面转移约束方程，其中ρ11为L轴方向的背应力，ρ33为T轴方向的背应力，△ρ11为L轴方向的背应力增量，△ρ33为T轴方向的背应力增量，防止增加的屈服面超过最终屈服面；
Figure FDA0002773316730000162
通过公式(46)可知，最初GTL-3和GTL-1均为1，对背应力增长无限制作用，当背应力总量达到最大值时，背应力停止增长，公式(46)均为0，而
Figure FDA0002773316730000163
和
Figure FDA0002773316730000164
表征背应力的增长方向，其中
Figure FDA0002773316730000165
和
Figure FDA0002773316730000166
为用材最终屈服强度，用试验获取，
Figure FDA0002773316730000167
为等效塑性应变增量，具体表示如下所示：
Figure FDA0002773316730000171
CN202011256618.6A 2020-11-11 2020-11-11 一种适于中国古建筑木结构用材本构关系的数值模拟方法 Active CN112541285B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011256618.6A CN112541285B (zh) 2020-11-11 2020-11-11 一种适于中国古建筑木结构用材本构关系的数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011256618.6A CN112541285B (zh) 2020-11-11 2020-11-11 一种适于中国古建筑木结构用材本构关系的数值模拟方法

Publications (2)

Publication Number Publication Date
CN112541285A true CN112541285A (zh) 2021-03-23
CN112541285B CN112541285B (zh) 2024-03-12

Family

ID=75015061

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011256618.6A Active CN112541285B (zh) 2020-11-11 2020-11-11 一种适于中国古建筑木结构用材本构关系的数值模拟方法

Country Status (1)

Country Link
CN (1) CN112541285B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113761729A (zh) * 2021-08-25 2021-12-07 中国林业科学研究院木材工业研究所 基于木材弱相结构的木材横纹承压本构关系模型构建方法、装置及存储介质
CN114112676A (zh) * 2021-12-03 2022-03-01 中国林业科学研究院木材工业研究所 一种木材横纹抗压全时程本构关系的构建方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003194686A (ja) * 2001-12-27 2003-07-09 Toyota Motor Corp 応力−ひずみ関係シミュレート方法および除荷過程における降伏点を求める方法
CN102364489A (zh) * 2011-10-25 2012-02-29 陈志勇 木材复杂各向异性本构关系模型的数值模拟方法
CN106227928A (zh) * 2016-07-20 2016-12-14 福州大学 木材各向异性塑性屈服本构模型的数值模拟方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003194686A (ja) * 2001-12-27 2003-07-09 Toyota Motor Corp 応力−ひずみ関係シミュレート方法および除荷過程における降伏点を求める方法
CN102364489A (zh) * 2011-10-25 2012-02-29 陈志勇 木材复杂各向异性本构关系模型的数值模拟方法
CN106227928A (zh) * 2016-07-20 2016-12-14 福州大学 木材各向异性塑性屈服本构模型的数值模拟方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
杨娜;张雷;秦术杰;: "一种描述木材受压的非线性本构模型及试验验证", 土木工程学报, no. 04, pages 84 - 92 *
陈志勇;祝恩淳;潘景龙;: "复杂应力状态下木材力学性能的数值模拟", 计算力学学报, no. 04, pages 142 - 147 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113761729A (zh) * 2021-08-25 2021-12-07 中国林业科学研究院木材工业研究所 基于木材弱相结构的木材横纹承压本构关系模型构建方法、装置及存储介质
CN113761729B (zh) * 2021-08-25 2024-01-02 中国林业科学研究院木材工业研究所 基于木材弱相结构的木材横纹承压本构关系模型构建方法、装置及存储介质
CN114112676A (zh) * 2021-12-03 2022-03-01 中国林业科学研究院木材工业研究所 一种木材横纹抗压全时程本构关系的构建方法
CN114112676B (zh) * 2021-12-03 2023-12-19 中国林业科学研究院木材工业研究所 一种木材横纹抗压全时程本构关系的构建方法

Also Published As

Publication number Publication date
CN112541285B (zh) 2024-03-12

Similar Documents

Publication Publication Date Title
Mazzotti et al. Nonlinear creep damage model for concrete under uniaxial compression
Barros et al. Near surface mounted CFRP strips for the flexural strengthening of RC columns: Experimental and numerical research
CN112541285A (zh) 一种适于中国古建筑木结构用材本构关系的数值模拟方法
CN103942441A (zh) 基于应力比影响的碳纤维复合材料疲劳寿命评估方法
Wang et al. Micromechanical modelling of the progressive failure in unidirectional composites reinforced with bamboo fibres
CN114065569A (zh) 一种确定玄武岩筋废旧钢纤维混凝土梁极限受弯承载力的方法
Lina et al. FEM analysis of spring-backs in age forming of aluminum alloy plates
Ju et al. Fixed-angle smeared-truss approach with direct tension force transfer model for torsional behavior of steel fiber-reinforced concrete members
Zhu et al. Numerical study on the influence of mesomechanical properties on macroscopic fracture of concrete
Khabaz Experimental and Numerical Investigation of Single Fiber Pull-Out Tests of Steel Macro-Fiber and Glass Micro-Fiber in a Cementitious Matrix
CN112287474A (zh) 用于模拟聚乙烯醇纤维混凝土力学行为的算法、设备及存储介质
Wen et al. Experimental study on the bond degradation of basalt fiber reinforced polymer grid-concrete interface under fatigue loading
Madhusoodanan et al. A physically based fatigue damage model for simulating three-dimensional stress states in composites under very high cycle fatigue loading
Khelifa et al. Fatigue study of E-glass fiber reinforced polyester composite under fully reversed loading and spectrum loading
CN114357564A (zh) 一种岩土材料本构模型建立方法
Diao et al. Simulation of fatigue performance of cross-ply composite laminates
Vienni et al. CRM reinforced brick masonry walls: Experimental and parametric numerical investigations
CN112710566A (zh) 一种界面ii型裂纹的临界能量释放率测试方法
Yang et al. Double pull specimen more suitable for measuring bond-slip relationship of FRP-to-concrete interface
Taqi et al. Numerical Analysis of Corrosion Reinforcements in Fibrous Concrete Beams
Hosseini-Toudeshky et al. Implementation of a micro-meso approach for progressive damage analysis of composite laminates
Dick-Nielsen et al. Simulation of strain-hardening in ECC uniaxial test specimen by use of a damage mechanics formulation
Moldovan et al. Analysing the Performance of Hollow Core Slabs Strengthened with CFRP
Thomson et al. Numerical prediction of the ballistic performance of hygrothermally aged CFRP laminates using a multi-scale modelling approach
Joffe et al. Damage evolution modeling in multidirectional laminates and the resulting nonlinear response

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant