Recently, there has been significant progress in the generation of ultrahigh pressures and temperatures through the use of new hard anvils and heater materials in large-volume presses (LVPs). Pressures of over 501 and 120 GPa2 have been achieved using hard tungsten carbide (WC) and sintered diamond (SD) cubes, respectively, as second-stage anvils in an LVP. For ultrahigh-temperature generation, boron-doped diamond (BDD) has been demonstrated to be a good heating material owing to its extremely high hardness, thermal conductivity, electrical conductivity, and other properties. With the use of BDD as a heater, temperatures of ∼3300 K at 30 GPa3 and ∼4000 K at 15 GPa4 have been reached using WC anvils in an LVP. Yoneda et al.5 generated a temperature of 2723 K at a pressure of around 50 GPa using harder SD anvils. More recently, an ultrahigh pressure and temperature of 50 GPa and above 3300 K have been achieved simultaneously using SD anvils and a BDD heater in combination with in situ synchrotron X-ray diffraction.6 The same study also demonstrated that thermal pressure could significantly enhance high-pressure generation inside a BDD chamber.6 Furthermore, the development of CVD-BDD and machinable BDD heaters,7,8 which are convenient for normal LVP users, has opened up new opportunities for melting or viscosity experiments under the high pressures and temperatures corresponding to those of the Earth’s lower mantle.

Unfortunately, the application of BDD heaters in ultrahigh-temperature–pressure studies has been limited owing to the lack of an appropriate method to estimate temperatures above the highest applicable temperature (∼2600 K) of thermocouples and the temperature distribution within the sample area. Currently, temperatures above 2600 K are usually estimated through extrapolation of a polynomial power–temperature relation below 2600 K or through melting points of refractory materials.9 The values estimated by polynomial extrapolations are strongly affected by the order of the polynomial, which causes a large uncertainty. Because of the significant thermal pressure inside a BDD heater, the uncertainties in the pressure in offline experiments reach 10 GPa or even larger, depending on the temperature,6 which in turn affects the melting points of materials. For example, the melting point of WC rises from 3140 K at ambient pressure to 3800 K at about 9 GPa.10 As a result, the method using melting points of refractory materials may induce additional uncertainties when applied to BDD heaters. For the temperature distribution below 2600 K within a BDD heater, an in situ synchrotron X-ray diffraction method has been used to measure the axial temperature gradient at high pressure.3,11 However, this method cannot be applied in offline experiments.

A finite element simulation (FES) model might be a good method to estimate temperatures beyond 2600 K for assemblies using BDD heaters. Yoneda et al.5 used a two-dimensional (2D) FES with direct current (DC) loading to estimate the temperature generation and distribution inside the chamber of a BDD heater.12 However, their model did not consider three-dimensional (3D) effects, the presence of electrodes, thermocouples, and second-stage anvils, or the commonly adopted alternative current (AC) loading method for heating. Therefore, there is a need for a more reliable 3D model that is based on a realistic cell assembly and is validated by real experiments.

In this study, we develop a 3D FES model as a virtual thermometer to estimate temperatures above 2600 K and assess the temperature homogeneity inside a BDD heater. The virtual thermometer is validated by ultrahigh-temperature experiments at loads of 2.8–7.9 MN (corresponding to pressures of 19–28 GPa). Using this virtual thermometer, we estimate the temperature generation and distribution above 2600 K inside BDD heaters for ultrahigh-temperature–pressure experiments.

Before performing the FES, we conducted experiments at high temperatures and pressures to obtain parameters for the simulation, such as the geometry of the cell assembly and the conductivity of the BDD. All the experiments were conducted using the JLUHC 1000-ton LVP (Walker-type, maximum load 10 MN) installed at the State Key Laboratory of Superhard Materials (SKLSHM), Jilin University. Figure 1 shows a schematic of the internal structure of the Walker-type LVP module. The specific experimental principles and details of the JLUHC 1000-ton LVP have been extensively discussed in our previous paper.13

