NUMERICAL STUDY OF THE MICROORGANISMS SWIMMING IN ANISOTROPIC VISCOELASTIC FLUIDS IN CHANNEL
-
摘要: 近年来, 微生物游动特性的研究逐渐成为流体力学的热点问题之一. 深入了解微生物在复杂流体中的游动特性可以为微型器件的设计制造提供帮助. 已有研究结果表明, 管道中微生物存在复杂的行为特征, 其内在的动力学机理尚未清楚. 本文采用浸没边界法, 引入泰勒波动板模型, 探讨槽道中微生物在各向异性黏弹性流体中的游动特性; 重点分析槽道壁面约束对微生物行为特征的影响规律. 数值模拟结果表明, 与自由游动相比, 槽道中微生物游动呈现加速现象, 而且槽道壁面附近加速效果更显著. 另外, 槽道中微生物游动效率比自由游动时更高. 进一步受力分析, 发现流体压力对微生物游动起推动作用, 黏性力起阻碍作用. 与这两者相比, 聚合物产生的弹性力则小很多, 几乎可以忽略. 槽道中微生物受到更大的流体压力, 认为加速现象与此相关.Abstract: Recently, microorganisms swimming in channel has gradually become one of the hot issues. The deep understanding of the microorganism swimming in complex fluids can provide guidance for the design and manufacture of micro-devices. There are complex behaviors when microorganisms swimming in channel and the hydrodynamic mechanism has not yet been clarified. In our work, immersed boundary method is adopted, and the Taylor’s swimming sheet model is introduced to investigate the swimming characteristics of micro-swimmer in anisotropic viscoelastic fluid in the channel. We focus on analyzing the effects of the wall on the micro-swimmer behaviors. Our results show that the micro-swimmer speeds up in channel, and it is significant when micro-swimmer close to the channel wall. In addition, the swimming efficiency in channel is greater than that of free swimming. We obverse that pressure plays a role in promoting the swimming of microorganisms and viscous force acts as a hindrance. Elastic force produced by the polymer is much smaller. The micro-swimmer gains greater pressure when swimming in channel, and we believe it is related to the acceleration phenomenon of microswimmer.
-
Keywords:
- microorganisms /
- viscoelastic fluid /
- anisotropy /
- channel /
- immersed boundary method
-
引 言
近年来, 在微纳米领域, 模仿微生物的微小元器件[1-3]越来越受到人们的关注, 它们在众多行业中具有极大的应用潜力. 掌握复杂流体中微生物的行为规律以及游动特性, 可以服务于微型元器件的生产制造. Taylor[4]在1951年对一类依靠鞭毛周期性摆动驱动自身推进的微生物进行了研究. 对于该类微生物, 鞭毛推动微生物在流体中波动运动 [5]. 依据微生物的运动形态特征, 以及摆动振幅$ A $远小于自身尺寸$ L $的假设, Taylor提出了经典的泰勒波动板模型; 并且采用摄动法, 基于小振幅进行渐进展开, 得到牛顿流体中微生物自由游动的推进速度. 结果表明, 牛顿流体中微生物的推进速度$ {U_N} $与自身振幅$ A $的平方成正比. 考虑到微生物生存的真实环境存在壁面约束, 如各种类型管道, Reynolds[6]基于Taylor的工作, 研究槽道中微生物的行为特征, 给出了微生物推进速度的近似公式. 其结果表明, 相比自由游动, 微生物在槽道中呈现加速现象; 且离壁面越近, 加速越显著. 其加速效果与微生物和壁面之间的距离$ {h_u} $相关. 该结论与实验研究工作[7-9]和数值模拟[10-12]的结果相符合. Reynolds的研究工作, 基于微生物和壁面间距离远大于微生物摆动振幅的假设, 并未考虑两者同一量级的情况. 有鉴于此, Katz[13]采用润滑力理论分析该问题, 并得到了狭窄槽道中微生物推进速度的理论公式. 其结果表明, 此时微生物的推进速度与$ {h_u}/A $的平方成反比例关系, 且存在最大值.
然而, 许多生物流体诸如血液、呼吸道黏液等, 会表现出明显的黏弹性流体特征[14-18]; 其中部分流体, 例如DNA分子溶液[19], 由于DNA分子特定的长链结构, 甚至表现出各向异性特征. 例如在生理温度下的各种生物膜通常表现为液晶态[20], 由双层脂质分子层和膜蛋白组成的红细胞膜在流动状态下也为液晶态. Gagnon等[21]研究在剪切变稀流体时发现, 微生物周围的流场发生明显的改变; 在其头部和尾部, 流场速度会有显著的增加. Lauga[22]对Oldroyd-B流体中微生物自由游动的行为特征进行了分析, 讨论了流体的黏弹性对微生物运动的影响规律. 在黏弹性流体中, 流场结构与牛顿流体差异较大, 微生物的推进速度与黏弹性流体的弛豫时间有关. 流体弹性应力起阻碍作用, 使得微生物呈现减速现象. 在此基础上, Ives等[23]探索槽道中微生物在黏弹性流体中的游动特性. 其结果表明, 与自由游动相比, 微生物在槽道中同样呈现加速现象, 与牛顿流体的结论相近.
进一步, Krieger等[24-25]采用Doi[26]的液晶模型研究各向异性黏弹性流体中微生物的行为特征. 对于自由游动问题, 在胆甾相液晶中, 微生物的游动速度总是低于牛顿流体中的游动速度. 在向列型液晶中, 微生物的游动速度在特定工况下高于牛顿流体中的游动速度, 并且在较高的旋转黏度下, 游动方向会发生逆转. 对于槽道中游动问题, Krieger等[24-25]首先对微生物小振幅游动问题进行理论分析, 相比自由游动的情形, 微生物推进速度变大. 其次在对微生物大振幅游动问题采用数值模拟进行研究中, 重点分析流体特性对微生物游动的影响, 与自由游动相似, 发现在较高的旋转黏度下, 微生物靠近壁面的过程中, 其推进速度的方向会发生逆转, 出现回游现象. Lin等[27-28]采用Q-tensor液晶模型[29-31], 结合理论分析和数值模拟, 对微生物在各向异性复杂流体中的游动问题进行了进一步的探索. 其结果表明, 与牛顿流体相比, 微生物游动方向与液晶取向平行时, 微生物呈现加速现象; 与液晶取向垂直时, 微生物呈现减速现象.
综上所述, 已有的相关研究工作大部分聚焦各向同性流体, 只有少部分的工作涉及各向异性流体. 因此, 本文采用数值模拟方法, 对槽道中波动形微生物在各向异性黏弹性流体中的游动问题进行初步探讨, 并尝试揭示其内在的动力学机制.
1. 模型和方法
本文研究内容属于典型的流固耦合问题, 采用浸没边界法(immersed boundary method)[32]进行求解, 将问题分解为流体子问题、黏弹性流体本构方程子问题以及微生物运动子问题进行求解. 模型方面, 本文针对波动形推进微生物, 采用经典的泰勒波动板模型, 而各向异性黏弹性流体则采用Doi[26]的Q-tensor模型. 黏弹性流体本构方程的数值求解方面, 则采用4阶龙格库塔格式结合乔列斯基分解进行求解. 下面详细介绍.
1.1 微生物模型
本文采用经典的泰勒波动板作为微生物模型. 其运动方程如下
$$ y = A\sin (ks - \omega t) $$ (1) 其中, $ A $表示行波振幅, $ k $表示行波波数, $ \omega $表示角频率. 行波初始形态$ {\boldsymbol{X}}(s,t = 0) = (s,A\sin (ks)) $, $ s \in [0,L] $.
1.2 黏弹性流体模型
本文关注各向异性黏弹性流体, 采用Doi[26]的Q-tensor液晶模型. 该模型采用棒状分子模拟流体黏弹性和各向异性的特征, 将高分子溶质简化为前后对称的线性棒状颗粒, 利用其统计分布形成的概率函数$ \varPsi ({\boldsymbol{x}},{\boldsymbol{p}},t) $来描述向列型液晶场. 描述液晶取向的二阶构型张量定义为
$$ {\boldsymbol{D}}({\boldsymbol{x}},t) = \int_{\boldsymbol{p}} {{\boldsymbol{pp}}\varPsi {\mathrm{d}}{\boldsymbol{p}}} $$ (2) 四阶张量为
$$ {\boldsymbol{S}}({\boldsymbol{x}},t) = \int_{\boldsymbol{p}} {{\boldsymbol{pppp}}\varPsi {\mathrm{d}}{\boldsymbol{p}}} $$ (3) $$ \begin{split} & {{\boldsymbol{\tau}} _p} = 2\nu {k_B}T\left({\boldsymbol{D}} - \frac{{\boldsymbol{I}}}{2}\right) - 2\nu {k_B}T\zeta ({\boldsymbol{D}} \cdot {\boldsymbol{D}} - {\boldsymbol{D}}:{\boldsymbol{S}}) +\\ &\qquad \frac{{\nu {k_B}T{\beta _0}}}{{2{d_r}{{(\nu {l^3})}^2}}}{\boldsymbol{E}}:{\boldsymbol{S}} \end{split} $$ (4) $$ \begin{split} & \frac{{\partial {\boldsymbol{D}}}}{{\partial t}} + {\boldsymbol{u}} \cdot \nabla {\boldsymbol{D}} - ({\boldsymbol{D}} \cdot \nabla {\boldsymbol{u}} + \nabla {{\boldsymbol{u}}^{\mathrm{T}}} \cdot {\boldsymbol{D}}) + 2{\boldsymbol{E}}:{\boldsymbol{S}} = \\ &\qquad -4{d_r}\left({\boldsymbol{D}} - \frac{{\boldsymbol{I}}}{2}\right) + 4\zeta {d_r}({\boldsymbol{D}} \cdot {\boldsymbol{D}} - {\boldsymbol{D}}:{\boldsymbol{S}}) + {d_t}\Delta {\boldsymbol{D}} \end{split} $$ (5) 其中, $ \nu $为有效体积分数, $ {k_B} $为玻尔兹曼常数, $ T $为绝对温度, $ {\beta _0} $是表征分子拥挤程度的经验系数, $ l $为棒状分子长度. $ \zeta $表示液晶分子之间相互作用的强度系数, 当$ 0 \leqslant \zeta < 4 $时, 液晶呈现各向同性性质, 当$ \zeta \geqslant 4 $时, 液晶呈现各向异性性质. $ {d_r} $表示旋转扩散系数, $ {d_t} $表示平动扩散系数. $ {\boldsymbol{E}} = (\nabla {\boldsymbol{u}} + \nabla {{\boldsymbol{u}}^{\mathrm{T}}})/2 $为应变率张量.
1.3 数值方法
选取行波波长$ {\lambda }_{c} = 2\text{π} /k $作为特征长度, 行波周期$ {T}_{c} = 2\text{π} /\omega $作为特征时间, 流体密度$ {\rho }_{f} $作为特征密度, 流体溶剂黏度$ {\mu }_{f} $为特征黏度. 无量纲方程组如下
$$ \nabla \cdot {\boldsymbol{u}} = 0 $$ (6) $$\left. \begin{aligned} & \frac{{\partial {\boldsymbol{u}}}}{{\partial t}} + {\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}} = \frac{1}{{Re}}{\nabla ^2}{\boldsymbol{u}} + \frac{{Er}}{{Re}}\nabla \cdot {{\boldsymbol{\tau}} _p} - \nabla p + {{\boldsymbol{f}}_e} \\ & {{\boldsymbol{\tau}} _p} = \left({\boldsymbol{D}} - \frac{{\boldsymbol{I}}}{2}\right) - \zeta ({\boldsymbol{D}} \cdot {\boldsymbol{D}} - {\boldsymbol{D}}:{\boldsymbol{S}}) + \beta {\boldsymbol{E}}:{\boldsymbol{S}}\end{aligned} \right\}$$ (7) $$ \begin{split} & \frac{{\partial {\boldsymbol{D}}}}{{\partial t}} + {\boldsymbol{u}} \cdot \nabla {\boldsymbol{D}} - ({\boldsymbol{D}} \cdot \nabla {\boldsymbol{u}} + \nabla {{\boldsymbol{u}}^{\mathrm{T}}} \cdot {\boldsymbol{D}}) + 2{\boldsymbol{E}}:{\boldsymbol{S}} = \\ &\qquad - \frac{1}{{Pe}}\left({\boldsymbol{D}} - \frac{{\boldsymbol{I}}}{2}\right) + \frac{\zeta }{{Pe}}({\boldsymbol{D}} \cdot {\boldsymbol{D}} - {\boldsymbol{D}}:{\boldsymbol{S}}) + \frac{1}{{P{e_t}}}\Delta {\boldsymbol{D}} \end{split} $$ (8) 其中, 雷诺数$ Re = 2\text{π} \omega {\rho _f}/({k^2}{\mu _f}) $, 佩克莱数$ Pe = \omega / (8\text{π} {d_r}) $, $ P{e_t} = 2\text{π} \omega /({d_t}{k^2}) $, 埃里克森数$ Er = 4\text{π} v{k_B}T/ ({\mu _f}\omega ) $表征溶质聚合物应力与溶剂黏性应力的比值. 无量纲经验系数$ \beta = {\beta _0}\omega /[8\text{π} {d_r}{(v{l^3})^2}] $.
式(7)中拟体力$ {{\boldsymbol{f}}_e} $采用Peskin[32]给出的方法进行计算. 对于泰勒波动板, 给定其运动方程的目标曲率$ {\kappa _0} $
$$ {\kappa _0}(s,t) = - {k^2}A\sin (ks - \omega t) $$ (9) 采用上述方式驱动波动板, 定义弹性能如下
$$ \begin{split} &E[{\boldsymbol{X}}(s,t)] = \frac{{k}_{s}}{2}{\displaystyle {\int }_{{\varOmega }_{L}}{\left(\left|\frac{\partial {\boldsymbol{X}}}{\partial s}\right|-1\right)}^{2}{\mathrm{d}}s + }\\ &\qquad \frac{{k}_{b}}{2}{\displaystyle {\int }_{{\varOmega }_{L}}{\left(\frac{{\partial }^{2}{\boldsymbol{X}}}{\partial {s}^{2}}\cdot {\boldsymbol{n}}-{\kappa }_{0}\right)}^{2}{\mathrm{d}}s}\end{split} $$ (10) 弹性能由两部分组成, 右端第1项表示由于拉伸产生的弹性能, $ {k_s} $表示拉伸模量; 第2项表示由于弯曲所产生的弹性能, $ {k_b} $表示弯曲模量. 拉格朗日坐标系下拟体力$ {{\boldsymbol{F}}_e} $由弹性能量的变分给出
$$ {{\boldsymbol{F}}_e}({\boldsymbol{X}},t) = - \frac{{\delta E[{\boldsymbol{X}}(s,t)]}}{{\delta {\boldsymbol{X}}}} $$ (11) 从而得到欧拉坐标系下的拟体力
$$ {{\boldsymbol{f}}_e}({\boldsymbol{x}},t) = \int_{{\varOmega _L}} {{{\boldsymbol{F}}_e}(s,t)\delta ({\boldsymbol{x}} - {\boldsymbol{X}}(s,t)){\mathrm{d}}s} $$ (12) 本文采用双线性插值函数作为离散$ \delta $函数实现物理量在欧拉网格节点$ x $和拉格朗日网格节点$ {\boldsymbol{X}} $之间的转换
$$ \delta ({\boldsymbol{x}} - {\boldsymbol{X}}) = \frac{1}{{{h^2}}}\varphi \left(\frac{{x - X}}{h}\right)\varphi \left(\frac{{y - Y}}{h}\right) $$ (13) 式中$ h $为网格步长, 函数$ \varphi $为线性函数
$$ \varphi (r) = \left\{ \begin{gathered} 1 - |r|,\quad |r| \leqslant 1 \\ {\text{0 }},\quad |r| > 1 \\ \end{gathered} \right.{\text{ }} $$ (14) 本文采用时间步分裂格式, 将方程式(6) ~ 式(8)解耦为流体运动以及黏弹性流体本构方程等两个子问题.
(1)流体运动子问题, 求解$ {{\boldsymbol{u}}^{n + 1}} $和$ {p^{n + 1}} $
$$ \begin{split} & \frac{{{{\boldsymbol{u}}^{n + 1}} - {{\boldsymbol{u}}^n}}}{{\Delta t}} - \frac{1}{{2Re}}{\nabla ^2}{{\boldsymbol{u}}^{n + 1}} = \frac{1}{{2Re}}{\nabla ^2}{{\boldsymbol{u}}^n}+ \\ &\qquad \frac{{Er}}{{Re}}\nabla \cdot {\boldsymbol{\tau}} _p^n - \nabla {p^{n + 1}} - \frac{1}{2}[3{({\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}})^n}- \\ &\qquad {({\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}})^{n - 1}}] + {\boldsymbol{f}}_e^{n + 1} \end{split} $$ (15) $$ \nabla \cdot {{\boldsymbol{u}}^{n + 1}} = 0 $$ (16) 采用投影格式将速度和压力解耦求解上述方程组, 得到速度亥姆霍兹方程
$$ \begin{split} & \frac{{{{\boldsymbol{u}}^\# } - {{\boldsymbol{u}}^n}}}{{\Delta t}} - \frac{1}{{2Re}}{\nabla ^2}{{\boldsymbol{u}}^\# } = \frac{1}{{2Re}}{\nabla ^2}{{\boldsymbol{u}}^n} + \\ &\qquad \frac{{Er}}{{Re}}\nabla \cdot {\boldsymbol{\tau}} _p^n - \nabla {p^n} - \frac{1}{2}[3{({\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}})^n}- \\ &\qquad {({\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}})^{n - 1}}] + {\boldsymbol{f}}_e^{n + 1} \end{split} $$ (17) 以及压力泊松方程
$$ {\nabla ^2}\phi = \frac{{\nabla \cdot {{\boldsymbol{u}}^\# }}}{{\Delta t}} $$ (18) 最后由修正方程求解速度以及压力
$$\qquad\qquad {{\boldsymbol{u}}^{n + 1}} = {{\boldsymbol{u}}^\# } - \Delta t\nabla \phi $$ (19) $$\qquad\qquad {p^{n + 1}} = {p^n} + \phi $$ (20) 流场采用半交错均匀网格, 空间导数离散采用二阶中心差分格式. 方程式(17) ~ 式(18), 本文采用交替方向隐式差分方法(ADI)[33] 进行求解, 更多细节可参考文献[34].
(2)黏弹性流体本构方程子问题, 求解$ {{\boldsymbol{D}}^{n + 1}} $.
为了保证构型张量$ {\boldsymbol{D}} $的正定对称性, 对张量进行乔列斯基分解[35]
$$ {\boldsymbol{D}} = {\boldsymbol{L}}\cdot{{\boldsymbol{L}}}^{{\mathrm{T}}} = \left(\begin{array}{cc}{l}_{11}& 0\\ {l}_{21}& {l}_{22}\end{array}\right)\left(\begin{array}{cc}{l}_{11}& {l}_{21}\\ 0& {l}_{22}\end{array}\right) $$ (21) 只要矩阵$ {\boldsymbol{L}} $对角元素为正数, 即可保证构型张量$ {\boldsymbol{D}} $的对称正定性质. 因此, 我们将矩阵$ {\boldsymbol{L}} $的分量进行变量替换
$$ \left.\begin{split} & {l_{11}} = \frac{1}{\text{π} }\arctan \eta + \frac{1}{2},\quad {l_{22}} = \frac{1}{\text{π} }\arctan \gamma + \frac{1}{2} \\ & {l_{12}} = \frac{2}{\text{π} }\arctan \chi \end{split} \right\}$$ (22) 将式(21) ~ 式(22)代入方程(8), 得到新的方程组. 利用已更新的速度$ {{\boldsymbol{u}}^{n + 1}} $, 采用显示格式求解新的方程组, 时间推进采用四阶龙格库塔方法, 详细细节可参考文献[27].
最后, 求解泰勒波动板的运动方程
$$ \frac{{\partial {\boldsymbol{X}}}}{{\partial t}} = {\boldsymbol{U}}(s,t) = \int_{{\varOmega _E}} {{\boldsymbol{u}}({\boldsymbol{x}},t)\delta ({\boldsymbol{X}}(s,t) - {\boldsymbol{x}}){\mathrm{d}}{\boldsymbol{x}}} $$ (23) 2. 算法验证
本节验证文中所采用算法及计算程序的精度以及可靠性. 文献[22]中Lauga采用摄动法, 得到Oldroyd-B流体中波动板小振幅情况下的推进速度
$$ {U_{{\mathrm{OB}}}} = \frac{1}{2}{A^2}k\omega \frac{{1 + D{e^2}{\eta _s}/\eta }}{{1 + D{e^2}}} $$ (24) 我们以此作为基准进行比较. 由于式(24)是基于斯托克斯流体得到的结果, 因此选取雷诺数$ Re = 0.01 $, 尽可能贴近文献参数. 流体参数$ De = 1 $, 黏度比值$ {\eta _s}/\eta = 2/3 $. 波动板参数$ L = 1 $, $ A = 0.001\sim 0.05 $, $ k = 2\text{π} $, $ \omega = 2\text{π} $, $ {k_b} = 80 $, $ {k_s} = 1000{k_b} $. 波动板游动方向选取为$ x $轴方向. 计算区域设置为$ [0,1] \;\times [0,1] $, 两个方向均为周期性边界条件. 网格步长为$ h = 1/128 $, 时间步长为$ \Delta t = 5 \times {10^{ - 6}} $.
得到波动板在稳定状态下推进速度如图1(a)所示. 当振幅$ A \leqslant 0.01 $时, 数值模拟结果与理论公式(24)基本吻合, 但仍存在大约30%左右的相对误差. 我们认为误差主要与波动板运动方程的处理方式有关. 文献[22]中Lauga推导理论公式时将波动板运动方程作为边界条件处理. 而数值模拟中, 我们采用Peskin[32]的处理方法, 引入能量函数, 结合浸没边界方法将其变分作为源项加入流体动量方程, 如上文式(9) ~ 式(11). 本质上能量函数类似于罚项, 当波动板实际运动与给定运动间有差异时, 施加一个“惩罚”, 使得实际运动趋向于给定运动. 该方法对运动方程进行近似处理, 我们认为这是误差的主要来源. 进一步, 将本文结果与文献[36]的结果进行比较. 文献[36]也采用Peskin[32]的处理方法, 但具体数值格式与本文不同. 如图1所示, 本文结果与文献[36]结果在小振幅工况下十分接近. 这说明本文的数值方法是可靠的. 有鉴于此, 本文数值模拟振幅参数设定为$ A = 0.01 $.
图1(b)给出网格收敛性验证的结果. 选取本文研究参数范围内的典型工况: $ Re = 0.01 $, $ Er = 1 $, $ Pe = 1 $, $ P{e_t} = 50 $, $ \zeta = 8 $, $ \beta = 0.000\;5 $. 泰勒波动板振幅选取为$ A = 0.01 $, 其余参数如上一部分所示. 计算区域为$ [0,1] \times [0,1] $, $ x $轴方向为周期性边界条件, $ y $轴方向为固壁边界条件, 如图2所示. 波动板中轴初始位置选取在槽道中心, 即$ {h_u} = 0.5 .$ 液晶的初始取向与波动板游动方向的夹角设置为$ \theta = \text{π} /4 .$ 图1(b)给出两种不同网格参数下波动板质心的瞬时速度. 由图中可以看出, 两种不同网格下, 波动板的瞬时速度$ x $轴方向分量$ {U_x} $差异极小, $ {U_y} $略有差异. 因此以下研究我们选取网格步长为$ h = 1/128 $, 时间步长为$ \Delta t = 1 \times {10^{ - 5}} $.
3. 结果与分析
本文主要关注槽道中微生物在各向异性黏弹性流体中的行为特征, 重点分析不同物性参数对微生物游动速度和效率的影响规律, 并进一步揭示其内在的水动力机理. 典型参数选取如下: $ Re = 0.01 $, $ \zeta = 8 $, $ Er = 0.1,1; $ $ P{e_t} = 50 ,$ $ Pe = 1,50; $ 相应 $ \beta = 0.000\;5, 0.005 ;$ $ \theta = 0\sim \text{π} /2 .$ 波动板参数 $ L = 1 ,$ $ A = 0.01, $ $ k = 2\text{π}, $ $ \omega = 2\text{π}, $ $ {k_b} = 80, $ $ {k_s} = 1000{k_b} ,$ 波动板位置$ {h_u} = 0.1\sim 0.5 $, 波动板和流体间为速度无滑移边界条件. 流场计算区域$ [0,1] \times [0,1] $, $ x $轴方向为周期性边界条件, $ y $轴方向为固壁边界条件. 网格步长为$ h = 1/128 $, 时间步长为$ \Delta t = 1 \times {10^{ - 5}} $. 重点考察微生物初始位置以及液晶初始取向对微生物游动特性的影响规律.
为了方便与自由游动工况相比较, 引入速度比$ {U_w}/{U_f} $[36]表征游动速度大小, 其中$ {U_w} $表示微生物在槽道中的推进速度, $ {U_f} $表示相同工况下微生物自由游动的推进速度. 如图3所示, 速度比值$ {U_w}/{U_f} $均大于1, 说明与自由游动相比, 槽道中微生物游动呈现加速现象. 而且加速效果主要与微生物和壁面间的距离相关, 壁面附近加速效果更明显, 微生物在槽道附近受到的约束增强. 这与已有文献[6, 13, 22]的结论相符合. 另外, 槽道中心附近速度比值几乎没有变化, 表明此时液晶取向对微生物游动速度的影响较小. 槽道壁面附近速度比值有所不同. $ Er = 0.1 $工况下, 速度比值几乎不受液晶取向的影响. $ Er = 1 $工况则不一样, $ Pe = 1 $时, 如图3(c)所示, 液晶取向与微生物游动方向接近垂直时, 即$ \theta = \text{π} /2 $, 速度比值出现极大值, 可以达到4倍以上. 对于$ Pe = 10 $工况, 如图3(d)所示, 加速效果在$ \theta = \text{π} /3 $附近最为显著, 速度比值达到5倍以上. 总体而言, 微生物与槽道壁面间的距离是主要影响因素, 而液晶取向则影响较小.
![]()
图 3 速度比值$ {U_w}/{U_f} $云图, 横轴为微生物与槽道上壁面间的距离, 纵轴为液晶取向; 图中虚线为等值线Figure 3. The contour of the speed ratio $ {U}_{w}/{U}_{f} $, the horizontal axis represents the distance between microswimmer and upper wall of the channel, and the vertical axis represents the liquid crystal orientation下面对微生物在游动过程中所受到的水动力进行分析, 探索其内在的动力学机理. 本文主要关注微生物在游动方向上, 即$ x $轴方向上的受力情况. 考虑到微生物做周期性运动, 我们考察微生物在一个周期内的平均受力情况, 即
$$\left. \begin{aligned} & {{\boldsymbol{F}}_p}(s,t) = \int_{{\varOmega _E}} {Re \cdot ( - \nabla p)\delta ({\boldsymbol{X}}(s,t) - {\boldsymbol{x}}){\mathrm{d}}{\boldsymbol{x}}} \\ & {F_{p,x}}(t) = \frac{1}{T}\int_t^{t + T} {\frac{1}{L}\int_{{\varOmega _L}} {{{\boldsymbol{F}}_p}(s,t') \cdot {\boldsymbol{e}}{}_x{\mathrm{d}}s} {\text{ }}{\mathrm{d}}t'} \end{aligned}\right\} $$ (25) $$ \left.\begin{aligned} & {{\boldsymbol{F}}_v}(s,t) = \int_{{\varOmega _E}} {({\nabla ^2}{\boldsymbol{u}})\delta ({\boldsymbol{X}}(s,t) - {\boldsymbol{x}}){\mathrm{d}}{\boldsymbol{x}}} \\ & {F_{v,x}}(t) = \frac{1}{T}\int_t^{t + T} {\frac{1}{L}\int_{{\varOmega _L}} {{{\boldsymbol{F}}_v}(s,t') \cdot {\boldsymbol{e}}{}_x{\mathrm{d}}s} {\text{ }}{\mathrm{d}}t'} \end{aligned}\right\}$$ (26) $$ \left.\begin{aligned} & {{\boldsymbol{F}}_e}(s,t) = \int_{{\varOmega _E}} {(Er\nabla \cdot {{\boldsymbol{\tau}} _p})\delta ({\boldsymbol{X}}(s,t) - {\boldsymbol{x}}){\mathrm{d}}{\boldsymbol{x}}} \\ & {F_{e,x}}(t) = \frac{1}{T}\int_t^{t + T} {\frac{1}{L}\int_{{\varOmega _L}} {{{\boldsymbol{F}}_E}(s,t') \cdot {\boldsymbol{e}}{}_x{\mathrm{d}}s} {\text{ }}{\mathrm{d}}t'} \end{aligned} \right\}$$ (27) 分别给出典型工况下的黏性力分量$ {F_{v,x}} $, 压力梯度分量$ {F_{p,x}} $以及弹性力分量$ {F_{e,x}} $, 如图4和图5所示. 压力$ {F_{p,x}} $始终为正值, 这表明压力与微生物游动速度方向一致, 起推力作用. 黏性力为负值, 这表明黏性力为阻力, 阻碍微生物游动. 弹性力$ {F_{e,x}} $相对于前两者则小很多, 几乎可以忽略. 上述结果表明, 在微生物游动过程中, 压力和黏性力起主导作用. 进一步分析, 相比自由游动, 微生物在槽道中受到更大的流体压力以及黏性力. 我们认为这与微生物在槽道中呈现加速现象有关. 由于微生物所处位置改变, 从而引起的微生物受力发生变化, 是导微生物游动速度改变的重要原因. 特别是壁面附近, 如$ {h_u} = 0.3 $工况, 微生物所受到的压力以及黏性力, 显著大于自由游动以及槽道中心$ {h_u} = 0.5 $工况. 比较$ Pe = 1,10 $工况下的结果, 可以发现微生物所受到的压力与黏性力差异不大. 对于弹性力, 液晶取向$ \theta = 0 $工况下, 微生物所受到的弹性力为正值, 对微生物游动起推动作用; 而$ \theta = \text{π} /3 $工况下, 弹性力为负值, 起阻碍作用. 这表明在微生物游动过程中, 弹性力性质与液晶取向有关. 上述结果与文献[27-28]中的结论一致. 由图3 ~ 图5的结果可以看出, 微生物在壁面附近位置的游动速度高于槽道中心处, 此时微生物受到更大的推力, 我们认为这与微生物受到更强的壁面约束有关. 最后, 我们给出微生物所受到的水动力的合力. 如图所示, 微生物所受到的合外力接近于0, 这表明研究的结果是可靠的.
![]()
图 4 ξ = 8, $ \theta = 0 $工况下, 微生物受到的水动力Figure 4. The hydrodynamic forces on microorganisms when ξ = 8, $ \theta = 0 $![]()
图 5 ξ = 8, $ \theta = \text{π} /3 $工况下, 微生物受到的水动力Figure 5. The hydrodynamic forces on microorganisms whenξ = 8, $ \theta = \text{π} /3 $接下来分析微生物的游动功率以及效率. 定义微生物在流体中的游动功率如下[37]
$$ Po = \frac{1}{L}\int_{{\varOmega _L}} {{\boldsymbol{u}} \cdot {{\boldsymbol{f}}_e}{\mathrm{d}}s} $$ (28) 引入功率比值, 定义为$ P{o_w}/P{o_f} $; 其中$ P{o_w} $表示微生物在槽道中游动的功率, $ P{o_f} $表示自由游动的功率. 图6中给出了微生物游动功率比值云图. 由图中结果可以看出, 游动功率比值总是大于1; 说明与自由游动相比, 微生物在槽道中游动功率增加. 相比槽道中心, 在近壁面处微生物游动功率更大. 与游动速度不同的是, 微生物游动功率对液晶取向$ \theta $的变化并不敏感, 几乎没有明显的变化.
![]()
图 6 游动功率比值$ P{o_w}/P{o_f} $云图, 横轴为微生物与槽道上壁面间的距离, 纵轴为液晶取向; 图中虚线为等值线Figure 6. The contour of the energy consumption ratio $ P{o_w}/P{o_f} $, the horizontal axis represents the distance between the microswimmer and the upper wall of the channel, and the vertical axis represents the liquid crystal orientation最后, 分析微生物在槽道中的游动效率. 基于游动速度和功率, 定义游动效率[37]为
$$ \eta = \frac{{{U^2}}}{{Po}} $$ (29) 游动效率比值定义为$ {\eta _w}/{\eta _f} $, $ {\eta _w} $表示微生物在槽道中的游动效率, $ {\eta _f} $表示自由游动时的游动效率. 如图7(a) ~ 图7(b)所示, 在$ Er = 0.1 $工况下, 微生物游动效率在距离槽道壁面0.1 ~ 0.2处达到极大值, 液晶取向对其影响较小. $ Er = 1 $工况下, 随着微生物所处位置逐渐接近壁面, 其游动效率表现出非单调性. $ Pe = 1 $时, 如图7(c)所示, 在距离壁面$ {h_u} \approx 0.2 $处, 游动效率达到极大值, 大约为自由游动效率的5倍以上. $ Pe = 10 $时, 微生物在壁面附近的游动效率仍显著高于槽道中心处. 但与$ Pe = 1 $工况不同, 当液晶取向$ \theta $接近$ \text{π} /3 $时, 游动效率达到极大值, 大约5倍左右.
![]()
图 7 游动效率比值$ {\eta _w}/{\eta _f} $云图, 横轴为微生物与槽道上壁面间的距离, 纵轴为液晶取向; 图中虚线为等值线 (续)Figure 7. The contour of the swimming efficiency ratio $ {\eta }_{w}/{\eta }_{f} $, the horizontal axis represents the distance between the microswimmer and the upper wall of the channel, and the vertical axis represents the liquid crystal orientation (continued)4. 结论
本文采用浸没边界法, 研究槽道中微生物在各向异性黏弹性流体中的游动特性以及内在的水动力制. 结果表明, 在本文研究参数范围内, 相比自由游动, 微生物在槽道中呈现加速现象, 而且加速效果主要与微生物和壁面间的距离相关. 在槽道壁面附近, 微生物加速效果更著. 在游动过程中, 流体压力为推力, 黏性力为阻力. 这两者起主导作用, 弹性力较小, 几乎可以忽略. 此外, 微生物在槽道中的游动功率和效率, 相较于自由游动工况均会增大, 在槽道壁面附近增加更显著. 游动效率呈现一定的非单调性, 极大值达到5倍以上.
未来将在本文基础上, 进一步研究复杂界面附近微生物的行为特征. 复杂界面的壁面效应一直是微生物游动问题的难点. 特别是壁面附近涉及微生物聚集[38-39], 使得问题更具挑战性. 研究复杂界面附近微生物的行为特征及其内在的动力学机理, 可以为微型元器件的生产制造提供理论支持, 在靶向药物的定向输送方面有巨大的应用潜力[40].
数据可用性声明
支撑本研究的科学数据已在中国科学院科学数据银行 ScienceDB 平台公开发布, 访问地址为: https://www.doi.org/10.57760/sciencedb.j00140.00043.
-
![]()
图 3 速度比值$ {U_w}/{U_f} $云图, 横轴为微生物与槽道上壁面间的距离, 纵轴为液晶取向; 图中虚线为等值线
Figure 3. The contour of the speed ratio $ {U}_{w}/{U}_{f} $, the horizontal axis represents the distance between microswimmer and upper wall of the channel, and the vertical axis represents the liquid crystal orientation
![]()
图 4 ξ = 8, $ \theta = 0 $工况下, 微生物受到的水动力
Figure 4. The hydrodynamic forces on microorganisms when ξ = 8, $ \theta = 0 $
![]()
图 5 ξ = 8, $ \theta = \text{π} /3 $工况下, 微生物受到的水动力
Figure 5. The hydrodynamic forces on microorganisms whenξ = 8, $ \theta = \text{π} /3 $
![]()
图 6 游动功率比值$ P{o_w}/P{o_f} $云图, 横轴为微生物与槽道上壁面间的距离, 纵轴为液晶取向; 图中虚线为等值线
Figure 6. The contour of the energy consumption ratio $ P{o_w}/P{o_f} $, the horizontal axis represents the distance between the microswimmer and the upper wall of the channel, and the vertical axis represents the liquid crystal orientation
![]()
图 7 游动效率比值$ {\eta _w}/{\eta _f} $云图, 横轴为微生物与槽道上壁面间的距离, 纵轴为液晶取向; 图中虚线为等值线 (续)
Figure 7. The contour of the swimming efficiency ratio $ {\eta }_{w}/{\eta }_{f} $, the horizontal axis represents the distance between the microswimmer and the upper wall of the channel, and the vertical axis represents the liquid crystal orientation (continued)
-
[1] Gutman E, Or Y. Symmetries and gaits for Purcell's three-link microswimmer model. IEEE Transactions on Robotics, 2015, 32(1): 53-69
[2] Cheang UK, Roy D, Lee JH, et al. Fabrication and magnetic control of bacteria-inspired robotic microswimmers. Applied Physics Letters, 2010, 97(21): 213704 doi: 10.1063/1.3518982
[3] Magdanz V, Medina-Sánchez M, Schwarz L, et al. Spermatozoa as functional components of robotic microswimmers. Advanced Materials, 2017, 29(24): 1606301 doi: 10.1002/adma.201606301
[4] Taylor GI. Analysis of the swimming of microscopic organisms. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 1951, 209(1099): 447-461
[5] Kumar MS, Philominathan P. The physics of flagellar motion of E. coli during chemotaxis. Biophys Rev, 2010, 2: 13-20
[6] Reynolds AJ. The swimming of minute organisms. Journal of Fluid Mechanics, 1965, 23(2): 241-260 doi: 10.1017/S0022112065001337
[7] Bianchi S, Saglimbeni F, Di Leonardo R. Holographic imaging reveals the mechanism of wall entrapment in swimming bacteria. Physical Review X, 2017, 7(1): 011010 doi: 10.1103/PhysRevX.7.011010
[8] Ezhilan B, Alonso-Matilla R, Saintillan D. On the distribution and swim pressure of run-and-tumble particles in confinement. Journal of Fluid Mechanics, 2015, 781: R4 doi: 10.1017/jfm.2015.520
[9] Bechinger C, Di Leonardo R, Löwen H, et al. Active Brownian particles in complex and crowded environments. Rev Mod Phys, 2016, 88(4): 045006 doi: 10.1103/RevModPhys.88.045006
[10] Galajda P, Keymer J, Chaikin P, et al. A wall of funnels concentrates swimming bacteria. Journal of Bacteriology, 2007, 189(23): 8704-8707 doi: 10.1128/JB.01033-07
[11] Frymier PD, Ford RM, Berg HC, et al. Three-dimensional tracking of motile bacteria near a solid planar surface. Proceedings of the National Academy of Sciences, 1995, 92(13): 6195-6199 doi: 10.1073/pnas.92.13.6195
[12] Li G, Tang JX. Accumulation of microswimmers near a surface mediated by collision and rotational Brownian motion. Physical Review Letters, 2009, 103(7): 078101 doi: 10.1103/PhysRevLett.103.078101
[13] Katz DF. On the propulsion of micro-organisms near solid boundaries. Journal of Fluid Mechanics, 1974, 64(1): 33-49 doi: 10.1017/S0022112074001984
[14] Rajendran RR, Banerjee A. Effect of non-Newtonian dynamics on the clearance of mucus from bifurcating lung airway models. Journal of Biomechanical Engineering, 2021, 143(2): 021011 doi: 10.1115/1.4048474
[15] Molla MM, Paul MC. LES of non-Newtonian physiological blood flow in a model of arterial stenosis. Medical Engineering & Physics, 2012, 34(8): 1079-1087
[16] Dharmaiah G, Prasad JL R, Balamurugan KS, et al. Performance of magnetic dipole contribution on ferromagnetic non-Newtonian radiative MHD blood flow: An application of biotechnology and medical sciences. Heliyon, 2023, 9(2): E13369 doi: 10.1016/j.heliyon.2023.e13369
[17] Bhatti MM, Sait SM, Ellahi R. Magnetic nanoparticles for drug delivery through tapered stenosed artery with blood based non-Newtonian fluid. Pharmaceuticals, 2022, 15(11): 1352 doi: 10.3390/ph15111352
[18] Faghiri S, Akbari S, Shafii MB, et al. Hydrothermal analysis of non-Newtonian fluid flow (blood) through the circular tube under prescribed non-uniform wall heat flux. Theoretical and Applied Mechanics Letters, 2022, 12(4): 100360 doi: 10.1016/j.taml.2022.100360
[19] Merchant K, Rill RL. DNA length and concentration dependencies of anisotropic phase transitions of DNA solutions. Biophysical Journal, 1997, 73(6): 3154-3163 doi: 10.1016/S0006-3495(97)78341-3
[20] Nayani K, Evans AA, Spagnolie SE, et al. Dynamic and reversible shape response of red blood cells in synthetic liquid crystals. Proceedings of the National Academy of Sciences, 2020, 117(42): 26083-26090 doi: 10.1073/pnas.2007753117
[21] Gagnon DA, Keim NC, Arratia PE. Undulatory swimming in shear-thinning fluids: experiments with Caenorhabditis elegans. Journal of Fluid Mechanics, 2014, 758: R3 doi: 10.1017/jfm.2014.539
[22] Lauga E. Propulsion in a viscoelastic fluid. Physics of Fluids, 2007, 19(8): 083104 doi: 10.1063/1.2751388
[23] Ives TR, Morozov A. The mechanism of propulsion of a model microswimmer in a viscoelastic fluid next to a solid boundary. Physics of Fluids, 2017, 29(12): 121612 doi: 10.1063/1.4996839
[24] Krieger MS, Spagnolie SE, Powers T. Microscale locomotion in a nematic liquid crystal. Soft Matter, 2015, 11(47): 9115-9125 doi: 10.1039/C5SM02194D
[25] Krieger MS, Spagnolie SE, Powers TR. Swimming with small and large amplitude waves in a confined liquid crystal. Journal of Non-Newtonian Fluid Mechanics, 2019, 273: 104185 doi: 10.1016/j.jnnfm.2019.104185
[26] Doi M. Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropic and liquid crystalline phases. Journal of Polymer Science: Polymer Physics Edition, 1981, 19(2): 229-243 doi: 10.1002/pol.1981.180190205
[27] Lin Z, Chen S, Gao T. Q-tensor model for undulatory swimming in lyotropic liquid crystal polymers. Journal of Fluid Mechanics, 2021, 921: A25 doi: 10.1017/jfm.2021.531
[28] Lin Z, Yu Z, Li J, et al. Anisotropic swimming and reorientation of an undulatory microswimmer in liquid-crystalline polymers. Journal of Fluid Mechanics, 2022, 946: A30 doi: 10.1017/jfm.2022.621
[29] Sonnet AM, Maffettone PL, Virga EG. Continuum theory for nematic liquid crystals with tensorial order. Journal of Non-Newtonian Fluid Mechanics, 2004, 119(1-3): 51-59 doi: 10.1016/j.jnnfm.2003.02.001
[30] Feng JJ, Sgalari G, Leal LG. A theory for flowing nematic polymers with orientational distortion. Journal of Rheology, 2000, 44(5): 1085-1101 doi: 10.1122/1.1289278
[31] Feng J, Chaubal CV, Leal LG. Closure approximations for theDoi theory: Which to use in simulating complex flows of liquid-crystalline polymers? Journal of Rheology, 1998, 42(5): 1095-1119
[32] Peskin CS. The immersed boundary method. Acta Numerica, 2002, 11: 479-517 doi: 10.1017/S0962492902000077
[33] Namiki T. A new FDTD algorithm based on alternating-direction implicit method. IEEE Transactions on Microwave Theory and Techniques, 1999, 47(10): 2003-2007 doi: 10.1109/22.795075
[34] Yu Z. A DLM/FD method for fluid/flexible-body interactions. Journal of Computational Physics, 2005, 207(1): 1-27 doi: 10.1016/j.jcp.2004.12.026
[35] Vaithianathan T, Collins LR. Numerical approach to simulating turbulent flow of a viscoelastic polymer solution. Journal of Computational Physics, 2003, 187(1): 1-21 doi: 10.1016/S0021-9991(03)00028-7
[36] Salazar D, Roma A, Ceniceros H. Numerical study of an inextensible, finite swimmer in Stokesian viscoelastic flow. Phys Fluids, 2016, 28(6): 063101
[37] Teran J, Fauci L, Shelley M. Viscoelastic fluid response can increase the speed and efficiency of a free swimmer. Physical Review Letters, 2010, 104(3): 038101 doi: 10.1103/PhysRevLett.104.038101
[38] Park Y, Kim Y, Lim S. Flagellated bacteria swim in circles near a rigid wall. Physical Review E, 2019, 100(6): 063112 doi: 10.1103/PhysRevE.100.063112
[39] Takaha Y, Nishiguchi D. Quasi-two-dimensional bacterial swimming around pillars: Enhanced trapping efficiency and curvature dependence. Physical Review E, 2023, 107(1): 014602 doi: 10.1103/PhysRevE.107.014602
[40] Park BW, Zhuang J, Yasa O, et al. Multifunctional bacteria-driven microswimmers for targeted active drug delivery. ACS Nano, 2017, 11(9): 8910-8923 doi: 10.1021/acsnano.7b03207
下载:
计量
- 文章访问数: 287
- HTML全文浏览量: 48
- PDF下载量: 55
