国产6061-T6铝合金的修正塑性本构模型

摘要

铝合金材料已广泛地应用于各种工程结构中。能准确描述铝合金在复杂应力状态下的塑性本构关系，对铝合金结构在复杂工况下的分析和设计具有重要意义。传统的von Mises屈服准则忽略了应力状态对金属塑性的影响，并不适用于复杂应力状态下的铝合金材料。通过试验证明了应力三轴度和Lode角均会对6061-T6铝合金的塑性行为产生影响，并在此基础上提出了适用于国产6061-T6铝合金在复杂应力状态下的塑性本构模型，并校准了全过程硬化规律和新屈服准则。通过试验和数值结果对比可知，该塑性本构模型可以精确地模拟国产6061-T6铝合金在复杂应力状态下的弹塑性行为。

关键词

铝合金; 塑性本构模型; 应力三轴度; Lode角

铝合金材料已经被广泛应用于体育馆、天文馆等重要的公共建筑中。在实际工程中，大多数关键构件通常处于复杂应力状态下，如铝合金单层网壳中最常用的板式节点。此外，铝结构连续性倒塌过程包含了材料在各种复杂应力状态下的大变形和断裂失效，采用精细的材料本构模型对模拟结构在局部损伤出现后的响应至关重要，这也是结构鲁棒性研究的发展趋势^［1］。因此，有必要充分研究结构用铝合金在各种应力状态下的塑性本构模型，为铝合金结构的分析和应用提供基础。

基于von Mises屈服准则的经典J₂塑性模型以其简洁性和适用性被广泛应用于描述金属材料的屈服行为，但其在多轴复杂应力状态下的局限性也逐渐被认识到^［2-7］。目前学界普遍认识到金属材料的屈服行为受到应力三轴度或Lode角的影响，或与二者均有关。应力三轴度和Lode角分别反映了静水压力和应力偏量第三不变量的影响。早在20世纪80年代，Spitzig和Richmond^［2］就发现铁基材料和铝合金的塑性行为受到静水压力的影响。Zhang^［3］给出了6组屈服函数和塑性流动势的表达式，能考虑静水压力和应力偏量第三不变量的影响，并将其应用于3种不同金属（铝合金、不锈钢、锆合金），认为5083铝合金的屈服准则与应力三轴度无关。Bai和Wierzbicki^［4］提出了同时考虑应力三轴度和Lode角效应的屈服准则，能较好地适用于2024-T351铝合金，但存在高应力三轴度时屈服函数不合理的问题。由此可见，不同研究对于金属材料（如铝合金）塑性行为的影响因素给出了不同的结论。

对于常用的结构钢，Kong等^［5］和Wang等^［6］针对国产普通低碳钢和高强钢进行了一系列试验，表明仅考虑Lode角效应的塑性模型足够精确地模拟结构钢的弹塑性行为。Cortese等^［7］对部分欧标钢材的研究也得到了相同的结论。然而，对于国产结构用铝合金，尚需通过更多的试验探究其塑性本构的影响因素及影响规律。

本文对常用的6061-T6挤压铝合金进行了不同应力状态下的拉伸或压缩试验，证明了经典J₂塑性本构模型在复杂应力状态下的局限性，且应力三轴度和Lode角均会影响铝合金的塑性行为。基于此，本文提出了适用于国产6061-T6铝合金材料的塑性本构模型，包括全过程硬化规律表达式和新的屈服准则，并通过有限元模拟结合试验结果验证了该塑性模型的准确性。

1 试验概况

1.1　试验设计

在三维应力空间中，任意一点P的应力状态可以由3个主应力（σ₁，σ₂，σ₃，即笛卡尔坐标）来表示，也可由一组柱坐标（r，θ，l）来表示，如图1所示。经数学推导，可得到柱面坐标系中r、θ、l 与3个应力不变量之间的关系，进而得到应力三轴度和Lode角的表达式。应力三轴度被定义为平均应力和von Mises应力（等效应力）的比值，Lode角则是柱坐标系中的坐标θ（0≤θ≤π/3），如式（1）—（2）所示。

图1  笛卡尔坐标系和柱面坐标系下的应力状态表征

Fig. 1  Characterization of stress states in Cartesian and cylindrical coordinate systems

$$T=\frac{I_{\mathrm{1}}}{\mathrm{3}q}=\frac{{\sigma }_{\mathrm{m}}}{q}$$

（1）

$$\theta =\frac{\mathrm{1}}{\mathrm{3}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{27}{J}_{\mathrm{3}}}{\mathrm{2}{q}^{\mathrm{3}}}\right)$$

（2）

式中：I₁为应力张量的第一不变量；J₃为偏应力张量的第三不变量；σ_m为平均应力（I₁/3）；q为von Mises应力，也称为等效应力。

