CN104729786A - 一种孔内灌浆压力的动态测量方法 - Google Patents

Abstract

一种孔内灌浆压力的动态测量方法，包括以下步骤，先是收集灌浆工艺的相关技术参数，然后确定灌浆压力测量的模型结构，基于格子玻尔兹曼方法建立注浆管道内浆液的流体微团的压力分布函数模型，再基于灌浆管道内特定孔深位置的浆液压力分布函数进行求和，得到特定灌浆部位的浆液宏观压力和浆液宏观速度；以注浆管道最下端的浆液宏观压力代表灌浆工艺中的灌浆压力，最后当注浆过程中流量的变化量达到10％或者浆液密度发生变化，则实时更新数据以指导灌浆作业。本发明的孔内灌浆压力的动态测量方法只要知道孔口的速度和浆液的密度值就能获取管道内的任意处压力，操作简单，且数据更加准确可靠。

Description

一种孔内灌浆压力的动态测量方法

技术领域

本发明涉及一种灌浆压力的测量方法，尤其涉及一种在线测量灌浆压力的动态测量方法。

背景技术

灌浆压力是影响灌浆质量的重要参数之一，是现行灌浆监测系统的一个主要监控参数。灌浆过程中，当作用在灌浆段内岩体上的浆液压力超出了地层的极限破坏值，就会引起地层抬动，出现施工事故。《水工建筑物水泥灌浆施工技术规范》(DL/T 5148—2012)中6.3.8中明确要求灌浆段在设计灌浆压力下的操作规则，实质是指灌浆部位的孔内灌浆压力要符合设计压力要求。因此，在灌浆施工中动态获取孔内灌浆压力是准确控制灌浆压力和准确了解施工状态的前提。

在深海采油和深海钻探领域，孔内压力的直接测量已取得初步成果。现在市场上出现了新型的无线智能传感器用于孔内油压的测量和高灵敏的光纤光栅传感器用于水下压力的感测。但由于灌浆工艺的特殊性，孔深有时深达70多米且钻孔的孔径小(仅为56mm-91mm左右)，且浆液粘性、流动性等物理性质与水、油也有很大差异，在灌浆部位安装传感器直接检测孔内灌浆压力，这种方法实施起来非常困难，且对施工效率有很大影响。由于孔内灌浆压力测量的瓶颈技术未解决，现行灌浆监控系统一般直接用孔口压力代替孔内灌浆压力。等的研究(参见L A,Scheid C M,Paraiso E C H,et al.Pressure Drop in CementSlurries Flow in Circular and Annular Regions in Primary Completion[J].Brazilian Journal ofPetroleum and Gas,2014,7(4):129-139.)表明：粘性浆液在小直径长管道中传输时，压力损失是不可以忽视的，直接采用孔口压力值代替孔内灌浆压力值存在较大的误差。

由于获取真实灌浆压力的重要性，部分工程中采用经验修正流体伯努利方程或简单的实验结果来修正灌浆压力，由于灌浆压力受其它因素的影响是复杂的,修正参数需要很多的经验。为了提高孔内灌浆压力模型的非线性预测能力，李凤玲在灌浆压力控制系统的关键技术研究论文中提出了一种智能建模方法[李凤玲.灌浆压力控制系统的关键技术研究[D],中南大学博士论文,2009,55-73.]，获取孔内灌浆压力，利用孔口压力、流量和密度等可测参数建立灌浆压力的自适应T-S(Takagi-Sugeno)模糊软测量模型。这种方法在理论上是可行的，但实际实施起来也存在一定困难。因为它属于一种“黑箱”建模法，前提是需要获取各种地层以及各种灌浆工况下的模型的输入输出数据样本，才能准确建立软测量模型来估计孔内灌浆压力。由于现行灌浆工程中未能获取孔内灌浆压力数据，需开展灌浆试验来获取各种灌浆工艺、浆液材料及地层下灌浆压力建模样本，实验工作量非常大，该方法的可行性较差。另外，该文献中采用Ansys商用软件的流体计算模块(CFD)对灌浆管道摩擦损失开展了数值仿真分析，这种方法只是分析了管道损失压力与管长、管道直径、流体流动速度、浆液密度及粘度的关系，这个模块的理论基础是宏观流体力学模型，需要很多的假设条件，且只针对孔内灌浆压力的影响因素开展离线分析，并没有建立孔内灌浆压力与其它参数的动态数学模型，不能应用于灌浆施工中监控孔内灌浆压力的动态变化，更不能反映复杂浆液对灌浆压力的局部影响。

