CN111709176B - 考虑塑性和损伤本构关系的有限元模拟方法及系统 - Google Patents
考虑塑性和损伤本构关系的有限元模拟方法及系统 Download PDFInfo
- Publication number
- CN111709176B CN111709176B CN202010732552.7A CN202010732552A CN111709176B CN 111709176 B CN111709176 B CN 111709176B CN 202010732552 A CN202010732552 A CN 202010732552A CN 111709176 B CN111709176 B CN 111709176B
- Authority
- CN
- China
- Prior art keywords
- stress
- damage
- representing
- tensor
- elastic
- 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.)
- Active
Links
- 230000006378 damage Effects 0.000 title claims abstract description 305
- 238000000034 method Methods 0.000 title claims abstract description 40
- 238000004088 simulation Methods 0.000 title claims abstract description 27
- 239000002184 metal Substances 0.000 claims abstract description 29
- 229910052751 metal Inorganic materials 0.000 claims abstract description 29
- 238000003754 machining Methods 0.000 claims abstract description 20
- 238000004458 analytical method Methods 0.000 claims abstract description 12
- 238000005381 potential energy Methods 0.000 claims description 77
- 239000000463 material Substances 0.000 claims description 76
- 208000027418 Wounds and injury Diseases 0.000 claims description 51
- 208000014674 injury Diseases 0.000 claims description 51
- 230000006399 behavior Effects 0.000 claims description 33
- 230000000694 effects Effects 0.000 claims description 28
- 238000004364 calculation method Methods 0.000 claims description 23
- 230000006870 function Effects 0.000 claims description 22
- 230000003902 lesion Effects 0.000 claims description 22
- 230000008878 coupling Effects 0.000 claims description 18
- 238000010168 coupling process Methods 0.000 claims description 18
- 238000005859 coupling reaction Methods 0.000 claims description 18
- 230000001186 cumulative effect Effects 0.000 claims description 14
- 230000006835 compression Effects 0.000 claims description 13
- 238000007906 compression Methods 0.000 claims description 13
- 230000006735 deficit Effects 0.000 claims description 13
- 230000004069 differentiation Effects 0.000 claims description 8
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 7
- 230000008859 change Effects 0.000 claims description 7
- 238000009795 derivation Methods 0.000 claims description 7
- 230000002706 hydrostatic effect Effects 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000002844 melting Methods 0.000 claims description 7
- 230000008018 melting Effects 0.000 claims description 7
- 230000009467 reduction Effects 0.000 claims description 7
- 238000003825 pressing Methods 0.000 claims description 6
- 238000010008 shearing Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 abstract description 12
- 238000012545 processing Methods 0.000 abstract description 10
- 238000013461 design Methods 0.000 description 3
- 102100021503 ATP-binding cassette sub-family B member 6 Human genes 0.000 description 2
- 101100000375 Homo sapiens ABCB6 gene Proteins 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000005489 elastic deformation Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000005336 cracking Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 238000005242 forging Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 150000002739 metals Chemical class 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
- 238000005482 strain hardening Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
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)
- Liquid Crystal (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
本发明提供了一种考虑塑性和损伤本构关系的有限元模拟方法及系统，包括：利用连续损伤力学和Johnson‑Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为。本发明通过应力三轴度和Lode系数表征的应变路径，将塑性损伤和Johnson‑Cook模型进行耦合，在不同温度和加载速度的环境下，对金属成形和加工过程中塑性力学行为和损伤演化进行仿真分析。
Description
技术领域
本发明涉及材料加工技术领域，具体地，涉及考虑塑性和损伤本构关系的有限元模拟方法及系统。尤其地，涉及一种考虑塑性和损伤本构关系的有限元模拟方法。
背景技术
金属成形和加工过程的有限元模拟对工业优化和设计高质量产品具有重大影响，它可以利用计算机虚拟地研究加工工艺参数、成形/加工工具，甚至材料力学行为。从而大量地降低生产成本，有效地缩短设计周期。目前，它已被应用于应力、应变、电、传热和流体流动等各个领域。它是一种求偏微分方程及其方程组近似解的数值技术，也可以求积分方程的近似解。简而言之，有限元分析是把一个非常复杂的问题分解成可以解决的小单元的一种方法。有限元方法在金属成形和加工过程的建模中是一种非常有效的工具，因为它为产品、工具、机器和工艺的设计提供了详细的信息。
然而，如何准确地预测材料的行为，特别是塑性损伤，对工程技术人员来说是一个巨大的挑战。建立合适的损伤演化模型，有利于避免金属成形过程中工件的损伤产生，有利于金属加工过程中工件和废料的分离。另一个挑战来自于金属成形和加工过程中复杂的应变路径和环境条件。例如，在不同的环境温度或加载速率下，耦合损伤行为的弹塑性响应变化也不同。同时，在不同的应变路径下损伤演化具有多样性。因此，准确描述损伤演化过程是提高金属成形和加工过程精确度的关键。
发明内容
针对现有技术中的缺陷，本发明的目的是提供一种考虑塑性和损伤本构关系的有限元模拟方法及系统。
根据本发明提供的一种考虑塑性和损伤本构关系的有限元模拟方法，包括：
利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为；
所述Johnson-Cook塑性流动准则：
其中，
表示应力；
表示应变率；
A表示屈服应力；
B表示硬化模量；
表示塑性应变；
C表示黏性参数；
ln表示数学运算符号，表示无量纲塑性应变率；
T表示当前温度；
T0表示室温；
Tm表示熔化温度；
n表示硬化指数；
m表示温度软化参数。
优选地，所述利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为，包括：
步骤S1：损伤的宏观表述；
步骤S2：耦合损伤与弹塑性；
步骤S3：耦合损伤与Johnson-Cook硬化模型。
优选地，所述步骤S1：
以损伤因子D表示材料实际承载面积的减少量，当材料处于初始状态，无损伤时，D＝0；当D＝1时表示材料完全破裂损伤；
基于等效应变假设其中/>将材料损伤中微孔洞的作用描述为材料弹性模量的变化；
通过引入单轴损伤演化参数h(0≤h≤1)来描述微裂纹在压缩加载状态下闭合效应，h＝1表示在没有裂纹闭合效应的情况下，拉伸状态下的损伤演化；h＝0表示裂纹闭合，纯压缩状态下无损伤演化；
根据弹性能量等效理论，在主坐标系下的各向同性损伤等效变量如式(2)：
其中，
表示损伤柯西应力张量；
D表示损伤变量；
σ表示柯西应力张量；
h表示单边应力状态参数；
表示损伤弹性张量；
ε e表示弹性张量；
We表示弹性势能；
E表示弹性模量张量；
表示损伤弹性模量张量。
优选地，所述步骤S2：
基于Helmholtz自由能方程，考虑弹性应变张量ε e和损伤变量D，弹性损伤势能方程可以表示为式(3)：
其中，
ρΨ(ε e,D)表示弹性损伤势能；
Ψ表示自由能；
tr表示运算函数；
λe,μe均为拉梅系数；
ρ表示材料密度；
Lame系数(λe,μe)可以用杨氏模量E和泊松比ν表示：
在该势能方程中，状态关系可以表示为：
其中，
ρ表示材料密度
表示自由能ψ关于弹性张量ε e的微分；
表示自由能ψ关于损伤D的微分；
Y(ε e,D)表示弹性损伤释放率；
自由能的微分形式可表示为式(8)：
其中，
为弹性损伤势能的微分形式；
为弹性张量的微分形式；
为损伤变量的微分形式。
优选地，所述步骤S3：
为了耦合损伤行为，控制应力空间内部变量的演化规律，定义势能方程如下：
F(σ,Y；D)＝fp+FY (9)
其中，
F(σ,Y；D)表示势能方程；
塑性势能fp和损伤势能FY计算方式如式(10)：
其中，
为应力张量第二不变量；
参数Y0,α,β和γ控制损伤势能演化；
σy为屈服应力，通过Johnson-Cook塑性流动准则进行求解；
根据增量理论，温度增量和应变速率/>可由式(11)进行计算：
其中，
η为自定义非弹性热系数；
为塑性应变率；
ρ为材料密度；
c为比热容；
应变速率为时间间隔Δt内的平均塑性应变；
根据上述应变速率的定义，Johnson-Cook硬化模型的模块可表示为:
其中，
σy表示屈服应力；
表示塑性应变；
表示初始应变率；
Johnson-Cook硬化模型中等效塑性应变相关变量的求导为：
其中，
η表示非弹性热系数；
R表示应力；
ρ表示材料密度；
c表示比热容；
对于特定的应力张量σ，选择的坐标轴不同，应力张量的矩阵形式也会有所不同。但是存在三个应力张量不变量：
其中，
J1表示应力偏张量第一不变量；
J2表示应力偏张量第二不变量；
J3表示应力偏张量第三不变量；
σH表示平均应力；
σ表示柯西应力张量；
σ1表示第一主应力；
σ2表示第二主应力；
σ3表示第三主应力；
σeq表示米塞斯等效应力；
det表示一个函数运算，类似于加减乘除，在数学计算中常用的一个函数；
S表示应力偏张量；
所以应力状态可用三个应力张量不变量来表示，为了简化，引入Lode角θ和应力三轴度ξ，在主应力空间中，通过坐标原点的偏应力平面为π平面；
Lode角表示应力张量与最大主应力轴夹角的投影，应力三轴度表示静水应力与等效应力的比值
应力状态用主应力空间中(σm,σeq,θ)来表示；
Lode角θ往往用Lode系数来表示使得应力三轴度ξ与Lode系数η均定义在[-1,1]内；
所有应力方向都可以通过上述定义的参数(ξ,θ)来表征；
为了考虑Lode角θ和应力三轴度ξ对损伤演化的影响，引入了参数h(ξ,θ)(0≤h(ξ,θ)≤1)；
h(ξ1,θ1)＝1表示在状态(ξ1,θ1)时拥有最快的损伤演化；h(ξ2,θ2)＝0表示在状态(ξ2,θ2)时的没有损伤演化；
应力状态参数h(ξ,θ)不仅可以考虑微裂纹闭合的效果，还可以区分剪切状态与拉压状态；
当用耦合了Lode角和应力三轴度的应力状态参数h(ξ,θ)替换了单轴损伤演化参数h后，一个更精确的本构模型被提出：
其中，
ρΨ(ε e,D)表示弹性损伤势能
等效应力张量和损伤释放率可以用式(16)进行更新：
其中，
h(ξ,θ)表示应力状态参数
Y(ε e,D)表示损伤释放率
通过耗散分析，累积等效塑性应变和损伤生长速率按式(17)更新：
其中，
表示累积等效塑性应变；
表示损伤生长速率；
J2(σ)表示应力偏张量第二不变量；
表示塑性系数；
Y0表示初始损伤参数；
β表示损伤相关系数。
应力状态参数明显耦合为塑性应变和损伤演化，损伤演化与应力状态参数有密切的关系，应力状态参数越高，损伤演化越快，低应力状态参数也导致损伤演化较慢，特别是当h(ξ,θ)＝0时，损伤没有演化。
根据本发明提供的一种考虑塑性和损伤本构关系的有限元模拟系统，包括：
利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为；
所述Johnson-Cook塑性流动准则：
其中，
表示应力；
表示应变率；
A表示屈服应力；
B表示硬化模量；
表示塑性应变；
C表示黏性参数；
ln表示数学运算符号，表示无量纲塑性应变率；
T表示当前温度；
T0表示室温；
Tm表示熔化温度；
n表示硬化指数；
m表示温度软化参数。
优选地，所述利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为，包括：
模块S1：损伤的宏观表述；
模块S2：耦合损伤与弹塑性；
模块S3：耦合损伤与Johnson-Cook硬化模型。
优选地，所述模块S1：
以损伤因子D表示材料实际承载面积的减少量，当材料处于初始状态，无损伤时，D＝0；当D＝1时表示材料完全破裂损伤；
基于等效应变假设其中/>将材料损伤中微孔洞的作用描述为材料弹性模量的变化；
通过引入单轴损伤演化参数h(0≤h≤1)来描述微裂纹在压缩加载状态下闭合效应，h＝1表示在没有裂纹闭合效应的情况下，拉伸状态下的损伤演化；h＝0表示裂纹闭合，纯压缩状态下无损伤演化；
根据弹性能量等效理论，在主坐标系下的各向同性损伤等效变量如式(2)：
其中，
表示损伤柯西应力张量；
D表示损伤变量；
σ表示柯西应力张量；
h表示单边应力状态参数；
表示损伤弹性张量；
ε e表示弹性张量；
We表示弹性势能；
E表示弹性模量张量；
表示损伤弹性模量张量。
优选地，所述模块S2：
基于Helmholtz自由能方程，考虑弹性应变张量ε e和损伤变量D，弹性损伤势能方程可以表示为式(3)：
其中，
ρΨ(ε e,D)表示弹性损伤势能；
Ψ表示自由能；
tr表示运算函数；
λe,μe均为拉梅系数；
ρ表示材料密度；
Lame系数(λe,μe)可以用杨氏模量E和泊松比ν表示：
在该势能方程中，状态关系可以表示为：
其中，
ρ表示材料密度
表示自由能ψ关于弹性张量ε e的微分；
表示自由能ψ关于损伤D的微分；
Y(ε e,D)表示弹性损伤释放率；
自由能的微分形式可表示为式(8)：
其中，
为弹性损伤势能的微分形式；
为弹性张量的微分形式；
为损伤变量的微分形式。
优选地，所述模块S3：
为了耦合损伤行为，控制应力空间内部变量的演化规律，定义势能方程如下：
F(σ,Y；D)＝fp+FY (9)
其中，
F(σ,Y；D)表示势能方程；
塑性势能fp和损伤势能FY计算方式如式(10)：
其中，
为应力张量第二不变量；
参数Y0,α,β和γ控制损伤势能演化；
σy为屈服应力，通过Johnson-Cook塑性流动准则进行求解；
根据增量理论，温度增量和应变速率/>可由式(11)进行计算：
其中，
η为自定义非弹性热系数；
为塑性应变率；
ρ为材料密度；
c为比热容；
应变速率为时间间隔Δt内的平均塑性应变；
根据上述应变速率的定义，Johnson-Cook硬化模型的模块可表示为:
其中，
σy表示屈服应力；
表示塑性应变；
表示初始应变率；
Johnson-Cook硬化模型中等效塑性应变相关变量的求导为：
其中，
η表示非弹性热系数；
R表示应力；
ρ表示材料密度；
c表示比热容；
对于特定的应力张量σ，选择的坐标轴不同，应力张量的矩阵形式也会有所不同。但是存在三个应力张量不变量：
其中，
J1表示应力偏张量第一不变量；
J2表示应力偏张量第二不变量；
J3表示应力偏张量第三不变量；
σH表示平均应力；
σ表示柯西应力张量；
σ1表示第一主应力；
σ2表示第二主应力；
σ3表示第三主应力；
σeq表示米塞斯等效应力；
det表示一个函数运算，类似于加减乘除，在数学计算中常用的一个函数；
S表示应力偏张量；
所以应力状态可用三个应力张量不变量来表示，为了简化，引入Lode角θ和应力三轴度ξ，在主应力空间中，通过坐标原点的偏应力平面为π平面；
Lode角表示应力张量与最大主应力轴夹角的投影，应力三轴度表示静水应力与等效应力的比值
应力状态用主应力空间中(σm,σeq,θ)来表示；
Lode角θ往往用Lode系数来表示使得应力三轴度ξ与Lode系数η均定义在[-1,1]内；
所有应力方向都可以通过上述定义的参数(ξ,θ)来表征；
为了考虑Lode角θ和应力三轴度ξ对损伤演化的影响，引入了参数h(ξ,θ)(0≤h(ξ,θ)≤1)；
h(ξ1,θ1)＝1表示在状态(ξ1,θ1)时拥有最快的损伤演化；h(ξ2,θ2)＝0表示在状态(ξ2,θ2)时的没有损伤演化；
应力状态参数h(ξ,θ)不仅可以考虑微裂纹闭合的效果，还可以区分剪切状态与拉压状态；
当用耦合了Lode角和应力三轴度的应力状态参数h(ξ,θ)替换了单轴损伤演化参数h后，一个更精确的本构模型被提出：
其中，
ρΨ(ε e,D)表示弹性损伤势能
等效应力张量和损伤释放率可以用式(16)进行更新：
其中，
h(ξ,θ)表示应力状态参数
Y(ε e,D)表示损伤释放率
通过耗散分析，累积等效塑性应变和损伤生长速率按式(17)更新：
其中，
表示累积等效塑性应变；
表示损伤生长速率；
J2(σ)表示应力偏张量第二不变量；
表示塑性系数；
Y0表示初始损伤参数；
β表示损伤相关系数。
应力状态参数明显耦合为塑性应变和损伤演化，损伤演化与应力状态参数有密切的关系，应力状态参数越高，损伤演化越快，低应力状态参数也导致损伤演化较慢，特别是当h(ξ,θ)＝0时，损伤没有演化。
与现有技术相比，本发明具有如下的有益效果：
本发明有限元法作为目前普遍应用的重要仿真手段，可以准确、高效的模拟分析材料在任意时刻和位置的变形行为。然而仿真结果的准确与否严重取决于仿真时所选用的材料本构模型。全耦合的损伤模型是目前广泛应用的，能够较为准确描述材料损伤行为的一种材料失效模型，强调材料损伤变量与其他状态变量的全耦合关系。同时，为预测不同应变状态下材料的塑性力学行为和损伤行为，应力三轴度和Lode系数是必须考虑的两个重要状态参数。本发明通过应力三轴度和Lode系数表征的应变路径，将塑性损伤和Johnson-Cook模型进行耦合。在不同温度和加载速度的环境下，对金属成形和加工过程中塑性力学行为和损伤演化进行仿真分析。
本发明还具有以下优点：
1.可以准确预测延性损伤的发生；
2.适用于多轴加载；
3.可以充分考虑应力三轴度和Lode角的影响；
4.可以对不同加载条件下的金属成形和加工过程进行仿真分析。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1为本发明提供的使用ABAQUS子程序仿真流程示意图。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
根据本发明提供的一种考虑塑性和损伤本构关系的有限元模拟方法，包括：
利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为；
所述Johnson-Cook塑性流动准则：
其中，
表示应力；
表示应变率；
A表示屈服应力；
B表示硬化模量；
表示塑性应变；
C表示黏性参数；
ln表示数学运算符号，表示无量纲塑性应变率；
T表示当前温度；
T0表示室温；
Tm表示熔化温度；
n表示硬化指数；
m表示温度软化参数。
具体地，所述利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为，包括：
步骤S1：损伤的宏观表述；
步骤S2：耦合损伤与弹塑性；
步骤S3：耦合损伤与Johnson-Cook硬化模型。
具体地，所述步骤S1：
以损伤因子D表示材料实际承载面积的减少量，当材料处于初始状态，无损伤时，D＝0；当D＝1时表示材料完全破裂损伤；
基于等效应变假设其中/>将材料损伤中微孔洞的作用描述为材料弹性模量的变化；
通过引入单轴损伤演化参数h(0≤h≤1)来描述微裂纹在压缩加载状态下闭合效应，h＝1表示在没有裂纹闭合效应的情况下，拉伸状态下的损伤演化；h＝0表示裂纹闭合，纯压缩状态下无损伤演化；
根据弹性能量等效理论，在主坐标系下的各向同性损伤等效变量如式(2)：
其中，
表示损伤柯西应力张量；
D表示损伤变量；
σ表示柯西应力张量；
h表示单边应力状态参数；
表示损伤弹性张量；
ε e表示弹性张量；
We表示弹性势能；
E表示弹性模量张量；
表示损伤弹性模量张量。
具体地，所述步骤S2：
基于Helmholtz自由能方程，考虑弹性应变张量ε e和损伤变量D，弹性损伤势能方程可以表示为式(3)：
其中，
ρΨ(ε e,D)表示弹性损伤势能；
Ψ表示自由能；
tr表示运算函数；
λe,μe均为拉梅系数；
ρ表示材料密度；
Lame系数(λe,μe)可以用杨氏模量E和泊松比ν表示：
在该势能方程中，状态关系可以表示为：
其中，
ρ表示材料密度
表示自由能ψ关于弹性张量ε e的微分；
表示自由能ψ关于损伤D的微分；
Y(ε e,D)表示弹性损伤释放率；
自由能的微分形式可表示为式(8)：
其中，
为弹性损伤势能的微分形式；
为弹性张量的微分形式；/>
为损伤变量的微分形式。
具体地，所述步骤S3：
为了耦合损伤行为，控制应力空间内部变量的演化规律，定义势能方程如下：
F(σ,Y；D)＝fp+FY (9)
其中，
F(σ,Y；D)表示势能方程；
塑性势能fp和损伤势能FY计算方式如式(10)：
其中，
为应力张量第二不变量；
参数Y0,α,β和γ控制损伤势能演化；
σy为屈服应力，通过Johnson-Cook塑性流动准则进行求解；
根据增量理论，温度增量和应变速率/>可由式(11)进行计算：
其中，
η为自定义非弹性热系数；
为塑性应变率；
ρ为材料密度；
c为比热容；
应变速率为时间间隔Δt内的平均塑性应变；
根据上述应变速率的定义，Johnson-Cook硬化模型的模块可表示为:
其中，
σy表示屈服应力；
表示塑性应变；
表示初始应变率；/>
Johnson-Cook硬化模型中等效塑性应变相关变量的求导为：
其中，
η表示非弹性热系数；
R表示应力；
ρ表示材料密度；
c表示比热容；
对于特定的应力张量σ，选择的坐标轴不同，应力张量的矩阵形式也会有所不同。但是存在三个应力张量不变量：
其中，
J1表示应力偏张量第一不变量；
J2表示应力偏张量第二不变量；
J3表示应力偏张量第三不变量；
σH表示平均应力；
σ表示柯西应力张量；
σ1表示第一主应力；
σ2表示第二主应力；
σ3表示第三主应力；
σeq表示米塞斯等效应力；
det表示一个函数运算，类似于加减乘除，在数学计算中常用的一个函数；
S表示应力偏张量；
所以应力状态可用三个应力张量不变量来表示，为了简化，引入Lode角θ和应力三轴度ξ，在主应力空间中，通过坐标原点的偏应力平面为π平面；
Lode角表示应力张量与最大主应力轴夹角的投影，应力三轴度表示静水应力与等效应力的比值
应力状态用主应力空间中(σm,σeq,θ)来表示；
Lode角θ往往用Lode系数来表示使得应力三轴度ξ与Lode系数η均定义在[-1,1]内；
所有应力方向都可以通过上述定义的参数(ξ,θ)来表征；
为了考虑Lode角θ和应力三轴度ξ对损伤演化的影响，引入了参数h(ξ,θ)(0≤h(ξ,θ)≤1)；
h(ξ1,θ1)＝1表示在状态(ξ1,θ1)时拥有最快的损伤演化；h(ξ2,θ2)＝0表示在状态(ξ2,θ2)时的没有损伤演化；
应力状态参数h(ξ,θ)不仅可以考虑微裂纹闭合的效果，还可以区分剪切状态与拉压状态；
当用耦合了Lode角和应力三轴度的应力状态参数h(ξ,θ)替换了单轴损伤演化参数h后，一个更精确的本构模型被提出：
其中，
ρΨ(ε e,D)表示弹性损伤势能
等效应力张量和损伤释放率可以用式(16)进行更新：
其中，
h(ξ,θ)表示应力状态参数
Y(ε e,D)表示损伤释放率
通过耗散分析，累积等效塑性应变和损伤生长速率按式(17)更新：
其中，
表示累积等效塑性应变；
表示损伤生长速率；
J2(σ)表示应力偏张量第二不变量；
表示塑性系数；
Y0表示初始损伤参数；
β表示损伤相关系数。
应力状态参数明显耦合为塑性应变和损伤演化，损伤演化与应力状态参数有密切的关系，应力状态参数越高，损伤演化越快，低应力状态参数也导致损伤演化较慢，特别是当h(ξ,θ)＝0时，损伤没有演化。
根据本发明提供的一种考虑塑性和损伤本构关系的有限元模拟系统，包括：
利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为；
所述Johnson-Cook塑性流动准则：
其中，
表示应力；
表示应变率；
A表示屈服应力；
B表示硬化模量；
表示塑性应变；
C表示黏性参数；
ln表示数学运算符号，表示无量纲塑性应变率；
T表示当前温度；
T0表示室温；
Tm表示熔化温度；
n表示硬化指数；
m表示温度软化参数。
具体地，所述利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为，包括：
模块S1：损伤的宏观表述；
模块S2：耦合损伤与弹塑性；
模块S3：耦合损伤与Johnson-Cook硬化模型。
具体地，所述模块S1：
以损伤因子D表示材料实际承载面积的减少量，当材料处于初始状态，无损伤时，D＝0；当D＝1时表示材料完全破裂损伤；
基于等效应变假设其中/>将材料损伤中微孔洞的作用描述为材料弹性模量的变化；
通过引入单轴损伤演化参数h(0≤h≤1)来描述微裂纹在压缩加载状态下闭合效应，h＝1表示在没有裂纹闭合效应的情况下，拉伸状态下的损伤演化；h＝0表示裂纹闭合，纯压缩状态下无损伤演化；
根据弹性能量等效理论，在主坐标系下的各向同性损伤等效变量如式(2)：
其中，
表示损伤柯西应力张量；
D表示损伤变量；
σ表示柯西应力张量；
h表示单边应力状态参数；
表示损伤弹性张量；
ε e表示弹性张量；
We表示弹性势能；
E表示弹性模量张量；
表示损伤弹性模量张量。
具体地，所述模块S2：
基于Helmholtz自由能方程，考虑弹性应变张量ε e和损伤变量D，弹性损伤势能方程可以表示为式(3)：
其中，
ρΨ(ε e,D)表示弹性损伤势能；
Ψ表示自由能；
tr表示运算函数；
λe,μe均为拉梅系数；
ρ表示材料密度；
Lame系数(λe,μe)可以用杨氏模量E和泊松比ν表示：
在该势能方程中，状态关系可以表示为：
其中，
ρ表示材料密度
表示自由能ψ关于弹性张量ε e的微分；
表示自由能ψ关于损伤D的微分；
Y(ε e,D)表示弹性损伤释放率；
自由能的微分形式可表示为式(8)：
其中，
为弹性损伤势能的微分形式；
为弹性张量的微分形式；
为损伤变量的微分形式。
具体地，所述模块S3：
为了耦合损伤行为，控制应力空间内部变量的演化规律，定义势能方程如下：
F(σ,Y；D)＝fp+FY (9)
其中，
F(σ,Y；D)表示势能方程；
塑性势能fp和损伤势能FY计算方式如式(10)：
其中，
为应力张量第二不变量；
参数Y0,α,β和γ控制损伤势能演化；
σy为屈服应力，通过Johnson-Cook塑性流动准则进行求解；
根据增量理论，温度增量和应变速率/>可由式(11)进行计算：
其中，
η为自定义非弹性热系数；
为塑性应变率；
ρ为材料密度；
c为比热容；
应变速率为时间间隔Δt内的平均塑性应变；/>
根据上述应变速率的定义，Johnson-Cook硬化模型的模块可表示为:
其中，
σy表示屈服应力；
表示塑性应变；
表示初始应变率；
Johnson-Cook硬化模型中等效塑性应变相关变量的求导为：
其中，
η表示非弹性热系数；
R表示应力；
ρ表示材料密度；
c表示比热容；
对于特定的应力张量σ，选择的坐标轴不同，应力张量的矩阵形式也会有所不同。但是存在三个应力张量不变量：
其中，
J1表示应力偏张量第一不变量；
J2表示应力偏张量第二不变量；
J3表示应力偏张量第三不变量；
σH表示平均应力；
σ表示柯西应力张量；
σ1表示第一主应力；
σ2表示第二主应力；
σ3表示第三主应力；
σeq表示米塞斯等效应力；
det表示一个函数运算，类似于加减乘除，在数学计算中常用的一个函数；
S表示应力偏张量；
所以应力状态可用三个应力张量不变量来表示，为了简化，引入Lode角θ和应力三轴度ξ，在主应力空间中，通过坐标原点的偏应力平面为π平面；
Lode角表示应力张量与最大主应力轴夹角的投影，应力三轴度表示静水应力与等效应力的比值
应力状态用主应力空间中(σm,σeq,θ)来表示；
Lode角θ往往用Lode系数来表示使得应力三轴度ξ与Lode系数η均定义在[-1,1]内；
所有应力方向都可以通过上述定义的参数(ξ,θ)来表征；
为了考虑Lode角θ和应力三轴度ξ对损伤演化的影响，引入了参数h(ξ,θ)(0≤h(ξ,θ)≤1)；
h(ξ1,θ1)＝1表示在状态(ξ1,θ1)时拥有最快的损伤演化；h(ξ2,θ2)＝0表示在状态(ξ2,θ2)时的没有损伤演化；
应力状态参数h(ξ,θ)不仅可以考虑微裂纹闭合的效果，还可以区分剪切状态与拉压状态；
当用耦合了Lode角和应力三轴度的应力状态参数h(ξ,θ)替换了单轴损伤演化参数h后，一个更精确的本构模型被提出：
其中，
ρΨ(ε e,D)表示弹性损伤势能
等效应力张量和损伤释放率可以用式(16)进行更新：
其中，
h(ξ,θ)表示应力状态参数
Y(ε e,D)表示损伤释放率
通过耗散分析，累积等效塑性应变和损伤生长速率按式(17)更新：
其中，
表示累积等效塑性应变；
表示损伤生长速率；
J2(σ)表示应力偏张量第二不变量；
表示塑性系数；
Y0表示初始损伤参数；
β表示损伤相关系数。
应力状态参数明显耦合为塑性应变和损伤演化，损伤演化与应力状态参数有密切的关系，应力状态参数越高，损伤演化越快，低应力状态参数也导致损伤演化较慢，特别是当h(ξ,θ)＝0时，损伤没有演化。
下面通过优选例，对本发明进行更为具体地说明。
优选例1：
本发明旨在准确地描述材料的变形行为，特别是塑性损伤演化行为，并考虑热软化和应变速率硬化的影响，提高金属成形和加工过程的仿真预测精度。
本发明建立了热-弹粘塑性与延性损伤的全耦合本构模型，对材料在大变形、温度、应变速率和延性损伤作用下，不同应变路径下的成形和加工过程进行仿真模拟。为了准确预测金属成形和加工模拟过程中材料行为，本发明在充分考虑应变硬化、应变速率和温度的影响基础之上，耦合材料的损伤演化行为。
Johnson-Cook模型作为一种在热力学领域广泛应用的经验型模型，已成功地用于模拟金属在不同温度和应变速率下的行为，但它不能解决由于微裂纹引起的材料刚度退化问题。因此，本发明的主要目的是利用连续损伤力学和Johnson-Cook塑性流动准则(式1)建立本构方程，以精确预测金属成形和加工过程中的损伤行为。
其中，
表示应力
表示应变率
A表示屈服应力
B表示硬化模量
表示塑性应变
C表示黏性参数
ln是一种数学运算符号，/>表示无量纲塑性应变率
T表示当前温度
T0表示室温
Tm表示熔化温度
n表示硬化指数
m表示温度软化参数
步骤1：损伤的宏观表述
本发明以损伤因子D表示材料实际承载面积的减少量，当材料处于初始状态，无损伤时，D＝0；当D＝1时表示材料完全破裂损伤。并基于等效应变假设其中将材料损伤中微孔洞的作用描述为材料弹性模量的变化。并通过引入单轴损伤演化参数h(0≤h≤1)来描述微裂纹在压缩加载状态下闭合效应，/>h＝1表示在没有裂纹闭合效应的情况下，拉伸状态下的损伤演化；h＝0表示裂纹闭合，纯压缩状态下无损伤演化。
根据弹性能量等效理论，在主坐标系下的各向同性损伤等效变量如式(2)：
其中，
表示损伤柯西应力张量
D表示损伤变量
σ表示柯西应力张量
h表示单边应力状态参数
表示损伤弹性张量
ε e表示弹性张量
We表示弹性势能
E表示弹性模量张量
表示损伤弹性模量张量
步骤2：耦合损伤与弹塑性
基于Helmholtz自由能方程，考虑弹性应变张量ε e和损伤变量D，弹性损伤势能方程可以表示为式(3)：
其中，
ρΨ(ε e,D)表示弹性损伤势能
Ψ表示自由能
tr是一种常用的运算函数
λe,μe均为拉梅系数
式中，ρ表示材料密度；Lame系数(λe,μe)可以用杨氏模量E和泊松比ν表示：
在该势能方程中，状态关系可以表示为：
其中，
ρ表示材料密度
在数学运算形式，与加减乘除一样，表示方程y关于x的微分求导，/>表示自由能ψ关于弹性张量ε e的微分；/>表示自由能ψ关于损伤D的微分；
Y(ε e,D)表示弹性损伤释放率
自由能的微分形式可表示为式(9)：
其中，
为弹性损伤势能的微分形式
为弹性张量的微分形式
为损伤变量的微分形式
步骤3：耦合损伤与Johnson-Cook硬化模型
为了耦合损伤行为，控制应力空间内部变量的演化规律，定义势能方程如下：
F(σ,Y；D)＝fp+FY (9)
其中，
F(σ,Y；D)表示势能方程
塑性势能fp和损伤势能FY计算方式如式(10)：
式中，为应力张量第二不变量；参数Y0,α,β和γ控制损伤势能演化；σy为屈服应力，通过Johnson-Cook塑性流动准则进行求解。
根据增量理论，温度增量和应变速率可由式(11)进行计算：
式中，η为自定义非弹性热系数，为塑性应变率，ρ为材料密度，c为比热容。应变率/>为时间间隔Δt内的平均塑性应变。根据上述绝热温度和应变速率的定义，Johnson-Cook硬化模型的模块可表示为:
其中，
σy表示屈服应力
表示塑性应变
表示初始应变率
Johnson-Cook硬化模型中等效塑性应变相关变量的求导为：
其中，
η表示非弹性热系数，R表示应力
ρ表示材料密度，c表示比热容
对于特定的应力张量σ，选择的坐标轴不同，应力张量的矩阵形式也会有所不同。但是存在三个应力张量不变量：
其中，
J1表示应力偏张量第一不变量
σH表示平均应力
σ表示柯西应力张量
σ1表示第一主应力
σ2表示第二主应力
σ3表示第三主应力
J2表示应力偏张量第二不变量
σeq表示米塞斯等效应力
J3表示应力偏张量第三不变量
det表示一个函数运算，类似于加减乘除，在数学计算中常用的一个函数
S表示应力偏张量
所以应力状态可用三个应力张量不变量来表示。为了简化，引入Lode角θ和应力三轴度ξ。在主应力空间中，通过坐标原点的偏应力平面为π平面。Lode角表示应力张量与最大主应力轴夹角的投影，应力三轴度表示静水应力与等效应力的比值所以应力状态可用主应力空间中(σm,σeq,θ)来表示。Lode角θ往往用Lode系数来表示/>使得应力三轴度ξ与Lode系数η均定义在[-1,1]内。所有应力方向(或称为加载条件或应力状态)都可以通过上述定义的参数(ξ,θ)来表征。
为了考虑Lode角θ和应力三轴度ξ对损伤演化的影响，引入了参数h(ξ,θ)(0≤h(ξ,θ)≤1)。h(ξ1,θ1)＝1表示在状态(ξ1,θ1)时拥有最快的损伤演化；h(ξ2,θ2)＝0表示在状态(ξ2,θ2)时的没有损伤演化。应力状态参数h(ξ,θ)不仅可以考虑微裂纹闭合的效果，还可以区分剪切状态与拉压状态。
因此，当用耦合了Lode角和应力三轴度的应力状态参数h(ξ,θ)替换了单轴损伤演化参数h后，一个更精确的本构模型被提出：
其中，
ρΨ(ε e,D)表示弹性损伤势能
等效应力张量和损伤释放率可以用式(16)进行更新：
/>
其中，
h(ξ,θ)表示应力状态参数
Y(ε e,D)表示损伤释放率
通过耗散分析，累积等效塑性应变和损伤生长速率按式(17)更新：
其中，
表示累积等效塑性应变；
表示损伤生长速率；
J2(σ)表示应力偏张量第二不变量
表示塑性系数
Y0表示初始损伤参数
β表示损伤相关系数
从这两个方程可以看出，应力状态参数明显耦合为塑性应变和损伤演化。可见，损伤演化与应力状态参数有密切的关系。应力状态参数越高，损伤演化越快。低应力状态参数也导致损伤演化较慢，特别是当h(ξ,θ)＝0时，损伤没有演化。
该方法可作为研究金属成形和加工过程损伤演化的有效工具，同时考虑了应力强度、应力三轴度和Lode角的影响。在金属成形过程中，例如在锻造和轧制中，工件压缩状态，损伤在负三轴加载下演化。在整个过程中，剪切损伤在损伤演化中占据主导地位。为了提高成形质量，必须避免在复杂的加载条件下产生的这种损伤。在复杂应力状态下，剪切带产生较大的剪切变形，最终产生裂纹。这些给准确预测复杂载荷条件下的损伤演化带来巨大的挑战，本发明通过定义应力状态参数h(ξ,θ)，使用耦合了应力三轴度和杆Lode角的损伤模型，对金属成形和加工过程进行有限元仿真分析。
优选例2：
结合附图对本发明进一步详细的说明，具体描述一个最佳实施例，这种描述的具体应使所属技术领域人员按照内容重现发明，而不必再创造性劳动。
本发明提出来一种考虑塑性和损伤本构关系的有限元模拟方法，下面将结合有限元仿真软件ABAQUS及其子程序UMAT对该发明进行说明。
用户子程序的主要作用就是对变量进行更新，在进行变量更新之前，他们的初始值首先被求解，如下：
式中，σJC为Johnson-Cook硬化准则，Y是损伤释放率。塑性流动法向张量用式(19)进行求解：
式(18)中的各变量以增量步Δt进行ti到ti+1分步求解，实现变量到/>的更新。本发明使用全隐式算法进行状态变量更新，如式(20)。使用牛顿迭代法求解塑性修正系数Δλ和损伤参数Di+1。
首先利用判定准则进行判断，若f＜0，则处于弹性变形阶段。在从ti到ti+1的时间增量内，发生的全是弹性变形，则
弹性应变张量可以表示为：
Lame系数(λe,μe)通过下式进行计算：
式中，E和ν分别表示杨氏模量和泊松比。其余变量可用式(24)进行更新。
式中，A,B,n表示硬化参数；C是粘性参数，m是温度软化参。应力三轴度和Lode角通过式(25)和(26)进行更新。
hi+1＝h(ξi+1,θi+1) (26)
塑性判定准则更新为：
如果塑性判定准则f＞0，则进入塑性变形阶段。该阶段的应变分为弹性应变和塑性应变两部分。因此，应变张量可表示为：
弹性应变和应力可由式(29)计算。
塑性流动方向表示为：
通过对式(20)和(30)观察，可以发现，所有非线性问题都集中于系数变量Δλ和Di+1上，转化为对下式的求解：
使用应力状态参数hi的损伤释放率Yi+1表示为：
利用一阶泰勒级数展开，得到每次迭代求解的方程如下:
M(Di+1,Δλ)和N(Di+1,Δλ)的偏导数定义如下：
式中，
对于的微分形式如下：
对于Johnson-Cook模型在ti+1时刻的变量求解如下：
偏微分的计算方法如下：
偏微分的计算方法如下：
/>
最终，对式(33)的计算结果如下：
结果δDi+1和δΔλ被用来对下一次迭代进行修正：
直至收敛，迭代结束，得到最终值Δλ和Dn+1，用于更新状态变量。
将上述过程，使用Fortan语言进行编译，并用于ABAQUS子程序的调用。使用有限元ABAQUS仿真软件及其子程序UMAT的流程如图1所示。
在本申请的描述中，需要理解的是，术语“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本申请和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本申请的限制。
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统、装置及其各个模块以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以，本发明提供的系统、装置及其各个模块可以被认为是一种硬件部件，而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构；也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。
Claims (6)
1.一种考虑塑性和损伤本构关系的有限元模拟方法，其特征在于，包括：
利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为；
所述Johnson-Cook塑性流动准则：
其中，
表示应力；
表示应变率；
A表示屈服应力；
B表示硬化模量；
表示塑性应变；
C表示黏性参数；
ln表示数学运算符号，表示无量纲塑性应变率；
T表示当前温度；
T0表示室温；
Tm表示熔化温度；
n表示硬化指数；
m表示温度软化参数；
所述利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为，包括：
步骤S1：损伤的宏观表述；
步骤S2：耦合损伤与弹塑性；
步骤S3：耦合损伤与Johnson-Cook硬化模型；
所述步骤S1：
以损伤因子D表示材料实际承载面积的减少量，当材料处于初始状态，无损伤时，D＝0；当D＝1时表示材料完全破裂损伤；
基于等效应变假设其中/>将材料损伤中微孔洞的作用描述为材料弹性模量的变化；
通过引入单轴损伤演化参数h(0≤h≤1)来描述微裂纹在压缩加载状态下闭合效应，h＝1表示在没有裂纹闭合效应的情况下，拉伸状态下的损伤演化；h＝0表示裂纹闭合，纯压缩状态下无损伤演化；
根据弹性能量等效理论，在主坐标系下的各向同性损伤等效变量如式(2)：
其中，
表示损伤柯西应力张量；
D表示损伤变量；
表示柯西应力张量；
表示单边应力状态参数；
表示损伤弹性张量；
表示弹性张量；
We表示弹性势能；
E表示弹性模量张量；
表示损伤弹性模量张量。
2.根据权利要求1所述的考虑塑性和损伤本构关系的有限元模拟方法，其特征在于，所述步骤S2：
基于Helmholtz自由能方程，考虑弹性应变张量εe和损伤变量D，弹性损伤势能方程可以表示为式(3)：
其中，
ρΨ(εe,)表示弹性损伤势能；
Ψ表示自由能；
tr表示运算函数；
λe,μe均为拉梅系数；
ρ表示材料密度；
Lame系数(λe,μe)可以用杨氏模量E和泊松比ν表示：
在该势能方程中，状态关系可以表示为：
其中，
ρ表示材料密度
表示自由能ψ关于弹性张量ε e的微分；
表示自由能ψ关于损伤D的微分；
Y(ε e,D)表示弹性损伤释放率；
自由能的微分形式可表示为式(8)：
其中，
为弹性损伤势能的微分形式；
为弹性张量的微分形式；
为损伤变量的微分形式。
3.根据权利要求1所述的考虑塑性和损伤本构关系的有限元模拟方法，其特征在于，所述步骤S3：
为了耦合损伤行为，控制应力空间内部变量的演化规律，定义势能方程如下：
F(σ，Y；D)＝fp+FY (9)
其中，
F(σ,Y；D)表示势能方程；
塑性势能fp和损伤势能FY计算方式如式(10)：
其中，
为应力张量第二不变量；
参数Y0,α,β和γ控制损伤势能演化；
σy为屈服应力，通过Johnson-Cook塑性流动准则进行求解；
根据增量理论，温度增量和应变速率/>可由式(11)进行计算：
其中，
η为自定义非弹性热系数；
为塑性应变率；
ρ为材料密度；
c为比热容；
应变速率为时间间隔Δt内的平均塑性应变；
根据上述应变速率的定义，Johnson-Cook硬化模型的模块可表示为:
其中，
σy表示屈服应力；
表示塑性应变；
表示初始应变率；
Johnson-Cook硬化模型中等效塑性应变相关变量的求导为：
其中，
η表示非弹性热系数；
R表示应力；
ρ表示材料密度；
c表示比热容；
对于特定的应力张量σ，选择的坐标轴不同，应力张量的矩阵形式也会有所不同。但是存在三个应力张量不变量：
其中，
J1表示应力偏张量第一不变量；
J2表示应力偏张量第二不变量；
J3表示应力偏张量第三不变量；
σH表示平均应力；
σ表示柯西应力张量；
σ1表示第一主应力；
σ2表示第二主应力；
σ3表示第三主应力；
σeq表示米塞斯等效应力；
det表示一个函数运算，类似于加减乘除，在数学计算中常用的一个函数；
S表示应力偏张量；
所以应力状态可用三个应力张量不变量来表示，为了简化，引入Lode角θ和应力三轴度ξ，在主应力空间中，通过坐标原点的偏应力平面为π平面；
Lode角表示应力张量与最大主应力轴夹角的投影，应力三轴度表示静水应力与等效应力的比值
应力状态用主应力空间中(σm,σeq,θ)来表示；
Lode角θ往往用Lode系数来表示使得应力三轴度ξ与Lode系数η均定义在[-1,1]内；
所有应力方向都可以通过上述定义的参数(ξ,θ)来表征；
为了考虑Lode角θ和应力三轴度ξ对损伤演化的影响，引入了参数h(ξ,θ)(0≤h(ξ,θ)≤1)；
h(ξ1,θ1)＝1表示在状态(ξ1,θ1)时拥有最快的损伤演化；h(ξ2,θ2)＝0表示在状态(ξ2,θ2)时的没有损伤演化；
应力状态参数h(ξ,θ)不仅可以考虑微裂纹闭合的效果，还可以区分剪切状态与拉压状态；
当用耦合了Lode角和应力三轴度的应力状态参数h(ξ,θ)替换了单轴损伤演化参数h后，一个更精确的本构模型被提出：
其中，
ρΨ(ε e,D)表示弹性损伤势能
等效应力张量和损伤释放率可以用式(16)进行更新：
其中，
(ξ,θ)表示应力状态参数
Y(ε e,D)表示损伤释放率
通过耗散分析，累积等效塑性应变和损伤生长速率按式(17)更新：
其中，
表示累积等效塑性应变；
表示损伤生长速率；
J2(σ)表示应力偏张量第二不变量；
表示塑性系数；
Y0表示初始损伤参数；
β表示损伤相关系数；
应力状态参数明显耦合为塑性应变和损伤演化，损伤演化与应力状态参数有密切的关系，应力状态参数越高，损伤演化越快，低应力状态参数也导致损伤演化较慢，特别是当h(ξ,θ)＝0时，损伤没有演化。
4.一种考虑塑性和损伤本构关系的有限元模拟系统，其特征在于，包括：
利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为；
所述Johnson-Cook塑性流动准则：
其中，
表示应力；
表示应变率；
A表示屈服应力；
B表示硬化模量；
表示塑性应变；
C表示黏性参数；
ln表示数学运算符号，表示无量纲塑性应变率；
T表示当前温度；
T0表示室温；
Tm表示熔化温度；
n表示硬化指数；
m表示温度软化参数；
所述利用连续损伤力学和Johnson-Cook塑性流动准则建立本构方程，以精确预测金属成形和加工过程中的损伤行为，包括：
模块S1：损伤的宏观表述；
模块S2：耦合损伤与弹塑性；
模块S3：耦合损伤与Johnson-Cook硬化模型；
所述模块S1：
以损伤因子D表示材料实际承载面积的减少量，当材料处于初始状态，无损伤时，D＝0；当D＝1时表示材料完全破裂损伤；
基于等效应变假设其中/>将材料损伤中微孔洞的作用描述为材料弹性模量的变化；
通过引入单轴损伤演化参数h(0≤h≤1)来描述微裂纹在压缩加载状态下闭合效应，h＝1表示在没有裂纹闭合效应的情况下，拉伸状态下的损伤演化；h＝0表示裂纹闭合，纯压缩状态下无损伤演化；
根据弹性能量等效理论，在主坐标系下的各向同性损伤等效变量如式(2)：
其中，
表示损伤柯西应力张量；
D表示损伤变量；
σ表示柯西应力张量；
表示单边应力状态参数；
表示损伤弹性张量；
ε e表示弹性张量；
We表示弹性势能；
E表示弹性模量张量；
表示损伤弹性模量张量。
5.根据权利要求4所述的考虑塑性和损伤本构关系的有限元模拟系统，其特征在于，所述模块S2：
基于Helmholtz自由能方程，考虑弹性应变张量εe和损伤变量D，弹性损伤势能方程可以表示为式(3)：
其中，
ρΨ(εe,)表示弹性损伤势能；
Ψ表示自由能；
tr表示运算函数；
λe,μe均为拉梅系数；
ρ表示材料密度；
Lame系数(λe,μe)可以用杨氏模量E和泊松比ν表示：
在该势能方程中，状态关系可以表示为：
其中，
ρ表示材料密度
表示自由能ψ关于弹性张量εe的微分；
表示自由能ψ关于损伤D的微分；
Y(εe,D)表示弹性损伤释放率；
自由能的微分形式可表示为式(8)：
其中，
为弹性损伤势能的微分形式；
为弹性张量的微分形式；
为损伤变量的微分形式。
6.根据权利要求4所述的考虑塑性和损伤本构关系的有限元模拟系统，其特征在于，所述模块S3：
为了耦合损伤行为，控制应力空间内部变量的演化规律，定义势能方程如下：
F(σ,Y；D)＝fp+FY (9)
其中，
F(σ,；)表示势能方程；
塑性势能fp和损伤势能FY计算方式如式(10)：
其中，
为应力张量第二不变量；
参数Y0,α,β和γ控制损伤势能演化；
σy为屈服应力，通过Johnson-Cook塑性流动准则进行求解；
根据增量理论，温度增量和应变速率/>可由式(11)进行计算：
其中，
η为自定义非弹性热系数；
为塑性应变率；
ρ为材料密度；
c为比热容；
应变速率为时间间隔Δt内的平均塑性应变；
根据上述应变速率的定义，Johnson-Cook硬化模型的模块可表示为:
其中，
σy表示屈服应力；
表示塑性应变；
表示初始应变率；
Johnson-Cook硬化模型中等效塑性应变相关变量的求导为：
其中，
η表示非弹性热系数；
R表示应力；
ρ表示材料密度；
c表示比热容；
对于特定的应力张量σ，选择的坐标轴不同，应力张量的矩阵形式也会有所不同。但是存在三个应力张量不变量：
其中，
J1表示应力偏张量第一不变量；
J2表示应力偏张量第二不变量；
J3表示应力偏张量第三不变量；
σH表示平均应力；
σ表示柯西应力张量；
σ1表示第一主应力；
σ2表示第二主应力；
σ3表示第三主应力；
σeq表示米塞斯等效应力；
det表示一个函数运算，类似于加减乘除，在数学计算中常用的一个函数；
S表示应力偏张量；
所以应力状态可用三个应力张量不变量来表示，为了简化，引入Lode角θ和应力三轴度ξ，在主应力空间中，通过坐标原点的偏应力平面为π平面；
Lode角表示应力张量与最大主应力轴夹角的投影，应力三轴度表示静水应力与等效应力的比值
应力状态用主应力空间中(σm,σeq,θ)来表示；
Lode角θ往往用Lode系数来表示使得应力三轴度ξ与Lode系数η均定义在[-1,1]内；
所有应力方向都可以通过上述定义的参数(ξ,θ)来表征；
为了考虑Lode角θ和应力三轴度ξ对损伤演化的影响，引入了参数h(ξ,θ)(0≤h(ξ,θ)≤1)；
h(ξ1,θ1)＝1表示在状态(ξ1,θ1)时拥有最快的损伤演化；h(ξ2,θ2)＝0表示在状态(ξ2,θ2)时的没有损伤演化；
应力状态参数h(ξ,θ)不仅可以考虑微裂纹闭合的效果，还可以区分剪切状态与拉压状态；
当用耦合了Lode角和应力三轴度的应力状态参数h(ξ,θ)替换了单轴损伤演化参数h后，一个更精确的本构模型被提出：
/>
其中，
ρΨ(ε e,D)表示弹性损伤势能
等效应力张量和损伤释放率可以用式(16)进行更新：
其中，
(ξ,θ)表示应力状态参数
Y(ε e,D)表示损伤释放率
通过耗散分析，累积等效塑性应变和损伤生长速率按式(17)更新：
其中，
表示累积等效塑性应变；
表示损伤生长速率；
J2(σ)表示应力偏张量第二不变量；
表示塑性系数；
Y0表示初始损伤参数；
β表示损伤相关系数；
应力状态参数明显耦合为塑性应变和损伤演化，损伤演化与应力状态参数有密切的关系，应力状态参数越高，损伤演化越快，低应力状态参数也导致损伤演化较慢，特别是当h(ξ,θ)＝0时，损伤没有演化。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010732552.7A CN111709176B (zh) | 2020-07-27 | 2020-07-27 | 考虑塑性和损伤本构关系的有限元模拟方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010732552.7A CN111709176B (zh) | 2020-07-27 | 2020-07-27 | 考虑塑性和损伤本构关系的有限元模拟方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111709176A CN111709176A (zh) | 2020-09-25 |
| CN111709176B true CN111709176B (zh) | 2023-09-12 |
Family
ID=72547796
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010732552.7A Active CN111709176B (zh) | 2020-07-27 | 2020-07-27 | 考虑塑性和损伤本构关系的有限元模拟方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111709176B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112347574B (zh) * | 2020-10-29 | 2022-05-24 | 中国石油大学(华东) | 一种增材制造金属材料的弹塑性损伤失效预测方法 |
| CN113642207A (zh) * | 2021-07-09 | 2021-11-12 | 恒大新能源技术(深圳)有限公司 | 金属失效模型构建方法、装置、终端设备及存储介质 |
| CN113764056B (zh) * | 2021-09-06 | 2023-04-07 | 北京理工大学重庆创新中心 | 一种获得材料多应变率下高精度硬化模型参数的方法 |
| CN114462146B (zh) * | 2022-01-28 | 2023-02-03 | 中国人民解放军陆军工程大学 | 考虑老化损伤的推进剂蠕变本构模型的构建与有限元应用方法 |
| CN115482890B (zh) * | 2022-09-14 | 2023-10-31 | 东南大学 | 基于弹性应变能密度分解的局部梯度损伤模型构建方法 |
| CN117497069B (zh) * | 2023-10-23 | 2024-05-24 | 华中科技大学 | 一种高聚物材料的超弹性本构模型的构建方法和装置 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407526A (zh) * | 2016-09-05 | 2017-02-15 | 大连理工大学 | 一种微铣削过程刀具后刀面磨损预测方法 |
| CN107655768A (zh) * | 2017-08-04 | 2018-02-02 | 中车青岛四方机车车辆股份有限公司 | 基于拉伸剪切和缺口拉伸试验的损伤模型参数标定方法 |
| CN111144044A (zh) * | 2019-12-16 | 2020-05-12 | 重庆大学 | 一种考虑温度影响的塑料齿轮接触疲劳寿命评估方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080078479A1 (en) * | 2006-09-29 | 2008-04-03 | Massachusetts Institute Of Technology | Determining constitutive models for fracture prediction in ductile materials undergoing plastic deformation |
-
2020
- 2020-07-27 CN CN202010732552.7A patent/CN111709176B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407526A (zh) * | 2016-09-05 | 2017-02-15 | 大连理工大学 | 一种微铣削过程刀具后刀面磨损预测方法 |
| CN107655768A (zh) * | 2017-08-04 | 2018-02-02 | 中车青岛四方机车车辆股份有限公司 | 基于拉伸剪切和缺口拉伸试验的损伤模型参数标定方法 |
| CN111144044A (zh) * | 2019-12-16 | 2020-05-12 | 重庆大学 | 一种考虑温度影响的塑料齿轮接触疲劳寿命评估方法 |
Non-Patent Citations (1)
| Title |
|---|
| 张杰.基于全耦合损伤模型的数值模拟研究.《电子技术与软件工程》.2020,摘要和第1-4节. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111709176A (zh) | 2020-09-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111709176B (zh) | 考虑塑性和损伤本构关系的有限元模拟方法及系统 | |
| Fei et al. | Springback prediction for incremental sheet forming based on FEM-PSONN technology | |
| Ponthot et al. | A cascade optimization methodology for automatic parameter identification and shape/process optimization in metal forming simulation | |
| Lee et al. | Shape optimization of the workpiece in the forging process using equivalent static loads | |
| CN110211645B (zh) | 微观-宏观尺度钣金成形工艺模型的损伤与疲劳寿命评估方法 | |
| CN112347574B (zh) | 一种增材制造金属材料的弹塑性损伤失效预测方法 | |
| Wang et al. | Evaluating interactions between the heavy forging process and the assisting manipulator combining FEM simulation and kinematics analysis | |
| Qin et al. | Position error compensation of the multi-purpose overload robot in nuclear power plants | |
| Wang et al. | Uncalibrated visual positioning using adaptive Kalman Filter with dual rate structure for wafer chip in LED packaging | |
| Mungi et al. | Residual stresses in cold axisymmetric forging | |
| Guo et al. | Analysis on quantifiable and controllable assembly technology for aeronautical thin-walled structures | |
| Glorieux et al. | Quality and productivity driven trajectory optimisation for robotic handling of compliant sheet metal parts in multi-press stamping lines | |
| Quan et al. | Numerical description of hot flow behaviors at Ti-6Al-2Zr-1Mo-1V alloy by GA-SVR and relative applications | |
| CN116702631A (zh) | 基于人工神经网络的电子束增材制造本构关系计算方法 | |
| Wan et al. | Non-geometric error compensation for long-stroke cartesian robot with semi-analytical beam deformation and gaussian process regression model | |
| Tran et al. | A computational design of robotic grasper by intelligence-based topology optimization for microassembly and micromanipulation | |
| Poulhaon et al. | Online prediction of machining distortion of aeronautical parts caused by re-equilibration of residual stresses | |
| Berthold et al. | Investigation of lateral compression effects in fiber reinforced soft pneumatic actuators | |
| Nguyen et al. | Dynamic Parameter Identification based on the Least Squares method for a 6-DOF Manipulator | |
| My et al. | Design analysis for a special serial-parallel manipulator | |
| Yeh et al. | Application of adaptive network fuzzy inference system to die shape optimal design in sheet metal bending process | |
| Gun | 3D boundary element analysis of creep continuum damage mechanics problems | |
| Nowak et al. | A plane stress elastic-plastic analysis of sheet metal cup deep drawing processes | |
| Esener et al. | A Sensitivity Analysis By Using Design Of Experiment And Its Application In Stamping | |
| Behrouzi et al. | A new approach for inverse analysis of springback in a sheet-bending process |
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 |