零件和涂层热膨胀匹配仿真研究报告

摘  要

结构在受热条件下会产生热膨胀效应。由于不同材料之间的热膨胀程度不同，导致不同材料的结构变形不一致，使得结构内部及不同材料交界面之间产生应力。当这种应力超过了材料的许用值时，就会产生局部破坏而失效。因此，在设计此类设备和产品的过程中，要充分考虑结构的热膨胀效应，从而避免失效的发生。

传统的设计手段是通过理论计算及经验公式的方法进行，可以处理简单结构及简单材料特性。但是对于复杂结构形式及需要局部精确计算时往往力不从心，面对参数设计时更是无从下手。因此参数化有限元分析及多场耦合技术解决此类问题的优势逐渐突显出来。本文采用ANSYS软件有限元分析方法及WB多场耦合技术及参数化分析手段建立了仿真的有限元模型，通过多场耦合技术计算了结构的热应力及交界面切向应力，同时利用参数优化分析技术确定了结构安全条件下的材料参数范围

1.项目计算内容及目的

涂料是现代制造业中大量采用的工业原料。其主要功能除了美观之外，还具有其他很多特殊的作用，例如隔热，防氧化，屏蔽电磁干扰等。因此对于许多重要的产品来说，所使用的涂料往往很大程度上决定了产品的质量和使用寿命。

电子设备作为技术含量极高的工业产品，对于涂料的性能指标和质量的优劣要求日趋严格。涂料质量的优劣除了体现在其自身产品的好坏之外，还体现在其与基体材料性质的匹配程度上。在热条件下的力学性能是其中的一个重要方面。在某种外界条件下当两者之间的热膨胀程度不匹配时，就会发生交界面的剪切应力超过粘接强度而发生涂料脱落的现象，从而失去作用。

本项目的主要目的在于找出零件与涂层之间的热匹配数据，通过改变涂料的热膨胀系数找到粘接面上剪切应力随膨胀系数变化的响应曲线，从而确定在安全条件下涂层热膨胀系数范围。同时研究在指定的温度载荷及边界条件下，零件和涂层的热传导过程及粘接面剪切应力的变化，从而深入了解整个系统变化过程。

为了达到上述目的，本项目中主要的计算内容包括：

Ø 瞬态温度场分析：计算25分钟内结构温度的变化过程

Ø 瞬态结构分析：计算25分钟内，温度变化导致的结构力学特征变化，粘接面剪切应力变化

Ø 耦合场分析：利用ANSYS多物理耦合场计算结构热应力

Ø 参数优化及响应面分析：利用ANSYS DX高级参数优化仿真模块计算各参数之间的函数关系，确定涂料设计参数的范围

2. 计算理论基础及失效评定准则

2.1. 热分析理论基础

有限单元法是将连续模型离散化后，通过数值计算模拟真实物理世界的数值分析方法。借助软件技术及计算机技术，可以求解多种物理场。在本次分析中主要是针对温度场及应力场进行分析。

热分析主要分为稳态热分析及瞬态热分析：

2.1.1. 稳态热分析

对于一个稳态热分析，温度矩阵{T}通过下面的矩阵方程解得：

式中：是传导矩阵，包含导热系数、对流系数及辐射率和形状系数

是节点温度向量；是节点热流率向量

在稳态分析中，不考虑时间因素的影响

如果系统的净热流量为0，流入系统的热量加上系统自身产生的热量等于流出系统的热量：，则系统处于热稳定状态。

2.1.2. 瞬态热分析

瞬态传热过程指一个系统的加热或冷却过程，其中系统的温度、热流率、热边界条件及系统内能随时间都有明显的变化。根据能量守恒原理，瞬态热平衡可以表达为：

其中是比热矩阵，考虑系统内能的增加；是温度对时间的导数

ANSYS Mechanical求解是基于热传导方程，即傅里叶定律。

2.2. 失效评定准则

通过热应力耦合仿真得到零件与涂层之间的热应力之后，可以与已知的粘接强度进行比较。ANSYS中可以输出粘接面剪切应力随时间变化的曲线，通过读取曲线上的最大值，以此来判断整个过程中是否会达到极限粘接强度。当超过极限粘接强度时，就会发生涂料脱落的现象，判定为涂层失效。

