1 Introduction

The goal of biological tissue artificial heating is to destroy the diseased fragments and to minimize the damage of surrounding healthy tissue. Artificial hyperthermia treatments are performed using, among others, the electromagnetic field or the laser action. In this paper, the influence of lasers heating on biological tissue is analyzed. In the planning of this type of treatment, the mathematical methods are used, which allow the laser intensity to be estimated, assuring the destruction of the target region of biological tissue. To analyze the problem discussed, the mathematical models describing the laser irradiation, the temperature distribution in the tissue domain and the tissue damage degree should be formulated (Jacques and Pogue 2008; Niemz 2007; Tuchin 2007; Welch 1984; Welch and van Gemert 2011). In the case of soft tissues, the scattering generally dominates over the absorption for wavelengths between 650 and 1300 nm and to determine the diffuse fluence rate the steady-state optical diffusion equation can be taken into account (Dombrovsky 2012; Farrel et al. 1992; Fasano et al. 2010; Guo et al. 2003). The temperature field in the irradiated biological tissue can be described by the different mathematical models. Three of them result from the theory of continuum mechanics, namely Pennes (1948), Vernotte (1958) and dual-phase lag (Xu et al. 2008) models. These models are widely used in numerical modeling of thermal processes occurring in the laser-irradiated biological tissues, e.g. Afrin et al. 2012; Fanjul-Vélez et al. 2009; He et al. 2006; Hooshmand et al. 2015; Jaunich et al. 2008; Kim and Guo 2007; Narasimhan and Sadasivam 2013; Zhang et al. 2017; Zhou et al. 2009a, b. The second group of models is based on the theory of porous media (Vafai 2010) and then the biological tissue is divided into two regions: vascular region (blood vessel) and extravascular region (tissue). Using this approach, the generalized dual-phase lag model (GDPL) is formulated (Zhang 2009). This model is starting to be used to determine the tissue and blood temperatures in the heated tissues (Afrin et al. 2012; Jasiński et al. 2016; Kałuża et al. 2017; Majchrzak et al. 2015). Knowledge of temperature distribution allows one to estimate the degree of tissue destruction using so-called Arrhenius integral (Jasiński 2013; Niemz 2007) or thermal dose parameter (Sapareto and Dewey 1984).

So far, the laser–tissue interactions have been modeled as the direct problems. For the assumed intensity of the laser, the source function related to the laser heating, the temperature distribution and finally the degree of tissue destruction have been determined. In this paper, the inverse problem connected with the estimation of laser intensity assuring the destruction of target region of biological tissue is considered. At first, the direct problem using the finite difference method is solved. To describe the temperature distribution the generalized dual-phase lag model supplemented by appropriate boundary and initial conditions is applied. The source function in the GDPL equation resulting from the laser heating is associated with the total fluence rate, where its diffuse part is determined using the steady-state optical diffusion equation. The degree of tissue destruction is estimated on the basis of Arrhenius integral. Next, taking into account these values, the inverse problem is solved using the gradient method. In the final part, the results of computations are shown.

2 Formulation of the Direct Problem

An axially symmetrical domain of biological tissue exposed to the laser beam, as shown in Fig. 1, is considered. Thermal processes can be described by generalized dual-phase lag model (Majchrzak et al. 2015; Zhang 2009). This model consists of two coupled equations:

$$\begin{aligned} C\left[ {\frac{\partial T(r,z,t)}{\partial t} + \tau_{q} \frac{{\partial^{2} T(r,z,t)}}{{\partial t^{2} }}} \right] = & \varLambda \nabla^{2} T(r,z,t) + \varLambda \tau_{T} \frac{\partial }{\partial t}\left[ {\nabla^{2} T(r,z,t)} \right] \\ & \quad + \;G\left[ {T_{b} (r,z,t) - T(r,z,t)} \right] + \varepsilon Q_{mb} + \left( {1 - \varepsilon } \right)Q_{{m{\kern 1pt} {\kern 1pt} }} + Q_{ex} (r,z,t) \\ & \quad + \;\frac{{\tau_{q} C}}{{\left( {1 - \varepsilon } \right)\rho c}}\left[ {\varepsilon \frac{{\partial Q_{mb} }}{\partial t} + \left( {1 - \varepsilon } \right)\frac{{\partial Q_{{m{\kern 1pt} {\kern 1pt} }} }}{\partial t} + \frac{{\partial Q_{ex} (r,z,t)}}{\partial t}} \right], \\ \end{aligned}$$
(1)

