Avizo用户使用手册-14

[TOC]

Chapter 14

14 Avizo XLabSuite Extension 用户指南

Avizo XLabSuite Extension提供数值仿真能力，用于从样品的3D图像（例如用CT、FIB/SEM、MRI等扫描得到）计算材料的物理属性。材料属性直接从已分割的3D图像计算得出。

所计算的材料属性有：

• 绝对渗透率(absolute permeability)，
• 分子扩散率(molecular diffusivity)，
• 形成因子与电导率(formation factor and electrical conductivity)，
• 热导率(thermal conductivity)。

注意：关于系统要求与硬件平台可用性，请参阅第1.4节 System Requirements。

对于其中每一种属性，都提供了两个不同的模块，对应两种不同的仿真思路。第一种思路是把某个给定属性估计为在实验室中所做实验的结果。为此，实验的外部条件由类似于实验室中存在的边界条件来模拟。这样就可以把模块结果与实际的实验结果作比较。第二种思路是有效属性计算。在这种情况下，该材料被视为它所取自的无限介质的代表。通过施加空间周期性边界条件，就可以获得宏观介质的有效属性。

实验仿真被设计为贴近真实的实验室实验。它考虑了非常强的扰动性边界条件：样品在四个面上被气密封闭，而在其余两个面上施加常数值。传输现象受到这类条件的高度约束。

张量计算基于一种数学方法，它把样品视为无限或宏观材料的代表。周期性是一种宽松得多的边界条件，传输现象更为自由。当达到代表性基元体(Representative Elementary Volume)时，张量计算通常是首选方法。

所提供的模块有：

• Absolute Permeability Experiment Simulation（绝对渗透率实验仿真）
  通过在四个面上气密封闭某个给定样品，同时在两个相对的面上添加实验装置以引导流动沿一个方向进行，从而仿真一个实验。
• Absolute Permeability Tensor Calculation（绝对渗透率张量计算）
  通过在一个代表性基元体上施加周期性边界条件，计算本征渗透率张量。
• Molecular Diffusivity Experiment Simulation（分子扩散率实验仿真）
  在实验室条件下求解菲克第二定律以仿真一个实验。
• Molecular Diffusivity Tensor Calculation（分子扩散率张量计算）
  用应用于菲克第二定律的体积平均法计算扩散率张量。所研究的样品被视为更大尺度材料的代表，从而允许施加周期性边界条件。
• Formation Factor Experiment Simulation（形成因子实验仿真）
  在实验室条件下求解欧姆定律以仿真一个实验。
• Effective Formation Factor Calculation（有效形成因子计算）
  用应用于欧姆定律的体积平均法计算电导率张量。所研究的样品被视为更大尺度材料的代表，从而允许施加周期性边界条件。形成因子由此推导得出。
• Thermal Conductivity Experiment Simulation（热导率实验仿真）
  在实验室条件下求解傅里叶定律以仿真一个实验。
• Thermal Conductivity Tensor Calculation（热导率张量计算）
  通过在一个代表性体上施加周期性边界条件，计算热导率张量。

以下各节针对每一种属性，概述计算模块所基于的理论，并给出这些模块的入门教程。

• 绝对渗透率计算入门
• 分子扩散率计算入门
• 形成因子与电导率计算入门
• 热导率计算入门

在这些教程中，有些步骤是必做的：你必须遵照执行才能成功完成教程。其他步骤是可选的：你至少应当阅读它们，因为它们的内容之后可能有用。某个步骤是可选的还是必做的，会在教程各节的开头指明。

致谢

Avizo XLabSuite Extension是与ICMCB-CNRS（法国佩萨克）研究主任Bernard博士合作开发的。

14.1 绝对渗透率计算入门

这个分步教程的目的，是让你更熟悉Avizo XLabSuite Extension为Avizo所提供的绝对渗透率计算模块的用法。本章将涉及以下主题：

• 关于绝对渗透率的理论基础
• 为仿真准备数据（步骤1到步骤6）：
  • 设定体素尺寸与单位
  • 分割孔隙空间
  • 用Avizo移除非渗流空间
  • 定义感兴趣区域
• 实验仿真（步骤7）：定义参数、运行仿真、解释并可视化结果
• 有效属性计算（步骤8）：定义参数、运行仿真、解释并可视化结果
• 用Kozeny-Carman方程验证结果（步骤9）。

还提供了一个演示脚本。该脚本会自动执行本教程中所详述的全部步骤。

14.1.1 理论要素

达西定律：绝对渗透率的定义

绝对渗透率被定义为对多孔材料传输单相流体能力的度量。它的SI单位是平方米($m^2$)，但平方微米($\mu m^2$)更为常用，因为它几乎等于一个达西(d)：$1d = 0.9869233\mu m^2$。它是材料的一种本征属性，与任何外部条件无关。

绝对渗透率在达西定律（见[1]）中作为一个联系流体、流动与材料参数的常数系数出现：

$$\frac{Q}{S} = -\frac{k}{\mu}\frac{\Delta P}{L}$$

其中：

• $Q$是通过该多孔介质的总流量（单位：$m^3 \cdot s^{-1}$）；
• $S$是流体所通过的样品横截面（单位：$m^2$）；
• $k$是绝对渗透率（单位：$m^2$）；
• $\mu$是流动流体的动力黏度（单位：$Pa \cdot s$）；
• $\Delta P$是施加在样品两端的压力差（单位：$Pa$）；
• $L$是样品在流动方向上的长度（单位：$m$）。

$\frac{Q}{S}$常记作$v$，代表通过该多孔介质的表观流速或平均流速，即达西速度。

绝对渗透率只考虑单相流体。多相流动涉及的是相对渗透率。

斯托克斯方程与流动条件

为了数值估计绝对渗透率，需要求解斯托克斯方程：

$$\begin{cases} \vec{\nabla} \cdot \vec{V} = 0 \\ \mu \nabla^2 \vec{V} - \vec{\nabla} P = \vec{0} \end{cases}$$

其中：

• $\vec{\nabla} \cdot$是散度算子；
• $\vec{\nabla}$是梯度算子；
• $\vec{V}$是材料流体相中流体的速度；
• $\mu$仍然是流动流体的动力黏度；
• $\nabla^2$是拉普拉斯算子；
• $P$是材料流体相中流体的压力。

这个方程组是纳维-斯托克斯方程的一种简化，它考虑了：

• 不可压缩流体，这意味着其密度是常数；
• 牛顿流体，这意味着其动力黏度是常数；
• 稳态流动，这意味着速度不随时间变化；
• 层流，这意味着所涉及的速度足够小，不会产生湍流。

最后一点等价于考虑低雷诺数下的流动（雷诺数的首次提出见[2]）。

一旦这个方程组被求解，估计渗透率系数就在于应用达西定律。该方程的所有值都可以从方程组的解($Q$、$\Delta P$)或者从外部条件($S$、$L$、$\mu$)推导出来。它由Absolute Permeability Experiment Simulation模块计算。

斯托克斯方程的体积平均形式

有效渗透率可以定义为固相对流体速度的影响。为了得到在整个体上有效的方程，需要进行一次尺度变换。体积平均法就是完成这种尺度变换的一种技术。它的主要目标是通过在一个体上对方程做平均，使方程在空间上变得平滑。感兴趣的读者可以在[3]中找到相关细节和参考文献。

这个非常一般的理论导出了一个闭合问题，它把斯托克斯方程变换为一个张量问题。尽管它是一个更高阶的问题，它仍与斯托克斯方程非常相似：

$$\begin{cases} \vec{\nabla} \cdot \overrightarrow{\overrightarrow{D}} = \vec{0} \\ \nabla^2 \overrightarrow{\overrightarrow{D}} - \vec{\nabla} \vec{d} = \overrightarrow{\overrightarrow{I}} \end{cases}$$

其中：

• $\overrightarrow{\overrightarrow{D}}$是一个张量，可视为速度空间偏差的源，我们称之为速度扰动场；
• $\vec{d}$是一个矢量，可视为压力空间偏差的源，我们称之为压力扰动场；
• $\overrightarrow{\overrightarrow{I}}$是单位张量。

渗透率张量通过在求解该系统所用的体$V$上计算$\overrightarrow{\overrightarrow{D}}$的平均值，从该问题的解中提取出来：

$$\overrightarrow{\overrightarrow{k}} = \frac{1}{V}\int_V \overrightarrow{\overrightarrow{D}}\,dV$$

这个渗透率张量给出了关于沿空间任意方向渗透率强度的额外信息。它可以体现多孔介质的各向异性，即渗透率强度对流动方向的依赖性。它由Absolute Permeability Tensor Calculation模块计算。

边界条件

Avizo XLabSuite Extension提供两种估计绝对渗透率的思路。

第一种是基于斯托克斯方程求解的实验仿真；这由Absolute Permeability Experiment Simulation模块完成。边界条件规定如下：

• 在流固界面处的无滑移条件。
• 在图像中不垂直于主流动方向的那些面上，添加一个一体素宽的固体平面（带无滑移条件）。这使样品与外界隔离，从而不允许流动从系统的输入面流出。
• 在图像中垂直于主流动方向的那些面上添加实验装置。它们的设计方式会创建出一个稳定区，在那里压力准静态，并且流体可以在样品的输入面上自由散布。
• 以下三个条件中的两个可由用户选择，第三个则从所选的两个估计得出：输入压力、输出压力、流量。

第二种思路求解由斯托克斯方程通过体积平均导出的闭合问题。这由Absolute Permeability Tensor Calculation模块完成。这种情况下所求解的张量问题，是通过对$\overrightarrow{\overrightarrow{D}}$、$\vec{d}$和几何形状施加周期性边界条件而闭合的。在流固界面处施加无滑移条件。该样品代表一种宏观的、无限的材料，因此必须是该多孔介质的代表。

人工压缩性

这些方程组无法用完全隐式方法（矩阵求逆）求解，因为这类系统的矩阵是奇异的。这就是为什么在系统中引入了一个人工压缩性系数和一些时间导数项。人工压缩性方法最早在[5]中被描述。

在方程组中引入这些项使得该问题可以被迭代求解。当时间导数趋于零时，就得到了唯一解。方程中所引入的时间没有物理意义。

方程组的离散化

Avizo XLabSuite Extension使用有限体积法来求解这些方程组。方程在一个交错网格（由[6]提出）上被离散化，从而能更好地估计无滑移边界条件。压力未知量位于体素中心，而速度未知量在体素的各个面上分解。

该离散化方案假定体素是各向同性的（立方体）。

参考文献

1. Darcy, H., Les fontaines publiques de la ville de Dijon, V. Dalmont, Paris, 1856
2. Reynolds, O., An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels, Philosophical Transactions of the Royal Society 174 (0): 935-982, 1883
3. Whitaker, S., The Method of Volume Averaging, Kluver Academic Publishers, 1999
4. Gray, W. G., A derivation of the equations for multiphase transport, Chemical Engineering Science, 30, 229-233, 1975
5. Chorin, A. J., A Numerical Method for Solving Incompressible Viscous Flow Problems, Journal of Computational Physics, 2, 12-26, 1967
6. Harlow, F. H., and Welch, J. E., Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Physics of Fluids, v.5, p.317, 1965

14.1.2 步骤1 - 激活Avizo单位管理

目的	激活Avizo中所提供的单位管理。
必做步骤	是。

长度单位是为了得到准确的渗透率结果而必须正确缩放的最重要参数之一。例如它被存储在Avizo格式中，但并非所有数据格式都保存这一信息。Avizo提供了一个长度单位管理工具，本教程必须激活它。

• 开始一个新项目 File > New Project
• 打开Edit > Preferences...对话框。
• 选择Units标签页。
• 选择Spatial information only。保持所有复选框都被勾选（见图14.1）。
• 单击OK。

图 14.1：(1) 在Edit > Preferences...对话框中选择Units标签页。(2) 选择Spatial information only，并 (3) 保持所有复选框都被勾选。

一旦被激活，每当加载一个数据集时，单位管理工具就会询问体素尺寸所要使用的长度单位。关于该功能的更多细节，请参阅Units in Avizo。

14.1.3 步骤2 - 加载数据集并为体素尺寸选择长度单位

目的	在Avizo中加载数据集并设定单位长度。
必做步骤	是。

本教程所使用的数据集是一堆随机堆积的玻璃球。它由ICMCB-CNRS研究主任Dominique Bernard扫描。这些球形颗粒经过筛分，直径在100–120微米范围内。该样品在扫描之前于700摄氏度下烧结了10分钟。图14.2展示了该数据集。

本教程所使用的完整3D数据集是一个200×200×200的立方体。