3. 结构有限元模型

3.1. 几何模型

本项目中模型如图所示，涂料覆盖在上表面，底面有用于固定的螺栓孔，孔径3.5mm。

图1：模型几何结构

3.2. 有限元模型建立及细节处理

几何模型在有限元分析之前需要进行网格划分，即进行结构离散。因此在生成分析模型时需要进行单元类型的选择、边界条件的等效、涂料与零件之间的接触控制等工作，以下将进行详尽叙述。

3.2.1. 单元类型

单元划分及类型选择是有限元分析的基础。几何模型离散就是将连续的求解域离散为网格状的分块区域的**体。网格化后每一个小格所包含的区域称为一个单元，有限个单元组成的网格被称为有限元模型。因此单元是组成有限元模型的基本元件，单元具有不同的类型以模拟不同的计算环境。

本项目中所使用的单元包括结构单元和热单元，结构单元使用SOLID185单元，用于计算结构的应力分布；热分析单元采用SOLID70单元，用于计算温度分布。粘接面采用ANSYS独触处理技术，接触面单元采用CONTA174，TARGE170单元计算粘接面应力。模型如下：

图2：有限元网格模型

3.2.2. 接触控制

实际过程中涂料与零件之间在未脱离前是不存在相对位移的，因此可以认为两者之间是固联的不存在相对滑动。在ANSYS中可以设置为绑定接触，同时为了能够在后处理中提取出粘接面的切向应力，边界协调方程采用增广拉格朗日协调方程。程序在处理过程中自动在零件与涂层之间形成接触对，单元分别为CONTA174，TARGE170单元，两者之间满足边界协调性。

接触控制中的关键参数除接触协调方程外，在热分析中需要接触面之间的导热系数。热量通过涂层向零件传递或从零件向涂层传递过程中，由于粘接面的结合程度不同存在热量的损失，不是完美的热传递过程。

图4；接触热阻

这种热阻效应会受到诸如接触压力、表面粗糙度、材料导热特性等方面的影响，因此在瞬态热传递过程中会影响温度分布。ANSYS中通过一个热导系数进行计算： 

其中：TCC——接触导热系数

      Ttarget——目标面温度

      Tcontact——接触面温度

本项目中由于没有确切的参数，因此作为一个参数研究其对温度变化的影响。

3.2.3. 边界条件处理

a) 热分析边界条件

零件由常温（20℃）上升至600℃，温度上升100℃所需时间为1分钟，此后保持25分钟。因此零件与环境的热交换主要是通过对流方式实现。ANSYS中可以直接模拟空气对流过程，所需要的参数为对流系数及环境温度；

分析中建立的对流方程为： 

q——对流热量

h——对流系数，计算中取常用的空气对流系数5e-6W/mm2℃

Tsuface——零件对流表面的温度

Tambient——环境温度（随时间变化）

本项目中环境温度随时间变化，可以通过ANSYS中提供的表数据加载，温度-时间梯度为1.67℃/s

图5：表数据

室温取20℃，按照梯度计算，348s后达到600℃，一直保持到25分钟后（1200s）。ANSYS可以在每个瞬态时间步内自动更新环境温度从而修正对流方程。

因此，整个热边界条件为零件、涂层与空气的对流，外界环境向结构中输入热量。

图6：对流边界条件

b) 结构分析边界条件

结构分析方面主要是零件的固定方式，本项目中主要的约束形式包括：螺栓孔固定约束，安装底面约束、螺栓孔对于零件安装法向约束。

ANSYS中提供了多种约束方式，本项目中主要涉及的有两种：位移约束及受压面约束。对应关系如下：

受压面约束——螺栓孔固定约束，安装面约束。受压面约束的特点是当两个面之间是相互挤压的关系时，约束其位移，当相互分离时，受压面自由移动不受约束。

位移约束——零件安装法向约束，约束零件沿安装面法向位移，防止其与安装面分离。 

图7：结构分析边界条件

由于在瞬态过程中，约束面的受压，分离状态是实时变化的，因此需要ANSYS进行非线性仿真修正约束状态，ANSYS能够自动修正边界条件进行计算，智能化程度极高。

