一种基于三维图像的有限差分法预测材料热导率的方法

摘要

一种基于三维图像的有限差分法预测材料热导率的方法，属于差分法预测材料热导率的方法，具体为一种利用材料的三维图像结合有限差分法预测材料热导率的方法。其基本原理为：一、利用图像分析法区分图像中不同的相或组分；二、利用计算机语言程序，将所有像素的位置和颜色信息读入计算机内存中；三、将每个像素构建为一个胞元，按照所属组分赋予其热导率；四、离散稳态导热方程，计算传热系数矩阵，构建计算方程组；五、计算方程组得出温度场；六、计算材料等效热导率。本发明采用三维模型，能够较准确的得出热导率结果；采用有限差分法，能够节省计算时间和内存消耗。全过程自动完成，对于节点数量巨大的模型非常实用。

说明书

技术领域

本发明涉一种数值模拟预测材料热导率的方法，特别是一种基于三维图像的有限差分 法预测材料热导率的方法。

背景技术

材料的性能对于其使用效能有决定性的作用，因此测试材料的性能是在使用它前必做 的一项功课。然而对于结构材料设计来说，常常要频繁的改变材料的结构以期能得到想要 的性能，这使得需要大量的测试工作。通过材料结构预测材料的性能，能够极大的缩短材 料设计的周期，因此大量应用在复合材料设计等领域。传统的预测方法一般都使用两种方 法：1)基于理想假设的理论计算和2)基于实验的经验公式，或者将两者结合起来。对 于第一种方法，实际情况往往比理想状态复杂的多，因此理论计算的结果和实际状态可能 会有较大偏差。而对于第二种方法，由于实验条件的不同，由经验公式计算的值也可能会 与真实结果有较大的差别。这使得无论是以上哪种方法都很难适用于普遍场合。

近些年随着计算机技术的发展，人们可以利用计算机程序对于简单的材料结构进行虚 拟构建，然后进行有限元等数值建模，再设置合理的边界条件和初始条件进行计算，能够 对材料的性能进行分析，这种方法即为计算机辅助设计CAD(ComputerAidedDesign)。 由于该方法是利用材料的真实结构进行建模的，因此计算结果有较好的准确率。目前，这 种方法广泛应用于复合材料结构设计领域。然而，对于一些没有周期结构的、非均质的复 合材料，例如热喷涂WC-Co涂层，由于其结构太复杂，很难进行构建，因此该方法也难 有用武之地。

基于图像的数值建模是指利用能够代表材料结构的二维或者三维数码照片，建立数值 计算模型的一种方法。该方法无需事先对材料结构进行构建，而是直接以材料的图像作为 其结构，因此能够适用于各种复杂的材料结构。在近十年来，有学者利用有限元方法对矿 物材料、建筑材料和涂层材料进行有限元建模，成功预测了这些材料的力学、热学等性能。 然而，当对象为一些较大的或者具有非常精细结构的三维图像时，由于节点和单元体数量 成指数增加，使得分析时间和所需计算机资源大大增加，这可能会失去数值模拟简便快捷 的优势。

发明内容

为了适应节点数量巨大的三维模型，改善其占用计算机资源多、计算效率低下的弊端， 保持数值模拟在材料性能预测周期上的优势，本发明提供一种基于三维图像的有限差分法 预测材料热导率的方法，能够节约计算机资源，提高运算速率。

本发明解决其技术问题所采用的技术方案是：该预测方法以数值图像的像素作为胞 元，像素的中心作为节点，传热基本方程为基础，再通过差分方程代替有限元法进行计算；

具体包括以下步骤：

(1)获取能代表材料结构的三维图像；

(2)以图像分析法将图像中不同的组分或相分割开；

(3)利用计算机语言程序，将处理过的图像信息读入内存中，以矩阵形式存储；

(4)以像素构建胞元，按照其所属组分或相，赋予热导率；

(5)计算每个节点向周边节点的热传导系数，将值储存进热传导系数矩阵；

(6)离散稳态传热方程，应用于每个节点，构建计算方程；

(7)设置边界条件和初始条件，求解整个域的温度场；

(8)换算整体热导率。

所述的步骤(2)中，分割图像中不同组分或相是依靠其不同的颜色或形貌来实现的； 具体步骤如下：