and

$$T_{b} (r,z,t) = T(r,z,t) - \frac{{\varepsilon \rho_{b} c_{b} }}{G}\frac{{\partial T_{b} (r,z,t)}}{\partial t}$$
(2)

where ε is the porosity defined as the ratio of blood volume to the total volume, \(C = \varepsilon \rho_{b} c_{b} + \left( {1{-}\varepsilon } \right)\rho c\) is the effective heat capacity (ρ, ρb, c, cb are the densities and specific heats of tissue and blood, respectively), \(\varLambda = \varepsilon \lambda_{b} + \left( {1{-}\varepsilon } \right)\lambda\) is the effective thermal conductivity (λ, λb are the thermal conductivities of tissue and blood, respectively), G is the coupling factor, τq is the relaxation time, τT is the thermalization time, Qm and Qmb are the constant metabolic heat sources, Qex(r, z, t) is the source function connected with the laser irradiation and T(r, z, t) and Tb(r, z, t) are the tissue and blood temperatures.

Fig. 1
figure 1

The domain considered

Full size image

The generalized dual-phase lag Eq. (1) is supplemented by boundary condition (Mochnacki and Majchrzak 2017):

$$(r,z) \in \varGamma_{0} \cup \varGamma :\quad - \varLambda \,\left[ {{\mathbf{n}} \cdot \nabla T(r,z,t) + \tau_{T} \frac{{\partial \left[ {{\mathbf{n}} \cdot \nabla T\left( {r,z,t} \right)} \right]}}{\partial t}} \right] = q_{b} \left( {r,z,t} \right) + \tau_{q} \frac{{\partial q_{b} \left( {r,z,t} \right)}}{\partial t}$$
(3)

where Γ0 corresponds to the axis of symmetry, Γ is the boundary of the cylinder, n is the normal outward vector, ∇T(r, z, t) is the temperature gradient and qb(r, z, t) is known boundary heat flux. Here, the no-flux boundary condition qb(r, z, t) = 0 is assumed. It should be noted that for strongly scattering tissues (like the soft tissues analyzed here), laser heating is considered as a body heat source but the irradiated surface is thermally insulated (Afrin et al. 2012). Thus, on the upper surface of the domain considered qb(r, z, t) = 0 is also assumed.

The initial conditions are also given:

$$t = 0:\quad T(r,z,t) = T_{p} ,\quad \left. {\frac{\partial T(r,z,t)}{\partial t}} \right|_{t = 0} = w(r,z)$$
(4)

where Tp is known temperature and w(r, z) is the initial heating rate.

The ordinary differential Eq. (2) is supplemented by condition

$$t = 0:\;\;\;\;\;T_{b} (r,z,t) = T_{p} .$$
(5)

Source function Qex(r, z, t) connected with the laser heating can be defined as follows (Kim et al. 1996):

$$Q_{ex} (r,z,t) = \mu_{a\,} \phi (r,z)p(t)$$
(6)

where μa [1/m] is the absorption coefficient, ϕ(r, z) [W/m2] is the total light fluence rate and p(t) is the function equal to 1 when the laser is on and equal to 0 when the laser is off.

The total light fluence rate ϕ(r, z) is the sum of collimated part ϕc and diffuse part ϕd (Gardner et al. 1996; Guo et al. 2003; Kim et al. 1996)

$$\phi (r,z) = \phi_{c} (r,z) + \phi_{d} (r,z).$$
(7)

In the case of soft tissues, to determine the diffuse fluence rate the steady-state optical diffusion equation (Farrel et al. 1992; Fasano et al. 2010; Welch and van Gemert 2011) should be solved

$$\begin{array}{*{20}c} {\left( {r,z} \right) \in \varOmega :\quad D\nabla^{2} \phi_{d} (r,z) - \mu_{a\,} \phi_{d} (r,z) + \mu_{s}^{{\prime }} \phi_{c} (r,z) = 0} \\ \end{array}$$
(8)

where

$$D = \frac{1}{{3\left[ {\mu_{a} + \left( {1 - g} \right)\mu_{s} } \right]}}$$
(9)

and \(\mu_{s}^{\prime } = (1 - g)\mu_{s}\) [1/m] is the effective scattering coefficient (μs is the scattering coefficient, g is the anisotropy factor).

Equation (8) is supplemented by boundary conditions