3.3. 计算单位制

采用如下单位制进行计算：长度单位-毫米（mm），质量单位-千克（kg），时间单位-秒（s），压力单位-兆帕（MPa），力的单位-牛（N）温度单位-摄氏度（℃）

3.4. 材料数据

本项目中零件材料分别为石英陶瓷，常温下涂层和石英陶瓷材料特性如下：

4. 瞬态热结构耦合分析

本项目中，涂层和零件在受热条件下，由于两者之间的热膨胀系数不同导致粘接面产生剪切应力；同时，受到零件支撑方式的影响，会产生内应力。因此结构实际是受到两个不同物理场的综合作用。ANSYS需要依次求解温度场和应力场，将温度场计算的瞬态结果作为应力场的输入条件，结合材料的热膨胀系数计算热应力和热应变，从而得到涂层和零件之间的接触应力。多物理场作为ANSYS最为领先的分析技术对于解决此类问题提供了全面的解决方案。

4.1. ANSYS多场耦合分析方法

ANSYS中多场耦合分析方法包括顺序耦合和直接耦合。顺序耦合是使用不同物理场对应的单元建立一套分析流程，将前一个物理场分析的结果作为后一个物理场分析的输入条件进行计算，在ANSYS中可以搭建这样的流程，如下图所示： 

图8：多场耦合分析流程

直接耦合通过ANSYS开发的耦合单元，在求解方程中同时包含多个物理场的自由度，例如耦合单元solid226可以一次计算同时求解温度和应力，完成热结构耦合分析。这种方法快速直接，但耗费计算机资源较大，对网格要求较高。本项目中，采取顺序耦合的方法，先计算温度的瞬态变化，再计算热应力分布。

4.2. 瞬态热结构耦合分析

4.2.1. 瞬态热分析结果

计算边界条件参数如下：

Ø 接触热导系数TCC=50 W/mm2℃

Ø 空气对流系数h=5e-6 W/mm2℃

Ø 环境温度从室温20℃开始，上升梯度为1.667℃/s，348秒后到达600℃保持到25分钟

结构在50s，150s，348s，600s温度分布如图10所示。

图10：温度分布

最大温度点变化曲线如图11所示。

图11：温度-时间曲线

从图中可以看出，在环境温度上升过程中，整个结构温度变化与时间成正比曲线。达到相同温度的时间，最大温度比环境温度滞后4秒左右。从图中可以看出，当348秒环境温度到达600℃时，整个结构已经接近环境温度，此后的时间整个结构保持不变，温度维持在600℃不变，可认为结构温升达到稳态。

4.2.2. 结构分析

将上述所有时刻点的温度分布结果输入到结构分析中即可进行热-结构耦合分析。分析结果如下选取50s，150s，348s，600s四个时刻的等效应力分布如图12所示。

 图12：应力分布

选取50s，150s，348s，600s四个时刻的粘接面剪切应力分布如图13所示。

 图13：剪切应力分布

粘接面最大剪切应力随时间变化曲线如图14所示。

图14：粘接面剪切应力变化

从以上分析结果可以看到，粘接面的剪切应力在整个时间过程中呈现线性变化，当环境温度到达600℃后，结构呈现稳态，应力分布不在变化。此时，结构达到最大应力状态，同时粘接面剪切应力达到最大值，涂层最容易出现开裂脱落现象。

4.3. 不同接触热阻条件下对比分析

由于接触热导是未知参数，因此本节将接触热阻作为研究对象，分析接触热阻对温度分布的影响。选取两组参数进行对比分析。

选取接触热导与空气对流系数一致5e-6 W/mm2℃，计算温度随时间的变化曲线：

 图15：温度分布

通过与4.2.1节中，接触热导为TCC=50 W/mm2℃相比，最高温度、最低温度相差一度，差别较小。同时温度随时间变化曲线如图16所示。

图16：应力变化曲线

可以看到，348秒后，环境温度达到600℃后，结构温升达到稳态，与环境温度保持一致。因此可以判断，涂层与零件之间的热量传递对各自的温度变化影响极小，结构温度的变化主要由外界空气对流热量决定。

