Verification of an Explicitly Coupled Thermal-Phase Field-Mechanical Electromagnetic (TPME) Framework by the Method of Manufactured Solutions

Travis S. Ramsay

Department of Mechanical Engineering, Iowa State University, Ames, IA, USA.

DOI: 10.4236/ojmsi.2021.91001                

Abstract

An explicitly coupled two-dimensional (2D) multiphysics finite element method (FEM) framework comprised of thermal, phase field, mechanical and electromagnetic (TPME) equations was developed to simulate the conversion of solid kerogen in oil shale to liquid oil through in-situ pyrolysis by radio frequency heating. Radio frequency heating as a method of in-situ pyrolysis represents a tenable enhanced oil recovery method, whereby an applied electrical potential difference across a target oil shale formation is converted to thermal energy, heating the oil shale and causing it to liquify to become liquid oil. A number of in-situ pyrolysis methods are reviewed but the focus of this work is on the verification of the TPME numerical framework to model radio frequency heating as a potential dielectric heating process for enhanced oil recovery. Very few studies exist which describe production from oil shale; furthermore, there are none that specifically address the verification of numerical models describing radio frequency heating. As a result, the Method of Manufactured Solutions (MMS) was used as an analytical verification method of the developed numerical code. Results show that the multiphysics finite element framework was adequately modeled enabling the simulation of kerogen conversion to oil as a part of the analysis of a TPME numerical model.

Keywords

Radio Frequency Heating, In-Situ Pyrolysis, Oil Shale, Multiphysics, Explicit Coupling, Finite Element Method, Method of Manufactured Solutions

1. Introduction

The term, “oil shale” itself is not geologically defined by a specific chemical formula but in general refers to fine-grained sedimentary rocks that yield shale oil upon pyrolysis or retort [1]. Kerogen, defined as an insoluble macromolecular organic matter, is also recognized as the most abundant form of organic matter on Earth [2]. The organic kerogen represents the premature form of petroleum that has not been exposed to sufficiently high temperature and pressure over a couple of millions of years to be converted to shale oil. Given the contemporaneous existence of kerogen in oil shale formations, it is understood that the natural hydrocarbon maturation process has yet to occur or is in the early stages of development. It is possible to imitate the natural maturation process through in-situ pyrolysis or surface retort so that oil and gas can be generated with the intent of production [3]. In-situ pyrolysis involves heating the oil shale directly in the subsurface, thus can be performed for shallow and deep formations. The applicability of in-situ pyrolysis across shallow and deep formations gives it a significant advantage [4]. A significant disadvantage of in-situ pyrolysis; however, is that its higher technical burden necessitates increased spending on advanced technologies which can only be achieved during times of increased oil prices [4]. Conversely, surface retort involves mining the oil shale, bringing it to the surface, performing surface pyrolysis followed by additional processing. Surface retort, as a result, is only economically feasible when the oil shale formations are shallow, with respect to Earth’s surface [4]. The advantage of surface retort is that since it can be conducted by mining and heating on the surface, the associated processing does not require highly sophisticated technology. In both cases once the temperatures reach ~700˚F (~370˚C), depending on the specific oil shale formation, chemical decomposition occurs which leads to the conversion of kerogen (solid) to oil (liquid) [4] [5].

Several technologies have been investigated in order to evaluate the production potential of oil shale. Given oil price uncertainty and increased U.S. desire for energy independence in the last couple of decades, the need for more reliable yet economic oil production techniques has increased. Challenges associated with production from oil shale include water use, net energy usage, carbon dioxide emissions, and commercial scalability [6], which if adequately addressed may lead to oil shale contributing to future supplies of shale oil. While radio frequency heating has been applied to heavy oil production [7] [8] [9], few have evaluated its usage for in-situ pyrolysis of oil shale [10] [11] [12]. Furthermore, the development and publication of numerical models that honor the key physical processes associated with the kerogen to oil conversion have been minimally disclosed in literature.

The goal of this study is to demonstrate the verification of an explicitly coupled 2D TPME code that was developed using a general-purpose finite element framework, leveraging the TalyFEM libraries, in order to analyze kerogen conversion to liquid oil. Upon successfully verifying the code it is intended that numerical modeling studies be undertaken to address parametric uncertainty analysis, subsurface formation description, solid-liquid conversion rates and mechanical formation response due to kerogen to oil conversion using the multiphysics TPME framework. While dimensionality and parametric description of TPME quantities are fundamental to the understanding of in-situ pyrolysis for a target formation, this work is focused on the mathematical verification of the underlying coupled equations. An explicit coupling scheme, analogous to the description provided by Dean et al. [13], was implemented where electric potential is solved with a quasi-static Maxwell equation, then the Allen-Cahn phase field method is used to characterize solid-liquid transition, the enthalpy equation is used to solve for temperature then the governing mechanical equilibrium equation is solved to describe the mean normal stress as a function of temperature. To the best of the authors’ knowledge, no commercial, or otherwise researched multiphysics TPME simulator solution, has been developed to evaluate the conversion of electromagnetic energy to thermal energy in a subsurface rock, such that the conversion of solid kerogen to liquid oil can be tracked as a moving interface.

Much of the work that has been conducted in developing in-situ pyrolysis has been proprietary thus details describing field trials have not been readily accessible in the public domain. As a result, there has been no attempt to verify the derived results with those obtained during field trials. Instead verification of the coupled governing equations describing TPME processes specific to radio frequency heating are undertaken by way of the Method of Manufactured Solutions (MMS). The executed study is primarily mathematical and only substantively related to the actual physical parametric description oil shale in-situ pyrolysis by radio frequency heating. Be that as it may, the main contribution of this work is the analytical verification of the developed multiphysics finite element method simulator, so that modeling the in-situ pyrolysis of oil shale by radio frequency heating can be achieved. The ability to model in-situ pyrolysis using radio frequency heating by a TPME approach is anticipated to enable environmental, economic, and technical analysis of production from oil shale formations; in this way, the associated advantages and disadvantages can be more accurately assessed.

2. Methods of Oil Shale In-Situ Pyrolysis

