**4. Analytical modeling**

In this section, the analytical solution of the DPL model-based bio-heat transfer equation is obtained using the FIT technique. Then, the algorithm for coupling this analytical solution with a solution of RTE is discussed to determine the laser-irradiated biological tissue.

To determine the temperature distribution inside the laser-irradiated biological tissue, the DPL model-based bio-heat transfer equation in terms of temperature is obtained using Eqs. (4) and (10), which can be expressed as [22]:

$$\begin{split} \rho c\_v \mathbf{r}\_q \frac{\partial^2 T}{\partial t^2} + \left(\rho c\_v + \tau\_q a \rho\_b \rho\_b c\_b\right) \frac{\partial T}{\partial t} &= k \left[\nabla \bullet (\nabla T) + \tau\_T \frac{\partial [\nabla \bullet (\nabla T)]}{\partial t}\right] + a \rho\_b \rho\_b c\_b (T\_b - T) \\ &+ \tau\_q a \rho\_b \rho\_b c\_b \frac{\partial T\_b}{\partial t} + Q\_m + \tau\_q \frac{\partial Q\_m}{\partial t} \end{split} \tag{24}$$

Under the assumptions of *Tb* = constant=37 *<sup>o</sup>* C, *Qm* = constant and *θ* defined as *T* � *Tb*, the governing equation (Eq. 24) in the two-dimensional Cartesian coordinate system can be expressed as [22]:

*Modeling of Laser-Irradiated Biological Tissue DOI: http://dx.doi.org/10.5772/intechopen.106794*

$$\rho c\_v \tau\_q \frac{\partial^2 \theta}{\partial t^2} + \left(\rho c\_v + \tau\_q a \rho\_b \rho\_b c\_b\right) \frac{\partial \theta}{\partial t} \quad = k \left[\frac{\partial^2}{\partial \mathbf{x}^2} \left(\theta + \tau\_T \frac{\partial \theta}{\partial t}\right) + \frac{\partial^2}{\partial \mathbf{y}^2} \left(\theta + \tau\_T \frac{\partial \theta}{\partial t}\right)\right] \tag{25}$$
 
$$-a \rho\_b \rho\_b c\_b \theta + Q\_m$$

Eq. (24) can also be expressed in a cylindrical coordinate system for axisymmetric biological tissue by Eq. (26) [23].