发明内容

本发明要解决的技术问题是克服现有技术的不足，提供一种可行性强、准确性高、操作方便的在线测量孔内灌浆压力的动态测量方法。

为解决上述技术问题，本发明提出的技术方案为一种孔内灌浆压力的动态测量方法，包括以下步骤，

第一步：根据灌浆工艺，分别通过流量计和密度计获取灌浆过程中实时的流量值Q和密度值ρ₀，并由流量值Q计算得到浆液在注浆管道入口处的速度U_y0；

第二步：确定灌浆压力测量的模型结构：基于格子玻尔兹曼方法建立注浆管道内浆液的流体微团的压力分布函数模型；

第三步：再根据不可压缩流体微团的压力分布函数模型构建灌浆工艺中浆液宏观压力Pm与浆液速度u、浆液密度ρ的泛函关系模型,即其中，f_k(u,ρ)为压力分布函数，且f_k(u,ρ)＝C_sN_k(u,ρ)，另外密度分函数N_k(u,ρ)＝m_kρ(1+3.0*e_k*u+4.5(e_k*u)²-1.5u)]，m_k＝{4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36}；同样根据流体微团的压力分布函数构建灌浆工艺中浆液宏观速度Um与浆液速度u、浆液密度ρ的关系模型,即其中，k＝1,2,3,4,5,6,7,8,9，e₁至e₉分别为{(0,0)，(1,0)，(0,1)，(-1,0)，(0，-1)，(1,1)，(-1,1)，(-1，-1)，(1，-1)}(表示一个节点附近的流体微团的九个离散粒子)；基于上述第二步中确定的灌浆管道内特定孔深位置的浆液压力分布函数进行求和，得到特定灌浆部位的浆液宏观压力和浆液宏观速度；

第四步，以注浆管道最下端的浆液宏观压力代表灌浆工艺中的灌浆压力，根据动态测量的灌浆压力控制灌浆过程中浆液注入流量大小，使得灌浆压力小于灌浆压力设计值；

第五步：当注浆过程中流量的变化量达到10％或者浆液密度发生变化时，则重复上述第一步至第四步，重新获取注浆管道的孔内灌浆压力。

上述的孔内灌浆压力的动态测量方法，优选的，所述构建浆液压力模型的第二步的具体操作包括以下步骤：

1)构建坐标系：测定注浆管道的直径为d，注浆管道伸入注浆孔内的长度为H，并以注浆管道的左下角为坐标原点建立坐标系，x方向为水平方向，y方向为重力反方向；离线确定x方向的步长Δx和y方向的步长Δy,并将坐标系构成的二维平面进行格子划分，x方向最大格子数记为N_X＝d/Δx，y方向最大格子数记为Ny＝H/Δy；

2)初始化：将整个网格N_X×Ny格子作为计算域，所有格子横竖相交的点设为节点，在整个格子场内将浆液速度、浆液密度和压力分布函数进行初始化；

2.1)浆液速度初始化，浆液在整个格子场内部x,y两个方向的初始速度场都是零,浆液流体流向为y的反方向；在注浆管道的最上端格子上浆液的初始速度为-U_y0；

2.2)浆液密度的初始化，在整个格子场内浆液密度ρ＝ρ₀；

2.3)压力分布函数的初始化，设定各格子内的压力分布函数f_k初始值等于各自的压力均匀分布函数所述压力均匀分布函数通过引入不可压缩流体模型(含义是采用二维平面模型，每个格子内采用9个粒子)来确立， $f_k^{\mathrm{eq}}=m_k{C}_S^2[\mathrm{\ρ}+{\mathrm{\ρ}}_0(1+3.0*e_k*u+4.5{(e_k*u)}^2-1.5u)],$ 其中m_k为权重系数、m₁至m₉分别为[4/9，1/9，1/9，1/9，1/9，1/36，1/36，1/36，1/36],u为格子内浆液速度，u为二维矢量；

3)注浆管道内边界格子的处理：

3.1)左边界处理：左边界格子的压力分布函数值＝本格子的压力均匀分布函数值+右边相邻格子的压力分布函数值-右边相邻格子的压力均匀分布函数值；