• 打开File > Open Data...对话框。
• 从AVIZO_ROOT目录打开data/tutorials/xlab/10mc3_200.vol.am文件（见图14.3）。
• 单击Open。

注意：文件10mc3_200.vol.am嵌入了单位信息，因此该数据的长度单位会被自动设定。不过，文件10mc3_400.vol.am没有这样的单位信息。由于Avizo单位管理工具在步骤1中已被激活，你需要在每次打开该文件时指定长度单位。每当加载一个没有长度单位信息的数据集时，图14.4所展示的对话框就会出现。例如，对于文件10mc3_400.vol.am，请选择micrometer [μm]作为本教程数据集的坐标单位。

注意：如果在选择长度单位之后保存了数据集，该信息会以Avizo格式存储，并在下次加载时被重用。之后可以使用Units Editor（见图14.5）以及一个看起来像图14.4那样的对话框来修改长度单位。

图 14.2：本教程所用球体堆积的3D可视化。

图 14.3：用于加载本教程数据集的Open Data...对话框。选择数据集的路径，然后选中并打开它。

图 14.4：Units Editor对话框。使用下拉菜单为体素尺寸选择长度单位。每当加载一个不含长度单位信息的数据集时，该对话框都会出现。

图 14.5：Units Editor按钮位于红框中（仅当在Project View中选中某个数据集时可见）。

14.1.4 步骤3 - 设定体素尺寸

目的	设定均匀数据集中体素的边长。
必做步骤	是。

我们先来显示所加载的数据集：

• 右键单击10mc3_200.vol.am并选择Bounding Box。
• 使用Ortho Slice模块在该体中浏览。

图 14.6：在步骤4开始时，连接好可视化模块之后，查看器和Project View应当看起来像这张图。

例如，体素尺寸被存储在Avizo格式中，但并非所有数据格式都保存这一信息。万一所存储的值有误，Avizo提供了一个编辑器，允许修改所加载数据集的体素尺寸。

• 通过单击图14.7中高亮显示的按钮打开Crop Editor。
• 可以在出现的对话框中修改数据集的体素尺寸（见图14.8）。

图 14.7：当数据集在Project View中被选中时，找到Crop Editor按钮。

图 14.8：在Crop Editor中修改体素尺寸所需填写的字段。

对于本教程的数据集，正确的体素尺寸是3.8，它已被默认设定。

如果Avizo单位管理未被激活，数据集中所指示的体素尺寸将被视为以微米为单位。否则，长度单位是在数据集加载时被设定的。

注意：Avizo XLabSuite Extension的计算必须使用各向同性（即立方体）体素。

14.1.5 步骤4 - 从数据集创建标签图像

目的	使用基本分割工具从数据集创建一个标签图像。
必做步骤	是。

本教程中的分割很简单，因为所使用的数据集从一开始就几乎是二值的。由于这不是本教程的重点，关于Avizo中分割工具的更多细节，请参阅用户指南的3D图像的分割一章。

应用一个快速平滑滤波器以移除图像中的噪声。

• 在Project View中右键单击10mc3_200.vol.am，并选择Image Processing > Smoothing and Denoising > Median Filter。
• 把Filter的interpretation设为3D（见图14.9）。
• 单击Apply。
• 在Project View中右键单击滤波后图像的结果(10mc3_200.vol.filtered)，并选择Display > Ortho Slice（见图14.10）。

图 14.9：Filter > Median Filter模块被修改的参数以红色箭头高亮显示。

图 14.10：滤波后数据集的可视化。查看器和Project View应当看起来像这张图。

图 14.11：Segment > Multi-Thresholding模块被修改的参数以红色箭头高亮显示。

图 14.12：阈值化后数据集的可视化。查看器和Project View应当看起来像这张图。

通过对平滑后的图像做阈值处理来创建一个标签图像。

• 在Project View中右键单击10mc3_200.vol.filtered，并选择Image Segmentation > Multi-Thresholding。
• 默认选项应当与图14.11相同。最重要的、需要核实的参数是Exterior-Range1滑块的值，它必须被设为114。
• 单击Apply。
• 在Project View中右键单击阈值操作的结果(10mc3_200.vol.labels)，并选择Display > Ortho Slice（见图14.12）。

14.1.6 步骤5 - 移除非渗流的孔隙空间

目的	使用高级分割工具移除阈值化图像中非渗流的孔隙空间。
必做步骤	否。仅在Avizo中可用。

在渗透率计算之前可以执行的一项初步测试是渗流测试(percolation test)。其目的是确保材料中存在的孔隙空间允许流体从一侧流到另一侧。如果这项测试不通过，就意味着从样品某个面注入的流体找不到任何路径从对面流出。这类材料的渗透率为零，但数值收敛到0可能相当漫长。从离散的角度看，所需的连通性是6连通：流体只能通过体素的面从一个体素流到它的邻居。

首先，球体的值为1，孔隙的值为0。为了进行渗流处理，这些值必须被取反。

• 在Project View中右键单击10mc3_200.vol.labels，并选择Compute > Logical Operations > Invert。
• 单击Apply。
• 在Project View中右键单击结果(10mc3_200.vol.invert)，并选择Display > Ortho Slice（见图14.13）。

图 14.13：取反后数据集的可视化。查看器和Project View应当看起来像这张图。

渗流测试可以应用到这个取反后的图像上。

• 在Project View中右键单击10mc3_200.vol.invert，并选择Image Processing > Propagation > Axis Connectivity。
• 在Neighborhood端口中选择6。
• 单击Apply。
• 在Project View中右键单击结果(10mc3_200.vol.Axis-Connectivity)，并选择Display > Ortho Slice（见图14.14）。

图 14.14：数据集中渗流孔隙率的可视化。查看器和Project View应当看起来像这张图。

这幅图像包含孔隙空间中所有这样的值为1的体素：从平面$z = 0$出发，只在值为1的体素中移动，一直到平面$z = 200$，都能到达它们。孤立的孔隙空间被移除，这使得后续计算稍快一些。当所得到的图像为空时，这项测试的用处就更大了。那意味着不存在连通的孔隙率，渗透率为零。

14.1.7 步骤6 - 子区域的选择

目的	为渗透率计算定义一个感兴趣区域。
必做步骤	是。

有一种简便方法可以在所加载数据集的某个子区域上计算渗透率，而不必裁剪它。可以定义一个感兴趣区域(ROI)并把它连接到Avizo XLabSuite Extension的模块上，这样计算就只在该ROI中进行。

• 在Project View中右键单击10mc3_200.vol.Axis-Connectivity，并选择Display > ROI Box。
• 把这个新模块Minimum端口的三个可编辑字段设为285，Maximum端口的三个可编辑字段设为473.1。见图14.15。

图 14.15：Display > ROI Box模块被修改的参数以红色矩形高亮显示。

该ROI定义了一个以完整体为中心、边长为50个体素的立方体（见图14.16）。本教程将用它在小体积上计算渗透率。在教程过程中，它可以随意在数据集中移动，不过为了便于说明，将使用它的初始位置。

图 14.16：代表感兴趣区域的立方体的可视化。查看器和Project View应当看起来像这张图。

14.1.8 步骤7 - 绝对渗透率实验仿真

目的	仿真一个绝对渗透率实验。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Axis-Connectivity（如果步骤5未完成，则为10mc3_200.vol.labels），并选择XLab Simulations > Absolute Permeability Experiment Simulation。
• 把该ROI Box连接到Absolute Permeability Experiment Simulation模块的ROI输入连接上。
• 如果步骤5已完成，请在Pore space端口中勾选带有1的复选框，因为10mc3_200.vol.Axis-Connectivity中孔隙空间的标签是1。如果步骤5未完成，请在Pore space中勾选带有0的复选框，因为10mc3_200.vol.labels中孔隙空间的标签是0。
• 如果步骤5已完成，模块参数应当看起来像图14.17。
• 单击Apply。

图 14.17：XLab Simulations > Absolute Permeability Experiment Simulation模块被修改的参数以红色箭头高亮显示。(1) 连接ROI Box模块以缩小计算域。(2) 选择标签1，因为它代表球体之间的孔隙空间。

默认参数仿真的是沿Z轴、输入压力为$1.3 \times 10^5$ Pa、输出为大气压的一个实验。默认黏度是水的黏度。这些选项可以被修改以仿真若干不同的实验。例如，主流动的方向可以被调整为X、Y或Z方向（默认为Z）。如果选择了若干个方向，计算会被依次执行。实验的边界条件也可以被修改，从而使速度场和压力场按这些值缩放。三者之中可以施加两个值：输入压力、输出压力、流量。修改这些值不会改变渗透率——它是该多孔介质的本征属性。它只会修改输出场。

14.1.8.1 检索并解释结果

Project View中出现了四个输出：

• 10mc3_200.vol.KExp.Spreadsheet：一个电子表格，包含计算最相关的结果：
  • 描述材料几何形状的数据集名称；
  • 计算所在的感兴趣区域；
  • 以$\mu m^2$为单位的渗透率值；
  • 以$d$（达西）为单位的渗透率值。
• 10mc3_200.vol.KExp.Error.Spreadsheet：一个电子表格，包含每次迭代时收敛判据的估计值。收敛判据相对于迭代次数的曲线由Plot 10mc3_200.vol.KExp.Error.Spreadsheet显示。
• 10mc3_200.vol.VelocityZ：一个矢量场，表示速度场——它是在由ROI所定义几何形状中求解的斯托克斯方程组解的一部分。它的单位是$\mu m \cdot s^{-1}$。
• 10mc3_200.vol.PressureZ：一个标量场，表示压力场——它是在由ROI所定义几何形状中求解的斯托克斯方程组解的第二部分。它的单位是$Pa$。

包含全部结果及相关信息的电子表格，可以通过在project view中选中它并单击Show按钮（在Properties区域中）来可视化。见图14.18。

图 14.18：包含Absolute Permeability Experiment Simulation计算主要结果的电子表格。

电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

14.1.8.2 可视化输出场

要可视化速度场：

• 通过单击它们的查看器开关，隐藏Bounding Box和那五个Ortho Slice。
• 右键单击10mc3_200.vol.VelocityZ并选择Compute > Magnitude。
• 右键单击10mc3_200.vol.VelocityZ并选择Display > Illuminated Streamlines，并按图14.19所述设定它的参数。
• 单击Apply。

所得到的可视化和Project View如图14.20所述。

图 14.19：表示实验仿真中速度场的Illuminated Streamlines被修改的参数以红色箭头高亮显示。

图 14.20：表示实验仿真中速度场的Illuminated Streamlines的可视化。查看器和Project View应当看起来像这张图。

要可视化压力场：

• 通过单击它的查看器开关隐藏Illuminated Streamlines模块。
• 右键单击10mc3_200.vol.PressureZ并选择Display > Height Map Slice。
• 按图14.21配置该Height Map Slice。

图 14.21：表示实验仿真中压力场的Height Map Slice被修改的参数以红色箭头高亮显示。

图 14.22：表示实验仿真中压力场的Height Map Slice的可视化。查看器和Project View应当看起来像这张图。

注意：实验装置是在计算启动之后自动设计的。图14.23描述了它们的形状。在垂直于流动方向的那些面上，会添加一个简单的形状。它由以下部分构成：

• 一个方形截面的通道，用于（从数值角度）稳定系统中的流动，以便于迭代算法的收敛；
• 一个渐扩部分，用于尽量减少添加到系统中的孔隙空间量，并把流动在输入面上散开；
• 一个完全为流体的区域，以确保流体能够通过样品的完整表面进入。

在样品的其他面上，会添加一个一体素宽的固体平面，以确保流体被约束在系统内部。

图 14.23：所创建实验装置的2D示例。该实验装置（系统的左侧和右侧）使流动能够在样品中散开。

14.1.9 步骤8 - 绝对渗透率张量计算

目的	计算本征渗透率张量。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Axis-Connectivity（如果步骤6未完成，则为10mc3_200.vol.Labels），并选择XLab Simulations > Absolute Permeability Tensor Calculation。
• 把该ROI Box连接到Absolute Permeability Tensor Calculation模块的ROI输入连接上。
• 如果步骤5已完成，请在Pore space端口中勾选带有1的复选框，因为10mc3_200.vol.Axis-Connectivity中孔隙空间的标签是1。如果步骤5未完成，请在Pore space中勾选带有0的复选框，因为10mc3_200.vol.Labels中孔隙空间的标签是0。
• 如果步骤5已完成，模块参数应当看起来像图14.24。
• 单击Apply。

使用这些参数，该模块将计算完整的本征渗透率张量。一次完整的张量计算需要三次在时间与内存消耗上等同于实验仿真的计算。

