Trans Pore的三维模型：模拟多孔介质的干燥在传热传质.doc

Trans Pore的三维模型：模拟多孔介质的干燥在传热传质Patrick Perr，Ian W. TurnerbE.N.G.R.E.F；森林实验室，试验产品队伍，Girardet 14街,F-54042南希Cedex,法国；数学科学学院，昆士兰科技大学；邮政信箱 2434 Q400布里斯班 澳大利亚东部城市摘要在运输方程支配干燥过程耦合和高度非线性的性质，强调了在这一领域的适用性数值模拟。这项工作主要是介绍关于最新版本的Trans Pore，这是个独特的三维干燥模型,能够处理多孔介质的传热传质。在目前的工作中使用的方程组，是为在宏观尺度模拟多孔介质干燥更全面的物理描述。一些模拟结果表明这样的一个工具能够提供新的可能性，尤其是该研究面临的是数效应研究了各向同性介质（一种轻质混凝土立方体）和各向异性介质（木板，其各向异性率和几何因素均较大）。这些模拟测试证明，只有三维空间都能够捕捉到了在高温干燥情况下中等宽度的效果。确实，在这种条件下，纵向方向需要突出一个强烈的各向异性介质中的超压的效果，而相应的厚度和宽度是热场之间的耦合和压力驱动流所不可缺的。1.引言在这个实验中我们将讲述一个新的三维的干燥模型。同时我们将从三个宏观状态变量（一个是能量守恒方程，一个是水分保持，另一个是气相的保持）来展开说明这个干燥过程。在这个构想中一些重要的物理现象包括毛细管和对流的流体流动，约束和气体扩散以及对流气体迁移，所需的导热性能及热蒸发。最后，等温曲线，热量和质量交换条件，在充分的饱和条件下形成这个数值模型。由此产生的代码很复杂、也很费时间，这与大型网络所遇到的三维空间网络中所消耗的CPU的时间不差上下。为了抵消这种结果，采用了最先进的数字技术，比如用来解决网络的非线性离散的准牛顿法以及用于解决稀疏矩阵系统的线性化的预条件共轭梯度。迄今为止，已有一些综合二维的模型可以用来研究干燥工艺。这些模型已被应用于木材烘干领域和建筑材料，比如砖。这些模型为干燥制度提供了一些基本的知识不仅使人们更了解这一程序同时，提出该更多具有创新的新的建议。 研发三维的干燥模型是根据以前一些专家有关二维模型的报告所得出的结论而引发的。首先，木材干燥应力在干燥过程中发展依靠强烈板的宽度，而纵向的方向是非常重要的代表传递现象。其次，在这种模式之前保证全面的无线电频率或模拟微波干燥可进行是绝对必要的。三维空间的模型在多孔介质中的使用并不陌生，但是据我们所知，综合性很强的计算模型目前在干燥研究这块领域还不存在，不过存在一些简单的用于研究木梁的吸收和解吸水的性能以及分析木料木块的结构的传输模型。高度耦合的非线性方程组的完整解决干燥模型似乎是一个新的有价值的文献。在地下水流领域，油藏模拟和污染建模可以提供有关的结构化和非结构化网格实现3 - D模型参考资源。虽然在理论上，从二维到三维是很顺其自然一气呵成的，但从以上所有研究中我们可以看出要实现模拟时间的高效计算，线性代数中一些微妙的数据分析和创新性概念是必不可少的。在这项研究中已经有两项研究来检验三维模型的性能。在实验1中描述的是在高温传递下轻质混凝土的烘干；实验2 针对高温与低温两种情况所产生的结果进行了分析。尤其需要强调的是在整个烘干过程中只有三维才能够正确地捕捉到多孔介质的物理行为。的确，高温干燥超压在一个重要的板材纵向的一个板盖上有高度各向异性的多孔介质。厚度和宽度的板必须解释热耦合机制、质量扩散和达西流。2.数学公式这个部分将着重介绍下干燥守恒方程以及相关通量和本构关系表达式。解析下离散化过程、简要概述下数值求解过程。还有讨论一些新的3 - D模型具体计算细节。2.1传输模型Whitaker详细研究了在多孔介质支配传热传质现象的宏观方程。这些公式后来也相应的被应用于软木干燥过程的模拟中。在那里关于它们有呈现几种不同干燥配置的能力也已得倒了证实。这些方程如下：液体守恒能量守恒气体守恒根据达西定律可以得出液体与气体的相对速度已知相潜力和是深度标就可以得出所求的数。其他符号还是原来的含义。2.2.边界和初始条件 有两种边界条件需要考虑，一个是外部条件，另一个是对称面。边界的外部干燥的样品表面提出的条件被认为是如下形式： Jw 和Je分别记录了在边界上总水的通量和总焓的通量。在干燥表面的外部压力是固定在大气的值。在计算时起重要作用的其中一个变量是大气的平均密度，Dirichlet边界条件Pg=P不得不修改了，以便形成一个适当的非线性方程为这个主要变量：在这期间每个外部边界在计算区内的控制卷的非线性迭代，而根据守恒定律，所以方程必须解决掉。为了减少总体所需的计算时间，在这个模型中我们采用了对称平面。把这个多孔性样本考虑在内的话，只有1/8的样品被计算过时间，同时所有流动的液体、蒸汽、空气和热气的通量在这个平面上检测到的是0。起初这些多空介质是有一定量的水分和热的分布。在这个过程中大气压力值一直都保持不变。2.3数值策略用于运输模式分辨率为了说明离散过程，同时也为了弄清楚这一过程，已决定从解析下液体守恒定律入手。这整个过程是和其他两个守恒定律是相同的。整个过程是把平衡定律与一个典型的立体的控制卷结合一起的，如图1 所示从文章中看在整合了计算机内的每一个控制卷后，一个非线性方程组的结果使本系统内每一个公式可以转化成离散的守恒方程，比如：（7）在方程7中表示的是总液体守恒量：。代表的是该点周围邻近节点p设置。这个方程式中的张量条件绝对磁导率；迁移率和扩散率在流量表达式是系统功能状态的变量。由于方程式7它是有必要近似于每一个梯度术语。对于正交网格，一定量的梯度可以和这个公式一样简单。但是需要注意的是当网格不是正角正交的时候，这个公式需要修改下，在这种情况下一个渐变的切向分量的修正项必须添加。然而，把这个公式带入7，就得到以下离散形式的光通量的结果：在这个公式中假设张量条件是对角线。如果该系统的各向异性放进笛卡尔坐标中的主惯性轴中，这个假设是合理的。这个复杂的配置已经在其他地方详细讨论过了，根据木板从树上切下了的位置可以来旋转材料轴。向量表示的是几何量，如例子：一个有效计算每个控制体积面标准的方法是运用式子（可以参考下图表1以及face ABCD)。注意下像术语表示的是at the th CV face对部件液体的流动张量评估。 离散模拟(7)液体的守恒定律是一个复杂的非线性函数的状态变量在(Sw, Pw, T, , , )节点p和周围邻居节点在控制卷考虑的范围内。表达式中的系数是把运输系数的非线性和控制体积面的相关几何因数集合在一起的。其中运输系数的非线性涉及到节点P和邻居节点nib的状态函数。注意，平流基本运输系数和扩散条件要区别对待，否则的话，方案可能不是单一的。 对于每一个变异系数来说，系统状态变量是被认为有一个确切的定义。事实上，六个变量中有三个是主要变量，而这三个主要变量对于每一个变异系数来说都是未知的自变量。剩下来的三个是辅助变量，它们是将根据主要变量的值来计算。所有变量的值在第N个时层中是假定已知的。一旦这些非线性函数在计算域为每个变异系数组合，其结果是一个形式的非线性函数方程组。对于每一个在网格内的变异系数，解向量U包含了主要的变量（三个一组）。为了及时加快所有的变量，这个系统每一次都将被计算。 这个解法是运用我们所熟知的内部与外部迭代阶段这个知识，在不同的阶段进行的。在外部循环期间，非线性方程组在拟牛顿方案中直线化了。在（n+1）段，解向量的估算是根据目前的解决方法，也就是记录： （9）和解出直线方程 （10）得出正确的向量来计算出n 阶段而得出的。在方程10中B代表的是一种逼近真实的雅可比矩阵。作者们最近正在忙着写一本以B的整合逼近为题材的书。然而，可以这样明确地说为建设所需要的繁生物是通过运用简单的有限差分近似进行数值运算的，来达到对过多而且昂贵的函数值需要与逼近接的平衡。为了说明这个方案，我们假设方程7可以写成以下形式：在这个式子中，比如，这个式子代表的是水分含量同时表示的是温度的无因次。关于 中的的派生物可以用以下方程式计算：在这个式子中和是经过挑选的适当的数值位移值。他们必须大到以避免舍入错误同时也要避免小到导数差逼近。假设那些主要的变量都是非因此为，然后有一个有用的经验法则- 根据该平方根的舍入的单位计算转移因子。我们已经利用，我们可以分析二次收敛的牛顿计划。雅可比矩阵（近似）的结构很稀疏，由于这一点，所以在解雅可比矩阵的时候要反复地进行演算。在这项研究中，在内在迭代相位，利用共轭梯度法方程式，和零水平的不完全因式分解，预处理技术可以解出方程式。另外一个对收敛迭代求解很重要而且同时对预优因子意义也非同一般的因数对于在解决方案领域中节点的排序有很大的关联。显然，如果镥的分解过程将被用来作为预优因子，那么必须尽量减少稀疏矩阵的宽带。这种小化可以通过使用反向的CuthillMcKee计算法。2.4 依据控制体积面评价张量条件从方程中可以看出，为了根据控制体积面来评价张量条件，我们有必要为空间平均技术下个定义。扩散条件往往均匀地被用于算术和调和平均的策略上。我们再次回到图表1中，同时给插值因数下定义：式中的扩展张量的组成部分和都对等的用于以下的算术平均法：该平流和对流的正确治疗对于数字计算格式的准确性和性能很重要。这项研究中上游加权已经被采用了。上游加权是一介空间平均法，这个空间平均法是已知物理的正确方法的收敛： 3.结果和讨论这一节我们将介绍下由上面讲到的计算机模型所得出的仿真结果。这项研究也主要是利用仿真结果来阐明在烘干过程中的三维耦合机制。希望通过由这个新的三维工具所提供的仿真的可能性将在整个随后的讨论中得到凸现出来。结果是在两个不同的配置中展现的，在高温烘干下从一块立体的各同向性材料（轻混泥土）形成一种强烈的各向异性材料（软质材料）的长形样本。相对于现存的TransPore两维，我们对于木块宽度的高度留意是一大提高。所以在选木头，我们是选三种不同大小的木板，而且随后我们还要对其进行相匹配： 1m3cm3cm 分别代表的是L-, R- and T-directions,); 1m6cm3cm 分别代表的是L-, R- and T-directions, 1m12cm3cm 分别代表的是L-, R- and T-directions,对于所有的配置来说，三个对称面的存在使得代码只能处理样品的八分一。此外，为了节省计算时间，网格不得不细化到接近交换表面，在这里变量梯度有可能变得最重要。不管样本的宽高比是多少，为了保证控制容量体积（xyz坐标是000）接近于一个立体，几何因素要考虑到内。三个网格的例子参考图2.在目前使用的网格是很粗糙的(12厘米厚的木板的L-, R- and T-direction 201510)，但是，他们是足够捕获所有在干燥过程所产生的物理现象。关于这个实验的证据可以在许多专家所进行的研究中找到。图表2 三个网格样品：留心下网格在交换表面细化部分： 一个444cm的混凝立方体的1/8，一个3350厘米的木板的1/8, 一个31250厘米的木板的1/8。3.1.在一个立方体的各向同性介质中的高温对流干燥TransPore三维版本的第一次试用是一个很简单的设置，一立方轻质混泥土，通过对流气流而烘干（干球温度是140C露点温度是80C）。图表3表示的是在2个小时干燥过程中计算变量的等值面。三个正交交换面的存在导致物体底部的水分含量很低。在当时的干燥期间，只有这个立方体的核心保持很高的水分含量。温度场和压力场使得更多微妙的现象得到观察。只有在水分迁徙的区域空间中的温度变化很大（水汽扩散，更为重要的是达西气体的流动。在这个领域，温度梯度是需要为液体转化成气体提供必要的热量。由于水蒸气和空气的交叉扩撒效应导致了气压梯度在这个领域的存在。根据达西法得出的总体流量运用到气相中，使得总空气量几乎为零。这个现象说明了为什么只有两个已经标记的气压（pressure isovalues）接近于交换表面。而在自由水领域几乎没有温度梯度的存在。图表3 混泥土干燥的数值模拟（干燥以及湿球温度计的温度分别是140/80C）。三维变量领域计算的是2个小时或烘干。在文章中看为了更清楚的展现三维空间转移的效果，另外两个模拟实验在同样的几何和相同的干燥条件下，但是不同的交换面下进行了。这个策略模式允许一维、二维或是三维的传动机制之间直接相互对比。平均含水量、中心温度以及中心压力在时间上的演化在图表4中得到了描述。从演化的曲线图看，在这个持续的过程中交换面数量所带来的影响可以被量化。有趣的是干燥过程所使用的时间和交换面的数量成正比，和一维方案中干燥所用时间相比，三维大约是一个半小时和两个交换面，如果是三个交换面则只需20分钟。很明显，这个简单的规则是由于每个方向的作用都是等效的所造成的：立体几何，同样的交换条件和各向同性介质。需要注意的是，除了干燥所花时间是一样外，所有的身体机制也是一样的。在干燥行为中并没有观察到根本性变化。尤其值得注意的是，超压的最大值是一样的，同时在干燥率保持在恒定不变的状态下，温度为8590C的稳定值是遵循所有的运行的。这个稳定水平是轻质混凝的一个特点，它源于一个动态的平衡，而这个动态平衡存在于热传递和蒸汽通量（扩散和体积流）中。在实验中它已经得到了解析，利用解析模型对其进行了解释说明。 3.2低温对流使软木干燥对流干燥软木是第二个选来测试代码的配置。和第一个测试相比较，现在的材料是高度各异性的，样品的形状涉及到一些重要的几何因数(1m12cm3cm in the L-, R- and T-directions)。干燥流体是湿空气，其中干球温度低于水的沸点（干球温度是80C露点温度是50C，热传系数是25 )鉴于三维结果很难使用图表表示和分析，所以我们通过画出在不同干燥时间的“篱笆图”（fence plots）说明三个状态变量之间的耦合。由于计算网格不是规律的进行排列，所以这些截面不得不通过线性插值来计算。图5 图6 图7描述的是在某些重要的干燥时期，三维领域的水分含量，温度和无量纲压力，换句话说就是；5分钟：短暂的加热期间。4小时：干燥率恒定时期到干燥率减少的转变30小时：干燥过程结束5分钟干燥后，温度场上的温度是不均匀的。木板的初始温度30C是远远低于露点温度。所以，是对流热传递和凝聚在表面的蒸汽而产生的潜伏热使得交换面开始变热。样本的核心类似于物体的初始温度场，但是样本最热的地方是靠近两个或三个毗邻交换面截距（木板的边距或是角）那里。图表4混泥土干燥的数值模拟（干燥以及湿球温度计的温度分别是140/80C）。空间的维数对于整个干燥曲线的影响。水分含量场导致了以下现象： 表面凝结； 在最热的部分开始了干燥； 纵向迁移有效地驱使缩合通量从终端件到样品的内部部分。由于温度的急剧上升，样品内的气压也随之上升。注意的是狄里克雷边界条件阻止了交换面附近的压力上升，特别由渗透率值很高所引起的纵向流变。 经过4小时的干燥过程，这个过程已经很好的确定好了。由于外部转化和内部转移的可能性决定了水分含量场。此外，三个转化面的存在以及材料各向性的效果可以从水分含量场观察得到。从一处可以辨别出两个不同的区域，木板的边缘已经在吸湿领域。结果是表面的水汽压力小于饱和蒸汽压。饱和蒸汽压可以使得温度得以升高，整个领域接近与半宽的表面还在自由水领域，木板的这部分一般是在干燥温度率恒定时期由机制设置驱动，而这时候的温度一般停止在湿球温度计的温度值，同时质量流量几乎保持不变。在压力场和温度场，相应的区域也可以看到。注意的是压力不足的地方温度还是湿球温度计的温度。这是和物理公式一样的（在恒温下气压体积增长），这个结果在之前的研究中报道过几次。有趣的是在第一个干燥期，干燥率很少是和实验曲线完全一致的。可以阐明三维计算干燥的第一阶段和第二阶段的顺利过渡期显示了更为准确的干燥曲线。在进行实验与理论的比较时，这也许是十分重要的。 最后，也就是在干燥过程进行到30小时的时候，这个过程已基本接近结束了。在样本的核心我们需要留心的是一些剩余的水分和一点点超压。在干球温度，温度场的温度接近于均强场。图表5 低温软质材料干燥(Tdryy=80C,Twet=50C) 。干燥过程5分钟中，水分含量，温度和无量纲压力的“篱笆图”。图表6 低温软木干燥(Tday=80C, Twet=50C) 干燥过程4小时中，水分含量，温度和无量纲压力的“篱笆图”。图表7 低温软木干燥(Tdry=80C, Twet=50C). 干燥过程30小时中，水分含量，温度和无量纲压力的“篱笆图”图8 低温软木干燥(Tdry=80C, Twet=50C). 整体干燥曲线和在不同设置中交换面的内压力。为了使一两个面更好地被绝缘，二维或是一维的设置已经被测算过了（无限宽度和无限长度）。同时最后的结果在全球含水量和中心压力的时间演化曲线中展现出来了。一些主要的备注可以从以下的曲线分析中获得 中心压力的演化只受出口端的状态影响（不管他是否被绝缘了还是没有）。这是在纵向流变中渗透率高的一个重要结论 全球干燥率曲线主要是受侧交换面的影响。出端口的转变可能性对于干燥曲线的影响有限的。 虽然没有二维的设置可以解释全三维几何，但在半身部分，纵向转换所带来的影响对于水分场的影响不大，因次在这个部分的应力场的作用也不大。总而言之，运用二维模型去计算半身部分的方法去测算低温配置还是很合理的。如果某个人在干燥过程对机械降解感兴趣的话，这是很正确的。3.3.软木的高温对流干燥使用和上个实验相同的木板，已经计算好高温配置的温度（干球温度140C露点温度80C,热传导系数25 ）12厘米宽的这块木板的仿真结果在干燥过程中的第2个小时，6个小时12个小时分别在图9图10图11中描述了。在第2个小时，几乎整个木板的温度都在干燥率恒定阶段。水分含量也渐渐地从内到外减少。尽管可以观察得到，但是纵向水分梯度一直保持相当低，这导致的直接后果是渗透性很强各向异性因素(103)。而且，木板的温度停滞在湿球温度计的温度，由于水分含量的降低而产生的气态阶段延伸使内压力减少到大气压。在之前的研究中，所有的机制都已经广泛的被分析过（见10例子）。在这里最有趣的是边缘效应，边缘效应是三维空间的直接人工制品。由于两个正交交换面的存在，边缘的每个区域的水分含量减少的很快。结果，木板的这部分已经到了干燥第2 阶段，这个阶段局部在升温。这一现象对压力的影响也可以在接近与角的整个木板的长度观察的到。湿区域的形成是轻微的超压所带来的结果，这个湿区域可以在交换面的出口端和接近边缘的地方看到。随着干燥程序继续，在进行到第2 阶段的时候这个区域朝半宽度（half-width）延伸。结果，一个超压的相应的区域也出现在木板的末端，同时在纵向方向产生了一个很重要的水迁移制。在木板的末端也就是与纵向投影区相应的地方中，这个现象导致了完全饱和区扩大，而这个纵向投影区有水和过多的压力存在。当中心超压充足的时候，该耦合迁移机制使得这个现象在中厚度显示出和后面的那个交换面很接近，并且逐步向mid-width延伸。图10对于一些错综复杂的机制解释的很明了。由于刚开始在拐角处的温度很高，所以在附近木板的边界出现了一个大的超压。注意的是末端出口的湿润出现在了中心部分。图9 高温软木干燥(Tdry=140C, Twet=80C). 在干燥进行的第2小时Fence plot的含水量温度和无量纲压力。另外一个有趣的现象发生在木板中部的。由于边缘部位温度和压力的上升，因此在区域内很长一段时间内出现了一个很重要的压力梯度，使得水分含量从边缘处传到了pressure gradient。在50厘米处，已经清晰地显示出相应的湿区。（水分场图10）由于含水量很高，在表面有自由水存在，使得温度接近于湿球温度。传动机制内的耦合是这样使得低温增强了压力梯度，反过来压力梯度有使得这一状态有一直持续下去。这个特殊的耦合是和模型内的物理现象和物理特性相一致的。事实上，这个现象并不是依靠于网格的大小，同时在RT区已经运用二维计算过。不过这个机制似乎对于物理参数值和干燥条件很敏感。据我们所了解的，这还从来没有在那个文献中报道过。要不然，值得留心下这个事实：末端口的再次饱和和它所给温度带来的影响已经详细的被研究过了，它的直接结果是 首先被TransPore二维版本预测的。因此，我们的观点是数值观察还需要在今后的试验中核实。图10高温软木干燥(Tdry=140C, Twet=80C). 在干燥进行的第4小时Fence plot的含水量温度和无量纲压力。图11高温软木干燥(Tdry=140C, Twet=80C). 在干燥进行的第8小时Fence plot的含水量温度和无量纲压力。随着程序的进展，这个液体区最后消失，整个末端口进入了一个吸湿性范围（图11图8干燥的第小时）。此后，这个耦合机制没有显现出来，变量梯度几乎和坐标轴的符号是一致的。温度从交换面向中心递减，无论方向如何。水分含量和压力场都是从中心向交换面递减。某些之前作用的迹象可以在温度和水分含量场中检测得出来。现在这个程序即将结束。在木板内还有能量在驱动，因为热传导。物体内液体在蒸发，湿气开始行动由于扩散（水和气和）和空气之水平流动（运用气相到达西法）。二维和一维设置也已经计算过（无线宽和无限长）。其结果在以水分含量和内压力的时间演化形式中报道过,图12:图12 软木高温干燥(Tdry=140C, Twet=80C). 在不同设置的交换面中全部的干燥曲线和内压力。从压力曲线中可以得出在mid-length and mid-width上，表面再次饱和的现象只有当交换是有可能从侧面开始时才存在。这是压力梯度驱使液体移向物体中心的一个条件。当纵向运移时，这个现象被增强了，并持续了200分钟（完整的三维模型）。一份所有的曲线图详细分析表证实了任何一个表面状态的改变都会引起全部程序运行的变动：1、内压力的演化经常受到交换面上任何一个改变影响.2、尽管整个干燥曲线主要受侧交换面的影响，但是通过达西法，我们可以得出在纵方向上的传递性对于高温设置有很大的作用。这个结果证实了在干燥时内部蒸发中我们所熟知的木板的运行状态。鉴于高温干燥这个例子，没有二维设置可以解释说明完整的立体几何。在谈及干燥率时，少不了侧面，然而末端口对于压力场的作用很大，此外也对贯穿于整个样本的机制有很大作用。总而言之，使用任何一个二维配置都会给干燥过程带来很大的粗的假设。从之前对于传递机制的分析来看，干燥的三维方面已经得到了很好的展示。特别是对于一个各向异性材料比如木块来说是尤其明显的。我们最后一个图表很好的阐明了在整个干燥过程中板料宽度的影响（图13）。在L- and T-directions上无线宽是和二维模拟相一致的