(a)对于灰度图像，采用适当的阀值将灰度图转化为与组分或相的种类相同数量的 色阶，每个色阶代表一个组分或相；

(b)对于彩色图像，采用先将图像转化为灰度图，或者直接根据RGB三原色通道， 采用阀值进行图像分割；

(c)或根据某些组分的特殊形态，将该组分所占的区域涂上一种颜色，与其他组分 或相相区别。

所述的步骤(4)中，构建胞元，不用预先建立实体模型，直接根据像素的归属和位 置信息构建胞元，节点的位置对应像素的中心。

所述的步骤(5)中，计算热传导系数，采用7点式模型，只计算胞元与共面的胞元 之间的直接传热，不考虑仅共线或仅共点的胞元之间的直接传热。

所述的步骤(6)中，构建计算方程，以有限差分的形式离散传热方程，采用直接求 解法或迭代法两种方法计算结果。

有益效果，由于采用了上述方案，本发明借助图像像素的离散特性，将像素直接构建 成六面体胞元，将材料的三维图像构建成有限差分模型，计算多组元或相结构材料的等效 热导率。与现有技术对比，本发明避免了建立实体模型再进行网格划分这两个步骤，节省 了系统资源和计算时间。

优点：采用三维模型，能够较准确的得出热导率结果；采用有限差分法，能够节省计 算时间和内存消耗。全过程自动完成，对于节点数量巨大的模型非常实用。

1)三维模型比二维模型更能接近材料内部真实传热行为，因而计算出的结果可信度 更高。

2)有限差分法采用差分方程对连续体进行离散，只计算节点的温度值，能够比有限 元法节省系统资源，极大的提高运算速度。

3)不用进行实体建模，以像素直接构建胞元，能够节省系统资源，避免不合理的网 格划分。

附图说明

图1是本发明基于图像的有限差分法预测材料热导率的流程图。

图2是本发明三维有限差分模型的胞元图。

图3是本发明等离子喷涂氧化钇稳定氧化锆涂层的三维重建图。

图4是以图3为对象的有限差分模型计算出的温度场。

具体实施方式

该预测方法：以数值图像的像素作为单元，像素的顶点作为节点，传热基本方程为基 础，再通过差分方程代替有限元法进行计算；具体操作步骤如下：

第一步，摄取能代表材料结构的三维图像；

第二步，以图像分析法将图像中不同的组分或相分割开；

第三步，利用计算机语言程序，将处理过的图像读入内存中，以矩阵形式存储：

第四步，以像素构建胞元，按照其所属组分或相，赋予热导率；

第五步，计算每个节点在直角坐标系三个方向上与相邻节点的热传导系数，将值储存 进热传导系数矩阵；

第六步，离散稳态传热方程，应用于每个节点，构建计算方程；

第七步，设置边界条件和初始条件，求解整个域的温度场；

第八步，换算整体热导率。

所述第一步，摄取的区域大小要适中，能够代表材料的结构；区域过小会造成计算结 果不准确的问题，过大会浪费系统资源，增加计算时间。

第二步，分割图像中不同组分或相是依靠其不同的颜色或形貌来实现的；具体的步骤 如下：

(1)对于灰度图像，采用适当的阀值将灰度图转化为与组分或相的种类相同数量的 色阶，每个色阶代表一个组分或相；

(2)对于彩色图像，采用先将图像转化为灰度图，或者直接根据RGB三原色通道， 采用阀值进行图像分割。

(3)根据某些组分的特殊形态，将该组分所占的区域涂上一种颜色，与其他组分或 相相互区别。

第三步中，认为组分或相是均匀体，性质同一。在处理过的图像中，同一组分或相具 有相同的颜色属性，也具有相同的热导率。孔隙和裂纹作为同一组分，其热导率由被其俘 获的气体性质决定。

第四步中，每个像素所在的域为一个传热的基本单元，称为胞元，其热导率取决于所 属的组分的材料属性。

第五步中，节点的位置为每个胞元的中心点，称为中心离散模型，节点的温度代表了 整个胞元的平均温度。

图2是三维模型的胞元和节点示意图，i、j、k代表三维矩阵中的位置角标，Kx、Ky 和Kz是直角坐标系三个方向的热传导系数；从节点(i,j,k)与相邻节点之间的热传导系数 Kx_(i,j,k)、Ky_(i,j,k)和Kz_(i,j,k)分别为：