图 14.24：XLab Simulations > Absolute Permeability Tensor Calculation模块被修改的参数以红色箭头高亮显示。(1) 连接ROI Box模块以缩小计算域。(2) 选择标签1，因为它代表球体之间的孔隙空间。

14.1.9.1 检索并解释结果

Project View中出现了五个输出：

• 10mc3_200.vol.KTensor.Spreadsheet：一个电子表格，包含计算最相关的结果：
  • 描述材料几何形状的数据集名称；
  • 计算所在的感兴趣区域；
  • 以3×3矩阵形式给出的完整渗透率张量，单位为$\mu m^2$；
  • 该张量的特征系统。特征值及其所关联的向量被描述在一行中。
• 10mc3_200.vol.KTensor.Error.Spreadsheet：一个电子表格，包含每次迭代、每个方向上收敛判据的估计值。收敛判据相对于迭代次数的曲线由Plot 10mc3_200.vol.KTensor.Error.Spreadsheet显示。
• 10mc3_200.vol.VelocityX、10mc3_200.vol.VelocityY、10mc3_2.vol.VelocityZ：三个矢量场，表示用体积平均方法从斯托克斯方程导出的张量问题的$\overrightarrow{\overrightarrow{D}}$解。

表示同一问题$\vec{d}$解的标量场结果没有在本教程中展示，因为它们的可视化通常不提供有价值的信息。

包含该张量及所有相关信息的电子表格，可以通过在Project View中选中它并单击Show按钮（在Properties区域中）来可视化。见图14.25。

图 14.25：包含Absolute Permeability Tensor Calculation计算主要结果的电子表格。

电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

14.1.9.2 可视化输出场

要可视化速度场：

• 通过单击它的查看器开关隐藏Height Map Slice。
• 右键单击10mc3_200.vol.VelocityX并选择Compute > Magnitude。
• 右键单击10mc3_200.vol.VelocityX并选择Display > Illuminated Streamlines，并按图14.26所述设定它的参数。
• 单击Apply。

所得到的可视化和Project View如图14.27所述。

图 14.26：用于绝对渗透率张量计算速度场结果的DisplayISL模块被修改的参数。

图 14.27：表示本征渗透率张量某一行计算中速度场的Illuminated Streamlines的可视化。查看器和Project View应当看起来像这张图。

对每一个速度场重复同样的流程，就可以把其他速度场也可视化出来。

14.1.10 步骤9 - 用Kozeny-Carman方程验证渗透率

目的	使用Kozeny-Carman方程检查结果的准确性。
必做步骤	否。

本教程所使用的数据集是随机堆积的球体。对于这类材料，可以使用Kozeny-Carman方程从孔隙率和球直径估计渗透率。Kozeny-Carman方程可以写成：

$$k_{KC} = \frac{1}{180}\frac{\epsilon^3}{(1-\epsilon)^2}d^2$$

其中：

• $\epsilon$是球体堆积的孔隙率；
• $d$是球的直径。

我们这里使用的是数据集的完整体、高分辨率版本。它可以在与本教程所用重采样数据集相同的目录中找到：data/tutorials/xlab/10mc3_400.vol.am。

球的平均直径可以使用Avizo工具集来确定，其流程类似于教程示例4：进一步的图像分析——泡沫中孔隙直径的分布。使用实验参数，球的直径在100与120$\mu m$之间。

孔隙率也可以用Avizo估计，其值为36.35%。

由这两个要素，从Kozeny-Carman方程计算出的渗透率值在$6.59d$与$8.00d$之间。

这些值可以与在完整几何体上（不使用本教程中所用的感兴趣区域）用绝对渗透率仿真模块所得到的渗透率作比较。每个方向上的实验仿真给出以下值：

$$\begin{aligned} k_X &= 7.80d \\ k_Y &= 7.85d \\ k_Z &= 7.75d \end{aligned}$$

该材料的本征渗透率张量为：

$$\begin{bmatrix} 7.57 & 0.25 & -0.23 \\ 0.25 & 7.54 & -0.18 \\ -0.23 & -0.18 & 7.15 \end{bmatrix}$$

仿真的实验值以及该张量的对角值，全都落在用Kozeny-Carman方程所估计的渗透率范围内。我们提醒用户：这些值是用数据集的完整体、高分辨率版本计算的：data/tutorials/xlab/10mc3_400.vol.am。还提供了第二个演示脚本，它在完整分辨率数据集上、限制在一个感兴趣区域内运行本教程的各个步骤。

14.2 分子扩散率计算入门

这个分步教程的目的，是帮助你更熟悉Avizo XLabSuite Extension为Avizo所提供的分子扩散计算模块的用法。本章将涉及以下主题：

• 关于分子扩散的理论基础
• 为仿真准备数据（步骤1到步骤5）
• 实验仿真（步骤6）：定义参数、运行仿真、解释并可视化结果
• 有效属性计算（步骤7）：定义参数、运行仿真、解释并可视化结果
• 用若干经验定律验证结果（步骤8）。

由于为仿真准备数据这一主题已在前一个关于绝对渗透率计算的教程中讲过，这里只作简要提及：

• 设定体素尺寸与单位，
• 分割孔隙空间，
• 用Avizo移除非渗流空间，
• 定义感兴趣区域。

还提供了一个演示脚本（仅适用于Microsoft Windows）。该脚本会自动执行本教程中所详述的全部步骤。

14.2.1 理论要素

菲克第一定律：分子扩散的定义

分子扩散是这样一个过程：溶解的质量通过随机的分子运动，从较高化学能状态被动地传输到较低化学能状态。自由溶液中某种化学物质的稳态扩散，可以用菲克第一定律来经验性地描述：

$$\vec{j} = -D \cdot \vec{\nabla}c$$

其中：

• $\vec{j}$是溶质质量通量（单位：$mol \cdot m^{-2} \cdot s^{-1}$），
• $D$是溶质在溶剂中的扩散系数（单位：$m^2 \cdot s^{-1}$），
• $c$是溶质在溶剂中的浓度（单位：$mol \cdot m^{-3}$）。

菲克第二定律

描述在均质（只有一个固相）、饱和（材料的孔隙空间充满溶剂）多孔介质中瞬态扩散的偏微分方程，可以从菲克第一定律与质量守恒推导出来。这个方程称为菲克第二定律，写作：

$$\frac{\partial c}{\partial t} - D \cdot \nabla^2 c = 0$$

在Molecular Diffusivity Experiment Simulation模块中，建议采用一个经典实验，它基于双储液室测试。两个具有相同体积$VR$的储液室被放置在样品沿所选方向的两侧。其他方向用气密平面封闭，因此没有扩散发生。两个储液室的初始浓度是不同的：$C{in}(t)$和$C{out}(t)$。样品初始时充满浓度为$C{in}(t_0)$的溶液，$t_0$是实验开始的瞬间。在时刻$t = t_0$，储液室与样品连通，扩散过程开始。重力的影响被忽略，只考虑被动扩散，不考虑对流。

考虑这些边界条件，菲克第二定律支配着扩散，并定义了样品中的浓度场。储液室的浓度也会演变，因为它们具有有限的体积$VR$。默认情况下，$V_R$被假定为样品中孔隙空间体积的100倍。我们记$\beta = \frac{V{voidspace}}{V_R}$为孔隙空间体积与储液室体积之比。

以下方程支配储液室中的浓度：

$$V_R\frac{\partial C_{in}(t)}{\partial t} = D\int_{S_{in}} \vec{\nabla}c \cdot \vec{n}\,dS$$ $$V_R\frac{\partial C_{out}(t)}{\partial t} = -D\int_{S_{out}} \vec{\nabla}c \cdot \vec{n}\,dS$$

其中$S{in}$和$S{out}$是样品上与储液室相连的那些面。

扩散过程一旦开始，样品中的浓度会快速演变，与储液室的交换是不对称的。随后，当与储液室的交换相等时，这个瞬态状态就被一个已建立的状态所取代。

这个已建立的状态由以下事实所刻画：$\frac{\partial C{in}(t)}{\partial t} = -\frac{\partial C{out}(t)}{\partial t}$。一旦达到该状态，储液室中的浓度会继续变化，直到它们达到平衡浓度$c\infty$。这些浓度之差$C{out}(t) - C_{in}(t)$遵循一个指数定律：

$$C_{out}(t) - C_{in}(t) = p \cdot exp(-\lambda^2 t)$$

其中$p$和$\lambda^2$是待确定的常数系数。

针对这个问题给出了一个解析解：

$$c(X,t) = A\left[\left[\frac{cos\left(\frac{\lambda}{\sqrt{D_{app}}}\right)-1}{sin\left(\frac{\lambda}{\sqrt{D_{app}}}\right)}\right]cos\left(\frac{\lambda}{\sqrt{D_{app}}}X\right) + sin\left(\frac{\lambda}{\sqrt{D_{app}}}X\right)\right]exp(-\lambda^2 t) + c_\infty$$

其中$c(X,t)$是位置$X$、时刻$t$处的局部浓度，$A$是一个常数系数。

已知该解必须满足前述通量相等的假设($\frac{\partial C{in}(t)}{\partial t} = -\frac{\partial C{out}(t)}{\partial t}$)，可推导出以下方程：

$$-\lambda^2\frac{cos\left(\frac{\lambda}{\sqrt{D_{app}}}\right)-1}{sin\left(\frac{\lambda}{\sqrt{D_{app}}}\right)} = D_{app}\beta\frac{\lambda}{\sqrt{D_{app}}}$$

它把$\lambda^2$系数与样品的表观扩散率$D_{app}$联系起来。

总结如下：

1. 在已建立的状态出现之前，必须先经历一个扩散过程启动的初始瞬态状态。
2. 一旦已建立的状态开始，储液室浓度之差遵循一个指数定律。因此，$ln(C{out}(t) - C{in}(t))$所遵循的线性曲线的斜率可以被轻松估计。
3. 这个斜率就是$\lambda^2$，即指数系数，它与表观扩散率$D_{app}$相关。

菲克定律的体积平均形式

有效分子扩散率张量给出关于材料扩散能力的全局信息，它由Molecular Diffusivity Tensor Calculation模块计算。

为了得到在整个体上有效的方程，需要进行一次尺度变换。体积平均法就是完成这种尺度变换的一种技术。它的主要目标是通过在一个体上对方程做平均，使方程在空间上变得平滑。感兴趣的读者可以在[Whitaker, S., The Method of Volume Averaging, Kluver Academic Publishers, 1999]中找到相关细节和参考文献。

该理论导出一个闭合问题，它把菲克方程变换为一个矢量问题；闭合变量$\vec{b}$被用于在一个新问题中表述浓度扰动：

$$\nabla^2\vec{b} = 0$$

当该问题被求解后，就可以计算如下定义的无量纲扩散率张量：

$$\epsilon\frac{\overrightarrow{\overrightarrow{D}}}{D_{solution}} = \epsilon\left(\overrightarrow{\overrightarrow{I}} + \frac{1}{V_f}\int_{S_{fs}} \overrightarrow{n_{fs}}\,\vec{b}\,ds\right)$$

其中：

• $\epsilon$是孔隙率，
• $D_{solution}$是本体溶液的扩散率，
• $V_f$是流体的体积，
• $S_{fs}$是流固界面的面积，
• $\overrightarrow{n_{fs}}$是流固界面的法线，方向从流体指向固相。

边界条件

Avizo XLabSuite Extension提供两种估计分子扩散率的思路。

第一种是基于菲克方程求解的实验仿真。如前所述，这由Molecular Diffusivity Experiment Simulation模块完成。固体的反应速率被假定为零：流固界面处不发生反应。于是流固界面处的边界条件为：

$$-\overrightarrow{n_{fs}} \cdot \vec{\nabla}c = 0$$

其中$\overrightarrow{n_{fs}}$是流固界面的法线，方向从流体指向固相。

除了这个流固界面条件之外，在图像中不垂直于主扩散方向的那些面上，还添加了一个一体素宽的固体平面。这使样品与外界隔离。

入口和出口处的边界条件需要知道储液室中的浓度。这些浓度随时间演变：

$$V_r\frac{\partial C_{in}}{\partial t} = \int_{S_{in}} \overrightarrow{n_{S_{in}}} \cdot \vec{\nabla}c\,ds$$ $$V_r\frac{\partial C_{out}}{\partial t} = -\int_{S_{out}} \overrightarrow{n_{S_{out}}} \cdot \vec{\nabla}c\,ds$$

其中$S{in}$和$S{out}$分别是样品的输入面和输出面；$\overrightarrow{n{S{in}}}$和$\overrightarrow{n{S{out}}}$分别是输入面和输出面的法线。