Figure 2 show schematics of the BDD cell assemblies used in this study, which were modified from those reported by Xie et al.4 We used ZK01F WC cubes with 3 mm truncated edge length as second-stage anvils and Cr₂O₃-doped MgO octahedra with 7 mm edge length as the pressure media. BDD tubes (Changsha 3 Better Ultra-Hard Materials Co., Ltd., China) were used as heaters. Photographs and microstructures of these BDD tubes are shown in Fig. S1 (supplementary material). The concentration of boron (B) in the BDD heater is ∼4692 ppm, which corresponds to a mass percentage of 0.422 wt. % from secondary ion mass spectrometry (SIMS) measurements [Fig. S1(d)]. The content of other elements such as H, O, and N is below 0.01 wt. %. We determined the electric conductivity of the BDD blocks using a multifunctional measurement system via the van der Pauw method at room temperature and atmospheric pressure. The measured conductivity (506 S/m) is slightly higher than that (464 S/m) in a previous study.5 Since the BDD heater in that study was converted from a boron-doped graphite heater in situ, a small amount of diamond whithout doping may have been present in their BDD heater owing to incomplete doping during graphite–diamond conversion,5 which may be the reason for the conductivity difference. Mo or TiC powder was selected as the electrode material, as suggested by Xie et al.4 Temperatures were measured by a D-type thermocouple (TC, W₉₇Re₃–W₇₅Re₂₅), whose hot junction was put at the center of the BDD tube. For most of the experiments, a tube of Ca-doped ZrO₂ was inserted between the BDD heater and the (Mg, Cr)O octahedron as a thermal insulator [Fig. 2(c)]. Powdered MgO or pyropic glass (En₇₅Cor₂₅, where En = MgSiO₃ enstatite, Cor = Al₂O₃ corundum, and the subscripts indicate mol%) samples were loaded into the BDD heater. Details of the glass preparation have been presented in our earlier paper.14 To remove moisture, all ceramic components were annealed at 1100 K for 2 h, except for the TiC powder, which was dried in a vacuum oven at 600 K for 2 h.

The cell assembly was first compressed to the target load at room temperature and then heated to the target temperature, which was maintained for around 3 min. After that, the experiment was quenched by turning off the electric power, and the pressure was released slowly over 10 h. The generated pressure reported in the present study is estimated based on the oil pressure-cell pressure calibration curve (Fig. S2, supplementary material).13 Owing to the extreme hardness of BDD, significant thermal pressure is expected inside a BDD heater.6 Therefore, the actual pressure might be higher than the reported pressure. An AC power source was adopted for heating (for details, see Fig. 3 in Sec. II B). After their recovery, the assemblies were embedded in epoxy resin and then ground and polished with sandpaper and polishing cloths using diamond paste. Phase assemblages were characterized by Raman spectroscopy (InVia Raman spectrometer equipped with a 514.5-nm laser, Renishaw Company) and a scanning electron microscope (SEM) combined with energy-dispersive spectroscopy (EDS, Magellan 400, Field Electron and Ion Company). The surfaces of the polished samples were coated with Au before SEM characterization.

We constructed digital twins of the high-temperature-pressure cell assemblies, which served as virtual thermometers for evaluating the sample temperature and its distribution. During Joule heating, on the one hand, the applied electrical current determines the heat source for the increasing temperature. On the other hand, the electrical conductivities of the parts of the cell assembly are influenced by the temperature as well. Thus, thermal–electrical coupling occurs during the Joule heating, which can be simulated by a finite element method. We simulated this coupling process using the commercial finite element software Abaqus (ver. 2021, SIMULIA, Dassault Systems).

The relationship between electric potential, material conductivity, and current density in the assemblies is expressed by the following equation:∫V∂δφ∂x⋅σE⋅∂φ∂xdV=∫SδφJ⋅dS+∫VδφrcdV,(1)where φ is the electric potential, δφ is a commonly used test function in the finite element method, σ is the electrical conductivity tensor, J is the electric current density (current per unit area), and r_c is the internal current source per unit volume.

The heating power P_ec generated by the current flow can be expressed as15Pec=J⋅E=∂φ∂x⋅σE⋅∂φ∂x,(2)where E is the electric field, defined as the negative of the gradient of the electric potential φ.