本文研究对象为常用的国产6061-T6铝合金。为考虑不同应力状态的影响，进行了光滑圆棒（SRB）、缺口圆棒（NRB）、开槽平板（FGP）、纯剪试件（PS）、压缩圆柱（CS）和缺口压缩圆柱（NCS）的单调拉伸或压缩试验。图2给出了所有试件的几何形状和尺寸，同一试件数量为3个。表1给出了各试件的理论应力三轴度和Lode角，计算方法参考文献［8］。

图2  试验件的几何形状及尺寸

Fig. 2  Geometries and dimensions of tested specimens

表1  各试件的理论应力三轴度和Lode角

Tab. 1  Theoretical stress triaxiality and Lode angle of specimens

试件类型	试件编号	T	θ	试件类型	试件编号	T	θ
光滑圆棒	SRB	1/3	0	开槽平板	FGP1（R=5 mm）	0.713	π/6
缺口圆棒	NRB1（R=2.5 mm）	1.314	0		FGP2（R=10 mm）	0.647	π/6
	NRB2（R=5 mm）	0.907	0	压缩圆柱	CS	-1/3	π/3
	NRB3（R=15 mm）	0.551	0	缺口压缩圆柱	NCS	-0.551	π/3
纯剪试件	PS	0	π/6	-	-	-	-

本试验在同济大学工程结构性能演化与控制重点实验室完成。试验加载照片如图3所示。试验采用位移加载，拉伸试验的加载速率控制在0.1~0.5 mm·min^-1，采用标距为50 mm的高精引伸计测量试件的变形；压缩试验的加载速率控制在0.2~1.0 mm·min^-1，采用位移计测量试件变形，标距为试件高度（30 mm），上下加载端各安装2个位移计以取其测量值的平均值。

图3  试验加载装置

Fig. 3  Test setup

1.2　试验现象

图4给出了部分代表性试件的破坏模式。由缺口圆棒和开槽平板试件的试验结果可知，应力三轴度越大，试件延性越差，断裂位移越小。纯剪试件的破坏模式是剪切断裂，裂纹贯通于整个临界截面。对于压缩试件，仅在个别试件中发现了与加载方向呈45°的斜裂纹；但铝合金常用于薄壁构件，通常在出现如此大的压缩应变之前就已发生失稳。

图4  试件破坏模式

Fig. 4  Failure modes of specimens

2 全过程真实应力‒真实塑性应变关系

2.1　全过程真实应力-真实塑性应变关系表达式

对于需考虑大塑性应变情况的数值分析，精确的全过程硬化规律即真实应力‒真实塑性应变关系曲线是必不可少的。常用模型有Swift模型、Voce模型、Ramberg-Osgood模型等^［9］，其表达式如式（3）：

$$\mathrm{S}\mathrm{w}\mathrm{i}\mathrm{f}\mathrm{t}{\sigma }_{\mathrm{t}}=K_{\mathrm{s}}{\left({\epsilon }_{\mathrm{s}}+{\epsilon }_{\mathrm{t}\mathrm{p}}\right)}^{m_{\mathrm{s}}}$$

（3a）

$$\mathrm{V}\mathrm{o}\mathrm{c}\mathrm{e}{\sigma }_{\mathrm{t}}=A_{\mathrm{v}}-K_{\mathrm{v}}\cdot e^{-m_{\mathrm{v}}{\epsilon }_{\mathrm{t}\mathrm{p}}}$$

（3b）

$$\mathrm{R}\mathrm{a}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{g}-\mathrm{O}\mathrm{s}\mathrm{g}\mathrm{o}\mathrm{o}\mathrm{d}\frac{{\epsilon }_{\mathrm{t}}}{{\epsilon }_{\mathrm{y}}}=\frac{{\sigma }_{\mathrm{t}}}{{\sigma }_{\mathrm{y}}}+B_{\mathrm{r}}{\left(\frac{{\sigma }_{\mathrm{t}}}{{\sigma }_{\mathrm{y}}}\right)}^{m_{\mathrm{r}}}$$

（3c）

式中：σ_t为真实应力；ε_t为真实应变；ε_tp为真实塑性应变；σ_y为屈服应力；ε_y为屈服应变；K_s、ε_s、m_s、A_v、K_v、m_v、B_r和m_r为根据试验数据校准得到的参数。

材料在单轴拉伸条件下的工程应力‒应变关系可通过光滑圆棒单调拉伸试验获得，本文通过以下方法将其转化为真实应力‒塑性应变关系。

在颈缩发生前，临界截面上的应力和应变均匀分布，采用式（4）实现从工程应力应变到真实应力应变的变换^［10］。

$$\begin{array}{l}{\epsilon }_{\mathrm{t}}=\mathrm{l}\mathrm{n}(\mathrm{1}+{\epsilon }_{\mathrm{e}})\\ {\sigma }_{\mathrm{t}}={\sigma }_{\mathrm{e}}(\mathrm{1}+{\epsilon }_{\mathrm{e}})\\ {\epsilon }_{\mathrm{t}\mathrm{p}}={\epsilon }_{\mathrm{t}}-{\sigma }_{\mathrm{t}}/E\end{array}$$