3.2)右边界处理：右边界格子的压力分布函数值＝本格子的压力均匀分布函数值+左边相邻格子的压力分布函数值–左边相邻格子的压力均匀分布函数值；

3.3)上边界处理：上边界格子的速度＝-U_y0，上边界格子的压力分布函数值＝本格子的压力均匀分布函数值+下边相邻格子的压力分布函数值–下边相邻格子的压力均匀分布函数值，

3.4)下边界处理：下边界格子的速度值＝上边相邻格子的速度值；下边界格子的压力分布函数值＝本格子的压力均匀分布函数值+上边相邻格子的压力分布函数值–上边相邻格子的压力均匀分布函数值；

4)确定压力分布函数演化过程：计算域内不同空间位置处每个格子内的压力分布函数f_k通过压力分布函数的演化方程来确定，压力分布函数的演化方程为 $f_k(x+\mathrm{\Δx},t+\mathrm{\Δt})=f_k(x,t)+(f_k^{\mathrm{eq}}(x,t)-f_k(x,t))/\mathrm{tau},$ 其中tau为压力分布函数的松弛因子，Δx为空间上的离散量，Δt为时间上的离散量，且Δx＝e_kΔt；在整个格子场开展压力分布函数演化方程的迭代运算，迭代次数为N_X×Ny×9，最终获取每个格子内的压力分布函数fk；

上述的孔内灌浆压力的动态测量方法，优选的，所述压力分布函数的演化方程主要通过引入不可压缩流体模型并基于以下步骤确立：由于时间和空间的离散并不是独立的，是通过粒子的离散速度联系起来，这表明粒子在两个时间步之间由一个节点迁移到相邻节点上，并在该节点上与其它粒子发生碰撞；利用该特性，我们将流体粒子的运动划分为碰撞和迁移两个过程；碰撞过程只发生在x本节点上，遵循公式 $f_k^{+}(x,t)=f_k(x,t)+(f_k^{\mathrm{eq}}(x,t)-f_k(x,t))/\mathrm{tau},$ 其重新分布了不同离散速度上的分布函数；不同离散速度方向上的粒子做迁移运动，运动到相邻节点上，遵循公式然后在该新的节点又发生碰撞，然后又迁移，进而确立得到压力分布函数的演化方程为 $f_k(x+\mathrm{\Δx},t+\mathrm{\Δt})=f_k(x,t)+(f_k^{\mathrm{eq}}(x,t)-f_k(x,t))/\mathrm{tau}.$

上述的孔内灌浆压力的动态测量方法，优选的，确定宏观压力和宏观速度的上述第三步的具体操作包括：基于上述最终获取的每个格子内的压力分布函数fk，再根据流体微团模型求出每个格子内浆液宏观压力值和浆液宏观速度值。

与现有技术相比，本发明的优点在于：本发明的注浆管道内孔内灌浆压力的动态测量方法能反映浆液在注浆管道内流动的动态变化本质，借助于现在的计算机数值计算技术很容易实现该方法的在线测量过程，能直接的测量孔内灌浆压力值；只要知道孔口的速度和浆液的密度值就能获取注浆管道内的任意处压力。当注浆管道内任意处的压力值达到或即将达到注浆管道的最大承受值时，现场施工人员则降低灌浆浆液的流量或者停止灌浆以避免事故的发生。

附图说明

图1为本发明实施例中灌浆装置的结构布置图。

图2为本发明实施例中注浆管道二维轴对称几何模型。

图3为本发明实施例中部分注浆管道内浆液速度分布云图。

图4为本发明实施例中沿着灌浆孔深度方向压力与孔深的关系图。

图例说明

1、灌浆液桶；2、高压泵；3、流量计；4、微机记录系统；5、阀门。

具体实施方式

为了便于理解本发明，下文将结合说明书附图和较佳的实施例对本发明作更全面、细致地描述，但本发明的保护范围并不限于以下具体的实施例。

除非另有定义，下文中所使用的所有专业术语与本领域技术人员通常理解的含义相同。本文中所使用的专业术语只是为了描述具体实施例的目的，并不是旨在限制本发明的保护范围。

实施例：

