雨水在沥青混凝土中饱和渗流的多松弛时间格子Boltzmann法模拟

杨大田,杨礼明,周 乾

(1.重庆交通大学 土木工程学院,重庆 400074;2.广西交通规划勘察设计研究院有限公司,南宁 530029)

沥青混凝土水损坏是沥青路面主要病害之一[1]。雨水在各种因素下渗入到沥青混凝土材料及路面结构内部[2-3]是沥青混凝土路面产生早期破坏的主要影响因素之一。沥青混凝土作为一种有孔多相复合材料,内部孔隙结构为雨水渗入和渗出以及滞留提供了可能。同时,车-水-路相互作用形成的路表有压水,使更多的雨水极易渗入到沥青混凝土内部及结构层层间。车辆荷载作用下,雨水在路面内部产生“瞬态脉动”水压[4],这种水压由静水压力和超孔隙水压力组成,超孔隙水压力是由孔隙结构变形导致孔隙中的水反耦合作用产生的[5]。水在多孔介质中的饱和渗流细观模拟常用格子Boltzmann方法[6-8]。Kutay[9]利用CT扫描技术构建沥青混凝土孔隙结构,用单松弛时间格子Boltzmann法模拟了水在沥青混凝土孔隙中渗流规律。

1 沥青混凝土二维数字图像构建

在160 ℃左右, 沥青混凝土被压路机压实, 但其内部存在一定孔隙, 其中一部分孔隙相互连通, 形成孔隙通道,与外界空气连通, 水能够自由流出; 另外一部分孔隙形成孤立的空穴。 压实沥青混凝土的孔隙结构,可用CT扫描技术获取。 CT扫描技术在国内刚处于起步阶段, 并且费用比较高。 国外对压实沥青混凝土的CT扫描, 发现压实沥青混凝土顶部和底部的孔隙率比中间部分的大[10], 形成盆状。 一般说来, 压实沥青混凝土,孔隙率越小, 孔隙的孔径越小, 相互连通的孔隙通道就越少。

本次模拟, 采用数字图像技术, 建立压实沥青混凝土的二维模型, 与外界大气连通的孔隙通道, 用“1”表示; 把孤立空穴和集料颗粒, 用“0”表示, 黑色代表固体颗粒, 白色代表孔隙。 一般来说,压实沥青混凝土的孔隙率范围分布较宽, 从密级配沥青混凝土的2%～6%, 到开级配沥青混凝土的18%～25%。 在本次模拟中,建立了3.9%、 6.0%、 9.1%、 12.5%、 15.8%、 18.6%和21.4%等7种孔隙率的压实沥青混凝土二维模型,在此仅列出了前3种孔隙率的二维模型，见图1。

图1 不同孔隙率下的压实沥青混凝土二维模型Fig.1 Two-dimensional model of compacted asphalt concrete indifferent porosity(593×1 256像素,1像素=8.466 7×10-5m)

2 多松弛时间格子Boltzmann 方法

2.1 多松弛时间格子Boltzmann 方法

格子Boltzmann方法利用一个粒子分布函数来描述流体分子集体运动行为。 在格子Boltzmann方法中,粒子被假定在一定规则的网格子上移动,并满足格子Boltzmann方程。 早期,最简单的格子Boltzmann方法是单松弛时间格子Boltzmann方法[11]。 但是,单松弛时间格子Boltzmann方法,数值不稳定、精准性差。

为了克服单松弛时间格子Boltzmann方法, 1992年d’Humières[12]提出了一种广义格子Boltzmann法(Generalized LBE, GLBE)。 2000年, Lallemand P等[13]对该模型作了详细的理论分析。

一个Q个离散速度的D维格子xi∈δxZD多松弛时间格子Boltzmann模型,一般表示为

fi(x+ciΔt,t+Δt)-fi(x,t)=-M-1·S·[m-meq],

(1)

式中:i=0,1,2,…,Q；M是一个Q×Q转换矩阵,它把粒子速度空间的粒子分布函数f∈V=RQ线性地传递到速度矩空间m∈M=RQ:

m=M·f,f=M-1·m,

(2)

S是一个非负Q×Q对角松弛系数矩阵。

多松弛时间模型演化为

碰撞:

(3)

迁移:

fi(x+ciΔt,t+Δt)=fi′(x,t)。

