混凝土结构中氯离子扩散分析的精细积分法

杨绿峰,洪 斌,高 钦,曾建聪

(广西大学工程防灾与结构安全教育部重点实验室,广西南宁 530004)

钢筋混凝土结构较长时间暴露于海洋等富含氯离子环境后,氯离子会侵入混凝土中并逐渐到达钢筋表面,造成钢筋锈蚀,从而影响钢筋混凝土结构的使用性能和耐久性。由于钢筋锈蚀的过程并不持续消耗氯离子,所以分析评估氯离子环境下钢筋混凝土结构耐久性和使用寿命的关键是确定氯离子在混凝土中的分布和扩散规律。吴庆令等[1]通过现场海洋暴露试验,研究了混凝土中氯离子的扩散特性。Collepardi等[2]用Fick第二扩散定律描述氯离子在混凝土中的扩散过程,得到在一定初始条件和边界条件下的数学解;Liang等[3]考虑了水泥水化物对氯离子扩散的影响,并对氯离子扩散系数进行了修正;宋子健等[4]研究了溶液成分对混凝土中氯离子扩散系数的影响;彭国军等[5]考虑骨料形状,对混凝土氯离子扩散系数进行了数值预测。余红发等[6]提出了考虑多因素扩散模型的一维解析解。杨绿峰等[7]研究了混凝土时变条件下氯离子扩散的封闭解。

由于实际的钢筋混凝土结构通常具有复杂的几何形状和边界条件,解析方法在应用中有较多困难,人们通常采用数值方法解决此类问题。在Fick第二定律及其数学模型的基础上,施养杭等[8]提出了氯离子侵入混凝土计算的有限差分法模型;Funahashi[9]也在氯离子扩散分析的有限差分法方面进行了研究;Sergi等[10]应用最小平方法(leastsquares methods)对氯离子在混凝土中的扩散规律进行了一维数值模拟;Han[11]采用有限元法分析氯离子扩散问题,并推导出具体计算格式;杨绿峰等[12]研究建立了氯离子扩散分析的边界元法。但是上述氯离子扩散分析的有限元法普遍存在两个问题,其一是忽视了基于半无限大扩散场的Fick第二定理及其边界条件对数值分析模型的影响,并造成计算误差;其二是在用有限元求解氯离子扩散微分方程时,通常采用显式的中央差分法或隐式的Newmark法、Wilson-θ法等,但计算结果的稳定性和计算精度时常不尽如人意。时域精细积分方法[13]可以将时域积分转换为矩阵相乘,能够取得非常高的计算精度;文献[14]利用三次样条插值函数模拟积分项中的被积函数,建立了求解非齐次动力方程特解的一种精细数值积分法。

由于氯离子侵入混凝土并发生扩散的过程非常缓慢,直接利用精细数值积分格式有时会带来明显的误差。本文利用误差函数建立了混凝土中氯离子扩散场的补偿长度和补偿系数,保证了有限元法在该类问题中的计算精度和计算效率;同时,在一致分布矩阵的基础上建立了集中分布矩阵,验证了集中分布矩阵具有更好的计算精度。在此基础上利用精细积分技术研究建立了混凝土中氯离子扩散分析的精细积分有限元法,克服了普通有限元法分析氯离子扩散时存在的问题。

1 氯离子一维扩散的有限元方程

氯离子在混凝土中的扩散过程可以用Fick第二定律来描述:

其边界条件和初始条件分别为

式中:D为混凝土中氯离子扩散系数;t为混凝土持续暴露于氯离子环境中的时间;C=C(x,t)表示t时刻混凝土试件内深度为x处氯离子浓度;C0为初始氯离子浓度;Cs为混凝土构件表面氯离子浓度。根据变分原理可知,上述控制方程是泛函的变分极值条件。

将扩散场沿 x方向离散为N个单元,可以建立典型单元(图1)上t时刻的浓度分布函数:式中:N为形函数矩阵;C为单元节点浓度参数列阵 ,且有 :

图1 扩散单元

C=(C1,C2)TN=(N1,N2)

式中:C1,C2分别为结点1和2上的氯离子浓度;Ni为插值函数 ,其中 N1=1-,N2=(=x/le为无量纲化单元局部坐标,le为单元长度)。

将式(4)代入式(3),可得单元泛函:式中:Ωe表示单元域;﹒C为C对时间t的一阶导数;N′为N对坐标x的一阶导数。

根据变分原理 δ Π=0可得混凝土中氯离子扩散分析的有限元方程:

其中

式中:M为氯离子分布矩阵;K为氯离子扩散矩阵。

当单元沿着与扩散方向相正交的另一个方向的尺寸取为单位值时,可求得式(7)为

2 一致分布矩阵和集中分布矩阵