For the heat conduction equations, a steady-state temperature field will be established for a DC loading mode. The steady-state temperature field satisfies the energy conservation equation, expressed in variational form as follows:15∫VρU̇δθdV+∫V∂δθ∂x⋅k⋅∂θ∂xdV=∫VδθrdV+∫Sδθq⋅dS,(3)where V is a volume of solid material with surface S, ρ is the density of the material, U is the internal energy, θ is the temperature, δθ is the test function, k is the thermal conductivity tensor, q is the heat flux per unit area of the body, flowing into the body, and r is the heat generated within the body.

For AC loading, the temperature field is governed by the transient heat transfer equation expressed as follows:ρcdθdt=∇⋅(k⋅∇θ)+q,(4)where c is the specific heat capacity of the material, t is the time, and θ is the temperature field as a function of time and space.

The phase angle control method is usually employed to regulate power within the AC voltage loading mode in a real experiment (Fig. 3). The voltage can be expressed as a piecewise function:V(t)=Vm⁡sin⁡ωt,nT+t0≤t≤nT+12T or nT+12T+t0≤t≤nT+T,0,otherwise,(5)where V_m is the amplitude of the AC voltage, ω is the angular frequency, t is the time, n is a natural number, T is the period of the AC power, and t₀ is the moment at which the thyristor conducts in one period. In a complete cycle, the thyristor is only conductive during the time intervals from t₀ to 12T and from t0+12T to T. The voltage V(t) is applied to the model as a boundary conditions in the AC loading mode. A frequency of 50 Hz is adopted in this study.

Instead of the real voltage, an equivalent voltage, defined as the voltage that does the same work within one cycle as the actual voltage, is commonly reported in real experiments. Thus, the equivalent voltage satisfies the following equation:∫0TVe2Rdt=∫t012TVm2⁡sin2⁡ωtRdt+∫12T+t0TVm2⁡sin2⁡ωtRdt,(6)where V_e is the equivalent voltage and R is the resistance value of the circuit.

Table I summarizes the material parameters used in this study. Both temperature and pressure may affect the material properties. The temperature effect is expressed by using polynomial functions as shown in Table I, which were obtained by fitting experimental data from the literature. Owing to a lack of high-pressure data on material properties, the effects of pressure on the material parameters was ignored in the simulations. In particular, although ZrO₂ undergoes multiple phase transitions under high-temperature and high-pressure conditions,16–19 the thermal conductivity of ZrO₂ under ambient conditions was adopted for all the pressures in the present study.

We obtained values of the electrical conductivity of BDD by directly measuring the resistance of the assembly in situ and the geometrical dimensions of the BDD heater in the recovered samples (Fig. S3, supplementary material). Given that the resistances of WC and TiC are more than three orders of magnitude lower than that of BDD in the assembly, the whole of the resistance of the assembly was attributed to the resistance of BDD in the calculation. The conductivities of BDD at various temperatures are provided in the heating logs (Tables S1–S5, supplementary material). To ensure the accuracy of the simulation, the unprocessed, raw electrical conductivity data were used directly.

Figure 4 shows the geometry and boundary conditions of the FES. We simplified the outer octahedral MgO pressure medium into a cylindrical shape with a height of 3.6 mm. To check the effect of the cylinder diameter on the computational results, we performed a series of testing simulations for various diameters. We found that the computational results were sensitive to the diameter size if the diameter was smaller than 2.075 mm, whereas the results converged if the diameter was larger than 2.075 mm. To determine the diameter for our model, we estimated the octahedron size under high temperature and pressure on the basis of recovered assemblage [Fig. S3(a), supplementary material]. The recovered assemblage had a complex shape, which was composed of 12 extruded fin-shaped MgO gaskets and an MgO octahedron truncated with edges and vertices. The thickness of the gasket near the octahedron was ∼0.8 mm [Fig. S3(a)]. The edge length of the octahedron without truncation should be 4.7 mm, the inscribed circle of whose triangular face has a diameter (2.7 mm) significantly larger than 2.075 mm. To save computational resources without sacrificing simulation accuracy, we adopted a cylinder diameter of 2.075 mm in our model. The WC anvils were also simplified as cylinders with a height of 1 mm and a radius of 2.075 mm at the top and bottom ends of the MgO pressure medium. For the model aiming to investigate the effect of the thermocouple on the temperature distribution, a D-type W–Re thermocouple with the same diameter (0.075 mm) as that used in the real experiments was set at the center of the assembly [as shown in Fig. 4(c)]. The shapes and dimensions of all the other internal components were the same as those used in the real experiments [as shown in Fig. S3(c)]. On the basis of the symmetry of the assembly, the 3D models with and without thermocouple were simplified to 1/8 and 1/4 assembly, respectively.

