混凝土中氯离子扩散的二维分数阶数值分析

Numerical Analysis on the Two-Dimensional Fractional Diffusion Model of Chloride Ion in Concrete

摘要: 氯离子对混凝土结构的侵蚀是影响其耐久性的重要原因，而扩散作用是氯离子侵入到混凝土中的最主要方式之一。对于混凝土中氯离子扩散的二维问题，首先基于Caputo分数阶导数建立了变扩散系数的时间分数阶扩散模型，随后结合L1算法建立了数值求解该模型的有限差分格式。接着，利用已有实验数据对分数阶导数的选取进行了分析，并对氯离子在混凝土中的扩散进行了预测。

Abstract: The erosion of chloride ions on concrete structures is an important reason to affect its durability, and diffusion is one of the most important ways for chloride ions to invade concrete. This paper is focusing on the two-dimensional problem of chloride ion diffusion in concrete. Firstly, the temporal fractional diffusion model of Chloride ion with variable diffusion coefficient in concrete is estab-lished based on the Caputo fractional derivative. Next, the finite difference scheme combined with L1 algorithm is established. At last, the value of fractional derivatives is analyzed using existing ex-perimental data, and the diffusion of chloride ions in concrete is predicted.

1. 引言

众所周知，钢筋混凝土具有抗压性能好、无收缩、剪切强度高、强度大、正常情况下不易变形等优点，这些优点使得钢筋混凝土被广泛应用于土木工程、水利建设、港口加固、桥梁搭建等领域中，造就了其在工程结构设计中的重要地位。然而氯离子对钢筋侵蚀进而导致混凝土结构提前失效已成为全世界普遍关注的问题 [1] 。氯离子能够进入到混凝土之中也不是凭借单一的途径，比如扩散作用、渗透作用、毛细作用、电迁移作用、吸附作用或者是以上两种或者多种方式的结合，而最主要的氯离子侵入混凝土的方式为前三种，而且扩散作用起主导作用 [2] 。

目前对于氯离子在混凝土中的扩散作用的描述大都是基于Fick第二扩散定律，常见的二维模型为：

$\frac{\partial C}{\partial t}=D\left(\frac{{\partial }^{2}C}{\partial x^2}+\frac{{\partial }^{2}C}{\partial y^2}\right)$ (1)

其中C为氯离子的浓度，t为时间，D为扩散系数，x和y分别代表扩散的两个方向。关于扩散系数D，Thomas等 [3] 提出氯离子的扩散系数是时变函数D_t的观点，即：

$D_t=D_0{\left(\frac{t_0}{t}\right)}^m$ (2)

式中D₀和D_t分别为扩散时间为t₀和t时混凝土的氯离子扩散系数，m为时间依赖性常数， $t\ge 0$ ，其中t₀时测定混凝土的氯离子扩散系数为D₀。

近年来，分数阶微积分已成为当前国际上的一个热点研究课题，因为分数阶微积分算子是一个具有记忆效应的非局部算子，它对时间上的历史依赖性和空间上的长程相关性。因而基于分数阶导数的扩散模型也倍受关注。Chen等 [4] 考虑了氯离子和混凝土之间的相互作用，使用朗缪尔等温线描述了氯离子浓度。Wei等 [5] 建立了Caputo型氯离子时间分数阶扩散模型，通过均方位移的方法近似得到扩散系数。杨帅 [6] 以孔隙介质中的非达西渗流和反常扩散现象为研究对象，借助分数阶微积分理论，系统地对孔隙介质中的低速、高速非达西渗流，溶质反常输运等进行了研究。王一晗等 [7] 建立了混凝土中变扩散系数的时间变分数阶的氯离子一维扩散模型。

本文将Caputo型分数阶导数引入模型(1)，并考虑时变扩散系数(2)，建立了变扩散系数的时间分数阶二维氯离子扩散模型，随后结合L1算法建立了数值求解该模型的有限差分格式。接着，利用已有实验数据对分数阶导数的选取进行了分析，并对氯离子在混凝土中的扩散进行了预测。

2. 扩散模型

基于Fick第二定律和Caputo型分数阶导数 [8] ，考虑氯离子在混凝土中的二维扩散问题，建立变系数时间分数阶扩散模型，描述为：

