微胶囊变形流动的三维数值模拟范文

《水动力学研究与进展A辑杂志》2015年第六期

摘要：

微胶囊是由一层极薄的弹性材料包裹微液滴而成的柔性细容器，在医药、食品、生物技术等领域应用广泛。该文考虑边界效应、流动黏性力和薄膜力学特性等因素，采用边界积分结合有限元方法的三维耦合数值模型，对矩形微流管中微胶囊的变形流动展开数值研究。结果表明，与正方形微流管相比，矩形微流管中微胶囊变形总长度下降，但宽度显著增加。长方形微流管中，微胶囊的变形程度与毛细数、尺寸比等参数密切相关。该结果有助于人们加深矩形微流管中微胶囊流动规律的认识，也为利用长方形微流管测量薄膜力学特性提供初步的理论依据。

关键词：

微流管；微胶囊；边界积分方法；流固耦合

充液微胶囊（liquid-filledmicrocapsule）简称微胶囊，是由一层极薄的固体弹性材料包裹微液滴形成的柔性微容器，其直径通常在1–1000μm之间[1]。自然界中类似结构随处可见，红细胞、囊泡等均是典型的天然微胶囊[2]。微胶囊结构将内容物与外界环境隔离，能在微米尺度上实现内容物运输、可控释放等功能[3]。自上世纪50年代首次工业应用以来，人造微胶囊技术迅猛发展，研究成果不断增加，相关产品遍及医药、生物技术、食品、新材料等与人们生活密不可分的领域[4,5]。

多数应用中，微胶囊功能的实现与其在微流管系统中的变形和流动行为密切相关。工业上，微流管中人造微胶囊的流动规律与其制造、运输、测量及应用等环节密不可分[6]。人体内，毛细血管中红细胞的变形流动是微循环的基础，也与镰状细胞贫血、疟疾、心血管疾病等联系紧密[7]。类似例子不胜枚举。由此可见，展开深入研究，揭示微流管中微胶囊的流动机理，具有重要的科学意义和现实背景。过去几十年里，关于微胶囊流动的数值研究迅速发展。微胶囊流动与变形的理论研究始于上世纪80年代。Barthès-Biesel[8]（1980）首次对无界剪切流中微胶囊的小变形问题给出了摄动解。Secomb&Skalak[9]（1982）假设微胶囊紧靠管壁作轴对称运动，利用润滑理论(lubricationtheory)描述微胶囊与管壁狭缝间的流动，建立起早期的理论模型。尽管上述理论模型并不完善，但其中许多建模思想和假设条件至今仍有借鉴意义。研究者也提出了大量描述微胶囊流动和变形的二维数值模型。Secomb等[10,11]（2007）提出了利用相互连接的黏弹性单元模拟微胶囊流动的二维有限元模型，随后被用于模拟二维分叉微流管中的多微胶囊流动。Kaoui等[12]（2011）将格子玻尔兹曼方法与浸入边界法相结合，对平行平板间剪切流作用下的微胶囊流动进行研究。

当微胶囊运动呈轴对称形式时，轴对称数值模型能在保证精度的同时提高计算效率。基于边界积分方法，Li等[13]（1988）提出了最早的描述无界拉伸流（elongationflow）中微胶囊变形的轴对称数值模型。值得一提的是，以轴对称模型为基础，Lefebvre等[6]（2008）提出了利用微圆管测量人造微胶囊薄膜力学特性的方法：在微圆管边界的限制下，微胶囊会产生稳定的轴对称形变；分析实验图像并与相应的轴对称数值结果进行对比，可逆分析（inverseanalysis）确定薄膜剪切弹性模量。该方法成功避免了微吸管、原子力显微镜等传统测量手段操作繁杂、设备昂贵等缺点，因而备受关注。但现实中，由于微流管形状多变，微胶囊的非轴对称运动十分普遍，因此轴对称数值模型的适用范围非常有限，三维数值模型不可或缺。

事实上，描述微流管中微胶囊流动的三维数值模型也在不断发展。Hsu&Secomb[14]（1989）最早利用润滑理论分析了微圆管中红细胞的非轴对称运动。但限于润滑理论的前提假设，该模型仅适用于微胶囊紧靠管壁的情形。基于边界元方法，Pozrikidis[15]（2005）对初始形态各异的单微胶囊在微圆管中的流动和变形展开研究。结果发现，在薄膜承受负应力时，其模型结果数值不稳定，存在改进的空间。Doddi&Bagchi[16]（2007）将界面跟踪法与浸入边界法相结合，对微胶囊在无限大平行平板间的迁移运动进行了三维数值研究。Walter等[17]（2010）提出了边界元和有限元相结合的三维耦合模型，显示出了较好的数值稳定性和计算精度，但其结果仅限于无界流场中单微胶囊的流动。本文作者等[18]（2012）曾改进边界元和有限元的三维耦合方法，对微圆管和矩形截面微流管中单微胶囊的流动和变形展开了系统研究。利用获得的数值结果，本文作者等[19]（2013）曾展开实验，使用正方形截面微流管成功测量了人造微胶囊表面薄膜的力学特性。但迄今为止，矩形微流管中微胶囊流动的三维数值研究并不多见。

