Numerical study of phase-change heat transfer in multilayer wellbores under moving heat sources

1. Introduction

Plug and Abandonment (P&A) is a critical stage in the decommissioning of oil and gas wells, aiming to ensure long-term wellbore integrity and prevent the migration of hydrocarbons or other contaminants. Previous studies have shown that abandoned wells may act as potential leakage pathways for natural gas or oil, posing environmental risks [1]. Conventional P&A operations typically rely on cement plugs to establish multiple sealing barriers. However, cement materials are prone to the formation of microcracks and interfacial micro-annuli during long-term service, which may compromise well integrity and lead to leakage risks [2]. Therefore, the development of more reliable and durable alternative sealing technologies is of significant engineering importance.

In recent years, downhole thermal plugging technologies based on aluminothermic reactions and low-melting-point alloys have attracted increasing attention. Thermite reactions serve as high-temperature heat sources capable of releasing substantial thermal energy downhole to heat and melt the alloy. The molten metal subsequently flows under gravity into fractures and pores within the wellbore and solidifies upon cooling, thereby forming a metallic sealing barrier.

Considerable efforts have been devoted to understanding heat transfer and phase change phenomena involved in thermal plugging. Existing studies have addressed combustion propagation mechanisms, heat conduction analysis, numerical simulation, and process optimization. Durães et al. [3] investigated radial combustion propagation in iron oxide/aluminum thermite mixtures. De Andrade et al. [4] proposed a hybrid analytical method for transient heat conduction in composite structures. da Silva et al. [5] employed an inverse heat conduction approach to estimate thermal input from thermite reactions, reducing modeling complexity. Abdelal et al. [6] developed a nonlinear finite element model to analyze the influence of bismuth solidification expansion on wellbore wall strength. Pena et al. [7-8] constructed thermal models of aluminothermic reactions, revealing the influence of pore structure on heat transfer efficiency and simulating heat conduction and phase change within well structures while accounting for thermal input to casing steel. Tobisawa and de Lemos [9] proposed a method for transporting thermite mixtures in fluid media with particle separation through filtration. Carpenter [10] reported a current-controlled thermal activation process that enhances heat transfer efficiency. dos Santos Magalhães and de Lemos [11] simulated temperature distributions during thermite-based P&A operations, incorporating the latent heat of pipe melting into the heat diffusion equation and modeling the problem in a three-dimensional domain.

These studies have provided valuable foundations for thermal field modeling and material design in downhole thermal plugging. However, systematic numerical investigations of the melting behavior of low-melting-point alloys under moving heat source conditions remain limited, particularly with respect to the quantitative prediction of phase interface evolution and complete melting time.

In this study, an axisymmetric transient heat transfer-phase change coupled finite element model is developed for a low-melting-point alloy subjected to a moving heat source. This work aims to characterize the transient temperature field, heat-transfer pathway, and phase-fraction evolution, and to quantitatively determine the complete melting time under the specified thermal plugging conditions. The results provide a theoretical basis for parameter design in downhole thermite-based P&A operations and for further investigation of the subsequent alloy solidification process.

2. Research method

3. Axisymmetric finite element model

An axisymmetric two-dimensional transient heat transfer finite element model is developed to investigate the melting behavior of low-melting-point alloys during downhole thermal plugging. Taking the 50 kg alloy case as representative, a multilayer composite wellbore structure with a total model height of 𝐻 = 1.5 m is constructed, including the thermite outer wall, phase-change alloy annulus, Tubing, annular water layer, casing and surrounding rock.

Under the axisymmetric assumption, the wellbore is simplified as an $r$-$z$ plane problem, with the heat source surface located 30 mm from the axis of symmetry. An axisymmetric boundary condition is imposed along the central axis, while the outer boundary is placed sufficiently far from the heat source to minimize the influence of the artificial truncated boundary. The model configuration is shown in Fig. 2.

Fig. 1Enthalpy-temperature curve of the alloy

Fig. 2Schematic diagram of the multilayer finite element model

