基于MPS方法的自密实混凝土流动数值模拟

李敬军， 田 雷， 邱流潮， 韩 宇

(中国农业大学 水利与土木工程学院， 北京 100083)

自密实混凝土(self-compacting concrete，SCC)[1]具有抗离析、高流动性和强填充黏结的特性，被广泛应用在建筑行业之中．SCC的流动性是其施工性能的重要量度，且最终的流动状态对硬化后的混凝土耐久性和强度有直接影响．室内试验一般很难真实再现施工环境的复杂性，因此数值模拟方法就成了一个较好的研究手段．

对于自密实混凝土的数值模拟来说，最重要的是对混凝土自由表面的计算．古川恭雄等[2]最早使用黏塑性有限单元法(VFEM)来描述均一连续介质，使用黏塑性悬浮单元法(VSEM)来描述非连续介质，从而通过数值模拟对混凝土的流动性进行了研究．Gram等[3]通过使用离散单元法(DEM)商业软件PFC3D对SCC的流动性进行模拟分析，获得了较高的精度．黄绵松[4]自主开发了离散元计算程序，并设计了正交仿真试验，对自密实砂浆性能进行了模拟分析，从而建立了使用DEM方法模拟自密实砂浆时参数确定的拟合公式．陈松贵[5]和张传虎等[6]通过使用基于Bingham流体的LBM-DEM数值模拟方法，对SCC在堆石孔隙中的流动进行了研究．

考虑到SCC流动的大变形性质和混凝土填充行为的复杂性，近年来发展起来的移动粒子半隐式法(moving-particle semi-implicit method，MPS)为混凝土流动性模拟提供了一个新的方法．MPS方法是日本东京大学的Koshizuka等[7]最早提出的一种用于计算不可压缩流体的粒子方法，在模拟流动和填充行为方面具有特殊的优势．MPS方法的提出，引起了国内外学者高度关注，并且已经被广泛应用于描述波-结构的相互作用[8]、溃坝问题[7]、液仓晃动[9-10]和入水砰击[11]等牛顿流体问题，也被应用于非牛顿流体[12-13]的研究当中．

本文使用MPS方法来研究SCC的流动性，采用Bingham流变模型对SCC在L形箱中的流动过程进行了三维数值模拟，对SCC的模拟流动状态和流动时间进行对比分析，以验证MPS方法的可行性与精度．在此基础上，本文还分析研究了塑性黏度对L形箱试验和V形漏斗试验流动状态的影响，建立了SCC在V形漏斗试验中下落时间与SCC密度、塑性黏度之间的关系公式，验证了所建立公式的精度，为SCC数值模拟过程中塑性黏度参数的确定提供了方法．

1 基本理论

1.1 MPS基础理论

1.1.1控制方程

MPS方法的控制方程包括连续性方程和N-S方程.对于不可压缩流体，可写成如下形式：

(1)

(2)

式中：ρ为流体密度,kg/m3；t为流动时间，s；u为速度向量，m/s；P为压力,Pa；v是流体的运动黏性系数；g是重力项．

1.1.2核函数

在MPS方法中，核函数用以描述粒子间的相互作用．常用的核函数W(r)由Koshizuka等[7]给出：

(3)

式中：re为有效半径，m；r为两粒子间的距离，m．

1.1.3梯度模型

通过将作用域内粒子间的位置矢量进行加权平均，可以获得MPS的梯度模型：

(4)

式中：φ为函数；D为空间维数，取D=3；n0为初始的粒子数密度；r为粒子的坐标矢量．

1.1.4Laplacian模型

N-S方程中的黏性项和压力Poisson方程的求解,都通过对拉普拉斯算子进行空间离散来进行计算：

(5)

(6)

式中：λ是与流体表面张力和接触角有关的参数，λ的引入是为了确保数值结果与扩散方程的解析相一致．

1.1.5压力Poisson方程

压力Poisson方程为：

(7)

式中：k为步数；n*为临时时刻的粒子数密度；Δt为时间步长，s．

1.1.6自由表面的判断

MPS中自由表面的判断采用以下粒子密度判别：

ni<βn0

(8)

式中：ni是粒子的粒子数密度；β是常数,本文取β= 0.97．

1.1.7固体边界条件