For the boundary conditions, the mode’s outer surface temperature was set to a constant room temperature of 298 K (yellow dashed lines in Fig. 4) and the potential of the bottom WC cylinder was set to zero (purple dashed lines). The loading voltage (red dashed line) was applied to the top WC cylinder for the 2D model, while half the loading voltage was applied to the upper surface of the BDD heater (areas within the pink dashed lines) for the 3D models. The initial temperature was set to 298 K for all the models.

A series of tests were conducted to refine the mesh size and increment step size. As shown in Fig. S8 (supplementary material), the computational results converged after the mesh size and increment step size were finer than 3 × 10⁻⁵ m and 1 × 10⁻⁴ s, respectively. To balance computational time and accuracy, the mesh size inside the heaters and the increment step size were set to 3 × 10⁻⁵ m and 1 × 10⁻⁴ s, respectively.

Table II summarizes the present experimental conditions and results. A series of experiments were conducted to assess the effects of electrode materials and thermal insulating materials on temperature generations. The first experiment [JLUC231, using the configuration shown in Fig. 2(a)] was conducted to determine the feasibility of transplanting the BDD heater into the Walker-type LVP. To simplify the construction procedure of the cell assembly, we adopted Mo rods as electrodes. In the subsequent experiment [JLUC263, using the configuration shown in Fig. 2(b)], we replaced the Mo electrodes with TiC powder, as recommended by Xie et al.,4,7,9 to increase the upper temperature limit of the cell assembly. In the experiment JLUC272 [using the configurations shown in Fig. 2(c)], we introduced a thermal insulation layer of ZrO₂ between the pressure medium and the BDD heater to increase the efficiency of temperature generation.

Figure 5 shows the relationship between temperature and heating power for some representative experiments under press loads of 2.8–7.9 MN, corresponding to pressures of 19–28 GPa on the basis of our recent pressure calibrations.26 Without the use of TiC electrodes and a ZrO₂ insulation layer, the power–temperature relationship is linear below 1500 K (JLUC231). As the temperature continues to rise, the heating power becomes unstable. Eventually, there is a sudden increase in resistance followed by a dramatic drop in power at 582 W, leading to termination of heating. The maximum temperature estimated from a cubic polynomial extrapolation (∼2850 K) is ∼600 K higher than that estimated from a linear extrapolation (∼2250 K). From cross-sectional SEM images of the recovered sample (Fig. S4, supplementary material), it can be seen that the Mo electrodes had undergone melting and entered the interior of the BDD cylinder, which may have caused the heating termination. Since the power–temperature relation became abnormal from ∼2200 K, far below the melting point of Mo (2900 K),32 it is likely that the temperature of the Mo electrode in the absence of a thermal insulator was sufficiently high for Mo to be oxidized into MoO₃. The MoO₃, with a melting point of 1100 K under ambient pressure, could then act as a melting flux and cause the abnormal melting of the Mo electrode. On replacing the Mo electrode with TiC powder in the run JLUC263, the heating efficiency became slightly higher than that in the run JLUC231. The maximum temperature estimated through a cubic polynomial extrapolation (3055 K) was ∼670 K higher than that estimated through a linear extrapolation (2380 K). We also observed frequent power fluctuations and blowouts in this cell assembly after heating to ∼2300 K. These may have been caused by the omission of thermal insulation in this experiment, which resulted in a rapid increase of temperature and softening of the pressure medium and gaskets with heating. After heating to ∼2300 K, the materials become too soft to seal the pressure inside, resulting in blowouts.6