第二种思路求解由菲克方程通过体积平均导出的闭合问题。这由Molecular Diffusivity Tensor Calculation模块完成。这种情况下所求解的矢量问题，是通过对$\vec{b}$和几何形状施加周期性边界条件而闭合的。流固界面条件具有以下类似形式：

$$-\overrightarrow{n_{fs}} \cdot \vec{\nabla}\vec{b} = \overrightarrow{n_{fs}}$$

方程组的离散化

Avizo XLabSuite Extension使用有限体积法来求解这些方程组。

该离散化方案假定体素是各向同性的（立方体）。

实验仿真的系统求解

实验仿真在入口和出口处的边界条件需要知道储液室中的浓度。这些浓度随时间演变，这使得必须对该问题进行显式求解，而不能像有效扩散率或表观/有效电导率那样采用直接求解。

周期性边界条件下的系统求解

一旦离散化，闭合方程组可以写成$Ax = b$，其中$A$是一个稀疏对称矩阵。

该方程组使用完全隐式方法（矩阵求逆）求解。线性系统的直接求解使用PETSc（Portable, Extensible Toolkit for Scientific Computation）库。

同时也执行一次带共轭梯度与ILU预条件子的迭代求解。所使用的收敛判据是残差$l_2$范数的相对下降。

14.2.2 步骤1 - 加载数据集

目的	在Avizo中加载数据集。
必做步骤	是。

本教程所使用的数据集是一堆随机堆积的玻璃球。它由ICMCB-CNRS研究主任Dominique Bernard扫描。这些球形颗粒经过筛分，直径在100–120微米范围内。该样品在扫描之前于700摄氏度下烧结了10分钟。图14.28展示了该数据集。

本教程所使用的完整3D数据集是一个200×200×200的立方体。

• 开始一个新项目 File > New Project
• 打开File > Open Data...对话框。
• 从AVIZO_ROOT目录打开data/tutorials/xlab/10mc3_200.vol.am文件。
• 单击Open。

图 14.28：本教程所用球体堆积的3D可视化。

如前一步所述，单位管理在这里不是必做的。不过，请参阅绝对渗透率仿真教程的步骤2以了解更多关于编辑长度单位的内容。

14.2.3 步骤2 - 设定体素尺寸

目的	设定均匀数据集中体素的边长。
必做步骤	否。

同样，如前面各步骤所述，正确设定体素的尺寸在这里不是必做的。不过，要了解更多关于使用Crop Editor设定体素尺寸的内容，请参阅绝对渗透率仿真教程的步骤3。

注意：Avizo XLabSuite Extension的计算必须使用各向同性（即立方体）体素。

14.2.4 步骤3 - 从数据集创建标签图像

目的	使用基本分割工具从数据集创建一个标签图像。
必做步骤	是。

这一步是使用Avizo所提供的一些分割工具，从被平滑并阈值化的初始数据集创建一个标签图像。

请遵照绝对渗透率仿真教程步骤4中所述的说明操作。

图 14.29：阈值化后数据集的可视化。

14.2.5 步骤4 - 移除非渗流的孔隙空间

目的	使用高级分割工具移除阈值化图像中非渗流的孔隙空间。
必做步骤	否。仅在Avizo中可用。

渗流测试可以在运行计算之前执行。这项测试的一个好处是，它会把所有孤立的孔隙空间从标签图像中移除（我们提醒用户：扩散系数以及固相的反应速率都被假定为零）。这样计算应当会稍快一些。从离散的角度看，所需的连通性是6连通：浓度只能在非零表面上被估计，也就是说在两个相邻体素之间的一个面上。

如果需要，你可以遵照绝对渗透率仿真教程步骤5中所述的说明来执行渗流测试。

图 14.30：数据集中渗流孔隙率的可视化。

14.2.6 步骤5 - 子区域的选择

目的	为分子扩散率计算定义一个感兴趣区域。
必做步骤	是。

有一种简便方法可以在所加载数据集的某个子区域上计算分子扩散率，而不必裁剪它。可以定义一个感兴趣区域(ROI)并把它连接到Avizo XLabSuite Extension的模块上，这样计算就只在该ROI中进行。

请遵照绝对渗透率仿真教程步骤6中给出的说明，定义一个以样品为中心、边长为50个体素的ROI立方体（见图14.31）。本教程将用该ROI在小体积上计算属性。在教程过程中，它可以随意在数据集中移动，不过为了便于说明，将使用它的初始位置。

为方便起见，同时也因为这会在本教程阶段节省时间，我们在样品的一个子区域上计算分子扩散率。不过，我们将在验证阶段看到重采样和使用ROI如何影响所计算结果的质量。

图 14.31：代表感兴趣区域的立方体的可视化。

14.2.7 步骤6 - 分子扩散率实验仿真

目的	仿真一次分子扩散率实验室测量。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Axis-Connectivity（如果步骤4未完成，则为10mc3_200.vol.Labels），并选择XLab Simulations > Molecular Diffusivity Experiment Simulation。
• 把该ROI Box连接到Molecular Diffusivity Experiment Simulation模块的ROI输入连接上。
• 如果步骤4已完成，请在Pore space端口中勾选指向标签1的复选框，因为10mc3_200.vol.Axis-Connectivity中孔隙空间的标签是1。如果步骤4未完成，请在Pore space中勾选指向标签0的复选框，因为10mc3_200.vol.Labels中孔隙空间的标签是0。
• 单击Apply。

默认参数仿真的是沿Z轴的一个实验：输入储液室在初始时刻的浓度为1711 $mol \cdot m^{-3}$，输出储液室在初始时刻的浓度为零。默认的溶液本体扩散率是1 $m^2 \cdot s^{-1}$。这些选项可以被修改以仿真若干不同的实验。例如，分子扩散的方向可以被调整为X、Y或Z方向（默认为Z）。如果选择了若干个方向，计算会被依次执行。用作实验边界条件的浓度值也可以被修改。修改这些值不会改变分子扩散率——它是该多孔介质的本征属性。它只会修改输出的浓度场。

14.2.7.1 检索并解释结果

Project View中出现了四个输出：

• 10mc3_200.DExp.Spreadsheet：一个电子表格，包含计算最相关的结果：
  • 数据集名称，
  • 计算所在的感兴趣区域，
  • 以$m^2 \cdot s^{-1}$为单位的表观分子扩散率，
  • 以$mol \cdot m^{-3}$为单位、在初始时刻施加于输入储液室的浓度，
  • 以$mol \cdot m^{-3}$为单位、在初始时刻施加于输出储液室的浓度，
  • 以$m^2 \cdot s^{-1}$为单位的溶液本体扩散率。
• 10mc3_200.vol.DExp.ResConcentration.Spreadsheet：一个电子表格，包含每次迭代时输入与输出储液室的均匀浓度。
• 10mc3_200.vol.DExp.Error.Spreadsheet：一个电子表格，包含每次迭代时收敛判据的估计值。收敛判据相对于迭代次数的曲线由Plot 10mc3_200.vol.DExp.Error.Spreadsheet显示。
• 10mc3_200.vol.ConcentrationZ：一个标量场，表示在受ROI限制的几何形状中所求解菲克方程组的解——浓度场（单位：$mol \cdot m^{-3}$）。

包含全部相关信息的电子表格10mc3_200.vol.DExp.Spreadsheet，可以通过在Project View中选中它并单击Show按钮（在Properties area中）来可视化。见图14.32。电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

图 14.32：包含Molecular Diffusivity Experiment Simulation计算主要结果的电子表格。

14.2.7.2 可视化储液室浓度的演变

要可视化储液室浓度的演变：

• 选择10mc3_200.vol.DExp.ResConcentration.Spreadsheet，并在宏按钮中选择Plot Spreadsheet。
• 在Plot Spreadsheet属性区域中，在Y端口中同时选择DExp.Z.Input concentration和DExp.Z.Output concentration。
• 按下Show按钮。

所得到的可视化应当看起来像图14.33。可以注意到，输入与输出浓度还远未达到终态浓度（在理论页中记作$C_\infty$）。正如理论页中所解释的，一旦输入浓度与终态浓度之差的对数相对于时间的行为变为线性，计算就会停止。

图 14.33：输入与输出储液室浓度演变的可视化。

14.2.7.3 可视化输出浓度场

要可视化浓度场：

• 通过单击它们的查看器开关，隐藏当前所有的显示模块，例如Bounding Box或Ortho Slice。
• 右键单击10mc3_200.vol.ConcentrationZ并选择Display > Ortho Slice。
• 在Ortho Slice属性中，从Edit菜单把所选颜色图改为temperature.icol。
• 再次选择颜色图的Edit菜单，然后选择Adjust range to > Data min-max。
• 在Orientation端口中选择yz。
• 放大该切片以可视化它。

所得到的可视化应当看起来像图14.34。你可以观察到浓度从输入储液室（在切片底部，黄色）到输出储液室（在切片顶部，浅蓝色）的下降。

图 14.34：Z方向分子扩散实验仿真中浓度场的可视化。

14.2.8 步骤7 - 分子扩散率张量计算

目的	计算本征分子扩散率张量。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Axis-Connectivity（如果步骤4未完成，则为10mc3_200.vol.Labels），并选择XLab Simulations > Molecular Diffusivity Tensor Calculation。
• 把该ROI Box连接到Molecular Diffusivity Tensor Calculation模块的ROI输入连接上。
• 如果步骤4已完成，请在Pore space端口中勾选指向标签1的复选框，因为10mc3_200.vol.Axis-Connectivity中孔隙空间的标签是1。如果步骤4未完成，请在Pore space中勾选指向标签0的复选框，因为10mc3_200.vol.Labels中孔隙空间的标签是0。
• 单击Apply。

使用这些参数，该模块将计算完整的本征扩散率张量。一次完整的张量计算需要三次计算，每一次在时间与内存消耗上都等同于一次实验仿真。

请注意，concentration field输出默认未被选中。该输出对应于$\vec{b}$矢量，即用体积平均方法从菲克方程导出的矢量问题的解（见分子扩散仿真理论页）。它被用于有效分子扩散率的计算，但它的可视化很难解释。

14.2.8.1 检索并解释结果

Project View中只生成了电子表格10mc3_200.vol.DTensor.Spreadsheet这一个输出。该电子表格可以通过在Project View中选中它并单击Show按钮（在Properties area中）来可视化。见图14.35。它包含关于计算结果的信息，汇集在两个表格中。

• 数据集名称，
• 计算所在的感兴趣区域，
• 以3×3矩阵形式给出的完整分子扩散率张量（无量纲），
• 该张量的特征系统解。特征值及其所关联的特征向量被描述在一行中。

电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

图 14.35：包含Molecular Diffusivity Tensor Calculation计算主要结果的电子表格中的各表。

14.2.9 步骤8 - 分子扩散率计算的验证

目的	使用实验结果与经验定律检查所计算分子扩散率的准确性。
必做步骤	否。

我们的验证基于若干项研究——它们报告了分子扩散率实验测量的结果，以及旨在计算材料分子扩散率的经验定律。

Vrettos等人(1989)报告了Carman(1956)在不同类型样品（玻璃珠堆、白砂、玻璃砂、石英砂、砂）上的测量结果。

Kim等人(1987)研究了玻璃球填充床的各向同性系统，并报告了Currie(1960)和Hoogschagen(1955)的测量结果，与他自己的测量作了比较。

Saez等人(1991)和Whitaker(1999)报告了分子扩散率相对于孔隙率$\epsilon$的以下解析估计：

• Maxwell (1881)：$\epsilon\frac{D}{D_{solution}} = \frac{2\epsilon}{3-\epsilon}$
• Weissberg (1963)：$\epsilon\frac{D}{D_{solution}} = \frac{\epsilon}{1-ln(\epsilon)/2}$
• Torquato (1985)：$\epsilon\frac{D}{D_{solution}} = \frac{\epsilon-0.5\epsilon\xi}{1.5-0.5\epsilon-0.5\epsilon\xi}$，其中$\xi = 0.21068(1-\epsilon) - 0.04693(1-\epsilon)^2 + 0.00247(1-\epsilon)^3$。

最后，Vrettos(1989)和Quintard(1993)在Voronoi网络上以及立方构型上得到了数值结果。

所有这些结果都显示在图14.36中。

我们把这些值与用Avizo XLabSuite Extension模块在本教程所用数据集（一堆随机堆积的球体）上计算的值作比较。这里我们使用数据集完整体、高分辨率版本的完整几何形状（也就是从data/tutorials/xlab/10mc3_400.vol.am得到的标签图像），不使用任何ROI Box。

孔隙率可以用Avizo估计（把Measure And Analyze > Global Measures > Intensity Integral in 3D interpretation连接到10mc3_400.vol.Axis-Connectivity上），其值为36.35%。

