Finite Difference Modeling of Ultrasound Tissue Heating
Pennes Bioheat Transfer in Cylindrical Coordinates
This is content was originally published as an appendex in my PhD 2003 dissertation “The Theoretical and Experimental Development of Interstitial Ultrasound Applicators for Coagulative Thermal Therapy.”
Heat Transfer in Biological Tissue
A working model of heat transfer in living tissue or ex vivo tissue has proven extremely useful for testing design parameters and estimating performance of interstitial ultrasound applicators. Several mathematical models have been proposed to predict heat transfer through tissue, but one of the most commonly used is the Pennes Bioheat Transfer equation. (Pennes 1948)
\[ \rho _{t}c_{t}\frac{\partial T}{\partial t}=\nabla \cdot k\nabla T+\dot{m} _{b}c_{b}\left( T_{b}-T\right) +\dot{q} \tag{1}\] The terms \(\rho _{t}\) and \(c_{t}\) are the density \([\mathrm{kg}\,\mathrm{m}^{-3}]\) and heat capacity \([\mathrm{J}\,\mathrm{kg}^{-1}\,\mathrm{°C}^{-1}]\) of the tissue, respectively. These are assumed to be independent of the partial derivative of time. The tissue temperature \([\mathrm{°C}]\) is given by \(T\), and the temperature of the perfusing blood is \(T_{b}\). The \(c_{b}\) term is the heat capacity of the blood \([\mathrm{J}\,\mathrm{kg}^{-1}\,\mathrm{°C}^{-1}]\), the \(\omega _{b}\) term is the mass flow rate of perfusing blood per unit tissue volume \([\mathrm{kg}\mathrm{m}^{-3}\mathrm{s}^{-1}]\), and the \(k\) term is the conductive heat transfer coefficient \([\mathrm{W}\mathrm{m}^{-1}\mathrm{K}^{-1}]\) of the tissue. Finally, \(\dot{q}\) represents the external power deposited in the tissue per unit volume \([\mathrm{W}\mathrm{m}^{-3}]\) from the ultrasound applicator.
For annular ultrasound transducers, Equation 1 is most conveniently studied in cylindrical coordinates. In the current form, the only term in Equation 1 affected by the change to cylindrical coordinates is the expression for three-dimensional thermal conduction:
\[ \begin{aligned} \nabla \cdot k\nabla T &=\nabla \cdot \left( k\frac{\partial T}{\partial r}+ \frac{k}{r}\frac{\partial T}{\partial \theta }+k\frac{\partial T}{\partial z}\right) \\ &=\frac{1}{r}\frac{\partial }{\partial r}\left( kr\frac{\partial T}{\partial r}\right) +\frac{1}{r^{2}}\frac{\partial }{\partial \theta }\left( k\frac{\partial T}{\partial \theta }\right) +\frac{\partial }{\partial z}\left( k\frac{\partial T}{\partial z}\right) \end{aligned} \tag{2}\]
Equation 2 gives the most general case in which \(k\) is a function of spatial position. Substituting Equation 2 into Equation 1 gives the Pennes Bioheat Transfer Equation in cylindrical coordinates.
\[ \rho _{t}c_{t}\frac{\partial T}{\partial t}=\frac{1}{r}\frac{\partial }{% \partial r}\left( kr\frac{\partial T}{\partial r}\right) +\frac{1}{r^{2}}% \frac{\partial }{\partial \theta }\left( k\frac{\partial T}{\partial \theta }% \right) +\frac{\partial }{\partial z}\left( k\frac{\partial T}{\partial z}% \right) +\dot{m}_{b}c_{b}\left( T_{b}-T\right) +\dot{q} \tag{3}\]
Ultrasound Power Deposition
When electrical power at the appropriate frequency is applied to a cylindrical ultrasound transducer, acoustic energy is generated. This energy propagates radially out from the transducer and deposits in the surrounding tissue. For an annular ultrasound source, the power deposition \(\dot{q}\) at a position \(r\) away from the source is given as (Diederich and Hynynen 1989).
\[ \dot{q}(r)=\frac{2\alpha I_{r}r_{o}}{r} \]
where $$ is the ultrasound absorption coefficient \([\mathrm{Np}\,\mathrm{m}^{-1}]\) , \(I(r)\) is the ultrasound intensity as a function of radial distance from the transducer \([\mathrm{W}\,\mathrm{m}^{-2}]\), \(r\) is the radial distance from the center of the transducer \([\mathrm{m}]\), and \(r_{o}\) is the radius of the transducer \([\mathrm{m}]\). The intensity term is given by
\[ I(r)=I_{o}e^{-2\mu \left( r-r_{transducer}\right) } \]
where \(I_{o}\) is the acoustic intensity at the transducer surface \([\mathrm{W}/ \mathrm{m}^{2}]\), and $$ is the ultrasound attenuation coefficient \([\mathrm{Np}\, \mathrm{m}^{-1}]\). The absorption coefficient $$ is actually a component of the attenuation coefficient $$ that represents only the absorbed energy. Frequently it is assumed that ultrasound energy that is scattered is absorbed sufficiently quickly that $$ can be approximated by $=$. However, it is possible to have effective attenuation mechanisms that do not result in a temperature rise, so the absorption coefficient could also be less than the attenuation coefficient.
The surface acoustic intensity is given by
\[ I_{o}=\frac{P\nu }{A}=\frac{P\nu }{2\pi r_{o}l} \]
where \(P\) is the applied electrical power \([\mathrm{W}]\), $$ is the acoustic efficiency, and \(A\) is the transducer outer surface area \([\mathrm{m}^{2}]\), and \(l\) is the transducer length \([\mathrm{m}]\). The final form of the power deposition term in cylindrical coordinates for the surrounding tissue is then given as follows:
\[ \dot{q}_{tissue}(r)=\frac{P\nu \alpha }{\pi lr}e^{-2\mu \left( r-r_{0}\right) } \tag{4}\]
Electrical energy applied to the transducer that is not converted to acoustic energy is converted to heat inside the transducer. Within the transducer, \(\dot{q}\) is given by (Diederich et al. 1996)
\[ \dot{q}_{transducer}(r)=\left( 1-\nu \right) \frac{P}{\pi \left( r_{0}^{2}-r_{inner}^{2}\right) } \tag{5}\]
where \(r_{inner}\) is the inner radius of the transducer.
For an annular transducer with air on its inner surface such as the ones frequently used for ultrasound thermal therapy, no acoustic power will radiated towards the center of the transducer. The impedence mismatch between the piezo-ceramic of the transducer and air cause all energy at this surface to be reflected back into the transducer. Therefore, in final form, energy deposition in the model can be written as
\[ \dot{q}(r) = \begin{cases} 0 & \text{for } r < r_{\text{inner}} \\ (1 - \nu) \frac{P}{\pi (r_o^2 - r_{\text{inner}}^2)} & \text{for } r_{\text{inner}} < r < r_o \\ \frac{P\nu\alpha}{\pi l r} e^{-2\mu (r - r_0)} & \text{for } r > r_o \end{cases} \tag{6}\]
Finite Differencing
Explicit Methods
Inserting Equation 6 into Equation 3 with appropriate boundary conditions gives a mathematical model for power deposition from a transducer in an attenuating medium. In theory it is possible to integrate this partial differential equation over a volume and in time to determine how the temperature of any point in that volume develops. In practice terms like attenuation, heat capacity, and thermal conductivity are functions of position, hence analytical solutions are difficult, but numerical techniques can be employed to successfully integrate the system.
Two techniques commonly used for integration of systems of partial differential equations are the method of finite differences and finite elements. Both techniques involve dividing the system into small nodal compartments collectively known as a mesh. Each technique has its advantages, but generally implementation of finite difference solutions is more straightforward when appropriate coordinate systems can be chosen to simplify the system geometry. The primary disadvantage of finite difference systems is that the spacing between nodes cannot vary arbitrarily at positions in the mesh. Generally, systems involving multiple complex geometric shapes require a large number of nodes and a long time to integrate in a finite difference implementation. Finite element solutions can utilize non-uniform meshes to overcome these problems, but generation of these meshes and their numerical solutions are often complicated.
The solution of the temperature profiles generated from a single cylindrical applicator in a volume is easily manageable by finite difference methods. To illustrate how finite differencing can be applied to this problem, it is convenient to start from a simplified version of the model. Let us assume that the only thermal conduction is in the radial direction, making the partial derivatives in $$ and \(z\) zero. Let us also assume that the perfusion and applied energy are zero. This leaves
\[ \rho _{t}c_{t}\frac{\partial T}{\partial t}=\frac{1}{r}\frac{\partial }{% \partial r}\left( kr\frac{\partial T}{\partial r}\right) =\frac{1}{r}\frac{% \partial kr}{\partial r}\frac{\partial T}{dr}+k\frac{\partial ^{2}T}{% \partial r^{2}} \tag{7}\]
Finally, we assume that $kk( r) $ , allowing Equation 7 to become
\[ \rho _{t}c_{t}\frac{\partial T}{\partial t}=\frac{k}{r}\frac{\partial T}{dr} +k\frac{\partial ^{2}T}{\partial r^{2}} \tag{8}\]
The essence of finite differencing is estimating the continuous partial derivatives as discrete differences between points very close together in space or time. In the case of Equation 8, a differenced system would consist of a number of discrete points along a line extending radially from the center of a cylindrical space. The general formulas for explicitly converting a partial differential equation into a finite difference equation are given below.(Press 1997)
\[ \begin{aligned} \frac{\partial u}{\partial t} &=\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t} \\ \frac{\partial u}{\partial r} &=\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2\Delta r} \\ \frac{\partial ^{2}u}{\partial r^{2}} &=\frac{ u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\left( \Delta r\right) ^{2}} \end{aligned} \tag{9}\] The subscript \(i\) refers to the \(i^{th}\) node along the system of radial points. The superscript \(n\) refers to the \(n^{th}\) time step in the simulation. The value \(\Delta t\) is the time between time step \(n\) and \(n+1\). The value \(\Delta r\) is the distance between nodes \(i\) and \(i+1\). These are assumed to be constant for between all \(n\) and \(i\) points, respectively. The value \(r_{i}\) is the radial distance from the center of the system to point \(i\). Rewriting Equation 8 using the rules in Equation 9 gives the following relationship:
\[ \rho _{t}c_{t}\left( \frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t}\right) =\frac{k}{ r_{i}}\left( \frac{T_{i+1}^{n}-T_{i-1}^{n}}{2\Delta r}\right) +k\left( \frac{ T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}}{\left( \Delta r\right) ^{2}}\right) \tag{10}\]
Equation 10 can be solved for the next time step temperature as follows:
\[ T_{i}^{n+1}=\frac{k\Delta t}{r_{i}\rho _{t}c_{t}}\left( \frac{% T_{i+1}^{n}-T_{i-1}^{n}}{2\Delta r}\right) +\frac{k\Delta t}{\rho _{t}c_{t}}% \left( \frac{T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}}{\left( \Delta r\right) ^{2}}% \right) +T_{i}^{n} \tag{11}\]
Equation 11 appears to be extremely useful since it gives an explicit value for the temperature at each node \(i\) at time \(n+1\) given simply the temperature at each node at time \(n\). Therefore, it should be possible to iterate through each node and compute its temperature value at the next time step. By deciding on the size of the time step and the distance between each node, the resolution can be controlled so that the model closely approximates a continuous system. However, all explicit differencing methods are limited by numerical stability. Even a correctly differenced explicit equation system will not converge to an answer unless certain criteria on the coefficients in the differenced equations are met. For example, in Equation 11 it is necessary to pick values for \(\Delta t\) and \(\Delta r\) to make the system numerically stable.
It can be shown that explicit finite differencing an equation consisting of only first order partial derivatives is unstable for any size time step or node spacing.(Press 1997) Explicit differencing of second order partial derivatives such as in Equation 8 is inherently more stable. In general, a general second order partial differential equation
\[ \frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left( D\frac{% \partial u}{\partial x}\right) \tag{12}\]
with constant \(D\) simplifies to
\[ \frac{\partial u}{\partial t}=D\frac{\partial ^{2}u}{\partial x^{2}} \tag{13}\]
The stability criterion for a differenced version of Equation 13 is given by (Press 1997)
\[ \frac{2D\Delta t}{\left( \Delta x\right) ^{2}}\leq 1 \]
This specific result leads to a general heuristic for explicit finite differencing: reducing the node spacing \((\Delta r)\) to more closely approximate a continuous system requires a significant reduction of the time spacing \((\Delta t)\) to maintain numerical stability. Simulations that are designed to resolve thermal effects within small regions such as a transducer a fraction of a millimeter thick would require a very small time step and subsequently a large amount of computing time.
Variable Terms in the Derivative
The model proposed in Equation 7 is different from Equation 13 in that the term \(D\) in the derivative is not constant. In Equation 7, \(D=kr\), so the equation must be differenced in an the appropriately centered manner shown as Equation 14
\[ \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=\frac{D_{i+1/2}\left( u_{i+1}^{n}-u_{i}^{n}\right) -D_{i-1/2}\left( u_{i}^{n}-u_{i-1}^{n}\right) }{% \left( \Delta x\right) ^{2}} \tag{14}\]
with the definitions
\[ \begin{aligned} D_{i+1/2} &=\frac{1}{2}\left( D_{i+1}+D_{i}\right) \\ D_{i-1/2} &=\frac{1}{2}\left( D_{i-1}+D_{i}\right) \end{aligned} \tag{15}\]
Following the scheme outlined in Equation 14 to Equation 15, Equation 7 becomes
\[ \frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t}=\left( \frac{1}{\rho _{t}c_{t}r_{i}}% \right) \frac{\left( kr\right) _{i+1/2}\left( T_{i+1}^{n}-T_{i}^{n}\right) -\left( kr\right) _{i-1/2}\left( T_{i}^{n}-T_{i-1}^{n}\right) }{\left( \Delta r\right) ^{2}} \]
with a stability criterion of
\[ \Delta t\leq \min_{i}\left[ \frac{\left( \Delta x\right) ^{2}}{2D_{i+1/2}}% \right] \tag{16}\]
Variable Node Spacing
With a time step determined by stability requirements and a node spacing determined by resolution requirements, simulations can become prohibitively memory and processor intensive because the number of nodes becomes large and the time step becomes small. One way to reduce the size of the mesh is to only use close grid spacing where it is needed for resolution and decrease the spacing where it is not. In the case of a transducer model in tissue, it is necessary to finely resolve power deposition and heat transfer in and around the transducer, but a courser grid spacing is acceptable far away from the transducer where gradients are less steep. For instance, assuming at least 10 nodes are required in the half millimeter transducer thickness, a node spacing of \(0.05\mathrm{mm}\) would be necessary. Further assuming the model should simulate to outer tissue boundary \(6\mathrm{cm}\) away, 1200 nodes would be required for just a 1-D simulation. However, a centimeter away from the transducer, a \(0.5\mathrm{mm}\) would probably have been sufficient, translating to only 120 nodes and a tenfold decrease in computing time.
The best solution to this problem is to use variable node spacing with tight node spacing close to the transducer and looser spacing some distance away. Given a function to determine the node spacing that is smooth and continuous the finite differencing scheme can be modified as follows: (Kálnay De Rivas 1972),(Quintenz 1978)
\[ \begin{aligned} \frac{\partial u}{\partial r} &=\frac{\left( \Delta r\right) _{+}^{2}\left( u_{i}^{n}-u_{i-1}^{n}\right) +\left( \Delta r\right) _{-}^{2}\left( u_{i+1}^{n}-u_{i}^{n}\right) }{\left( \Delta r\right) _{+}\left( \Delta r\right) _{-}\left[ \left( \Delta r\right) _{+}+\left( \Delta r\right) _{-} \right] } \\ \frac{\partial ^{2}u}{\partial r^{2}} &=\frac{2\left( \Delta r\right) _{+}\left( u_{i-1}^{n}-u_{i}^{n}\right) +2\left( \Delta r\right) _{-}\left( u_{i+1}^{n}-u_{i}^{n}\right) }{\left( \Delta r\right) _{+}\left( \Delta r\right) _{-}\left[ \left( \Delta r\right) _{+}+\left( \Delta r\right) _{-} \right] } \end{aligned} \tag{17}\]
Here the subscripts to the \(r\) spacings refer to the spacing from the positive and negative \(i\) direction. Variable node spacing helps to speed up model computation time by reducing the number of nodes while still maintaining the necessary spatial resolution. Nevertheless, the tight node spacing near the transducer still produces rigid stability constraints on the time step as governed by Equation 16. To further speed up the simulation, a method to increase the time step is required in addition to reducing the average node spacing.
Implicit Differencing
While the differencing techniques discussed so far have been forward in time or explicit, it is also possible to difference in backward time or implicitly. This can be accomplished by modifying the first and second-order portions of Equation 9 as follows:
\[ \begin{aligned} \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t} &=\left[ \frac{ u_{i+1}^{n+1}-u_{i-1}^{n+1}}{2\left( \Delta x\right) }\right] \\ \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t} &=\left[ \frac{ u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1}}{\left( \Delta x\right) ^{2}}\right] \end{aligned} \]
Since each \(u_{i}^{n+1}\) is directly or inderectly a function of all other \(% u^{n+1}\) at all other values of \(i\), the entire solution at each time step is a system of linear simultaneus equations. Even though solving the simultaneous equations at each step takes additional time, the important feature of fully implicit differencing schemes is that they are unconditionally stable. The correct equilibrium solution is always reached, regardless of the spatial or time step, which can be large to speed up simulation. The sacrifice is that the small-scale fluctuations that can drive explicit methods to instability are ignored in order to converge to the large-scale equilibrium solution. In a dynamic simulation where small-scale variations in temperature can effect parameters for the next time step, this is unacceptable.
Mixed Implicit and Explicit
One method that takes advantage of the stability of implicit methods and the accuracy of explicit methods involves taking some arbitrary weighted average of the estimate from each technique. For example, the Crank-Nicholson method takes a simple average of the explicit and implicit estimates to get the next time step. The Crank-Nicholson method of differencing Equation 13 would appear as follows:
\[ \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=\frac{1}{2}D\left( \left[ \frac{ u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1}}{\left( \Delta x\right) ^{2}}\right] +\left[ \frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\left( \Delta x\right) ^{2} }\right] \right) \]
By incorporating variable node spacing and no longer assuming a constant \(D\) , a very robust form for differencing Equation 12 can be written as
\[ \begin{aligned} \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t} &= \frac{(\Delta x_{+})(D_{-1/2}) \left[ u_{i-1}^{n+1}-u_{i}^{n+1}+u_{i-1}^{n}-u_{i}^{n} \right] }{(\Delta x_{+})(\Delta x_{-})(\Delta x_{+} + \Delta x_{-})} \\ &\quad + \frac{(\Delta x_{-})(D_{+1/2}) \left[ u_{i+1}^{n+1}-u_{i}^{n+1}+u_{i+1}^{n}-u_{i}^{n} \right] }{(\Delta x_{+})(\Delta x_{-})(\Delta x_{+} + \Delta x_{-})} \end{aligned} \tag{18}\] Like with the fully implicit method, Equation 18 are stable for any timestep \(\Delta t\), but also have improved spatial and temporal accuracy.
Differenced Bioheat Transfer Equations
The techniques discussed above can easily be implemented to use finite differences to integrate an ordinary differential equation in one dimension. However, the full cylindrical model shown in Equation 3 is a partial differential equation in which the temperature of any finite volume depends on the temperature in all surrounding volumes. This would require solving a huge number of simultaneous equations at each time step even for just explicit differencing. A better solution is to apply the Crank-Nicholson method of finite differencing to each axis in the cylindrical geometry independently and use then numerical tricks to recombine the system.
\(R\)-axis of thermal model
By ignoring all generation terms and thermal conduction along all axes except the \(r\)-axis, the thermal conduction portion of the model can be written as follows:
\[ \begin{aligned} \rho c\frac{\partial T}{\partial t} &=\frac{1}{r}\frac{\partial }{\partial r% }\left( rk\frac{\partial T}{\partial r}\right) \\ \frac{\partial T}{\partial t} &=\frac{1}{\rho cr}\frac{\partial }{\partial r% }\left( rk\frac{\partial T}{\partial r}\right) \\ \frac{\partial T}{\partial t} &=\kappa \frac{\partial }{\partial r}\left( \beta \frac{\partial T}{\partial r}\right) \end{aligned} \tag{19}\] where
\[ \begin{aligned} \kappa &=\frac{1}{\rho cr} \\ \beta &=rk \end{aligned} \]
The Crank-Nicholson differencing scheme for Equation 19 using variable spacing along the \(r\)-axis is given as
\[ \begin{aligned} \frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t} &= \kappa_{i} \frac{(\Delta r_{i-})(\beta_{i+1/2}) \left[ T_{i+1}^{n+1}-T_{i}^{n+1}+T_{i+1}^{n}-T_{i}^{n} \right]}{\Delta r_{i+} \Delta r_{i-} (\Delta r_{i+} + \Delta r_{i-})} \\ &\quad + \kappa_{i} \frac{(\Delta r_{i+})(\beta_{i-1/2}) \left[ T_{i-1}^{n+1}-T_{i}^{n+1}+T_{i-1}^{n}-T_{i}^{n} \right]}{\Delta r_{i+} \Delta r_{i-} (\Delta r_{i+} + \Delta r_{i-})} \end{aligned} \tag{20}\]
Using
\[ \begin{aligned} \delta _{i-} &=\frac{1}{\left( \Delta r_{i-}\right) \left( \Delta r_{i+}+\Delta r_{i-}\right) } \\ \delta _{i+} &=\frac{1}{\left( \Delta r_{i+}\right) \left( \Delta r_{i+}+\Delta r_{i-}\right) } \end{aligned} \]
and simplifying gives the following expression:
\[ \begin{aligned} &(\Delta t \kappa_i \delta_{i-} \beta_{i-1/2}) T_{i-1}^{n+1} + (1 + \Delta t \kappa_i \delta_{i+} \beta_{i+1/2} + \Delta t \kappa_i \delta_{i-} \beta_{i-1/2}) T_{i}^{n+1} \\ &\quad - (\Delta t \kappa_i \delta_{i+} \beta_{i+1/2}) T_{i+1}^{n+1} = (\Delta t \kappa_i \delta_{i-} \beta_{i-1/2}) T_{i-1}^{n} \\ &\quad + (1 - \Delta t \kappa_i \delta_{i+} \beta_{i+1/2} - \Delta t \kappa_i \delta_{i-} \beta_{i-1/2}) T_{i}^{n} + (\Delta t \kappa_i \delta_{i+} \beta_{i+1/2}) T_{i+1}^{n} \end{aligned} \tag{21}\]
Writing the equation in this form has a number of advantages for rapid numerical solution. These are more obvious through the following constant substitutions
\[ \begin{aligned} a_{i} &=\Delta t\kappa _{i}\delta _{i-}\beta _{i-1/2} \\ b_{i} &=-\Delta t\kappa _{i}\delta _{i+}\beta _{i+1/2} \\ c_{i} &=1+a_{i}-b_{i} \\ d_{i} &=a_{i}T_{i-1}^{n}+\left( 1+b_{i}-a_{i}\right) T_{i}^{n}-b_{i}T_{i+1}^{n} \end{aligned} \tag{22}\]
to give \[ \begin{aligned} a_{i}T_{i-1}^{n+1}+\left( 1-b_{i}+a_{i}\right) T_{i}^{n+1}+b_{i}T_{i+1}^{n+1} &=a_{i}T_{i-1}^{n}+\left( 1+b_{i}-a_{i}\right) T_{i}^{n}+b_{i}T_{i+1}^{n} \\ a_{i}T_{i-1}^{n+1}+c_{i}T_{i}^{n+1}+b_{i}T_{i+1}^{n+1} &=d_{i} \end{aligned} \tag{23}\]
Equation 23 represents a system of simultaneous equations that can neatly be written in matrix format as follows:
\[ \begin{bmatrix} c_{1} & b_{1} & 0 & 0 & \cdots & 0 & 0 & 0 \\ a_{2} & c_{2} & b_{2} & 0 & \cdots & 0 & 0 & 0 \\ 0 & a_{3} & c_{3} & b_{3} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & a_{n-1} & c_{n-1} & b_{n-1} \\ 0 & 0 & 0 & 0 & \cdots & 0 & a_{n} & c_{n} \end{bmatrix} \left[ \begin{array}{c} T_{1} \\ T_{2} \\ T_{3} \\ \vdots \\ T_{n-1} \\ T_{n}% \end{array}% \right] =% \begin{bmatrix} d_{1} \\ d_{2} \\ d_{3} \\ \vdots \\ d_{n-1} \\ d_{n}% \end{bmatrix} \tag{24}\] The matrix show as Equation 24 is a special form of sparse matrix known as tridiagonal which can be inverted rapidly using known numerical linear algebra methods. By posing the system in this manner it can be integrated rapidly despite the implicit component of the differencing scheme. Additional development must be done to impose realistic boundary conditions on the model. As written above, the model has a constant zero-degree boundary on both ends of the radial axis.
\(Z\)-axis of thermal model
Similar to the case for the \(r\)-axis, ignoring all generation terms and thermal conduction along all axes except the \(z\)-axis gives the following model:
\[ \rho c\frac{\partial T}{\partial t}=\frac{\partial }{\partial z}\left( k% \frac{\partial T}{\partial z}\right) \]
Crank-Nicholson differentiation with variable node spacing gives
\[ \begin{aligned} \frac{T_{j}^{n+1}-T_{j}^{n}}{\Delta t} &= \frac{1}{\rho_{j}c_{j}} \frac{(\Delta z_{j-})(k_{j+1/2}) \left[ T_{j+1}^{n+1}-T_{j}^{n+1}+T_{j+1}^{n}-T_{j}^{n} \right]}{\Delta z_{j+} \Delta z_{j-} (\Delta z_{j+} + \Delta z_{j-})} \\ &\quad + \frac{1}{\rho_{j}c_{j}} \frac{(\Delta z_{j+})(k_{j-1/2}) \left[ T_{j-1}^{n+1}-T_{j}^{n+1}+T_{j-1}^{n}-T_{j}^{n} \right]}{\Delta z_{j+} \Delta z_{j-} (\Delta z_{j+} + \Delta z_{j-})} \end{aligned} \]
Using
\[ \begin{aligned} \delta _{j-} &=\frac{1}{\left( \Delta z_{j-}\right) \left( \Delta z_{j+}+\Delta z_{j-}\right) } \\ \delta _{j+} &=\frac{1}{\left( \Delta z_{j+}\right) \left( \Delta z_{j+}+\Delta z_{j-}\right) } \end{aligned} \]
and simplifying gives the following expression:
\[ \begin{aligned} &(\Delta t \delta_{j-} k_{j-1/2}) T_{j-1}^{n+1} + (\rho_k c_j + \Delta t \delta_{j+} k_{j+1/2} + \Delta t \delta_{j-} k_{j-1/2}) T_{j}^{n+1} \\[1ex] &\quad - (\Delta t \delta_{j+} k_{j+1/2}) T_{j+1}^{n+1} = (\Delta t \delta_{j-} k_{j-1/2}) T_{j-1}^{n} \\[1ex] &\quad + (\rho_j c_j - \Delta t \delta_{j+} k_{j+1/2} - \Delta t \delta_{j-} k_{j-1/2}) T_{j}^{n} + (\Delta t \delta_{j+} k_{j+1/2}) T_{j+1}^{n} \end{aligned} \]
Like Equation 23 this can be written as a tridiagonal matrix and solved accordingly.
\(\theta\)-axis of thermal model
Carrying out the same procedure for the $$-axis as for the \(r\)- and $% z $-axes gives the following form of the thermal model:
\[ \rho c\frac{\partial T}{\partial t}=\frac{k}{r^{2}}\frac{\partial ^{2}T}{\partial \theta ^{2}} \tag{25}\]
Equation 25 shows a simpler case than the \(r\)-axis because there is no term inside the derivative. In addition, variable spacing along the $% $-axis will not be used, although this could easily be done. Crank-Nicholson differencing gives
\[ \frac{T_{^{k}}^{n+1}-T_{k}^{n}}{\Delta t}=\frac{k_{k}}{2r_{k}^{2}\rho _{k}c_{k}}\left[ \frac{% T_{k-1}^{n}-2T_{k}^{n}+T_{k+1}^{n}+T_{k-1}^{n+1}-2T_{k}^{n+1}+T_{k+1}^{n+1}}{% \left( \Delta \theta \right) ^{2}}\right] \]
Simplification using
\[ \gamma _{k}=\frac{\Delta t}{2\left( \Delta \theta \right) ^{2}}\left( \frac{k}{r^{2}\rho c}\right) _{k} \] gives
\[ \left( -\gamma _{k}\right) T_{k-1}^{n+1}+\left( 1+2\gamma _{k}\right) T_{k}^{n+1}-\left( \gamma _{k}\right) T_{k+1}^{n+1}=\left( \gamma _{k}\right) T_{k-1}^{n}+\left( 1-2\gamma _{k}\right) T_{k}^{n}+\left( \gamma _{k}\right) T_{k+1}^{n} \]
or
\[ e_{k}T_{k-1}^{n+1}+f_{k}T_{k}^{n+1}+e_{k}T_{k+1}^{n+1}=g_{k} \]
with
\[ \begin{aligned} e_{k} &=-\gamma _{k} \\ f_{k} &=1+2\gamma _{k} \\ g_{k} &=\left( \gamma _{k}\right) T_{k-1}^{n}+\left( 1-2\gamma _{k}\right) T_{k}^{n}+\left( \gamma _{k}\right) T_{k+1}^{n} \end{aligned} \tag{26}\]
Since this axis represents a circle, the matrix formulation of the system is slightly different. The first and last nodes overlap, resulting in the following set of simultaneous equations:
\[ \begin{bmatrix} f_{1} & e_{1} & 0 & 0 & \cdots & 0 & 0 & a_{n} \\ e_{2} & f_{2} & e_{2} & 0 & \cdots & 0 & 0 & 0 \\ 0 & e_{3} & f_{3} & e_{3} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & e_{n-1} & f_{n-1} & e_{n-1} \\ e_{1} & 0 & 0 & 0 & \cdots & 0 & e_{n} & f_{n} \end{bmatrix} \left[ \begin{array}{c} T_{1} \\ T_{2} \\ T_{3} \\ \vdots \\ T_{n-1} \\ T_{n}% \end{array} \right] = \begin{bmatrix} g_{1} \\ g_{2} \\ g_{3} \\ \vdots \\ g_{n-1} \\ g_{n} \end{bmatrix} \]
This system, known as cyclically tridiagonal, also has quick inversion methods.
Boundary Conditions
While the Pennes Bioheat Equation models the system behavior in the bulk of the tissue, it is also necessary to describe appropriate boundary conditions for the outer edges of the model and for interfaces within the model. The differencing equations described above are modeled based on conductive boundaries between each node. We will also need to describe constant temperature and convective boundaries for some of the nodes.
Constant Temperature
A simple constant temperature boundary described by \(\frac{\partial T}{\partial t}=0\) can be differenced to give \(T_{i}^{n+1}=T_{i}^{n}\). All energy transferred to this boundary is lost to an infinite constant temperature heat sink.
Insulating
At an insulating boundary, no energy is transferred through. For instance, the case of the \(r\)-axis Equation 19 differenced as Equation 20 with an insulating boundary on the positive direction side of node \(i\), $$\(_{i+1/2}=0\) leaves
\[ \frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t}=\kappa _{i}\frac{\left( \Delta r_{i+}\right) \left( \beta _{i-1/2}\right) \left[ T_{i-1}^{n+1}-T_{i}^{n+1}+T_{i-1}^{n}-T_{i}^{n}\right] }{\left( \Delta r_{i+}\right) \left( \Delta r_{i-}\right) \left( \Delta r_{i+}+\Delta r_{i-}\right) } \]
or
\[ \begin{aligned} &(\Delta t \kappa_{i} \delta_{i-} \beta_{i-1/2}) T_{i-1}^{n+1} + (1 + \Delta t \kappa_{i} \delta_{i-} \beta_{i-1/2}) T_{i}^{n+1} \\ &\quad = (\Delta t \kappa_{i} \delta_{i-} \beta_{i-1/2}) T_{i-1}^{n} + (1 - \Delta t \kappa_{i} \delta_{i-} \beta_{i-1/2}) T_{i}^{n} \end{aligned} \]
Convective
A convective boundary describes the energy lost to a moving fluid along a surface. The total amount of energy moving across a convective boundary is given by
\[ q=hA\left( T-T_{\infty }\right) \]
where \(q\) is energy \([\mathrm{W}]\), \(h\) is the convective heat transfer coefficient \([\mathrm{W}\mathrm{m}^{-2}\mathrm{°C}^{-1}]\), \(T\) is the boundary temperature \([\mathrm{°C}]\), and \(T_{\infty }\) is the temperature of the convecting fluid \([\mathrm{°C}]\). To convert this to cylindrical coordinates, consider the case of water cooling the inner surface of an annulus with inner radius \(r\). The boundary has thickness \(\Delta r\), volume \(V\), surface area \(A\), and length \(l\).
\[ \begin{aligned} \dot{q}& =\frac{q}{V}=\frac{hA\left( T-T_{\infty }\right) }{V}=\frac{2\pi rlh\left( T-T_{\infty }\right) }{\pi l\left( 2r\Delta r+\left( \Delta r\right) ^{2}\right) } \\ & =\frac{2rh\left( T-T_{\infty }\right) }{\Delta r\left( 2r+\Delta r\right) }% \approx \frac{h\left( T-T_{\infty }\right) }{\Delta r} \notag \end{aligned} \tag{27}\]
The equation for a node with energy deposition (or removal) due to conduction on one side and convection on the other along the \(r\)-axis only in cylindrical coordinates is generally described as an extension to Equation 19:
\[ \frac{\partial T}{\partial t}=\kappa \frac{\partial }{\partial r}\left( \beta \frac{\partial T}{\partial r}\right) -\frac{1}{\rho c}\frac{h\left( T-T_{\infty }\right) }{V} \tag{28}\]
Using the results from Equation 27 and differencing like in Equation 20 with $$\(_{i-1/2}=0\) (no conduction in the negative \(r\)-direction due to the convective boundary) gives
\[ \begin{aligned} \frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t} &= \kappa_{i} \frac{\beta_{i+1/2} [T_{i+1}^{n+1}-T_{i}^{n+1}+T_{i+1}^{n}-T_{i}^{n}]}{\Delta r_{i+} (\Delta r_{i+}+\Delta r_{i-})} \\[2ex] &\quad + \frac{h}{2\rho_{i}c_{i}} \frac{[T_{i-1}^{n+1}-T_{i}^{n+1}+T_{i-1}^{n}-T_{i}^{n}]}{\Delta r_{i-}} \end{aligned} \]
Simplification yields
\[ \begin{aligned} &-(\Delta t \kappa_i \delta_{i+} \beta_{i+1/2}) T_{i+1}^{n+1} + \left( 1 + \Delta t \kappa_i \delta_{i+} \beta_{i+1/2} + \frac{h \Delta t}{2\rho_i c_i \Delta r_{i-}} \right) T_{i}^{n+1} \\[1.5ex] &\quad - \left( \frac{h \Delta t}{2\rho_i c_i \Delta r_{i-}} \right) T_{i-1}^{n+1} = (\Delta t \kappa_i \delta_{i+} \beta_{i+1/2}) T_{i+1}^{n} \\[1.5ex] &\quad + \left( 1 - \Delta t \kappa_i \delta_{i+} \beta_{i+1/2} - \frac{h \Delta t}{2\rho_i c_i \Delta r_{i-}} \right) T_{i}^{n} + \left( \frac{h \Delta t}{2\rho_i c_i \Delta r_{i-}} \right) T_{i-1}^{n} \end{aligned} \tag{29}\]
Equation 29 describes a convective boundary node in a format that can easily be inserted into a matrix representing the entire system.
Matrix Formulation of Boundaries
As an example, we will formulate the matrix description of 1-D heat propagation through a cylindrical annulus in which the inner surface is cooled convectively by a fluid at a constant temperature. In this formulation, node 1 will represent the constant temperature fluid boundary, node 2 will represent the convective boundary, and node \(n\) will represent the outer constant temperature boundary.
\[ \begin{bmatrix} 1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ l_{2} & m_{2} & n_{2} & 0 & \cdots & 0 & 0 & 0 \\ 0 & a_{3} & c_{3} & b_{3} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & a_{n-1} & c_{n-1} & b_{n-1} \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1 \end{bmatrix} \left[ \begin{array}{c} T_{1} \\ T_{2} \\ T_{3} \\ \vdots \\ T_{n-1} \\ T_{n} \end{array} \right] = \begin{bmatrix} T_{1} \\ p_{2} \\ d_{3} \\ \vdots \\ d_{n-1} \\ T_{n} \end{bmatrix} \]
where
\[ \begin{aligned} l_{2} &= -(\Delta t \kappa_{2} \delta_{2+} \beta_{2+1/2}) \\ m_{2} &= \left( 1 + \Delta t \kappa_{2} \delta_{2+} \beta_{2+1/2} + \frac{h \Delta t}{2\rho_{2}c_{2}\Delta r_{2-}} \right) \\ n_{2} &= -\left( \frac{h \Delta t}{2\rho_{2}c_{2}\Delta r_{2-}} \right) \\ p_{2} &= (\Delta t \kappa_{2} \delta_{2+} \beta_{2+1/2}) T_{3}^{n} + \left( 1 - \Delta t \kappa_{2} \delta_{2+} \beta_{2+1/2} - \frac{h \Delta t}{2\rho_{2}c_{2}\Delta r_{2-}} \right) T_{2}^{n} \\ &\quad + \left( \frac{h \Delta t}{2\rho_{2}c_{2}\Delta r_{2-}} \right) T_{1}^{n} \end{aligned} \tag{30}\]
and the terms \(a_{i}\), \(b_{i}\), \(c_{i}\), and \(d_{i}\) are given as in Equation 22.
Example: 1-D Model of US Applicator in Tissue
A worthwhile example is to derive the matrix formulation of a model that fully describes heat deposition from an internally cooled, tubular, interstitial ultrasound applicator in homogenous, perfused tissue with conduction in the radial direction only. This approximates the temperature distribution perpendicular to the surface of the middle of a very long, axially symmetric, ultrasound transducer. This simple model also assumes that the properties of tissue such as ultrasound attenuation and perfusion remain constant with temperature.
The model has a convective boundary Equation 28 on the inner surface of the transducer, which is being cooled by constant-temperature water. The outer boundary is held constant at body temperature. Energy deposition in the system is described by Equation 6. The temperature evolution is modeled by Equation 3. Since the energy removal by tissue perfusion and deposition by ultrasound are linear terms in the bioheat equation, so they can simply be added to the right hand side of the difference equation.
The full model is shown below with the constant terms \(a\) through \(p\) above described in equations Equation 22 and Equation 30. The matrix can be divided into 5 distinct layers. Row 1 represents the cooling water held at a constant temperature \(T_{cooling}\). Row 2 describes the conductive boundary on the inner surface of the transducer. Rows 3 through \(s\) describe heat conduction and generation within the transducer. Rows \(s+1\) through \(t-1\) model the heat generation and conduction in the tissue. Finally, row \(t\) describes the outer constant temperature boundary. The interface between the transducer and the tissue (rows \(s\) and \(s+1\)) is modeled as conductively coupled.
\[ \begin{aligned} &\begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0 & \dots & 0 \\ l_{2} & m_{2} & n_{2} & \dots & 0 & 0 & \dots & 0 \\ 0 & a_{3} & c_{3} & \dots & 0 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & a_s & c_s & b_s & 0 & 0 \\ 0 & 0 & 0 & 0 & a_{s+1} & c_{s+1} & b_{s+1} & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & 0 & a_{t-1} & c_{t-1} & b_{t-1} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} T_{1} \\ T_{2} \\ T_{3} \\ \vdots \\ T_{s} \\ T_{s+1} \\ \vdots \\ T_{t-1} \\ T_{t} \end{bmatrix} \\ &= \begin{bmatrix} T_{\text{cooling}} \\ p_{2} + Q_{\text{trans}} \\ d_{3} + Q_{\text{trans}} \\ \vdots \\ d_{s} + Q_{\text{trans}} \\ d_{s+1} + Q_{\text{tissue}}(r_{s+1}) + \text{perf}_{s+1} \\ \vdots \\ d_{t-1} + Q_{\text{tissue}}(r_{t-1}) + \text{perf}_{t-1} \\ T_{\text{bodytemp}} \end{bmatrix} \end{aligned} \tag{31}\]
Time Step Splitting
In multiple dimensions, each temperature point is a function of every other point in the mesh, and it becomes very difficult to solve all of the equations simultaneously. A well established method known as time step splitting can be used to solve each axis individually. In each time step or iteration, the equation set for each coordinate direction is solved independently. In the example of two dimensions, the result after calculating the solution to the columns after a half time step and then calculating the solution to the rows after a half time step using the updated temperature values is a close approximation to the result of calculating the rows and columns simultaneously for a single time step. Specialized methods known as alternating direction implicit (ADI) make this approximation even more accurate.
Considering an equation similar to Equation 12 but in two dimensions
\[ \frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left( D\frac{% \partial u}{\partial x}\right) +\frac{\partial }{\partial y}\left( D\frac{% \partial u}{\partial y}\right) \]
we could employ the notation
\[ \zeta _{x}^{2}u_{i,j,n}=\frac{D_{i+1/2,j,n}\left( u_{i+1,j,n}-u_{i,j,n}\right) -D_{i-1/2,j,n}\left( u_{i,j,n}-u_{i-1,j,n}\right) }{\left( \Delta x\right) ^{2}} \tag{32}\]
by using the results of Equation 14. Therefore, the Crank-Nicholson scheme in two dimension could be written as
\[ \frac{u_{i,j,n+1}-u_{i,j,n}}{\Delta t}=\frac{1}{2}\left( \zeta _{x}^{2}+\zeta _{y}^{2}\right) \left( u_{i,j,n+1}+u_{i,j,n}\right) \tag{33}\]
Imposing an ADI modification on Equation 33 gives
\[ \begin{aligned} \frac{1}{2}\zeta _{x}^{2}\left( u_{i,j,n+1}^{\ast }+u_{i,j,n}\right) +\zeta _{y}^{2}u_{i,j,n} &=\frac{u_{i,j,n+1}^{\ast }-u_{i,j,n}}{\Delta t} \\ \frac{1}{2}\zeta _{x}^{2}\left( u_{i,j,n+1}^{\ast }+u_{i,j,n}\right) +\frac{1% }{2}\zeta _{y}^{2}\left( u_{i,j,n+1}+u_{i,j,n}\right) &=\frac{% u_{i,j,n+1}-u_{i,j,n}}{\Delta t} \end{aligned} \]
where \(u_{n+1}^{\ast }\) is a first approximation of \(u_{i,j,n+1}\) from evaluating the \(x\) derivative only.
In three dimensions, the result is generalized as follows with spatial indices omitted:
\[ \begin{aligned} \frac{1}{2}\zeta _{x}^{2}\left( u_{n+1}^{\ast }+u_{n}\right) +\zeta _{y}^{2}u_{n}+\zeta _{z}^{2}u_{n} &=\frac{u_{n+1}^{\ast }-u_{n}}{\Delta t} \\ \frac{1}{2}\zeta _{x}^{2}\left( u_{n+1}^{\ast }+u_{n}\right) +\frac{1}{2}% \zeta _{y}^{2}\left( u_{n+1}^{\ast \ast }+u_{n}\right) +\zeta _{z}^{2}u_{n} &=\frac{u_{n+1}^{\ast \ast }-u_{n}}{\Delta t} \\ \frac{1}{2}\zeta _{x}^{2}\left( u_{n+1}^{\ast }+u_{n}\right) +\frac{1}{2}% \zeta _{y}^{2}\left( u_{n+1}^{\ast \ast }+u_{n}\right) +\frac{1}{2}\zeta _{z}^{2}\left( u_{n+1}+u_{n}\right) &=\frac{u_{n+1}-u_{n}}{\Delta t} \end{aligned} \]
Finally, variable grid spacing can be used by applying the result of Equation 17 to Equation 32.