1. 引言

PEEK材料自从研制出来，一直被作为一种重要的战略性国防军工材料，许多国家限制出口[1]。从上世纪八十年代的“七五”计划开始，国家特别重视高性能聚合物的研究与开发，为有关研究投入了大量的人力和物力。在完全使用国产原料的情况下，吉林大学特种工程塑料研究中心成功研制出聚醚醚酮树脂的实验小样，使中国成为了继英国之后的第二个系统掌握聚醚醚酮树脂制备技术的国家[2]-[4]。文周等[5]采用正交试验法，研究了FDM成型技术中打印速度、打印温度、填充比等打印参数对PEEK试样抗拉性能的影响。在拉伸实验的基础上，通过有限元分析软件ABAQUS对PEEK材料紧凑拉伸试样进行了结构特性评价，得出了最优打印条件：打印速度为40 mm/s，打印温度为350℃，填充比为20%。赵丙辉等[6]研究了运用3D打印制造技术加工PEEK材料骨科植入样件的力学性能及氨基改性对样件成骨性能的影响。S. Berretta等[7]通过FDM技术制得了PEEK碳纳米聚合物，分析了碳纳米在PEEK基体中的分布以及不同喷头温度下PEEK碳纳米聚合物中的纤维质量对FDM成型PEEK零件的力学性能影响。本文主要研究在FDM成型过程中，喷头温度、打印速度和打印层厚对温度场的影响各有不同，本文主要讨论喷头温度的影响；通过有限元的方法模拟仿真打印过程，分析模型不同时刻的温度场分布、不同节点的温度–时间变化曲线和温度梯度的变化趋势，得出不同打印条件的影响规律，用以指导实际生产生活，对产品的加工有着重要的意义。

2. FDM成型传热模型的建立

在FDM成型过程中，PEEK材料在密闭的打印腔中逐层堆叠在基板平台上，虽然成型腔温度能维持在80℃，基板平台温度能维持在150℃，但是对于熔融态的PEEK (350℃左右)，它们之间存在一定程度的温度差，因此，热能的传递不仅发生在材料的内部，在材料表面与周围介质之间也存在。根据传热学理论，热量的传递有三种基本形式，即热传导、热对流和热辐射。

2.1. 热传导

在固体、液体和气体中均有热传导的发生。在FDM成型过程中，熔融态的PEEK将热量传递到基板平面上。热传导遵循傅利叶定律[8]：

$q^{″}=-k\frac{\text{d}T}{\text{d}x}$ (1)

式中， $q^{″}$ 为热流密度(W*m⁻²)， $k$ 为导热系数(W/m∙℃)，负号表示热量向低温方向流动。

2.2. 热对流

热对流现象仅仅存在于流体中。由于PEEK材料从熔融态转变为固态的时间短暂，因此不考虑熔融态PEEK的对流换热，本文重点考虑已成型表面与周围空气之间的热对流。

热对流可以用牛顿冷却方程来表示：

$q^{″}=-h\left(T_S-T_B\right)$ (2)

式中， $q^{″}$ 为热流密度(W*m⁻²)， $h$ 为对流换热系数， $T_S$ 为固体表面温度(℃)， $T_B$ 为流体温度(℃)。

2.3. 热辐射

热辐射与其他两种传热方式的区别在于，热传导和热对流均需要介质才能传递热量，而热辐射不需要任何媒介就可直接向外传递。因此，在进行有限元模型的建立时必须考虑热辐射的边界条件。

热辐射可以用斯蒂芬–波尔兹曼方程来表示：

$q=\epsilon \sigma T^4$ (3)

式中， $q$ 为物体热辐射能量密度， $\epsilon$ 为辐射率， $\sigma$ 为斯蒂芬–波尔兹曼常数 $\sigma =5.67\times {10}^{-8}{\text{W/m}}^2{\text{k}}^4$ ， $T$ 为绝对温度(℃)。

3. FDM成型热力场数值模拟模型的建立

3.1. 单元类型

Figure 1. Scan trajectories

图1. 扫描轨迹