每个方向上的实验仿真给出以下结果：

• $\epsilon\frac{DX}{D{solution}} = 0.249$
• $\epsilon\frac{DY}{D{solution}} = 0.255$
• $\epsilon\frac{DZ}{D{solution}} = 0.257$

所计算的有效分子扩散率张量为：

$$\epsilon\frac{\overrightarrow{\overrightarrow{D}}}{D_{solution}} = \begin{bmatrix} 0.213 & 0.002 & -0.003 \\ 0.002 & 0.213 & -0.001 \\ -0.003 & -0.001 & 0.209 \end{bmatrix}$$

这些值也显示在图14.36中。可以看到Avizo XLabSuite Extension的结果与各项研究的结果吻合良好。

图 14.36：实验测量、经验定律与数值仿真（其中包括Avizo XLabSuite Extension模块）在分子扩散率随孔隙率变化的测定上的比较。

还提供了第二个演示脚本（仅适用于Microsoft Windows），它在完整分辨率数据集上、限制在一个感兴趣区域内运行本教程的各个步骤。

参考文献：

1. Carman P.C., Flow of Gases through Porous Media, Butterworths, London, 1956
2. Currie J.A., Gaseous diffusion in porous media. Part I - a non-steady state method, Brit. J. Appl. Phys., II, 314-324, 1960
3. Hoogschagen J., Diffusion in porous catalysts and absorbents, Ind. Eng. Chem., 47, 906-913, 1955
4. Kim J.-H., Ochoa J.A., Whitaker S., Diffusion in Anisotropic Porous Media, Transport in Porous Media, 2, 327-356, 1987
5. Maxwell J.C., Treatise on Electricity and Magnetism, Vol. I, 2nd edn., Clarendon Press, Oxford, 1881
6. Quintard M., Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations, Transport in Porous Media, 11, 187-199, 1993
7. Saez A.E., Perfetti J.C., Rusinek I., Prediction of Effective Diffusivities in Porous Media using Spatially Periodic Models, Transport in Porous Media, 6, 143-157, 1991
8. Torquato S., Effective electrical conductivity of two-phase disordered composite media, J. Appl. Phys., 58, 3790-3797, 1985
9. Vrettos N.A., Imakoma H., Okazaki M., Transport Properties of Porous Media from the Micro-geometry of a Three-dimensional Voronoi Network, Chem. Eng. Process, 26, 237-246, 1989
10. Weissberg H.L., Effective diffusion coefficients in porous media, J. Appl. Phys., 34, 2636-2639, 1963
11. Whitaker S., The method of volume averaging, Theory and applications of transport in porous media Vol. 13, Kluwer Acad. Pub., Dordrecht, 1999

14.3 形成因子与电导率计算入门

这个分步教程的目的，是更熟悉Avizo XLabSuite Extension为Avizo所提供的电导率与形成因子计算模块的用法。本章将涉及以下主题：

• 关于电导率的理论基础
• 为仿真准备数据（步骤1到步骤5）
• 实验仿真（步骤6）：定义参数、运行仿真、解释并可视化结果
• 有效属性计算（步骤7）：定义参数、运行仿真、解释并可视化结果
• 用若干经验定律验证结果（步骤8）。

由于为仿真准备数据这一主题已在前一个关于绝对渗透率计算的教程中讲过，这里只作简要提及：

• 设定体素尺寸与单位；
• 分割孔隙空间；
• 用Avizo移除非渗流空间；
• 定义感兴趣区域。

还提供了一个演示脚本（仅适用于Microsoft Windows）。该脚本会自动执行本教程中所详述的全部步骤。

14.3.1 理论要素

欧姆定律：电传导的定义

电传导是这样一个过程：电荷在电场的作用下被传输。在电解质中，电能由自由离子传递。均质材料中的电传导由欧姆定律描述：

$$\vec{j} = -\sigma\vec{\nabla}v$$

其中：

• $\vec{j}$是电流密度（单位：$A \cdot m^{-2}$），
• $\sigma$是材料的电导率（单位：$S \cdot m^{-1}$或$A \cdot V^{-1} \cdot m^{-1}$），
• $v$是电势（单位：$V$）。

多孔材料中的欧姆定律

在多孔材料中，电传导会因电解质周围材料的存在而被改变。一般来说，在形成因子最常见的应用中，固相被视为绝缘体，因为它的电导率比电解质电导率低若干个数量级。我们考虑这样一种多孔介质：

• 均质的：只有一个固相，
• 饱和的：材料的孔隙空间充满溶剂。

由于我们考虑的是导电溶液，不存在电荷积累。欧姆定律与电荷守恒导出以下方程：

$$\vec{\nabla} \cdot \vec{j} = 0$$

又由于我们只考虑一个均质流体相，电导率$\sigma$在空间中不变化，这导出：

$$\nabla^2 v = 0$$

在Formation Factor Experiment Simulation模块中所仿真的实验，是在材料样品的两个相对面之间施加一个恒定的电势差（使用直流电）。样品的其他面被电绝缘体包围。当达到最终状态时，输入与输出电流通量相等；通过在整个体上应用欧姆定律，表观电导率被估计为：

$$\frac{j_{total}}{S} = \sigma\frac{V_{in} - V_{out}}{L}$$

其中

• $j_{total}$是通过输入面的总电通量，
• $S$是输入面的面积，
• $\sigma$是材料的电导率，
• $V{in}$和$V{out}$是所施加的输入与输出电势，
• $L$是材料样品的长度。

通过输入面的总电通量$j_{total}$可以通过局部应用欧姆定律来计算：

$$j_{total} = \int_S -\sigma_{solution}\vec{\nabla}v \cdot \vec{ds}$$

其中$\sigma_{solution}$是自由溶液的电导率。

重要说明：在该实验中，我们考虑的是直流电（没有瞬态现象），并且我们忽略了皮肤效应和浓度变化。实验者通常必须使用校正模型来消除那些实验现象，而我们在这里用理想条件下的实验仿真规避了它们。

形成因子通过电导率的倒数——电阻率——与电导率直接相关。确实，它是”充满水的岩石的电阻率与该水的电阻率之比”（Schlumberger Oilfield Glossary: formation factor, 2009）。

欧姆定律的体积平均形式

有效电导率张量给出关于材料电传导能力的全局信息。为了得到在整个体上有效的方程，需要进行一次尺度变换。体积平均法就是完成这种尺度变换的一种技术。它的主要目标是通过在一个体上对方程做平均，使方程在空间上变得平滑。感兴趣的读者可以在[Whitaker, S., The Method of Volume Averaging, Kluver Academic Publishers, 1999]中找到相关细节和参考文献。

这个非常一般的理论导出一个闭合问题，它把欧姆方程变换为一个矢量问题。闭合变量$\vec{b}$被用于在一个新问题中表达电势扰动：

$$\nabla^2\vec{b} = 0$$

当该问题被求解后，就可以计算如下定义的无量纲电导率张量：

$$\epsilon\frac{\overrightarrow{\overrightarrow{\sigma}}}{\sigma_{solution}} = \epsilon\left(\overrightarrow{\overrightarrow{I}} + \frac{1}{V_f}\int_{S_{fs}} \overrightarrow{n_{fs}}\,\vec{b}\,ds\right)$$

其中：

• $\epsilon$是孔隙率，
• $\sigma_{solution}$是溶液的电导率，
• $V_f$是流体的体积，
• $S_{fs}$是流固界面的面积，
• $\overrightarrow{n_{fs}}$是流固界面的法线，方向从流体指向固相。

电导率张量的逆就是形成因子张量，而形成因子标量是后一个张量各特征值的平均值。所有这些都由Effective Formation Factor Calculation模块计算。

边界条件

Avizo XLabSuite Extension提出两种估计电导率的思路。

第一种是基于欧姆方程求解的实验仿真。如前所述，这由Formation Factor Experiment Simulation模块完成。材料一般由岩石构成，可以被视为绝缘体。于是流固界面处的边界条件为：

$$-\overrightarrow{n_{fs}} \cdot \vec{\nabla}v = 0$$

其中$\overrightarrow{n_{fs}}$是流固界面的法线，方向从流体指向固相。除了这个流固界面条件之外，边界条件还有：

• 在图像中不垂直于电通量主方向的那些面上，添加一个一体素宽的电绝缘体平面。这使样品与外界隔离。
• 输入与输出（即垂直于通量主方向的那些面）被设计为施加电势的一体素宽平面。

第二种思路求解由欧姆方程通过体积平均导出的闭合问题。这由Effective Formation Factor Calculation模块完成。这种情况下所求解的矢量问题，是通过对$\vec{b}$和几何形状施加周期性边界条件而闭合的。流固界面条件具有以下类似形式：

$$-\overrightarrow{n_{fs}} \cdot \vec{\nabla}\vec{b} = \overrightarrow{n_{fs}}$$

方程组的离散化

Avizo XLabSuite Extension使用有限体积法来求解这些方程组。

该离散化方案假定体素是各向同性的（立方体）。

系统求解

一旦离散化，这些方程组可以写成$Ax = b$，其中$A$是一个稀疏对称矩阵。

这些方程组使用完全隐式方法（矩阵求逆）求解。线性系统的直接求解使用PETSc（Portable, Extensible Toolkit for Scientific Computation）库。

同时也执行一次带共轭梯度与ILU预条件子的迭代求解。所使用的收敛判据是残差$l_2$范数的相对下降。

14.3.2 步骤1 - 加载数据集

目的	在Avizo中加载数据集。
必做步骤	是。

本教程所使用的数据集是一堆随机堆积的玻璃球。它由ICMCB-CNRS研究主任Dominique Bernard扫描。这些球形颗粒经过筛分，直径在100–120微米范围内。该样品在扫描之前于700摄氏度下烧结了10分钟。图14.37展示了该数据集。

本教程所使用的完整3D数据集是一个200×200×200的立方体。

• 开始一个新项目 File > New Project
• 打开File > Open Data...对话框。
• 从AVIZO_ROOT目录打开data/tutorials/xlab/10mc3_200.vol.am文件。
• 单击Open。

图 14.37：本教程所用球体堆积的3D可视化。

如前一步所述，单位管理在这里不是必做的。不过，请参阅绝对渗透率仿真教程的步骤2以了解更多关于编辑长度单位的内容。

14.3.3 步骤2 - 设定体素尺寸

目的	设定均匀数据集中体素的边长。
必做步骤	否。

同样，如前面各步骤所述，正确设定体素的尺寸在这里不是必做的。不过，要了解更多关于使用Crop Editor设定体素尺寸的内容，请参阅绝对渗透率仿真教程的步骤3。

注意：Avizo XLabSuite Extension的计算必须使用各向同性（即立方体）体素。

14.3.4 步骤3 - 从数据集创建标签图像

目的	使用基本分割工具从数据集创建一个标签图像。
必做步骤	是。

这一步是使用Avizo所提供的一些分割工具，从被平滑并阈值化的初始数据集创建一个标签图像。

请遵照绝对渗透率仿真教程步骤4中所述的说明操作。

图 14.38：阈值化后数据集的可视化。

14.3.5 步骤4 - 移除非渗流的孔隙空间

目的	使用高级分割工具移除阈值化图像中非渗流的孔隙空间。
必做步骤	否。仅在Avizo中可用。

渗流测试可以在运行计算之前执行。这项测试的一个好处是，它会把所有孤立的孔隙空间从标签图像中移除（我们提醒用户：固体材料被视为电绝缘体，因此电流不会到达孤立空间中的溶液）。这样计算应当会稍快一些。从离散的角度看，所需的连通性是6连通：电流密度只能在非零表面上被估计，也就是说在两个相邻体素之间的一个面上。

你可以遵照绝对渗透率仿真教程步骤5中所述的说明来执行渗流测试。

图 14.39：数据集中渗流孔隙率的可视化。

14.3.6 步骤5 - 子区域的选择

目的	为电导率计算定义一个感兴趣区域。
必做步骤	是。

有一种简便方法可以在所加载数据集的某个子区域上计算电导率，而不必裁剪它。可以定义一个感兴趣区域(ROI)并把它连接到Avizo XLabSuite Extension的模块上，这样计算就只在该ROI中进行。

请遵照绝对渗透率仿真教程步骤6中给出的说明，定义一个以样品为中心、边长为50个体素的ROI立方体（见图14.40）。本教程将用该ROI在小体积上计算属性。在教程过程中，它可以随意在数据集中移动，不过为了便于说明，将使用它的初始位置。

为方便起见，同时也因为这会在本教程阶段节省时间，我们在样品的一个子区域上计算电导率和形成因子。不过，我们将在验证阶段看到重采样和使用ROI如何影响所计算结果的质量。