为防止流体粒子在边界处产生缺失或者穿透，本文在边界处设置了多层边界粒子．为获得流体粒子在边界处正确的压力梯度，将靠近流体粒子的第1层边界粒子进行压力Poissin方程的计算，其他边界粒子按下式进行计算：

(9)

1.2 SCC的Bingham流变本构模型

在模拟过程中，本文将SCC简化为各向同性的均质流体，并且用Bingham流变模型对其进行描述[14]．塑性黏度(ηp)和屈服强度(τy)是Bingham流变模型的2个主要参数，只和混凝土的自身特性有关．在试验过程中可通过坍落扩展度试验和相应公式进行测定计算．其中Roussel等[15-16]给出了坍落扩展度与混凝土屈服应力τy之间的换算关系：

(10)

(11)

式中：V为材料总体积，m3；R为扩展度半径，m；对于一般自密实混凝土λ可取0.005，如果忽略表面张力，式(10)中后一项可以去掉；S为坍落度，cm．公式适用在坍落度为5～25cm时的试验结果．

在坍落扩展度试验中，材料到达某一规定扩展度的时间，可以间接反映材料塑性黏度的大小，时间越长，则材料塑性黏度越小.流动时间与塑性黏度ηp的关系为[17-18]：

(12)

式中：H为坍落度测量桶的高度，m；T为从流动到静止的总时间，s．

为防止计算的不稳定性，本文采用双黏度模型来对流体进行求解，模型将流体视为处于流体状态的黏塑性流体，以及处于刚性状态的高黏性流体．

(13)

2 数值模型的验证

2.1 试验概况

采用文献[19]中的L形箱法来测定SCC的流动性，同时记录SCC流到固定距离(lx=200、300、400、500、600、700mm)的时间．L形箱的尺寸及形状如图1所示.

图1 L形箱尺寸、形状示意图Fig.1 Size and shape of L-box(size:mm)

2.2 计算参数及模型

三维MPS计算模型如图2所示，其中MPS颗粒直径为7mm，粒子总数为31527个，混凝土密度为2300kg/m3,塑性黏度ηp=100Pa·s，屈服应力τy=1600Pa．

图2 L形箱试验的三维模型Fig.2 3-D model of L-box test

2.3 模拟结果与讨论

图3为ηp=100Pa·s时L形箱试验的MPS三维数值模拟流动过程图．模拟从挡板一次性打开开始，SCC在自重的作用下从垂直段向水平段进行扩展，L形箱的垂直部分可清晰模拟出SCC黏附L箱壁面的特征，反映了三维模拟过程中可以考虑壁面与底面的黏性效应．

图3 L形箱试验模拟流动过程图Fig.3 Snap shots of 3D L-box at different instants of time

在重力作用下，SCC从L形箱垂直段向水平段流动，图4为MPS模拟与试验结果对比．由图4可见：模拟值与试验值吻合较好，但在流动前期，试验值要比模拟值偏小，这是因为试验状态下打开隔板活动门的时间延迟，而数值模拟可以达到隔板活动门瞬间打开的理想状态．

图4 试验结果与模拟结果的对比Fig.4 Comparison between simulation results and experimental results

Bingham流变模型中塑性黏度系数反映了混凝土的流动性能．塑性黏度越大，SCC流动越慢，反之亦然．因此，通过调整塑性黏度值进行L形箱流动性能试验的模拟，可以研究Bingham流变模型下塑性黏性对于SCC在L形箱中的流动性能的影响．本文分别设置了4种工况来探究SCC在L形箱中的流动性能，分别为ηp=10、40、50、100Pa·s．

图5为不同塑性黏度条件下，t=2.00s时SCC在L形箱水平段的自由表面线和水平段流动距离与流动时间的关系．由图5(a)可见：流动停止之前的相同时刻下，随着塑性黏度的增大，SCC的流动距离越来越短．图5(b)对比了到达相同距离处不同塑性黏度下SCC的流动时间，SCC的流动速度也随着SCC塑性黏度的减小而增大，成负相关的趋势．也直接说明了塑性黏度是SCC流动速度的重要影响因素．

图5 塑性黏度对SCC在L形箱中流动状态的影响Fig.5 Effect of plastic viscosity on flow state of SCC in L-box

3 V形漏斗试验塑性黏度影响数值模拟分析

3.1 V形漏斗试验模型及计算参数