计算式(7)中分布矩阵M所用到的形函数与计算扩散矩阵K所用到的形函数一样时,称 M为一致分布矩阵。由于求解方程(6)时采用一致分布矩阵M有时会导致计算结果不稳定,这里尝试通过近似数值积分方法重新计算氯离子分布矩阵M。

考虑数值积分近似计算的梯形公式:式中:l为积分域的长度;k为积分点,取积分域的两个端点;Fk(Ni,Nj)为被积函数F(Ni,Nj)在积分点k上的值,结合式(7)可以定义被积函数:

根据形函数的性质可以求得分布矩阵:

式(12)中只有对角元素不为零,且等于式(9)中相应行的全部元素的叠加。称式(12)的分布矩阵为集中分布矩阵。

按照普通有限元法标准步骤,可以将单元分布矩阵和单元扩散矩阵集成为总体分布矩阵和总体扩散矩阵,相应地可将单元扩散方程集成为总体扩散方程。以下第3节、第4节的计算格式都是针对混凝土中氯离子扩散分析的总体有限元方程。

3 精细积分法

由于集中分布矩阵 M存在逆矩阵,所以由式(6)集成的总体有限元方程可以改写为式中A为定常矩阵,且A=-M-1K。

按照微分方程求解理论,齐次方程(6)的通解为

式中eAt为矩阵At的指数函数。当将时间域[t0,ta]离散为p个等步长的子域,则时间步长Δt和时间域内的第i个离散结点分别为

与之相对应,在ti时刻,结点氯离子浓度向量表示为 Ci,且有:

式中指数矩阵 T=eAΔt。

根据矩阵加法定理有:

其中m为正整数,可选择m=2N,N为正整数,因此m通常是非常大的正整数。由于Δt是有限大的时间步长,则τ=Δt/m通常非常小。由此可得:

式中I为单元矩阵。

将式(18)代入式(17),有:

将上述分解计算过程重复进行至第 n次时,有:

其中

由此可以看出,当 T0确定后,可以按照迭代计算公式(21)求出 Tn,当 n=N时,代入式(20)求得指数矩阵T,从而避免在计算过程中发生大数吃小数的情况;进而根据迭代计算公式(16),可以求解时间域上各个离散结点的氯离子浓度参数列阵Ci。

以上迭代求解方程的过程将时程积分问题转化为矩阵相乘问题,简化了计算,称之为混凝土中氯离子扩散分析的精细积分有限元法。

4 混凝土中氯离子扩散场的补偿长度

式(1)的解析解为

式中erf(s)为误差函数,且:

由于Fick第二定理所导出的解析模型,要求远端边界取在无穷远处,并满足远端边界条件C(x=+∞)=C0,如式(2)所示。但有限元、边界元等数值方法所建立的计算模型通常要求求解域为有限大,远端边界不可能取在无穷远处,因此对于氯离子一维扩散分析的有限元、边界元等数值方法[12],需要定义氯离子扩散场的补偿长度L,如图2所示,以此确定远端边界的位置,并保证远端边界条件C(x=L)=C0仍然能够得到满足。

图2 混凝土试件长度及补偿长度

将式(23)的误差函数用图表示,见图3。假定远端边界位于x=y=L处,如图2中虚线所示,它能够保证扩散方程的远端边界条件成立,由式(24)

图3 误差函数y=erf(s)

可得

根据式(25)可得

式中,k=2erf-1(1),即有:

式(23)定义的误差函数如图3所示,可以看出,当1.5≤s≤2时误差函数接近于1。因此,通过比较式(23)和式(27),容易看出式(27)成立时必然要求3≤k≤4。这里,k称为混凝土的氯离子扩散长度补偿系数,简称为补偿系数。通常情况下建议取k=3.5。

根据式(26)可以确定计算扩散场的补偿长度L,从而能够保证远端边界条件的成立。补偿长度L表示了混凝土的氯离子扩散问题中,正确计算混凝土内部各点的氯离子浓度所需的扩散场的最小长度,如图2所示。图2实线 ABCD表示试件的实际边界,如果试件的实际长度 l≥L,则按照试件的实际长度建立有限元计算模型。反之,L≥l,那么在建立有限元法等数值方法的计算模型时,需要将图2所示的混凝土试件的远端边界BC虚拟地移动至虚线B′C′处,从而在保证远端边界条件仍能成立的同时,方便有限元等数值方法计算模型的建立。

5 算例分析

算例1混凝土试件的尺寸为200mm,混凝土表面的氯离子浓度(氯离子与混凝土的质量百分比,下同)Cs=1.0%,初始氯离子浓度C0=0,氯离子扩散系数D=31.536mm2/a。分别采用基于集中分布矩阵的精细积分有限元法(LCM)和基于协调分布矩阵的精细积分有限元法(CCM)计算扩散时间分别为50a和80a时试件不同深度的氯离子浓度值,并将计算结果和封闭解(CFS)相比较。图4给出了扩散时间50a时的LCM,CCM和CFS的结果。

