在回转黏度仪中，水泥净浆的流变特性参数通常是通过雷诺-里符林公式^[1-3]计算. 通过黏度仪所测得的扭矩-转速曲线，可以计算出水泥净浆的黏度和屈服应力.

水泥净浆通常被认为是一种宾汉姆流体^[4-6]；同时，水泥净浆具有触变性^[7-11]，是一种剪切变稀的流体，在剪切作用下，其表观黏度会下降. Hattori等^[12-13]认为水泥净浆的表观黏度与单位体积的可逆键连接数量有关，并将单位体积的可逆键连接数表示成了剪切速率和絮结速率以及时间的函数，建立了服部-泉理论（Hattori-Izumi theory，HI理论）. Wallevik等^[14-15]对上述函数进行了修正，用2个记忆模量代替了剪切速率和絮结速率，同时对屈服应力和表观黏度的计算方式进行了修正，建立了修正的HI理论，即颗粒流交互理论（particle flow interaction theory，PFI理论）.

回转黏度仪中水泥净浆的黏度并非均匀分布，每个位置的黏度和屈服应力均与其所受的剪切速率有关，因此，回转黏度仪中水泥净浆的流场为非线性流场，雷诺-里符林公式无法准确计算不同位置的黏度和屈服应力. Wallevik等^[15-16]基于PFI理论，采用有限差分方法^[17-20]，对黏度仪内的净浆特性进行了仿真，并进行了相应的实验，仿真结果和对应的实验结果显示了高度的一致性. 由于有限差分方法的基础为计算流体力学，需要对整个流场区域进行网格划分，且每个网格点都需要计算多个变量的差分，仿真时间较长，很难满足实验时标定仿真参数的要求.

基于理论流体力学，本文提出一种迭代算法，用于快速计算回转黏度仪流场中水泥净浆的流变特性参数，以及水泥净浆PFI参数的标定. 针对此类非线性流场，本文根据回转黏度仪中水泥净浆流场中心对称的性质，对流场进行径向分层，假设层内的流体具有均匀性和连续性，对每层运用雷诺-里符林公式计算剪切速率，进而计算出整个流场中水泥净浆的流变特性参数.

水泥净浆的黏度主要受时间、剪切和温度等因素的影响. 本文研究当温度为20~25°C时，在回转黏度仪叶片剪切作用下，水泥净浆的表观黏度在流场中的分布和变化情况.

流变性是指物质在外力作用下的变形和流动性质，对于流体来说，就是施加于流体上的剪切应力 $\tau $与剪切变形率（剪切速率） $\dot \gamma $之间的定量关系.

按照是否符合牛顿内摩擦定律，可以将流体分为牛顿流体和非牛顿流体. 牛顿流体是一种流变性与时间无关的纯黏性流体，其流变方程为

(1) $ \tau = {\eta _{{}_{\rm N}}}\dot \gamma. $

式中： ${\eta _{{}_{\rm N}}}$为牛顿流体的动力黏度，简称黏度，是剪切应力与剪切速率的比值，是一个恒定值.

在非牛顿流体中，按照流变性是否与剪切持续时间相关，可以把非牛顿流体分成2类：非时变性非牛顿流体和时变性非牛顿流体. 静态下水泥净浆的近似流变模型—宾汉姆流体即为一种非时变性非牛顿流体，其流变方程如下：

(2) $ \tau = {\tau _0} + \mu \dot \gamma . $

式中： $\mu $为宾汉姆黏度，又称为塑性黏度， ${\tau _0}$为屈服应力.

在非牛顿流体中，剪切应力与剪切速率的比值被定义为黏度 $\eta $，称为表观黏度，宾汉姆流体的表观黏度为

(3) $ \eta = \frac{\tau }{{\dot \gamma }}{\rm{ = }}\frac{{{\tau _0}}}{{\dot \gamma }}{\rm{ + }}\mu. $

由式（3）可知，表观黏度是一个与剪切速率 $\dot \gamma $有关的非恒定值.

宾汉姆流体的流变特性参数可通过回转黏度仪测定，对于回转黏度仪存在以下公式：

(4) $ M = 2{\text{π}} {r^2}h\tau = 2{\text{π}} {r^2}h\left({\tau _0}{\rm{ + }}\mu r\frac{{{\rm d}\omega }}{{{\rm d}r}}\right). $