5. 涂层热膨胀系数参数化研究及结果分析

此项分析中需要确定安全条件下，涂层的热膨胀系数范围，因此需要将涂层的热膨胀系数参数化，确定其与粘接面剪切应力之间的函数关系。找出当粘接面剪切应力在1MPa（粘接强度）时，热膨胀系数的取值范围。

5.1. ANSYS参数化分析方法

ANSYS参数化方法有很多，常用的包括ANSYS Workbench参数化管理平台及ANSYS  APDL参数化建模语言。

5.1.1. ANSYS Workbench参数化管理平台

ANSYS Workbench参数化平台是ANSYS各种分析功能的高度集成化的数据管理平台。平台可以集成ANSYS旗下所有的分析功能，分析模块之间可以迅速实现数据连接及共享。因此可以实现快速的多物理场耦合仿真。同时，每个分析流程的相关参数可以通过一个统一的平台进行管理，包括几何设计参数、材料参数、载荷参数、分析控制参数等。本项目中所涉及的所有参数均可通过此平台进行管理。

5.1.2. ANSYS  APDL参数化设置

ANSYS  APDL是ANSYS结构分析的参数化设计语言。建模及分析的所有过程都可以通过参数控制，本项目中需要控制的参数可以通过“=”实现，命令流用文本编辑器即可，部分命令流如下：

!************建立参数*************

thickness=5e-3  !隔热层厚度

width=100e-3

length=120e-3

h1=31.5e-3

h2=40.3e-3

rad=6e-3

dens1=410      !涂层密度

dens2=1650     !石英陶瓷密度

alpx1=1.46e-6   !涂层热膨胀系数

alpx2=1.5e-6    !石英陶瓷热膨胀系数

c1=0.94        !涂层比热

c2=0.91        !石英陶瓷比热

ex1=5.3e8      !涂层弹模

ex2=1.8e10     !石英陶瓷弹模

!************定义材料参数*************

/prep7

mp,ex,1,ex1       !调用参数建立材料模型

mp,prxy,1,prxy

mp,dens,1,dens1

mp,kxx,1,kxx1

mp,c,1,c1

mp,alpx,1,alpx1

通过“=”赋值符号定义的变量，可以在后续的APDL语句中进行调用，变量名代表的数值可以用于定义几何参数，材料参数，结果参数等。

5.2. 瞬态热结构耦合参数优化研究

本项目中需要建立涂层热膨胀系数与粘接面剪切力之间的函数关系，利用ANSYS DesignXplorer参数优化分析分析工具，可以建立任意多个输入参数与输出参数之间的响应关系。

本项目中首先将涂层热膨胀系数设置为输入参数：

图17：输入参数设置

将整个导热过程中粘接面最大剪切应力作为输出参数。设定涂层热膨胀系数的范围为1×10-8~4.5×10-6，程序提供了几种样本点的确定方法，通常采用中心复合差分法确定在设定的区域内选取哪些样本点。程序需要通过这些样本点计算输出参数的取值。

图18：样本点结果

得到这些样本点后，ANSYS DesignXplorer响应面分析可以得到输入参数和输出参数之间的函数曲线，采用完全二阶多项式法拟合函数关系：

图19：输入输出参数曲线

5.3. 分析结论

从上图中可以看到，当涂层的热膨胀系数小于3.7×10-6/℃时，粘接面的最大剪切应力小于1MPa，属于安全范围。同时，由于零件材料的热膨胀系数为1.5×10-6/℃，因此当涂层的热膨胀系数接近零件材料时，粘接面剪切应力较小。越远离则应力越大，两者的匹配性越差，在使用过程中越容易出现开裂失效。

6. 结束语

本次项目技术内容烦琐复杂，时间周期要求短。在整个项目工作过程中，不可避免地涉及到大量的技术交流与讨论工作。通力凯顿公司成都分公司的技术团队对模型细节、计算结果进行了仔细校验和讨论，保证了该项目工作的顺利完成，在此一并致谢！

7. 参考资料

【1】王勖成 有限单元法

【2】尚晓江 ANSYS结构有限元高级分析方法与范例应用（第2版）