超声波在粘弹性混凝土介质中传播机理研究

2019-07-12鲁光银张志勇朱自强付渝渝

物探化探计算技术 2019年3期

鲁光银，张志勇.2，朱自强，付渝渝

(1.中南大学地球科学与信息物理学院，长沙 410083; 2.河南省交通规划设计研究陈院股份有限公司，郑州 450000)

0 引言

近年来，随着各级建设工程要求地提高，混凝土成为各种工程建筑最常用的一种结构材料[1]。在混凝土施工过程中，由于混凝土材质、外界因素的影响，在内部结构中往往会有缺陷存在[2]，但混凝土的质量对整个结构工程至关重要，关乎到工程的使用安全性和使用期限[3]。因此，在混凝土浇筑完成之后，需要对混凝土的质量进行相应的检测，确保其安全性。同时，在对混凝土检测时，又要尽量保证混凝土原有的质量尽可能地不受到破坏[4]。因此，常用无损检测进行混凝土质量的判断。而超声波具有很强的穿透能力，方向性强，检测设备简单，仪器操作方便，节省人力和检测成本，分辨率较高，因此在混凝土无损检测领域应用最为广泛[5]。

为了更确切地研究超声在混凝土中的传播规律，不同的数值计算方法逐渐被引用到超声数值模拟中来。模拟超声波传播常用的数值方法有：有限元法、弹性动力学有限积分技术、有限差分法等[6]。

笔者主要介绍有限差分法在混凝土正演中的应用进展。有限差分法在弹性波数值模拟方面应用最为广泛，因此逐渐地被应用到超声波的数值模拟中。Day等[7]使用Pade近似法进行了二维时间域的数值模拟；Emmerich等[8]提出了广义标准线性固体，推导了粘弹性模量的有理式;Robertsson等[9]提出了速度一应力有限交错网格法，时间上二阶精度，空间上四阶精度，用这种方法研究了二维空间中横波和纵波的传播性质与规律，并从二维推广到三维，取得了较好的效果；H.Yamawaki等[10]在常规的有限差分交错网格的基础上，提出了改进交错网格节点的计算方案，使得有限差分法更能适应于各项异性介质；魏东等[11]应用时域有限差分法，研究了弹性固体中的二维声场特性，采用PML(完全匹配层)吸正收边界条件，解决了有限区域的计算问题；张俊等[12]基于有限差分法，设置建立了含有多个探头的相控阵超声检测模型；密士文[13]在混凝土超声数值模拟波动方程中加入粘弹性介质的拉梅常数项。但其只是推导出速度—应力公式，未对超声正演模拟效果做出相应的研究。

相对于弹性介质，其受到应力和产生应变是瞬发的，且满足Hook定律中应变与应力成线性关系[14]。并且，当失去外部应力后，弹性体会迅速恢复到原来的状态。而对于粘弹性介质，其物理性质介于弹性与粘性之间。在受到外力作用时，粘弹性介质会永久发生形变，且在外力消失后，不会恢复到原来的状态。这种性质不仅跟粘弹性介质自身因素有关，还与压力大小、温度等外界因素有关。粘弹性介质对混凝土中能量的衰减机理比较复杂。

混凝土作为一种多项非均匀复合材料，超声在其内部传播，在做数值模拟时，弹性动力方程已经不能更确切地模拟超声声波的传播方式。笔者分析总结前人工作，把有限差分法应用到混凝土超声波数值模拟上来。将混凝土视为粘弹性介质，基于Kelvin粘弹性介质[15]，详细推导出了超声在粘弹性介质中的声波方程，在时间二阶、空间任意偶数精度下的交错网格差分格式，利用数值模拟的方法来模拟混凝土中超声波的传播过程。

1 混凝土中超声波传播的数值模拟

粘弹性体常见的基本模型有：Maxwell 模型、Kelvin模型、标准线性模型和达朗贝尔模型。