（4）

式中：σ_e为工程应力；ε_e为工程应变；E为弹性模量。

在实际应用式（4）时，可通过拟合试验数据得到具体的表达式。通常，单阶段曲线难以准确描述铝合金颈缩前的真实应力‒塑性应变曲线。因此，本文采用双阶段曲线进行拟合，即式（5）。通过比对拟合曲线和试验数据，确定分段点所对应的真实塑性应变为ε_tp0=0.2 %。

$${\sigma }_{\mathrm{t}}=\left\{\begin{array}{ll}{{\sigma }_{\mathrm{t},}}_{\mathrm{0.2}}{\left(\mathrm{1}+\frac{{\epsilon }_{\mathrm{t}\mathrm{p}}-{\epsilon }_{\mathrm{t}\mathrm{p}\mathrm{0}}}{{\epsilon }_{\mathrm{0}}}\right)}^{n_{\mathrm{0}}}& {\epsilon }_{\mathrm{t}\mathrm{p}}\le {\epsilon }_{\mathrm{t}\mathrm{p}\mathrm{0}}\\ {{\sigma }_{\mathrm{t},}}_{\mathrm{0.2}}{\left(\mathrm{1}+\frac{{\epsilon }_{\mathrm{t}\mathrm{p}}-{\epsilon }_{\mathrm{t}\mathrm{p}\mathrm{0}}}{{\epsilon }_{\mathrm{1}}}\right)}^{n_{\mathrm{1}}}& {\epsilon }_{\mathrm{t}\mathrm{p}\mathrm{0}}<{\epsilon }_{\mathrm{t}\mathrm{p}}\le {\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}\end{array}\right.$$

（5）

$${\epsilon }_{\mathrm{1}}=\left({\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}-{\epsilon }_{\mathrm{t}\mathrm{p}\mathrm{0}}\right)/\left[{\left({\sigma }_{\mathrm{t},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}/{\sigma }_{\mathrm{t},\mathrm{0.2}}\right)}^{\mathrm{1}/n_{\mathrm{1}}}-\mathrm{1}\right]$$

（6）

式（5）—（6）中：ε₀、n₀、n₁为根据试验曲线拟合的参数；σ_t，0.2为ε_tp等于0.2%时对应的真实应力。为保证曲线的连续性，ε₁由式（6）计算，其中σ_t，neck和ε_tp，neck分别为颈缩起始点的真实应力和真实塑性应变。

颈缩发生后，由于临界截面上的应力和应变分布极不均匀，不能由式（4）获得真实应力应变。通常单个模型难以精确描述颈缩后材料的力学行为，本文采用组合Swift-Voce模型，引入权重因子w对Swift和Voce模型组合，得到半经验的修正模型。综上所述，铝合金全过程真实应力‒真实塑性应变关系的表达式可由式（5）—（7）表述。

$$\begin{array}{cc}{\sigma }_{\mathrm{t}}=w{K}_{\mathrm{s}}{\left({\epsilon }_{\mathrm{s}}+{\epsilon }_{\mathrm{t}\mathrm{p}}\right)}^{m_{\mathrm{s}}}+\left(\mathrm{1}-w\right)\left(A_{\mathrm{v}}-K_{\mathrm{v}}\cdot {\mathrm{e}}^{-m_{\mathrm{v}}{\epsilon }_{\mathrm{t}\mathrm{p}}}\right)& {\epsilon }_{\mathrm{t}\mathrm{p}}>{\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}\end{array}$$

（7）

通常将真实应力‒应变曲线中斜率等于当前应力的点视为颈缩起始点^［9］；由于颈缩点附近材料的塑性已充分发展，弹性应变增量可忽略不计。故颈缩起始点可由式（8）确定：

$$\frac{\mathrm{d}{\sigma }_{\mathrm{t}}}{\mathrm{d}{\epsilon }_{\mathrm{t}}}{\left|\right.}_{\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{g}}={\sigma }_{\mathrm{t},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}\Rightarrow \frac{\mathrm{d}{\sigma }_{\mathrm{t}}}{\mathrm{d}{\epsilon }_{\mathrm{t}\mathrm{p}}}{\left|\right.}_{\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{g}}={\sigma }_{\mathrm{t},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}$$

（8）

为保证颈缩起始处曲线的连续性和光滑性，将颈缩起始处的真实应力和真实塑性应变代入式（3）和（8），利用颈缩前试验数据得到Swift和Voce模型的拟合参数，则颈缩后真实应力‒真实塑性应变曲线可变换为

$${\sigma }_{\mathrm{t}}\left({\epsilon }_{\mathrm{t}\mathrm{p}}\right)=\begin{array}{cc}w{\sigma }_{\mathrm{t},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}{\left(\mathrm{1}+\frac{{\epsilon }_{\mathrm{t}\mathrm{p}}-{\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}}{{\epsilon }_{\mathrm{s}}+{\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}}\right)}^{{\epsilon }_{\mathrm{s}}+{\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}}+\left(\mathrm{1}-w\right){\sigma }_{\mathrm{t},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}\left(\mathrm{1}+\frac{\mathrm{1}}{m_{\mathrm{v}}}\left(\mathrm{1}-{\mathrm{e}}^{-m_{\mathrm{v}}\left({\epsilon }_{\mathrm{t}\mathrm{p}}-{\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}\right)}\right)\right)& {\epsilon }_{\mathrm{t}\mathrm{p}}>{\epsilon }_{\mathrm{t}\mathrm{p},\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{k}}\end{array}$$

（9）

2.2　全过程真实应力-真实塑性应变关系的参数确定

2.2.1　颈缩前关系式参数确定

表2给出了6061-T6铝合金的工程力学性能参数，包括弹性模量E、名义屈服强度f_0.2、抗拉强度f_u以及抗拉强度对应的名义应变ε_u。颈缩前关系式的各参数可以根据试验数据获得。σ_t，0.2可通过名义屈服强度f_0.2由式（4）计算得到；σ_t，neck和ε_tp，neck可由式（8）确定，也可用抗拉强度f_u及其对应的名义应变ε_u采用式（4）计算得到；ε₁可由式（6）计算；ε₀、n₀、n₁可由试验数据拟合得到。各参数见表3。

表2  6061-T6铝合金的工程力学性能参数

Tab. 2  Engineering mechanical properties of 6061-T6 aluminum alloy

试件	弹性模量E /MPa	名义屈服强度f0.2 /MPa	抗拉强度fu /MPa	抗拉强度对应的名义应变εu
SRB-1	69 156	298.84	319.12	0.065 37
SRB-2	68 555	297.79	315.65	0.068 91
SRB-3	68 905	298.95	317.02	0.065 48
平均值	68 872	298.53	317.26	0.066 58

表3  6061-T6铝合金颈缩前真实应力-真实塑性应变关系的各参数

Tab. 3  Parameters of true stress-true plastic strain relationship of 6061-T6 aluminum alloy before necking

试件	σt，0.2 /MPa	σt，neck /MPa	εtp，neck	ε0	n0	n1	ε1
SRB-1	300.73	339.98	0.058 404	0.001 991	0.033 473	0.146 407	0.043 005
SRB-2	299.68	337.40	0.061 715	0.002 005	0.035 069	0.166 624	0.057 575
SRB-3	300.84	337.78	0.058 524	0.002 004	0.034 298	0.144 440	0.045 986
平均值	300.42	338.39	0.059 547	0.002 000	0.034 280	0.152 490	0.048 855

2.2.2　颈缩后关系式参数确定

（1） 采用数值迭代法确定权重因子w

2.1节给出了颈缩后真实应力‒真实塑性应变关系表达式，但权重因子w无法直接确定。本节通过数值迭代法来校准w，即对光滑圆棒试件进行有限元分析得到数值解，通过更新迭代w，直至有限元和试验结果吻合。采用ABAQUS 6.14^［11］建立模型，根据试件对称性建立1/8模型，采用C3D8R单元，如图5所示。经网格敏感性分析，确定关键区域网格尺寸为0.25 mm，其他区域网格尺寸为0.5 mm。

图5  光滑圆棒有限元模型

Fig. 5  Numerical model of smooth round bar

通过数值迭代法得到了6061-T6铝合金颈缩后的真实应力‒真实塑性应变关系，各拟合参数为w=0.25，ε_s=0.006 786，m_v=23.08。采用该曲线得到的有限元结果（组合Swift-Voce）和试验结果对比见图6，二者非常吻合，表明颈缩后曲线的有效性。图7给出了全过程真实应力‒真实塑性应变曲线。图6中还给出了颈缩后曲线采用三个常用模型的模拟结果，这三者的准确性均低于组合Swift-Voce模型。

图6  采用4种模型的数值结果和试验结果对比

Fig. 6  Comparison of experimental results and numerical results based on four models

图7  6061-T6铝合金全过程真实应力‒真实塑性应变曲线

Fig. 7  Full range true stress-true plastic strain of 6061-T6 aluminum alloy

（2） 采用瞬时截面积迭代方法验证

数值迭代法操作简单，但需要提前确定颈缩后真应力‒应变曲线形式。Ho等^［12］提出的瞬时截面积迭代法克服了该缺点，本文对其进行改进，使用塑性应变量代替总应变量。通过该方法来验证数值迭代法的准确性。其流程如图8所示。图中式（a）即为式（4），式（b）和（c）中参数说明如下：N_i为实测荷载，A_i为实测临界截面面积，A₀为初始截面面积；应变校正系数η_ε为临界截面上最大等效塑性应变ε^p_max和平均塑性应变ε^p_av的比值，应力校正系数η_σ为临界截面上最大等效应力q_max和平均纵向应力σ_av的比值。首先，采用高清摄像机获得试验过程中试件的变形^［12］，得到颈缩处的截面直径变化，假定临界截面保持为圆形，可得到实测截面面积；其次，大应变下真应力和真应变由图中式（b）得出，η_ε、η_σ的初值η_ε，0、η_σ，0设为1.0；再次，通过数值模型得到临界截面应力应变不均匀分布情况η_ε，j、η_σ，j；最后，将η_ε，j、η_σ，j回代到式（b），得到j+1次修正的真应力σ_t，j+1和真塑性应变ε_tp，j+1。重复以上步骤，直至校正系数收敛。

图8  瞬时截面积迭代法流程图

Fig. 8  Flowchart of instantaneous cross-sectional area iteration method

图9给出了光滑圆棒试件的迭代结果曲线。到第5步迭代（it5）结束，有限元荷载‒变形曲线和试验结果吻合良好。由图7可知，通过数值迭代法和瞬时截面积迭代法得到的硬化曲线整体较为一致，二者互相验证，证明本节得到的颈缩后真实应力‒塑性应变关系表达式是准确的。下文可直接采用数值迭代法得到铝合金的全过程真实应力‒塑性应变关系。

图9  瞬时截面积迭代方法模拟结果

Fig. 9  Numerical results based on instantaneous cross-sectional area iteration method

3 屈服准则的确定

3.1　von Mises屈服准则的适用性

本文建立了各类试件的有限元模型，见图10。为减少计算时间，除纯剪试件采用1/2模型外，其他试件均采用1/8模型。材料特性取自上一节确定的全过程硬化规律，采用von Mises屈服准则。通过网格敏感性分析确定网格尺寸，在保证精度的情况下提高计算效率。

图10  试件的有限元模型

Fig. 10  Numerical models of specimens

图11为部分试件试验和模拟结果的对比。由缺口圆棒结果可知，应力三轴度会明显影响铝合金的屈服强度，而开槽平板和纯剪试件则反映出Lode角对屈服强度也存在较大的影响。可见，von Mises屈服准则难以准确预测铝合金在各种应力状态下的响应。

图11  试验荷载-变形曲线和采用von Mises屈服准则的有限元结果对比

Fig. 11  Comparison of experimental load-deformation curves and numerical results using von Mises criterion

3.2　修正屈服准则

本文试验证明了应力三轴度和Lode角均对铝合金材料的塑性行为产生影响。能同时考虑二者影响且常用的屈服准则是Bai和Wierzbicki^［4］提出的屈服准则F（见式（10）），Bai准则将应力三轴度和Lode角的影响分别考虑，即式（10）中的C（T）和C（θ）。

$$F=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)C(T)C(\theta )=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)\left[\mathrm{1}-c_T\left(T-T_{\mathrm{0}}\right)\right]\left[c_{\mathrm{s}}+\left(c_{\mathrm{a}\mathrm{x}}-c_{\mathrm{s}}\right)\left(\gamma -\frac{{\gamma }^{m+\mathrm{1}}}{m+\mathrm{1}}\right)\right]$$

