聚碳酸酯挤压形变力学行为有限元模拟

王博伦，王 韬*，霍钟祺，孙琦伟，李茂源，张 云，陈宇宏，颜 悦*

（1 中国航发北京航空材料研究院，北京100095；2 北京市先进运载系统结构透明件工程技术研究中心，北京100095；3 华中科技大学 材料成形与模具技术国家重点实验室，武汉430074）

随着航空材料技术的发展，聚碳酸酯（polycarbonate，PC）因其轻质、高韧性等优点，加之良好的力学与光学性能，成为飞机座舱透明件的重要选材之一［1-2］。在飞机结构中，座舱透明件的边缘需要与金属骨架相连接以保证整件的坚固稳定，因为材料匹配性和后期机械加工等潜在影响，该边缘连接区域成为透明件整体结构的薄弱位置［3］，这种边缘区域的“软弱”现象可能演变为失效破坏，增大飞行事故风险，是透明件服役中不容忽视的问题。因此，需要对强度低于航空有机玻璃、加工应力敏感的聚碳酸酯透明件在边缘区域进行补强处理，提高材料的拉伸、冲击等力学性能，以适应服役过程中动态/准静态承载的要求，保证飞行安全。

聚合物自增强是一种低成本、高效率，且对透明件边缘补强适用性好的力学强化手段，该方法的原理在于：利用聚碳酸酯自身易受外力作用而使分子链显著取向的结构及性能特点，对固态聚碳酸酯在一定应变率下进行挤压，使其产生塑性形变，进而引起分子链的有序取向，有效提升材料力学性能。Bouaksa 等［4］针对聚碳酸酯等通道角挤压过程进行了实验与模拟研究，2D-WAXS 表征结果显示产生大形变的材料内部存在一定程度的分子链取向。本课题组［5］开展了聚碳酸酯挤压形变的力学强化效果研究，结果显示形变后材料取向度明显增大，拉伸屈服强度、缺口冲击强度等力学性能显著提高。

事实上，基于自增强的有效性，认识聚碳酸酯受挤压产生大形变时的力学特性，可为聚碳酸酯形变实验的工艺选用，实现材料力学强化的目标提供重要理论支撑。近年以来，关于聚碳酸酯的形变力学行为受到广泛关注，研究者们基于不同角度提出了一些较为完善的本构模型，例如Wang 等［6］在不可逆热力学和连续损伤力学（CDM）的基础上，提出了一种基于损伤的黏弹-塑性本构模型，并将模型应用于到有限元分析中，在低于0.2 s-1的应变率下模拟了聚碳酸酯单轴压缩大变形过程，观察到模型可在一定应变率和温度范围内描述材料的损伤演变；Miehe 等［7］提出了一个基于塑性度量理论的无定形聚合物本构模型，发现在温度23~75 ℃以及应变率低于0.1 s-1的范围内，该模型可准确描述聚碳酸酯的压缩与拉伸大变形力学响应。值得注意的是，用以表征黏弹-塑性变形行为的DSGZ模型［8］因可描述聚碳酸酯在多个温度与应变率下发生形变时的非线性弹性、屈服、应变软化以及应变硬化过程而受到关注，Dar 等［9］验证了DSGZ 模型适用于表征聚碳酸酯在高应变率不同温度下的压缩力学行为，李玉梅［10］使用改进的DSGZ 模型预测聚碳酸酯的低温和中温拉伸实验，结果表明实验曲线和模型曲线的一致性较好。

因此，为认识材料大形变力学行为，并进一步模拟聚碳酸酯受挤压后弹-塑性形变过程，推动聚合物自增强在聚碳酸酯透明件力学强化中的应用，本工作关注目前研究较少且与应用需求相符合的中应变率形变（由于低应变率下慢速、低效，高应变率下易造成材料破坏，中应变率代表有应用价值、且容易实现的形变速度范围），首先基于聚碳酸酯中应变率压缩实验结果，选择DSGZ 本构模型，依据材料的应变率和温度响应进行数值拟合，再以DSGZ 为材料属性的输入条件，在ABAQUS 有限元软件中模拟单轴压缩变形过程，分析多物理场特征，验证模型的预测效果；第二，利用有限元模拟进一步研究聚碳酸酯平板制件的挤压形变，模拟受限挤压后材料沿一维长度方向产生固态流动的弹-塑性形变过程，预测制件不同时刻的变形状态，分析形变过程中应力、应变随温度与应变率的变化规律，为材料大形变的工程实现提供理论支持。