图 14.40：代表感兴趣区域的立方体的可视化。

14.3.7 步骤6 - 形成因子实验仿真

目的	仿真一次形成因子实验室测量。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Axis-Connectivity（如果步骤4未完成，则为10mc3_200.vol.Labels），并选择XLab Simulations > Formation Factor Experiment Simulation。
• 把该ROI Box连接到Formation Factor Experiment Simulation模块的ROI输入连接上。
• 如果步骤4已完成，请在Pore space端口中勾选指向标签1的复选框，因为10mc3_200.vol.Axis-Connectivity中孔隙空间的标签是1。如果步骤4未完成，请在Pore space中勾选指向标签0的复选框，因为10mc3_200.vol.Labels中孔隙空间的标签是0。
• 单击Apply。

默认参数仿真的是沿Z轴、输入电势为1 V、输出电势为零的一个实验。默认的溶液电导率是0.0001 $S \cdot m^{-1}$。这些选项可以被修改以仿真若干不同的实验。例如，电通量的方向可以被调整为X、Y或Z方向（默认为Z）。如果选择了若干个方向，计算会被依次执行。用作实验边界条件的电势值也可以被修改。修改这些值不会改变电导率或形成因子——它们是该多孔介质的本征属性。它只会修改输出的电势场。

14.3.7.1 检索并解释结果

Project View中出现了两个输出：

• 10mc3_200.vol.FFExp.Spreadsheet：一个电子表格，包含计算最相关的结果：
  • 数据集名称；
  • 计算所在的感兴趣区域；
  • 以$S \cdot m^{-1}$为单位的表观电导率；
  • 表观形成因子；
  • 以V为单位、施加在实验装置输入端的电势；
  • 以V为单位、施加在实验装置输出端的电势；
  • 以$S \cdot m^{-1}$为单位的溶液电导率。
• 10mc3_200.vol.PotentialZ：一个标量场，表示在受ROI限制的几何形状中所求解欧姆方程组的解——电势场（单位：V）。

包含全部相关信息的电子表格，可以通过在Project View中选中它并单击Show按钮（在Properties area中）来可视化。见图14.41。电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

图 14.41：包含Formation Factor Experiment Simulation计算主要结果的电子表格。

14.3.7.2 可视化输出电势场

要可视化电势场：

• 通过单击它们的查看器开关，隐藏当前所有的显示模块，例如Bounding Box或Ortho Slice。
• 右键单击10mc3_200.vol.PotentialZ并选择Display > Ortho Slice。
• 在Ortho Slice属性中，从Edit菜单把所选颜色图改为temperature.icol。
• 再次选择颜色图的Edit菜单，然后选择Adjust range to > Data min-max。
• 在Orientation端口中选择yz。
• 放大该切片以可视化它。

所得到的可视化应当看起来像图14.42。你可以观察到电势从装置的输入端（在切片底部，黄色）到输出端（在切片顶部，浅蓝色）的下降。

图 14.42：Z方向电通量实验仿真中电势场的可视化。

14.3.8 步骤7 - 有效形成因子计算

目的	计算本征电导率张量。计算本征形成因子。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Axis-Connectivity（如果步骤4未完成，则为10mc3_200.vol.Labels），并选择XLab Simulations > Effective Formation Factor Calculation。
• 把该ROI Box连接到Effective Formation Factor Calculation模块的ROI输入连接上。
• 如果步骤4已完成，请在Pore space端口中勾选指向标签1的复选框，因为10mc3_200.vol.Axis-Connectivity中孔隙空间的标签是1。如果步骤4未完成，请在Pore space中勾选指向标签0的复选框，因为10mc3_200.vol.Labels中孔隙空间的标签是0。
• 单击Apply。

使用这些参数，该模块将计算完整的本征电导率张量。一次完整的张量计算需要三次计算，每一次在时间与内存消耗上都等同于一次实验仿真。

请注意，potential field输出默认未被选中。该输出对应于$\vec{b}$矢量，即用体积平均方法从欧姆方程导出的矢量问题的解（见电导率仿真理论页）。它被用于有效电导率的计算，但它的可视化很难解释。

14.3.8.1 检索并解释结果

Project View中只生成了电子表格10mc3_200.vol.FFTensor.Spreadsheet这一个输出。该电子表格可以通过在Project View中选中它并单击Show按钮（在Properties area中）来可视化。见图14.43。它包含关于计算结果的信息，汇集在两个表格中。

• 电导率张量表：
  • 数据集名称；
  • 计算所在的感兴趣区域；
  • 以3×3矩阵形式给出的完整电导率张量（无量纲）；
  • 该张量的特征系统解。特征值及其所关联的特征向量被描述在一行中。
• 形成因子张量表：
  • 数据集名称；
  • 计算所在的感兴趣区域；
  • 以3×3矩阵形式给出的完整形成因子张量（无量纲）；
  • 该张量的特征系统解。特征值及其所关联的特征向量被描述在一行中。

由于形成因子张量是作为电导率张量的逆来计算的，因此只有在三个方向全部被选中、完整电导率张量已被计算的情况下，第二个表格才会被生成。

电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

图 14.43：包含Effective Formation Factor Calculation计算主要结果的电子表格中的各表。

14.3.9 步骤8 - 形成因子计算的验证

目的	使用若干经验定律检查所计算形成因子的准确性。
必做步骤	否。

我们的验证基于Lemaitre等人(1988)的研究。他们测量了由密集堆积的二元球体混合物构成、并充满导电流体的多孔介质的电导率（并由此推导出形成因子）。该实验是在具有两种球直径比和若干孔隙率值的混合物上进行的。

结果表明，形成因子几乎只依赖于孔隙率，而不依赖于堆积结构。

他们还回顾了众多试图把形成因子F与孔隙率$\epsilon$联系起来的作者所得到的主要结果：

• Archie定律 (1942)：$F = A\epsilon^{-m}$，其中A通常被假定为1，而m取决于多孔介质的几何形状（m通常属于范围$[1,3]$），
• Sen等人 (1981)：$F = \epsilon^{-3/2}$，对于$\epsilon > 0.2$的烧结玻璃球堆积是一个良好的近似。
• Berryman (1983)：$F = 2/(\epsilon(1+\epsilon))$，它很好地适用于孔隙率很大的材料。

结果显示，Lemaitre等人在玻璃球混合物上所测得的实验值，落在由Sen等人和Berryman的理论所得值定义的范围内，并且能被m=1.46的Archie定律很好地拟合。

本教程所使用的数据集是随机堆积的球体。这里我们使用数据集完整体、高分辨率版本的完整几何形状（也就是从data/tutorials/xlab/10mc3_400.vol.am得到的标签图像），不使用任何ROI Box。

孔隙率可以用Avizo估计（把Measure And Analyze > Global Measures > Intensity Integral in 3D interpretation连接到10mc3_400.vol.Axis-Connectivity上），其值为36.35%。

由不同经验定律计算出的形成因子值为：

• Sen等人：F = 4.563
• Berryman：F = 4.035

我们把这些值与在data/tutorials/xlab/10mc3_400.vol.am上用电导率仿真模块所计算的值作了比较。

每个方向上的实验仿真给出以下结果：

• $F_X = 4.599$
• $F_Y = 4.602$
• $F_Z = 4.623$

所计算的有效形成因子为F = 4.714。

注意：我们可以注意到，与在降采样情形10mc3_200.vol.am的某个ROI上所得到的值相比，存在相当大的差别。这显示了重采样和使用ROI如何影响所计算结果的质量。

仿真的实验值和有效值全都与Sen等人及Berryman的理论所得的值处于同一数量级。如果我们应用Archie定律，所计算的值对应于在1.51到1.53之间变化的m系数，这是非常合理的。

还提供了第二个演示脚本（仅适用于Microsoft Windows），它在完整分辨率数据集上、限制在一个感兴趣区域内运行本教程的各个步骤。

参考文献：

1. Archie G.E., The electrical resistivity log as an aid in determining some reservoir characteristics, Petroleum Transactions of AIME, 146, 54-62, 1942
2. Berryman J.G., Effective conductivity by fluid analogy for a porous insulator filled with a conductor, Phys. Rev. B 27, 7789-7792, 1983
3. Lemaitre J., Troadec J.P., Bideau D., Gervois A. and Bougault E., The formation factor of the pores space of binary mixtures of spheres, J. Phys. D: Appl. Phys., 21, 1589-1592, 1988
4. Sen P.N., Scala C. and Cohen M.H., A self-similar model from sedimentary rocks with application to dielectric constant of fused glass beads, Geophysics, 46, 781-795, 1981

14.4 热导率计算入门

这个分步教程的目的，是更熟悉Avizo XLabSuite Extension为Avizo所提供的热导率计算模块的用法。本章将涉及以下主题：

• 关于热导率的理论基础
• 为仿真准备数据（步骤1到步骤5）
• 实验仿真（步骤6）：定义参数、运行仿真、解释并可视化结果
• 有效属性计算（步骤7）：定义参数、运行仿真、解释并可视化结果
• 用若干经验定律验证结果（步骤8）。

由于为仿真准备数据这一主题已在前一个关于绝对渗透率计算的教程中讲过，这里只作简要提及：

• 设定体素尺寸与单位，
• 分割孔隙空间，
• 用Avizo移除非渗流空间，
• 定义感兴趣区域。

还提供了一个演示脚本（仅适用于Microsoft Windows）。该脚本会自动执行本教程中所详述的全部步骤。

14.4.1 理论要素

傅里叶定律：热传导的定义

热传导是材料把热量从高温区域传导到低温区域的能力。在稳态条件下，均质材料中的热传导由傅里叶定律描述：

$$\vec{\varphi} = -\lambda\vec{\nabla}T$$

其中：

• $\vec{\varphi}$是热通量（单位：$W \cdot m^{-2}$），
• $\lambda$是材料的热导率（单位：$W \cdot m^{-1} \cdot K^{-1}$），
• $T$是温度（单位：$K$）。

非均质材料中的傅里叶定律

大多数材料是非均质的，并且包含多于一个固相。因此，这类材料中的热传导比傅里叶定律所表达的更为复杂。

在某个名为$\alpha$的给定相中的瞬态热传导由以下偏微分方程描述：

$$(\rho c_p)_\alpha\frac{\partial T_\alpha}{\partial t} - \lambda_\alpha\nabla^2 T_\alpha = 0$$

其中：

• $(\rho cp)\alpha$是$\alpha$相的热容（单位：$J \cdot m^{-3} \cdot K^{-1}$）；
• $\rho_\alpha$是$\alpha$相的密度（单位：$kg \cdot m^{-3}$）；
• $c{p\alpha}$是$\alpha$相的比热（单位：$J \cdot kg^{-1} \cdot K^{-1}$）；
• $\lambda_\alpha$是$\alpha$相的热导率（单位：$W \cdot m^{-1} \cdot K^{-1}$）；
• $T_\alpha$是$\alpha$相的温度（单位：$K$）。

我们认为已达到稳态。那么我们现在要求解的方程是：

$$\lambda_\alpha\nabla^2 T_\alpha = 0$$

在Thermal Conductivity Experiment Simulation模块中所仿真的实验，是在材料样品的两个相对面之间施加一个恒定的热通量。例如，输入与输出温度可以通过输入端的一个加热电阻和输出端的一个冷冻浴来保持恒定。样品的其他面是理想的热绝缘平面。当达到最终状态时，输入与输出热通量相等，傅里叶定律在整个体上写作：

$$\frac{\varphi_{total}}{S_{in}} = \lambda\frac{T_{in} - T_{out}}{L}$$

其中

• $\varphi_{total}$是通过输入面的总热通量，
• $S_{in}$是输入面的面积，
• $\lambda$是材料的表观热导率，
• $T{in}$和$T{out}$是所施加的输入与输出温度，
• $L$是材料样品的长度。

只要我们知道其他各项，表观热导率$\lambda$就可以借助前面的表达式计算出来。实验者控制温度$T{in}$和$T{out}$，而通过输入面的总热通量很容易用局部的傅里叶定律确定：

$$\varphi_{total} = \int_{S_{in}} -\lambda_\alpha\vec{\nabla}T_\alpha \cdot \vec{ds}$$

其中$\lambda\alpha$和$T\alpha$是与输入表面$S_{in}$接触的材料任意相$\alpha$的热导率和温度。

应用于傅里叶定律的多尺度均质化理论

有效热导率张量给出关于材料热传导能力的全局信息。

均质化理论在于同时在一个宏观域和一个微观域上考虑该问题。这里，宏观域是周期性域，具有特征长度$L$；而微观域是那个可被视为代表性基元体(R.E.V.)的周期。特征局部长度$l$是该周期的长度，且$L \gg l$。引入特征无量纲空间变量$x^$和$y^$，使得$y^ = X/l$、$x^ = X/L$，其中$X$是物理空间变量。