$$\begin{aligned} \begin{array}{*{20}c} {\left( {r,z} \right) \in \varGamma :} & { - \,D\,{\mathbf{n}} \cdot \nabla \phi_{d} (r,z) = \frac{{\phi_{d} (r,z)}}{2}} \\ \end{array} \hfill \\ \begin{array}{*{20}c} {\left( {r,z} \right) \in \varGamma_{0} :} & { - D\,{\mathbf{n}} \cdot \nabla \phi_{d} (r,z) = 0} \\ \end{array} . \hfill \\ \end{aligned}$$
(10)

The collimated fluence rate is given as (Fasano et al. 2010)

$$\phi_{c} (r,z) = I\exp ( - \mu_{t}^{\prime } \,z)\exp \left( { - \frac{{r^{2} }}{{r_{D}^{2} }}} \right)$$
(11)

where I [W/m2] is the surface irradiance of laser, rD is the radius of laser beam and \(\mu_{t}^{\prime }\) [1/m] is the attenuation coefficient defined as

$$\mu_{t}^{\prime } = \mu_{a} + \mu_{s}^{\prime } .$$
(12)

Solution of the above-presented problems allows one to determine the Arrhenius integral (Jasiński 2013; Niemz 2007)

$$A\left( {r,z,t^{f} } \right) = P\,\int\limits_{0}^{{t^{f} }} {\exp \left( { - \frac{E}{R\,T(r,z,t)}} \right){\text{d}}t}$$
(13)

where P[1/s] is the pre-exponential factor, E[J/mol] is the activation energy, R[J/(mol K)] is the universal gas constant and [0, tf] is the time interval under consideration.

A value of damage integral A(r, z, tf) = 1 corresponds to a 63% probability of cell death at a specific point (r, z), while A(r, z, tf) = 4.6 corresponds to 99% probability of cell death at this point (Chang and Nguyen 2004).

3 Inverse Problem

The inverse problem formulated here concerns the estimate of laser intensity I which ensures the destruction of target region of biological tissue. Thus, the following criterion can be formulated:

$$S(I) = \sum\limits_{f = 1}^{F} {\sum\limits_{i = 1}^{M} {\left[ {A(r_{i} ,z_{i} ,t^{f} ,I) - A_{m} (r_{i} ,z_{i} ,t^{f} )} \right]^{2} } }$$
(14)

where \(A_{m} (r_{i} ,z_{i} ,t^{f} )\) is the ‘measured’ Arrhenius integral. \(A(r_{i} ,z_{i} ,t^{f} ,I)\) is the calculated Arrhenius integral obtained from the direct problem solution with current estimate of the unknown parameter I, while M is the number of points and F is the number of time steps.

This inverse problem is solved using the gradient method (Kurpisz and Nowak 1995). Thus, using the necessary condition of optimum one has

$$\frac{{{\text{d}}S(I)}}{{{\text{d}}I}} = 2\sum\limits_{f = 1}^{F} {\sum\limits_{i = 1}^{M} {\left[ {A_{i}^{f} - A_{m} (r_{i} ,z_{i} ,t^{f} )} \right]} } \frac{{\partial A(r_{i} ,z_{i} ,t^{f} ,I)}}{\partial I} = 0$$
(15)

where (c.f. Eq. 13)

$$V_{i}^{f} = \frac{{\partial A(r_{i} ,z_{i} ,t^{f} ,I)}}{{\partial {\kern 1pt} I}} = P\int\limits_{0}^{{t^{f} }} {\frac{E}{{RT^{2} (r_{i} ,z_{i} ,t)}}\exp \left( { - \frac{E}{{RT(r_{i} ,z_{i} ,t)}}} \right)U(r_{i} ,z_{i} ,t){\text{d}}t}$$
(16)

while

$$U(r_{i} ,z_{i} ,t) = \frac{{\partial T(r_{i} ,z_{i} ,t)}}{\partial I}$$
(17)

are the sensitivity coefficients and \(A_{i}^{f} = A(r_{i} ,z_{i} ,t^{f} ,I)\).

Function \(A_{i}^{f}\) is expanded into a Taylor series about known value of I k, this means

$$A_{i}^{f} = \left( {A_{i}^{f} } \right)^{k} + \left( {V_{i}^{f} } \right)^{k} \left( {I_{{}}^{k + 1} - I_{{}}^{k} } \right)$$
(18)

where k is the number of iteration and Ik for k = 0 is the arbitrary assumed value of I, while Ik for k > 0 results from the previous iteration.

Introducing formula (18) to Eq. (15) one obtains