$\{\begin{array}{l}\frac{{\partial }^{\alpha }C}{\partial t^{\alpha }}=D_0{\left(\frac{t_0}{t}\right)}^m\left(\frac{{\partial }^{2}C}{\partial x^2}+\frac{{\partial }^{2}C}{\partial y^2}\right)\\ C\left(x,y,0\right)=C_0\\ C\left(0,y,t\right)=C_s,C\left(+\infty ,y,t\right)=C_0\\ C\left(x,0,t\right)=C_s,C\left(x,+\infty ,t\right)=C_0\end{array}$ , (3)

式中 ${\partial }^{\alpha }/\partial t^{\alpha }$ 为Caputo分数阶算子；x和y分别为两个方向侵蚀深度；t为混凝土结构暴露于氯离子环境中的时间， $t\ge 0$ ； $C\left(x,y,t\right)$ 表示t时刻 $\left(x,y\right)$ 位置处的氯离子浓度；D₀和D_t分别为扩散时间为t₀和t时的扩散系数；α为阶数， $0<\alpha \le 1$ ；C_s为混凝土表面氯离子浓度；C₀为氯离子内部初始浓度；C₀和C_s是常数；m为时间依赖性常数。

3. 数值格式

本节使用有限差分格式和L1算法相结合的方法，求解上面所建立的变系数时间分数阶扩散模型(3)。

首先定义 $x_i=i\Delta x$ ， $i=0,1,\cdots ,M$ ， $y_j=j\Delta y$ ， $j=0,1,\cdots ,N$ ， $t_k=k\Delta t$ ， $k=0,1,\cdots ,R$ ，其中 $\Delta x$ 和 $\Delta y$ 为空间步长， $\Delta t$ 为时间步长。模型(3)中的整数阶导数项在点( $x=x_i$ , $y=y_j$ , $t=t_{k-}{}_{1/2}$ )的离散格式为：

$\begin{array}{c}\frac{{\partial }^{2}C}{\partial x^2}|{}_{{}_{{}_{\left(x_i,y_j,t_{k-1/2}\right)}}}=\frac{C\left(x_{i+1},y_j,t_k\right)-2C\left(x_i,y_j,t_k\right)+C\left(x_{i-1},y_j,t_k\right)+C\left(x_{i+1},y_j,t_{k-1}\right)}{2\Delta x^2}\\ +\frac{C\left(x_{i-1},y_j,t_{k-1}\right)-2C\left(x_i,y_j,t_{k-1}\right)}{2\Delta x^2}+\text{Ο}\left(\Delta x^2+\Delta t^2\right),\end{array}$ (4)

$\begin{array}{c}\frac{{\partial }^{2}C}{\partial y^2}|{}_{{}_{{}_{\left(x_i,y_j,t_{k-1/2}\right)}}}=\frac{C\left(x_i,y_{j+1},t_k\right)-2C\left(x_i,y_j,t_k\right)+C\left(x_i,y_{j-1},t_k\right)+C\left(x_i,y_{j+1},t_{k-1}\right)}{2\Delta y^2}\\ +\frac{C\left(x_i,y_{j-1},t_{k-1}\right)-2C\left(x_i,y_j,t_{k-1}\right)}{2\Delta y^2}+\text{Ο}\left(\Delta y^2+\Delta t^2\right).\end{array}$ (5)

分数阶导数项借助L1算法进行数值离散，L1算法的定义见 [8] ，模型(3)中的分数阶导数项在点( $x=x_i$ , $y=y_j$ , $t=t_{k-}{}_{1/2}$ )的离散格式为：

$\begin{array}{l}\frac{{\partial }^{\alpha }C}{\partial t^{\alpha }}|{}_{{}_{{}_{(x_i,y_j,t_{k-1/2})}}}=\frac{\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}[\frac{C\left(x_i,y_j,t_k\right)+C\left(x_i,y_j,t_{k-1}\right)}{2}-{\alpha }_{k-1}\frac{C\left(x_i,y_j,t_0\right)+C\left(x_i,y_j,t_1\right)}{2}\\ -{\displaystyle \underset{n=1}{\overset{k-1}{\sum }}\left({\alpha }_{k-n-1}-{\alpha }_{k-n}\right)\frac{C\left(x_i,y_j,t_n\right)+C\left(x_i,y_j,t_{n-1}\right)}{2}}]+O\left(\Delta t^{2-\alpha }\right),\end{array}$ (6)

其中 ${\alpha }_s={\left(s+1\right)}^{1-\alpha }-s^{1-\alpha },s=0,1,2,\cdots R$ 。

最后，在( $x=x_i$ , $y=y_j$ , $t=t_{k-}{}_{1/2}$ )处离散扩散模型(3)，将(4) (5) (6)代入化简差分格式整理可得：

$\begin{array}{c}\left(r_1+2r_2^k+2r_3^k\right)C_{i,j}^k-r_2^k{C}_{i+1,j}^k-r_2^k{C}_{i-1,j}^k-r_3^k{C}_{i,j+1}^k-r_3^k{C}_{i,j-1}^k\\ =-\left(r_1+2r_2^k+2r_3^k\right)C_{i,j}^{k-1}+r_2^k{C}_{i+1,j}^{k-1}+r_2^k{C}_{i-1,j}^{k-1}+r_3^k{C}_{i,j+1}^{k-1}+r_3^k{C}_{i,j-1}^{k-1}\\ +r_1{\alpha }_{k-1}\left(C_{i,j}^0+C_{i,j}^1\right)+r_1{\displaystyle \underset{n=1}{\overset{k-1}{\sum }}\left({\alpha }_{k-n-1}-{\alpha }_{k-n}\right)\left(C_{i,j}^n+C_{i,j}^{n-1}\right)}\end{array}$ (7)

其中，

$r_1=\frac{\Delta t^{-\alpha }\Delta x^2\Delta y^2}{\Gamma \left(2-\alpha \right)},r_2^k=D_0{t}_0^m{\left(k-\frac{3}{2}\right)}^{-m}{\left(\Delta t\right)}^{-m}\Delta y^2,r_3^k=D_0{t}_0^m{\left(k-\frac{3}{2}\right)}^{-m}{\left(\Delta t\right)}^{-m}\Delta x^2.$ (8)

$C_{ki,j}$ 表示点( $x=x_i$ , $y=y_j$ , $t=t_k$ )处氯离子浓度的数值解。结合模型(3)中的初始条件和边界条件，求解方程组(7)，就可求得 $C_{ki,j}$ 。

4. 工程算例

张伟等 [9] 对某码头的构件的氯离子含量进行了测定。构件保护层厚度50 mm，全长为400 m，初始氯离子浓度C₀为0，表面氯离子浓度C_s为0.29%，扩散系数D₀为3.48 × 10⁻⁶ mm²/s。全部工程完工15年后，对150 mm处混凝土中不同深度x (mm)的氯离子浓度进行了三次测量，结果如表1所示。

Table 1. Measured values of chloride ion concentration C (%) at different depths x (mm) in concrete at y = 150 mm [9]

表1. y = 150 mm处混凝土中不同深度x (mm)的氯离子浓度C (%)的三组测量值 [9]

4.1. 参数分析

对于模型(3)，工程上t₀一般取为28天，时间依赖性常数m和分数阶导数的阶数α对氯离子在混凝土中的扩散有很大影响。选取时间为15年，将不同的m和α代入数值格式(7)得到图1。

Figure 1. Curve of chloride ion concentration under different m and α at y = 150 mm

图1. 不同m和α下y = 150 mm处混凝土中氯离子浓度曲线

从图1(a)可以看出，当时间依赖性常数m固定时，随着分数阶导数的阶数α的增大，扩散越来越快，恰好说明了分数阶微分对时间上的历史依赖性。而当分数阶导数的阶数α固定时，随着时间依赖性常数m的增大，导致扩散系数减小，从而减缓了扩散速度，如图1(b)所示。其中m = 0时，就是扩散系数为常数的情况。

4.2. 参数选取

为了进一步预测氯离子扩散趋势，对于模型(3)中的时间依赖性常数m和分数阶导数的阶数α，根据表1中的值，我们定义

$均方误差E=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{10}{\sum }}{\left(计算值-测量平均值\right)}^2}}{10}}$ ， (9)