图4 扩散50a时氯离子浓度分布

从图4可以看出,LCM和CCM的结果都与CFS比较接近,能够给出正确的氯离子浓度分布规律,其中LCM的结果更为准确。为了便于比较LCM和CCM,这里通过和CFS比较计算两种方法的误差,见图5,可以看出,LCM的计算结果精度优于CCM。

图5 LCM和CCM的相对误差比较

算例2制作一批混凝土试件,尺寸为150mm×150mm×200mm,养护28d后将5个面用环氧树脂封闭,只留下1个150mm×150mm的面暴露于海水中,混凝土的氯离子扩散系数D=9.252×10-12m2/s=291.77mm2/a,试件中氯离子初始浓度C0=0,试件表面氯离子浓度Cs=0.565%,利用精细积分有限元法(PIFEM)计算暴露时间分别为50a和80a时试件内部的氯离子分布浓度,并与封闭解进行比较。在精细积分有限元法的计算模型中,取时间步长Δt=0.1a。

首先按照试件实际长度建立精细积分有限元法离散模型,计算试件暴露50 a和80a时的氯离子浓度分布,并将计算结果与CFS相比较,如图6(a)所示,可以看出此时PIFEM计算结果与CFS之间的误差相当大。

图6 算例2计算结果

按照本文提出的扩散场补偿长度公式(26)确定有限元法离散模型(图2),并重新利用PIFEM计算。当暴露时间分别取50a和80a时,根据式(26),可以求得扩散场补偿长度分别为544mm和688mm,二者都大于试件的实际长度,所以应根据扩散场补偿长度建立有限元离散模型。将计算结果同CFS相比较,如图6(b)所示,可以看出,按照本文提出的扩散场补偿长度理论建立有限元法计算模型,计算结果和封闭解吻合很好,能够取得很高的计算精度。

6 结 论

a.基于混凝土中氯离子扩散场的补偿长度及其表达式,建立了有限元法离散模型,克服了传统数值方法离散模型中半无限大假设导致的计算误差,从而保证计算结果的正确性。

b.将一致分布矩阵转换为集中分布矩阵,并与精细积分技术相结合,研究了混凝土中氯离子扩散分析的有限元法控制方程的求解方法,不仅保证了计算结果的稳定性,而且能够取得很好的计算精度。

[1]吴庆令,余红发,王甲春,等.现场海洋区域环境中混凝土的Cl-扩散特性[J].河海大学学报:自然科学版,2009,37(4):410-414.

[2]COLLEPARDI M,MARCIALIS A,TURRIZZANIR.The kinetics of penetration of chloride ions into the concrete[J].Cemento,1970(4):157-164.

[3]LIANG M T,WANG K L,LIANG C H.Service life prediction of reinforced concrete structures[J].Cement and Concrete Research,1999,29:1411-1418.

[4]宋子健,蒋林华.溶液成分对混凝土中氯离子扩散系数影响的研究综述[J].水利水电科技进展,2011,31(5):82-87.

[5]彭国军,郑建军,周颖琼.考虑骨料形状时混凝土氯离子扩散系数预测的数值方法[J].水利水电科技进展,2009,29(6):13-16.

[6]余红发,孙伟,鄢良慧,等.混凝土使用寿命预测方法的研究[J].硅酸盐学报,2002,30(6):686-701.

[7]杨绿峰,李冉.混凝土中氯离子二维扩散规律的解析研究[J].水利水电科技进展,2009,29(3):20-23.

[8]施养杭,罗刚.含多种因素的氯离子侵入混凝土的有限差分计算模型[J].工业建筑,2004,34(5):7-10.

[9]FUNANSHI M.Prediction corrosion-free service life of a concrete structure in a chloride environment[J].ACI Materials Journal,1990,87(6):581-587.

[10]SERGI H,PAGE C.Diffusion of chloride and hydroxyl ions in cementitious materials exposed to a saline environment[J].Magazine of Concrete Research,1992,44:63-69.

[11]HAN S H.Influence of diffusion coefficient on chloride ion penetration of concrete structure[J].Construction and Building Materials,2007,21:370-378.

[12]杨绿峰,陈正,王炎炎,等.混凝土中氯离子二维扩散分析的边界元法[J].硅酸盐学报,2009(7):1110-1117.

[13]ZHONG W X,WILLIAMS F W.A precise time step integration method[J].Journal of Mechanical Science,1994,208:427-430.

[14]富明慧,刘祚秋,林敬华.一种广义精细积分法[J].力学学报,2007,39(5):672-677.