Disparate methods have been undertaken by oil and gas operators to evaluate the technical feasibility of commercial oil shale production. The more prominent processes that have been considered for commercial scale in-situ pyrolysis have included: Shell In-situ Conversion Process (ICP) [4] [14] [15] [16]; ExxonMobil Electrofrac™ [16] [17] [18] [19]; Chevron CRUSH technology [20]; and Raytheon-CF Technologies partnership for radio frequency heating with Critical Fluids Technology [20] to list a few.

The ExxonMobil Electrofrac™ process, as outlined in Hoda et al. [18], Kelkar et al. [19] and Martemyanov et al. [16] highlights a typical in-situ conversion of oil shale by heating the oil shale formation to pyrolyze the kerogen. The result is the generation of liquid hydrocarbon that can be conventionally produced. This process can be extended by creating hydraulic fractures in the rock formation and filling those fractures with conductive material, described as a mixture of calcined coke and cement slurry as well as graphite. Electricity is then conducted across the length of the fracture converting the filled fracture to a resistive heating element for in-situ pyrolysis. The main advantage of this process is that it produces gas of high calorific value. Rock displacement induced by Electrofrac heating was measured by multi-point extensometers as the kerogen conversion would weaken the structural integrity of the oil shale formation. This displacement can be considered as a disadvantage since it can negatively impact surface facilities during in-situ pyrolysis operations. Additionally, this method requires horizontal drilling and horizontal fracturing before conducting electricity through the conductive material. Published work illustrates numerically modeled independent electrical and thermal models [18], coupled mass transport and heat flow through combined porous and fractured media including equations of state [21], as well as coupled thermal-hydrological-mechanical-chemical (THMC) [19] numerical models.

The Shell ICP process involves heating the oil shale formation with resistive heating elements surrounded by a hexagonal pattern of heating wells with the converted oil and gas produced through a production well in the center of the ring [16]. This process has been shown to add value by producing good quality gas and oil. In order to protect the environment during the heating process Shell developed a freeze wall technology to isolate groundwater from the production zone [22]. The freeze wall is located concentrically about the heater and production wells thus impeding the groundwater flux into the productive zone. This freeze wall is formed by injecting an ammonia-water solution of −45˚C to freeze the proximal groundwater at shallow depths then deepening the injection to create a containment volume 224 ft across in the planer direction [15] [20]. The implementation of a freeze wall in this method represents added value by protecting proximal groundwater. The main disadvantages associated with this method include the requirement of many wells for production as well as heating, in addition to the need for long heating times to stimulate conversion.

The Chevron CRUSH process involves using pressurized and heated carbon dioxide to raise the temperature of the oil shale; however, it depends on large quantities of water and leads to negative impact on the environment [20] [23]. The result of the outlined process is fracturing of the oil shale formation to increase the contact area of kerogen in the matrix with the injected carbon dioxide. The carbon dioxide is then recycled and reheated so that it can be used to further extend fractures into the formation [24].

The coupled TPME computational framework developed for this investigation reflects the Raytheon Radio frequency with Critical Fluids Technology process described by Allix et al. [15] and Pan et al. [20]; which represents a dielectric heating method. The description of the physical system is shown in Figure 1 for reference highlighting a homogeneous oil shale as kerogen-rich according to a

Figure 1. Generalized subsurface description of the target oil shale zone of interest including the idealized model domain.

numerical description without consideration of grade or constituent mineral components. Here the electromagnetic energy from the electrode is converted to thermal energy so that it heats the subsurface oil shale, enabling in-situ pyrolysis of kerogen to take place; then critical CO₂ is injected to displace the oil to the production well. The CO₂ is separated and recycled for further usage as a part of the process. While the injection of CO₂ is not considered in this investigation, the focus of this work is instead the numerical modeling of radio frequency dielectric heating in the subsurface given specific mathematical considerations. Of particular interest in the use of the multiphysics framework is the conversion of solid kerogen to liquid oil for production. The advantage of the radio frequency heating considered in the Raytheon process is that oil and gas can be produced in months compared to other methods which have a longer conversion and production times, order of years. One of the main disadvantages of radio frequency heating that has been identified is the reduction in efficiency that occurs as half of the generated electromagnetic energy is lost to the formation surrounding the installed electrodes in the conversion process [15]. Given the lower heating requirements; however, radio frequency heating has recently been considered a suitable enhanced recovery technique in tar sands and heavy oil assets [25].

Radio frequency heating represents a dielectric heating process where electromagnetic energy is converted to thermal energy in the oil shale. Molecules in the kerogen of the shale formation experience a change in electrical polarity as they are introduced to an alternating electromagnetic field from electrodes that are placed in the subsurface rock; this leads to increased thermal energy. Once conversion temperatures are reached, the solid kerogen in the oil shale is converted to liquid oil. As a result, kerogen conversion can be characterized as a moving interface, or phase field, defined by the conversion of a solid to a liquid. As the solid kerogen is converted to liquid, the converted zone is anticipated to become mechanically unstable as it may no longer support the stress of the overburden rock. Given this outcome, the TPME was developed to specifically model thermal, phase field, mechanical, and electromagnetic processes that constitute the radio frequency heating process.

3. Numerical Theory

A complex electrostatic potential solution is obtained by solving a quasi-static Maxwell equation for “lossy” dielectric material [26] [27]. An approximation was applied considering frequencies in the range of 100 kHz to 1 MHz had wavelengths that were large compared to the separation distance of the in-place electrodes. This allows the electromagnetic field to be approximated as a static, linear in electrical potential. A power conversion term (P) is constructed as a function of the applied alternating electromagnetic field. The description of the quasi-static Maxwell equation is given according to:

$-\nabla \cdot \left[\left(\sigma +j2\pi f{\epsilon }_0{\epsilon }^{\prime }\right)\nabla V\right]=0$ (1)

where ${\epsilon }^{″}$ is the loss factor of oil shale (imaginary part) describing the material’s ability to convert the electromagnetic field energy to heat, ${\epsilon }^{\prime }$ is the relative dielectric constant of oil shale (real part) describing the lossless energy interaction of the material, j is equal to (−1)^0.5 and V is the electric potential difference. The power conversion term P is expressed as:

$P=2\pi f{\epsilon }_0{\epsilon }^{″}{\left|E\right|}^2=\sigma {\left|-\nabla V\right|}^2$ (2)

which describes the conversion of electromagnetic to thermal energy.

The power conversion term P is used to couple the quasi-static Maxwell equation to the enthalpy and Allen-Cahn phase field equations in order to compute the temperature field. The radio frequency power conversion term causes the temperature of the kerogen to increase and eventually leads to kerogen upgrading to oil. The physics of the actual TPME numerical framework is defined such that interpolation of all intrinsic rock type properties -i.e. thermal conductivity, occurs as a function of the Allen-Cahn phase field and temperature solutions. This method is typically applied in capturing parametric transitions across an interface within phase field models.

The characterization of the fluid conversion interfaces during the modeled pyrolysis of kerogen to shale oil is performed using an Allen-Cahn phase field description. In mathematical physics, the Allen-Cahn equation represents a reaction-diffusion equation which can describe physical processes like liquefaction, such as what takes place as a result of dielectric heating of oil shale. The results illustrate a moving interface (spatio-temporal) expressed through a highly localized area of the domain. The Allen-Cahn phase field equation of Dyja et al. [28] was used for the basis of modeling the evolution of the phase field interface. In order to account for the phase transition during electromagnetic energy conversion to thermal energy, the power conversion term is included in the Allen-Cahn equation explicitly and scaled by the bulk modulus (K_m) of the reference material. The Allen-Cahn equation is described as:

$\frac{\partial \varphi }{\partial t}-D{C}_n^2{\nabla }^2\varphi +2DA\varphi \left(1-3\varphi +2{\varphi }^2\right)-Dk-\frac{P}{K_m}=0$ (3)

By convention, the phase field term, $\varphi$, and D, a scalar term, are dimensionless. As for the remaining terms; C_n is a diffusion coefficient, A and k are respectively frequency parameters for the solid-liquid interface, with A being a scaling factor and k being a bulk adjusting term.

The enthalpy equation described in Belhamadia et al. [29] was used as a starting point for computing system temperature in this study. This formulation included latent heat effects, negated the effects of the velocity field in the solid, and assumed negligible flow rates in the liquid. This definition also included the assumption that the fluid is incompressible and Newtonian. It is well understood that phase change does not instantaneously occur in a given rock type. This is due to variations in the distribution of spatial rock type characteristics that cause phase change to manifest over a small temperature range about a defined melting temperature $\left[T_m-ϵ,T_m+ϵ\right]$ –where $T_m$ is the melting temperature, and $ϵ$ is a relatively small numerical value defining quantitative variation. Thus, phase temperature is treated in a similar manner in its interpretation for this study. Even though the chemical kinematic model for the upgrading of oil shale does involve the decomposition of kerogen into shale oil (liquid) and gas as a result of the high pyrolysis temperatures, a limited two-phase description was employed which only considered solid (oil shale) conversion to a liquid (oil) phase.

The modified enthalpy equation for this multiphysics model was required to include the power conversion term, describing electromagnetic to thermal energy conversion, and the solid-liquid phase transition during pyrolysis. Belhamadia et al. [29] utilized the enthalpy-porosity model to perform numerical simulation of water solidification and pure gallium melting thus it was believed to similarly describe the thermophysical process under consideration in oil shale in-situ pyrolysis. The modified enthalpy equation which was used in the study is described as:

$\alpha \left(\varphi \right)\frac{\partial T}{\partial t}+{\rho }_l{L}_s\frac{\partial \varphi }{\partial t}-\kappa \left(\varphi \right){\nabla }^{2}T=P$ (4)

Here $\alpha \left(\varphi \right)$ is the volumetric heat capacity defined by density times the specific heat capacity as a function of the phase variable $\varphi$. It should be noted that phase related physical properties across the phase transition interface are interpolated between the respective pure solid and liquid phases. As for the remaining terms, the density of the liquid phase is ${\rho }_l$, the latent heat of fusion is given by L_s and the thermal conductivity, as a function of the phase variable, is expressed as $\kappa \left(\varphi \right)$.

It is anticipated that with suitably weak formation properties the upgraded zone becomes unable to support the stress of the overburden once the in-place solid kerogen upgrades to liquid oil. In this structurally unstable situation, as the liquid can no longer support the overburden stress, compaction of the target formation ensues. The complexity of such poro-elastic and plastic deformation exhibited by oil shale as a result of in-situ pyrolysis has been extensively discussed [19] [30] [31] [32]. The mean normal stress is computed as a primary variable while the equation is modified so that stress is a function of temperature. The corresponding mechanical equilibrium equation from Fakcharoenphol et al. [33] is described as:

$\frac{3\left(1-\upsilon \right)}{\left(1+\upsilon \right)}{\nabla }^2{\sigma }_m+\nabla \cdot F-\frac{2\left(1-2\upsilon \right)}{1+\upsilon }\left(3\beta K_m{\nabla }^{2}T\right)=0$ (5)

Expressed within this equation is the Poisson’s Ratio $\upsilon$, a force density term $F$ which is taken to be the force density due to the overburden and the linear thermal expansion coefficient is $\beta$. Within this equation the Poisson’s Ratio, linear thermal expansion coefficient and bulk modulus are functions of the phase variable.

The description of the explicitly coupled TPME formulation is shown in Figure 2. The terms that are used to couple each of the equations are highlighted and labeled as “coupling parameters”. According to the outlined process, the power conversion term is first computed based on the resulting electric potential difference. This solution is then passed to the Allen-Cahn equation that computes the phase field solution. The power conversion term and phase field parameter are passed to the enthalpy equation so that the system temperature can be computed. Lastly, once the temperature solution is obtained, the result is included in the mechanical equilibrium equation so that the mean normal stress can be computed. Next, MMS is used in order to analytically verify the results of the developed TPME code.

4. Results of the Method of Manufactured Solutions

4.1. Verification Method Background and Model Application

Few studies have been published about radio frequency heating as an enhanced