式中： $M$为黏度仪旋转时的推动力矩，h为内筒在物料中的浸没深度， $r$为半径， $\omega $为半径r处的转速. 进一步推导可以得到

(5) $ \varOmega = \frac{M}{{4{\text{π}} h\mu }}\left(\frac{1}{{{R_1}^2}} - \frac{1}{{{R_2}^2}}\right) - \frac{{{\tau _0}}}{\mu }\ln \frac{{{R_2}}}{{{R_1}}}. $

式中：R₁、R₂分别为黏度仪内、外筒半径， $\varOmega $为黏度仪内筒转速，式（5）被称为雷诺-里符林公式.

雷诺-里符林公式为给定内筒转速下回转黏度仪水泥净浆流场中不同半径处的剪切速率计算提供了方法.

水泥净浆具有触变性，即水泥净浆的黏度会在持续的剪切作用下渐渐变小，剪切作用消失后，黏度会逐渐恢复.

Hattori^[12]从微观结构出发，对水泥净浆的触变性作出了如下解释：在水泥与水的水合反应过程中，相互接触的颗粒之间会出现絮结，形成连接键，由范德华力形成的连接键可被剪切破坏，称为可逆键，由于水合反应形成的连接键不可被剪切破坏，被称为不可逆键. 在受到剪切作用时，水泥净浆中的可逆键被破坏，水泥颗粒重新成为自由粒子，在宏观上导致水泥净浆的黏度减小，呈现剪切变稀的性质. 剪切作用消失后，由于范德华力的作用，自由的水泥颗粒之间会形成可逆连接键，宏观上导致水泥净浆的黏度逐渐恢复.

HI理论认为水泥净浆的表观黏度与单位体积的可逆键连接数量有关，即

(6) $ {\eta _{\rm{HI}}} = {B_3}J_{\rm t}^{2/3} + {\text{其他可忽略的相关因素}} \approx {B_3}J_{\rm t}^{2/3}. $

式中： ${B_3}$为水泥净浆中基本粒子之间的摩擦系数， ${J_{\rm t}}$为单位体积内水泥颗粒之间的可逆键连接数量. 定义单位体积的水泥净浆中的水泥粒子数为 ${n_3}$，则 ${J_{\rm t}}$可通过下式计算：

(7) $ {J_{\rm t}} = \frac{{{n_3}\left[{U_0}(\dot \gamma H{t^2}{\rm{ + }}1){\rm{ + }}Ht\right]}}{{{\rm{(}}Ht{\rm{ + 1)(}}\dot \gamma t{\rm{ + }}1{\rm{)}}}}. $

式中： ${U_0}$为 $t = 0$时刻与单个水泥颗粒相连的平均可逆键的数量，H为絮结速率常数，t为时间.

由式（6）和式（7）得到表观黏度为

(8) $ {\eta _{{}_{\rm{HI}}}} = {B_3}{J_{\rm t}}^{2/3} = {B_3}{\left\{ {\frac{{{n_3}\left[{U_0}(\dot \gamma H{t^2}{\rm{ + }}1){\rm{ + }}Ht\right]}}{{{\rm{(}}Ht{\rm{ + 1)(}}\dot \gamma t{\rm{ + }}1{\rm{)}}}}} \right\}^{2/3}}. $

由于剪切速率 $\dot \gamma $和絮结速率常数大多数情况下是一个变化的值，式（8）在计算中应用价值不大. Wallevik^[14]对式（8）进行了修正，使其计算更方便，计算结果拟合效果更好，修正后的HI理论被称为PFI理论. 修正过程如下.

1）用2个记忆模量 $\tilde \varGamma $和 $\tilde \varTheta $分别代替式（8）中的 $\dot \gamma t$和Ht. $\tilde \varGamma $和 $\tilde \varTheta $的计算公式分别为

(9) $ \tilde \varGamma = \int_0^t {\alpha (t - t')} \dot \gamma (t'){\rm d}t', $

(10) $ \tilde \varTheta = \int_0^t {\beta (t - t')} H(t'){\rm d}t', $

(11) $ {J_{\rm t}} = \frac{{{n_3}\left[{U_0}(\tilde \varGamma \tilde \varTheta {\rm{ + }}1){\rm{ + }}\tilde \varTheta \right]}}{{{\rm{(}}\tilde \varTheta {\rm{ + 1)(}}\tilde \varGamma {\rm{ + }}1{\rm{)}}}}. $

式中：

(12) $ \alpha (t - t') = {\exp\;[{ - (t - t')/{m_{\rm a}}}]}, $

