一种基于物质浓度的时空梯度分布确定微尺度下变截面扁平微通道高度方向平均流速的方法

本发明属于微尺度下流体流速测量技术领域，涉及一种利用物质的时空浓度梯度分布确定微尺度下变截面扁平微通道高度方向平均流速的方法，是基于流体力学和物质传输原理、光学成像技术、图像处理技术、优化算法以及神经网络算法的检测微流控芯片中变截面扁平微通道内高度方向平均流速的方法。

背景技术：

微流控(microfluidics)技术是近年来发展起来的一种广受关注的多学科交叉领域之一，其通过构建具有不同功能的微纳米尺度的微通道网络单元，实现对微量多相流体的精准操控，复杂生物力学微环境的离体构建，生化反应过程的控制等，已成功应用于分析化学、生物学、医学等领域。微流控芯片系统中各网络单元的功能实现与流体的运动息息相关，微通道中微尺度流体流速分布的测定对于在微流控芯片上进行定量检测、精准操控等都至关重要。

目前，微流控芯片中微尺度流体的测量方法有热膜传感器测量方法，这是一种基于热扩散理论，通过对微通道壁面进行加热，测量热传递来检测微通道壁面附近流速分布的方法。其需要将高精度的传感器植入微流控芯片内，大大提高了微流控芯片的制作成本和制作难度。微粒子图像测速技术(micro-piv)则是基于传统piv技术发展而来用于微尺度流体测速的技术，通过记录示踪微粒在流体中的位移，计算流体的速度。由于其示踪粒子大小与微通道几何尺寸量级相同，因此粒子跟随性和遮挡等问题都会使测量结果产生误差。为克服上述缺陷，申请人提出了一种基于动态荧光粉浓度确定均匀扁平微通道平均流速和剪切力的方法(详见发明专利zl201610139388.2)，然而该方法仅能检测均匀截面扁平微通道纵向平均流速，无法应用于变截面扁平微通道得到二维平面上的平均速度场分布。

本发明针对微流控芯片上具有高度远小于横向和纵向尺寸的变截面扁平微流通道提出一种测量微通道高度方向平均流速分布的方法。

技术实现要素：

本发明是一种利用纳米尺度标记物溶液浓度的时空梯度分布确定微流控芯片中变截面扁平微通道内高度方向流体平均流速分布的方法。在入口处设计能够产生浓度时空梯度分布的生成装置，根据流体力学和物质传输原理，建立描述物质在变截面扁平微通道中对流-扩散现象的对流-扩散方程，设计优化算法利用标记物浓度的时空梯度分布进行逆向求解，得到速度分布。并基于上述工作构建具有物理约束的神经网络框架，通过优化算法和数值仿真的数据对神经网络进行预训练，进一步提高本方法在实施过程中的速度，使其可以应用于对实时性有要求的测量中。

本发明的技术方案：

一种基于物质浓度的时空梯度分布确定微尺度下变截面扁平微通道高度方向平均流速的方法，步骤如下：

设待检测的微流控芯片是满足高度h远小于横向尺寸w和纵向尺寸l的扁平的微通道；向扁平的微通道中加载流量和浓度随时间变化的具有动态浓度时空梯度分布的标记物溶液；微通道中的不可压缩粘性流体运动满足navier-stokes方程和连续性方程：

其中，v表示速度，p表示压力，ρ是流体密度，μ是粘度。

在微通道内，具有动态浓度时空梯度分布的标记物溶液不仅沿着流速方向作对流运动，还会沿着浓度梯度方向进行扩散。微通道内的物质传输满足对流-扩散方程：

其中，c是浓度，d为扩散系数。

忽略侧面通道对流动的影响以及重力作用，仅考虑顶部和底部的表面粘性力对流体运动的影响，则微通道中的navier-stokes方程简化为：

其中t为时间，x,y,z分别为长度、宽度、高度方向的坐标，u，v分别表示x方向和y方向上的速度。扁平的微通道满足低宽高比的几何特征，且微通道中的流体运动为小雷诺数运动，womesley数很小，满足准定常假设。因此微通道中的流速满足：

其中分别表示x和y方向上的高度方向平均流速，所述的高度方向即为z方向，满足：

由于扁平微通道中h＜＜w、h＜＜l，认为标记物溶液在高度方向上会快速形成均匀浓度。因此，对浓度在高度方向上取平均，得到平均浓度：

则扁平的微通道中的taylor-aris弥散方程为：

令其中i、j为笛卡尔坐标下的基向量，则上述taylor-aris弥散方程简记为：

并满足约束方程：

使用有限差分法近似弥散方程，以空间步长δx、δy分别沿着x方向和y方向将待检测的矩形区域离散化，网格点分别为xm，yn，其中m，n分别满足m＝0,1,2,…,m，n＝0,1,2,…,n。以时间步长δt将时间离散化，离散时间点记为tk，其中k＝0,1,2,…,k。使用如下的差分公式近似弥散方程中的各项：

其中表示tk时刻坐标(xm,yn)处的标记物溶液，和分别表示tk时刻坐标(xm,yn)处x方向和y方向上的高度方向平均速度。根据实际测速要求和检测仪器的精度，可以使用更高精度的差分格式。