Spatial discretization is performed using the two-dimensional thermal element PLANE55, with mesh size determined based on the outer pipe wall thickness. A mesh-independence check confirmed that the adopted mesh is sufficient for resolving the phase-fraction evolution and temperature distribution. Thermophysical properties are assigned by region. The alloy is implemented via the enthalpy-temperature relationship to account for latent heat, whereas constant properties are used for water, steel, and rock (Table 1). The initial temperature is set to 50 °C to represent downhole conditions.

The heat source is described by an idealized three-stage time-dependent function: linear increase to the peak heat flux during 0-100 s (ignition stage), linear decrease to zero during 100-200 s (attenuation stage), and zero heat flux during 200-500 s (cooling stage). The peak heat flux density is $q_{peak}$ = 1×10⁶W·m^-2. The heat source moves axially at a constant velocity of 0.01 m/s, simulating combustion front propagation along the wellbore. The time scale and piecewise-linear variation were inspired by the heat-flux characteristics reported by da Silva et al. [5] for aluminothermic heat release, while simplified for the present numerical analysis.

4. Results and discussion

To characterize the melting behavior of the alloy under thermal plugging conditions, a 50 kg alloy configuration was selected as the baseline case. Under an environmental temperature of 50 °C, the heat-transfer characteristics and phase-change evolution were analyzed. Due to the presence of the moving heat source, heat transfer within the confined wellbore exhibits a characteristic pattern of “upper-region priority” followed by axial downward propagation.

Fig. 3Heat flux vector distributions at 100 s and 200 s

The heat flux vector distributions at 100 s and 200 s (Fig. 3) show that thermal energy is concentrated near the heat source, forming a high-temperature-gradient region. Heat propagates downward along the axial direction while simultaneously conducting outward through the tool wall into the annular water layer and further into the casing and surrounding rock. This heat transfer pathway leads to a spatially non-uniform melting process, where phase change does not occur simultaneously but progresses along the axial direction following the movement of the heat source.

The evolution of the temperature field further confirms this behavior. Temperature contours within the phase-change alloy at representative time instants (10 s, 50 s, 150 s, and 225 s), shown in Fig. 4, illustrate the staged melting process. During the early stage (10 s), the upper portion of the alloy first enters the phase transition range and forms a localized high-temperature zone. At 50 s and 150 s, the isotherms shift downward and the phase-change region progressively expands, indicating axial propagation of the melting front. At 225 s, the temperature distribution becomes relatively uniform, and the alloy temperature exceeds the liquidus temperature throughout the domain.

The time evolution of the phase fraction is presented in Fig. 5. During the initial 0-50 s, the phase fraction increases slowly, indicating that this stage is dominated by sensible heat accumulation and temperature gradient development. From 50 s to 200 s, the phase fraction increases approximately linearly, suggesting that the melting front propagates at a nearly constant rate under conduction-controlled conditions. After 200 s, the increase in phase fraction slows down, and the melting rate gradually decreases until complete melting is achieved at approximately 225 s.

Fig. 4Temperature contours of the phase-change alloy at different time instants (10 s, 50 s, 150 s, and 225 s)

Fig. 5Time evolution of phase fraction

5. Conclusions

An axisymmetric two-dimensional transient heat transfer-phase change coupled finite element model was developed to investigate the melting behavior of a low-melting-point alloy under a moving heat source.

The simulation results show that the moving heat source induces axial downward propagation of the melting region within the confined wellbore, resulting in a spatially non-uniform evolution of the phase interface. For the 50 kg alloy case under the specified thermal plugging conditions, the melting process can be divided into three distinct stages: an initial stage dominated by sensible heating before 50 s, a quasi-steady melting stage from 50 s to 200 s, and a final slowing stage after 200 s. Complete melting was achieved at approximately 225 s.

These results clarify the alloy melting behavior in the plugging tool and provide a numerical basis for the design of thermite-based P&A operations. In practical P&A applications, the subsequent cooling and solidification process is also critical for forming a stable metallic barrier. Therefore, future work will extend the present model to the cooling and solidification stages to further evaluate the formation and integrity of the final metallic barrier.