在SCC性能的研究中，SCC流过V形漏斗的时间是反映其流动性能的重要参数，它可以很好地反映出SCC在流动中的表观黏性．一般来讲，相同密度的SCC流过V形漏斗的时间越长，其表观黏性越大．不同工况下V形漏斗试验数值模拟计算结果可以反映Bingham流变模型下塑性黏度与SCC流动性能的关系．V形漏斗试验装置及尺寸如图6所示[4]．

图6 V形漏斗尺寸、形状示意图Fig.6 Size and shape of V-funnel(size:mm)

三维MPS计算模型如图7所示，其中MPS颗粒直径为5mm，粒子总数为70056个，混凝土密度为2300kg/m3．

图7 V形漏斗试验的三维模型Fig.7 3-D model of V-funnel test

3.2 模拟结果与讨论

图8给出了ηp=20Pa·s时V形漏斗试验的MPS三维数值模拟流动过程图．在自重作用下，SCC从V形漏斗下端出口流出．由图8可见：由于边壁效应，V形漏斗内SCC在边壁处流速较小，中间处流速较大，从而出现液面中间低两边高的形态，与试验结果相同，说明MPS可以较好地模拟SCC在V形漏斗中的掉落过程，并再次验证了基于Bingham流变模型的MPS方法对于SCC的流态可以进行较好地模拟．

本文设置密度和塑性黏度2个参数影响的正交试验，来探究SCC模拟参数对于V形漏斗试验的影响，并将ρ=2300kg/m3,t=3s时的SCC在V形漏斗的流动状态绘制成图9．由图9可见：相同时刻下，相同密度的SCC塑性黏度越高，流动性能越差，流出V形漏斗的速度越小，相同时间内的流出量越少，V形漏斗内部中心液面高度值也越大，这与理论分析趋势相同．

不同塑性黏度下SCC的V形漏斗下落试验模拟的各工况下落时间如表1所示．

图8 V型漏斗试验模拟流动过程图Fig.8 Snap shots of 3D V-funnel at different instants of time

图9 塑性黏度对SCC自由表面的影响Fig.9 Relationship between plastic viscosity and flow state of SCC

表1 不同密度与不同塑性黏度下SCC的V形漏斗模拟下落时间

对模拟结果进行线性拟合，得到下落时间T与SCC密度、塑性黏度的关系为：

T=-0.0054ρ+0.3086ηp+13.6303

(14)

相关度R2为0.984．式(14)即为Bingham塑性黏度与SCC流动性能的对应关系．分析可知，相同密度的SCC在V形漏斗中的下落时间随塑性黏度增大而增大，并基本呈线性关系．

为验证公式精度，本文对文献[4]中的3种V形漏斗的试验工况分别进行数值模拟，试验参数如表2所示，其中Te为V形漏斗下落时间的试验值．

表2 SCC配合比、试验值及模拟参数值

对V形漏斗不同试验工况进行数值模拟，将其试验下落时间与模拟下落时间进行对比分析，两者误差如表3所示，其中Ts为V形漏斗模拟值，E为误差率．

表3 模拟精度分析

通过式(14)对试验中塑性黏度值进行计算，并使用MPS方法对其进行数值模拟．对比分析结果可知，模拟结果最大误差为5.9%，说明本文建立的方程具有一定的精度，且可以为以Bingham流变模型描述SCC的流动状态时，对塑性黏度参数的标定提供参考依据．

从上述分析可知，本文提出的基于Bingham流变模型的MPS方法可以较好地模拟SCC的流动过程并可以对SCC的最终流动状态进行预测．通过使用文中方法将流变参数分析与实际工程相结合，可以为SCC的配合比计算、施工方案设计及施工技术指导提供一定的参考依据．

4 结论

(1)基于Bingham流变模型的MPS方法可以比较准确地模拟SCC的流动状态，且可以考虑到SCC与试验装置壁面之间的黏性效应．

(2)塑性黏度对于SCC在L形箱中的流动时间有较大的影响，且随着塑性黏度的增加，SCC流动到指定长度的时间也随之延长．

(3)建立了V形漏斗下落时间与SCC密度、塑性黏度之间的关系方程，为以Bingham流变模型对SCC进行数值模拟计算时的塑性黏度参数选取提供了参考依据．