The hyperbolic heat transfer equation and the ablation problem: Theory and experiment

TL;DR

This paper investigates the hyperbolic heat transfer equation in an ablation context, providing an analytic solution, experimental measurements of thermal relaxation time, and emphasizing the importance of hyperbolic modeling over parabolic equations in electrocardiology ablation.

Contribution

It offers a new analytic solution for the hyperbolic heat equation in ablation, with experimental validation and implications for medical applications.

Findings

Thermal relaxation time τ is about 7 minutes for 0.5% NaCl in water.

Hyperbolic heat equation is necessary over parabolic models in ablation.

Analytic solution shows approach to steady state in ablation process.

Abstract

We study the ablation problem for the hyperbolic heat equation in an axisymmetrical geometry which can be conveniently realized in the lab. We determine an analytic solution which shows the approach to steady state. The thermal relaxation time $\tau$ is best obtained from the small time behavior. The measurements give a surprisingly large $\tau$ of about 7 minutes for 0.5 % NaCl in water. This shows that the hyperbolic equation must certainly be used instead of the parabolic heat equation in the ablation problem of electrocardiology.

Equations

$$\frac{{\partial }^{2}T}{\partial t^2}+\frac{1}{\tau }\frac{\partial T}{\partial t}-\frac{\kappa }{\tau \varrho c}△T=\frac{1}{\varrho c}\frac{\partial Q}{\partial t}+\frac{Q}{\tau \varrho c}.$$

$$V(r)=C_0+C_1\log r,$$

$$V(r_0)=V_0,V(r_1)=0$$

$$V(r)=V_0\frac{\log r-\log r_1}{\log r_0-\log r_1}.$$

$$E_r=-\frac{\partial V}{\partial r}=-\frac{V_0}{\log r_0-\log r_1}\frac{1}{r}$$

$$Q=\sigma E_r^2=\frac{\beta }{r^2}\mathrm{with}\beta =\sigma \frac{V_0^2}{(\log r_0-\log r_1{)}^2}=\frac{P}{2\pi L{\log }^2(r_0/r_1)}$$

$$T(t,r)=T_1(r)+T_2(t,r),$$

$${d\over dr}\Bigl{(}r{dT_{1}\over dr}\Bigl{)}=-{\beta\over\kappa r}.$$

$$T_1(r)=-\frac{b}{2}(\log r{)}^2+C_2\log r+C_3$$

$${dT_{1}\over dr}\Bigl{|}_{r_{0}}=0,$$

$$T_1(r_1)=T_0.$$

$$T_1(r)=T_0+b\log \frac{r_0}{r_1}\log \frac{r}{r_1}-\frac{b}{2r^2}{\log }^2\frac{r}{r_1}.$$

$$T_1(r_0)=T_0+\frac{b}{2{\log }^2\frac{r_1}{r_0}}.$$

$$T_2(t,r)=T_3(t)T_4(r).$$

$$\frac{1}{a}\frac{{\dot{T}}_3}{T_3}+a\tau \frac{{\ddot{T}}_3}{T_3}=\frac{T_4^{\prime \prime }}{T_4}+\frac{1}{r}\frac{T_4^{\prime }}{T_4}=const.=-k^2$$

$$a=\frac{\kappa }{\varrho c}.$$

$$a\tau {\ddot{T}}_3+\frac{1}{a}{\dot{T}}_3=-k^2{T}_3$$

$$r^2{T}_4^{\prime \prime }+r{T}_4^{\prime }=-k^2{r}^2{T}_4=0.$$

$$y(r)=\sqrt{r}{T}_4(r)$$

$$-y^{\prime \prime }-\frac{y}{4r^2}=k^{2}y.$$

$$y(r)=\sqrt{r}(D_1{J}_0(kr)+D_2{Y}_0(kr))$$

$$T_{4}(t,r_{1})=0,\quad{\partial T_{4}\over\partial r}\Bigl{|}_{r_{0}}=0$$

$$y(r_1)=0,r_0{y}^{\prime }(r_0)-\frac{1}{2}y(r_0)=0.$$