${\mathrm{Kx}}_{(i,j,k)}=\frac{2}{\frac{1}{{\lambda }_{(i,j,k)}}+\frac{1}{{\lambda }_{(i+1,j,k)}}}---\left(1.1\right)$

${\mathrm{Ky}}_{(i,j,k)}=\frac{2}{\frac{1}{{\lambda }_{(i,j,k)}}+\frac{1}{{\lambda }_{(i,j+1,k)}}}---\left(1.2\right)$

${\mathrm{Kz}}_{(i,j,k)}=\frac{2}{\frac{1}{{\lambda }_{(i,j,k)}}+\frac{1}{{\lambda }_{(i,j,k+1)}}}---\left(1.3\right)$

第六步中，本发明使用稳态传热方程计算材料的热导率，该方程的表达式为：

$-\overrightarrow{▿}.\left(\lambda \overrightarrow{▿}T\right)=0---\left(1.4\right)$

其中λ为热导率，T为温度。在离散该传热方程时，本发明只考虑相邻像素所在域之间直 接传热，即七点式三维传热(图2)。对于三维传热，方程(1.4)的离散过程为：

$\begin{array}{c}-\overrightarrow{▿}·\left(\lambda \overrightarrow{▿}T\right)=-\frac{\partial }{\partial x}\left(\lambda \frac{\partial T}{\partial x}\right)-\frac{\partial }{\partial y}\left(\lambda \frac{\partial T}{\partial y}\right)-\frac{\partial }{\partial z}\left(\lambda \frac{\partial T}{\partial z}\right)\\ =-\frac{{\mathrm{Kx}}_{\left(i,j,k\right)}\left(T_{\left(i+1,j,k\right)}-T_{\left(i,j,k\right)}\right)-{\mathrm{Kx}}_{\left(i-1,j,k\right)}\left(T_{\left(i,j,k\right)}-T_{\left(i-1,j,k\right)}\right)}{{\mathrm{\Delta x}}^2}\\ -\frac{{\mathrm{Ky}}_{\left(i,j,k\right)}\left(T_{\left(i,j+1,k\right)}-T_{\left(i,j,k\right)}\right)-{\mathrm{Ky}}_{\left(i,j+1,k\right)}\left(T_{\left(i,j,k\right)}-T_{\left(i,j+1,k\right)}\right)}{{\mathrm{\Delta y}}^2}\\ -\frac{{\mathrm{Kz}}_{\left(i,j,k\right)}\left(T_{\left(i,j,k+1\right)}-T_{\left(i,j,k\right)}\right)-{\mathrm{Kz}}_{\left(i,j,k+1\right)}\left(T_{\left(i,j,k\right)}-T_{\left(i,j,k+1\right)}\right)}{{\mathrm{\Delta z}}^2}=0\end{array}---\left(1.5\right)$

由于像素的各边边长相等(Δx＝Δy＝Δz)，可消去边长Δx、Δy、Δz，化简后，方程(1.5)就 改写成：

Kx_(i,j,k)T_(i+1,j,k)+Kx_(i-1,j,k)T_(i-1,j,k)+Ky_(i,j,k)T_(i,j+1,k)+Ky_(i,j-1,k)T_(i,j-1,k)(1.6)

+Kz_(i,j,k)T_(i,j,k+1)+Kz_(i,j,k-1)T_(i,j,k-1)-K_(i,j,k)T_(i,j,k)＝0

其中

K_(i,j,k)＝Kx_(i,j,k)+Kx_(i-1,j,k)+Ky_(i,j,k)+Ky_(i,j-1,k)+Kz_(i,j,k)+Kz_(i,j,k-1)(1.7)

第七步中，给待计算方向的两边界设置固定的温度T₁、T₂(T₁≠T₂)，给内部节点设置 固定的温度梯度，使温度沿传热方向线性分布。其他边界设置绝热边界条件。

第七步中，所有的节点方程构成线性齐次方程组：

[K]·[T]＝[C](1.8)

其系数矩阵[K]为稀疏对角阵。

1.对于较小的域，可以采用以下四种方式求解[T]：

(1)消元法；

(2)矩阵对角化法；

(3)克莱姆法：