$$\sum\limits_{f = 1}^{F} {\sum\limits_{i = 1}^{M} {\left[ {\left( {A_{i}^{f} } \right)^{k} + \left( {V_{i}^{f} } \right)^{k} \left( {I^{k + 1} - I^{k} } \right) - A_{m} (x_{i} ,t^{f} )} \right]\left( {V_{i}^{f} } \right)^{k} = 0} } .$$
(19)

After the mathematical manipulations one has

$$I^{k + 1} = I^{k} + \frac{{\sum\limits_{f = 1}^{F} {\sum\limits_{i = 1}^{M} {\left[ {A_{m} (r_{i} ,z_{i} ,t^{f} ) - \left( {A_{i}^{f} } \right)^{k} } \right]} } \left( {V_{i}^{f} } \right)^{k} }}{{\sum\limits_{f = 1}^{F} {\sum\limits_{i = 1}^{M} {\left[ {\left( {V_{i}^{f} } \right)^{k} } \right]}^{2} } }},\quad k = 0,1,2, \ldots ,K$$
(20)

where K is the assumed number of iterations.

To determine the sensitivity coefficients (17), the governing equations should be differentiated with respect to the unknown parameter I (Kleiber 1997; Majchrzak and Mochnacki 2014). At first, the source function Qex(r, z, t) connected with the laser irradiation is differentiated

$$\frac{{\partial Q_{ex} (r,z,t)}}{\partial I} = \mu_{a\,} \left[ {\frac{{\partial \phi_{d} (r,z)}}{\partial I} + \frac{{\partial \phi_{c} (r,z)}}{\partial I}} \right]p(t)$$
(21)

or

$$\frac{{\partial Q_{ex} (r,z,t)}}{\partial I} = \mu_{a\,} \left[ {W_{d} (r,z) + W_{c} (r,z)} \right]p(t)$$
(22)

where

$$\begin{array}{*{20}c} {W_{d} (r,z) = \frac{{\partial \phi_{d} (r,z)}}{\partial I},\;\;W_{c} (r,z) = \frac{{\partial \phi_{c} (r,z)}}{\partial I}} \\ \end{array} .$$
(23)

Taking into account the dependence (11) one has

$$W_{c} (r,z) = \exp ( - \mu_{t}^{\prime } z)\exp \left( { - \frac{{r^{2} }}{{r_{D}^{2} }}} \right)$$
(24)

while the function Wd(r, z) results from the differentiation of Eq. (8) and boundary conditions (10). Thus

$$\begin{array}{*{20}c} {\left( {r,z} \right) \in \varOmega :\quad D\nabla^{2} W_{d} (r,z) - \mu_{a\,} W_{d} (r,z) + \mu_{s}^{\prime } W_{c} (r,z) = 0} \\ \end{array}$$
(25)

and

$$\begin{array}{*{20}c} {\left( {r,z} \right) \in \varGamma :} & { - D\,{\mathbf{n}} \cdot \nabla W_{d} (r,z) = \frac{{W_{d} (r,z)}}{2}} \\ \end{array}$$
(26)

Next, Eqs. (1) and (2) are differentiated with respect to the parameter I and then

$$\begin{aligned} C\left[ {\frac{\partial U(r,z,t)}{\partial t} + \tau_{q} \frac{{\partial^{2} U(r,z,t)}}{{\partial t^{2} }}} \right] & = \varLambda \nabla^{2} U(r,z,t) + \varLambda \tau_{T} \frac{\partial }{\partial t}\left[ {\nabla^{2} U(r,z,t)} \right] \\ & \quad + G\left[ {U_{b} (r,z,t) - U(r,z,t)} \right] + \frac{{\partial Q_{ex} (r,z,t)}}{\partial I} \\ \end{aligned}$$
(27)
$$U_{b} (r,z,t) = U(r,z,t) - \frac{{\varepsilon \rho_{b} c_{b} }}{G}\frac{{\partial U_{b} (r,z,t)}}{\partial t}$$
(28)

where

$$\begin{array}{*{20}c} {U(r,z,t) = \frac{\partial T(r,z,t)}{\partial I},\;\;U_{b} (r,z,t) = \frac{{\partial T_{b} (r,z,t)}}{\partial I}.} \\ \end{array}$$
(29)

The boundary condition (3) and the initial conditions (4) (5) are also differentiated