基于分析得到的差分方程，可以构建求解流场速度分布的优化问题。

由于微流控芯片的入口处的浓度分布随时间变化，可以使得流动方向上的浓度梯度不为0。求解tk时刻的浓度分布对应的tk时刻的速度分布，构建如下的优化问题：

其中，实际运算时采用的是上述方程的差分形式。对于流速分布不随时间变化的定常流情况，可以通过求解不同时刻的速度分布取平均降低速度求解的误差。对于流速分布随时间变化的脉动流情况，代入不同时刻的浓度分布数据重复上述优化问题求解，即代入tk时刻的浓度分布求解tk时刻的速度分布，其中k＝0,1,2,…,k,…,k，即可得到脉动流情况下，速度分布随时间的变化情况。所使用的优化算法应根据实际实验数据进行选取，包括但不限于常见的单纯形法、karmarkar算法、lagrange算法、lemke算法等，通常而言采用经典的单纯形法即可获得较好的效果。

进一步可以通过构建神经网络减少浓度场反演速度场的计算时间，满足对实时性有要求的测量中，可以实现速度分布的实时观测。所使用的神经网络包括输入层、隐藏层与输出层。输入层为时间tk，横向坐标xi，纵向坐标yj以及该时刻对应位置的浓度输出层为相应的浓度横向平均速度纵向平均速度采用全连接的方式连接输入层与输出层，输入层与输出层之间的网络被称为隐藏层，表示为zl，其中l＝1，……，l-1。两个相邻层之间连接如下：

zl＝σl(wl·zl-1+bl)(28)

其中σl为第l层激活函数，wl表示l层的神经网络权重向量，bl表示l层的偏差值。

神经网络的预测结果表示为l层的输出，即：

构建神经网络的损失函数，除了包含神经网络预测值与真实值的差异之外，还通过对输出层求微分，将上述推导的taylor-aris弥散方程(13)与约束方程(14)(15)作为惩罚项加入损失函数中，使之具有物理条件约束。损失函数l(w,b)定义如下：

其中，f是真实值，神经网络的预测值。

通过优化算法与数值计算的获取数据，完成神经网络的训练。在实施过程中，使用光学成像技术得到标记物溶液浓度时空分布的图像，经过图像处理后将其输入训练好的神经网络，即可在实施过程中实时的得到该区域的速度分布

进一步，该方法采用的装置包括动态流量生成装置、动态浓度生成装置、可产生空间浓度梯度分布的圣诞树结构、工控计算机、待检测的微流控芯片、显微镜和废液回收容器。其中，动态流量生成装置由可编程气压泵与密闭储液容器连接构成，通过工控计算机控制气压波形产生动态可控的流体流量；动态浓度生成装置包括可编程控制注射泵、注射器和三通接口，两个注射器中为不同浓度的纳米尺度标记物溶液，通过控制两注射器的流量比即可产生浓度随时间动态变化的标记物溶液。动态流量生成装置和动态浓度生成装置分别接入圣诞树结构的三个入口，在圣诞树的出口形成可控的具有浓度时空梯度分布的标记物溶液，将其通入待检测微流控芯片的微通道中，利用荧光标记光谱成像或无标记光谱成像，如反斯托克斯拉曼光谱cars、红外、拉曼等光学成像技术，实时记录检测区域并得到随时间变化的浓度分布图像。

本发明的有益效果：本方法将物质传输理论与流体力学原理结合起来，建立了变截面扁平微通道中二维浓度场与二维速度场的关系。通过在微通道中形成具有浓度时空梯度分布的标记物溶液，使用优化算法即可基于浓度场反演速度场。进一步利用神经网络算法，可实现实施过程中对所检测区域高度方向平均速度分布的实时观测。

附图说明

图1是变截面扁平微通道示意图。

图2是本发明的装置结构的示意图。图中：1为动态浓度标记物溶液生成装置，其中，1-1是可编程控制注射泵，1-2为注射器，1-3为三通接口；2是动态流量生成装置，其中2-1是一个可编程控制气压泵，2-2是一个密闭储液容器；3是用于控制注射泵与气压泵的工控计算机；4是用于产生空间浓度梯度分布的圣诞树结构(或金字塔结构)；5是待检测的微流控芯片；6是用于光学成像的显微镜；7是废液回收容器。

图3是具有物理约束的全连接人工神经网络模型框架示意图。

图4是具体实施方式中所使用待检测微通道二维结构在极坐标下的示意图。

图5是扁平微通道中浓度分布数据观测示意图。图中8是用于光学成像的显微镜；9是待检测的变截面扁平微通道；10是离散化处理后的浓度分布图像。

图6是具体实施方式中计算速度的结果图，其中(a)为仿真计算的速度矢量，(b)为使用本专利方法基于浓度场反演速度场求解得到的速度矢量，(c)为使用神经网络方法计算的速度矢量。

图7体现了神经网络训练过程中训练集与验证集的损失值随迭代次数的变化，使用圆圈标注了完成训练的位置。

具体实施方式