（10）

$$\begin{array}{c}{c}_{\mathrm{a}\mathrm{x}}^{}=\left\{\begin{array}{cc}{c}_{\mathrm{t}}^{}& \mathrm{0}\le \theta \le \mathrm{\pi }/\mathrm{6}\\ c_{\mathrm{c}}^{}& \mathrm{\pi }/\mathrm{6}<\theta \le \mathrm{\pi }/\mathrm{3}\end{array}\right.,\end{array}$$

$$\gamma =\frac{\sqrt[]{\mathrm{3}}}{\mathrm{2}-\sqrt[]{\mathrm{3}}}\left[\frac{\mathrm{1}}{\mathrm{c}\mathrm{o}\mathrm{s}(\theta -\mathrm{\pi }/\mathrm{6})}-\mathrm{1}\right]$$

（11）

式中：F为屈服函数；$\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)$为对应于等效塑性应变${\overline{\epsilon }}_{\mathrm{p}}$的等效应力$\overline{\sigma }$，即2.2节中的σ_t（ε_tp）；γ为关于Lode角的函数；$c_T$为关于应力三轴度的拟合参数；$c_{\mathrm{a}\mathrm{x}}$、$c_{\mathrm{s}}$、$c_{\mathrm{t}}$、$c_{\mathrm{c}}$和m为关于Lode角的拟合参数。