The introduction of the ZrO₂ thermal insulation tube significantly enhances the heating efficiency (JLUC272). The temperature and power show an excellent linear relationship below 2600 K, indicating an improved stability of the assembly. Note that phase transitions of ZrO₂ should occur above ∼1300 K.33 Considering the temperature gradient, the sample temperature could be several hundred kelvin higher than the temperature of ZrO₂. The excellent linear relationship without any abnormal power–temperature relationship between 1300 and 2000 K suggests that the phase transitions of ZrO₂ should have little effect on its thermal conductivity. The maximum temperature possibly reached ∼3500 K at a heating power of 470 W from a linear extrapolation. The runs of JLUC300 and JLUC308 adopted the same assembly as JLUC272 to verify the repeatability of the experiment.

Figure 6 displays back-scattered electron (BSE) images of the recovered samples in the runs JLUC308 [Fig. 6(a)] and JLUC272 [Fig. 6(b)], both of which were subjected to temperatures above 3400 K. Both samples exhibit a hemispherical mixed phase at the bottom of the BDD heater and dendrites texture throughout the heater chamber [Figs. 6(c) and 6(d)]. Chemical analysis showed that these phases contain not only elements of Mg, Al, and Si from the starting materials, but also elements of Zr, Ti, Ca, W and Re from surrounding materials (Table S6, supplementary material). These observations indicate that eutectic melting occurred among the starting material, ZrO₂, TiC, and thermocouple. The melts moved to the bottom of the BDD heater owing to the melt–solid density contrast and formed a hemispherical mixed phase at the bottom of the heater. The BSE images of the other recovered samples are shown in Fig. S5 (supplementary material). The distribution of elements in these samples is similar to that shown in Figs. 6(c) and 6(d). Their run products were divided mainly into two types of Mg/Si/Al-rich (melt 1) and Zr-rich (melt 2) quenched crystals from melts with different compositions (Table S6, supplementary material).

After we had obtained the detailed parameters of heating experiments under various pressures using BDD heaters, we constructed an FES model to serve as a virtual thermometer to verify the temperature reached from our experiments. First, we examined the FES model within the measurement range of the thermocouple (<2600 K). Two model dimensions (2D and 3D) and two voltage-loading methods (DC and AC) were considered. In total, four models (2D&DC, 2D&AC, 3D&DC, and 3D&AC) were tested. Let us take the temperature generation of 2473 K at a voltage of 3.69 V in run JLUC306 as an example. Figure 7 shows the temperature evolution of the sample center in the four models. Temperature stabilizes within 2 s in all four models [Fig. 7(a)]. Figure 7(b) shows that the temperature fluctuation caused by AC loading is ±2 K, which is much smaller than the observed experimental temperature fluctuation (±14 K, Fig. S6, supplementary material). The remaining temperature variations will be due to other random factors such as electrical noise in the thermocouple caused by the heater. The 2D models systematically yield a higher temperature (∼10 K) than the 3D models, which may be due to mathematical differences between the models or the smaller heat-conducting surface in the 2D models compared with the 3D models. The AC-loading models systematically yield a temperature higher by ∼40 K compared with the DC-loading models, which may be caused by the difference in the resistance of the assembly between the AC- and DC-loading models. Figure 8(a) summarizes the temperature difference (experimental value minus simulated value) of the present four models for run JLUC306. The 3D&AC model gives the closest values to the experimental results. Therefore, this model was selected as a virtual thermometer in this study.