$$(r,z) \in \varGamma :\quad - \varLambda \left[ {{\mathbf{n}} \cdot \nabla U(r,z,t) + \tau_{T} \frac{{\partial \left[ {{\mathbf{n}} \cdot \nabla U\left( {r,z,t} \right)} \right]}}{\partial t}} \right] = 0$$
(30)
$$t = 0:\;\;\;\;\;U(r,z,t) = 0,\;\;\;\;\left. {\frac{\partial \,U(r,z,t)}{\partial \,t}} \right|_{t = 0} = \frac{\partial \,w(r,z)}{\partial \,I}$$
(31)
$$t = 0:\;\;\;\;\;U_{b} (r,z,t) = 0.$$
(32)

In Fig. 2, the summary of the identification algorithm is presented.

Fig. 2
figure 2

Flowchart of the identification algorithm

Full size image

4 Method of Solution

The direct problem and additional ones connected with the sensitivity functions are solved using the finite difference method. For this purpose, the uniform spatial grid of dimensions n × n (h is the grid spacing) shown in Fig. 3 is introduced.

Fig. 3
figure 3

Differential grid

Full size image

At first, the optical diffusion Eq. (8) is approximated

$$D\left( {\frac{1}{{r_{i,j} }}\frac{{\phi_{{d\,{\kern 1pt} i,j + 1}} - \phi_{{d\,{\kern 1pt} i,j - 1}} }}{2h} + \frac{{\phi_{{d\,{\kern 1pt} i,j + 1}} - 2\phi_{{d\,{\kern 1pt} i,j}} + \phi_{{d\,{\kern 1pt} i,j - 1}} }}{{h^{2} }} + \frac{{\phi_{{d\,{\kern 1pt} i - 1,j}} - 2\phi_{{d\,{\kern 1pt} i,j}} + \phi_{{d\,{\kern 1pt} i + 1,j}} }}{{h^{2} }}} \right) - \mu_{a\,} \phi_{d\,i,j} + \mu_{s}^{{\prime }} \phi_{c\,i,j} = 0$$
(33)

and then (i = 1, 2,…, n − 1, j = 1, 2,…, n − 1)

$$\begin{aligned} \begin{array}{*{20}c} {\begin{array}{*{20}c} {\phi_{d\,i,j} = \frac{{D\left( {2r_{i,j} - h} \right)}}{{2r_{i,j} \left( {\mu_{a\,} h^{2} + 4D} \right)}}\phi_{{d\,{\kern 1pt} i,j - 1}} + \frac{{D\left( {2r_{i,j} + h} \right)}}{{2r_{i,j} \left( {\mu_{a\,} h^{2} + 4D} \right)}}\phi_{{d\,{\kern 1pt} i,j + 1}} } \\ \end{array} } \\ { + \frac{D}{{\mu_{a\,} h^{2} + 4D}}\left( {\phi_{{d\,{\kern 1pt} i - 1,j}} + \phi_{{d\,{\kern 1pt} i + 1,j}} } \right) + \frac{{\mu^{\prime}_{s} {\kern 1pt} h^{2} }}{{\mu_{a\,} h^{2} + 4D}}\,\phi_{c\,i,j} } \\ \end{array} . \hfill \\ \hfill \\ \end{aligned}$$
(34)

The boundary conditions (10) are also approximated

$$\begin{array}{*{20}c} { - D\frac{{\phi_{d\,i,n} - \phi_{d\,i,n - 1} }}{h} = \frac{{\phi_{d\,i,n} }}{2}} & \to & {\phi_{d\,i,n} = \frac{2D}{2D + h}\phi_{d\,i,n - 1} ,} & {i = 1,2, \ldots n - 1} \\ {D\frac{{\phi_{d\,1,j} - \phi_{d\,0,j} }}{h} = \frac{{\phi_{d\,0,j} }}{2}} & \to & {\phi_{d\,0,j} = \frac{2D}{2D + h}\phi_{d\,1,j} ,} & {j = 1,2, \ldots ,n - 1} \\ { - D\frac{{\phi_{d\,n,j} - \phi_{d\,n - 1,j} }}{h} = \frac{{\phi_{d\,n,j} }}{2}} & \to & {\phi_{d\,n,j} = \frac{2D}{2D + h}\phi_{d\,n - 1,j} ,} & {j = 1,2, \ldots ,n - 1} \\ {D\frac{{\phi_{d\,i,1} - \phi_{d\,i,n} }}{h} = 0} & \to & {\phi_{d\,i,0} = \phi_{d\,i,1} ,} & {i = 1,2, \ldots ,n - 1} \\ \end{array} .$$
(35)