本文利用边界元（流体流动）结合有限元（薄膜变形）的三维耦合数值模型，对矩形微流管中微胶囊的变形流动展开数值模拟。首先基于正方形微流管，分析微胶囊流动特性与变形特征，理论分析其机理；其次对比正方形和长方形微流管，研究微流管几何外形对微胶囊变形的影响，探讨皱折等现象的生成机理；最后着眼于长方形微流管，对不同流动参数下微胶囊流动展开数值模拟，分析毛细数、尺寸比等参数变化对微胶囊变形规律的影响。

1数值模型

1.1问题描述如图1所示，矩形微流管高度为2l，其截面分别为正方形或长方形（宽高比为2）。初始悬浮的圆球形微胶囊半径为a，内外均为同种不可压牛顿流体（密度为ρ，黏性为μ）。微胶囊表面薄膜为一层极薄的不可透材料，内外流体间不存在物质交换。受到微胶囊扰动前，微流管内平均流速为V。微胶囊进入微流管后，在流体黏性力、边界限制、薄膜力学特性等交叉作用下，常呈现较大变形。

1.2薄膜定律薄膜力学特性与微胶囊的变形程度和流动特性密切相关。数值模拟过程中，合理选择薄膜定律十分重要。作者在前期工作中曾成功测量一类卵清蛋白（ovalbumin）微胶囊薄膜的力学特性[19]，结果发现neo-Hookean定律能很好地描述该类薄膜的力学行为。本文继续采用neo-Hookean定律，假定薄膜是一层三维的各向同性不可压材料。薄膜厚度极薄，其抗弯刚度可忽略，因此形变仅在平面内发生。薄膜的剪切模量和面积扩张模量分别为SG和SK（neo-Hookean定律中假定3SSK=G）。

1.3数值过程微流管内微胶囊流动涉及大变形的流固耦合问题。本文采用边界积分（内外流动）结合有限元（薄膜变形）的三维耦合数值模型，对矩形微流管中微胶囊的变形流动展开数值模拟。该模型曾成功用于求解无界剪切流、圆形或正方形微流管环境中的微胶囊流动，其优势在于仅需在流场边界上划分网格。该数值模型的详细描述见文献[17,18]，本文仅简要介绍数值过程。从球形二十面体（20个三角单元围成的正多面体）出发，依次细分各三角单元，可获得满足所需精度的微胶囊薄膜网格[18]。本文采用的薄膜网格含1280个2P单元和2562个节点。矩形微流管网格可由开源软件Modulef（INRIARocquencourt,France）生成。该数值模型的主要输入参数为毛细数Ca与尺寸比a/l。毛细数/SCa=μVG，主要衡量流体黏性力与薄膜弹性力之比。当流体和薄膜材质不变时，Ca可视为流动强度的度量。尺寸比a/l为微胶囊初始大小与微流管截面高度的比值。当微流管高度恒定时，a/l可视为微胶囊初始大小的度量。该模型的主要输出结果为微胶囊形变、质心速度等。数值计算中，数值迭代采用无量纲时间步5tV/l510−Δ=×。当薄膜表面积达到稳定值时，我们认为微胶囊的变形流动达到稳定状态。本文假定，若无量纲时间段Vt/l=1内，薄膜表面积变化小于5×104（24πa）时，即认为数值结果稳定从而结束数值模拟。

2结果与讨论

2.1正方形微流管本文首先考察正方形微流管中neo-Hookean薄膜微胶囊的变形流动特性。考察毛细数Ca=0.05时正方形微流管中尺寸比a/l=1.1的微胶囊，数值计算获得其稳定变形如图2。稳定状态下，微胶囊形状和表面积维持不变。因为内部流体质点间无相对运动，所以内部流场压强均匀稳定。在流管边界限制下，微胶囊高度下降但长度沿流向伸展。众所周知，管流中压强沿流向下降，因此微胶囊头尾承受不同的环境压强。在内部流体压强为常数的情况下，微胶囊头尾部曲率形成显著差异，即下游头部曲率远大于上游尾部曲率。当流动强度较大时（见图2），微胶囊尾部向内凹陷形成空腔，呈现常见的“降落伞状”（parachuteshape）形变。因neo-Hookean薄膜富于弹性，该类微胶囊容易变形。在流体粘性力、边界限制等交叉作用下，微胶囊逐渐趋于流管截面外形，见图2。微胶囊能延伸至正方形微流管的边角区域，呈现出典型的非轴对称形变。由此可见，前人的轴对称数值模型并不适求解该类非轴对称问题。本文利用三维数值模型求解矩形微流管中的微胶囊流动问题具备合理性。