如图1至图4所示的一种孔内灌浆压力的动态测量方法，注浆管道为圆形钢管，注浆管道直径20毫米，如图1所示，高压泵2将灌浆液桶1内的浆液抽送到注浆管道的灌浆段处，在输送管道上安装有与微机记录系统4连接的流量计3和阀门5。本实施例中，注浆管道下端的压力代表孔内压力值；灌浆管道孔深50米，则注浆管道下端距离地表距离为50米处的压力值即为灌浆孔内压力值。

本实施例的孔内灌浆压力的动态测量方法，包括以下步骤：

第一步：收集灌浆过程相关技术参数，本实施例中，经过测量灌浆孔的长度为50米，注浆管道的直径为20毫米，浆液的初始密度ρ₀为1.8g/cm³，浆液流量为100L/min，根据圆形管道的流量与速度计算公式Q＝AU_y0，求出浆液初始速度为U_y0＝5.308m/s，即为上端浆液入口格子的y方向的初始速度。其中A为圆形管道横截面积。

第二步：选取模型变量，确定灌浆压力测量的模型结构：基于格子玻尔兹曼方法建立注浆管道内浆液的流体微团的压力分布函数模型(因为格子玻尔兹曼也是建立管道内浆液的流体压力的一种模型；这个模型要通过数值求解算法得到模型的数值解)；压力分布函数的确立具体包括以下步骤：

1)构建坐标系：以注浆管道的左下角为坐标原点建立坐标系，x方向为水平方向，y方向为重力反方向；离线确定x方向的步长Δx和y方向的步长Δy,并将坐标系构成的二维平面进行格子划分，x方向最大格子数记为N_X＝d/Δx，y方向最大格子数记为Ny＝H/Δy，d为注浆管道的直径为20毫米，H为注浆管道伸入注浆孔内的长度为50米；如图2所示，本实施例中对坐标系构成的二维平面进行格子划分，x方向每个格子的步长为1毫米，y方向每个格子的步长为50毫米，x方向的步长和y方向的步长为经验参数，根据灌浆过程密度和速度区间值反复离线计算确定。x方向最大格子数记为N_X＝20；y方向最大格子数记为Ny＝1000；整个计算域格子数则产生20×1000的格子数目。

2)初始化：将整个网格20×1000格子作为计算域，在整个格子场内将浆液速度、浆液密度和压力分布函数进行初始化；

2.1)浆液速度初始化，浆液在整个格子场内部x,y两个方向的初始速度场都是零,浆液流向为y的反方向；在注浆管道的最上端格子上浆液的初始速度的大小为U_y0＝5.308m/s；

2.2)浆液密度的初始化，在整个格子场内浆液密度ρ＝ρ₀＝1.8g/cm³；

2.3)压力分布函数的初始化，设定各格子内的压力分布函数f_k初始值等于各自的压力均匀分布函数所述压力均匀分布函数通过格子玻尔兹曼模型来确立， $f_k^{\mathrm{eq}}=m_k{\mathrm{\ρ}}_0[1+3.0*e_k*u+4.5{(e_k*u)}^2-1.5u],$ 其中m_k为权重系数、m₁至m₉分别为[4/9，1/9，1/9，1/9，1/9，1/36，1/36，1/36，1/36],u为格子内浆液速度，u为二维矢量；

3)注浆管道内边界格子的处理：

3.1)左边界处理：左边界格子的压力分布函数值＝本格子的压力均匀分布函数值+右边相邻格子的压力分布函数值-右边相邻格子的压力均匀分布函数值；

3.2)右边界处理：右边界格子的压力分布函数值＝本格子的压力均匀分布函数值+左边相邻格子的压力分布函数值–左边相邻格子的压力均匀分布函数值；

3.3)上边界处理：上边界格子的速度＝-5.308m/s，上边界格子的压力分布函数值＝本格子的压力均匀分布函数值+下边相邻格子的压力分布函数值–下边相邻格子的压力均匀分布函数值；

3.4)下边界处理：下边界格子的速度值＝上边相邻格子的速度值；下边界格子的压力分布函数值＝本格子的压力均匀分布函数值+上边相邻格子的压力分布函数值–上边相邻格子的压力均匀分布函数值；