$$\left(\rho c\_v \tau\_q \frac{\partial^2 \theta}{\partial t^2} + \left(\rho c\_v + \tau\_q \alpha b\_b \rho\_b c\_b\right) \frac{\partial \theta}{\partial t} = k \left[\frac{\partial^2 \theta}{\partial r^2} + \frac{1}{r} \frac{\partial \theta}{\partial r} + \frac{\partial^2 \theta}{\partial z^2}\right.\tag{26}$$

$$+ \tau\_T \frac{\partial}{\partial t} \left(\frac{\partial^2 \theta}{\partial r^2} + \frac{1}{r} \frac{\partial \theta}{\partial r} + \frac{\partial^2 \theta}{\partial z^2}\right)\right] - \alpha\_b \rho\_b c\_b \theta + Q\_m$$

The FIT technique [24] converts the partial differential equation into the ordinary differential equation by removing the space derivative using the integral transformation pair. Thus, the ordinary differential equation is only the function of time. This equation with the transformed initial conditions is solved. The transformed temperature can be inverted back into the temperature using the inversion formula.

The boundary conditions Eqs. (27)–(30) and initial conditions Eqs. (31) and (32) required for solving Eq. (25) are given below:

$$\theta(\mathbf{0}, \mathbf{y}, t) = \mathbf{0} \tag{27}$$

$$\theta(L, \mathcal{y}, t) = \mathbf{0} \tag{28}$$

$$\theta(x, \mathbf{0}, t) = \mathbf{0} \tag{29}$$

$$\frac{\partial\theta(\mathbf{x},L\_{\mathcal{V}},t)}{\partial\mathbf{y}} + H\theta(\mathbf{x},L\_{\mathcal{V}},t) = H\theta\_{\infty} \tag{30}$$

$$\theta(\mathfrak{x}, \mathfrak{y}, \mathfrak{0}) = \mathfrak{0} \tag{31}$$

$$\frac{\partial \theta(x, y, \mathbf{0})}{\partial t} = \mathbf{0} \tag{32}$$

$$\text{Here } H = \frac{h}{k} \dots$$

The integral transform pair for a given function *θ*ð Þ *x*, *y*, *t* concerning the *x* variable can be expressed as below

Inversion formula:

$$\theta(\mathbf{x}, y, t) = \sum\_{p=1}^{\infty} \frac{X(\nu\_p, \mathbf{x})}{N(\nu\_p)} \overline{\theta}(\nu\_p, y, t) \tag{33}$$

Integral transform:

$$\overline{\theta}(\nu\_p, \mathbf{y}, t) = \int\_{x'=0}^{L} X(\nu\_p, \mathbf{x'}) \theta(\mathbf{x'}, \mathbf{y}, t) d\mathbf{x'} \tag{34}$$

Similarly, the integral transform pair for the function *θ νp*, *y*, *t* � � concerning the *y* variable can be expressed as

Inversion formula:

$$\overline{\theta}(\nu\_p, \jmath, t) = \sum\_{m=1}^{\infty} \frac{Y(\eta\_m, \jmath)}{N(\eta\_m)} \overline{\theta}(\nu\_p, \eta\_m, t) \tag{35}$$

Integral transform:

$$\bar{\overline{\theta}}(\nu\_p, \eta\_m, t) = \int\_{\mathcal{Y}'=0}^{W} Y(\eta\_m, \mathcal{Y}) \overline{\theta}(\nu\_p, \mathcal{Y}, t) d\mathcal{Y}' \tag{36}$$

Here, *<sup>ν</sup>p*,*η<sup>m</sup>* are eigenvalues, *<sup>X</sup> <sup>ν</sup>p*, *<sup>x</sup>* � � , *<sup>Y</sup>*ð Þ *<sup>η</sup>m*, *<sup>y</sup>* are Eigen function and *<sup>N</sup> <sup>ν</sup><sup>p</sup>* � �, *N*ð Þ *η<sup>m</sup>* are the norm. The eigenvalues, Eigen function, and norm for the givens set boundary conditions Eqs. (27)–(30) can be found elsewhere [24].

Multiplying both sides of Eq. (25) by the operatorÐ *<sup>L</sup> x*0 ¼<sup>0</sup>*<sup>X</sup> <sup>ν</sup>p*, *<sup>x</sup>*<sup>0</sup> � �*dx*<sup>0</sup> and using the Eq. (34) and Eqs. (27) and (28), the derivative concerning *x* variables is removed from Eq. (25), and the resulting equation is in terms of *θ νp*, *y*, *t* � �. Similarly, the derivative concerning *y* variables is removed from this transformed governing equation by multiplying both sides by Ð *<sup>W</sup> y*0 ¼<sup>0</sup>*<sup>Y</sup> <sup>η</sup>m*, *<sup>y</sup>*<sup>0</sup> ð Þ*dy*<sup>0</sup> and using the Eq. (36) and the transformed boundary conditions of Eqs. (29) and (30), we get [22]

$$\frac{d^2\bar{\theta}(\nu\_p, \eta\_m, t)}{dt^2} + a\frac{d\bar{\theta}(\nu\_p, \eta\_m, t)}{dt} + b\tilde{\theta}(\nu\_p, \eta\_m, t) = \frac{1}{\rho c\_v \tau\_q} \left( kH\theta\_\infty \sin \eta\_m L\_\gamma \frac{[1 - (-1)^p]}{\nu\_p} \right) \tag{37}$$

$$\text{where } a = \frac{k\_{\text{TT}} \left(\nu\_p^2 + \eta\_m^2\right) + \left(\rho c\_v + \tau\_q \alpha y\_b \rho\_b c\_b\right)}{\rho c\_v \tau\_q}, b = \frac{k \left(\nu\_p^2 + \eta\_m^2\right) + \alpha y\_b \rho\_b c\_b}{\rho c\_v \tau\_q} \text{ .}$$

Eq. (37) is the second-order non-homogenous ordinary differential equation, and its solution using the transformed initial conditions has been obtained, which is given below [22]:

*θ*ð Þ¼ *x*, *y*, *tn*þ<sup>1</sup> X∞ *p*¼1 *Hθ*<sup>∞</sup> <sup>1</sup> � �ð Þ<sup>1</sup> *<sup>p</sup>* ½ � *νp sinh b*1ð Þ *p y b*1ð Þ *p cosh b*1ð Þ *p W* þ *Hsinhhb*1ð Þ *p W sin νpx N ν<sup>p</sup>* � � <sup>þ</sup> <sup>X</sup><sup>∞</sup> *p*¼1 X∞ *m*¼1 ~ *θ νp*, *ηm*, *tn* � � � *kHθ*<sup>∞</sup> sin *<sup>η</sup>mW ρcvτqb* � � <sup>1</sup> � �ð Þ<sup>1</sup> *<sup>p</sup> νp* � � � � *<sup>γ</sup>* <sup>þ</sup> *<sup>ω</sup>*<sup>1</sup> 2*ω*<sup>1</sup> � � *exp* f g �ð Þ *<sup>γ</sup>* � *<sup>ω</sup>*<sup>1</sup> <sup>Δ</sup>*<sup>t</sup>* � � *<sup>γ</sup>* � *<sup>ω</sup>*<sup>1</sup> 2*ω*<sup>1</sup> � � *exp* f g �ð Þ *<sup>γ</sup>* <sup>þ</sup> *<sup>ω</sup>*<sup>1</sup> <sup>Δ</sup>*<sup>t</sup>* � *sin νpx N ν<sup>p</sup>* � � *sin ηmy N*ð Þ *η<sup>m</sup>* <sup>þ</sup> <sup>X</sup><sup>∞</sup> *p*¼1 X∞ *m*¼1 *d*~ *θ νp*, *ηm*, *tn* � � *dt exp* f g �ð Þ *γ* � *ω*<sup>1</sup> Δ*t* � *exp* f g �ð Þ *γ* þ *ω*<sup>1</sup> Δ*t* 2*ω*<sup>1</sup> � � *sin νpx N ν<sup>p</sup>* � � *sin ηmy N*ð Þ *η<sup>m</sup>* , *a* 2 � �<sup>2</sup> >*b*

*Modeling of Laser-Irradiated Biological Tissue DOI: http://dx.doi.org/10.5772/intechopen.106794*

$$\begin{split} \theta(\mathbf{x},\mathbf{y},t\_{n+1}) &= \sum\_{p=1}^{\infty} H\theta\_{\infty} \frac{[1-(-1)^{p}]}{\nu\_{p}} \frac{\sinh b\_{1}(p)y}{b\_{1}(p)\cosh b\_{1}(p)W + H\sinh b h\_{1}(p)W} \frac{\sin \nu\_{p}\mathbf{x}}{N(\nu\_{p})} \\ &+ \sum\_{p=1}^{\infty} \sum\_{m=1}^{\infty} \left[ \tilde{\theta}(\nu\_{p}, \eta\_{m}, t\_{n}) - \left(\frac{kH\theta\_{\infty}\sin\eta\_{m}W}{\rho c\_{1}x\_{0}b}\right) \left[\frac{1-(-1)^{p}}{\nu\_{p}}\right] \right] \left[\frac{\exp\left(-\gamma\Delta t\right)}{o\_{2}}\right] |o\_{2}\cos\omega\_{2}\Delta t| \\ &+ \gamma\sin\alpha\omega\_{2}\Delta t\Big] \frac{\sin\nu\_{p}\mathbf{x}}{N(\nu\_{p})} \frac{\sin\eta\_{m}\mathbf{y}}{N(\eta\_{m})} + \sum\_{p=1}^{\infty} \sum\_{m=1}^{\infty} \frac{\tilde{d\tilde{\theta}}(\nu\_{p}, \eta\_{m}, t\_{n})}{dt} \left[\frac{\sin\alpha\_{2}\Delta t \exp\left(-\gamma\Delta t\right)}{o\_{2}}\right] \frac{\sin\nu\_{p}\mathbf{x}}{N(\nu\_{p})} \frac{\sin\eta\_{m}\mathbf{y}}{N(\eta\_{m})}, \\ & \left(\frac{a}{2}\right)^{2} < b \end{split} \tag{39}$$

where *ω*<sup>1</sup> ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *a* 2 � �<sup>2</sup> � *<sup>b</sup>* q , and *ω*<sup>2</sup> ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *<sup>b</sup>* � *<sup>a</sup>* 2 � �<sup>2</sup> <sup>q</sup>

It is to be noted here that the subscript "*n*" denotes the previous time instant while subscript "*n*+1" represents the current time instant. Here, Δ*t* ¼ *tn*þ<sup>1</sup> � *tn* and ~ *θ νp*, *ηm*, *tn* � �is the temperature distribution at the previous time instant *tn*. If *tn* <sup>¼</sup> 0, Eqs. (31) and (32) is used as the initial conditions.

When multiplying both sides of Eq. (37) by *τ<sup>q</sup>* and then substituting *τ<sup>q</sup>* ¼ *τ<sup>T</sup>* ¼ 0 into Eq. (37), the DPL model-based bio-heat transfer equation becomes the Fourier model-based bio-heat transfer equation, and its solution is expressed by Eq. (40).

$$\begin{aligned} \theta(\mathbf{x}, \mathbf{y}, t\_{n+1}) &= \sum\_{p=1}^{\infty} H \theta\_{\infty} \frac{[1 - (-1)^p]}{\nu\_p} \frac{\sinh b\_1(p) y}{b\_1(p) \cosh b\_1(p) W + H \sinh b h\_1(p) W} \frac{\sin \nu\_p \pi}{N(\nu\_p)} \\ &+ \sum\_{p=1}^{\infty} \sum\_{m=1}^{\infty} \left[ \tilde{\theta} \begin{pmatrix} \nu\_p & \eta\_m & t\_n \end{pmatrix} \right. \\\\ &- \left( \frac{k H \theta\_{\infty} \sin \eta\_m W}{\rho c\_p b} \right) \left[ \frac{1 - (-1)^p}{\nu\_p} \right] \exp \left( -b \Delta t \right) \frac{\sin \nu\_p \pi \sin \eta\_m y}{N(\nu\_p)} \end{aligned} \tag{40}$$

Following the approach discussed in Section 2, the temperature rise due to a single pulse (Δ*T*) was first obtained by solving Eq. (9). After that, this temperature rise has been added to the temperature distribution already obtained for the previous time instant to determine the temperature distribution at the current time instant:

$$
\tilde{\overline{\theta}}(\nu\_p, \eta\_m, t\_n) = \tilde{\overline{\theta}}(\nu\_p, \eta\_m, t\_n) + \Delta \theta(\nu\_p, \eta\_m) \tag{41}
$$

where

$$
\Delta\theta\left(\nu\_p, \eta\_m\right) = \int\_{\mathbf{x'}=0}^{L} \int\_{\mathbf{y'}=0}^{W} \Delta T(\mathbf{x'}, \mathbf{y'}) \sin \nu\_p \mathbf{x'} \sin \eta\_m \mathbf{y'} d\mathbf{x'} d\mathbf{y'} \tag{42}
$$

The trapezoidal rule has been employed for evaluating the double integral appearing on the right-hand side of Eq. (42).

Similarly, the analytical solution of the DPL model-based bio-heat transfer equation in the cylindrical coordinate system Eq. (26) can be obtained using FIT, and the details can be given in elsewhere [23].