由于存在这两个特征变量，空间导数变为$\nabla{x^*} + \epsilon^{-1}\nabla{y^}$，其中$\epsilon = x^/y^*$，因此$\epsilon \ll 1$。

未知量$T$被写为相对于$\epsilon = x^/y^$的渐近展开。

考虑新的空间导数和$T$的渐近展开，傅里叶偏微分方程被重写。最后，把$\epsilon$的同次幂各项一一对应，导出对若干个相继问题的求解，其中在R.E.V.上有以下正则方程：

$$\nabla(\lambda_\alpha(\vec{\nabla}\vec{b}_\alpha + \overrightarrow{\overrightarrow{I}})) = 0$$

其中$\vec{b}$可被视为温度场的一个扰动。$\vec{b}$还满足：

$$\overrightarrow{\overrightarrow{\lambda_{eff}}} = \frac{1}{V}\sum_\alpha \int_{V_\alpha} \lambda_\alpha\left(\vec{\nabla}\vec{b}_\alpha + \overrightarrow{\overrightarrow{I}}\right)dv$$

其中：

• $\overrightarrow{\overrightarrow{\lambda_{eff}}}$是有效热导率张量，
• $V$是样品的总体积，
• $\alpha$是一个导热相，
• $V_\alpha$是每个相$\alpha$所占据的体积。

感兴趣的读者可以在[Auriault, J.-L., Boutin, C. and Geindreau, C., Homogénéisation de phénomènes couplés en milieux hétérogènes 1, Lavoisier, 2009]中找到相关细节和参考文献。

边界条件

如前所述，大多数材料不是均质的，并且包含多于一个相。每个相的特性可以完全不同，材料各组分之间的界面在整体热导率中起着主要作用。

我们来考虑一种两相材料，由一个$\alpha$相和一个$\beta$相构成。这可以推广到任意数量的相。必须为各相之间的每一个界面定义边界条件。于是待求解的系统为：

$$\begin{cases} \lambda_\xi\nabla^2 T_\xi = 0 & \text{对于}\ \xi = \alpha, \beta \\ T_\alpha = T_\beta & \text{在}\ \alpha - \beta\ \text{界面处} \\ -\overrightarrow{n_{\alpha\beta}} \cdot \lambda_\alpha\vec{\nabla}T_\alpha = -\overrightarrow{n_{\alpha\beta}} \cdot \lambda_\beta\vec{\nabla}T_\beta & \text{在}\ \alpha - \beta\ \text{界面处} \end{cases}$$

其中$\overrightarrow{n_{\alpha\beta}}$是垂直于两相之间界面、方向从$\alpha$指向$\beta$的单位矢量。

这些边界条件规定了温度和热通量的法向分量在各相之间的界面处是连续的。

Avizo XLabSuite Extension提出两种估计热导率的思路。

第一种是基于傅里叶方程求解的实验仿真。如前所述，这由Thermal Conductivity Experiment Simulation模块完成。除了相界面条件之外，边界条件还有：

• 在图像中不垂直于热通量主方向的那些面上，添加一个一体素宽的热绝缘体平面。这使样品与外界隔离。
• 输入与输出（即垂直于通量主方向的那些面）被设计为施加温度的一体素宽平面。
• 以下三个条件中的两个可由用户施加，第三个则从所选的两个估计得出：输入温度、输出温度、热通量。

第二种思路求解在一个无限周期性域上、由傅里叶方程通过均质化导出的正则问题。这由Thermal Conductivity Tensor Calculation模块完成。两个边界条件源自两相界面处温度以及热通量法向分量的连续性条件：

$$\begin{cases} \vec{b}_\alpha = \vec{b}_\beta \\ -\lambda_\alpha\overrightarrow{n_{\alpha\beta}} \cdot (\vec{\nabla}\vec{b}_\alpha + \overrightarrow{\overrightarrow{I}}) = -\lambda_\beta\overrightarrow{n_{\alpha\beta}} \cdot (\vec{\nabla}\vec{b}_\beta + \overrightarrow{\overrightarrow{I}}) \end{cases}$$

对$\vec{b}_\alpha$和几何形状施加周期性边界条件。

方程组的离散化

Avizo XLabSuite Extension使用有限体积法来求解这些方程组。

该离散化方案假定体素是各向同性的（立方体）。

系统求解

一旦离散化，这些方程组可以写成$Ax = b$，其中$A$是一个稀疏对称矩阵。

这些方程组使用完全隐式方法（矩阵求逆）求解。线性系统的直接求解使用PETSc（Portable, Extensible Toolkit for Scientific Computation）库。

同时也执行一次带共轭梯度与ILU预条件子的迭代求解。所使用的收敛判据是残差$l_2$范数的相对下降。

14.4.2 步骤1 - 加载数据集

目的	在Avizo中加载数据集。
必做步骤	是。

本教程所使用的数据集是一堆随机堆积的玻璃球。它由ICMCB-CNRS研究主任Dominique Bernard扫描。这些球形颗粒经过筛分，直径在100–120微米范围内。该样品在扫描之前于700摄氏度下烧结了10分钟。图14.44展示了该数据集。

本教程所使用的完整3D数据集是一个200×200×200的立方体。

• 开始一个新项目 File > New Project
• 打开File > Open Data...对话框。
• 从AVIZO_ROOT目录打开data/tutorials/xlab/10mc3_200.vol.am文件。
• 单击Open。

图 14.44：本教程所用球体堆积的3D可视化。

如前一步所述，单位管理在这里不是必做的。不过，请参阅绝对渗透率仿真教程的步骤2以了解更多关于编辑长度单位的内容。

14.4.3 步骤2 - 设定体素尺寸

目的	设定均匀数据集中体素的边长。
必做步骤	否。

同样，如前面各步骤所述，正确设定体素的尺寸在这里不是必做的。不过，要了解更多关于使用Crop Editor设定体素尺寸的内容，请参阅绝对渗透率仿真教程的步骤3。

注意：Avizo XLabSuite Extension的计算必须使用各向同性（即立方体）体素。

14.4.4 步骤3 - 从数据集创建标签图像

目的	使用基本分割工具从数据集创建一个标签图像。
必做步骤	是。

这一步是使用Avizo所提供的一些分割工具，从被平滑并阈值化的初始数据集创建一个标签图像。

请遵照绝对渗透率仿真教程步骤4中所述的说明操作。所得到的可视化应当看起来像图14.45。

图 14.45：阈值化后数据集的可视化。

14.4.5 步骤4 - 移除孤立的导热材料部分

目的	使用高级分割工具移除被绝缘材料包围的导热材料部分。
必做步骤	否。仅在Avizo中可用。

如果材料的某一个或若干个相被视为热绝缘体（且仅在这种情况下），强烈建议运行一项测试，以检测并移除完全被绝缘体包围的导热材料部分。被移除的部分将被视为绝缘体，如图14.46所示。

这会避免为传输方程求解构建出奇异矩阵的风险。该测试可以在运行计算之前执行。从离散的角度看，所需的连通性是6连通：热量只能通过非零表面传输，也就是说通过两个相邻体素之间的一个面。

这项测试背后的思想与绝对渗透率仿真教程步骤5中所执行的非渗流测试的思想相同。因此，同样的工具可以用于执行这两项测试。

在本教程中，我们将认为材料的所有相都导热；但如果你遇到必须同时考虑绝缘相与导热相的情形，你可以遵照绝对渗透率仿真教程步骤5中所述的说明操作。绝对渗透率仿真教程中流体流动的孔隙空间，在这里对应于导热材料；而多孔材料在这里对应于绝缘体。

图 14.46：移除属于导热材料但被绝缘材料包围的体素。

14.4.6 步骤5 - 子区域的选择

目的	为热导率计算定义一个感兴趣区域。
必做步骤	是。

有一种简便方法可以在所加载数据集的某个子区域上计算热导率，而不必裁剪它。可以定义一个感兴趣区域(ROI)并把它连接到Avizo XLabSuite Extension的模块上，这样计算就只在该ROI中进行。

请遵照绝对渗透率仿真教程步骤6中给出的说明，定义一个以样品为中心、边长为50个体素的ROI立方体（见图14.47）。本教程将用该ROI在小体积上计算属性。在教程过程中，它可以随意在数据集中移动，不过为了便于说明，将使用它的初始位置。

为方便起见，同时也因为这会在本教程阶段节省时间，我们在样品的一个子区域上计算热导率。不过，我们将在验证阶段看到重采样和使用ROI如何影响所计算结果的质量。

图 14.47：代表感兴趣区域的立方体的可视化。

14.4.7 步骤6 - 热导率实验仿真

目的	仿真一次热导率实验室测量。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Labels，并选择XLab Simulations > Thermal Conductivity Experiment Simulation。
• 把该ROI Box连接到Thermal Conductivity Experiment Simulation模块的ROI输入连接上。
• 在Conducting Materials端口中，勾选这两个相的复选框。
• 在Thermal Conductivity端口中把Inside相的热导率设为1000 $W \cdot m^{-1} \cdot K^{-1}$，并把Exterior相（对应于孔隙空间）的热导率保留为默认值1 $W \cdot m^{-1} \cdot K^{-1}$。
• 单击Apply。

默认参数仿真的是沿Z轴、输入温度为298 K、输出温度为273 K的一个实验。这些选项可以被修改以仿真若干不同的实验。例如，热通量的方向可以被调整为X、Y或Z方向（默认为Z）。如果选择了若干个方向，计算会被依次执行。实验的边界条件也可以被修改。三者之中必须施加两个值：输入温度、输出温度、热通量。修改这些值不会改变热导率——它是该材料的本征属性。它只会修改输出的温度场。

14.4.7.1 检索并解释结果

Project View中出现了两个输出：

• 10mc3_200.vol.TCExp.Spreadsheet：一个电子表格，包含计算最相关的结果：
  • 数据集名称；
  • 计算所在的感兴趣区域；
  • 以$W \cdot m^{-1} \cdot K^{-1}$为单位、限制在该ROI内的材料表观热导率；
  • 以K为单位、实验装置输入端的温度；
  • 以K为单位、实验装置输出端的温度；
  • 以$W \cdot m^{-2}$为单位、通过样品的热通量；
  • 以$W \cdot m^{-1} \cdot K^{-1}$为单位、材料各相的热导率。
• 10mc3_200.vol.TemperatureZ：一个标量场，表示在受ROI限制的几何形状中所求解傅里叶方程组的解——温度场（单位：K）。

包含全部相关信息的电子表格，可以通过在Project View中选中它并单击Show按钮（在Properties area中）来可视化。见图14.48。电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

图 14.48：包含Thermal Conductivity Experiment Simulation计算主要结果的电子表格。

14.4.7.2 可视化输出温度场

要可视化温度场：

• 通过单击它们的查看器开关，隐藏当前所有的显示模块，例如Bounding Box或Ortho Slice。
• 右键单击10mc3_200.vol.TemperatureZ并选择Display > Ortho Slice。
• 在Ortho Slice属性中，从Edit菜单把所选颜色图改为temperature.icol。
• 再次选择颜色图的Edit菜单，然后选择Adjust range to > Data min-max。
• 在Orientation端口中选择yz。
• 放大该切片以可视化它。

所得到的可视化应当看起来像图14.49。你可以观察到温度从装置的输入端（在切片底部，黄色）到输出端（在切片顶部，浅蓝色）的下降。你还可以猜到，热锋面在热导率较大的相中推进得更快。

图 14.49：Z方向热通量实验仿真中温度场的可视化。

为了强调这一现象，我们将使用另外一些显示工具。

• 把一个Color Wash模块附加到该Ortho Slice上。
• 在Color Wash模块的Data端口中选择10mc3_200.vol.Labels。
• 把Weight Factor设为0.75。

现在更容易可视化Inside相的颗粒，以及它们较高的热导率对热锋面的影响。

• 把一个Isosurface模块附加到10mc3_200.vol.TemperatureZ上(Display > Isosurface)。
• 勾选Properties area底部的auto-refresh复选框。
• 移动Threshold游标，以可视化不同温度值的等值面。

等值面轮廓的形态突出了这样一个事实：热锋面在Inside相的颗粒中推进得更快。（如果这一点不明显，你可以使用一个Volume Rendering来帮助在3D中可视化这些颗粒。）

• 把一个Isocontour Slice模块附加到10mc3_200.vol.TemperatureZ上(Display > Isocontour Slice)。
• 把Clipping Plane Orientation设为yz，并在Translate端口中把位置设为51。
• 在Isocontour Slice的Properties area中，把Values范围设为273, 298，线的数量设为30。
• 从Parameters端口把line width设为1。

