水下火山岩似材料凝固过程数值模拟：基于内聚力区模型（CZM）与有限元法（FEM）的热-水-力耦合及柱状节理形成机制研究

《Journal of Volcanology and Geothermal Research》：The solidification of volcanic rock-like materials using cohesive zone models and finite element method

为解决水下火山熔岩冷却凝固过程中热开裂、水渗透与传热耦合机制不清的问题，研究者开展了基于内聚力区模型（CZM）与有限元法（FEM）的2D热-水-力耦合数值模拟，成功捕捉了结壳演化的两阶段特征及柱状节理形成，结果与Kīlauea Iki和Grímsv?tn野外数据吻合良好，为火山岩凝固与相关地质过程研究提供了有效框架。

　　

当地表熔岩流入海洋、湖泊或被降水覆盖时，一场极其复杂的“冰与火之歌”便在水下上演。与空气中缓慢的冷却不同，水的存在会急剧加速熔岩的散热，引发剧烈的热梯度。这种不均匀收缩会在刚刚形成的固态结壳（crust）中产生巨大的热应力，导致岩石破裂，形成独特的网状裂隙（joints）。更关键的是，这些裂隙并非仅仅是破坏，它们会成为水深入高温岩体的通道，引发沸腾换热，进一步改变凝固前沿（solidification front）的推进速度与最终岩石的力学性能。然而，传统纯粹导热的解析模型难以刻画这种热-水-力（Thermo-Hydro-Mechanical, THM）多场耦合的复杂过程，实验研究也因成本高昂和极端环境难以复现而受限。为了揭开这一地质过程的神秘面纱，Léo Falquet等人开展了一项前沿的数值模拟研究，相关成果发表在《Journal of Volcanology and Geothermal Research》。

为了准确回答上述问题，研究人员构建了一个二维（2D）数值模型，巧妙地将内聚力区模型（Cohesive Zone Models, CZM）与有限元法（Finite Element Method, FEM）相结合。该模型不仅模拟了热传导与熔岩相变（液-固），还通过经验定律量化了裂隙几何形态对对流热通量的影响，并最终通过与夏威夷Kīlauea Iki火山口湖和冰岛Grímsv?tn冰盖下火山的实际观测数据进行对比验证。研究成功再现了结壳演化的两个关键阶段，并揭示了水供给强度是控制凝固速率和裂隙间距的核心因素。

作者开展研究用到的主要关键技术方法包括：1. 内聚力-体积有限元法（cohesive-volumetric FEM），在单元间嵌入牵引-分离定律（traction–separation law）以模拟准脆性材料的热开裂与损伤；2. 非光滑接触动力学（NSCD）求解器与θ-方法，处理接触力与内聚力交互；3. 视在热容法（apparent heat capacity method）处理凝固潜热（latent heat），以及有效导热系数法（乘以努塞尔数 Nu）模拟熔岩内的自然对流；4. 基于Jones等（1984）关系的经验沸腾曲线，估算裂隙内水渗入（water ingression）带来的对流换热通量；5. 模拟验证基于Kīlauea Iki和Grímsv?tn的野外观测数据。

研究结果如下：

1. 1.
   结壳因热梯度诱导的开裂（Crust cracking induced by thermal gradients）
   在假设结壳为各向同性弹性脆体、上表面为水饱和温度（100 °C）、下表面为凝固温度（T_sol）的条件下，模拟显示初始热梯度最大，热应力首先达到抗拉强度（σ_max）和临界能量释放率（G_c），触发Mode I开裂。裂隙垂直于冷却面（水平方向）垂直向下扩展。研究发现，需约100 °C的过冷度（相对于T_sol）才会始发开裂，这与Lamur等的实验观察一致。
2. 2.
   水渗入的热通量耗散（Heat flux dissipation by water-ingression）
   通过CZM计算裂隙开度（crack opening, d），结合Jones等（1984）的渗流介质临界热通量关系，推导出单位裂隙长度局部热通量（φ_CHF^local）与d的立方成正比。研究还引入了基于裂隙温度的修正系数k，划分了三个区间：裂隙尖端高温区（T_crack> T_dryout，800 °C）无水接触（k=0）；中间剧烈沸腾区（T_dryout≥ T_crack> T_threshold，130 °C）最大换热（k=1）；低温区（≤130 °C）线性衰减。这表明水渗入换热高度依赖于裂隙开度与壁面温度。
3. 3.
   熔岩凝固（Lava solidification）
   模型将域分为液态（约束自由度，仅热存储）和固态结壳（激活CZM）。当单元温度降至T_sol以下即“活化”为固体。利用视在热容法处理潜热（在T_sol±ΔT内增大C_p），并用Nu=3模拟熔岩内自然对流。模拟显示：早期（约1小时）结壳厚~5 cm，裂隙间距几厘米，但开度仅0.3 mm，不足以显著引水；中期（2小时）结壳厚~10 cm，部分裂隙继续扩展，开度达0.4 mm，水渗入对流与表面传导换热相当；后期（4小时+）仅少数裂隙持续延伸，开度达0.7 mm，对流换热占据主导，最终形成约15 cm厚的传导区（conduction zone），其下为线性推进的凝固前沿。

研究结论与讨论部分指出，该数值框架成功捕捉了水下火山岩似材料凝固的两个阶段：初期以表面热传导主导（结壳厚度随√t演化），后期进入对流与传导平衡态（结壳厚度线性增长），速度与Kīlauea Iki及Grímsv?tn数据吻合（后者因水源充足、~45000 cm/yr，速度~1×10^-5m/s，与模拟的6×10^-6m/s更接近）。模拟所得裂隙间距约7 cm，与Grímsv?tn观测的10-20 cm接近。参数分析表明，热膨胀系数（α）和凝固温度是关键控制因子，低于~2×10^-6K^-1则无裂隙，仅为传导冷却。该工作建立的验证框架可扩展至多种几何与条件，对理解火山灾害、地热储层（裂隙网络渗透性）及地质工程稳定性具有重要意义。未来方向包括向3D扩展以还原真实六方柱状节理、引入温度依赖材料属性及近T_sol的玻璃化转变流变行为。