在有限元计算中，模型的尺寸和划分网格的大小对仿真分析的计算精度和效率有着很大的影响，为了提高分析精度和计算效率，选择合适的模型尺寸和网格大小非常关键。在实际加工过程中，如果样件尺寸较大，导致加工时间较长，根据有关文献的查阅和实际加工经验，以长20 mm，宽3 mm，高5 mm的立方体作为加工尺寸。由于实际打印过程中，扫描轨迹如图1所示，则需要将仿真模型的宽划分为7份。考虑计算效率的前提下，将仿真模型的长划分为40份，而仿真模型的高则需要根据不同实验打印层厚的大小进行划分。

ANSYS 19.0中可用于非线性瞬态温度场分析的单元类型有Solid 70和Solid 90，其中Solid 70为8节点，Solid 90为20节点，由Solid 70高阶演变而来。在一定的误差允许范围内，为了提高计算效率，减少计算时间，选择Solid 70单元进行建模，结果如图2。

Figure 2. Schematic diagram of the unit division

图2. 单元划分示意图

3.2. PEEK材料属性

在FDM成型过程中，PEEK材料的密度会随着结晶程度的不同而发生变化，这使得其材料热物理特性表现出随温度的变化而发生改变的特性，并且与温度的关系是非线性的。根据目前现有的文献资料表明，PEEK材料的热物理性能参数并不齐全，在熔融温度及以上的参数是空白，因此本研究采用了一定的简化。表1为温度场分析需要的PEEK材料各项物理性能参数表。

Table 1. The thermophysical performance parameters of PEEK

表1. PEEK热物理性能参数

3.3. 初始条件和边界条件

在FDM成型过程中，PEEK材料会经历两次相变：从固态到熔融态，再到固态。整个成型过程的理论分析非常复杂，其中有扫描轨迹、打印条件和PEEK材料性能等多种因素的共同影响，全面考虑多种因素作用下的成型过程是非常困难的。因此，本研究需要在一定的假设基础上建立简化模型，假设如下[9]：

(1) 在FDM成型过程中，PEEK材料在挤出时的温度为喷头温度，且因为层厚较小，忽略时间的影响，认为材料落在基板上的温度与喷头温度一致。

(2) 打印过程中，已成型的部分与空气之间的热量传递方式为热对流与热辐射，在有限元分析中需要考虑综合换热系数。

(3) 忽略PEEK材料不同结晶率对密度的影响，采用成型温度区间的平局密度值计算。

(4) 采用ANSYS分析方法中的生死单元法，认为单元格的厚度为材料的分层厚度，每个单元格的激活时间与实际扫描填充速度有关。

在打印之前，需要对打印腔进行足够长时间的预热以达到预先设定的温度，假设经过预热后打印腔内部的温度一致，初始条件为：

${T\left(x,y,z,t\right)|}_{t=0}=T_0$ (4)

FDM打印过程中，PEEK材料逐层堆积，层与层之间存在热量的传递，其主要形式是热传递。同时，打印样件与空气之间存在温度差，有热量交换现象。在本研究中，PEEK材料的熔点为334℃，而打印腔的温度为80℃，整个打印过程均在高温环境中进行。当温度较高时，热量传递的主要形式为热辐射，因此，样件表面的散热方式主要为热对流和热辐射的综合作用，根据资料，其折算换热系数 $h$ 为[10]：

$h=\sigma \epsilon \left(T_p^2+T_c^2\right)\left(T_p+T_c\right)+h_c$ (5)

式中， $h_c$ 为自然对流换热系数(W/m²∙℃)， $\sigma$ 为斯蒂芬–波尔兹曼常数 $\sigma =5.67\times {10}^{-8}{\text{W/m}}^2{\text{k}}^4$ ， $\epsilon$ 为样件材料黑体温度， $T_p$ 为样件温度， $T_c$ 为环境温度。

在实际打印过程中，样件在基板平面上成型，基板平面始终保持较高温度，在打印结束后快速将样件从基板平面上剥离下来，然后在空气中冷却。打印过程中，样件与基板平面之间热量传递的方式主要为热传导，打印结束后主要为热对流和热辐射。仿真分析中，采用上述中样件表面的换热系数来模拟整个过程与基板接触面的热量传递。

