-
水在土壤中的流动是自然界中一种常见的渗流现象。达西定律常被用来描述宏观上统计均匀的土壤多孔介质中的单向低速流动, 在实际问题的求解中发挥了重要作用。但仅应用达西定律建立的渗流方程, 只能对流体的运动特征进行宏观表征, 无法在孔隙尺度(微米级)上精细地描述流体的运动特征, 具有研究尺度的限制。近年来, 人们逐渐意识到在孔隙尺度上探究多孔介质中的流体运动过程对于自然现象的认识和工程设备的设计有着极其重要的意义[1-5]。然而, 多孔介质包含复杂的孔隙结构, 在其内部的渗流、溶质运移过程比一般均匀介质中的相应过程要复杂[6]。流体在通过孔隙时, 流速、方向、流态随着流-固边界不断变化, 土柱试验难以详细探究[7-8], 即传统试验手段无法在孔隙尺度对土壤多孔介质进行研究。因此, 有学者提出采用模拟手段探究流体运动过程。一方面, 可将多孔介质抽象概化为简单的结构(孔网模型、毛细管模型等), 通过得到相近的解析解来获取等效的宏观参数, 但这种方式仍无法得到实际土壤结构的精确参数值[9-10]。另一方面, 可利用传统的数值模拟手段(有限元法、有限体积法、有限差分法等)进行微观尺度描述, 即采取离散非线性偏微分的Euler方程或Navier-Stocks方程等宏观方程, 再由计算机求解代数方程组的方式; 这种方式虽然可以从孔隙尺度进行流体运动的分析, 却在处理多孔介质的复杂边界时遇到困难, 模拟结果难以收敛[11-12]。
在此背景下, 格子Boltzmann法(Lattice Boltzmann Method, LBM)作为一种联系宏观与微观、具有介观粒子背景的数值模拟方法, 因易于处理边界、界面问题, 在近年逐渐发展起来[13]。它起源于传统物理学描述气体运动的连续方程Boltzmann方程和理论计算机领域的元胞自动机模型[14]。独特的动力学背景和计算机编码优势使得它在连续方程离散化、并行计算等方面有较好的表现。多项研究表明, 格子Boltzmann法在应对多孔介质孔隙流模拟中涉及的复杂边界有显著的优势[15-17]。Vogel等[11]比较了不同模型用于评估玻璃珠结构水力学特性的优劣, 发现LBM在流体运动实时输出上有更可信的结果。Chen等[18]利用LBM动态模拟研究了胶体沉积对多孔介质结构的影响, 并描述了相应流体运动和溶质运移的变化情况。Zhang等[5]用LBM模拟了分辨率为30 μm的土样结构中的溶质运移过程, 并与有限体积法和解析解的结果进行比较, 指出了在未详细研究孔隙尺度流体行为的情况下解释宏观试验结果是非常不准确的。基于分子动理论的Boltzmann方程在微观层面通过对结构内部的分子体系采样, 获得分子的详细轨道图景及系统的各物理量分布, 在给定条件下可以回溯宏观的Navier-Stokes方程, 是一种自底而上的建模方式, 它在计算中考虑了流体孔隙内的微观流动, 相比传统的体积平均方法有着物理本质刻画上的优越性, 同时由于空间和时间离散方式简单, 并行性优良, 可以有效减少庞大节点数量带来的计算耗时。目前已知的前人所做的相关研究主要集中于多孔介质中的流体流动模拟及并行, 较少考虑溶质运移。尤其受限于计算能力, 已有的研究大多考虑二维结构或简单三维结构中的流动和运移, 极少考虑真实土壤结构中的三维流动与溶质迁移。
本研究采用格子Boltzmann法在孔隙尺度分别对土壤等多孔介质进行水流运动和溶质运移的模拟, 从而得出渗透率、迂曲度、纵向弥散系数等宏观参数以描述多孔介质的水力学特性；同时利用GPU并行技术缩短计算时间。本研究有助于进行孔隙尺度土壤水力学特性快速评估, 以期对土壤改良中的结构优化提供新的机理认识。
-
本研究采用三维格子Boltzmann法的D3Q19模型(三维空间、19个速度方向, 如图 1所示), 用于模拟孔隙流流动和溶质运移过程, 模型的构建基础为单弛豫时间的BGK模型, 其采用Bhatnagar、Gross和Krook提出的BGK近似, 用一个简单算子代替Boltzmann方程中的碰撞项[19]。相比于D3Q13、D3Q15、D3Q17等三维格子Boltzmann模型[20], 本研究采用的D3Q19模型可以更精准地描绘复杂结构中的流体运动, 其演化方程包括水流模拟和溶质运移两部分。
图 1 格子Boltzmann法的D3Q19模型
Figure 1. The D3Q19 lattice structure in the LBM
-
水流流场模拟的演化方程[14]为
$$ \underbrace {{f_i}\left( {\mathit{\boldsymbol{x}} + {\mathit{\boldsymbol{e}}_i}\Delta t, t + \Delta t} \right) - {f_i}\left( {\mathit{\boldsymbol{x}}, t} \right)}_{迁移} = \underbrace { - \frac{1}{\tau }\left[ {{f_i}\left( {\mathit{\boldsymbol{x}}, t} \right) - f_i^{{\rm{eq}}}\left( {\rho , \mathit{\boldsymbol{u}}} \right)} \right]}_{碰撞}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {i = 1, 0, 1, 2, \cdots , 18} \right) $$ (1) 包括碰撞和迁移两部分, 其中fi表示i方向的分布函数, 平衡态分布函数feq则由式(2)表示:
$$ f_i^{{\rm{eq}}}\left( {\rho , \mathit{\boldsymbol{u}}} \right) = {w_i}\rho \left( \mathit{\boldsymbol{x}} \right)\left[ {1 + \frac{{3{\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{u}}}}{{{c^2}}} + \frac{{{9}{{\left( {{\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{u}}} \right)}^2}}}{{2{c^4}}} - \frac{{3{\mathit{\boldsymbol{u}}^2}}}{{2{c^2}}}} \right] $$ (2) 式中: i代表方向:当i=0时, ei=(0, 0, 0)c；当i=1, 2时, ei=(±1, 0, 0)c；当i=3, 4时, ei= (0, ±1, 0)c；当i=5, 6时, ei=(0, 0, ±1)c；当i=7, …, 10时, ei=(±1, ±1, 0)c；当i=11, …, 14时, ei=(0, ±1, ±1)c；当i=15, …, 18时, ei=(±1, 0, ±1)c。c为单位格子长度与时间步的比值; τ为量纲一的松弛时间; wi为权重系数, 对于D3Q19模型来说, 在i=0时, wi=1/3, 对于i=1, …, 6, wi=1/18, 而当i=7, …, 18时, wi=1/36。宏观流体密度ρ和速度u可以由式(3)和式(4)得到:
$$ \rho = \sum\limits_{i = 0}^{18} {{f_i}} $$ (3) $$ \mathit{\boldsymbol{u = }}\frac{1}{\rho }\sum\limits_{i = 0}^{18} {{f_i}{\mathit{\boldsymbol{e}}_i}} $$ (4) 该模型能回溯纳维-斯托克斯方程(Navier-Stokes equations)。
系统边界采用了在进出口处压力差恒定的周期边界[19]。对于不可压缩流体来说, 流体质点的密度ρ在运动过程中保持不变, 某一方向上相邻节点的压差Δp可以通过p=cs2ρ计算得到, cs为音速。在D3Q19模型中, cs2=c2/3=1/3, 则p=ρ/3[26-27]。为使宏观流体符合达西定律的应用前提, 确保雷诺数(Re)远小于1, 同时流体密度变异小, 即马赫数(Ma)处在一个较小的值(说明:一般认为Ma≤0.3的流体是不可压缩的[28])。本文中给定的压差是一个较小的值, 流体流动符合达西定律, 同时在进口和出口两侧增加一定数量的格子作为缓冲区。
通过以上水流模拟得到的参数可进一步得出水力学特性参数渗透率及迂曲度。渗透率的大小表达了流体通过多孔介质的能力, 是工程应用中非常重要的参数。当通过多孔介质流体的流速较低或者雷诺数较小(< 1)时, 经典的达西定律被广泛用来描述多孔介质内的流体流动问题, 并可以通过该定律来计算多孔介质的渗透率。本研究中的渗透率可以通过下式计算得到:
$$ {k_j} = \frac{{\overline {\rho q\upsilon } L}}{{\Delta P}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {j = x, y, z} \right) $$ (5) 式中: ρq为通过样品的平均流量；截面流量q为通过平均孔隙流速u得到；υ为运动黏度；L为样品尺度；ΔP为进出口之间的压差。
迂曲度(Tj)表示流体流经的实际轨迹长度与直线距离的比值, 在一定程度上可以表征多孔介质孔隙结构的复杂性。但是, 这个比值难以直接测得, 通常通过图像分析和溶质运移模拟间接得到。本研究中的迂曲度采用了Muljadi等[24]提出的方法, 通过孔隙流的流速计算得到:
$$ {T_j} = \frac{{\mathit{\boldsymbol{\overline u}} }}{{\overline {{u_j}} }} $$ (6) 式中: u代表流体的平均流速；uj为j方向上的速度分量, 在这项研究中, j为x、y或z方向。
-
溶质运移的分布函数用gi定义如下:
$$ {g_i}\left( {\mathit{\boldsymbol{x}} + {\mathit{\boldsymbol{e}}_i}\Delta t, t + \Delta t} \right) = {g_i}\left( {\mathit{\boldsymbol{x}}, t} \right) - \frac{{{g_i}\left( {\mathit{\boldsymbol{x}}, t} \right) - g_i^{{\rm{eq}}}\left( {C, \mathit{\boldsymbol{u}}} \right)}}{{{\tau _{\rm{s}}}}}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {i = 1, 0, 1, 2, \cdots , 18} \right) $$ (7) 与流场模拟类似, τs为溶质羽的量纲一弛豫时间, gieq(C, u)为x处、t时刻、i方向上的平衡态, 可以用下式表征:
$$ g_i^{{\rm{eq}}}\left( {C, \mathit{\boldsymbol{u}}} \right) = {w_i}C\left[ {1 + \frac{{3{\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{u}}}}{{{c^2}}} + \frac{{9{{\left( {{\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{u}}} \right)}^2}}}{{2{c^4}}} - \frac{{3{\mathit{\boldsymbol{u}}^2}}}{{2{c^2}}}} \right] $$ (8) 式中: C为宏观溶质浓度, 计算公式为
$$ C = \sum\limits_{i = 0}^{18} {{g_i}} $$ (9) 通过式(9)可回溯为对流弥散方程:
$$ \frac{{\partial C}}{{\partial t}} + \left( {\mathit{\boldsymbol{u}} \cdot \nabla } \right)C = \nabla \cdot \left( {D{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \nabla C} \right) $$ (10) 其中格子单位下的分子扩散系数D可以通过D=(τs－0.5)Δx2/3Δt计算得到。对于溶质场的边界, 为在z方向上的进出口面采用∂C/∂z=0, 其他4个面及内部边界则采用无流量格式。
通过以上溶质运移模拟得到的溶质场可得出多孔介质中溶质运移的弥散系数及穿透曲线。本文的研究在稳态流场的基础上, 进行了瞬时面源溶质运移的模拟, 主要的评估参数纵向弥散系数(DL)的计算由矩估计得到:
$$ {D_{\rm{L}}} = {D_z} = {\rm{d}}\sigma _z^2\left( t \right)/2{\rm{d}}t $$ (11) 式中: σz2(t)为污染羽在z方向上的空间二阶矩。为了直观表征分子扩散之外的因素对纵向弥散的影响, 本文采用纵向弥散系数(DL)与分子扩散系数(Dm)来归一评估溶质运移。
通过实时监测介质出口处的溶质浓度变化, 可绘制得到穿透曲线。为确保LBM代码的准确性, 在之前的工作中, 将其与多个明确定义的几何结构中的流动状态解析解进行对比, 并将数值模拟结果与简单圆球结构的渗透率土柱实验结果进行比较[25]。
-
从水流和溶质运移的计算过程可以看出, 对某三维格子n而言, 其碰撞和迁移计算所需的数据仅和与之相邻的18个格子有关, 并且各个格子之间的计算相互独立。因此, 基于D3Q19模型的水流和溶质运移模拟具有良好的可并行性。与CPU相比, GPU具有极高的加速比、擅长单指令多线程并发、更新换代速度快、天然的并行架构等优势[26]。本文基于GPU架构实现水流和溶质运移模拟的并行。
-
本文使用了3种具有不同空间异质性的多孔介质结构, 分别是粒径一致规则排列的圆球结构、粒径不等随机排列的圆球结构(3种结构的基本信息见表 1)以及同步辐射X射线显微CT扫描的土样(图 2)。为验证方法的有效性, 通过人工构建的简单结构做模拟的探索, 进而应用到真实的多孔介质中, 用于测试的结构异质性亦逐步增加。
表 1 多孔介质的基本信息
Table 1. Physical properties of porous structure
样品序号 计算区域 真实大小/mm3 总孔隙度/% Sphere300 300×300×300 1.11×1.11×1.11 47.61 Sphere160 160×160×160 0.592×0.592×0.592 47.89 Isphere160 160×160×160 0.592×0.592×0.592 16.48 Soil400 400×400×400 1.48×1.48×1.48 31.12 
图 2 多孔介质结构示意
Figure 2. Three-dimensional pore geometry of glass beads
(1) 粒径一致规则排列的圆球结构, 包括5行5列5排粒径均为32像素的圆球堆叠(Sphere160)以及10行10列10排粒径均为30像素的圆球堆叠(Sphere300)(图 2(a))。
(2) 粒径不等随机排列的圆球结构(Isphere160)则由粒径范围为6~32像素的圆球随机排列生成, 为了和真实土壤结构更为接近, 圆球球体允许重叠, 置放位置完全随机, 先大球后小球, 为了增加固相占比, 允许圆球超出边界, 最后再修剪(图 2(b))。
(3) 土壤样品(Soil400)是孔隙度为31.12%的改良黏土团聚体, 其CT扫描在上海同步辐射光源BL13W光束线站进行, 各角度的扫描投影图由像素分辨率为3.7 μm的CCD探测器记录(图 2(c))。CT图在去除了伪影和进行了图像分割之后, 进行三维重构, 最终截取了400×400×400像素大小的三维图像作为计算区域(对应的真实大小为1.48 mm×1.48 mm×1.48 mm)。为了便于不同样品之间的对比, 构造的多孔介质分辨率也设为3.7 μm(一个像素为一个节点, 分辨率即对应节点的边长)。
在完成多孔介质图像的三维重构之后, 多孔介质的结构信息转成二进制文件供Boltzmann模拟器使用, 用于定义流-固边界, 随后分别进行流场和溶质场的模拟。利用Thorne[27]的研究成果, 本文设定流体和溶质的松弛时间为1, 其中流体的松弛时间$\tau = \frac{{6\upsilon + 1}}{2}, \upsilon = \frac{1}{6}; $溶质的松弛时间$ {\tau _{\rm{s}}} = \frac{{6D + 1}}{2}, D = \frac{1}{6}$。流场的模拟通过在样品某个方向(x、y或z方向)的进口处和出口处施加恒定的压差ΔP, 边界处理方式为周期性边界。而其他4个面及结构内部的流-固交界面则采用无滑移反弹边界。ΔP根据雷诺数的大小来确定(ΔP=0.05), 以确保宏观流场符合达西定律, 模拟过程中雷诺数维持在1以内。两侧各增加10个格子的缓冲区, 以确保入口和出口处的流体平滑。
在经历了一定时间步的流体粒子碰撞和迁移过程之后, 本项研究在得到稳态流场的基础上, 叠加溶质场模拟的模块, 在样品中央注入一个瞬时面源溶质, 并在污染羽超出边界之前进行纵向弥散系数的计算, 稳态流场的标准是当两个时间步之间的平均流速符合|[u (t)－ u (t－Δt)]/ u (t)| ≤10－4, 选择样品中心线位置注射是为了确保统计弥散系数的时刻溶质羽没有溢出边界。最后通过矩估计得到纵向弥散系数。
本研究的整个计算流程均实现了每个时间步内的并行, 由于t+1时刻的计算依赖于t时刻的结果, 因此对每个时间步内各个格子的计算进行并行。在完成数据分区、计算数组建立、计算线程分配之后, 进入对应的时间步计算各个格子新的分布函数, 进行各计算节点之间的通信和同步。
在GPU架构上实现本文算法的并行为步骤(1)—(3)；同时, 为扩大计算尺度, 采用MPI的多GPU并行需增加步骤(4)—(6)。
(1) 数据存储由于GPU显存和内存之间的数据传输较为耗时, 为尽可能提高计算速度, 本文将计算数据(流体速度、密度、分布概率函数等)全部存储于GPU显存。计算过程均在GPU内完成, CPU负责控制计算的整体流程, 仅在提取结果时将数据拷贝至内存。
(2) 线程分配为充分利用GPU的高度并行性, 本文对每一个格子均分配一个计算线程。每个时间步的循环中, 计算是按线程块轮流执行的, 当一个线程块正在执行数据读取命令时, GPU将自动切换到其他线程块进行计算, 从而掩盖数据读取的延迟。
(3) 并行计算在每个线程中, 由于寄存器的读取速度远快于显存的读取速度, 程序将从GPU显存中读取计算所需数据到每个线程的寄存器中, 进一步加快计算速度。
(4) 数据分区将尺寸为m×n×h的样品在h方向根据MPI计算节点数量等分为多个区域, 计算节点间通过InfiniBand等高速低延迟的通信链路进行数据传输。
(5) 各GPU内部的计算各个GPU对分配到的计算任务同步计算, 实现GPU之间的并行。
(6) 计算节点间的通信与同步完成各GPU内部的计算之后, 对相邻分区共用格子的数据进行通信和同步, 从而保证计算结果的正确性。
-
多孔介质内部的水流场分布十分复杂, 受孔隙连通情况、孔隙形态、颗粒分布和大小的影响。如图 3所示, 规则排列的圆球结构流线与有恒定压差的z方向大致保持平行, 流速均匀；而真实土样的流线则曲折多变, 分布随机, 在部分孔隙处呈现较快的流速(图 3(b)由于真实土样的三维固相过于复杂, 故只显示流线)。
图 3 样品流线模拟
Figure 3. Simulated streamlines of water flow
-
渗透率是多孔介质地下水流系统中最为重要的水文地质参数[28]。表 2显示规则堆放圆球结构(Sphere300及Sphere160)的渗透率远大于随机排列圆球结构(Isphere160)和真实土壤样品(Soil400), 考虑到Sphere300及Sphere160的孔隙率较高(表 1), 说明多孔介质孔隙率的增加对渗透率的提升有一定的积极作用。相对Isphere160样品, Soil400孔隙率高, 其渗透率反而更低, 表明结构的空间异质性(如连通性、迂曲度、各向异性)对多孔介质的渗透率会有明显影响。另外, Sphere300及Sphere160为规则堆放结构, 因此, 它们沿坐标轴3个方向的渗透率值各自相同, 表现出各向同性, 且该性质与样品尺寸无关。随机排列圆球结构(Isphere160)和真实土样(Soil400)结构均展现出一定的各向异性[29], 同一样品3个方向的渗透率值不同。
表 2 多孔介质的水力学特性
Table 2. Hydraulic properties of porous structure
样品序号 kx/mD ky/mD kz/mD Tx Ty Tz DL/Dm Sphere300 28 772 28 772 28 772 1.02 1.02 1.02 0.696 8 Sphere160 25 188 25 188 25 188 1.03 1.03 1.03 0.656 1 Isphere160 261 782 936 1.66 1.45 1.41 0.397 8 Soil400 696 854 604 1.46 1.40 1.59 0.234 3 注:渗透率单位mD为毫达西, 1mD=10-15 m2。 -
通过格子Boltzmann法进行的流体运动模拟, 分别得到了4个样品在不同方向上的水流流速, 并对其进行了迂曲度计算。从计算结果(表 2)可以发现规则排列的圆球结构(Sphere300及Sphere160)在3个方向上的迂曲度相同, 而其余两个样品的迂曲度呈现各向异性, 这与结构的各向异性一致。从整体来看, 渗透率大的样品, 迂曲度较小。
-
在相等压差下的(DL/Dm)计算结果详见表 2, (DL/Dm)对应z方向(即施加压力梯度的方向), 可以看出, 在压力梯度一定的前提下, 渗透率较大的样品, 纵向弥散系数也较大。
在均质的多孔介质中, 弥散系数是一个恒定的值。然而在对真实土壤样品的测试中, 发现随着时间的推进, 该土壤介质中的纵向弥散系数与迁移时间呈线性递增的关系, 也就是通常说的反常运移[30]。在土壤等复杂结构中, 对于溶质运移具体机制的研究, 经常使用到一个参量——佩克莱数(Péclet number, 简称Pe数)。佩克莱数是一个量纲一数。其物理意义为对流速率与扩散速率之比, 其中扩散速率是指在一定浓度梯度驱使下的分子扩散速率。其定义式为Pe=uLc/D, 其中u为流速, Lc为特征长度, D为分子扩散系数, 均采用格子单位。对真实土样在不同Pe数下的纵向弥散系数的时间演变进行了刻画, 发现随着Pe数的增加, 曲线斜率更加陡峭, 虽然由于样品尺寸的限制, 无法进行更多时间步的测试, 但是根据图 4的趋势可以预计, 越大Pe数的样品, 达到稳定值所需的时间步越多。同时表 2中(DL/Dm)的计算基于相等压差的条件下, Pe数在10-2水平上。虽然规则排列的两个圆球样品(Sphere300、Sphere160)的平均流速高于Isphere160和Soil400两个数量级, 但在Pe数非常小的情况下, 溶质运移主要由分子扩散主导, 对流对弥散系数的影响有限, 因此两者的纵向弥散系数没有出现量级差异。
图 4 真实土样的纵向弥散系数的时间演变
Figure 4. Temporal evolution of the longitudinal dispersion coefficient, DL/Dm of Soil400
-
一般认为在溶质运移过程中观察到的反常弥散是不同尺度下的介质异质性导致, 而均匀介质中的溶质运移则符合菲克定律[31-32]。为了进一步探究溶质运移的特性, 在溶质运移的模拟过程中, 实时监测了出口处的溶质浓度变化, 并绘制了渗出液的浓度穿透曲线。由于Soil400样品的渗透性很弱, 溶质场中的溶质羽历经每一万个时间步大约在z方向上迁移90格, 得到较为完整的穿透曲线需要的时间步数超过80万步。图 5展示了计算域尺度均为160的Sphere160、Isphere160和尺度为400的Soil400的穿透曲线, 分别展示为红线(粒径一致规则排列的圆球结构), 蓝线(粒径不等随机排列的圆球结构)和黑线(真实土样)。溶质羽(Cout为出口处浓度, Co为初始浓度)在流场稳定后的t*=50 000时刻注入(流场达到稳定大概在t*=20 000~30 000), 可以发现规则排列的Sphere160在溶质羽注入后渗出液迅速达到一个峰值, 而后骤降, 渗出液浓度峰型集中, 而结构异质性较强的Isphere160的渗出液浓度上升较缓慢, 完成穿透曲线需要更多的时间步, 且渗出液达到波峰后的拖尾较长, 展现出较为强烈的非高斯性质, 即反常运移[34]。由于异质性的增加, 相较于Sphere160样品, Isphere160的曲线波峰右移对应峰值较低, 达到波峰后下降速度减缓, 这表明溶质在结构中迁移速度减缓, 且同时间内残留物质更多。与此相对, Sphere160的穿透曲线上升下降过程更快, 由于对流作用是影响该结构中溶质运移的主要因素, 流场中的水流速度大, 使溶质运移加快。由于本文不考虑溶质的化学反应和吸附, 认为Isphere160的拖尾现象是介质孔隙形态、颗粒随机分布的不均匀性造成的, 孔隙间的连通性较差以及狭小的孔隙通道使得在溶质运移后期仍有较大残留。对于真实样品(图 5(b), 起始时间从注射溶质开始), 在该样品内观察到了较长的拖尾, 其空间异质性远高于圆球样品, 非高斯现象较之圆球样品更为明显。
图 5 样品穿透曲线
Figure 5. Breakthrough curve of Sphere160, Isphere160 and Soil400
本文使用了浙江大学计算机辅助设计与图形学(CAD & CG)国家重点实验室的GPU集群, GPU卡为NIVIDIA Tesla K20M(频率500 MHz, 核心数2 496个)。每个主机上有4个GPU卡, 采用MPI实现联机多GPU并行。表 3进行了GPU并行加速效果的展示, 以真实土壤样品(Soil400)为例, Soil400利用双CPU工作站进行5万步的流场模拟耗时30 h, 而使用GPU集群后, 5万时间步的耗时缩短至55 min。Isphere160在CPU计算5万个时间步需28 h, 实施GPU并行之后在单显卡上耗时7 min, 在GPU集群上仅需4 min。大小相同、结构不同多孔介质样品模拟得到相对完整的穿透曲线耗时不同, 短的需要5万个时间步, 长的则需要超过80万个时间步, 说明计算耗时不仅与结构大小有关, 更与结构中孔隙所占的节点数量以及孔隙的连通性(图 6)有关。
表 3 GPU并行的加速效果
Table 3. Effect of GPU parallel computation
样品序号 孔隙节点数 CPU时间(tCPU)/h GPU时间(tGPU)/min tCPU/tGPU Sphere160 1 961 750 50.70 9.57 318 Isphere160 675 041 27.96 6.67 252 Soil400 19 914 848 29.59 56.05 32 注:Sphere160和Isphere160计算条件为CPU(i5-7300, 2.50 GHz)和计算显卡NVIDIA GeForce GTX TITAN X(GM200, 1 000 MHz)；Soil 400计算条件为双CPU(Intel Xeon E5-2680 2.8 GHz)和GPU集群, 4张NIVIDIA Tesla K20M(频率500 MHz, 核心数2 496个)。 
图 6 样品孔隙骨骼化示意
Figure 6. Skeletons of Sphere160 and Isphere160
-
(1) 固有渗透率与介质的空间异质性相关, 同时, 本研究发展的方法能准确量化出多孔结构非均匀性造成的固有渗透率的各向异性。
(2) 计算结果表明渗透率大的样品, 迂曲度较小；在压力梯度一定的前提下, 纵向弥散系数较大。
(3) 通过对溶质运移穿透曲线的研究, 发现即使在结构较为简单的多孔介质中, 也会出现非高斯对称形式的反常运移。
(4) GPU并行显著减少了计算时间, 并将前人200左右的研究尺度提升到400。
此方法为三维饱和格子Boltzmann模型, 可用于土壤中流体运动及溶质运移过程的研究, 有助于加深对于土壤结构改良效应的机理认识。此外, 也将进一步开发非饱和模型, 应用于污染物运移、微生物运动、植物养分吸收等方向的研究。
Pore-scale simulations of fluid flow and solute transport in porous media by high-performance Lattice Boltzmann Method
-
摘要: 深入探究孔隙尺度下的流体流动特性和溶质运移规律对石油开采、农田养分管理、地下水污染修复有着重要意义。以人工构建的多孔介质结构和同步辐射X射线显微CT扫描的土壤团聚体（分辨率3.7 μm）为研究对象，在空间节点数多达64 000 000的情况下，基于格子Boltzmann模型和GPU并行技术计算得到多孔介质流体运动和溶质运移过程的关键参数，并据此探究多孔介质空间异质性对水力学特性的影响。通过对3组不同结构的多孔介质比较发现，结构复杂程度最高的土壤样品和不规则堆叠的圆球结构的渗透率在100 mD（即10-13m2）量级，远低于规则堆叠的圆球结构（>20 000 mD）；土壤的迂曲度为1.40~1.60，明显高于规则堆叠的圆球结构。研究结果表明，渗透率大的样品具有较小的迂曲度，这与结构的空间异质性有较强的关系；土壤的渗透率和迂曲度呈现各向异性；在水力梯度一定的前提下，渗透率较大的样品，纵向弥散系数也较大；同时，结构的异质性也会影响溶质的穿透曲线。本研究提出的模拟方法可在土壤结构中进行高效的水流运动和溶质运移模拟，可用于土壤多孔介质在孔隙尺度下的水力学特性研究。
-
关键词:
- 孔隙尺度 /
- 格子Boltzmann法 /
- 流体流动 /
- 溶质运移
Abstract: Understanding the mechanism of fluid flow and solute transport at the pore scale is of great importance for oil recovery, crop nutrient management and groundwater pollution restoration. This study employed lattice Boltzmann model combining with GPU parallel technology to investigate the porous media of computer-generated structures and synchrotron-based X-ray micro-CT scans of soil aggregates (resolution 3.7 μm). The key parameters of fluid flow and solute transport in the porous media were obtained, and the influence of spatial heterogeneity of porous media on hydraulic properties was explored by high-performance simulation (spatial nodes up to 64 000 000). By comparing the three groups of porous media with different structures, it was found that the permeabilities of the soil sample with the highest structural complexity and beads irregularly stacked are on the order of 100 mD (i.e. 10-13m2), which is much lower than that of the regularly stacked beads (>20 000 mD); The soil sample has a tortuosity of 1.40~1.60, which is significantly higher than that of the regularly stacked beads. Our results show that the porous media with high permeabilities have small degree of tortuosity, indicating that the permeabilities of porous media are related to the spatial heterogeneity of the structure. The permeability and tortuosity of soil aggregate are anisotropic. At given pressure gradient, the longitudinal diffusion coefficient is greater for a sample with higher permeability. The heterogeneity of the pore structure also affects the breakthrough curve. The method established in this work can simulate water flow and solute migration in real soil structure, and can be used to study the hydraulic characteristics of porous media at the pore scale.-
Key words:
- pore-scale modeling /
- Lattice Boltzmann Method /
- fluid flow /
- solute transport
-
图 2 多孔介质结构示意
(a)规则排列的圆球结构(b)随机排列的圆球结构(c)土样扫描
(左列三维重构图, 右列二维扫描图)Figure 2. Three-dimensional pore geometry of glass beads
图 4 真实土样的纵向弥散系数的时间演变
Figure 4. Temporal evolution of the longitudinal dispersion coefficient, DL/Dm of Soil400
图 5 样品穿透曲线
(a)规则圆球结构、随机排列圆球结构和土壤样品(b)土壤样品(Soil400)
Figure 5. Breakthrough curve of Sphere160, Isphere160 and Soil400
表 1 多孔介质的基本信息
Table 1. Physical properties of porous structure
样品序号 计算区域 真实大小/mm3 总孔隙度/% Sphere300 300×300×300 1.11×1.11×1.11 47.61 Sphere160 160×160×160 0.592×0.592×0.592 47.89 Isphere160 160×160×160 0.592×0.592×0.592 16.48 Soil400 400×400×400 1.48×1.48×1.48 31.12 
下载: 导出CSV
表 2 多孔介质的水力学特性
Table 2. Hydraulic properties of porous structure
样品序号 kx/mD ky/mD kz/mD Tx Ty Tz DL/Dm Sphere300 28 772 28 772 28 772 1.02 1.02 1.02 0.696 8 Sphere160 25 188 25 188 25 188 1.03 1.03 1.03 0.656 1 Isphere160 261 782 936 1.66 1.45 1.41 0.397 8 Soil400 696 854 604 1.46 1.40 1.59 0.234 3 注:渗透率单位mD为毫达西, 1mD=10-15 m2。
下载: 导出CSV
表 3 GPU并行的加速效果
Table 3. Effect of GPU parallel computation
样品序号 孔隙节点数 CPU时间(tCPU)/h GPU时间(tGPU)/min tCPU/tGPU Sphere160 1 961 750 50.70 9.57 318 Isphere160 675 041 27.96 6.67 252 Soil400 19 914 848 29.59 56.05 32 注:Sphere160和Isphere160计算条件为CPU(i5-7300, 2.50 GHz)和计算显卡NVIDIA GeForce GTX TITAN X(GM200, 1 000 MHz)；Soil 400计算条件为双CPU(Intel Xeon E5-2680 2.8 GHz)和GPU集群, 4张NIVIDIA Tesla K20M(频率500 MHz, 核心数2 496个)。
下载: 导出CSV
-
[1] NING Y, JIANG Y, LIU H L, et al. Numerical modeling of slippage and adsorption effects on gas transport in shale formations using the lattice Boltzmann method[J]. Journal of Natural Gas Science and Engineering, 2015, 26: 345-355. doi: 10.1016/j.jngse.2015.06.015 [2] 姚军, 赵建林, 张敏, 等.基于格子Boltzmann方法的页岩气微观流动模拟[J].石油学报, 2015, 36(10): 1280-1289. http://d.old.wanfangdata.com.cn/Periodical/syxb201510011 YAO J, ZHAO J L, ZHANG M, et al. Microscale shale gas flow simulation based on Lattice Boltzmann Method[J]. Acta Petrolei Sinica, 2015, 36(10): 1280-1289. (in Chinese) http://d.old.wanfangdata.com.cn/Periodical/syxb201510011 [3] XU H, DANG Z. Lattice Boltzmann modeling of carbon deposition in porous anode of a solid oxide fuel cell with internal reforming[J]. Applied Energy, 2016, 178: 294-307. doi: 10.1016/j.apenergy.2016.06.007 [4] KANG Q J, CHEN L, VALOCCHI A J, et al. Pore-scale study of dissolution-induced changes in permeability and porosity of porous media[J]. Journal of Hydrology, 2014, 517: 1049-1055. doi: 10.1016/j.jhydrol.2014.06.045 [5] ZHANG X X, CRAWFORD J W, FLAVEL R J, et al. A multi-scale Lattice Boltzmann model for simulating solute transport in 3-D X-ray micro-tomography images of aggregated porous materials[J]. Journal of Hydrology, 2016, 541: 1020-1029. doi: 10.1016/j.jhydrol.2016.08.013 [6] 康学远, 施小清, 邓亚平, 等.基于EnKF融合地球物理数据刻画含水层非均质性[J].水科学进展, 2018, 29(1): 40-49. http://journal16.magtechjournal.com/Jweb_skxjz/CN/abstract/abstract2795.shtml KANG X Y, SHI X Q, DENG Y P, et al. Assimilation of hydrogeophysical data for the characterization of subsurface heterogeneity using Ensemble Kalman Filter (EnKF)[J]. Advances in Water Science, 2018, 29(1): 679-684. (in Chinese) http://journal16.magtechjournal.com/Jweb_skxjz/CN/abstract/abstract2795.shtml [7] 周昊, 芮淼, 岑可法.多孔介质内流体流动的大涡格子Boltzmann方法研究[J].浙江大学学报(工学版), 2012, 46(9): 1660-1665. doi: 10.3785/j.issn.1008-973X.2012.09.017 ZHOU H, RUI M, CEN K F. Study of flow in porous media by LES-LBM coupling method[J]. Journal of Zhejiang University (Engineering Science), 2012, 46(9): 1660-1665. (in Chinese) doi: 10.3785/j.issn.1008-973X.2012.09.017 [8] 蒋昌波, 刘洋, 邓斌, 等.出渗对均匀大粒径泥沙附近水动力特性影响[J].水科学进展, 2017, 28(1): 67-75. http://journal16.magtechjournal.com/Jweb_skxjz/CN/abstract/abstract2694.shtml JIANG C B, LIU Y, DENG B, et al. Effect of upward seepage on hydrodynamic characteristics around uniform coarse grains[J]. Advances in Water Science, 2017, 28(1): 67-75. (in Chinese) http://journal16.magtechjournal.com/Jweb_skxjz/CN/abstract/abstract2694.shtml [9] AHRENHOLZ B, TOLKE J, LEHMANN P, et al. Prediction of capillary hysteresis in a porous material using lattice-Boltzmann methods and comparison to experimental data and a morphological pore network model[J]. Advances in Water Resources, 2008, 31(9): 1151-1173. doi: 10.1016/j.advwatres.2008.03.009 [10] MEHMANI Y, BALHOFF M T. Mesoscale and hybrid models of fluid flow and solute transport[J]. Reviews in Mineralogy and Geochemistry, 2015, 80(1): 433-459. doi: 10.2138/rmg.2015.80.13 [11] VOGEL H J, TOLKE J, SCHULZ V P, et al. Comparison of a Lattice-Boltzmann model, a full-morphology model, and a pore network model for determining capillary pressure-saturation relationships[J]. Vadose Zone Journal, 2005, 4(2): 380-388. doi: 10.2136/vzj2004.0114 [12] 徐鹏, 李翠红, 柳海成, 等.多尺度多孔介质有效气体输运参数的分形特征[J].中国地质大学学报, 2017, 42(8): 1373-1378. http://d.old.wanfangdata.com.cn/Periodical/dqkx201708015 XU P, LI C H, LIU H C, et al. Fractal features of the effective gas transport coefficient for multiscale porous media[J]. Journal of China University of Geosciences, 2017, 42(8): 1373-1378. (in Chinese) http://d.old.wanfangdata.com.cn/Periodical/dqkx201708015 [13] CHEN S Y, CHEN H D, MARTNEZ D, et al. Lattice Boltzmann model for simulation of magnetohydrodynamics[J]. Physical Review Letters, 1991, 67(27): 3776. doi: 10.1103/PhysRevLett.67.3776 [14] 何雅玲, 王勇, 李庆.格子Boltzmann方法的理论及应用[M].北京:科学出版社, 2009. HE Y L, WANG Y, LI Q. Lattice Boltzmann method: theory and applications[M]. Beijing: Science Press, 2009. (in Chinese) [15] CHEN C, HU D D, MARTYSEVICH V N. Applications of high-resolution imaging and high-performance parallel computing in unconventional energy recovery[C]//Proceedings of the Abu Dhabi International Petroleum Exhibition and Conference. Abu Dhabi: Society of Petroleum Engineers, 2014. [16] PAN C X, LUO L S, MILLER C T. An evaluation of lattice Boltzmann schemes for porous medium flow simulation[J]. Computers & fluids, 2006, 35(8): 898-909. doi: 10.1016-j.compfluid.2005.03.008/ [17] ZHOU H X, YU X L, CHEN C, et al. Evaluating hydraulic properties of biochar-amended soil aggregates by high-performance pore-scale simulations[J]. Soil Science Society of America Journal, 2018, 82(1): 1-9. doi: 10.2136/sssaj2017.02.0053 [18] CHEN C, LAU B L T, GAILLARD J F, et al. Temporal evolution of pore geometry, fluid flow, and solute transport resulting from colloid deposition[J]. Water Resources Research, 2009, 45(6): W06416.1-W06416.12. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.1029/2008WR007252 [19] BHATNAGAR P L, GROSS E P, KROOK M. A model for collision processes in gases: I: small amplitude processes in charged and neutral one-component systems[J]. Physical Review, 1954, 94(3): 511. doi: 10.1103/PhysRev.94.511 [20] 杨艳霞, 李静. 3-D格子Boltzmann传质模型模拟生物膜降解有机污水[J].农业工程学报, 2018, 34(10): 225-230. doi: 10.11975/j.issn.1002-6819.2018.10.028 YANG Y X, LI J. Simulation of biomembrane degrading organic wastewater by 3-D lattice Boltzmann mass transfer model[J]. Transactions of the Chinese Society of Agricultural Engineering, 2018, 34(10): 225-230. (in Chinese) doi: 10.11975/j.issn.1002-6819.2018.10.028 [21] TAKAJI I, MASATO Y, OGINO F. Lattice Boltzmann simulation of flows in a three-dimensional porous structure[J]. International Journal for Numerical Mehods in Fluids, 1999, 29: 737-748. doi: 10.1002/(SICI)1097-0363(19990415)29:7<737::AID-FLD813>3.0.CO;2-H [22] INAMURO T, YOSHINO M, OGINO F. Accuracy of the lattice Boltzmann method for small Knudsen number with finite Reynolds number[J]. Physics of Fluids, 1997, 9(11): 3535-3542. doi: 10.1063/1.869426 [23] THOMPSON P A. Compressible-fluid dynamic[M]. New York:McGraw-Hill, 1971. [24] MULJADI B P, BLUNT M J, RAEINI A Q, et al. The impact of porous media heterogeneity on non-Darcy flow behaviour from pore-scale simulation[J]. Advances in Water Resources, 2015.[doi: 10.1016/j.advwatres.2015.05.019] [25] CHEN C, PACKMAN A I, GAILLARD J F. Pore-scale analysis of permeability reduction resulting from colloid deposition[J]. Geophysical Research Letters, 2008, 35: L07404. doi: 10.1029-2007GL033077/ [26] 蔡勇, 李光耀, 王琥. GPU通用计算平台上中心差分格式显式有限元并行计算[J].计算机研究与发展, 2013, 50(2): 412-419. http://d.old.wanfangdata.com.cn/Periodical/jsjyjyfz201302021 CAI Y, LI G Y, WANG H. Parallel computing of central difference explicit finite element based on GPU general computing platform[J]. Journal of Computer Research and Development, 2013, 50(2): 412-419. (in Chinese) http://d.old.wanfangdata.com.cn/Periodical/jsjyjyfz201302021 [27] THORNE D T. Lattice Boltzmann modeling: an introduction for geoscientists and engineers[M]. Berlin:Springer, 2006. [28] 施小清, 吴吉春, 袁永生, 等.含水介质各向异性对渗透系数空间变异性统计的影响[J].水科学进展, 2005, 16(5): 679-684. doi: 10.3321/j.issn:1001-6791.2005.05.011 SHI X Q, WU J C, YUAN Y S, et al. Effect of the anisotropy in porous media on the spatial variability of the hydraulic conductivity[J]. Advances in Water Science, 2005, 16(5): 679-684. (in Chinese) doi: 10.3321/j.issn:1001-6791.2005.05.011 [29] 周志芳, 庄超, 戴云峰, 等.单孔振荡式微水试验确定裂隙岩体各向异性渗透参数[J].岩石力学与工程学报, 2015, 34(2): 271-278. http://d.old.wanfangdata.com.cn/Periodical/yslxygcxb201502006 ZHOU Z F, ZHUANG C, DAI Y F, et al. Determing anisotropic hydraulic conductivity in fractured rocks based on single-borehole slug tests[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(2): 271-278. (in Chinese) http://d.old.wanfangdata.com.cn/Periodical/yslxygcxb201502006 [30] 夏源, 吴吉春, 张勇.改进时间分数阶模型模拟非Fick溶质运移[J].水科学进展, 2013, 24(3): 349-357. http://journal16.magtechjournal.com/Jweb_skxjz/CN/abstract/abstract2297.shtml XIA Y, WU J C, ZHANG Y. Tempered time-fractional advection-dispersion equation for modeling non-Fickian transport[J]. Advances in Water Science, 2013, 24(3): 349-357. (in Chinese) http://journal16.magtechjournal.com/Jweb_skxjz/CN/abstract/abstract2297.shtml [31] ZHANG X X, LYU M C. Persistence of anomalous dispersion in uniform porous media demonstrated by pore-scale simulations[J]. Water Resources Research, 2007, 43(7).[doi: 10.1029/2006WR005557] [32] HOU Y S, JIANG J G, WU J C. Anomalous solute transport in cemented porous media: pore-scale simulations[J]. Soil Science Society of America Journal, 2018.[doi: 10.2136/sssaj2017.04.0125] -
点击查看大图
图(6) / 表 (3)
计量
- 文章访问数: 672
- HTML全文浏览量: 148
- PDF下载量: 128
- 被引次数: 0
