ECAH理论模型及数值计算探析

4.1.1.1　ECAH理论模型

该理论模型基于守恒型水动力学方程，首先做如下假设：

（1）采用Navier-Stokes形式的动量守恒方程，热应力由于在本问题中影响甚微，可以忽略。

（2）整个系统为一个准稳态过程，即不计温度和压力随时间的变化。

（3）将黏性系数和热传导系数当作常数处理。

（4）应力应变张量满足以下关系：

对颗粒悬浊液中的声衰减进行理论分析，可以通过质量、动量和能量守恒定律，应力应变关系以及热力学状态方程求得压缩波、剪切波、热波在弹性、各向同性、导热的固体介质以及连续相介质中的波动方程。波动方程在球坐标下求解，按照球贝塞尔函数和球谐函数的级数展开，其中包含了待定散射系数。在颗粒与介质界面运用边界条件，可以得到一个6阶的线性方程组。求解这一方程组即可得到与复波数有关的散射系数。

一般情况下，向量v可以写成标量势ψ的梯度和矢量势A的旋度的和的形式，即

对于压缩纵波、热波和剪切波，分别得到三个关于各自复波数的亥母霍兹波动方程：

其中，kc、kT和ks为复波数，它们的定义如下：

式中　ω——超声角频率；

c——压缩波声速；

αL——纵波声衰减系数；

σ——热扩散系数；

ρ——密度；

η——剪切黏度。

参照第3章中的图3.1，在球体内部和球面向外散射的压缩波、热波和剪切波都分别满足前面的波动方程，对它们求解可以得到下面的通解形式。

连续相：

颗粒相：

式中　jn——球贝塞尔函数；

hn——球汉克尔函数；

Pn——勒让德多项式；

带“'”的参数——颗粒内部参数；

An、Bn、Cn、A'n、B'n、C'n——六个待定系数。

在颗粒与连续相界面r=a，其中a为颗粒半径，根据速度、应力分量以及温度、热流连续的边界条件：

对压缩波与热波有

其中，热相关系数Γ由下式给出：

式中　γ——比热比，γ=Cp/Cv。

β——热扩散系数。

式（4.2）与式（4.3）中的近似结果在式（4.4）的条件下成立：

对剪切波有

将通解形式代入边界条件，得到六个线性方程组：

式中　ac、as、aT——分别为颗粒半径R与压缩波、黏性波和热波波数的乘积；

jn、hn——分别为球贝塞尔函数和球汉克尔函数；

带“'”的参数——球贝塞尔函数和球汉克尔函数的导数（即j'n和h'n），或

表示颗粒内部参数（如a'c）。由于剪切波和热波在液体中迅速衰减而不能被远场探测器检测到，因此最终颗粒两相系超声衰减系数α仅仅与压缩波散射系数An有关，即(www.xing528.com)

式中　R——颗粒半径；

ψ——颗粒相体积浓度。

由式（4.6）可以看出，超声衰减系数的获得就归结为对系数An的求解，但其很难有一个显式的表达式。将式（4.5）变为矩阵形式：

其中

式（4.7）即为待求线性方程组，其中待求向量中An、Bn、Cn分别指代压缩波、热波和剪切波散射系数，上述系数的“'”指代球体内部值。等号右边变量αc=kcR为无因次量。

4.1.1.2　数值计算

表4.1给出水与玻璃颗粒的物理参数，对于颗粒相与连续相均分别需要七个物理参数。

表4.1　水与玻璃颗粒的物理参数（25℃）

对于玻璃微珠的水悬浊液，取超声频率1 MHz，颗粒半径由1～100μm递增做数值计算。从图4.1结果看出，随着颗粒半径增加，矩阵条件数相应剧烈增加。为分析原因，进一步观察矩阵M中的元素变化情况。

图4.1　矩阵条件数（数值越大表明矩阵越接近于奇异）