同样，由这些等值线所表示的热锋面，显然在Inside相的颗粒中移动得更快。

所得到的可视化应当看起来像图14.50。

图 14.50：温度等值面T = 284 K与30条温度等值线在带Color Wash渲染的Ortho Slice上的可视化。

14.4.8 步骤7 - 有效热导率计算

目的	计算本征热导率张量。可视化并解释计算结果。
必做步骤	是。

• 右键单击10mc3_200.vol.Labels，并选择XLab Simulations > Thermal Conductivity Tensor Calculation。
• 把该ROI Box连接到Thermal Conductivity Tensor Calculation模块的ROI输入连接上。
• 在Conducting Materials端口中，勾选这两个相的复选框。
• 在Thermal Conductivity端口中把Inside相的热导率设为1000 $W \cdot m^{-1} \cdot K^{-1}$，并把Exterior相（对应于孔隙空间）的热导率保留为默认值1 $W \cdot m^{-1} \cdot K^{-1}$。
• 单击Apply。

使用这些参数，该模块将计算完整的本征热导率张量。一次完整的张量计算需要三次计算，每一次在时间与内存消耗上都等同于一次实验仿真。

请注意，temperature field输出默认未被选中。该输出对应于$\vec{b}$矢量，即用均质化方法从傅里叶方程导出的矢量问题的解（见热导率仿真理论页）。它被用于有效热导率的计算，但它的可视化很难解释。

14.4.8.1 检索并解释结果

Project View中只生成了电子表格10mc3_200.vol.TCTensor.Spreadsheet这一个输出。该电子表格可以通过在Project View中选中它并单击Show按钮（在Properties area中）来可视化。见图14.51。它包含关于计算结果的信息，汇集在两个表格中。

• 数据集名称，
• 计算所在的感兴趣区域，
• 以3×3矩阵形式给出的完整热导率张量（单位：$W \cdot m^{-1} \cdot K^{-1}$），
• 该张量的特征系统解：特征值及其所关联的特征向量被描述在一行中，
• 以$W \cdot m^{-1} \cdot K^{-1}$为单位、材料各相的热导率。

电子表格可以被导出为若干种格式（例如CSV、XML、txt）。

图 14.51：包含Thermal Conductivity Tensor Calculation计算主要结果的电子表格中的各表。

14.4.9 步骤8 - 热导率计算的验证

目的	使用实验结果与经验定律检查所计算热导率的准确性。
必做步骤	否。

我们的验证基于若干旨在计算具有几种非常特定几何形状的材料热导率的经验定律。”内部”材料（或几何体）被称为$\beta$相，具有热导率$\lambda\beta$和体积分数$\epsilon\beta$。包围它的材料被称为$\alpha$相，具有热导率$\lambda\alpha$和体积分数$\epsilon\alpha$。我们定义比值$\Lambda = \lambda\beta/\lambda\alpha$以及$\Lambda{eff} = \lambda{eff}/\lambda\alpha$，其中$\lambda{eff}$是整个材料的有效热导率。

14.4.9.1 不相接触3D圆柱的周期性阵列

我们研究的第一个模型是不相接触3D圆柱的周期性阵列，如图14.52所示。这些圆柱被称为$\beta$相，而包围圆柱的材料被称为$\alpha$相。

图 14.52：不相接触3D圆柱周期性阵列的示例。

在Ochoa-Tapia等人(1994)中，给出了不相接触圆柱阵列在串联分布（即两相相对于热通量方向在热学上串联）情形下有效热导率的以下解析估计：

$$\Lambda_{eff} = \frac{2\Lambda - \epsilon_\alpha(\Lambda-1)}{2 + \epsilon_\alpha(\Lambda-1)}$$

在Perrins等人(1979)中，有效热导率的解析估计被汇集在一张表中，给出了具有各种电导率比$\Lambda$和体积分数的方形圆柱阵列的导率值。

我们把Ochoa-Tapia等人的公式应用到我们的模型上（其中$\epsilon\alpha = 0.49$），并把它与Avizo实验仿真和张量计算所得到的结果，以及Perrins等人对$\epsilon\alpha = 0.5$所给出的热导率值（因为表中没有0.49）作比较：

$\Lambda$	Ochoa-Tapia等人 $\Lambda_{eff}$	Perrins等人 $\Lambda_{eff}$	Avizo TCExp $\Lambda_{eff}$	Avizo TCTensor $\Lambda_{eff}$
2	1.406	1.401	1.397	1.397
3.5	1.783	1.776	1.765	1.765
5	2.020	2.013	1.999	1.999
10	2.416	2.416	2.395	2.395
20	2.692	2.701	2.677	2.677
50	2.897	2.915	2.890	2.890
100	2.972	na	2.970	2.970
1000	3.045	na	3.046	3.045
10000	3.053	na	3.054	3.054
$\infty$	3.053	3.080	na	na

我们观察到解析估计与仿真之间吻合良好。

在Meredith and Tobias (1960)中，给出了不相接触圆柱阵列在并联分布（即两相相对于热通量方向在热学上并联）情形下有效热导率的以下解析估计：

$$\Lambda_{eff} = \epsilon_\alpha + (1-\epsilon_\alpha)\Lambda$$

我们把该公式应用到我们的模型上（其中$\epsilon_\alpha = 0.49$），并把它与Avizo实验仿真和张量计算所得到的结果作比较：

$\Lambda$	Meredith and Tobias $\Lambda_{eff}$	Avizo TCExp $\Lambda_{eff}$	Avizo TCTensor $\Lambda_{eff}$
10	5.559	5.474	5.474
100	51.152	50.218	50.217
1000	507.08	497.65	497.65
10000	5066.4	4972.0	4971.9

同样，我们观察到解析估计与仿真之间吻合良好。

Whitaker (1999)报告了Nozad等人(1985)针对一个简单二维方形单元周期性阵列所计算的结果，并表明所得到的结果与Perrins等人(1979)的结果相似。在图14.53中，我们把Nozad等人对各种$\epsilon\alpha$值的$\Lambda{eff}$相对于$\Lambda$的结果与Avizo仿真作比较，可以看到后一条曲线位置良好，行为正确。由于Avizo实验仿真与张量计算的结果几乎相等，因此只显示了一条曲线。

图 14.53：Avizo仿真与Nozad等人热导率计算的比较。

14.4.9.2 不相接触3D球体的周期性阵列

我们研究的第二个模型是不相接触3D球体的周期性阵列，如图14.54所示。这些球体被称为$\beta$相，而包围球体的材料被称为$\alpha$相。

图 14.54：不相接触3D球体周期性阵列的示例。

Maxwell (1881)对球体充分分离以致互不相互作用的球体填充床，得到了热导率的以下表达式：

$$\Lambda_{eff} = \frac{3\Lambda - 2\epsilon_\alpha(\Lambda-1)}{3 + \epsilon_\alpha(\Lambda-1)}$$

Meredith and Tobias (1960)对同一问题给出了一个解析表达式：

$$\Lambda_{eff} = \frac{\frac{2+\Lambda}{1-\Lambda} - 2\epsilon_\beta + 0.409\frac{6+3\Lambda}{4+3\Lambda}\epsilon_\beta^{7/3} - 2.133\frac{3-3\epsilon_\beta}{4+3\epsilon_\beta}\epsilon_\beta^{10/3}}{\frac{2+\Lambda}{1-\Lambda} + \epsilon_\beta + 0.409\frac{6+3\Lambda}{4+3\Lambda}\epsilon_\beta^{7/3} - 0.906\frac{3-3\Lambda}{4+3\Lambda}\epsilon_\beta^{10/3}}$$

我们把这两个公式应用到我们的模型上（其中$\epsilon_\alpha = 0.64$），并把它们与Avizo实验仿真和张量计算所得到的结果作比较：

$\Lambda$	Maxwell $\Lambda_{eff}$	Meredith and Tobias $\Lambda_{eff}$	Avizo TCExp $\Lambda_{eff}$	Avizo TCTensor $\Lambda_{eff}$
10	2.110	2.243	2.138	2.138
100	2.611	2.884	2.732	2.732
1000	2.680	2.976	2.820	2.820
10000	2.687	2.986	2.830	2.830

我们观察到解析估计与仿真之间吻合良好。

14.4.9.3 具有颗粒-颗粒接触的3D球体周期性阵列

我们研究的第三个模型是具有颗粒-颗粒接触的3D球体周期性阵列，如图14.55所示。接触面由接触尺寸与球半径之比控制。这些球体被称为$\beta$相，而包围球体的材料被称为$\alpha$相。

图 14.55：从周期性阵列中提取出的两个处于颗粒-颗粒接触的球体。

Whitaker (1999)报告了Nozad等人(1985)在具有颗粒-颗粒接触的二维方形阵列上所计算的结果。在那里，Nozad等人把这些结果与在一个较窄孔隙率范围（$\epsilon\alpha$在$[0.39, 0.41]$内）上的实验测量作了比较。尽管两者的特性存在重要差异，理论与实验吻合良好。基于这一观察，我们把Nozad等人的计算结果与在我们模型上（其中$\epsilon\alpha = 0.51$）的Avizo仿真作比较。由于Avizo实验仿真与张量计算的结果几乎相等，因此只显示了一条曲线。

图 14.56：Avizo仿真与Nozad等人热导率计算的比较。

我们观察到Avizo曲线的行为是正确的，并且它相对于其他曲线的位置大体上是可接受的。

14.4.9.4 周期性双层复合材料

我们研究的最后一个模型是周期性双层复合材料，如图14.57所示。最薄的层材料被称为$\beta$相，另一种材料被称为$\alpha$相。我们定义$r$为各层之间的厚度比。层状结构假定在Z方向上。

图 14.57：周期性双层复合材料的示例。

Auriault (2009)针对这种双层材料给出了有效热导率张量的一个解析表达式：

$$\overrightarrow{\overrightarrow{\lambda}}_{eff} = \begin{pmatrix} r\lambda_\alpha + (1-r)\lambda_\beta & 0 & 0 \\ 0 & r\lambda_\alpha + (1-r)\lambda_\beta & 0 \\ 0 & 0 & \frac{\lambda_\alpha\lambda_\beta}{r\lambda_\alpha + (1-r)\lambda_\beta} \end{pmatrix}$$

我们针对两组厚度比与热导率比，把该解析表达式与Avizo张量计算的结果作比较：

$r$	$\Lambda$	Auriault $\overrightarrow{\overrightarrow{\lambda}}_{eff}$	Avizo TCTensor $\overrightarrow{\overrightarrow{\lambda}}_{eff}$
$\frac{1}{3}$	0.02	$\begin{pmatrix} 33.667 & 0 & 0 \ 0 & 33.667 & 0 \ 0 & 0 & 2.885 \end{pmatrix}$	$\begin{pmatrix} 33.667 & 0 & 0 \ 0 & 33.667 & 0 \ 0 & 0 & 2.885 \end{pmatrix}$
$\frac{1}{4}$	0.01	$\begin{pmatrix} 75.25 & 0 & 0 \ 0 & 75.25 & 0 \ 0 & 0 & 3.883 \end{pmatrix}$	$\begin{pmatrix} 75.25 & 0 & 0 \ 0 & 75.25 & 0 \ 0 & 0 & 3.884 \end{pmatrix}$

通过Avizo张量计算所得到的结果与解析表达式完美吻合。

参考文献：

1. Maxwell J.C., Treatise on Electricity and Magnetism, Vol. I, 2nd edition, Clarendon Press, Oxford, 1881
2. Meredith R.E., Tobias C.W., Resistance to potential flow through a cubical array of spheres, J. Applied Phys., Vol. 31, 1270-1273, 1960
3. Nozad I., Carbonell R.G., Whitaker S., Heat conduction in multiphase systems I: Theory and experiment for two-phase systems, Chem. Engng. Sci. 40, 843-855, 1985
4. Ochoa-Tapia J.A., Stroeve P., Whitaker S., Diffusive transport in two-phase media: Spatially periodic models and Maxwell’s theory for isotropic and anisotropic systems, Chem. Engng. Sci. 49, 709-726, 1994
5. Perrins W.T., McKenzie D.R., McPhedran R.C., Transport properties of regular arrays of cylinders, Proc. Roy. Soc. Lond. A369, 207-225, 1979
6. Whitaker S., The method of volume averaging, Theory and applications of transport in porous media Vol. 13, Kluwer Acad. Pub., Dordrecht, 1999
7. Auriault, J.-L., Boutin, C., Geindreau, C., Homogénéisation de phénomènes couplés en milieux hétérogènes 1, Lavoisier, 2009