1 实验材料与方法

1.1 单轴压缩实验

实验原材料为TEIJIN 公司的聚碳酸酯，牌号L1250，密度为1200 kg/m3（ASTM D792），熔融指数为8.7 g/10 min（ASTM D1238，300 ℃/11.76 N）。首先对粒料进行干燥除湿处理，然后在KraussMaffei 公司的CX130-750 注射成型机上制备尺寸为ϕ20 mm×30 mm 的圆柱形压缩试样，采用INSTRON VHS 160高速液压伺服实验机完成中应变率条件下聚碳酸酯的单轴压缩实验。

1.2 本构模型参数计算

DSGZ 模型是由Duan 等［8］在Johnson-Cook 模型、G’Sell-Jonas 模型、Matsuoka 模型和Brooks 模型的基础上提出的用于描述无定形聚合物变形过程的唯象本构模型，模型表达式如（1）~（5）所示。模型中含有8 个材料参数：K，C1，C2，C3，C4，a，m，α，这些参数无实际物理意义，均需通过1.1 节中所述单轴压缩实验得出的应力-应变曲线确定。其中，T为温度，ε为应变，为应变率，f(ε)反映材料发生大变形时的应变硬化过程，g(˙，T)是h(˙，T)的无量纲形式，j(，T)为形变前期应力-应变的函数关系。

为了计算模型中各参数，首先将聚碳酸酯发生大变形时的模型近似处理：

简化后得：

选择两条相同应变率、不同温度（T1，T2）下的压缩应力应变曲线，在同一大应变的位置分别选择点（ε0，σ1）和点（ε0，σ2），代入式（6）后对应力作比值，可得：

简化后得：

为了计算C3的值，选用压缩实验中一条应力-应变曲线，读出屈服应力所对应的应变εb，可以求得：

最后，在获得参数m，a和C3的基础上，根据压缩应力-应变实验数据进行粒子群优化算法（particle swarm optimization，PSO）拟合，得到K，C1，C2，C4和α。依据上述参数求解方法，聚碳酸酯压缩过程的DSGZ 模型参数拟合结果如表1 所示。

表1 DSGZ 模型拟合参数值Table1 Fitted parameter values of DSGZ model

1.3 有限元建模

1.3.1 本构模型算法实现

ABAQUS/Explicit 是ABAQUS 软件中的显式分析模块，需要的时间增量步很小，可对短时间发生的动力学事件给予准确解答［11-12］，因此，本工作采用这一模块完成有限元模拟的求解过程。结合DSGZ 唯象模型以及Dunne 等［13］阐述的材料形变算法基础理论，使用一种唯象模型应力-应变更新算法［14-15］，在此基础上，利用用户自定义子程序功能，编写唯象模型的子程序，实现ABAQUS 对压缩过程中应力、应变的迭代计算，完成仿真模拟。

1.3.2 挤压形变模型构建

运行ABAQUS/Explicit，针对单轴压缩与平板受限挤压两种过程，分别创建一个底面直径20 mm、高30 mm 的可变形圆柱体以及一个长100 mm、宽60 mm、高15 mm 的可变形立方体。在材料属性中设定聚碳酸酯的密度为1200 kg/m3，泊松比为0.417，并将表1 中DSGZ 本构模型各参数输入用户材料界面列表之中。为实现圆柱在厚度方向的单轴压缩形变，在载荷模块中对其上表面施加恒定加载速度，下表面设置为完全固定；为实现平板在挤压形变时，制件因侧面受限而仅在一维长度方向发生变形伸长的效果，在自定义坐标（X，Y，Z轴分别位于立方体的宽度、高度、长度方向）的基础上，对平板模型边界的约束条件设置如下：立方体的上表面设置Y轴负方向的恒定加载速度，下表面在Y轴方向固定，两个100 mm×15 mm的侧面在X轴方向固定。两模型的网格划分方式均选用三维八节点减缩积分单元（C3D8R），近似全局尺寸为2，圆柱与平板的划分网格单元总数分别为1740，720，如图1 所示。最后，将针对DSGZ 模型所编写好的UMAT 子程序代码（. FOR）文件载入ABAQUS中，排除流程设置中的错误，运行程序以完成有限元模拟。挤压形变温度与应变率设定如下：单轴压缩温度选用0，25，50，70 ℃，应变率选用0.05，0.5，5，50 s-1；平板受限挤压温度选用20，50，90，120，150，160 ℃，应变率选用0.007，0.03，0.3，0.7，1.3 s-1。