3.4. 生死单元法

在有限元仿真模拟过程中，模型单元无法在求解过程中建立，因此需要在建模时创建所有实体单元。使用生死单元法可以通过“EKILL”命令将某一时刻不存在的单元杀死，使其不参与计算，如果需要该单元时则使用“EALIVE”命令将其激活。杀死单元的方法是通过在其刚度矩阵上乘以一个非常小的系数，使单元的各种物理特性在仿真分析的过程中约等于零，对整个计算结果不产生影响，此系数在ANSYS软件中的默认值为1.0 × 10⁻⁶，也可以根据实际需要设置为其他非零值，而激活单元只需要将单元的属性值恢复到之前的设定值[11]。

仿真开始之前，在ANSYS软件中建立模型，并划分单元，并将所有单元全部杀死，然后按照预设的扫描轨迹激活相应的单元，然后施加温度载荷，从而模拟材料逐层堆积的过程，过程如图3所示。

Figure 3. Schematic diagram of the life and death unit method

图3. 生死单元法示意图

3.5. 计算流程及简化

APDL (ANSYS Parametric Design Language)参数化编程语言是ANSYS中通过命令流的方式实现变量参数化并且建立分析模型求解的语言。通过APDL的方式，可以实现智能化的有限元分析，对于需要多次重复分析的复杂模型，可以通过修改已经定位好的参数化变量，实现不同模型的分析。

本文编写了一套APDL命令流，采用有限元分析方法中的生死单元法，模拟不同打印条件下样件温度场和热应力耦合场，实现打印过程中材料按照预设的扫描轨迹堆积成型，在不同时刻激活相应的单元的过程。在实际打印过程中，喷头连续加工每一层单元，有四次扫描轨迹。在有限元仿真模拟过程中，本研究采用简化加载，将每一层按照打印轨迹分为四次加载，在误差允许的范围内，极大地减少了计算量，大大提高了计算效率。有限元仿真计算流程如图4所示。

Figure 4. Calculation flow of temperature field simulation

图4. 温度场仿真计算流程

4. 温度场仿真结果分析

在FDM成型中，模型的温度场分布情况与喷头温度、打印速度和打印层厚等打印条件都有一定的关系。本文主要讨论喷头温度的影响；为准确地分析不同喷头温度对样件温度场分布情况的影响规律，得到影响样件变形的主要因素，为实际生产加工提供理论基础，本研究对喷头温度进行了实验，如表2所示，并分别对每个实验组编写APDL命令流，进行温度场有限元分析，通过ANSYS软件对仿真进行后处理，得到温度场分布和节点的温度变化，并对这些数据进行分析。

Table 2. Process parameters of FDM molding at different nozzle temperatures

表2. 不同喷头温度下FDM成型工艺参数

4.1. 温度场分布

(1) 打印过程中不同时刻温度场分布

通过ANSYS软件后处理功能，可以查看不同时刻模型的温度场求解结果。选用FDM成型中常用的打印条件：喷头温度350℃，打印速度30 mm/s，打印层厚0.2 mm。打印时间为112.75 s，冷却时间为200 s。图5为模型在打印过程中各个时间段温度场的分布情况。

(a) t = 14.1 s时的温度场 (b) t = 46.89 s时的温度场

(c) t = 75.26 s时的温度场 (d) t = 94.71 s时的温度场

(e) t = 112.75 s打印完成时的温度场 (f) t = 162.75 s冷却50 s时的温度场

(g) t = 212.75 s冷却100 s时的温度场 (h) t = 312.75 s冷却200 s时的温度场

Figure 5. Temperature field distribution at different printing times

图5. 不同打印时刻温度场分布情况

(a) T = 340℃ (b) T = 350℃

(c) T = 360℃ (d) T = 370℃

Figure 6. Temperature field distribution during model cooling for 200 s at different nozzle temperatures

图6. 不同喷头温度下模型冷却200 s时温度场分布