$$D_2(k)=-\frac{J_0(k{r}_1)}{Y_0(k{r}_1)}{D}_1(k)$$

$$D_1(k)J_0^{\prime }(k{r}_0)+D_2(k)Y_0^{\prime }(k{r}_0)=0$$

$$\frac{Y_0(k_n{r}_1)}{J_0(k_n{r}_1)}=\frac{Y_1(k_n{r}_0)}{J_1(k_n{r}_0)}.$$

$$y_n(r)=\sqrt{r}[Y_0(k_n{r}_1)J_0(k_{n}r)-J_0(k_n{r}_1)Y_0(k_{n}r)]\equiv \sqrt{r}{Z}_0(k_{n}r).$$

$$\underset{r_0}{\overset{r_1}{\int }}{y}_n^2(r)dr=\frac{1}{2}[r_1^2{Z}_1^2(k_n{r}_1)-r_0^2{Z}_0\text{texttwosuperior}(k_n{r}_0)].$$

$$k_n\approx \frac{{\alpha }_n}{r_1},$$

$$a\tau {\omega }^2+\frac{\omega }{a}+k^2=0$$

Full text

The hyperbolic heat transfer equation and the ablation problem: Theory and experiment

Günter Scharf111e-mail: scharf@physik.uzh.ch

Physics Institute, University of Zürich

Lam Dang222e-mail: lam.dang@gmail.com

HerzGefässZentrum, Klinik im Park

8022 Zürich

We study the ablation problem for the hyperbolic heat equation in an axisymmetric geometry which can be conveniently realized in the lab. We determine an analytic solution which shows the approach to steady state. The thermal relaxation time $\tau$ is best obtained from the small time behavior. The measurements give a surprisingly large $\tau$ of about 7 minutes for 0.5 % NaCl in water. This shows that the hyperbolic equation must certainly be used instead of the parabolic heat equation in the ablation problem of electrocardiology.

1 Introduction

The hyperbolic heat transfer equation

[TABLE]

has been applied to the ablation problem in previous studies ([1] and references given there). In (1.1) $T$ is temperature, $t$ is time, $\triangle$ is the Laplace operator and $Q$ describes the heat generation. The constant $\kappa$ is the thermal conductivity, $\varrho$ the mass density and $c$ the specific heat of the medium. However the value of the thermal relaxation time $\tau$ is badly known and consequently, it is unclear whether it is necessary to use (1) instead of the simpler parabolic heat equation. In [2] a value $\tau=16$ seconds has been measured in “processed meat” by methods using thermal conduction only, without heat generation. This value then has been used also in simulations of the ablation problem [1]. But if Joule’s heat is generated by an electric current the physics is completely different from pure heat conduction in [2]. The electric energy is first transferred to kinetic energy of moving charged ions (Na+ and Cl- in blood or water). Then it is dissipated to the water molecules where it is measured by the thermometers. This total relaxation process is much slower. In fact, we shall see that $\tau$ is now measured in minutes instead of seconds as in [2]. This means that we must certainly use equation (1.1) for the ablation problem in electrocardiology.

2 The axisymmetric ablation problem

Let us consider a cylindrical electrode of radius $r_{0}$ and length $L$ with potential $V_{0}$ in a medium with electrical conductivity $\sigma$. A second dispersive cylindrical electrode with radius $r_{1}$ ($r_{1}\gg r_{0}$) is grounded with potential $V=0$. The potential $V(r)$ in the medium is given by the simple solution of the two-dimensional Laplace’s equation

[TABLE]

if boundary effects at the end of the electrodes are neglected. Using the boundary conditions

[TABLE]

we have

[TABLE]

The corresponding electric field strength is

[TABLE]

and the heat generation is equal to

[TABLE]

where $P$ is the electric power in watts. With this heating we want to calculate the transient temperature $T(t,r)$ as solution of equation (1.1) in two dimensions.

We write the solution as

[TABLE]

where $T_{1}(r)$ is the steady state solution satisfying

[TABLE]

This is easily integrated

[TABLE]