Figure 2. TPME numerical framework description highlighting electromagnetic, phase field, enthalpy and mechanical equilibrium components of the coupled equations.

oil recovery method or about numerical modeling of the same. Furthermore, results of field trials, recoverable volumes of hydrocarbon and economic viability of in-situ recovery methods remain closely guarded. The intellectual property that has been developed in conjunction with these methods are overwhelmingly reserved to maintain strategic advantage by the respective corporate entities which developed or acquired the specific upgrading technology. This lack of published production and numerical modeling results was the main motivation for at least verifying the mathematical and algorithmic implementation of the multiphysics framework.

The employed technique for code verification in this study was the MMS. While this technique was developed some time in the 1980’s or 1990’s there is no exact attribution as many have claimed to use it without specifically referencing it. Be that as it may, a description of the technique may be reviewed in Salari and Knupp [34].

According to Salari and Knupp [34], the following guidelines exist for devising manufactured solutions for computational code verification:

1) The manufactured solutions should be comprised of smooth analytical functions so that solutions may be easily computed. These would include trigonometric or polynomial functions. This guideline is designed to ensure that theoretical order-of-accuracy is attainable.

2) The solution should sample every term in the governing equation being evaluated.

3) An appropriate number of non-trivial derivatives for the solution should exist.

4) The solution derivatives should be constrained by a small constant and not contain singularities.

5) Successful execution of the code should not be impeded by the imposed solution.

6) The solution domain should be defined in a connected subset of the modeled space.

7) The differential operators in the partial differential equations should be formed in a manner that is logical.

In each case the MMS was used to verify the coupled equations that were incorporated into the multiphysics TPME solution. The verification was performed on the coupled code so that coupling terms, and non-solution variables were modified to be constants in this way the primary variables are independent of implemented coupling terms. In each case the MMS was applied to a 2D rectilinear FEM mesh with dimensions [0:1] × [0:1] and 20 elements in the X-direction and Y-direction, respectively. Essential boundary conditions were enforced on each of the boundaries in the quasi-static Maxwell and mechanical equilibrium equations. Conversely, in the Allen-Cahn and enthalpy equations, essential boundary conditions were applied to the left and right boundaries while natural boundary conditions were applied on the top and bottom boundaries. Appendix A contains the list of the parameters used in the underlying code verification performed by MMS. It is important to note that the listed parameters have no physical significance with respect to oil shale in-situ pyrolysis but are considered mathematically relevant to the verification of the code using MMS.

4.2. Quasi-Static Maxwell Equation

Given the aforementioned guidelines a solution to V for the quasi-static Maxwell equation in Equation (1) was set to the following analytical spatial function for $m\left(x,y\right)$ :

$m\left(x,y\right)=\sin \pi x\sin \pi y$ (6)

An illustration of Equation (6) is shown in Figure 3 as a reference. In order to use MMS, Equation (6) is substituted into the V term of Equation (1), leading to Equation (7). Next, the right-hand side of Equation (1) is set to $\Gamma \left(x,y\right)$ in Equation (7):

$\Gamma \left(x,y\right)=\left[\sigma +j2\pi f{\epsilon }_0{\epsilon }^{\prime }\right]\left(2{\pi }^2\sin \pi x\sin \pi y\right)$ (7)

Equation (1) is modified to include Equation (7) in its right-hand side. It is then solved for $V\left(x,y\right)$ using the scalable linear equation Krylov solver (KSP) in the PETSc framework [35] [36] [37] that was called from the TalyFEM libraries. The resulting solutions, real and imaginary parts of $V\left(x,y\right)$, were found to match the assumed solution from Equation (7). Since the solution for $V\left(x,y\right)$ followed that of Equation (7) the implementation of the code corresponding to the solution of the quasi-static Maxwell equation was verified within the coupled framework. The similitude of the solution is verified by comparing the computed

Figure 3. Function $m\left(x,y\right)=\sin \pi x\sin \pi y$ used in MMS as an analytical solution to static equations of V and ${\sigma }_m$.

solution to $m\left(x,y\right)$ shown in Figure 3 to that of the resulting finite element solution to $V\left(x,y\right)$ in Figure 4 for the real (Vreal) and imaginary (Vimg) parts of electric potential difference in Equation (1). A direct comparison of the real and imaginary parts of the electric potential difference to $m\left(x,y\right)$ is performed on a node by node basis in order to evaluate solution continuity. Cross-plots of these comparisons are shown in Figure 5(a) & Figure 5(b), for the real and imaginary parts, respectively; the solution to $m\left(x,y\right)$ is shown as “v_a” in the plot. Figure 5(a) includes a linear trend line that is described by an R² coefficient equal to unity; this suggests that the solutions at each node are consistent with the analytical solution. This is further evidenced by the slope and y-intercept of the linear trend line equation that has values of 0.9938 and 4e−07, respectively. As a corollary, the anticipated slope of an exact match in solution is unity and the y-intercept is zero. The identical solution is obtained for the imaginary part of the electric potential difference in Figure 5(b). The deviations of $V\left(x,y\right)$ from the idealized solution $m\left(x,y\right)$ are illustrated as a solution map of absolute node error, computed as $\left|V\left(x,y\right)-m\left(x,y\right)\right|$ in Figure 6(a) & Figure 6(b). The results in Figure 6(a) & Figure 6(b) show the best match in solution occurs near the boundaries while largest deviation occurs

Figure 4. Verification by MMS for the quasi-static Maxwell equation; (a) real part and (b) imaginary parts.

Figure 5. Cross-plot including Linear trend line and R² coefficient for (a) real and (b) imaginary components of the electric potential difference.

Figure 6. Absolute node error of the difference between the $m\left(x,y\right)$ solution with the computed (a) real and (b) imaginary electric potential difference components.

near the center of the domain.

4.3. Allen-Cahn Phase Field Equation

The code for the Allen-Cahn equation, from Equation (3), was verified by assuming the phase variable $\varphi$ as equal to the following spatio-temporal analytical function for $g\left(x,t\right)$ :

$g\left(x,t\right)=\sin \pi x\sin \pi t$ (8)

A contour plot of Equation (8) is shown in Figure 7 for reference. MMS was performed by substituting Equation (8) into Equation (3) as a solution to the phase field variable yields:

$\begin{array}{c}\chi \left(x,t\right)=\pi \sin \pi x\cos \pi t+D{C}_n^2{\pi }^2\sin \pi x\sin \pi t\\ +2DA\sin \pi x\sin \pi t\left[1-3\sin \pi x\sin \pi t+2{\left(\sin \pi x\sin \pi t\right)}^2\right]\\ -D\kappa -1\times {10}^{-1}\end{array}$ (9)

Equation (3) was then revised so that Equation (9) became the right-hand side of the equation and then solved for the phase field variable $\varphi \left(x,t\right)$. The solution for the phase field variable was then obtained in the numerical model by using the Newton-Raphson method of the SNES solver in the PETSc framework [35] [36] [37]; which was called through the TalyFEM libraries. By comparing the plot of Equation (8), represented in Figure 7, with the result of the finite element solution of $\varphi \left(x,t\right)$ at selected times (t), as shown in Figure 8, it is clear that the Allen-Cahn equation, Equation (3), has been properly implemented in the code. This is confirmed by comparing the solutions of $g\left(x,t\right)$ and $\varphi \left(x,t\right)$ between Figure 7 and Figure 8. It should be noted that the power conversion

term of $\frac{P}{K_m}$ in the coupling between the quasi-static Maxwell equation and the

Allen-Cahn equation was fixed to a constant value of 1 × 10⁻¹ in order to independently verify the implementation of the Allen-Cahn equation in the coupled

Figure 7. Function $g\left(x,t\right)=\sin \pi x\sin \pi t$ used in MMS for spatio-temporal Allen-Cahn phase field ( $\varphi$ ) and temperature (T): (a) t = 6; (b) t = 25; (c) t = 50 and (d)t = 100.

framework in the absence of complimentary coupling terms. The cross-plot of node values for $\varphi \left(x,t\right)$ -as the variable “u” in the x-axis, and the function $g\left(x,t\right)$ -as the variable “u_a”, is shown in Figures 9(a)-(d) for select times. An alternating positive-negative response is observed in these results and is attributed to the time-dependent sinusoidal analytical function $g\left(x,t\right)$ used in the MMS. In addition to Figure 9, this behavior is also observed in Figure 8 where the finite element solution to MMS is shown. The results of Figure 9 demonstrate that while the solutions were visually congruent there are discernible quantitative differences. In each case the y-intercept of the added trend line is practically zero; however, the slope of the trend line varies from 0.0054 to 0.159. The slope of the trend line tends to be consistent when t = 50 and t = 68, when the phase of the solution is the same (positive value/phase). When t = 63 and t = 75 the phase of the solution differs (negative value/phase); however, the slope of the trend line increases as a function of time. Be that as it may, there is greater solution accuracy when t = 50 and t = 68 as R² in both cases equals to 0.9997. This is compared to t = 63 where R² equals 0.9983 andt = 75 where R² equals 0.9661. Even though the t = 75 solution has a trend line slope that is closer to unity it is important to note that R² value is more indicative of the relationship between $\varphi \left(x,t\right)$ and $g\left(x,t\right)$ than is the slope of the trend line. Similar results are illustrated in the absolute error maps of $\left|\varphi \left(x,t\right)-g\left(x,t\right)\right|$ in Figures 10(a)-(d) where changes in the solution at the node are observed over time but appear to

Figure 8. Verification of the Allen-Cahn equation by the MMS: a) t = 6; b) t = 25; c) t = 50 and d) t = 100.

Figure 9. Cross-plot including Linear trend line and R² coefficient for solutions of the phase field equation at (a) t = 50, (b) t = 63, (c) t = 68, and (d)t = 75 for dt= 1.0.

Figure 10. Absolute node error of the difference between the $g\left(x,t\right)$ solution with the solution to the Allen-Cahn phase field at (a) t = 50, (b) t = 63, (c) t = 68, and (d)t = 75 with dt = 1.0.

be consistent att = 50 and t = 68. Solutions with the greatest accuracy occur near the left and right-side boundaries where essential boundary conditions were defined in contrast to the top and bottom side boundaries of the domain where natural boundary conditions were defined.

4.4. Enthalpy Equation

The code for the enthalpy equation, Equation (4), remained coupled, by definition, to Equation (3) through the phase field variable. As a result, this code was not verified in the absence of or the explicit fixing of a constant phase variable value using MMS. The assumption was made that the temperature variable T was set to Equation (8) to provide a solution to the enthalpy equation.

Substitution was performed for Equation (8) into Equation (4) and upon rearranging the following result was derived:

$\gamma \left(x,t\right)=\alpha \left(\varphi \right)\pi \sin \pi x\cos \pi t+{\rho }_l{L}_s\frac{\partial \varphi }{\partial t}+\kappa {\pi }^2\sin \pi x\sin \pi t-1.0$ (10)