基于Bai屈服准则，本文提出了一种新的适用于铝合金的屈服准则，下面将从应力三轴度和Lode角两个角度分别进行阐述。

3.2.1　应力三轴度对屈服准则的影响

Bai准则^［4］中假定屈服函数为应力三轴度的线性函数；但当T为正时，屈服强度随T增加迅速小于0，如图12所示，这是不合理的。因此本文参考Karr等^［13］提出的指数型各向同性屈服准则（式（12）），拟克服应力三轴度为正时屈服强度迅速小于0的问题。

$$F=\sqrt[]{J_{\mathrm{2}}}+k_{\mathrm{2}}{\mathrm{e}}^{k_{\mathrm{3}}{I}_{\mathrm{1}}}-k_{\mathrm{1}}$$

（12）

式中：J₂为应力偏量第二不变量；k₁、k₂、k₃为表征材料属性的常数。

图12  Bai准则和新准则中C（T）曲线对比示意图

Fig. 12  Comparison of C(T) curve in Bai criterion and new criterion proposed

首先，将I₁量纲一化，转化为应力三轴度。将式（12）改写为等效应力q、应力三轴度T和材料参数（K₁、K₂、K₃）的表达式：

$$F=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)C(T)=q+K_{\mathrm{2}}{\mathrm{e}}^{K_{\mathrm{3}}T}-K_{\mathrm{1}}$$