图1 圆柱模型（a）和平板模型（b）的约束条件与网格划分示意图Fig.1 Schematic diagram of constraint conditions and meshing of cylindrical （a） and plate model（b）

2 结果与分析

2.1 单轴压缩形变模拟

通过中应变率下不同温度的有限元仿真，获得了一系列聚碳酸酯单轴压缩形变的力学模拟结果，图2显示单轴压缩后的模拟试样已发生显著变形，模拟试样与实验试样的尺寸基本相符，受网格精度影响，二者轮廓形状略有差别。

图2 单轴压缩后的实验试样（a）与模拟试样（b）Fig.2 Experimental sample（a） and simulated sample（b）after uniaxial compression

为了获知模拟结果的有效性，提取并整理了聚碳酸酯在温度从0~70 ℃，应变率从0.05~50 s-1下的单轴压缩应力-应变曲线的实验与模拟结果，如图3 所示。可以看到，DSGZ 本构模型可有效描述聚碳酸酯的大变形压缩力学特征对温度和应变率的相关性，并可完整模拟变形全过程，即弹性形变、屈服、应变软化以及应变硬化的应力-应变关系。对比图3（a）~（d），可以发现0 ℃和25 ℃下模拟与实验结果的曲线重合度较高，而50 ℃和70 ℃下模型对压缩屈服强度与应变软化部分的模拟精度有所下降。总体来看，模型对聚碳酸酯中应变率压缩大变形的模拟效果较好。

图3 不同温度下各应变率单轴压缩应力-应变曲线模拟与实验值（a）0 ℃；（b）25 ℃；（c）50 ℃；（d）70 ℃Fig.3 Experimental and simulated true stress-strain curves of PC compressed in various strain rates（a）0 ℃；（b）25 ℃；（c）50 ℃；（d）70 ℃

为进一步分析DSGZ 本构模型在有限元模拟中预测压缩力学行为的准确性，本工作选用统计学方法以量化聚碳酸酯单轴压缩完整弹-塑性形变过程的模拟精度，利用公式（12），（13）分别对相关系数R与均方根误差R1进行计算［16］：

式中：Xi和Yi分别表示实验值和预测值；和分别为各实验值与各预测值的平均数；n为数据总个数。相关系数代表两组数据的相关程度，而均方根误差反映两组数据的离散程度，因此，R值越高或R1值越低时，实验与预测值的相关程度越大。表2 中列举了不同温度下，中应变率单轴压缩形变有限元模拟与实验值的R与R1值，可以看出，16 组相关系数R均在0.96以上，其中，0 ℃和25 ℃的R值超过0.98，略高于50 ℃和70 ℃下的结果，对应于图3 中模拟与实验曲线较高的重合度，表明理论计算与实验数据的高度相关。同时，除了0 ℃，0.5 s-1组以外，其他组的R1值均小于5，表明了模拟结果的准确性较高。综合比较各组结果，DSGZ 本构模型对聚碳酸酯单轴压缩形变力学行为有较高的预测精度，其中25 ℃，0.05 s-1和25 ℃，0.5 s-1下的预测效果最佳。

表2 聚碳酸酯单轴压缩模拟精度统计Table2 Simulation accuracy of PC uniaxial compression

基于模型的准确性，本工作进一步关注了聚碳酸酯在单轴压缩过程中不同时刻的形变状态。图4 显示了聚碳酸酯在25 ℃，0.05 s-1下单轴压缩时，圆柱各区域形变量随时间的变化。由于模拟设定总压缩行程为15 mm，总压缩时间为10 s，因此图4（c）中可见受压缩的顶面最终在厚度方向产生15 mm 位移，而底面位移接近于零，形变量从顶面到底面逐渐减小，同时，图4（b）反映压缩过程中圆柱形变沿厚度方向的梯度变化均匀，表明整个制件经历了完全、充分的单轴压缩。