The updated Equation (4) was then modified so that P was subtracted from each side of the equation, then Equation (10) was set to be the right-hand side of the equation then solved for temperature. Here the power conversion term, used to couple the quasi-static Maxwell equation to the enthalpy equation, was set to a constant value of unity solely for the purpose of code verification. This minimized the impact of the solution to the electromagnetic equation by isolating the temperature term. The phase field term was maintained as it was a component of the enthalpy equation that was used as a starting point the numerical model. The temperature solution, obtained using the MMS, is illustrated in Figure 11. Once juxtaposed to Figure 7 it is clear that the solution matches the bulk description of the function definition for temperature following Equation (8). The largest deviation in response is roughly a factor of 200; however, the accuracy of the match is critiqued based on the overall trend as well as consideration that the magnitude of respective solutions is very small. Specifically, the magnitude of the temperature and the $g\left(x,t\right)$ solutions has very small values. The cross-plot of node values for $T\left(x,t\right)$ as the variable “h” in the x-axis, and the function $g\left(x,t\right)$ as the variable “h_a”, is shown in Figure 12 at select times. The results show a consistent linear trend between the computed function $T\left(x,t\right)$ and the assumed $g\left(x,t\right)$ function for the listed times in Figure 12 as each case illustrates an R² coefficient of unity. While the y-intercept in each plot of Figure 12 is approximately zero there is a change in magnitude of the slope between the computed and assumed solutions that illustrates that the computed solution is larger in absolute magnitude than the analytical MMS solution. An observation of systematic time dependent variation in the solutions is observed in Figure 12(a), Figure 12(c) for t = 50 and 68, then Figure 12(b), Figure 12(d) for t = 63 and 75 as two distinct groups. The differences in these two groups are attributed to the time-dependent sinusoidal behavior of $g\left(x,t\right)$ which is illustrated in the positive-negative value alternating solution response in Figure 11 as a function of time and quantified in the identical trend line slopes computed for the cross-plots of similar phase. Recall, a similar observation was made in Figure 8 and Figure 9 for the Allen-Cahn phase field solution. The trend line slopes are identical for Figure 12(a), Figure 12(c) collectively then Figure 12(b), Figure 12(d) as a separate group. A depiction of the absolute node error $\left|T\left(x,t\right)-g\left(x,t\right)\right|$ is shown in Figures 13(a)-(d). Again, the results show the most accurate solution along the left and right-side boundaries where essential boundary conditions were set. Compared to the absolute node error of the phase field solution in Figures 10(a)-(d), the results of $\left|T\left(x,t\right)-g\left(x,t\right)\right|$ show an

Figure 11. Temperature solution from the enthalpy equation solved by the MMS: (a) t = 6; (b) t = 25; (c) t = 50 and (d)t = 100.

Figure 12. Cross-plot including Linear trend line and R² coefficient for solutions of the temperature at (a) t = 50, (b)t = 63, (c) t = 68, and (d) t = 75 for dt= 1.0.

Figure 13. Absolute node error of the difference between the $g\left(x,t\right)$ solution with the solution for temperature at (a) t = 50, (b)t = 63, (c)t = 68, and (d) t = 75 with dt= 1.0.

increase in absolute node error with time fromt = 50 until t = 68 but then the absolute node error decreases oncet = 75. While additional times were not analyzed it is anticipated that due to the sinusoidal behavior of $g\left(x,t\right)$ the response of the absolute node error would continue to alternate with time.

4.5. Mechanical Equilibrium Equation

The governing mechanical equation, Equation (5), was verified while maintaining the coupling with the temperature solution obtained by solving the enthalpy equation, Equation (4). The assumption was made that the mean normal stress ${\sigma }_m$ was equal to Equation (6) which is shown in Figure 3. Equation (6) was then substituted into ${\sigma }_m$ of Equation (5) and the terms were rearranged leading to the following equation:

$\zeta \left(x,y\right)=\frac{-6\left(1-\upsilon \right){\pi }^2\sin \pi x\sin \pi y}{1+\upsilon }-\frac{2\left(1-2\upsilon \right)}{1+\upsilon }\left(3\beta K_m{\nabla }^{2}T\right)$ ( 1 )

In the next step, Equation (5) was revised so that Equation (11) was set as the right-hand side of Equation (5) then solved for the mean normal stress. In this case, the enthalpy equation was coupled to the mechanical equilibrium equation through an enthalpy-stress coupling parameter that was set to a value of 1.6e1. The result of the executed mechanical simulation code by MMS led to a solution for the mean normal stress which is shown in Figure 14. Comparing the results for the mean normal stress in Figure 14 with the analytical mean normal stress described in Figure 3, demonstrates that the model and analytical descriptions are coincident. This observation is further maintained by analyzing cross-plot of the analytical mean normal stress solution “s_a” versus the mean normal stress “s” computed by the MMS equation in Figure 15. The results in Figure 15 show an R² coefficient as well as a trend line that is described by a slope of unity and a y-intercept in the trend line equation that is approximately zero. The cross-plot highlights the accurate match of the MMS computed solution to that of the initial analytical function. The spatial description of the deviation between the

Figure 14. Mean normal stress solution from the mechanical equilibrium equation solved by MMS.

MMS solution and the analytical solution at each node, $\left|{\sigma }_m\left(x,y\right)-m\left(x,y\right)\right|$ is shown in Figure 16. In Figure 16, there is greater deviation in the solutions towards the center of the domain compared to significantly less at the boundaries. Again, due to the essential boundary conditions that are applied to the numerical model there is less deviation between solutions at the boundaries. Even though the solution deviates toward the center of the domain, the maximum deviation is rather small in magnitude and does not exceed 6.3e−3 in absolute value. These results lead to the determination that the mechanical equilibrium equation is accurately coded in the coupled TPME computational framework.

Figure 15. Cross-plot including Linear trend line and R² coefficient for the mean normal stress.

Figure 16. Absolute node error of the difference between the $m\left(x,y\right)$ solution with the computed mean normal stress.

5. Conclusion

A multiphysics finite element computational framework was developed which explicitly coupled thermal, phase field, mechanical and electromagnetic equations for the purpose of numerically simulating the in-situ pyrolysis of oil shale by radio frequency heating. Specific processes for in-situ pyrolysis were discussed but the multiphysics solution proposed in this work was specifically designed to address radio frequency heating as an enhanced oil recovery method. There have been limited publications that have addressed either the development of multiphysics numerical models for simulating the in-situ pyrolysis process or production from oil shale using radio frequency heating. Furthermore, very few publications have demonstrated the development or use of multiphysics finite element solutions to model in-situ radio frequency heating. As a result, verification of the developed multiphysics finite element model necessitated the use of MMS. The results from the use of MMS in this work show that the electromagnetic, phase field, enthalpy and mechanical equilibrium solutions are respectively and collectively implemented in an accurate manner in the computational TPME framework. The observed matches between the computed and analytical solutions of the coupled equations using MMS demonstrated combinations of negligible or numerically insignificant differences between the numerical solutions. As a corollary, the TPME numerical framework has been verified for use in numerically modeling oil shale in-situ pyrolysis by radio frequency heating.

Acknowledgements

The author wishes to thank Dr. Baskar Ganapathysubramanian for fruitful discussions and supervision in the development of the multiphysics TPME finite element framework.

Nomenclature