其中计算值是求解数值格式(7)得到的结果，测量平均值为三次测量值的平均。不同m和α取值时计算均方误差E，部分如表2所示。结合上面参数分析的结果，为了使均方误差E尽可能小，固定m和α中的一个取值，可以确定另外一个参数的取值。比如，m = 0.10时，选取α = 0.94。

Table 2. Mean squared error E at different m and α

表2. 不同m和α得到的均方误差E

Figure 2. Comparing computing value with measruing value of chloride ion concentration under different m and α at y = 150 mm

图2. 不同m和α下y = 150 mm处混凝土中氯离子测量值和计算值比较

按照这种选取方法确定m和α的取值，图2中绘制了在y = 150 mm处混凝土中氯离子测量值和计算值比较曲线。比如图2(b)所示，m = 0.05, α = 0.91时，黑色曲线为15年后在y = 150 mm处不同深度氯离子的浓度分布。不难看出，计算值在一定程度上较好地吻合了测量值。

4.3. 扩散预测

根据模型(3)和数值格式(7)，以及上述方法确定的参数，就可以对氯离子的扩散情况加以预测。图3预测了两组参数分别在20年后和30年后氯离子浓度在构件内分布的等高线图，(a)和(b)为m = 0, α = 0.89时的分布情况，(c)和(d)为m = 0.15, α = 0.96时的分布情况。显然，随着时间的推移，氯离子扩散深度在不断地增加。氯离子含量为0.25的位置在30年时比20年时的情况距左边界更近, 而且差别较为明显，前后十年混凝土构件的被腐蚀程度区别已经较为显著。

Figure 3. Contour line distribution of predicted chloride ion diffusion in concrete

图3. 混凝土中氯离子扩散预测值的等高线分布

5. 结论

本文建立了混凝土中变系数时间分数阶二维氯离子扩散模型(3)，以及求解该模型的有限差分数值格式(7)。结果表明：时间依赖性常数m的增大会减缓扩散；增大分数阶导数的阶数α会加快扩散。数值预测说明：氯离子在混凝土的扩散作用十分明显，为保证建筑物的安全必须引起足够的重视。

基金项目

北京建筑大学教育科学研究项目(Y2205)。