4)宏观压力和宏观速度的计算，包括压力分布函数演化方程的计算和宏观压力和宏观速度的获取，压力分布函数的演化方程通过不可压浆液流体的格子玻尔兹曼模型来确立，压力分布函数的演化方程为 $f_k(x+\mathrm{\Δx},t+\mathrm{\Δt})=f_k(x,t)+(f_k^{\mathrm{eq}}(x,t)-f_k(x,t))/\mathrm{tau},$ 其中tau为压力分布函数的松弛因子；在整个格子场开展压力分布函数演化方程的迭代运算，迭代次数为20×1000×9，最终获取每个格子内的压力分布函数f_k。

第三步：再根据不可压缩流体微团的压力分布函数构建灌浆工艺中浆液宏观压力Pm与浆液速度u、浆液密度ρ的关系模型,即同样根据流体微团的压力分布函数构建灌浆工艺中浆液宏观速度Um与浆液速度u、浆液密度ρ的关系模型,即其中，k＝1,2,3,4,5,6,7,8,9，e₁至e₉分别为{(0,0)，(1,0)，(0,1)，(-1,0)，(0，-1)，(1,1)，(-1,1)，(-1，-1)，(1，-1)}；基于上述确定的灌浆管道内特定孔深位置的浆液压力分布函数进行求和，求出每个格子内浆液宏观压力值和浆液宏观速度值，如图3所示为本实施例中部分注浆管道内浆液速度分布云图。

第四步，以注浆管道最下端的浆液宏观压力代表灌浆工艺中的灌浆压力，本实施例中，沿着灌浆孔深度方向压力与孔深的关系如图4所示，在50米处的压力 根据获得的灌浆压力指导灌浆过程中浆液流量大小的控制，使得灌浆压力小于灌浆压力设计值。

第五步：当注浆过程中流量的变化量达到10％或者浆液密度发生变化，则重复第一步至第四步，重新获取注浆管道的灌浆压力。当注浆管道孔口流量变化为80L/min时，即U_y0为4.246m/s，注浆管道50米处的压力变化1.33KPa，变为5.246MPa。

本实施例中，将第一步至第五步以C语言的方式嵌入灌浆监控系统4中，根据实际的灌浆孔长度和监测流量和浆液密度，获取注浆管道内不同部位的压力值；同时节约场地空间。

本实施例的注浆管道内孔内灌浆压力的动态测量方法能反映浆液在灌浆管道内流动的动态变化本质，借助于现在的计算机数值计算技术很容易实现该方法的在线测量过程，能直接的测量孔内灌浆压力值；只要知道孔口的速度和浆液的密度值就能获取管道内的任意处压力。当注浆管道内任意处的压力值达到或即将达到注浆管道的最大承受值时，现场施工人员则降低灌浆浆液的流量或者停止灌浆以避免事故的发生。

Claims (4)

1.一种孔内灌浆压力的动态测量方法，其特征在于，包括以下步骤：

第一步：根据灌浆工艺，分别通过流量计和密度计获取灌浆过程中的流量值Q和密度值ρ₀，并由流量值Q计算得到浆液在注浆管道入口处的速度U_y0；

第二步：确定灌浆压力测量的模型结构：基于格子玻尔兹曼方法建立注浆管道内浆液的流体微团的压力分布函数模型；

第三步：再根据不可压缩流体微团的压力分布函数模型构建灌浆工艺中浆液宏观压力Pm与浆液速度u、浆液密度ρ的泛函关系模型,即其中，f_k(u,ρ)为压力分布函数；同样根据流体微团的压力分布函数构建灌浆工艺中浆液宏观速度Um与浆液速度u、浆液密度ρ的关系模型，即其中，e₁至e₉分别为{(0,0)，(1,0)，(0,1)，(-1,0)，(0，-1)，(1,1)，(-1,1)，(-1，-1)，(1，-1)}；基于上述第二步中确定的灌浆管道内特定孔深位置的浆液压力分布函数进行求和，得到特定灌浆部位的浆液宏观压力和浆液宏观速度；

第四步，以注浆管道最下端的浆液宏观压力代表灌浆工艺中的灌浆压力，根据动态测量的灌浆压力控制灌浆过程中浆液注入流量大小，使得灌浆压力小于灌浆压力设计值；

第五步：当注浆过程中流量的变化量达到10％或者浆液密度发生变化时，则重复上述第一步至第四步，重新获取注浆管道的孔内灌浆压力。