A = Allen-Cahn frequency parameter, t⁻¹, s⁻¹

$C_n$ = Allen-Cahn diffusion coefficient, L²/t, m²/s

D = Allen-Cahn scaling factor, dimensionless

E = quasistatic electric field, V/L, V/m

f = frequency, t⁻¹, s⁻¹

$F$ = force density, F/L³, N/m³

g = general spatio-temporal function

j = unit imaginary number, $\sqrt{-1}$

k = Allen-Cahn frequency parameter, t⁻¹, s⁻¹

$K_m$ = bulk modulus, m/Lt², Pa

$L_s$ = latent heat of fusion (solid), E/m, J/kg

m = general spatial function

P = power conversion term, t³I²V²/L³m³, s³∙A²∙V²/m³∙kg³

t = time, t, s

T = temperature, T, K

V = voltage, V, V

x = distance, L, m

y = distance, L, m

$\alpha$ = volumetric heat capacity, E/L³T, J/m³∙K

$\beta$ = linear thermal expansion coefficient, T⁻¹, K⁻¹

$\gamma$ = general spatio-temporal function

$\Gamma$ = general spatial function

$ϵ$ = relatively small numerical value

${\epsilon }_0$ = electromagnetic permittivity of free space, t⁴I²/L³m, F/m

${\epsilon }^{\prime }$ = electromagnetic relative dielectric constant, dimensionless

${\epsilon }^{″}$ = electromagnetic loss factor, dimensionless

$\zeta$ = general spatial function

$\kappa$ = thermal conductivity, P/LT, W/m∙K

$\rho$ = density, m/L³, kg/m³

${\rho }_l$ = density (liquid), m/L³, kg/m³

$\sigma$ = electrical conductivity, t³I²/L³m, s³∙A²/m³∙kg

${\sigma }_m$ = mean normal stress, m/Lt², kg/m∙s²

$\upsilon$ = Poisson’s ratio, dimensionless

$\varphi$ = Allen-Cahn phase field, dimensionless

$\chi$ = general spatio-temporal function