打印开始时，未成型部分的单元处于杀死的状态，其温度与打印腔温度一致，待加工区域喷头按轨迹扫描加工，其温度瞬间增高，已成型部分通过热传导、热对流和热辐射的方式进行换热，使温度逐渐降低。随着打印进行，热影响区域不断扩大，当前打印层的等温线分布密集，下层等温线较为稀疏，同一时刻模型的温度场分布不均匀，温度变化大。打印完成后模型的温度下降较快，最后趋于环境温度。

(2) 喷头温度对温度场分布的影响

PEEK材料的熔点为334℃，熔融区间的温度范围较大，当PEEK为熔融态时可以从喷头顺利挤出，图6为340℃~370℃下模型不同时刻冷却200 s温度场的分布情况。打印完成时，不同喷头温度下模型的温度场分布规律基本相同，都是上层温度高，越往下温度越低，底层温度最低。在不同冷却时间，模型中心温度高，四周温度低，模型整体温度随着冷却时间降低。随着喷头温度升高，打印不同时刻的最高温度升高，高温范围扩大，模型最低温度也有所提高，但是整体变化较小。喷头温度的变化对打印完成时模型的整体温度场分布影响较小。

4.2. 节点温度变化分析

分析不同时刻模型节点温度随时间变化情况，对分析温度场的变化规律有着重要的作用。选用常用打印条件作为节点温度变化分析模型，喷头温度T = 350℃，打印速度V = 30 mm/s，打印层厚H = 0.2 mm。所取节点位置如图7所示，图8为相应节点的温度–时间变化曲线。

Figure 7. Schematic representation of the model sampling nodes

图7. 模型取样节点示意图

Figure 8. Model the temperature-time change curves of different nodes

图8. 模型不同节点的温度–时间变化曲线

各个节点的温度均是在短暂的时间里到达最高值，斜率较大；到达最高温度之后，由于热量传递的原因，节点的温度开始下降，但是下降的速率相比较上升时较小，最后趋于平缓，直至变为水平线。其中，在温度下降的过程中，由于材料在不断的堆积，刚成型的部分温度高于已成型部分，会有热量传递下来，所以曲线呈现波动下降。在打印结束后，模型处于空气中，其换热系数增大，温度曲线下降更加迅速，之后曲线趋于平滑。FDM成型是一个高度瞬态加热过程，在加热过程和冷却过程都具有温度骤变的特点，但是冷却过程的斜率变化比加热过程缓慢，并不是突然变化的，具有一定斜率梯度的平滑下降。

4.3. 温度梯度分析

温度梯度是有限元温度场分析中非常重要的因素，是指在特定区域环境内，描述温度变化最迅速的方向以及温度是以何种速率变化的物理量。分析温度梯度对研究样件热变形具有重要意义。选取不同打印条件下，模型节点4的温度梯度进行分析，图9为不同喷头温度下节点4的温度梯度变化曲线。结果表明喷头温度的变化对节点温度梯度影响较小，随着喷头温度的降低，温度梯度值降低。根据仿真结果分析，在完成加工要求的前提下，降低喷头温度，可以降低温度梯度。

Figure 9. Temperature gradient change curve of the characteristic point at different nozzle temperatures

图9. 不同喷头温度下特征点温度梯度变化曲线

5. 结论

本文基于瞬态非线性热分析理论，在ANSYS平台上运用生死单元法编写APDL命令，按照实际打印轨迹简化加载，模拟打印过程中温度场的变化。通过分析不同打印条件下温度场的分布情况、各节点温度随时间变化规律及其冷却速率，得出了以下结论：

(1) 同一时刻，正处于打印层的温度场等温线分布密集，已打印完成的下层区域等温线分布较为稀疏，并且模型整体温差较大，温度分布不均匀。

(2) 对模型不同打印条件下的温度场进行分析，结果表明：升高喷头温度，模型在打印完成时的整体温度升高，温度分布更加均匀。

(3) 对不同节点的温度–时间变化曲线进行分析，结果表明：节点温度随时间快速变化，能够在短时间内达到峰值，然后随时间波动下降，最后趋于平缓，直至温度冷却到打印腔温度。

(4) 对模型不同喷头温度下节点温度梯度变化曲线进行分析，结果表明：降低喷头温度，可以降低温度梯度。

参考文献