图4 不同时刻圆柱单轴压缩厚度方向变形量分布图（a）0 s；（b）5.01 s；（c）10 sFig.4 Distribution diagram of deformation degree of cylinder in thickness direction under uniaxial compression at different moments（a）0 s；（b）5.01 s；（c）10 s

在观察材料形变的基础上，本工作提取了25 ℃，0.05 s-1下聚碳酸酯单轴压缩应变随时间的变化结果，如图5 所示。由图可见，虽然在3 s 与10 s 时，制件受压缩顶面附近区域的应变值略小于其他区域，即压缩过程中圆柱在水平方向出现了一定程度的扭转与滑移，各时间圆柱不同区域的应变大小一致程度较高，DSGZ 模型可以较好地反映材料单轴压缩时的弹-塑性形变过程。

图5 不同时刻单轴压缩圆柱的应变分布图 （a）0 s；（b）3 s；（c）6 s；（d）10 sFig.5 Strain distribution diagram of cylinder under uniaxial compression at different moments （a）0 s；（b）3 s；（c）6 s；（d）10 s

聚碳酸酯单轴压缩时，自身力学行为对温度及应变率的敏感性也是本工作关注的一项重要内容。整理图3 中应变率为0.5 s-1时，不同温度下圆柱单轴压缩应力-应变模拟曲线如图6 所示，可以看出，随着温度升高，材料的弹性模量降低，屈服强度以及后屈服阶段的压缩应力均显著下降。图6（a）中屈服强度从0 ℃时的86.9 MPa 降至70 ℃时的58 MPa，降幅达到33.2%。随着温度升高，应力下降幅度有所减小，模拟过程有效反映了实验中所体现的单轴压缩温度敏感性。图6（b）展示了圆柱在25 ℃，各应变率下模拟单轴压缩时的应力-应变曲线，可以看出，在同一温度下，应变率从50 s-1下降至0.05 s-1时，材料弹性模量、屈服强度以及屈服点后的应力均有降低，屈服强度从84.9 MPa 降至68.6 MPa，降幅19.2%。对比图6（a），（b）可见，模拟结果充分展现了聚碳酸酯单轴压缩响应温度比应变率更敏感的这一特点。

图6 0.5 s-1（a）与25 ℃（b）下的单轴压缩应力-应变模拟曲线Fig.6 Simulated stress-strain curves of cylinder under uniaxial compression at 0.5 s-1（a） and 25 ℃（b）

2.2 平板受限挤压形变模拟

基于聚碳酸酯圆柱试样单轴压缩模拟结果，利用准确度较好的DSGZ 模型模拟平板受限挤压形变，重点关注平板多物理场的分布与变化规律。

首先，图7 展示了90 ℃，0.03 s-1下挤压形变后，平板各个区域在自身厚度与长度方向上的形变量分布。模拟中设定完整挤压行程为5 mm，因此在图7（a）中，受到挤压后的顶面Y1最终在厚度方向上发生5 mm 的位移，同时底面Y2的位移接近于零，形变量从Y1到Y2逐渐减小。图7（b）中，平板长度方向的中点Z0在该方向上的形变量基本为零，而平板的两端Z1和Z2均已产生约22.6 mm 的形变，且形变量从Z1到Z0、Z2到Z0均逐渐减小，有效体现了聚碳酸酯平板在宽度受限的状态下受到单轴挤压后，以中心位置为基准，沿一维长度方向向两端扩展的流动方式。

图7 挤压形变后平板变形量分布图 （a）厚度方向；（b）长度方向Fig.7 Distribution diagrams of deformation degree in plates after press-induced deformation （a）thickness direction；（b）length direction

进一步观察图8 所示平板在90 ℃，0.03 s-1下挤压形变过程中不同时刻的弹、塑性总应变分布，可以明显看出，尽管平板与圆柱均经历从弹性到塑性形变的多个阶段，但在侧边受限的情况下，顶面受到挤压后的平板各位置所产生的应变大小的一致性较圆柱更好，未出现局部偏差较大的塑性变形区域，印证了DSGZ 模型的可靠性。宏观来看，相比单轴压缩试样，较大体积的聚碳酸酯制件在承受外力后也可产生均匀的大变形。