2.根据权利要求1所述的孔内灌浆压力的动态测量方法，其特征在于，上述第二步的具体操作包括以下步骤：

1)构建坐标系：测定注浆管道的直径为d，注浆管道伸入注浆孔内的长度为H，并以注浆管道的左下角为坐标原点建立坐标系，x方向为水平方向，y方向为重力反方向；离线确定x方向的步长Δx和y方向的步长Δy,并将坐标系构成的二维平面进行格子划分，x方向最大格子数记为N_X＝d/Δx，y方向最大格子数记为Ny＝H/Δy；

2)初始化：将整个网格N_X×Ny格子作为计算域，所有格子横竖相交的点设为节点，在整个格子场内将浆液速度、浆液密度和压力分布函数进行初始化；

2.1)浆液速度初始化，浆液在整个格子场内部x,y两个方向的初始速度场都是零,浆液流体流向为y的反方向；在注浆管道的最上端格子上浆液的初始速度设为-U_y0；

2.2)浆液密度的初始化，在整个格子场内浆液密度ρ＝ρ₀；

2.3)压力分布函数的初始化，设定各格子内的压力分布函数f_k初始值等于各自的压力均匀分布函数所述压力均匀分布函数通过引入不可压缩流体模型来确立， $f_k^{\mathrm{eq}}=m_k{C}_S^2[\mathrm{\ρ}+{\mathrm{\ρ}}_0(1+3.0*e_k*u+4.5{(e_k*u)}^2-1.5u)],$ 其中m_k为权重系数、m₁至m₉分别为[4/9，1/9，1/9，1/9，1/9，1/36，1/36，1/36，1/36],u为格子内浆液速度；

3)注浆管道内边界格子的处理：

3.1)左边界处理：左边界格子的压力分布函数值＝本格子的压力均匀分布函数值+右边相邻格子的压力分布函数值-右边相邻格子的压力均匀分布函数值；

3.2)右边界处理：右边界格子的压力分布函数值＝本格子的压力均匀分布函数值+左边相邻格子的压力分布函数值–左边相邻格子的压力均匀分布函数值；

3.3)上边界处理：上边界格子的速度＝-U_y0；上边界格子的压力分布函数值＝本格子的压力均匀分布函数值+下边相邻格子的压力分布函数值–下边相邻格子的压力均匀分布函数值；

3.4)下边界处理：下边界格子的速度值＝上边相邻格子的速度值；下边界格子的压力分布函数值＝本格子的压力均匀分布函数值+上边相邻格子的压力分布函数值–上边相邻格子的压力均匀分布函数值；

4)确定压力分布函数演化过程：计算域内不同空间位置处每个格子内的压力分布函数f_k通过压力分布函数的演化方程来确定，压力分布函数的演化方程为 $f_k(x+\mathrm{\Δx},t+\mathrm{\Δt})=f_k(x,t)+(f_k^{\mathrm{eq}}(x,t)-f_k(x,t))/\mathrm{tau},$ 其中tau为压力分布函数的松弛因子，Δx为空间上的离散量，Δt为时间上的离散量，且Δx＝e_kΔt；在整个格子场开展压力分布函数演化方程的迭代运算，迭代次数为N_X×Ny×9，最终获取每个格子内的压力分布函数f_k。

3.根据权利要求2所述的孔内灌浆压力的动态测量方法，其特征在于，所述压力分布函数的演化方程主要通过引入不可压缩流体模型并基于以下步骤确立：将流体粒子的运动划分为碰撞和迁移两个过程；碰撞过程只发生在x本节点上，遵循公式其重新分布了不同离散速度上的分布函数；不同离散速度方向上的粒子做迁移运动，运动到相邻节点上，遵循公式然后在该新的节点又发生碰撞，然后又迁移，进而确立得到压力分布函数的演化方程为 $f_k(x+\mathrm{\Δx},t+\mathrm{\Δt})=f_k(x,t)+(f_k^{\mathrm{eq}}(x,t)-f_k(x,t))/\mathrm{tau},$

4.根据权利要求1、2或3所述的孔内灌浆压力的动态测量方法，其特征在于，上述第三步的具体操作包括以下步骤：基于上述最终获取的每个格子内的压力分布函数f_k，再根据流体微团模型求出每个格子内浆液宏观压力值和浆液宏观速度值。