（13）

其次，当T=1/3时，材料处于单轴受拉状态，屈服准则将退化为von Mises屈服准则，则可得到：

$$F=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)=q+K_{\mathrm{2}}{\mathrm{e}}^{K_{\mathrm{3}}/\mathrm{3}}-K_{\mathrm{1}}$$

（14）

将式（14）代入式（13），并将参数合并成c_i（i=1~3），整理可得C（T）的表达式：

$$C\left(T\right)=c_{\mathrm{1}}+(\mathrm{1}-c_{\mathrm{1}}){\mathrm{e}}^{K_{\mathrm{3}}\left(T-\mathrm{1}/\mathrm{3}\right)}$$

（15）

为避免屈服应力随应力三轴度增大时迅速小于零，将${\mathrm{e}}^{K_{\mathrm{3}}\left(T-\mathrm{1}/\mathrm{3}\right)}$变换为$\mathrm{1}/(\mathrm{1}+{\mathrm{e}}^{-c_{\mathrm{3}}\left(T-\mathrm{1}/\mathrm{3}\right)})$：

$$C\left(T\right)=c_{\mathrm{1}}+(\mathrm{1}-c_{\mathrm{1}})\frac{\mathrm{1}}{\mathrm{1}+{\mathrm{e}}^{-c_{\mathrm{3}}(T-\mathrm{1}/\mathrm{3})}}$$

（16）

式（16）是关于单轴拉伸应力状态中心对称的表达形式。对其进行扩展，增强其适用性，参考Bai准则引入T₀对应力三轴度T进行修正，使新准则处于标定硬化规律试件的应力状态时，可退化为von Mises屈服准则，即此时C（T）=1。

$$C\left(T\right)=c_{\mathrm{1}}-\frac{c_{\mathrm{2}}}{\mathrm{1}+{\mathrm{e}}^{-c_{\mathrm{3}}(T-T_{\mathrm{0}})}}$$

（17）

式中：c₁、c₂、c₃、T₀是需确定的参数， c₁、c₂、c₃、T₀不完全独立，可由其中3个参数得到第4个参数的值。

图12给出了Bai准则和新准则关于应力三轴度的变化规律对比的示意图，新准则避免了应力三轴度增大情况下C（T）值迅速小于0的情况，更为合理。

3.2.2　Lode角对屈服准则的影响

对于Lode角对屈服准则的影响，Bai准则^［4］中对此已进行了充分的论证。但当应力状态为单轴拉伸时， T=1/3、θ=0，则C（T）=1、C（θ）= m/（m+1），Bai准则退化为

$$F=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)C(T)C(\theta )=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)\frac{m}{m+\mathrm{1}}$$

（18）

即此时接近却不等于von Mises屈服准则。因此，本文将m取为正无穷大，则m/（m+1）会趋向于1，式（18）退化成von Mises屈服准则，则C（θ）的表达式如下所示：

$$C(\theta )=c_{\mathrm{s}}+\left(c_{\mathrm{a}\mathrm{x}}-c_{\mathrm{s}}\right)\gamma$$

（19）

3.2.3　修正屈服准则表达式

综上所述，基于现有的屈服准则形式和试验数据，本文提出了一种新的屈服准则：

$$F=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)C(T)C(\theta )=q-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)\left[c_{\mathrm{1}}-\frac{c_{\mathrm{2}}}{\mathrm{1}+{\mathrm{e}}^{-c_{\mathrm{3}}(T-T_{\mathrm{0}})}}\right]\left[c_{\mathrm{s}}+\left(c_{\mathrm{a}\mathrm{x}}-c_{\mathrm{s}}\right)\gamma \right]$$

（20）