where $b=\beta/\kappa$. The integration constants are fixed by the assumption of no heat flux at the electrode

[TABLE]

which is reasonable for a thin electrode, and by the second boundary condition

[TABLE]

This gives the following steady state solution

[TABLE]

The maximal temperature is found at the electrode, of course

[TABLE]

The remaining homogeneous equation for $T_{2}$ is solved by separation of variables

[TABLE]

Then we have

[TABLE]

where the dot means time derivative and the prime radial derivative and

[TABLE]

This yields

[TABLE]

and

[TABLE]

This last equation is a special case of Bessel’s equation [3]. We transform it to a self-adjoint form

[TABLE]

so that

[TABLE]

A fundamental system of solutions is given by

[TABLE]

where $J_{0}$ and $Y_{0}$ are the Bessel functions of first and second kind.

The integration constants in (2.20) are determined by the two boundary conditions

[TABLE]

which lead to the self-adjoint boundary conditions

[TABLE]

The first one yields

[TABLE]

and the second one gives

[TABLE]

Using $J^{\prime}_{0}=J_{1}$ and the same for $Y^{\prime}_{0}$ we get a transcendental equation for eigenvalues $k_{n}$

[TABLE]

The corresponding eigenfunctions are equal to

[TABLE]

By general theorems about self-adjoint eigenvalue problems [3] the eigenfunctions form a complete orthogonal system in the Hilbert space $L^{2}([r_{0},r_{1}])$. The normalization integral is easily calculated ([2], p.485)

[TABLE]

We denote the normalized eigenfunctions by $\varphi_{n}(r)$. The solutions $k_{n}$ of (2.23) are approximately given by

[TABLE]

where $\alpha_{n}$,$n=0,1,2\ldots$ are the zeros of the Bessel function $J_{0}(z)$. This follows from the fact that the left-hand side of (2.23) assumes arbitrary values in the neighborhood of $\alpha_{n}$, and then can be made equal to the right-hand side. By the oscillation theorem [4] $\varphi_{n}(r)$ has $n$ zeros in the interval $[r_{0},r_{1}]$.

Now we are ready to write down the general solution of our ablation problem. The time dependent factor $T_{3}(t)$ follows from (2.16): $\sim\exp(\omega_{n}^{\pm}t)$ where $\omega_{n}^{\pm}$ are the solutions of the quadratic equation

[TABLE]

namely

[TABLE]

Then the total solution is equal to

[TABLE]

[TABLE]

The unknown coefficients $a_{n}^{\pm}$ are determined by the two initial conditions

[TABLE]

In previous studies [1] the simple initial condition $f(r)=0$ has been assumed without any justification. However, since the heat generation $Q(r)$ (2.5) is $r$-dependent we expect $f(r)\neq 0$. Our experiments clearly show this (see next section). Using the two conditions (2.28) in (2.27) we obtain the two equations

[TABLE]

[TABLE]

Since the eigenfunctions $\varphi(r)$ form a complete orthonormal system the functions on the right side of (2.29) (2.30) can be uniquely expanded in the usual way by calculating scalar products, for example

[TABLE]

and similarly for (2.29). This finally allows to determine the coefficients $a_{n}^{+}$ and $a_{n}^{-}$ separately.

The form (2.27) of the solution is usefull to study the behavior for large and intermediate times $t$. Since the lowest eigenvalue $k_{0}$ is positive, $\omega_{0}^{+}$ in (2.29) is negative. Consequently, for $t\to\infty$ the solution goes exponentially towards the steady state. In the application to ablation in medicine one is interested in small $t$, as we will see in the next section. The analytic solution above may also be useful for testing numerical codes. But for small $t$ a different treatment must be used.

3 Experimental study of the axisymmetrical ablation problem

We have taken a stainless steel pot with radius $r_{1}=11$ cm with an isolating plate at the bottom. At the center $r=0$ we have placed a bar electrode of length $L=8.5$ cm and small radius $r_{0}=0.45$ cm. Filling the pot with water plus $0.3\%$ NaCl and applying a 13 Volt AC voltage of 50 Hz we have obtained a fairly accurate axisymmetric electric field as assumed in the last section. We measure the temperature in steps of 0.1 degree Celsius at various places $r$ as a function of time.