The system of Eqs. (34) and (35) is solved using the iterative method. In a similar way, the problem described by Eqs. (25) and (26) is solved.

Next, Eqs. (1) and (2) are approximated

$$\begin{aligned} C\left[ {\frac{{T_{i,j}^{f} - T_{i,j}^{f - 1} }}{\Delta t} + \frac{{T_{i,j}^{f} - 2T_{i,j}^{f - 1} + T_{i,j}^{f - 2} }}{{\left( {\Delta t} \right)^{2} }}} \right] & = \varLambda \nabla^{2} T_{i,j}^{f - 1} + \frac{{\varLambda \tau_{T} }}{\Delta t}\left( {\nabla^{2} T_{i,j}^{f - 1} - \nabla^{2} T_{i,j}^{f - 2} } \right) \\ & \quad + \;G\left( {T_{b\,i,j}^{f - 1} - T_{i,j}^{f - 1} } \right) + \varepsilon Q_{{m{\kern 1pt} {\kern 1pt} b}} + \left( {1 - \varepsilon } \right)Q_{{m{\kern 1pt} {\kern 1pt} }} + Q_{ex\,\,i,j}^{f - 1} \\ \end{aligned}$$
(36)

and

$$T_{b\,i,j}^{f} = T_{i,j}^{f - 1} - \frac{{\varepsilon \rho_{b} c_{b} }}{G}\frac{{T_{b\,i,j}^{f} - T_{b\,i,j}^{f - 1} }}{\Delta t}$$
(37)

where

$$\nabla^{2} T_{i,j}^{s} = \frac{{T_{i - 1,j}^{s} - 2T_{i,j}^{s} + T_{i + 1,j}^{s} }}{{h^{2} }} + \frac{1}{{r_{i,j} }}\frac{{T_{i,j + 1}^{s} - T_{i,j - 1}^{s} }}{2h} + \frac{{T_{i,j - 1}^{s} - 2T_{i,j}^{s} + T_{i,j + 1}^{s} }}{{h^{2} }}$$
(38)

while s = f − 1 or s = f − 2 and ∆t is the time step.

After the mathematical manipulations one has

$$\begin{aligned} & T_{i,j}^{f} = \frac{{C\left( {\Delta t + 2\tau_{q} } \right) - G\left( {\Delta t} \right)^{2} }}{{C\left( {\Delta t + \tau_{q} } \right)}}T_{i,j}^{f - 1} - \frac{{\tau_{q} }}{{\Delta t + \tau_{q} }}T_{i,j}^{f - 2} + \frac{{\varLambda {\kern 1pt} \Delta t\left( {\Delta t + \tau_{T} } \right)}}{{C\left( {\Delta t + \tau_{q} } \right)}}\left[ {\left( {\frac{1}{{h^{2} }} - \frac{1}{{2h{\kern 1pt} r_{i,j} }}} \right)T_{i,j - 1}^{f - 1} + \left( {\frac{1}{{h^{2} }} + \frac{1}{{2h{\kern 1pt} r_{i,j} }}} \right)T_{i,j + 1}^{f - 1} + \frac{1}{{h^{2} }}\left( {T_{i - 1,j}^{f - 1} + T_{i + 1,j}^{f - 1} } \right) - \frac{4}{{h^{2} }}T_{i,j}^{f - 1} } \right] \\ & - \frac{{\varLambda {\kern 1pt} \Delta t\tau_{T} }}{{C\left( {\Delta t + \tau_{q} } \right)}}\left[ {\left( {\frac{1}{{h^{2} }} - \frac{1}{{2h{\kern 1pt} r_{i,j} }}} \right)T_{i,j - 1}^{f - 2} + \left( {\frac{1}{{h^{2} }} + \frac{1}{{2h{\kern 1pt} r_{i,j} }}} \right)T_{i,j + 1}^{f - 2} + \frac{1}{{h^{2} }}\left( {T_{i - 1,j}^{f - 2} + T_{i + 1,j}^{f - 2} } \right) - \frac{4}{{h^{2} }}T_{i,j}^{f - 2} } \right] + \frac{{\left( {\Delta t} \right)^{2} }}{{C\left( {\Delta t + \tau_{q} } \right)}}\left[ {GT_{b\;i,j}^{f} + \varepsilon Q_{{m{\kern 1pt} {\kern 1pt} b}} + \left( {1 - \varepsilon } \right)Q_{{m{\kern 1pt} {\kern 1pt} }} + Q_{ex\,\,i,j}^{f - 1} } \right] \\ \end{aligned}$$
(39)