式中：c₁、c₂、c₃、T₀、c_s、c_c是需校准的参数，同时$\sqrt[]{\mathrm{3}}/\mathrm{2}\le c_{\mathrm{s}}^{}/c_{\mathrm{a}\mathrm{x}}^{}\le \mathrm{1}$以满足屈服面的外凸性^［14］；c_ax和γ由式（11）计算。

4 流动法则的确定

对于金属材料，常用的流动法则是关联流动法则。根据第3节确定的屈服准则，其表达式如下所示：

$$\begin{array}{l}\mathrm{d}{\mathit{\epsilon }}_{}^{\mathbf{p}\mathbf{l}}=\mathrm{d}\lambda \frac{\partial G}{\partial \mathit{\sigma }}=\mathrm{d}\lambda \left\{\frac{\mathrm{3}}{\mathrm{2}q}\mathit{s}-\overline{\sigma }\left({\overline{\epsilon }}_{\mathrm{p}}\right)\frac{-c_{\mathrm{2}}{c}_{\mathrm{3}}{\mathrm{e}}^{-c_{\mathrm{3}}(T-T_{\mathrm{0}})}}{{\left(\mathrm{1}+{\mathrm{e}}^{-c_{\mathrm{3}}(T-T_{\mathrm{0}})}\right)}^{\mathrm{2}}}\left[c_{\mathrm{s}}+\left(c_{\mathrm{a}\mathrm{x}}-c_{\mathrm{s}}\right)\gamma \right]\left(\frac{\mathrm{1}}{\mathrm{3}q}\mathit{I}-\frac{\mathrm{3}T}{\mathrm{2}{q}^{\mathrm{2}}}\mathit{s}\right)\right.\\ \left.-\overline{\sigma }\left({\overline{\epsilon }}^{\mathrm{p}}\right)\left[c_{\mathrm{1}}-\frac{c_{\mathrm{2}}}{\mathrm{1}+{\mathrm{e}}^{-c_{\mathrm{3}}(T-T_{\mathrm{0}})}}\right]\left(c_{\mathrm{a}\mathrm{x}}-c_{\mathrm{s}}\right)\left(\frac{\mathrm{3}\sqrt[]{\mathrm{3}}}{\mathrm{2}-\sqrt[]{\mathrm{3}}}\right)\left(\frac{\mathrm{t}\mathrm{a}\mathrm{n}(\theta -\mathrm{\pi }/\mathrm{6})}{\mathrm{c}\mathrm{o}\mathrm{s}(\theta -\mathrm{\pi }/\mathrm{6})}\right)\frac{\mathrm{1}}{q\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{3}\theta }\left(\frac{\mathrm{1}}{\mathrm{3}}\mathit{I}+\frac{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{3}\theta }{\mathrm{2}q}\mathit{s}-\frac{\mathrm{3}}{\mathrm{2}{q}^{\mathrm{2}}}\mathit{s}\cdot \mathit{s}\right)\right\}\end{array}$$

（21）

式中：G为塑性势函数；$\mathrm{d}{\mathit{\epsilon }}_{}^{\mathbf{p}\mathbf{l}}$为塑性应变张量；σ为应力张量；s为偏应力张量；I为二阶单位张量；$\mathrm{d}\lambda$为塑性乘子，亦为等效塑性应变增量。

已有研究证明，尽管金属材料的屈服准则依赖于平均应力，但是其塑性体积变形远远小于关联流动法则所预测的，仍可忽略^［15］。由式（21）发现$\mathrm{d}{\epsilon }_{ii}^{\mathrm{p}\mathrm{l}}$不等于0，违背了塑性变形体积不变的假设；若要满足该假设，需要将$\partial T/\partial \mathit{\sigma }$中的$\mathit{I}/\mathrm{3}q$删去，又与$\mathrm{d}{\overline{\epsilon }}_{\mathrm{p}}=\mathrm{d}\lambda$相冲突，说明关联流动法则不再适用。因此，本文仍采用很多研究的普遍做法^［6，16］，采用非关联流动法则，选择von Mises屈服函数作为塑性势函数，则流动法则如式（22）所示：

$$\mathrm{d}{\mathit{\epsilon }}_{}^{\mathbf{p}\mathbf{l}}=\mathrm{d}\lambda \frac{\partial G}{\partial \mathit{\sigma }}=\mathrm{d}\lambda \frac{\mathrm{3}}{\mathrm{2}q}\mathit{s}$$

（22）

5 数值模型分析

5.1　有限元模型

本节通过有限元模型验证所提出的修正塑性本构模型，本构模型基于UMAT子程序实现，实现过程中联合使用了图形返回算法与隐式欧拉向后算法。除材料特性外，模型参数与3.1节一致，见图10。

5.2　屈服准则参数校准

由试验结果反推本构模型参数时，需采用有限元模型迭代计算。本节采用遗传算法来得到较优的反演结果。该算法的适应度函数f如下所示：

