Finite Element Analysis of Dynamic Behavior of Microbubble with Different Wall Shapes
-
摘要: 为研究刚性平面、凸面和凹面附近的微泡动力学行为的差异，建立超声激励下3种刚性壁面附近微泡的有限元模型. 结果显示，刚性凹面附近的微泡形变较明显，易产生瞬态空化导致破裂. 同时，微泡具有偏离初始位置向壁面振荡靠近的动力学行为. 在声学参数以及距离壁面底部距离相等的情况下，近刚性凹面下的微泡重心振荡更剧烈，偏离初始位置距离最大，且凹面受到的压力较大，凸面受到的压力相对较小，壁面受到的压力与入射声压呈正相关，平面所受压力和附近微泡重心偏离程度介于凸面和凹面之间. 本模型可为靶向药物治疗等研究提供理论参考.Abstract: In order to study the difference in dynamic behavior of microbubble near rigid plane, convex and concave surfaces, a finite element model of three kinds of microbubble near rigid walls under ultrasonic excitation was established. Results show that the microbubble deformation near the rigid concave surface is more obvious, and it is easy to cause transient cavitation and rupture. At the same time, the microbubble has a dynamic behavior that was deviating from the initial position and oscillating toward the wall. When the acoustic parameters and the distance from the bottom of the wall are equal, the center of gravity of the microbubble under the nearly rigid concave surface oscillates more violently, and the deviation from the initial position is the largest. The pressure on the concave surface is relatively large, the pressure on the convex surface is relatively small, and the pressure on the wall surface is positively correlated with the incident sound pressure, and the deviation between the pressure on the plane and the center of gravity of the nearby microbubble is between the convex and concave surfaces. The proposed model can provide a theoretical reference for targeted drug therapy and other aspects.
-
Key words:
- ultrasound /
- microbubble /
- rigid wall /
- finite element
-
随着对微泡空化效应等特性研究的深入，其在医学临床诊断与治疗上得到广泛应用，尤其是利用微泡在生物壁面产生的生物力学效应. 2020年，国内学者刘晓晖等[1]对肿瘤患者采取低频脉冲超声联合微泡增强肿瘤化疗，其效果较为显著，能够增加阿霉素在肿瘤组织中的转运效率，提高治疗效果. 这与超声场中微泡的非线性动力学行为有关.
目前关于超声场中的微泡空化有两种假设：一是瞬态空化，即超声微泡表现剧烈的非线性振荡和不可控制，并在局部产生高温、高压现象，对细胞形成不可逆性损伤；二是稳态空化，即超声微泡在平衡半径附近呈现近周期性振荡，因而细胞膜的功能性损伤较小[2]. 早在1917年Rayleigh提出了无限大液体中的微泡动力学方程，之后Plesset考虑了液体的黏滞力和表面张力对上述方程进行补充，得到Rayleigh-Plesset方程（RP方程）[3]. 此后各国研究者在RP方程基础上进行了不断修正，但大多数都是在无限大的液体区域中. 1991年，Zudin等[4]推导出有限圆柱形液体区域中的微泡动力学方程，但是此理论研究建立在球形空化泡的假设基础上，实际上微泡在有限区域中振荡并不是呈现球形. 2008年，国内同济大学姜学平等[5]利用镜像法推导出描述不可压缩理想流体中无限大刚性平面附近微泡运动特性的类RP方程并对其进行数值求解,模拟微泡压缩至极小半径前形状的剖面. 2011年，国内邱晓晖等[6]建立了超声激励下微血管内微泡的有限元轴对称模型，得出微管内微泡动力学行为不同于无限大液体中，微泡振荡的频率和有效半径膨胀率随着微管半径减小而抑制. 目前研究大部分基于平面附近和微管内. 事实上，凸形壁面和凹形壁面附近微泡空化在人体组织中并不少见，医学影像显示了日益增长的动脉粥样硬化斑块表面呈现凸形和凹形的壁面[7]，在冲击波聚焦肾结石和膀胱结石的切除过程中，空化微泡常在靶物表面产生许多空洞或凹坑[8]. 因此研究刚性平面、凸面和凹面附近的微泡非线性动力学行为的差异对声学影像和超声溶栓方面的研究有一定的价值[9].
实验研究微泡动力学行为所需相关仪器精度要求较为苛刻，且具有不易观测和成本高等特点. 有限元法计算具有快速、灵活模拟不同条件下微泡的动力学行为的优势. 本文基于生物力学原理，创建了一个气泡—流体耦合模型. 重点研究了不同形状的刚性壁面对微泡动力学响应的影响，并对计算结果进行了分析.
1. 理论模型
本研究在低强度超声压下，研究微泡在不同形状壁面附近的振荡响应. 建立3种情况下的气泡—流体耦合的二维轴对称模型，如图1所示. 其中z轴为模型的对称轴. 模型中相关材料尺度物理量：a为微泡横向半径；
b 为微泡纵向半径；d为微泡重心距离壁面的距离；Rw为凹、凸形壁面半径. 本文研究中d和Rw分别取4 μm和5 μm，且假设初始状态下气泡为球形，即a=b=R0 ，血液为不可压缩流体.1.1 微泡模型
血液和微泡的交界面的初始压力
pg0 为pg0=p0+2σR0−pv (1) 其中：
p0 为微泡内初始气压，文中与初始血液压力相等；pv 为饱和蒸汽压；σ 为气液表面的表面张力；R0 为微泡初始半径.气液界面的初始压力与瞬时压力的关系为
pg−pvpg0=(V0V)γ (2) 其中：
pg 为气液界面瞬时压力；V0 为微泡初始体积；V 为微泡瞬时体积，通过对微泡区域进行积分获得.血液与微泡的耦合边界的关系为
\left( {\left( {{{{p}}_{\rm{g}}} - p} \right)\boldsymbol{I} + \mu \left( {\nabla \cdot\boldsymbol{u} + \nabla\ {\boldsymbol{u}^{\rm{T}}}} \right)} \right) \cdot {\boldsymbol{n}} = {\rm{2}}\sigma k\boldsymbol{n} (3) 式中：
{{p}} 为血液压力；{ \boldsymbol{I}} 为单位张量；n为法向量；\boldsymbol{ u } 为血流速度，气液边界的气泡速度与血液速度相等；\nabla 为哈密顿算符；k 为微泡表面曲率，其值通过弱形式边界偏微分方程（WB）模块对微泡边界进行差分求得.1.2 血液模型
假设血液是不可压缩流体，血流动力学满足Navier-Stokes方程[10]，公式为
\nabla \cdot \boldsymbol{u} = 0 (4) \begin{split} &\rho \left( {\frac{{\partial \boldsymbol{u} }}{{\partial t}} + \left( {\boldsymbol{u} \cdot \nabla } \right)\boldsymbol{u}} \right) =\\ &\qquad\qquad\qquad \nabla \cdot \left( { - p\boldsymbol{I} + \mu \left( {\nabla \boldsymbol{u} + \nabla {\boldsymbol{u} ^{\rm{T}}}} \right)} \right) \end{split} (5) 式中：
\rho 、\mu 分别为血液密度和动力黏度；壁面条件选择为无滑移类型壁面.1.3 超声激励模型
在血液域的边界处设置超声激励条件，表示无限远处的声压. 入射超声满足的关系为
{p_{\rm{u}}} = {p_{{\rm{in}}}}\sin \left( {2{\text{π}} ft} \right) (6) 式中：pu为入射超声；t为时间；
{p}_{{\rm{in}}} 为超声入射声压；f为超声入射频率. 凹面中的微泡振荡较不稳定，在微泡破裂时，考虑有限元计算的收敛性和网格畸变等局限性导致计算终止，故模拟时长为5 μs.1.4 有限元方法
本研究借助于Comsol Multiphysics 5.3a软件中的层流（CFD）模块、动网格（ALE）模块和WB模块建立有限元模型，采用二维轴对称板块建立三维模型来简化计算. 根据有限元法，建立有限元矩阵方程，对整个区域进行网格划分. 使用单位大小的自由三角形网格划分模型：血液域包含2 071个顶点和3 936个三角形单元. 采用MUMPS求解器进行瞬态求解，时间进步选择向后差分法，相对容差0.01，绝对容差0.05，时间步长设置为0.001 μs. 计算并获得时速和时压分布的二维云图，通过选取微泡和壁面边界，计算出每个时刻微泡的形状、位移、微泡附近产生的微射流速度以及附近壁面受到的压力.
2. 模型验证
无限大液体空间中，微泡动力学行为一般用RP方程表示[11]，方程为
\begin{split} & \rho \left( {R\ddot R + 1.5{{\dot R}^2}} \right) = \left( {{P_0} + \dfrac{{2\sigma }}{{{R_0}}} + {P_{\rm{v}}}} \right){\left( {\dfrac{{{R_0}}}{R}} \right)^{3\gamma }} + \\ &\quad {p_{\rm{v}}} - \dfrac{{2\sigma }}{R} - \dfrac{{4\mu \dot R}}{R} - {P_{\rm{u}}} \end{split} (7) 式中：
\dot R 为微泡半径对时间的一阶求导；\ddot R 为微泡半径对时间的二阶求导. 模型计算涉及参数见表1，参数数值的大小与文献[12]一致. 将RP方程计算结果与有限元模型的结果进行对比，如图2所示.表 1 模型参数Table 1. Parameters of geometry model名称 参数 数值 气体多方指数 \gamma 1.070 0 饱和蒸汽压 {P_{\rm{v}}} / Pa 2 330 0 初始半径 {R_0} / μm 2.000 0 气液表面张力系数 \sigma / （N·m−1） 0.056 0 血液密度 \rho / （kg·m−3） 1.059 0 血液动力黏度 \mu / （Pa·s） 0.003 5 超声声压 {P_{{\rm{in}}} } / MPa 0.100 0 超声频率 f / MHz 1.000 0 初始血液压力 { {{p} }_0} / MPa 0.101 3 为模拟无限大液体环境，液体域设置为微泡体积的1 000倍，其计算结果如图2（a）所示，使用RP方程的计算结果如图2（b）所示. 由图2可知，两种方法的曲线趋势相近，验证了模型具有一定可行性. 但也存在差别，主要原因是计算方法和物理场环境的差异. 根据图2（a），微泡横向半径、纵向半径以及有效半径曲线重合，表示微泡在无限血液区域中呈球形振荡.
3. 结果与讨论
超声场下刚性平面、凸面和凹面附近微泡在几个时刻的形状如图3所示. 由图3（a）和图3 (b)可知，微泡在刚性平面和凸面的形变以及产生的微射流无明显区别. 根据图3（c）可知，刚性凹面附近微泡的状态非常不稳定，在t=3.643 μs之后微泡发生了破裂，产生的微射流高达120 m/s. 此外，根据图3所示，3种壁面下的微泡重心位置都竖直向下发生了不同程度的偏离.
超声场下刚性平面、凸面和刚性凹面附近微泡的横向半径和纵向半径随时间的变化曲线如图4所示. 根据图4可知，凹形壁面附近微泡的横向半径的膨胀率明显大于另外两种情况，纵向半径的膨胀率和压缩率也相对较大，在3.643 μs时微泡纵向半径为0，此时已发生了破裂. 图4显示，相对于平面和凸形壁面下的情况，凹形壁面附近微泡的振荡频率较小.
图 4 不同壁面下微泡半径的变化曲线Figure 4. Radius change curves of microbubbles radius under different wall surfaces超声场下刚性平面、凸面和凹面附近微泡的重心位置随时间的变化曲线如图5所示. 由图5可知，3种壁面情况下，微泡在5个正弦周期内发生的最大位移分别为1.07 μm、0.71 μm和3.61 μm，而此现象与文献[13]中得到的两微泡产生的振荡靠近的现象类似，产生的原因与气泡周围的流体速度场变化相关. 声辐射力导致气泡与壁面之间产生交替吸引和排斥的现象[13]. 微泡重心呈振荡靠近壁面，根据图5，在微泡初始位置距离壁面底部相等的情况下，凹形壁面附近微泡重心振幅最大，平面壁面附近微泡次之，凸形壁面附近微泡重心的振幅最小.
图 5 不同壁面下微泡重心位置随时间的变化曲线Figure 5. Chang carves of microbubbles center position with time under different walls超声场下刚性平面、凸面和凹面附近所受最大压力柱状图如图6所示. 根据图6可知，当入射声压为0.1 MPa时，刚性平面、凸面和凹面附近所受最大压力分别为0.65、0.50和1.0 MPa，由于凹形壁面附近的微泡发生了破裂，其产生的微射流对壁面的压力几乎是凸形壁面压力的3倍，是平面所受压力的2倍；当入射声压为0.15 MPa时，3种壁面所受压力几乎是Pin = 0.1 MPa情况的3倍. 由图6所示，凹形壁面对微泡的动力学行为影响较大，凸形壁面对微泡的动力学行为较小.
4. 结 语
本文建立了超声激励下刚性壁面附近微泡的有限元模型. 计算结果表明刚性凹面附近微泡的形变更明显，且凹面附近微泡的膨胀率和压缩率相对较大，微泡的振荡频率较小. 同时，微泡重心具有向壁面振荡靠近的行为. 在声学参数以及距离壁面底部距离相等的情况下，5个正弦周期内近刚性凹形壁面下的微泡重心位置振荡更剧烈并发生破裂，偏离初始位置相对较远，近刚性凸面下的微泡重心位置相对初始位置距离相对较小，微泡状态比较稳定；凹面受到的最大压力最大，凸面受到的压力相对较小. 分析结果表明：3种形状壁面下，凹形壁面对微泡的影响较大，易产生瞬态空化. 然而，本模型未考虑壁面的声反射以及柔性壁面的影响. 因此可在本研究的基础上耦合组织进行生物组织的力学行为分析，为超声治疗方面提供理论指导.
-
图 4 不同壁面下微泡半径的变化曲线
Figure 4. Radius change curves of microbubbles radius under different wall surfaces
图 5 不同壁面下微泡重心位置随时间的变化曲线
Figure 5. Chang carves of microbubbles center position with time under different walls
表 1 模型参数
Table 1. Parameters of geometry model
名称 参数 数值 气体多方指数 \gamma 1.070 0 饱和蒸汽压 {P_{\rm{v}}} / Pa 2 330 0 初始半径 {R_0} / μm 2.000 0 气液表面张力系数 \sigma / （N·m−1） 0.056 0 血液密度 \rho / （kg·m−3） 1.059 0 血液动力黏度 \mu / （Pa·s） 0.003 5 超声声压 {P_{{\rm{in}}} } / MPa 0.100 0 超声频率 f / MHz 1.000 0 初始血液压力 { {{p} }_0} / MPa 0.101 3 
下载: 导出CSV
-
[1] 刘晓晖, 任艳, 韩婷婷, 等. 超声空化联合微泡造影剂增强肿瘤化疗的初步应用[J] . 临床医药文献电子杂志,2020,7(37):54. [2] MOVAHED P, KREIDER W, MAXWELL A D, et al. Cavitation-induced damage of soft materials by focused ultrasound bursts: A fracture-based bubble dynamics model[J] . Journal of the Acoustical Society of America,2016,140(2):1374 − 1386. doi: 10.1121/1.4961364 [3] 江行军, 牛传筱, 吴宇鹏, 等. 低频超声场中微血管内微泡动力学仿真研究[J] . 中国医学物理学杂志,2017,34(2):182 − 187. [4] ZUDIN Y B. Analog of the Rayleigh equation for the problem of bubble dynamics in a tube[J] . Journal of Engineering Physics and Thermophysics,1992,63:672 − 675. doi: 10.1007/BF00853959 [5] 姜学平, 程茜, 钱梦騄. 刚性边界附近微泡运动特性的计算及数值模拟[C]//2008年全国声学学术会议论文集, 上海: 中国声学学会, 2008. [6] 邱晓晖, 沈圆圆, 钱建庭, 等. 刚性微管内微泡动力学行为的有限元数值分析[J] . 生物医学工程学杂志,2011,28(5):911 − 915. [7] JEREMY E, PÁLFI KATALIN, LUISE D F, et al. 3D imaging and quantitative analysis of vascular networks: A comparison of ultramicroscopy and micro-computed tomography[J] . Theranostics,2018,8(8):2117 − 2133. doi: 10.7150/thno.22610 [8] JOHNSEN E, COLONIUS T. Shock-induced collapse of a gas bubble in shockwave lithotripsy[J] . Journal of the Acoustical Society of America,2008,124(4):2011 − 2020. doi: 10.1121/1.2973229 [9] 盛常睿, 陈赛君, 严利明, 等. 注射用六氟化硫微泡造影剂剂量与机械指数对孕鼠胎盘超声造影成像的影响[J] . 中华妇幼临床医学杂志(电子版),2020,16(3):329 − 334. [10] QIN S P, FERRARA K W. Acoustic response of compliable microvessels containing ultrasound contrast agents[J] . Physics in Medicine & Biology,2006,51(20):5065 − 5088. [11] WEISS H L. Mechanical damage from cavitation in high intensity focused ultrasound accelerated thrombolysis[D]. Berkeley: University of California, 2012. [12] 姚文瑛, 许松林, 吴云, 等. 基于计算流体力学的血液和血栓通过静脉瓣时流动分析[J] . 大连理工大学学报,2020,60(4):339 − 348. [13] 蔡晨亮, 屠娟, 郭霞生, 等. 包膜黏弹特性及声驱动参数对相互作用微泡动力学行为的影响[J] . 声学学报,2019,44(4):772 − 779. -
下载:
下载:
点击查看大图
计量
- 文章访问数: 195
- HTML全文浏览量: 263
- PDF下载量: 93
- 被引次数: 0