and

$$T_{b\,i,j}^{f} = \frac{G\Delta t}{{G\Delta t + \varepsilon \rho_{b} c_{b} }}\left( {T_{i,j}^{f - 1} + \frac{{\varepsilon \rho_{b} c_{b} }}{G\Delta t}T_{b\,i,j}^{f - 1} } \right)$$
(40)

The approximate form of boundary condition (3) is as follows (i = 1, 2,…, n − 1, j = 1, 2,…, n − 1):

$$\left[ {{\mathbf{n}} \cdot \nabla T\left( {r,\,\,z,\,\,t} \right)} \right]_{i,j}^{f} + \frac{{\tau_{T} }}{\Delta t}\left\{ {\left[ {{\mathbf{n}} \cdot \nabla T\left( {r,\,\,z,\,\,t} \right)} \right]_{i,j}^{f} - \left[ {{\mathbf{n}} \cdot \nabla T\left( {r,\,\,z,\,\,t} \right)} \right]_{i,j}^{f - 1} } \right\} = 0.$$
(41)

It means

$$\begin{array}{*{20}c} {\frac{{\Delta t + \tau_{T} }}{\Delta t}\frac{{T_{i,n}^{f} - T_{i,n - 1}^{f} }}{h} - \frac{{\tau_{T} }}{\Delta t}\frac{{T_{i,n}^{f - 1} - T_{i,n - 1}^{f - 1} }}{h} = 0} & \to & {T_{i,n}^{f} = T_{i,n - 1}^{f} + \frac{{\tau_{T} }}{{\Delta t + \tau_{T} }}\left( {T_{i,n}^{f - 1} - T_{i,n - 1}^{f - 1} } \right)} \\ { - \frac{{\Delta t + \tau_{T} }}{\Delta t}\frac{{T_{i,1}^{f} - T_{i,0}^{f} }}{h} + \frac{{\tau_{T} }}{\Delta t}\frac{{T_{i,1}^{f - 1} - T_{i,1}^{f - 1} }}{h} = 0} & \to & {T_{i,0}^{f} = T_{i,1}^{f} - \frac{{\tau_{T} }}{{\Delta t + \tau_{T} }}\left( {T_{i,0}^{f - 1} - T_{i,1}^{f - 1} } \right)} \\ { - \frac{{\Delta t + \tau_{T} }}{\Delta t}\frac{{T_{1,j}^{f} - T_{0,j}^{f} }}{h} + \frac{{\tau_{T} }}{\Delta t}\frac{{T_{1,j}^{f - 1} - T_{0,j}^{f - 1} }}{h} = 0} & \to & {T_{0,j}^{f} = T_{1,j}^{f} - \frac{{\tau_{T} }}{{\Delta t + \tau_{T} }}\left( {T_{1,j}^{f - 1} - T_{0,j}^{f - 1} } \right)} \\ {\frac{{\Delta t + \tau_{T} }}{\Delta t}\frac{{T_{n,j}^{f} - T_{n - 1,j}^{f} }}{h} - \frac{{\tau_{T} }}{\Delta t}\frac{{T_{n,j}^{f - 1} - T_{n - 1,j}^{f - 1} }}{h} = 0} & \to & {T_{n,j}^{f} = T_{n - 1,j}^{f} + \frac{{\tau_{T} }}{{\Delta t + \tau_{T} }}\left( {T_{n,j}^{f - 1} - T_{n - 1,j}^{f - 1} } \right)} \\ \end{array} .$$
(42)

Because the explicit scheme of finite difference method is considered, the stability criterion should be determined (Majchrzak and Mochnacki 2016).

In a similar way, Eqs. (27) and (28) related to the sensitivity functions, supplemented by boundary and initial conditions (30), (31), (32) are solved.

Finally, the Arrhenius integral (13) and the sensitivity function (16) are calculated

$$A_{i}^{f} = P\,\sum\limits_{k = 1}^{f} {\exp \left( { - \frac{E}{{R\,T_{i}^{k} }}} \right)\Delta t}$$
(43)

and

$$V_{i}^{f} = P\,\sum\limits_{k = 1}^{f} {\frac{E}{{R\,\left( {T_{i}^{k} } \right)^{2} }}\exp \left( { - \frac{E}{{R\,T_{i}^{k} }}} \right)U_{i}^{k} \Delta t}$$
(44)

5 Results of Computations