Appendix: A Summary of Parameters for the Coupled TPME Code Verification

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Altun, N.E., Hicyilmaz, C., Hwang, J.Y., Suat Bagci, A. and Kok, M.V. (2006) Oil Shales in the World and Turkey; Reserves, Current Situation and Future Prospects: A Review. Oil Shale, 23, 211-227.
[2]	Vandenbroucke, M. and Largeau, C. (2007) Kerogen Origin, Evolution and Structure. Organic Geochemistry, 38, 719-833. https://doi.org/10.1016/j.orggeochem.2007.01.001
[3]	Peters, K.E., Walters, C.C. and Moldowan, J.M. (2005) The Biomarker Guide, Volume 2: Biomarkers and Isotopes in Petroleum Exploration and Earth History. 2nd Edition, Cambridge University Press, Cambridge.
[4]	Fan, Y., Durlofsky, L. and Tchelepi, H. (2010) Numerical Simulation of the In-Situ Upgrading of Oil Shale. SPE Journal, 15, 368-381. https://doi.org/10.2118/118958-PA
[5]	Potter, J., Craddock, P.R., Kleinberg, R.L. and Pomerantz, A.E. (2017) Downhole Estimate of the Enthalpy Required to Heat Oil Shale and Heavy Oil Formations. Energy and Fuels, 31, 362-373. https://doi.org/10.1021/acs.energyfuels.6b02495
[6]	Brandt, A.R. (2008) Converting Oil Shale to Liquid Fuels: Energy Inputs and Greenhouse Gas Emissions of the Shell in Situ Conversion Process. Environmental Science & Technology, 42, 7489-7495. https://doi.org/10.1021/es800531f
[7]	Sahni, A., Kumar, M. and Knapp, R. (2000) Electromagnetic Heating Methods for Heavy Oil Reservoirs. Proceedings of the SPE/AAPG Western Regional Meeting, Long Beach, 19-23 June 2000, SPE-62550. https://doi.org/10.2118/62550-MS
[8]	Davletbaev, A., Kovaleva, L. and Babadagli, T. (2011) Mathematical Modeling and Field Application of Heavy Oil Recovery by Radio frequency Electromagnetic Stimulation. Journal of Petroleum Science and Engineering, 78, 646-653. https://doi.org/10.1016/j.petrol.2011.07.006
[9]	Mukhametshina, A. and Martynova, E. (2013) Electromagnetic Heating of Heavy Oil and Bitumen: A Review of Experimental Studies and Field Applications. Journal of Petroleum Engineering, 2013, Article ID: 476519. https://doi.org/10.1155/2013/476519
[10]	Kinzer, D. (2008) Past, Present and Pending Intellectual Property for Electromagnetic Heating of Oil Shale. Proceedings of the 28th Oil Shale Symposium, Golden, 13-15 October 2008, 1-18.
[11]	Pan, Y., Chen, C. and Yang, S. (2012) Development of Radio Frequency Heating Technology for Shale Oil Extraction. Open Journal of Applied Sciences, 2, 66-69. https://doi.org/10.4236/ojapps.2012.22008
[12]	Han, H., et al. (2016) Numerical Simulation of in Situ Conversion of Continental Oil Shale in Northeast China. Oil Shale, 33, 45-57. https://doi.org/10.3176/oil.2016.1.04
[13]	Dean, R.H., Gai, X., Stone, C.M. and Minkoff, S.E. (2006) A Comparison of Techniques for Coupling Porous Flow and Geomechanics. SPE Journal, 11, 132-140. https://doi.org/10.2118/79709-PA
[14]	Crawford, P.M., Biglarbigi, K., Dammer, A.R. and Knaus, E. (2008) Advances in World Oil Shale Technologies. Proceedings of the SPE Annual Technical Conference and Exhibition, Denver, 21-24 September 2008, SPE-116570-MS. https://doi.org/10.2118/116570-MS
[15]	Allix, P., Burnham, A., Fowler, T., Herron, M., Kleinberg, R. and Symington, B. (2011) Coaxing Oil from Shale. Oilfield Review, 22, 4-15.
[16]	Martemyanov, S., Bukharki, A., Koryashov, I. and Ivanov, A. (2016) Analysis of Applicability of Oil Shale for in Situ Conversion. Proceedings of the AIP XIII International Conference of Students and Young Scientists, Tomsk, 26-29 April 2016, 1-6.
[17]	Symington, W.A., Olgaard, D.L., Otten, G.A., Phillips, T.C., Thomas, M.M. and Yeakel, J.D. (2008) ExxonMobil’s Electrofrac Process for in Situ Oil Shale Conversion. Proceedings of the AAPG Annual Convention, San Antonio, 20-23 April 2008, Search and Discovery #40316.
[18]	Hoda, N., Fang, C., Lin, W.M., Symington, W.A. and Stone, M.T. (2010) Numerical Modeling of ExxonMobil’s Electrofrac TM Field Experiment at Colony Mine. Proceedings of the 30th Oil Shale Symposium, Golden, 18-20 October 2010, 1-13.
[19]	Kelkar, S., Pawar, R. and Hoda, N. (2011) Numerical Simulation of Coupled Thermal-Hydrological-Mechanical-Chemical Processes During in Situ Conversion and Production of Oil Shale. Proceedings of the 31st Oil Shale Symposium, Golden, 17-19 October 2011, 1-8.
[20]	Pan, Y., Mu, J., Ning, J. and Yang, S. (2012) Research on In-Situ Oil Shale Mining Technology. Journal of Pharmaceutical Science Invention, 1, 1-7.
[21]	Lee, K., Moridis, G. and Ehlig-Economides, C. (2016) A Comprehensive Simulation Model of Kerogen Pyrolysis for the In-Situ Upgrading of Oil Shales. SPE Journal, 21, 1612-1630. https://doi.org/10.2118/173299-PA
[22]	Vinegar, H. (2006) Shell’s In-Situ Conversion Process. The 26th Oil Shale Symposium, Golden, 16-18 October 2006, 1-23. http://www.ceri-mines.org/documents/R05a-HaroldVinegar.pdf
[23]	Sun, Y.-H., Bai, F.T., Lu, X.-S., Li, Q., Liu, Y.-M., Guo, M.-Y., Guo, W. and Liu, B.-C. (2015) A Novel Energy-Efficient Pyrolysis Process: Self-Pyrolysis of Oil Shale Triggered by Topochemical Heat in Horizontal Fixed Bed. Scientific Reports, 5, Article No. 8290. https://doi.org/10.1038/srep08290
[24]	Biglarbigi, K., Dammer, A., Cusimano, J. and Mohan, H. (2007) Potential for Oil Shale Development in the United States. Proceedings of the SPE Annual Technical Conference and Exhibition, Anaheim, 11-14 November 2007, SPE-110590-MS, 1-15. https://doi.org/10.2118/110590-MS
[25]	Hazra, K. (2014) Comparison of Heating Methods for In-Situ Oil Shale Extraction. Master’s Thesis, Texas A&M University, College Station. https://doi.org/10.3997/2214-4609.20142631
[26]	Choi, C. and Konrad, A. (1991) Finite Element Modeling of the RF Heating Process. IEEE Transactions on Magnetics, 27, 4227-4230. https://doi.org/10.1109/20.105034
[27]	Huang, Z., Zhu, H. and Wang, S. (2015) Finite Element Modeling and Analysis of Frequency Heating Rate in Mung Beans. Transactions of the American Society of Agricultural and Biological Engineers, 58, 149-160.
[28]	Dyja, R., Ganapathysubramanian, B. and Van Der Zee, K.G. (2018) Parallel-in-Space-Time Adaptive Finite Element Framework for Nonlinear Parabolic Equations. SIAM Journal on Scientific Computing, 40, C283-C304. https://doi.org/10.1137/16M108985X
[29]	Belhamadia, Y., Kane, A.S. and Fortin, A. (2012) An Enhanced Mathematical Model for Phase Change Problems with Natural Convection. International Journal of Numerical Analysis and Modeling Series B, 3, 192-206.
[30]	Burnham, A.K. (2003) Slow Radio Frequency Processing of Large Oil Shale Volumes to Produce Petroleum-Like Shale Oil. US Department of Energy Report No. UCRL-ID-155045. https://doi.org/10.2172/15004663
[31]	Burnham, A.K. and McConaghy, J.R. (2006) Comparison of the Acceptability of Various Oil Shale Processes. Proceedings of the 26th Oil Shale Symposium, Golden, 16-18 October 2006, 1-11.
[32]	Tran, T.Q. and McLennan, J.D. (2016) Geomechanical and Fluid Transport Properties. In: Spinti, J., Ed., Utah Oil Shale: Science, Technology, and Policy Perspective, Taylor & Francis Group, Boca Raton, 209-240. https://doi.org/10.1201/9781315370873-10
[33]	Fakcharoenphol, P., Xiong, Y., Hu, L., Winterfeld, P.H., Xu, T. and Wu, Y.-S. (2013) User’s Guide of TOUGH2-EGS: A Coupled Geomechanical and Reactive Geochemical Simulator for Fluid and Heat Flow in Enhanced Geothermal Systems Version 1.0. Golden, U.S. Department of Energy Contract No. DE-EE0002762. https://doi.org/10.2172/1136243
[34]	Salari, K. and Knupp, P. (2000) Code Verification by the Method of Manufactured Solution. U.S. Department of Energy Contract No. DE-AC04-94AL85000. https://doi.org/10.2172/759450
[35]	Balay, S., Gropp, W.D., McInnes, L.C. and Smith, B.F. (1997) Efficient Management of Parallelism in Object Oriented Numerical Software Libraries. In: Arge, E., Bruaset, M. and Langtengen, H.P., Eds., Modern Software Tools in Scientific Computing, Birkhauser Press, Boston, 163-202.
[36]	Balay, S., Abhyankar, S., Adams, M.F., et al. (2018) PETSc Webpage. http://www.mcs.anl.gov/petsc
[37]	Balay, S., Abhyankar, S., Adams, M.F., et al. (2019) PETSc User’s Manual. Argonne National Laboratory ANL-95/11—Revision 3.11.