下面的实施例将对本发明予以进一步说明，但并不因此限制本发明的保护范围。

本实施例使用仿真的方式说明本发明所述方法的有效性，仿真参数设置与实际情况相吻合。仿真模拟了的一种微尺度下变截面扁平微通道高度方向平均流速分布的测量装置，如图2所示，包括可编程气压泵、密闭储液容器、可编程注射泵、注射器、三通接口、导管、圣诞树结构(或金字塔结构)的微通道、显微镜、工控计算机、废液池。

在本实施例中，选择纳米荧光粉溶液作为标记物溶液。如图2所示，在注射器1-1中充满不含荧光粉的缓冲液，注射器1-2中充满浓度为200umol/ml的荧光粉溶液。通过工控计算机3控制可编程注射泵使注射器1-1和注射器1-2的流量比随时间按照一定规律变化，产生浓度随时间正弦变化的荧光粉溶液，同时调整气压泵的输出压力以产生脉动流或定常流。浓度与流量随时间变化的波形表达式如下：

φ(t)＝φ0[1+δsin(2πfφt)](31)

q(t)＝q0[1+λsin(2πfqt)](32)

φ(t)表示动态浓度生成装置产生的动态浓度变化，q(t)则表示微通道中总的流量变化。式(31)中的φ0为平均浓度，fφ为浓度变化的频率，式(32)中的q0为平均流量，fq为流量变化的频率，δ与λ则为相应的无量纲比例因子。当fq为零时即对应定常流情况，反之则对应脉动流情况。

将上述的动态流量生成装置、动态浓度生成装置分别接入圣诞树结构的不同入口，通过圣诞树结构即可在出口处生成具有浓度空间梯度分布的荧光粉溶液，在仿真中设置圣诞树出口处浓度纵向分布为线性梯度。所使用的待检测微流控芯片二维结构如图4中所示，o为极坐标系的原点，a为该微通道的长度，b为入口处的宽度，其边界可由极坐标方程描述如下：

其中θ0与r0分别为：

其中r为极坐标系的极径，θ为极坐标系的极角，n为无量纲的形状参数。该微通道的上下边界与流线重合，出口边界与等势线重合，该微通道具有在底面形成线性切应力梯度分布的特征。该微通道的高度为30μm。

通过数值仿真得到上述微通道中的速度分布与浓度分布情况，仿真参数设置如下表所示：

表1：仿真所使用参数的默认值

仿真模拟基于荧光标记成像技术，使用连接ccd相机的荧光显微镜记录待检测微通道中的荧光粉溶液浓度变化的过程。如图5所示，使用ccd相机记录荧光显微镜视野内某一固定矩形位置不同时刻的荧光强度分布图像，得到时间间隔为△t的荧光强度分布图像。通过图像处理技术，对采集的图像进行去噪与归一化处理，空间离散化后即可得到荧光粉溶液浓度的时空分布数据。

在扁平微通道中，荧光粉溶液浓度的时空分布数据满足说明书中的taylor-aris弥散方程公式(12)。建立求解流场分布的优化问题，如公式(26)(27)所示。

在微通道中建立如图5所示的笛卡尔坐标系，使用专利所提出的方法计算x坐标在300-330μm，y坐标在90-110μm的矩形区域内的速度矢量。所计算出的速度矢量的结果如图6所示，图6中(a)为仿真计算出的速度矢量，(b)为使用本专利方法基于浓度场求解的速度矢量，(c)为使用神经网络方法计算的速度矢量。图6中各点的速度如下表所示：

表2：仿真的速度矢量与基于浓度场求解的速度矢量

定义如下式所示的相对误差衡量本方法求解速度分布的准确性：

上式中表示该点速度的准确值，在本用例中即为仿真计算出的速度，则表示使用本专利方法基于浓度场反演速度场计算得到的该点处的速度。分别计算所选取矩形区域中各离散点的re值，取平均后得到该矩形区域中本专利方法求解速度的平均相对误差为0.18％。

下面验证说明书中所述的神经网络的作用。通过数值方法获取上述矩形区域内在不同时刻的不同浓度分布与不同速度分布的数据，该数据在实际实施中需要通过优化算法计算得到，本实施例中则直接使用了仿真数据来验证神经网络的有效性。

神经网络的结构如图3所示，包括3个隐藏层，每个隐藏层有6个神经元。取10000组数据，将数据集分为三类，80％用于神经网络的训练的训练集，10％用于验证网络是否正在泛化并在过拟合前停止训练的验证集，最后10％将用作神经网络的独立测试的测试集。定义如公式(30)所示的损失函数，使之具有物理条件约束，调用adam优化器训练神经网络，图7为训练过程中训练集与验证集损失函数值随迭代次数的变化，最后使用测试集验证训练后的神经网络求解速度的准确性。使用同等条件的计算机与编译器，对于一个3×3共9点的浓度数据，输入训练好的神经网络得到其速度分布仅耗时0.09s，而使用优化算法的方法计算速度分布则需要20.30s，可见使用训练后的神经网络可以极大缩短计算速度场的时间，在实施中实现对于流场的实时观测。