目前在处理粘弹性介质中波动方程问题主要使用两种方法：①复模量法(粘性参数法)；②添加粘弹性介质品质因子Q值的方法。在全波列模拟中，通常采用的是粘性参数法，而在波场、波形分析时，则常采用Q值法[16]。笔者主要进行声波传播某一时刻波场快照分析和提取某一道声波信号进行分析，因此，采用的是添加Q值的方法。

1.1 Kelvin粘弹性介质

二维的Kelvin粘弹性声波方程的一阶速度—应力方程表示形式为：

(1)

式(1)中的拉梅系数计算式为：

λ+2μ=ρv2

(2)

(3)

其中:λ和μ表示弹性介质的拉梅系数；λ′和μ′表示粘性介质的拉梅系数；ν表示模型介质的纵波速度；ω为圆频率。

在粘弹性介质中，用Q表示纵波的品质因子，反映模型介质对于超声波能量的吸收情况。式(4)中，f表示超声波的频率，C为超声波的相速度，a表示超声波振幅随距离的衰减系数。品质因子Q与频率成函数关系，不同频率的超声波的能量衰减情况不同[16]。频率越高，能量衰减越大，相对应的Q值越大。但品质因子Q的值表示的是粘弹性介质本身的能量衰减属性，只与粘弹性介质物理属性有关，与频率的大小无关，因此，在进行数值模拟时，通常是采用Q值固定的方法。品质因子Q值越大，表示介质的粘滞性越小，即介质属性越趋近于弹性。当Q值为无穷大时，式(3)可表示为λ′+2μ′=0，即变成了完全弹性介质声波方程的表达形式。

1.2 高阶精度差分粘弹性声波方程

采用交错网格法计算粘弹性声波方程，分别取质点在和的应力和速度的值进行计算。得到速度—应力在时间一阶导数上的二阶精度差分格式：

对式(4)进行交错网络剖分，i表示网格计算节点的深度方向序号，j表示节点水平方向序号，k表示时间方向序号，i+1/2、j+1/2、k+1/2分别表示相对应的半程节点序号。分别在时间节点k和空间节点j+1/2位置处进行计算水平速度分量vx和垂直速度分量vz。