图8 挤压形变不同时刻平板的应变分布图 （a）0 s；（b）2.56 s；（c）5.03 s；（d）7.55 s；（e）10 sFig.8 Strain distribution diagrams of plates at different moments of press-induced deformation （a）0 s；（b）2.56 s；（c）5.03 s；（d）7.55 s；（e）10 s

基于宏观现象，本工作着眼于局部区域，考察平板发生大变形后不同位置的力学状态随时间变化的差异，在模型中心截面上选取了7 个特征单元，以分析90 ℃，0.03 s-1下形变过程中制件不同位置的应力与应变，如图9 所示。其中1~5 号特征单元位于平板长度方向，6 号、3 号、7 号特征单元位于厚度方向。提取7 个特征单元的力学数据，绘制出图10 所示各特征单元的应变-时间曲线与应力-应变曲线。从图10（a），（b）中可以明显看到，各个特征单元的应变与时间均近似于线性关系，制件应变随时间延长匀速增加。1~5 号特征单元的应变-时间曲线十分贴近，6 号、3 号和7号特征单元的曲线同样表现为高度重合。可以认为，挤压形变过程中，平板在长度方向与厚度方向上的不同区域应变的变化规律具有一致性。同时，图10（c），（d）中分别显示的1~5 号以及6 号、3 号和7 号特征单元的应力-应变曲线，均可完整体现聚碳酸酯的弹-塑性形变特征，且各曲线基本重合，有力证明了整个制件变形的充分与均匀。

图9 平板特征单元位置示意图Fig.9 Schematic diagram of characteristic unit positions in a plate

图10 90 ℃，0.03 s-1形变时平板特征点的应变-时间曲线（a），（b）与应力-应变曲线（c），（d）Fig.10 Strain-time curves（a），（b） and stress-strain curves（c），（d） of the characteristic points in plates deformed at 90 ℃ and 0.03 s-1

聚碳酸酯平板的挤压形变力学行为同样表现出温度与应变率敏感性。从图11（a）所示应变率0.3 s-1下，各温度挤压形变的制件应力-应变曲线模拟值可以看出，各曲线均完整反馈弹-塑性变形特征，弹性模量、屈服强度随温度升高而显著降低，屈服强度从20 ℃下的73.4 MPa 降至160 ℃下的41.6 MPa，降幅达到43.3%。随着温度升高，屈服强度的下降逐步减缓，150 ℃和160 ℃下的曲线几乎重合。由此可见，与前文所述圆柱相近，模拟过程有效反映了聚碳酸酯平板挤压形变的温度敏感性。图11（b）展示了平板在20 ℃，各应变率下模拟挤压形变时的应力-应变曲线。不难看出，同一温度下，应变率在0.007~1.3 s-1范围中变化时，材料弹性模量与屈服强度变化幅度较小，随着形变速度降低，屈服强度仅从最高78.3 MPa 降至最低66 MPa，降幅15.7%，这一结果同样与圆柱单轴压缩模拟的率相关性一致。因此，聚碳酸酯平板挤压形变模拟体现了材料响应温度变化比形变速度变化更敏感的这一特征，可为较大尺寸制件大形变力学行为提供理论预测。

图11 0.3 s-1（a）与20 ℃（b）下的平板应力-应变模拟曲线Fig.11 Simulated stress-strain curves of plates at 0.3 s-1（a） and 20 ℃（b）

3 结论

（1）基于DSGZ 唯象本构模型，对聚碳酸酯中应变率单轴压缩形变过程进行了有限元模拟，模拟与实验应力-应变曲线吻合度高。量化计算结果进一步表明，DSGZ 模型对聚碳酸酯压缩弹-塑性形变力学行为有较高的预测精度，可成为聚碳酸酯挤压形变有限元模拟中材料属性的有效输入。

（2）基于模型预测的准确性，分析聚碳酸酯圆柱单轴压缩形变过程的多物理场分布，圆柱形变沿厚度方向的梯度变化均匀，不同区域应变大小一致度较高，圆柱整体经历了完全、充分的单轴压缩。相比应变率变化，温度变化使材料力学行为出现更为敏感的响应。

（3）聚碳酸酯平板受限挤压形变模拟结果显示，不同时刻平板制件各位置应变一致，整体变形均匀，应力响应温度变化比响应形变速度变化更为敏感。