2.2不同截面微流管（正方形vs长方形）本文采用高度一致的正方形和长方形（宽高比为2）微流管展开数值模拟，研究流管几何外形对微胶囊变形流动的影响规律。采用固定的毛细数Ca=0.05和尺寸比a/l=1.1，即流动强度、微胶囊初始大小等均一致。数值模拟获得不同管道中微胶囊稳定变形，微胶囊在各剖面上的轮廓对比见图3。图3（a）中管流从左至右。微胶囊在流管边界挤压下沿流向伸展，不同微流管内微胶囊的形变特征差异明显。与正方形微流管相比，长方形微流管中微胶囊变形程度较小，不仅总长度更短，而且尾部空腔更浅。上述差异在图3（c）所示的俯视图中尤为明显。微流管截面的几何差异是造成微胶囊形变差异的主要原因。与正方形相比，长方形截面能为微胶囊提供更多的变形空间。长方形微流管中微胶囊能沿展向延伸，从而形成更长的变形宽度（x方向）和较矮的变形高度（y方向）。相应地，正方形微流管具有更强的边界约束效应。在边界约束下，微胶囊薄膜内可能形成负应力，造成薄膜质点间相互挤压，因此在宏观上引起薄膜皱折，见图3（b）。在逆分析实验结果测量薄膜力学特性时，需要测量微胶囊剖面面积[6,19]，薄膜皱折会引起测量困难或结果误差。利用长方形微流管减少皱折生成，可能是一种值得考虑的方法。

2.3长方形微流管本文关注长方形微流管中微胶囊的变形流动规律，着重研究毛细数Ca、尺寸比a/l对微胶囊变形规律的影响。首先采用固定的尺寸比a/l=1.1，令毛细数Ca从0.05、0.10到0.15依次增加，考察不同流动强度下同等大小微胶囊的变形规律，数值模拟获得微胶囊稳定变形的剖面轮廓如图4。随着毛细数Ca扩大，微胶囊沿流向的总长度不断增加，展向宽度不断增大。Ca越大，微胶囊尾部空腔越深，而且头部曲率增加。由此可见，毛细数Ca变化对微胶囊形变的影响十分显著。进而采用固定的毛细数Ca=0.05，使尺寸比a/l从0.9、1.0到1.1依次增加，考察相同流动强度下微胶囊初始大小对稳定变形的影响规律，数值模拟获得微胶囊稳定变形的剖面轮廓见图5。随着尺寸比a/l增大，微胶囊沿流的总变形长度逐渐增加，尾部空腔不断加深，而且微胶囊的变形宽度也逐渐扩张，见图5（a）、图5（c）。当毛细数Ca等参数固定时，尺寸比a/l增加使微胶囊与流管边界间的流动润滑层（lubricationfilm）更薄，流管边界的限制作用更加明显，从而引起微胶囊变形程度增加。值得注意的是，不同尺寸比a/l微胶囊的变形高度非常接近（见图5（b）），似乎高度方向的形变受尺寸比a/l影响较小。流管的几何特性是引起上述现象的主要原因。长方形微流管的高度小于宽度，即使对尺寸比较小(a/l=0.9)的微胶囊而言，高度方向的变形空间也十分有限，因此高度形变程度较早地达到极限。模型参数变化能造成结果显著差异，是利用数值模型对实验结果展开逆分析的前提。上述结果表明，模型输入参数（Ca与a/l）对长方形微流管中微胶囊的变形程度影响显著，因此基于长方形微流管进行逆分析以测量薄膜力学特性具备一定可行性。

3结论

本文利用边界元（内外流动）耦合有限元（薄膜变形）的三维数值模型，对矩形微流管中微胶囊的变形流动展开数值模拟。基于正方形和长方形微流管，较深入地探讨了不同参数下微胶囊的变形规律。结果表明，微流管几何外形对微胶囊变形特性影响显著。正方形微流管中，微胶囊变形程度较大，其薄膜易出现皱折。长方形微流管中，微胶囊的变形程度明显降低，薄膜皱折不易生成，毛细数Ca和尺寸比a/l等参数对微胶囊变形特性影响显著。由此可见，利用长方形微流管进行逆分析，不仅能规避薄膜皱折等难解现象，而且具备理论可行性。本文结果有助于人们深入了解矩形微流管中微胶囊的变形流动规律，也为利用长方形微流管进行逆分析以测量薄膜力学特性提供了初步的理论依据。

作者：胡徐趣 吴排青 王滔 单位：湖南大学 汽车车身先进设计制造国家重点实验室 湖南大学 机械与运载工程学院