{px(t+Δt2)=px(t-Δt2)+(λ+2μ)∂vx∂x+(λ′+2μ′)∂vx∂x

pz(t+Δt2)=pz(t-Δt2)+(λ+2μ)∂vz∂z+(λ′+2μ′)∂vz∂z

vx(t)=vx(t-Δt)+Δtρ∂p∂x

vz(t)=vz(t-Δt)+Δtρ∂p∂z

(4)

为了便于说明，在差分格式中，将水平速度分量vx设为U，垂直速度分量vz设为W；水平加速度分量∂U∂t=∂vx∂t设为U1，垂直加速度分量∂W∂t=∂vz∂t设为W1；设P、S、T分别表示应力p、水平应力分量px和垂直应力分量pz。得出的时间2阶，空间2N阶精度的交错网络粘弹性质声波方程的差分格式表示为：

水平应力分量为：

S|k+12i,j=S|k-12i,j+Δt(λ+2μ)(i,j)

1Δx∑Nn=1CNn[Uk(i,j+2n-12)-Uk(i,j-2n-12)]+

Δt(λ' + 2μ')(i,j)1Δx∑Nn=1CNn

[Uk1(i,j+2n-12)-Uk1(i,j-2n-12)]

(5)

垂直应力分量：

T|k + 12i,j=T|k-12i,j+ Δt(λ+ 2μ)(i,j)

1Δz∑Nn=1CNn[Wk(i+2n-12,j)-Wk(i-2n-12,j)]+

Δt(λ' + 2μ')(i,j)1Δz∑Nn=1CNn

[Wk1(i + 2n-12,j)-Wk1(i-2n-12,j)]

(6)

水平速度分量：

U|ki,j+12=U|k-1i,j+12+1ρi,j+12

ΔtΔx∑Nn=1CNn(Pki,j+n-Pki,j-n+1)

(7)

垂直速度分量：

W|ki+12,j=W|k-1i+12,j+1ρi+12,j

ΔtΔz∑Nn=1CNn(Pki+n,j-Pki-n+1,j)

(8)

水平加速度分量：

U1|ki,j+12=1ρi,j+12ΔtΔx∑Nn=1CNn(Pki,j+n-Pki,j-n+1)

(9)

垂直加速度分量：

W1|ki+12,j=1ρi+12,jΔtΔz∑Nn=1CNn(Pki+n,j-Pki-n+1,j)

(10)

1.3 完全匹配层吸收边界

图1 稳定系数与空间差分精度关系曲线Fig.1 Relationship between stability coefficient and spatial differential accuracy

采用完全匹配层吸收边界(PML),在计算区域四周区域选取50层差分网格厚度的吸收层，用来吸收和衰减向模型之外传播的波，假设在建立的模型区域与匹配层的交界面处产生最小的虚假反射[17]。

1.4 稳定性条件

在时间2M阶、空间2N阶精度下，二维声波交错网格形式下的有限差分方程的稳定性条件为：

0≤∑Mm=1{(-1)m-1(2m-1)!(Δt2S)2m·

[∑Nn=1CNn(-1)n-1]2m}≤1

(11)

差分精度选取的阶数不同，稳定性也不同，精度要求越高，稳定性条件选取的越严格[18]。当空间差分精度达到12阶时，稳定系数变化值已经很小，再增加空间差分精度对稳定系数的影响不大，因此，本文在做有限差分计算时，选取空间12阶差分格式。

2 数值模拟分析

通过对混凝土模型进行定量化和数字化，对混凝土中的超声波场进行计算，模拟声波传播特性。通过数值模拟可以真实有效的重复计算超声在混凝土中的传播过程及变化情况，分析超声波频率、相位、振幅及波的类型、能量等物性参数的变化，为推断混凝土内部结构可能存在的缺陷提供借鉴。

2.1 混凝土骨料模型

骨料是混凝土中除水泥浆外所占体积比最大的成分，骨料形状不规则，分布无规律。模型尺寸设置为400 mm×200 mm规格，骨料最大颗粒直径为Amax=14 mm(图2)。

从图3可以看出：在均匀水泥砂浆介质模型中，振幅衰减曲线平滑，振幅衰减值与传播距离基本成正比关系，可以推测，均匀介质中能量衰减应该只与波前扩散的距离有关。在含骨料混凝土介质中，衰减曲线不平滑，衰减值无规律，较均匀水泥砂浆介质衰减更快，这是由于骨料分布无规律、界面不规则，超声在骨料界面发生反射、透射及多重散射，能量衰减程度也无法确切判断。图4为两种介质模型(200,50)处质点振动对比图，从图4中可以看出：两种介质中质点振动速度不同，且混凝土中的质点相位相对比较滞后，这是由于混凝土自身的粘弹性，使得介质中质点在震源波消失后，仍会处于振动状态，不会立即静止，在这一时间段内存在剩余的应变，不仅会造成能量的衰减，还会对信号产生干扰。

图2 混凝土模型Fig.2 Concrete model

图3 水泥浆、混凝土衰减对比Fig.3 The attenuation comparison between cement slurry and concrete

图4 质点振动图Fig.4 Particle vibration map

图5 波场快照对比图Fig.5 The wavefield snapshot(a)t=24 μs时刻波场快照；(b)t=36 μs时刻波场快照

图5为超声波在含有骨料的混凝土模型中传播不同时刻的波场快照。如图5(a)所示，是声波纵波、横波的初至波波前面的位置，此刻超声波包络面相对比较规则连续；图5(b)所示，声波传播范围内，由于受不规则骨料的影响，超声在骨料界面发生多重散射、发射等现象，造成波形包络面不连续，且包络面内波形出现紊乱。为了更加清晰地分析超声传播过程中能量的衰减情况，抽取了混凝土模型横向中心位置(200,50)处的单道波形，如图6所示：超声信号整体上振幅逐渐减小，与波前扩散的衰减有关；但振幅衰减不规则变化，这是由于在混凝土介质中受骨料的多重散射、反射的影响，造成的能量衰减。相对于均匀水泥浆介质，混凝土中振幅不再线性衰减。

2.2 空洞模型

混凝土砌筑过程中，由于工艺技术、人为失误及环境等因素的影响，常会有空洞缺陷的存在，为了研究混凝土中空洞对超声传播特性的影响，建立含空洞的混凝土模型。模型尺寸设置为400 mm×200 mm规格，骨料最大颗粒直径为Amax=14 mm。空洞半径为R0=20 mm(如图7)。

图6 抽取单道声波波形记录Fig.6 The single channel waveform record

图7 含空洞混凝土模型Fig.7 Concrete model including round hole

图8 波场快照对比图Fig.8 The wavefield snapshot(a) t=28 μs时刻波场快照；(a) t=62 μs时刻波场快照

图9 抽取单道声波波形记录Fig.9 The single channel waveform record

图10 含裂隙混凝土模型Fig.10 Concrete model including horizontal cracks

图8为超声声波在含空洞混泥土模型中传播时不同时刻的波场快照。从图8(a)可以看出，t=28 μs时刻，超声纵波已经完全传播到空洞下空间，但因为混凝土与空洞两者之间波阻抗相差很大，声波在空洞界面发生全反射，在空洞周围形成了明显的反射波包络面。同时根据惠更斯原理，超声声波在空洞界面上，以新的点源的形式继续传播，形成新的波前圆弧包络面，但因为混凝土的衰减作用，这一部分的声波能量降低，使得圆弧包络面不连续，特别是空洞正下方的波形，由于全反射出现部分波形缺失。如图8(b)所示：t=62 μs时刻，纵波初至波波前面传播到模型下界面，被吸收层吸收，基本无反射波的形成；横波初至波波前面传播到空洞下界面，同样，空洞正下方的波形由于受全反射的影响出现缺失；此时形成的横波的全反射包络面较明显，而t=28 μs时刻形成的纵波的全反射波，由于能量的衰减，纵波的全反射包络面不再明显。为了更加清晰地分析波形传播过程中能量的衰减情况，抽取了在含空洞混凝土模型横向中心位置(200,50)处的单道波形(图9)：超声信号在t=22 μs时刻纵波初至波被换能器接收，波形振幅达到最大，在t=39 μs时刻在空洞上界面发生全反射，在t=47 μs时刻在空洞下界面形成反射波，但由于上界面的全反射作用，下界面反射波振幅较小；在t=54 μs时刻形成的波峰，是声波横波在空洞上界面形成的全反射形成的，此后，横波在空洞中无法继续向下传播。相对于图6混凝土模型中的单道信号：整体上超声波能量也是随着传播距离逐渐衰减，但前期振幅减小较小，且振幅变化紊乱，是由于声波在空洞上界面发生全反射，能量主要聚集在空洞上方，也造成波形传播到空洞下方时，振幅整体表现较小。空洞垂直方向上的范围大致可以通过两波峰间的时间差确定。

图11 波场快照对比图Fig.11 The wavefield snapshot(a)t=26 μs时刻波场快照；(b)t=56 μs时刻波场快照

图12 抽取单道声波波形记录Fig.12 The single channel waveform record

图13 含倾斜裂隙混凝土模型Fig.13 Concrete model including tilted cracks

2.3 裂隙模型

在混凝土浇筑过程中，由于水泥砂浆凝固过程中的涨缩特性，混凝土在横向和纵向上会存在应力的作用，为了防止涨缩力过大，通常会留有横向和纵向的缝隙来缓冲涨缩力的作用。为了研究横向裂隙和空洞对混凝土中声波传播影响的差异，模型尺寸设置为400 mm×200 mm规格，骨料最大颗粒直径为Amax=14 mm。裂隙长度为50 mm，宽度为10 mm。裂隙的最小尺寸不能小于声波最小波长(图10)。

图11是超声声波在含裂隙混凝土模型中传播不同时刻的波场快照。从图11(a)可与看出：超声纵波传播到裂隙位置处，由于裂隙与混凝土之间波阻抗差异较大，声波会在裂隙上界面发生全反射，纵波形成的强反射波与横波的初至波相互干扰，造成横波波前的包络面紊乱，同时由于裂隙上界面对纵波的强反射，造成裂隙正下方的纵波初至波包络面缺失。图11(b)显示，超声横波传播到裂隙后，在裂隙上界面同样会发生反射，但横波在裂隙中无法传播，造成裂隙正下方的横波包络面波形缺失；同时因为裂隙横向尺寸较大，在裂隙上界面发生全反射的范围较广，传播到裂隙下方的声波较少，形成一个能量较弱的梯形区域。为了更加清晰地分析波形传播过程中能量的衰减情况，抽取了在含裂隙混凝土模型横向中心位置(200,50)处的单道波形，如图12所示：波形整体上能量随传播距离逐渐衰减，超声信号在t=22 μs时刻的振幅波峰为纵波初至波，在t=32 μs时刻的振幅波峰为纵波在裂隙上界面全反射形成的反射波，由于裂隙宽度较小，纵波在裂隙中传播的时间较短，在t=38 μs时刻形成的波峰即为纵波在裂隙下界面形成的反射波，随后纵波振幅逐渐减小；在t=55 μs时刻出现的振幅较大，是横波在裂隙上界面形成的强反射波。相对于图9含空洞混凝土模型中抽取的单道信号，整体上振幅在前期衰减较慢，但在t=55 μs时刻之后，振幅迅速衰减，且减小值较大，这是由于纵波、横波在裂隙界面形成全反射，且裂隙横向截面较大，全反射范围较大，使得裂隙上、下空间能量分布存在较大差异。

图14 波场快照Fig.14 The wavefield snapshot(a)t=36 μs时刻波场快照；(b)t=72 μs时刻波场快照

图15 抽取单道声波波形记录Fig.15 The single channel waveform record

2.4 倾斜裂隙模型

为了对比研究横向裂隙和空洞对混凝土中声波传播影响的差异，设置倾斜裂隙模型，模型尺寸设置为400 mm×200 mm规格，骨料最大颗粒直径为Amax=14 mm。倾斜裂隙长度为99 mm，宽度为9.66 mm。裂隙的最小尺寸不能小于声波波长(图13)。

图14是超声声波在含倾斜裂隙混凝土模型中传播不同时刻的波场快照。由图14(a)可知：在t=36 μs时刻，超声声波传播到倾斜裂隙位置处，超声波在倾斜裂隙上界面发生强反射，反射波向背倾角的方向传播，形成明显的反射弧形包络面，在倾斜裂隙下方，出现初至波包络面的错位、缺失，形成一处能量分布较小的区域，随着波形的传播，与裂隙的接触面积越大，该区域面积越大。由图14(b)可知：根据惠更斯原理，超声在倾斜裂隙界面以新的点源的形式向外继续传播，t=36 μs时刻在倾斜裂隙下方形成的能量分布较小的区域，由于新的点源形成的波形的传播，能量逐渐增大，逐渐形成比较完整、连续的包络面。为了更加清晰地分析波形传播过程中能量的衰减情况，抽取了在含倾斜裂隙混凝土模型横向中心位置(200，50)处的单道波形，由图15可知：超声信号在t=22 μs时刻为纵波初至波，显示的波形振幅最大；t=37 μs时刻为纵波在裂隙上界面发生强反射形成的反射波，t=44 μs时刻的波峰，推测是横波初至波引起的。对比图12中含横向裂隙混凝土模型抽取的单道信号：倾斜裂隙由于倾角的存在，反射波会向背倾角的方向传播，能量向背斜方向发散，在裂隙正上方(200,50)处接收到的能量较少，使得振幅衰减较大，振幅值较小。