(4)

对于多松弛时间D2Q9模型,在式(4)中ci(i=0,1,2,…,8)分别为

c0=0;c1=(1,1);c2=(0,1);c3=(-1,0);

c4=(0,-1);c5=(1,1);c6=(-1,1);

c7=(-1,-1);c8=(1,-1)。

宏观密度和速度

(5)

(6)

对于多松弛时间D2Q9模型,速度矩

m=(ρ,e,ε,jx,qx,jy,qy,pxx,pxy),

(7)

对应的平衡态

meq=(1,-2ρ+3u2,ρ-3ρu2,jx,-jx,jy,

(8)

式中:ρ是密度;jx=ρux和jy=ρuy是流体动量J在x和y轴的分量, 在流体系统中这3个矩守恒, 其他几个矩不守恒, 它们的平衡矩是前面3个矩的函数。

对称松弛时间矩阵{si}

S=diag(sρ,se,sε,sx,sq,sy,sq,sυ,sυ)。

(9)

因为ρ、jx和jy是守恒变量,因此sρ和sx不影响系统,因而可取任何值。因此,松弛时间矩阵{si}取值为

S=diag(0,1.64,1.54,0,1.9,0,1.9,sυ,sυ)。

(10)

相应的转换矩阵M为

(11)

把作用力加入多松弛时间D2Q9模型中, 采用Guo等[14]作用力模型,因此,式(6)整理成

(12)

同时,式(1)重新整理成

(13)

在低马赫数下,格子Boltzmann方法实际上是一个弱压缩流方法,因此压力与密度有关,这个关系方程为状态方程

(14)

2.2 边界条件及体力

所取沥青混凝土模型厚约5 cm,宽10.6 cm,只是整个沥青混凝土路面层的一部分,因此用周期边界模拟雨水在模型中自由流进流出。 体力为重力加速度,取值9.81 N/kg或9.81 m/s2。

3 渗透系数计算

一般形式的二维达西定律方程可表示为

qx=-(KxPx/μ);

(15)

qy=-(KyPx/μ)。

(16)

式中:qx是x方向达西流量;qy是y方向达西流量;Kx是x方向渗透系数;Ky是y方向渗透系数;Px是x方向压力梯度;μ是动力粘度, 单位N·s/m2。

达西流量q与流体速度U通过有效介质孔隙率neff联系起来,即

q=neffU。

(17)

根据式(16)和(17),得到渗透系数

(18)

(19)

由式(18)和式(19)计算得到的渗透系数单位是m2, 还需通过式(20)转换为m/s:

k=Kg/v,

(20)

式中,g是重力加速度, m/s2;v是运动粘度, m2/s。

4 曲折度计算

水在沥青混凝土中流动路径是十分复杂和曲折,反映这种曲折和复杂,常用曲折度或迂曲度[7]

τ=〈|v|〉/〈vx〉,

(21)

式中: |v|是某点流体速度的绝对值;vx是速度在x方向的分量, 也即是压力梯度的方向, 在这里是指重力加速度的方向; 〈〉是在空间上的平均。

5 模拟结果与分析

5.1 孔隙率对平均渗流速度的影响分析

从图2可以看出,随着孔隙率增大,x方向平均渗透速度增大(x方向是重力加速度的方向, 下同);约在3%～14.9%范围,x方向平均渗透速度变化比较缓慢;在孔隙率14.9%后,x方向平均渗透速度急剧增大。这充分说明14.9%是x方向平均渗透速度变化的一个临界孔隙率。

图2 x方向渗透速度与孔隙率关系曲线Fig.2 Relationship between porosity and seepage velocity at x-direction

图3是不同孔隙率条件下的x方向平均渗透速度。 从图3a、 b可以看出, 在x方向, 最大速度出现在孔径比较狭窄的地方; 对于孔隙率比较大的模型(图3c、 d), 在汇集孔道中出现最大速度。

从图4可以看出,y方向平均渗透速度与孔隙率相关性较小, 但总体规律表现为, 随着孔隙率增大,先增大后减小, 在孔隙率14%左右出现最大值, 这与x方向平均渗透速度的临界孔隙率较为一致。

从图5可以看出,迂曲度与孔隙率相关系数R2只有0.685,呈线性变化,孔隙率越大,迂曲度越小。

5.2 孔隙率对平均渗透系数的影响分析

从图6可以看出,孔隙率与x方向平均渗透系数存在较强的相关性,随着孔隙率增大,x方向平均渗透系数增大,当孔隙率超过15.6%时,沥青混凝土模型出现超渗透性,对应的渗透系数为6.374×10-3m/s。这与实际沥青混凝土渗透试验结果存在一定差距,但与文献[15]相比较,是合理的。

图3 不同孔隙率下模型中某点x方向的平均渗透速度Fig.3 Average seepage velocity of a certain point at x-direction in different porosity

图4 y方向渗透速度与孔隙率关系曲线Fig.4 Relationship between porosity and seepage velocity at y-direction

图5 迂曲度与孔隙率关系曲线Fig.5 Relationship between porosity and tortuosity

图6 x方向渗透系数与孔隙率关系曲线Fig.6 Relationship between porosity and permeability at x-direction

从图7可以看出,孔隙率与y方向平均渗透系数也存在较强的相关性,随着孔隙率增大,y方向平均渗透系数增大,当孔隙率超过15.6%时,沥青混凝土出现超渗透性,对应的渗透系数为3.325×10-3m/s。可以发现y方向的渗透系数比x方向的渗透系数约小2倍。

5.3 孔隙壁剪应力分析

在重力作用下, 雨水不仅对集料表面上的沥青膜有压力作用, 此这压力等于大气压力, 即静态压力,一般不会造成沥青膜的脱落;同时还有剪应力作用,其中剪应力是造成沥青膜在集料表面上剥离脱落主要因素之一。因此,分析不同孔隙率模型的孔隙壁上最大剪应力,可揭示沥青混凝土的水损坏的机理。

图7 y方向渗透系数与孔隙率关系曲线Fig.7 Relationship between porosity and permeability at y-direction

从图8可以看出,在重力作用下,孔隙率越大,孔隙壁剪应力越大。实际上,孔隙壁剪应力与孔隙通道中速度梯度有关。这可从计算剪应力公式[9]可以得到解释：

(22)

式中，μ是运动粘度；τ是剪应力。

以孔隙率9.1%模型为例, 观察剪应力在孔隙通道中分布情况。 从图9可以看出,在孔隙通道内的剪应力都主要集中在0～0.025 Pa, 只有少数孔隙通道的剪应力突然增大如A、B和C点。A点是两个孔隙通道汇成一个孔隙通道上的一个点,A点的孔径约4像素, 在A点前段, 孔径大小约9像素。 同样地,B点处于两个孔隙通道汇成一个孔隙通道上,B点孔径约3像素, 在B点前部分的孔径约8像素。C点也是一样的。 由于孔隙的孔径由大变小, 引起孔隙中渗流速度发生变化, 这个变化幅度越大,产生剪应力越大。 因此, 在这3个点处的剪应力突然增大, 图中还隐隐约约地观察到其他几个类似点。 因此, 可以说沥青混凝土路面水损坏首先发生在孔隙的孔径比较小的地方。

图10是A、B、C3点的剪应力放大图。3处的剪应力在0.20～0.27 Pa,小于沥青膜与石灰石集料表面在6 h浸水后的粘结强度1.0 MPa[16]。因此,在重力作用下雨水在沥青混凝土孔隙中流动,不会引起沥青膜从集料表面脱落。

图8 孔隙率与孔隙壁剪应力关系曲线Fig.8 Relationship between porosity and shear stress of pore wall

图10 不同点处的剪应力分布Fig.10 Shear stress distribution in different sites

6 结 论

利用数字图像方法,建立了7种沥青混凝土的数字模型,代表了不同沥青混凝土孔隙。利用多松弛时间格子Boltzmann法模拟了雨水在重力作用下在沥青混凝土渗流过程,发现在重力作用下,孔隙率约为14.9%时雨水在x方向渗透速度突然增大;y方向渗透速度与孔隙率相关性较差;随着孔隙率增大,模型孔隙通道复杂性越小;同时,发现孔隙率约为15.6%时x方向渗透系数和y方向渗透系数突增,x方向渗透系数比y方向渗透系数约大2倍。还发现剪应力突然增大地方一般是孔隙孔径突然变狭窄的地方,这可能是沥青混凝土路面发生水损坏开始地方。