(13) $ \beta (t - t') = {\exp\;[{ - (t - t')/{m_{\rm b}}}]}. $

其中， $\alpha $和 $\beta $为记忆函数，用以表征过去的状态对当前剪切黏度的影响， ${m_{\rm a}}$和 ${m_{\rm b}}$为对应的记忆参数，用以表征剪切黏度能记忆过去状态的程度.

2）修正屈服应力和表观黏度的计算方式. 塑性黏度和屈服应力的计算公式如下：

(14) $ \mu = {a_1}{\eta _{{}_{\rm{HI}}}} + \widetilde {{\eta _0}} = \widetilde {{\eta _0}} + {a_1}{B_3}{J_{\rm t}}^{2/3}, $

(15) $ {\tau _0} = {a_2}{\eta _{{}_{\rm{HI}}}} + \widetilde {{\tau _0}} = \widetilde {{\tau _0}} + {a_2}{B_3}{J_{\rm t}}^{2/3}. $

式中： $\widetilde {{\eta _0}}$ 为初始黏度， $\widetilde {{\tau _0}}$为初始屈服应力，a₁、a₂均为与水泥颗粒之间形成可逆键连接时产生的与动量交换有关的常数.

3）修正絮结速率的计算方法，计算公式如下：

(16) $ \left. \begin{array}{*{20}{l}} H(\dot \gamma ,t) = \dfrac{{K(t)}}{{{{\dot \gamma }^2} + l}},\;\;\;\;\;\;\;\;\;&t > 0;\\ H(\dot \gamma ,0) =\dfrac{{{k_1}(1 - {U_0})}}{{4l}},&t = 0. \end{array} \right\} $

式中：l为一个经验常数，恒等于1 s⁻². K(t)为经验函数，计算方法如下：