The cylinder of dimensions 0.02 m × 0.02 m is considered. The following values of parameters are assumed: ρ = 1000 kg/m3, ρb = 1060 kg/m3, c = 4000 J/(kg K), cb = 3770 J/(kg K), λ = λb = 0.5 W/(m K), ε = 0.0041, G = 34,785 W/(m3 K), τq = 0.46772 s, τT = 0.46771 s (Zhang 2009), Qm = Qmb = 250 W/m3, μa = 50 1/m, μs = 16,900 1/m, g = 0.952.

At first, the direct problem is solved under the assumption that the radius of laser beam is equal to rD = 0.001 m and laser intensity equals I = 0.3 MW/m2. The initial tissue and blood temperatures are equal to Tp = 37 °C and initial heating rate w(r, z) = 0. In the Arrhenius integral (13): P = 3.1 × 1098 1/s, E = 6.28 × 105 J/mol, R = 8.314 J/(mol K). The problem is solved using explicit scheme of the finite difference method (100 × 100 nodes, time step ∆t = 0.02 s).

In Figs. 4 and 5 the distributions of collimated fluence rate ϕc(r, z) and diffuse fluence rate ϕd(r, z) are shown, while Fig. 6 illustrates the distribution of source function Qex(r, z) related to the laser irradiation. As can be seen, only in the small part of the tissue (0 ≤ r ≤ 0.004 m, 0 ≤ z ≤ 0.006 m) the source function associated with the laser irradiation is non-zero.

Fig. 4
figure 4

Distribution of function ϕc(r, z) [kW/m2]

Full size image
Fig. 5
figure 5

Distribution of function ϕd(r, z) [kW/m2]

Full size image
Fig. 6
figure 6

Distribution of source function Qex(r, z) [kW/m3]

Full size image

Figures 7 and 8 illustrate the courses of the same functions, as previously in the z direction for r = 0. Maximal values of these functions are clearly visible.

Fig. 7
figure 7

Distribution of function ϕc(r, z) and ϕd(r, z) for r = 0

Full size image
Fig. 8
figure 8

Distribution of function Qex(r, z) for r = 0

Full size image

The tissue temperature distributions after the time 60 and 120 s are shown in Figs. 9 and 10, while in Figs. 11 and 12 the distributions of Arrhenius integral for the same moments of time are presented. The subdomain marked in black (Figs. 11, 12) corresponds to this part of the tissue which is destroyed (A(r, z, tf) ≥ 1).

Fig. 9
figure 9

Tissue temperature distribution after 60 s

Full size image
Fig. 10
figure 10

Tissue temperature distribution after 120 s

Full size image
Fig. 11
figure 11

Arrhenius integral distribution after 60 s

Full size image
Fig. 12
figure 12

Arrhenius integral distribution after 120 s

Full size image

It should be noted that to better illustrate the results of computations, in the Figs. 4, 5, 6, 11 and 12 only the fragment of the domain corresponding to 0 ≤ r ≤ 0.006 m, 0 ≤ z ≤ 0.006 m is shown.

Next, the inverse problem is solved. Thus, on the basis of the values of Arrhenius integrals, at the nodes located in the subdomain 0 ≤ r ≤ 0.001 m and 0 ≤ z ≤ 0.001 m, obtained from the direct problem solution the laser intensity is identified. In Fig. 13 the results for two initial values I0 = 0.25 MW/m2 and I0 = 0.5 MW/m2 are shown. For these values, the iteration process is convergent and the final value of I is obtained after 20 and 25 iterations, respectively.

Fig. 13
figure 13

Identification of laser intensity I (for different initial values of I)

Full size image

In conclusion, it is possible to estimate the intensity of the laser which will ensure the appropriate distribution of Arrhenius integral—it means the assumed tissue destruction.

6 Conclusions

The solution of the inverse problem consisting in estimating the laser intensity which ensures the destruction of assumed tissue subdomain is presented. The temperature field in the biological tissue is described by the generalized dual-phase lag model, while the degree of tissue destruction is determined using the Arrhenius integral. The direct problem and additional one related to the sensitivity analysis are solved using the finite difference method, while to solve the inverse problem the gradient method is applied. However, it should be noted that the gradient algorithm is not always convergent. Thus, appropriate selection of the starting point I0, which ensures convergence, is, therefore, extremely important.

Presented approach can be used, among others, in the case of electromagnetic field interaction on biological tissue, e.g. (Paruch 2017). For the postulated degree of tissue destruction it is possible to identify, for example, the voltage on the electrodes.