表4.2中列出了半径增至100μm时矩阵中元素数值。可以看到，矩阵第二、三列元素在颗粒粒径较大时迅速减小，而第五列元素则迅速增大。这导致了矩阵元素间的数值上差异非常大，其条件数也很大。

分析原因，以矩阵M第二列为例，各元素包含了球汉克尔函数hn（αT）及其1阶、2阶导数项。当粒径增大，其自变量aT=kTR相应增加且具有很大的虚部，对应表4.2的参数，aT=468.93+468.93i。此时球汉克尔函数值将随着其自变量线性增加呈指数规律递变，如aT继续增加，该列元素则会很快超过双精度数值值域2.2×10-308～1.8×10308，导致数值溢出，计算失败。

图4.2　ECAH模型计算中数值溢出范围（双精度计算条件）

图4.2给出了ECAH模型计算时发生数值溢出的范围。可以看出，即使对于超声频率1 MHz，对应的可计算颗粒半径不超过150μm，且随着超声频率增加，对应粒径减小。而在超声波颗粒粒径检测中，为了获得更全面的衰减谱信息，有时使用的超声频率可以达到50 MHz甚至100 MHz，显然模型当前计算域无法满足要求。

表4.2　矩阵M中元素数值

4.1.1.3　计算方法改进

由于颗粒粒径较大或者频率较高（即kR较大），数值溢出问题会限制ECAH模型的可能计算范围，降低模型实用性。本研究提出通过构造与球贝塞尔函数和球汉克尔函数有关的变量，改善矩阵M进而从根本上拓展ECAH模型计算上限。

考虑贝塞尔函数（汉克尔函数也适用），其具有如下递推关系：

式中　z——函数自变量；

n——阶数。

为此引入几个与之有关的变量及其递推公式。

变量Mn（z）：

并有

变量Ln（z）：

变量Qn（z）：

图4.3为不同变量随自变量z递增的变化趋势，可以看出随着复自变量的增加，球汉克尔函数及其1阶、2阶导数均快速递减并趋于零，但是新构造的变量Ln（z）、Qn（z）却趋于数值1，因此用其替代原有函数可以避免数值溢出问题。采用新变量还使得系数矩阵的条件数始终稳定在一定数量级上，不会随着自变量递增。如对于f=1 MHz，R=100μm，矩阵条件数3.7×1026与原矩阵M在粒径1～10μm时条件数相当。

图4.3　不同变量随自变量增加的变化[n=1，自变量z=w（1+i）]

在引入了上述新的变量后，对矩阵M进行优化，即对其第二、三、五列的元素分别除以hn（αT）、hn（αs）和jn（α'T）项，由于这里贝塞尔函数自变量均为复数，不会出现除零点情况，数值上稳定。相应地，新构造矩阵中各列元素由数值1替代了hn（z），新变量Ln（z）和Qn（z）替代原有的h'n（z）和h"n（z）。应注意，这里的改进方法因为重构了矩阵M，形如

但只要系数An结果不变，即可确保最终超声衰减系数的求解是正确的。

图4.4　超声衰减系数随粒径变化情况

图4.5　不同粒径时的超声衰减谱

如图4.4所示，以1%体积浓度的玻璃微珠悬浊液为例计算，对于颗粒粒径较小的情况，改进算法和式（4.7）的直接求解方法结果是一致的，但是当颗粒粒径增加且频率较高（如100 MHz曲线），原矩阵条件数过高，结果发生偏差。粒径继续增加，直接求解法发生数值溢出，图中数据点中断。而采用本研究改进方法，计算范围拓展到了超声频率100 MHz和粒径1 000μm，即原来范围10余倍以上。图4.5中对于粒径10μm颗粒，两方法的超声衰减谱直到50 MHz完全吻合；但是当粒径增加至50μm时，直接求解法只能计算得9 MHz内频谱；当粒径为100μm，其频谱范围进一步减小至约2 MHz；对于粒径为1 000μm时，直接求解完全失效，采用改进方法仍获得频率至50 MHz的频谱。