(17) $ K(t) = \left\{\begin{array}{*{20}{l}} {k_1} = 0.005\;\,{{\rm s}^{ - 3}},& {\rm d}\omega (t)/{\rm d}t \geqslant 0; \\ {k_2} = 0.1\;\,{{\rm s}^{ - 3}},& {\rm d}\omega (t)/{\rm d}t \leqslant 0 . \\ \end{array} \right. $

PFI理论将水泥净浆的塑性黏度和屈服应力表示成了剪切速率和时间的函数，为受剪切作用时水泥净浆的黏度和屈服应力的计算提供了方法，参数 $\widetilde {{\tau _0}}$、 $\widetilde {\eta_0} $、 ${a_1}$、 ${a_2}$、 ${U_0}$、 ${m_{\rm a}}$、 ${m_{\rm b}}$、 ${B_3}$、 ${n_3}$由水泥净浆的自身特性决定，称为PFI参数. PFI参数可通过水泥净浆的试验数据进行标定，本文中PFI参数的确定方法与文献[8]中所用方法相同.

基于以上理论，本文所述迭代算法先通过雷诺-里符林公式计算流场中不同半径处的剪切速率，再通过PFI理论计算流场中不同半径处的塑性黏度和屈服应力.

在采用回转黏度仪测量水泥净浆的流变特性时，转子运动带动附着于其表面的水泥净浆运动，在整个测量区域内形成流场. 通过雷诺-里符林公式可以推导出整个流场区域内的剪切速率的计算公式如下：

(18) $ \dot \gamma = \frac{2}{{{r^2}}}\left(\frac{1}{{{R_1}^2}} - \frac{1}{{{R_2}^2}}\right)\left[\varOmega + \frac{{{\tau _0}}}{\mu }\ln \frac{{{R_1}}}{{{R_2}}}\right] - \frac{{{\tau _0}}}{\mu }. $

式中： $\dot \gamma $为半径r处的剪切速率.

由式（18）可见，测量区域中的剪切速率 $\dot \gamma $随着半径r、屈服应力 ${\tau _0}$、塑性黏度 $\mu $的变化而变化. 因此，整个流场区域内水泥净浆的特性是不均匀的，黏度仪内的流场为非线性流场，不能通过雷诺-里符林公式来对整个流场区域的剪切速率进行计算.

Wallevik^[16]利用有限差分法，结合PFI理论对回转黏度仪中的水泥净浆流场进行了仿真，利用Taylor展开，取一阶截断，得到如下公式：

(19) $ \frac{{\partial f}}{{\partial x}} = {\frac{1}{{\Delta x}} }\left[{{f\left({x_0} + \dfrac{{\Delta x}}{2}\right) - f\left({x_0} - \dfrac{{\Delta x}}{2}\right)}}\right]+ O(\Delta {x^2}). $

式（19）被称为空间中心差分，对整个流场区域进行网格划分，在格点上应用式（19）得到速度的更新公式（式（20））. 如图1所示，v为每个网格点的速度，i为黏度仪转子半径方向的单位向量，j为黏度仪转子轴向.

(20) $\begin{gathered} { v}_{i,j}^{k + {1 / 2}} - { v}_{i,j}^k = \frac{\,\beta\, }{\,\rho\, }\left[T_{i + {0.5},j}^{k + {0.5}} - T_{i - {0.5},j}^{k + {0.5}}\right] + \frac{{{\theta _i}}}{\rho }T_{i,j}^{k + {0.5}} + \\ \frac{\,\chi \,}{\,\rho\, }\left[Z_{i,j + {0.5}}^k - Z_{i,j - {0.5}}^k\right]. \\ \end{gathered} $

式中：

(21) ${\theta _i}{\rm{ = }}\Delta t/{r_i},\;\;\beta {\rm{ = }}\Delta t/\left( {2\Delta r} \right),\;\; \chi {\rm{ = }}\Delta t/\left( {2\Delta z} \right),$

${Z_{i,j +{0.5}}} = {\eta _{i,j + {0.5}}}\frac{{{{ v}_{i,j + 1}} - {{ v}_{i,j}}}}{{\Delta z}} = \frac{{{\eta _{i,j + {0.5}}}}}{{\Delta z}}{{ v}_{i,j + 1}} - \frac{{{\eta _{i,j + {0.5}}}}}{{\Delta z}}{{ v}_{i,j}},$

${Z_{i,j - {0.5}}} = {\eta _{i,j - {0.5}}}\frac{{{{ v}_{i,j}} - {{ v}_{i,j - 1}}}}{{\Delta z}} = \frac{{{\eta _{i,j - {0.5}}}}}{{\Delta z}}{{ v}_{i,j}} - \frac{{{\eta _{i,j - {0.5}}}}}{{\Delta z}}{{ v}_{i,j - 1}}.$

式中： $\Delta t$为时间步长， $\Delta r$为径向步长， $\Delta z$为轴向步长， ${ v}_{i,j}^k$为半径为 $i\Delta r$、轴向高度为 $j\Delta z$位置网格点在 $t = k\vartriangle t$时刻的速度， ${\eta _{i,j - {0.5}}}$为半径为 $i\Delta r$、轴向高度为 $(j - {0.5})\Delta z$位置网格点水泥净浆的黏度.

Wallevik^[16]比较了有限差分方仿真结果与试验结果，发现仿真曲线与试验曲线相似度极高，验证了有限差分方法仿真的准确性.

在每个时间步长内，需要对流场区域内的速度、转矩、黏度等参数进行计算，以得到一个收敛的解. 如果得到的计算结果不收敛，则本次仿真失败，需要使用更小的时间步长进行仿真. 因此整个仿真的耗时较长，效率较低，很难满足实验时标定水泥净浆的PFI参数的要求.

针对三维仿真计算存在的问题，根据回转黏度仪内部的水泥净浆运动特性，提出二维状态下回转黏度仪内部各点黏度及剪切速率计算方法.

首先对回转黏度仪的结构及其测量区域内水泥净浆的流动状态进行分析. 为了提高黏度仪的测量精度，黏度仪转子通常采用上、下部分离结构，如图2所示.

回转黏度仪在采集数据时，只采集转子上部的数据. 由于机械结构的分离，转子上、下部分由不同的驱动结构带动，转子下部复杂的流场分布对转子数据采集的影响被消除. 在这种状况下，可以认为对于转子上部而言，从上至下，流体对转子的作用力是均匀的，因此可以通过流场区域的径向分层来简化计算过程.

如图3所示，基于回转黏度仪流场中心对称的特性，在黏度仪半径方向上，将流场划分为许多很薄的剪切层，并认为划分之后的剪切层符合以下条件：

1）可认为层内流体的流体特性（黏度和屈服应力等）相同；

2）层与层之间接触部分的转速相同；

3）整个流场区域内所有的剪切层受到的转矩为同一值（转子转动的驱动转矩）.

4）转子表面的流体与转子转速相同.

根据以上条件可以构建简单的迭代算法，在层内运用雷诺-里符林公式计算剪切速率，当整个流场符合以上4个条件时，即可认为迭代算法收敛，计算得到的转速符合实际流场的剪切运动. 具体计算步骤如图4所示.

在仿真过程中，流场中计算得到的转速分布状况有4种，不同转速分布状况的处理方式如下.

1）若转速由内至外不断减小且均为正值，最外层的转速小于阈值，则可认为本次选取的转矩值与真实流场的转矩值误差在可接受范围内，停止迭代.

2）若转速由内至外不断减小且均为正值，最外层的转速大于阈值，则可认为本次选取的转矩值小于真实流场的转矩值，须根据误差增大转矩值.

3）若转速由内至外不断减小，出现转速为负的状况，则可认为本次选取的转矩值大于真实流场的转矩值，须根据误差减小转矩值.

4）若转速由内至外增加，则说明转速太小，靠近外部的水泥净浆由于所受剪切应力小于屈服应力而静止，此时认为出现阻塞流，须将出现阻塞流位置的半径作为外筒半径代入算法，重新计算转矩.

单个时间步的仿真结束后，可以通过PFI理论计算这个时间步的水泥净浆特性，然后根据当前的水泥净浆特性，计算下一个时间步的剪切速率，直至仿真结束，具体步骤如图5所示.

本文将迭代算法的仿真结果分别与文献[16]中有限差分法的仿真结果以及试验结果进行对比，以验证迭代算法仿真结果的准确性.

根据文献[16]中所提供的参数，设置仿真参数和黏度仪转速. 如图6所示为仿真中黏度仪的转速曲线，仿真参数设置如表1所示.

图7为迭代算法仿真结果与有限差分法仿真结果的对比，图中，T为转矩. 由图7可知，迭代算法仿真得到的结果与文 献[16]利用有限差分法仿真得到的结果极为吻合. 经计算，2条曲线的平均误差率为0.89%. 该迭代算法是从理论流体力学出发，通过对流场的径向分层，假设每层内轴向的水泥净浆的特性相同，将三维问题转化为二维问题，避免了网格划分所带来的巨大计算量，迭代算法的计算速度大大提高. 本次仿真时长仅为2.55 s，相较于有限差分法(仿真时长为5 d)，大大减小了仿真时长，提高了运算效率，因此该算法可用于水泥净浆PFI参数的标定试验，对水泥净浆的PFI参数进行实时标定.

为进一步验证迭代算法仿真结果的准确性，本文采用Brookfield-DV2TLV型号回转黏度仪对牌号为PO52.5、水灰比为0.4的水泥净浆进行变转速剪切试验，试验黏度仪转速曲线如图8所示，采用图8中的曲线作为仿真转速曲线，取表2的仿真参数作为迭代算法仿真参数，对变转速试验进行仿真.

图9为试验结果与迭代算法仿真结果的对比. 由图9可见，试验结果与仿真结果吻合度较高，经计算，仿真转矩与试验测试转矩的平均误差率为4.23%，整体误差较小，可认为迭代算法仿真结果准确，可用于回转黏度仪中水泥净浆流场的高效计算.

本文基于PFI理论，提出了一种用于回转黏度仪的水泥净浆流变性分布迭代算法. 该算法采用对流场径向分层的方法，在每层运用雷诺-里符林公式计算剪切速率，并结合PFI理论计算每层的流变性参数，将复杂的三维仿真问题转化为了二维问题，简化了计算过程. 所提算法计算结果与有限差分法仿真结果高度相似，与试验结果的平均误差率为4.23%，计算结果较准确. 相较于传统的有限差分法，该迭代算法将计算时间由5 d缩短到了2.55 s，计算效率较高，可用于水泥净浆PFI参数的实时标定，并对相似的对称流场的宾汉姆流体仿真有一定借鉴价值. 该算法仍存在一定的局限性，在变转速试验转速上升阶段有一定的仿真误差，因此，在后期的工作中须进一步完善转速上升时的计算方式.