With help of a voltmeter we have first checked the radial dependence (2.3) of the AC potential. We found exact agreement with (2.3) for $r$ between 4 and 5 cm. For smaller and larger $r$ there are deviations of unknown origin. So we have restricted the temperature measurements to $r=4$ cm and 5 cm.

For each $r$ we have fitted the rise of the temperature $T(t)$ by means of a quadratic polynomial

[TABLE]

Let us try to describe the measured $T(t,r)$ by the usual first order heat equation

[TABLE]

In this case we have at $t=0$

[TABLE]

Calculating $Q$ from equation (2.5) there is a serious problem. The applied electric power $P_{0}$ which can be measured is only partly used for heating. Most of it goes into the electrolysis of Na Cl. In addition there appear bubbles of gas which also costs energy or power. In electrochemistry the heating power $P$ in (2.5) is written as

[TABLE]

where $V_{0}$ is called overpotential and $I$ is the total current [5]. Since $V_{0}$ cannot be directly measured we treat it as a parameter to be determined.

Differentiating (3.2) we get for $t=0$

[TABLE]

[TABLE]

where (2.5) has been inserted. Since this is equal to $2a_{2}$ we can eliminate the unknown $\beta$

[TABLE]

This relation is strongly violated by the measurements which shows that the simple heat equation (3.2) cannot be used for small times.

We proceed similarly with the hyperbolic heat equation

[TABLE]

For $t=0$ we have

[TABLE]

or

[TABLE]

The heat generation $Q(r)$ is given by (2.5)

[TABLE]

where $L$ and $r$ are measured in cm and $P$ in watts (4.185 J = 1 cal). The relevant power $P=V_{0}I$ (3.4) cannot be measured, therefore, we write equ. (3.8) as

[TABLE]

with

[TABLE]

where $P$ is the unknown heating power in watts. The factor 60 is inserted because $a_{1}$ in (3.10) is measured in degree Celsius/min. Now if we calculate $a_{1}(r)$ and $a_{2}(r)$ for two radial distances $r=\varrho_{1}$ and $\varrho_{2}$ from the temperature measurements, then we can determine $\tau$ and $y$.

The determination of $a_{1}$ and $a_{2}$ by a least square fit of the polynomial (3.1) to measured temperature values $T(t_{j})$ is uncertain, depending on the number $N$ of values $t_{j}$ taken with. With our simple equipment (not very good thermometer) we have found for $N=4$

[TABLE]

[TABLE]

for $N=5$:

[TABLE]

[TABLE]

and for $N=6$:

[TABLE]

[TABLE]

We see that $r=4$ cm is a critical distance where $a_{2}$ changes sign and more accurate measurements are necessary.

From the above results, equation (3.10) taken with $r=\varrho_{1}$ and $\varrho_{2}$ gives the following values for the interesting quantities $\tau$ and $y$:

[TABLE]

[TABLE]

[TABLE]

We see that the results obtained are quite uncertain due to the low quality thermometers which did not have the required sensitivity of 0.1 degree. The most reliable values are those of $N=6$ because they contain most measurements. From (3.12) we get $y=0.269P$. Then the electric heating power is $P=y/0.269$ watts. For $y=1.31$ this gives 4.9 watts compared to 13 watts total applied power. This shows that most energy goes into the electrolysis and the bubble generation of gas at the central electrode.

For the medical application the value of the thermal relaxation time $\tau=6.5$ min is most important. In radiofrequency ablation one works with much smaller times $t=0.3-0.4$ min. On the one hand this makes life complicated because the hyperbolic heat equation must certainly be used. On the other hand this equation can simply be solved by a power series (3.1) in $t$. As we have seen the first two terms only depend on the heating power $Q(r)$. The heat conduction contributes to the cubic term

[TABLE]

If one has determined an approximate value of

[TABLE]

one is able to calculate a correction

[TABLE]

for heat conduction.