$T_{(i,j,k)}=\frac{\left[K_{i,j,k}\right]}{\left[K\right]}---\left(1.9\right)$

(4)则求出稀疏矩阵的逆矩阵，然后直接求解:

[T]＝[K]^-1·[C](1.10)

2.对于较大的域，以迭代算法求解，如高斯-赛德尔迭代，其三维问题的迭代形式可以写 成：

$T_{(i,j,k)}^{m+1}=\frac{{\mathrm{Kx}}_{(i-1,j,k)}{T}_{\left(i-1,j,k\right)}^{m+1}+{\mathrm{Ky}}_{(i,j-1,k)}{T}_{(i,j-1,k)}^{m+1}+{\mathrm{Kz}}_{\left(i,j,k-1\right)}{T}_{\left(i,j,k-1\right)}^{m+1}+{\mathrm{Kx}}_{\left(i,j,k\right)}{T}_{(i+1,j,k)}^m+{\mathrm{Ky}}_{\left(i,j,k\right)}{T}_{\left(i,j+1,k\right)}^m+{\mathrm{Kz}}_{(i,j,k)}{T}_{\left(i,j,k+1\right)}^m}{K_{(i,j,k)}}---\left(1.11\right)$

其中，m为迭代的次数；但高斯-赛德尔迭代收敛比较慢，因此可采用连续超松弛法，其 迭代形式为：

$T_{(i,j,k)}^{m+1}=T_{(i,j,k)}^m+\omega (T_{(i,j,k)}^{m+1,GS}-T_{(i,j,k)}^m)---\left(1.12\right)$

其中指的是由高斯-赛德尔法，即公式(1.11)，计算的第m+1次迭代后节点(i,j,k)的温 度值；ω为松弛因子，其取值范围为(1,2)；ω的值会影响该方法的收敛速度，使迭代收敛 速度最快的ω_opt值可取：

${\omega }_{opt}=2/\left(1+sin\left(\frac{\pi }{N+1}\right)\right)---\left(1.13\right)$

其中N为温度梯度方向的节点数。

第八步中，整个域等效热传导系数为通过边界的热流率q和温度梯度之比。对于 三维热传导问题，温度载荷在x方向，x方向的热导率计算公式为：

$K_{eff}=\frac{q}{▿T}=\frac{[\underset{k=1}{\overset{Nz}{\Sigma }}\underset{j=1}{\overset{Ny}{\Sigma }}{\mathrm{Kx}}_{(1,j,k)}(T_{2,j,k}-T_2)+\underset{k=1}{\overset{Nz}{\Sigma }}\underset{j=1}{\overset{Ny}{\Sigma }}{\mathrm{Kx}}_{(Nx-1,j,k)}(T_1-T_{Nx-1,j,k})]/(2Ny·Nz)}{(T_1-T_2)/Nx}---\left(1.14\right)$

Nx、Ny、Nz分别为x、y、z方向的胞元数。

温度载荷在y方向，y方向的热导率计算公式为：

$K_{eff}=\frac{q}{▿T}=\frac{[\underset{k=1}{\overset{Nz}{\Sigma }}\underset{i=1}{\overset{Nx}{\Sigma }}{\mathrm{Ky}}_{(i,1,k)}(T_{i,2,k}-T_2)+\underset{k=1}{\overset{Nz}{\Sigma }}\underset{i=1}{\overset{Nx}{\Sigma }}{\mathrm{Ky}}_{(i,Ny-1,k)}(T_1-T_{i,Ny-1,k})]/(2Nx·Nz)}{(T_1-T_2)/Ny}---\left(1.15\right)$

温度载荷在z方向，z方向的热导率计算公式为：

$K_{eff}=\frac{q}{▿T}=\frac{[\underset{i=1}{\overset{Nx}{\Sigma }}\underset{j=1}{\overset{Ny}{\Sigma }}{\mathrm{Kz}}_{(i,j,1)}(T_{i,j,2}-T_2)+\underset{i=1}{\overset{Nx}{\Sigma }}\underset{j=1}{\overset{Ny}{\Sigma }}{\mathrm{Kz}}_{(i,j,Nz-1)}(T_1-T_{i,j,Nz-1})]/\left(2Nx·Ny\right)}{(T_1-T_2)/Nz}---\left(1.16\right).$