$$\mathrm{m}\mathrm{i}\mathrm{n}f=\mathrm{m}\mathrm{i}\mathrm{n}\left[{\sum _j\sum _i\left(\frac{{\widehat{F}}_{i,j}-F_{i,j}}{F_{i,j}}\right)}^{\mathrm{2}}\right]$$

（23）

式中：$F_{i,j}$为第j个试件的第i个试验值；${\widehat{F}}_{i,j}$为第j个试件的第i个数值解。

参数校准具体流程见图13。通过该方法得到了新屈服准则的各参数：c₂=567，c₃=0.019 84，T₀=-260，c_s=0.875，c_t=1.0，c_c=1.008。屈服面形状如图14所示。

图13  基于遗传算法的参数校准流程图

Fig. 13  Flowchart of genetic algorithm-based parameter calibration

图14  屈服面形状示意图

Fig. 14  Diagram of yield surface shape

5.3　有限元计算结果

本节基于试验和有限元模拟结果对新提出的塑性本构模型进行验证，结果如下。

（1） 缺口圆棒试件

缺口圆棒试件关键点（最小截面中心）的Lode角均为0，应力三轴度随缺口半径增大而减小。图15为缺口圆棒试件的荷载‒变形曲线、采用新模型的模拟结果（new model）和采用von Mises屈服准则的模拟结果（von Mises）。Lode角不变时，随着应力三轴度增加，von Mises模型结果与试验结果差异增大，说明屈服强度随应力三轴度的增大而降低。对比各曲线可知，von Mises屈服准则高估了高应力三轴度下的屈服强度，而新塑性模型则能更好预测高应力三轴度下屈服强度降低的情况。

图15  缺口圆棒试件荷载‒变形曲线对比

Fig. 15  Comparison of load-deformation curves of notched round bar specimens

（2） 开槽平板试件

图16给出了开槽平板试件的试验和有限元荷载‒变形曲线对比结果。开槽平板试件考察的是位于平面应变状态的试件中心区域，其Lode角为π/6，应力三轴度则随缺口半径增大而减小。与von Mises屈服准则相比，新模型中影响屈服强度的因素既有应力三轴度，也有Lode角，二者均会降低开槽平板试件的屈服荷载，模拟结果更精确，表明新塑性模型比von Mises准则更适用于平面应变状态。

图16  开槽平板试件荷载‒变形曲线对比

Fig. 16  Comparison of load-deformation curves of grooved flat plate specimens

（3） 纯剪试件

理论上纯剪试件的应力三轴度应为0，Lode角应为π/6，但是受限于加工精度和试验条件，只能使得应力状态表征参数接近于期望值。在纯剪应力状态下，Lode角会降低屈服强度，而应力三轴度会略微增大屈服强度。图17对比结果表明新塑性模型更适用于纯剪状态下的铝合金。

图17  纯剪试件荷载变形曲线对比

Fig. 17  Comparison of load-deformation of pure shear specimens

（4） 压缩圆柱和缺口压缩圆柱

在轴对称压缩状态下，圆柱中心的应力三轴度为负，Lode角为π/3。此时，新屈服准则会从应力三轴度和Lode角两个角度提高铝合金的屈服强度。从图18

图18  压缩圆柱和缺口压缩圆柱的荷载变形曲线对比

Fig. 18  Comparison of load-deformation curves of compression cylinders and notched compression cylinders

对比结果可以看出，新塑性模型在受压状态下具有更好的适用性和准确性。

5 结论

本文完成了国产结构用铝合金在不同应力状态下的拉伸或压缩试验，对其塑性本构模型进行了研究，主要结论如下：

（1） 基于一系列拉伸或压缩试验，得到了国产6061-T6铝合金在不同应力状态下的塑性响应。结果表明，铝合金的塑性行为强烈地依赖于其应力三轴度和Lode角， von Mises屈服准则难以准确预测多种应力状态下铝合金的塑性响应。

（2） 提出了铝合金的全过程真实应力‒真实塑性应变关系的表达式，给出了颈缩前表达式中参数的确定方法，且通过数值迭代法和瞬时截面积迭代法对颈缩后表达式进行了相互验证。采用该方法得到的硬化曲线模拟光滑圆棒拉伸试件，其数值模拟结果与试验结果高度吻合。

（3） 针对6061-T6铝合金，提出了一种考虑应力三轴度和Lode角影响的屈服准则，并使用遗传算法对屈服准则中各参数进行校准。试验数据和有限元结果对比表明，在各种复杂应力状态下，新屈服准则比von Mises屈服准则具有更高的准确性。

作者贡献声明

郭小农：论文研究内容的整体规划以及审阅修改，提供技术及经费支持。

宗绍晗：调研及整理文献，完成试验和数值计算分析，撰写论文。

张超众：参与试验和数值计算分析，审阅并修订论文。

成张佳宁：参与试验研究。