The uncertainty of the virtual thermometer was assessed by comparing the temperature difference between our predictions and the measured experimental values. The possible error sources are the boundary conditions, the mesh method, the increment step of the simulations, and variations in the thermal–electrical behaviors of cell components due to phase transitions, pressure, and temperature changes. As described in Sec. II B 3, we minimized the effects of boundary conditions, mesh density and increment step. Figure 8(b) shows the temperature difference as a function of temperature. The temperature difference increases roughly linearly with increasing temperature. Since most of the material parameters at high temperatures were extrapolated from low-temperature data, this trend is probably caused by the increased errors in the material parameters with increasing temperature. Note that no significant increases in temperature difference were observed within the temperature range (1300–2000 K) of the ZrO₂ phase transitions, which indicates the small effect of these phase transitions on the simulation results and is consistent with the experimental results. The average simulation error at 3500 K is ∼80 K, which is significantly smaller than that from the heating power–temperature extrapolation method (∼600 K). Figure 8(c) shows the temperature difference as a function of the press load. No significant relationship can be observed, indicating that the thermometer is insensitive to pressure. This insensitivity to pressure is confirmed by the good match between the simulated and experimental temperature–power relations at different pressures (Fig. 9). Because the pressure inside the BDD heater increases significantly with temperature, the insensitivity to pressure makes the virtual thermometer suitable as a thermometer for this heater. Given that significant changes occur in material properties during melting, the lowest melting temperature of the cell assembly components represents the highest temperature at which our virtual thermometer can be used.

Figure 9 compares the temperature–power relationship obtained from high-pressure-temperature experiments and FES. The 3D&AC FES reproduces the temperature–power relationships well for high-pressure–temperature experiments below 2600 K, confirming the validity of the thermometer. The temperature–power relationship can be linear or nonlinear, depending on the cell assembly. Special caution is required when estimating temperature using a temperature–power relationship. According to the experiments, melting occurs at ∼3500 K owing to the eutectic melting of the assembly at a load of 7.9 MN. The lowest estimated temperature (3500 K) of the melting experiments represents the highest temperature that can be measured by our thermometer.

Using the thermometer, we also examined the temperature distribution within the BDD heaters in our assemblies. Figure 10 presents simulated temperature contour maps of run JLUC306 at ∼3500 K. Figure 11 shows the temperature variations along the radial (X) and longitudinal (Y) directions. The central position of the sample was set as 0 on both the X and Y axes. When no thermocouple is used in the assembly, the center temperature of the sample is 6–8 K lower than the inner surface of the BDD tube. The temperature gradient of the sample along the X direction is 19–26 K/mm, and that along the Y (longitudinal) direction is 83 K/mm. The temperature gradient inside the BDD heaters along both the radial and longitudinal directions is lower than those inside traditional heaters such as graphite or LaCrO₃ (100–200 K/mm).34,35 The larger temperature gradient in the longitudinal direction than the radial direction indicates that heat is primarily dissipated in the longitudinal direction, which will cause the temperature of the electrode anvil top to be higher than those of the other anvils’ tops. This could be why the top of the electrode anvil, rather than the tops of the other anvils, sometimes partially melts during the generation of ultrahigh temperatures (Fig. S7, supplementary material). When a thermocouple is used, the temperature at the sample center is reduced by ∼80 K [Fig. 10(c)], and the temperature gradient is also reduced to 6.55 K/mm within the heater (Fig. 11), owing to the high thermal conductivity of the thermocouple [180 and 80 W/(m K)].24 In contrast, the thermocouple transfers part of the Joule heat outside the heater and raises the temperature there.

We have established a 3D AC-loading FES, serving as a virtual thermometer to estimate the temperature and analyze the temperature field above 2600 K in a high-temperature–pressure assembly. This virtual thermometer gives results consistent with those of a thermocouple for temperatures below 2600 K, confirming its validity, and gives temperature predictions with an error smaller than 80 K. It is also insensitive to pressure, making it a good thermometer for experiments with a BDD heater. Simulation results indicate that the extremely high thermal conductivity of the BDD heater contributes to a highly uniform temperature distribution within the heater (19–26 K/mm in the radial direction and 83 K/mm in the longitudinal direction). Our virtual thermometer provides a good temperature reference above 2600 K for experiments in LVPs, enabling investigations of the melting properties of minerals and materials in both Earth and material sciences.

Acknowledgment. The authors thank Dr. Donghui Yue for his guidance in the finite element simulation work. This work is supported financially by the National Key R&D Program of China (Grant No. 2022YFB3706602), the National Natural Science Foundation of China (Grant Nos. 42272041, 41902034, and 12011530063), and the Jilin University High-Level Innovation Team Foundation, China (Grant No. 2021TD-05).