实施例1：对等离子喷涂氧化钇稳定氧化锆(YSZ)涂层的沿喷涂方向热导率进行预测。

该实例以等离子喷涂YSZ涂层为例，展示了本发明对于多孔材料热导率的预测过程。

采用本发明对离子喷涂氧化钇稳定氧化锆(YSZ)涂层的热导率进行预测，包括以下 步骤：

(1)重构等离子喷涂YSZ涂层结构的三维图像，如图3所示；图3是本发明等离子喷涂 氧化钇稳定氧化锆(YSZ)涂层的三维重建图像(300×300×300像素)。图中黑色是 孔隙和裂纹，白色为YSZ。

(2)利用计算机语言程序，将图3所示图像读入内存中，以三维矩阵形式存储[A_(i,j,k)]；

(3)以像素构建胞元，按照其所属材质，赋予热导率[λ_(i,j,k)]；

(4)计算每个节点向在直角坐标系下三个方向的热传导系数，将值储存进热传导系数矩 阵[Kx_(i,j,k)]、[Ky_(i,j,k)]、[Kz_(i,j,k)]；

(5)离散稳态传热方程，应用于每个节点，构建计算方程；

(6)边界条件设置为：上表面所有胞元的温度设置为100℃，下表面胞元的温度设置为 0℃，四个侧表面设置为绝热条件；初始条件设置为：从上到下温度线性分布。

(7)以连续超松弛迭代算法求解整个域的温度场；

(8)用公式(1.15)计算涂层y方向的等效热导率。

等离子喷涂YSZ涂层是一种多孔的涂层，因此可以认为涂层中包含两种组分：YSZ实 体材料和孔隙，两种组分的热导率如表1所示。由于该涂层孔隙的结构十分复杂，包括了 各种空间形状的封闭孔裂纹和未结合界面等，因此采用二维模型很难模拟其内部的传热行 为，需采用三维模型进行计算。对于图3所示的等离子喷涂YSZ涂层的三维重建图像 (300×300×300像素)，其模型中含有27,000,000个胞元。表2中展示了采用有限差分法 和有限元法(ANSYS12.1)对进行图3进行建模计算所消耗的内存和时间比较。如果采用 有限元软件ANSYS计算，估算约需要超过70GB的内存，因此一般计算机根本无法完成计 算，而采用有限差分模型则可以解决这个问题。为了适应计算机的计算能力，将图3所示 图像分别分割为150×150×150像素和100×100×100像素。可以发现，对于150×150×150 像素有限差分模型仅需要ANSYS计算时间的1.2％和消耗内存的1.6％；而对于100×100×100 像素，有限差分模型仅需要ANSYS计算时间的1.5％和消耗内存的1.5％。这充分的说明了 本发明的三维模型计算效率高于有限元方法。

表3展示了采用有限差分法和有限元法(ANSYS12.1)对图3所示图像进行模拟计算 的热导率的比较(括号里为标准差)，同样地，150×150×150像素和100×100×100像素是 由原图像分割所得。可以发现，有限差分法计算得到的热导率值小于有限元法，对于 150×150×150像素的图像小17.9％，100×100×100像素的图像小16.9％。这是由于二者的插 值方程不同造成的。而参考该涂层实测的热导率0.99W·m^-1·K^-1，有限差分法计算结果更符 合实测值。

表4展示了采用有限差分法对图3所示图像分别进行二维和三维模拟计算的热导率结 果比较(括号里为标准差)，二维模拟的对象为该三维图像沿z轴不同层的截面。可以发 现，二维维模拟结果远低于三维结果，仅为后者的65.2％，也远低于实测值。

表1等离子喷涂YSZ涂层各组分热导率

表2：采用有限差分法和有限元法(ANSYS12.1)对进行图3进行建模计算所消耗的内存 和时间比较。

表3：采用有限差分法和有限元法(ANSYS12.1)对图3所示图像进行三维模拟计算的热 导率结果比较(括号里为标准差)，小图像由原图像分割所得。

表4：采用有限差分法对图3所示图像分别进行二维和三维模拟计算的热导率结果比较(括 号里为标准差)，二维模拟的对象为该三维图像沿z轴不同深度的截面。

通过上述实例可以证明，本发明能够通过三维图像预测材料的热导率，预测过程较为 简便迅速，消耗系统资源少，结果较二维模型更为准确。