板材冲压成形的晶体塑性有限元法模拟
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
目前在晶体塑性有限元模拟中都采用实体单元，而壳单元模型在板材冲压成形模拟中效 率更高；计算中大都采用随机选取离散集合中大量晶体方位的方法实现晶体取向分布规律向 有限元模型的映射，计算量大。为了提高晶体塑性有限元模拟的效率，将晶体塑性本构模型 引入 Mindlin 曲壳单元与动力显式有限元法，根据织构组分的正态分布特性在单元积分点处 建立细观晶体模型与宏观有限元模型的映射关系。开发晶体塑性有限元程序，对板材冲压成 形过程进行数值模拟，预示成形中板材的变形行为与织构演化特性。
5. 结论
本文基于率相关多晶体塑性动力显式有限元模型，对退火铝板的方盒件拉深成形过程进 行了数值模拟，较好的预示了制耳的形成，并和实验结果相一致。数值模拟结果表明了本文 所建立模型的有效性，能够很好地反映板材的初始织构特性对成形性的影响。
A
A
(a) 计算结果
(b) 试验结果
图 4 盒形件成形构形模拟结果与试验的比较
设变形前第s个滑移系滑移方向的单位矢量以n(s)表示，滑移面的单位法矢量为m(s)。晶
体塑性变形后产生晶格畸变，滑移方向及滑移面法向将变为n*(s)和 m*(s)
n*(s) = F e ⋅ n(s) m*(s) = m (s) ⋅ (F e )−1
(4)
设变形中在滑移系s中由于滑移产生的塑性剪切应变率为 γ&(s) ，塑性应变率张量 ε& p' 及塑
=
x,ξ
/
x,ξ
⎨eˆ3 = ( x,ξ ×x,η ) / x,ξ ×x,η
(17)
⎪ ⎩
eˆ2
=
eˆ3
× eˆ1
在单元高斯积分点局部坐标系下的应变率 ε&ˆ 可以由整体坐标系下的节点速度 u& 来计算
[7]
εˆ& = Bˆ u&
(18)
在高斯点局部坐标系下，本构关系方程为
∇
σˆ ij
=
De ijkl
2
2.2 滑移系的确定及塑性应变增量计算
对于剪切应变速率 γ&(s) ，采用率相关形式[10]：
γ&
(s
)＝a&
(
s
)
⎡ ⎢ ⎣
τ g
(s) (s)
⎤⎡ ⎥⎢ ⎦⎣
τ (s) g (s)
⎤ (1/ m)−1 ⎥ ⎦
(9)
式中 τ (s)源自文库—滑移系 s 中的切应力
a&(s) —滑移系上的参考剪切应变速率
=
De ijkl
(∆εˆkl
− ∆εˆkpl ) − Γ ijkl∆εˆkl
(22)
Γ ijkl
=
1 2
(σˆ
lj
δik
+ σˆ kj δlj
+ σˆli δkj
+ σˆ ki δlj )
(23)
于是整体坐标系下的第二类 Kirchihoff 应力增量为
则整体坐标系下的 Cauchy 应力可以更新为 式中 ｜J｜—雅可比矩阵行列式
( 0 ≤ θ ≤ 1 ) (10)
-2-
http://www.paper.edu.cn
将式（12）进行台劳展开并代入（11）式，可以得到
∑ ( ) N sj∆γ(s) = γ&ts + Q(s) : ε ∆t
(11)
( j)
Q(s)
=
⎜⎜⎝⎛
θ∆tγ&t( mτ (s
s )
)
⎟⎟⎠⎞ R ( s )
∆Sij = Θil ∆Sˆlk Θ jk
(24)
( ) σ n+1 ij
=
1 J
∂xin+1
∂x
n+1 j
∂x
n k
∂xln
Sinj + ∆Sij
(25)
板材弹塑性变形的动力学方程为：
∂σ ij ∂x j
+ Fi
− ρu&&i
− γu&i
=
0
(26)
式中 σij —Cauchy应力张量
Fi —外力 ρ —金属板料的质量密度
γ —阻尼系数
u&i , u&&i —板料速度和加速度
根据散度定理以及边界条件，系统的虚功方程可以表示为
∫ ∫ ∫ ∫ V ρu&&i δu& i dV + V γu& i δu& i dV = V piδu& i dV + V σ ij δεij dV
(27)
进行有限元离散之后，得到有限元方程
mu&& + cu& = Fout − Fint
uiN
和绕
v
1N i
、
v i2 N
的二个角位移
α1N
、
α
N 2
插值表示为
∑ ∑ ui
=
4
ψ N uiN
N =1
+ζ 2
4
ψ Nδ
N =1
N (vi1N α1N
−
v
2 i
N
α
N 2
)
(16)
为描述应力应变关系，在单元高斯积分点处建立局部坐标系 xˆ i （i=1,2,3），其基矢量
eˆi 为
⎧ ⎪
eˆ1
1. 引 言
金属板材在成形过程中通常会显示出各向异性，这种各向异性对板材的成形性能有重要 的影响。作为多晶体材料，金属板材产生塑性各向异性的原因，主要是材料在加工过程中形 成晶粒的择优取向，即织构，随着后续塑性加工的进行，织构还将发生演化。
为了反映板材的塑性各向异性，Hill[1]、Gotoh[2]、Barlat[3]等人提出了多个宏观现象学 的屈服函数及相应的材料本构关系模型，并应用于板材冲压成形模拟。在应用这些宏观现象 学的各向异性本构模型过程中，模拟的准确性往往与屈服函数的形式密切相关。模型中的参 数通过简单实验得到，一般保持不变，因而无法考虑塑性变形后的后继屈服轨迹与塑性变形 的相关性。
2. 晶体塑性有限元模型
1 本课题得到教育部博士点基金(20030248029)的资助。
-1-
http://www.paper.edu.cn
2.1 晶体塑性本构方程
在每个晶体上沿晶体方向[100], [010]和[001]建立晶体坐标系 xi' (i=1,2,3），宏观整体坐
标系 xi(i=1,2,3）与晶体坐标系 xi' (i=1,2,3）之间的关系可以表示为
性旋转率张量 ω& p 可以表示为
∑ ε&p' = P (s)γ&(s)
(5)
s
∑ ω& p = W (s)γ&(s)
(6)
s
P (s) = 1 (m *(s) n *(s) +n *(s) m *(s) )
(7)
2
W (s) = 1 (m *(s) n *(s) −n *(s) m *(s) )
(8)
晶体塑性理论的发展为精确描述材料塑性变形中的各向异性提供了理论基础，它从塑性 变形的物理本质出发，能够在细观尺度下描述板材的塑性各向异性。晶体塑性理论的研究起 源于二十世纪二十年代 Taylor 以及随后 Bishop 与 Hill 的工作，近三十年来随着细观力学、 实验测试技术、材料科学的发展，使得人们能够从材料微结构的变化对材料的塑性变形行为 进行理论分析和实验验证。基于晶体塑性理论的本构模型以位错理论为基础，可以直接引入 材料的晶体学织构并考虑其在变形中的演化，近几年来已被逐渐应用到板材冲压成形的模拟 中。Anand 等[4] 采用单晶体面心金属模拟冲压中的制耳现象。Nakamachi 等人[5-9] 建立 了弹塑性多晶体有限元模型，并应用于金属板成形分析。
(28)
式中 m —一致质量矩阵
c —阻尼矩阵
Fout —节点外力列矢量 Fint —为节点内力列矢量。
采用中心差分时间积分
-4-
http://www.paper.edu.cn
4. 方盒件冲压成形模拟
基于以上晶体塑性有限元模型编制了有限元计算程序，对方盒件冲压成形过程进行了数 值模拟。
图 2 所示为方盒件冲压成形示意图，板料为轧制铝板，材料参数为：弹性模量为 69GPa, 泊松比为 0.3, 板厚为 0.8mm，应变硬化指数为 0.2，滑移系初始临界剪切应力为 69MPa,速 率敏感指数为 0.001, 参考剪切应变率为 0.001, 初始硬化率为 614MPa。冲头行程为 14mm。 图 3 所示为板料的初始晶体取向分布函数（ODF）图，退火铝板主要由再结晶、立方织构及 铜织构和 S 织构组成。
(12)
R(s) = C e : P(s) +W (s) ⋅ σ − σ ⋅W (s)
(13)
N sj
= δsj
+
⎜⎜⎝⎛
θ∆tγ& t( mτ (s
s )
)
⎟⎟⎠⎞⎢⎣⎡
R
(
s)
τ
:P
(s)
(
j
)
+ sgn(τ ( j) )
H sj ⎤
g
(s)
⎥ ⎦
(14)
式中 H st —硬化模量。
3. Mindlin 曲壳单元
图 2 方盒件成形示意图
图 3 板料初始晶体取向分布
图 4 所示为模拟得到的盒形件最终构形及与试验的对比，图 5 为成形后 ODF 的计算结 果及试验对比，二者趋势基本吻合；由于在初始织构离散与应力计算中忽略了随机取向，模 拟结果具有较强的密度峰值。在方盒件冲压成形过程中，铜织构和 S 织构为不稳定取向， 变形后逐渐转到其他取向； 而立方织构为稳定取向，变形后取向密度值增加。
http://www.paper.edu.cn
板材冲压成形的晶体塑性有限元法模拟1
彭颖红，张少睿，李大永
上海交通大学机械与动力学院，上海（200240）
E-mail：shaoruiz@sjtu.edu.cn
摘 要：金属板材的织构特性及其演化是板材各向异性的主要原因。本文将率相关晶体塑性 本构理论引入 Mindlin 曲壳单元模型与动力显式有限元法，以潜在硬化模型描述应变硬化， 根据取向分布函数在取向空间中的分布规律将晶体取向分配给各个单元的积分点，建立了弹 塑性大变形条件下的率相关多晶体塑性模型，并将其引入动力显式有限元法。对退火铝板方 盒件拉深成形过程进行了模拟，并与实验结果进行了对比分析，计算结果与试验结果显示了 较好的一致性。 关键词：晶体塑性模型，织构，动力显式有限元，Mindlin 壳单元 中图分类号：TG302
局部坐标系下的塑性应变率与整体坐标系下的塑性应变率之间可由下式转换
εˆ&ipj = Θli ε&lpkΘkj
(20) (21)
Θ 为局部坐标系与整体坐标系之间的变换张量。整体坐标系下塑性应变率 ε& p 通过晶体塑性
本构模型计算得到。
由式（19）可以得到第二类 Kirchihoff 应力增量为
∆Sˆ ij
ε&ˆkel
=
De ijkl
(εˆ&kl
− εˆ&kpl )
(19)
-3-
http://www.paper.edu.cn
∇
式中 σˆ ij —局部坐标系下 Cauchy 应力 σˆij 的 Jaumann 率张量
De ijkl
=
De ijkl
−
D D e e ij33 kl 33
/
De 3333
式中 De —四阶弹性本构张量
g (s) —滑移系上的临界切应力
m—速率敏感指数
为提高计算的稳定性，采用 PIERCE 等[11]提出的切线系数法。在时间步长∆t 内增量 ∆γ(s)
可以由 t 和 t+∆t 时刻的 γ&(s) 线性插值得到，即 ∆γ(s) = [θγ&(s)(t + ∆t) + (1− θ)γ&(s)(t)]∆t
δ N —节点 N 处的厚度 vi3N —节点 N 处中面的单位法矢量。
图 1 四节点 Mindlin 曲壳单元示意图
设
v
1N i
和
v i2 N
分别为与
v
3N i
相互正交的单位矢量，假定壳体的中面法线变形之后仍保持
为直线，但不一定再垂直于变形后的壳中面，则单元内任一点的位移可由各中面节点的三个
线位移分量
xi' = Λij x j
(1)
式中 Λij—包含欧拉角(ϕ1,φ,ϕ2 )坐标转换张量[8]
两坐标系下塑性应变率之间的转换关系为
ε&ipj' = Λik ε&kpl Λjl
(2)
晶体的变形梯度可通过乘法分解表示为
F = Fe ⋅Fp
(3)
式中 FP—对应于滑移面上的塑性剪切变形
Fe—对应于弹性变形与晶格的转动
图 1 所示四节点 Mindlin 曲壳单元由上下两个 曲面及周边以壳体厚度方向的直线为母线的曲面 所围成。设 xi(i=1,2,3）为宏观整体坐标系，单元坐 标插值分别为
∑ ∑ xi
=
4
ψ N xiN
N =1
+
4 N =1
ζ 2
ψ
N
δ
N
v
3N i
(15)
式中 ψ N —四节点线性单元的形函数
η
ζ
ξ
-5-
http://www.paper.edu.cn
(a) 计算结果
(b) 试验结果
图 5 盒形件成形后 ODF 计算结果及试验对比
-6-
http://www.paper.edu.cn
参考文献
[1] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Royal Society of London Proceedings, 1948, A: 193-281 [2] Gotoh M. A finite element analysis of the rigid plastic deformation of the flange in a deep drawing process based on a fourth-degree yield function, International Journal of Mechanical Science, 1980, 22: 367-377 [3] Barlat F, Lian J. Plastic behavior and stretchability of sheet metals. Part I: a yield function for orthotropic sheet under plane stress conditions. International Journal of Plasticity, 1989, 51～66 [4] Anand L, Balasubramanian S. Polycrystal Plasticity: Application to Earing in Cup Drawing, Annals of the CIRP, 1996, Vol. 45(1):263-268 [5] Xie C L, Nakamachi E. Investigations of the formability of BCC steel sheets by using crystalline plasticity finite element analysis. Materials and Design, 2002, 23: 59-68 [6] Nakamachi E, Xie C L. Harimato M., Drawability assessment of Bcc steel by using elastic/crystalline viscoplastic finite element analysis. Int. J. Mech. Sci., 2001, 43: 631-652 [7] Nakamachi E, Hiraiwa K, Morimoto H, Harimoto M. Elastic/crystalline viscoplastic finite element analyses of single- and poly-crystal sheet deformations and their experimental verification, Int. J. of Plasticity, 2000, 16:1419-1441 [8] Nakamachi E, Dong X. Study of texture effect on sheet failure in a limit dome height test by using elastic/crystalline viscoplastic finite element analysis, J. Appl. Mech. Trans. ASME(E), 1997, 64: 519-524 [9] Nakamachi E, Huo T. Dynamic-explicit elastic plastic finite-element simulation of hemispherical punch-drawing of sheet metal. Engineering Computations, 1996, 13: 327-338 [10] Pierce D, Asaro R J and Needleman A. Material rate dependent and localized deformation in crystalline solids[J]. Acta Metall., 1983, 31:1951 ~1972 [11] Pierce D, Shih C F and Needleman A. A tangent modulus method for rate dependent solids[J]. Computers & Structures, 1984, 18: 875-887.
5. 结论
本文基于率相关多晶体塑性动力显式有限元模型，对退火铝板的方盒件拉深成形过程进 行了数值模拟，较好的预示了制耳的形成，并和实验结果相一致。数值模拟结果表明了本文 所建立模型的有效性，能够很好地反映板材的初始织构特性对成形性的影响。
A
A
(a) 计算结果
(b) 试验结果
图 4 盒形件成形构形模拟结果与试验的比较
设变形前第s个滑移系滑移方向的单位矢量以n(s)表示，滑移面的单位法矢量为m(s)。晶
体塑性变形后产生晶格畸变，滑移方向及滑移面法向将变为n*(s)和 m*(s)
n*(s) = F e ⋅ n(s) m*(s) = m (s) ⋅ (F e )−1
(4)
设变形中在滑移系s中由于滑移产生的塑性剪切应变率为 γ&(s) ，塑性应变率张量 ε& p' 及塑
=
x,ξ
/
x,ξ
⎨eˆ3 = ( x,ξ ×x,η ) / x,ξ ×x,η
(17)
⎪ ⎩
eˆ2
=
eˆ3
× eˆ1
在单元高斯积分点局部坐标系下的应变率 ε&ˆ 可以由整体坐标系下的节点速度 u& 来计算
[7]
εˆ& = Bˆ u&
(18)
在高斯点局部坐标系下，本构关系方程为
∇
σˆ ij
=
De ijkl
2
2.2 滑移系的确定及塑性应变增量计算
对于剪切应变速率 γ&(s) ，采用率相关形式[10]：
γ&
(s
)＝a&
(
s
)
⎡ ⎢ ⎣
τ g
(s) (s)
⎤⎡ ⎥⎢ ⎦⎣
τ (s) g (s)
⎤ (1/ m)−1 ⎥ ⎦
(9)
式中 τ (s)源自文库—滑移系 s 中的切应力
a&(s) —滑移系上的参考剪切应变速率
=
De ijkl
(∆εˆkl
− ∆εˆkpl ) − Γ ijkl∆εˆkl
(22)
Γ ijkl
=
1 2
(σˆ
lj
δik
+ σˆ kj δlj
+ σˆli δkj
+ σˆ ki δlj )
(23)
于是整体坐标系下的第二类 Kirchihoff 应力增量为
则整体坐标系下的 Cauchy 应力可以更新为 式中 ｜J｜—雅可比矩阵行列式
( 0 ≤ θ ≤ 1 ) (10)
-2-
http://www.paper.edu.cn
将式（12）进行台劳展开并代入（11）式，可以得到
∑ ( ) N sj∆γ(s) = γ&ts + Q(s) : ε ∆t
(11)
( j)
Q(s)
=
⎜⎜⎝⎛
θ∆tγ&t( mτ (s
s )
)
⎟⎟⎠⎞ R ( s )
∆Sij = Θil ∆Sˆlk Θ jk
(24)
( ) σ n+1 ij
=
1 J
∂xin+1
∂x
n+1 j
∂x
n k
∂xln
Sinj + ∆Sij
(25)
板材弹塑性变形的动力学方程为：
∂σ ij ∂x j
+ Fi
− ρu&&i
− γu&i
=
0
(26)
式中 σij —Cauchy应力张量
Fi —外力 ρ —金属板料的质量密度
γ —阻尼系数
u&i , u&&i —板料速度和加速度
根据散度定理以及边界条件，系统的虚功方程可以表示为
∫ ∫ ∫ ∫ V ρu&&i δu& i dV + V γu& i δu& i dV = V piδu& i dV + V σ ij δεij dV
(27)
进行有限元离散之后，得到有限元方程
mu&& + cu& = Fout − Fint
uiN
和绕
v
1N i
、
v i2 N
的二个角位移
α1N
、
α
N 2
插值表示为
∑ ∑ ui
=
4
ψ N uiN
N =1
+ζ 2
4
ψ Nδ
N =1
N (vi1N α1N
−
v
2 i
N
α
N 2
)
(16)
为描述应力应变关系，在单元高斯积分点处建立局部坐标系 xˆ i （i=1,2,3），其基矢量
eˆi 为
⎧ ⎪
eˆ1
1. 引 言
金属板材在成形过程中通常会显示出各向异性，这种各向异性对板材的成形性能有重要 的影响。作为多晶体材料，金属板材产生塑性各向异性的原因，主要是材料在加工过程中形 成晶粒的择优取向，即织构，随着后续塑性加工的进行，织构还将发生演化。
为了反映板材的塑性各向异性，Hill[1]、Gotoh[2]、Barlat[3]等人提出了多个宏观现象学 的屈服函数及相应的材料本构关系模型，并应用于板材冲压成形模拟。在应用这些宏观现象 学的各向异性本构模型过程中，模拟的准确性往往与屈服函数的形式密切相关。模型中的参 数通过简单实验得到，一般保持不变，因而无法考虑塑性变形后的后继屈服轨迹与塑性变形 的相关性。
2. 晶体塑性有限元模型
1 本课题得到教育部博士点基金(20030248029)的资助。
-1-
http://www.paper.edu.cn
2.1 晶体塑性本构方程
在每个晶体上沿晶体方向[100], [010]和[001]建立晶体坐标系 xi' (i=1,2,3），宏观整体坐
标系 xi(i=1,2,3）与晶体坐标系 xi' (i=1,2,3）之间的关系可以表示为
性旋转率张量 ω& p 可以表示为
∑ ε&p' = P (s)γ&(s)
(5)
s
∑ ω& p = W (s)γ&(s)
(6)
s
P (s) = 1 (m *(s) n *(s) +n *(s) m *(s) )
(7)
2
W (s) = 1 (m *(s) n *(s) −n *(s) m *(s) )
(8)
晶体塑性理论的发展为精确描述材料塑性变形中的各向异性提供了理论基础，它从塑性 变形的物理本质出发，能够在细观尺度下描述板材的塑性各向异性。晶体塑性理论的研究起 源于二十世纪二十年代 Taylor 以及随后 Bishop 与 Hill 的工作，近三十年来随着细观力学、 实验测试技术、材料科学的发展，使得人们能够从材料微结构的变化对材料的塑性变形行为 进行理论分析和实验验证。基于晶体塑性理论的本构模型以位错理论为基础，可以直接引入 材料的晶体学织构并考虑其在变形中的演化，近几年来已被逐渐应用到板材冲压成形的模拟 中。Anand 等[4] 采用单晶体面心金属模拟冲压中的制耳现象。Nakamachi 等人[5-9] 建立 了弹塑性多晶体有限元模型，并应用于金属板成形分析。
(28)
式中 m —一致质量矩阵
c —阻尼矩阵
Fout —节点外力列矢量 Fint —为节点内力列矢量。
采用中心差分时间积分
-4-
http://www.paper.edu.cn
4. 方盒件冲压成形模拟
基于以上晶体塑性有限元模型编制了有限元计算程序，对方盒件冲压成形过程进行了数 值模拟。
图 2 所示为方盒件冲压成形示意图，板料为轧制铝板，材料参数为：弹性模量为 69GPa, 泊松比为 0.3, 板厚为 0.8mm，应变硬化指数为 0.2，滑移系初始临界剪切应力为 69MPa,速 率敏感指数为 0.001, 参考剪切应变率为 0.001, 初始硬化率为 614MPa。冲头行程为 14mm。 图 3 所示为板料的初始晶体取向分布函数（ODF）图，退火铝板主要由再结晶、立方织构及 铜织构和 S 织构组成。
(12)
R(s) = C e : P(s) +W (s) ⋅ σ − σ ⋅W (s)
(13)
N sj
= δsj
+
⎜⎜⎝⎛
θ∆tγ& t( mτ (s
s )
)
⎟⎟⎠⎞⎢⎣⎡
R
(
s)
τ
:P
(s)
(
j
)
+ sgn(τ ( j) )
H sj ⎤
g
(s)
⎥ ⎦
(14)
式中 H st —硬化模量。
3. Mindlin 曲壳单元
图 2 方盒件成形示意图
图 3 板料初始晶体取向分布
图 4 所示为模拟得到的盒形件最终构形及与试验的对比，图 5 为成形后 ODF 的计算结 果及试验对比，二者趋势基本吻合；由于在初始织构离散与应力计算中忽略了随机取向，模 拟结果具有较强的密度峰值。在方盒件冲压成形过程中，铜织构和 S 织构为不稳定取向， 变形后逐渐转到其他取向； 而立方织构为稳定取向，变形后取向密度值增加。
http://www.paper.edu.cn
板材冲压成形的晶体塑性有限元法模拟1
彭颖红，张少睿，李大永
上海交通大学机械与动力学院，上海（200240）
E-mail：shaoruiz@sjtu.edu.cn
摘 要：金属板材的织构特性及其演化是板材各向异性的主要原因。本文将率相关晶体塑性 本构理论引入 Mindlin 曲壳单元模型与动力显式有限元法，以潜在硬化模型描述应变硬化， 根据取向分布函数在取向空间中的分布规律将晶体取向分配给各个单元的积分点，建立了弹 塑性大变形条件下的率相关多晶体塑性模型，并将其引入动力显式有限元法。对退火铝板方 盒件拉深成形过程进行了模拟，并与实验结果进行了对比分析，计算结果与试验结果显示了 较好的一致性。 关键词：晶体塑性模型，织构，动力显式有限元，Mindlin 壳单元 中图分类号：TG302
局部坐标系下的塑性应变率与整体坐标系下的塑性应变率之间可由下式转换
εˆ&ipj = Θli ε&lpkΘkj
(20) (21)
Θ 为局部坐标系与整体坐标系之间的变换张量。整体坐标系下塑性应变率 ε& p 通过晶体塑性
本构模型计算得到。
由式（19）可以得到第二类 Kirchihoff 应力增量为
∆Sˆ ij
ε&ˆkel
=
De ijkl
(εˆ&kl
− εˆ&kpl )
(19)
-3-
http://www.paper.edu.cn
∇
式中 σˆ ij —局部坐标系下 Cauchy 应力 σˆij 的 Jaumann 率张量
De ijkl
=
De ijkl
−
D D e e ij33 kl 33
/
De 3333
式中 De —四阶弹性本构张量
g (s) —滑移系上的临界切应力
m—速率敏感指数
为提高计算的稳定性，采用 PIERCE 等[11]提出的切线系数法。在时间步长∆t 内增量 ∆γ(s)
可以由 t 和 t+∆t 时刻的 γ&(s) 线性插值得到，即 ∆γ(s) = [θγ&(s)(t + ∆t) + (1− θ)γ&(s)(t)]∆t
δ N —节点 N 处的厚度 vi3N —节点 N 处中面的单位法矢量。
图 1 四节点 Mindlin 曲壳单元示意图
设
v
1N i
和
v i2 N
分别为与
v
3N i
相互正交的单位矢量，假定壳体的中面法线变形之后仍保持
为直线，但不一定再垂直于变形后的壳中面，则单元内任一点的位移可由各中面节点的三个
线位移分量
xi' = Λij x j
(1)
式中 Λij—包含欧拉角(ϕ1,φ,ϕ2 )坐标转换张量[8]
两坐标系下塑性应变率之间的转换关系为
ε&ipj' = Λik ε&kpl Λjl
(2)
晶体的变形梯度可通过乘法分解表示为
F = Fe ⋅Fp
(3)
式中 FP—对应于滑移面上的塑性剪切变形
Fe—对应于弹性变形与晶格的转动
图 1 所示四节点 Mindlin 曲壳单元由上下两个 曲面及周边以壳体厚度方向的直线为母线的曲面 所围成。设 xi(i=1,2,3）为宏观整体坐标系，单元坐 标插值分别为
∑ ∑ xi
=
4
ψ N xiN
N =1
+
4 N =1
ζ 2
ψ
N
δ
N
v
3N i
(15)
式中 ψ N —四节点线性单元的形函数
η
ζ
ξ
-5-
http://www.paper.edu.cn
(a) 计算结果
(b) 试验结果
图 5 盒形件成形后 ODF 计算结果及试验对比
-6-
http://www.paper.edu.cn
参考文献
[1] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Royal Society of London Proceedings, 1948, A: 193-281 [2] Gotoh M. A finite element analysis of the rigid plastic deformation of the flange in a deep drawing process based on a fourth-degree yield function, International Journal of Mechanical Science, 1980, 22: 367-377 [3] Barlat F, Lian J. Plastic behavior and stretchability of sheet metals. Part I: a yield function for orthotropic sheet under plane stress conditions. International Journal of Plasticity, 1989, 51～66 [4] Anand L, Balasubramanian S. Polycrystal Plasticity: Application to Earing in Cup Drawing, Annals of the CIRP, 1996, Vol. 45(1):263-268 [5] Xie C L, Nakamachi E. Investigations of the formability of BCC steel sheets by using crystalline plasticity finite element analysis. Materials and Design, 2002, 23: 59-68 [6] Nakamachi E, Xie C L. Harimato M., Drawability assessment of Bcc steel by using elastic/crystalline viscoplastic finite element analysis. Int. J. Mech. Sci., 2001, 43: 631-652 [7] Nakamachi E, Hiraiwa K, Morimoto H, Harimoto M. Elastic/crystalline viscoplastic finite element analyses of single- and poly-crystal sheet deformations and their experimental verification, Int. J. of Plasticity, 2000, 16:1419-1441 [8] Nakamachi E, Dong X. Study of texture effect on sheet failure in a limit dome height test by using elastic/crystalline viscoplastic finite element analysis, J. Appl. Mech. Trans. ASME(E), 1997, 64: 519-524 [9] Nakamachi E, Huo T. Dynamic-explicit elastic plastic finite-element simulation of hemispherical punch-drawing of sheet metal. Engineering Computations, 1996, 13: 327-338 [10] Pierce D, Asaro R J and Needleman A. Material rate dependent and localized deformation in crystalline solids[J]. Acta Metall., 1983, 31:1951 ~1972 [11] Pierce D, Shih C F and Needleman A. A tangent modulus method for rate dependent solids[J]. Computers & Structures, 1984, 18: 875-887.