Ecuación hiperbólica de transmisión del calor para el estudio de la ablación corneal

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Ecuaci´on hiperb´olica de transmisi´on del calor para el estudio de la ablaci´on corneal E. J. Berjano1, J. A. L´ opez Molina2, M. J. Rivera2, M. Trujillo2 1Dpto. de Ingenier´ıa Electr´onica, Centro de Investigaci´on e Innovaci´on en Bioingenier´ıa, Universidad Polit´ecnica de Valencia, Camino de Vera s/n, 46022, Valencia. E-mail: [email protected]. 2Dpto. Matem´atica Aplicada, Instituto de Matem´atica Pura y Aplicada, Universidad Polit´ecnica de Valencia, Camino de Vera s/n, 46022, Valencia. E-mails: [email protected], [email protected], [email protected] . Palabras clave: ecuaci´on hiperb´olica del calor, c´ornea, ablaci´on Resumen En este trabajo presentamos el estudio del problema de transmisi´on del calor que se produce en el proceso de ablaci´on de la c´ornea humana por radiofrecuencia, tal es el caso de la queratoplastia conductiva para la correcci´on de la hipermetrop´ıa y/o el astigmatismo mediante el empleo de electrodos intracorneales de muy peque˜nas dimensiones. En estos casos, las distancias de inter´es son lo suficientemente peque˜nas y los tiempos de excitaci´on lo suficientemente cortos, como para que el estudio de modelado de la transmisi´on de calor deba contemplar la ecuaci´on hiperb´olica de transmisi´on del calor en lugar de la parab´olica o cl´asica. 1. Introducci´on Es cada vez mayor el n´umero de trabajos en los que para abordar el estudio de un problema de transmisi´on del calor se utiliza el modelo hiperb´olico, en lugar del cl´asico o parab´olico. El modelo cl´asico de Fourier supone una velocidad infinita de transmisi´on del calor y la existencia de flujos de valor infinito. Aunque estos inconvenientes f´ısicos no han supuesto una barrera para utilizar el modelo en multitud de aplicaciones ordinarias a la Ingenier´ıa, s´ı que se encuentran disparidades entre los resultados te´oricos obtenidos con el modelo y la experiencia en situaciones f´ısicas en las que se aplican grandes cantidades de calor en peque˜nos intervalos de tiempo. Con el creciente desarrollo de nuevas tecnolog´ıas es cada vez mayor el n´umero de situaciones f´ısicas de estas caracter´ısticas (como por ejemplo el procesado de materiales 1 E. J. Berjano, J. A. L´opez Molina, M. J. Rivera, M. Trujillo mediante pulsos de l´aser o la irradiaci´on electromagn´etica de s´olidos). En estos casos el empleo del modelo hiperb´olico subsana los errores del cl´asico. El modelo hiperb´olico supone una velocidad finita de transmisi´on del calor y la existencia de flujos de valor finito. La ecuaci´on hiperb´olica del calor en el caso de materiales is´otropos es (ver [3]) ∂T ∂t (x, t) + τ∂2T ∂t2(x, t) = α∆T(x, t) + 1 ρc µQ(x, t) + τ∂Q ∂t (x, t)¶,(1) donde Tes la temperatura, xytla variables espacial y temporal, respectivamente, τel par´ametro de relajaci´on, ρla densidad, cel calor espec´ıfico, αla difusividad y Qlas fuentes internas de calor. Estamos interesados en el estudio del problema de transmisi´on del calor que se produce en el proceso de ablaci´on de la c´ornea humana mediante radiofrecuencia, tal es el caso de la queratoplastia conductiva para la correcci´on de la hipermetrop´ıa y/o el astigmatismo mediante el empleo de electrodos intracorneales de dimensiones muy peque˜nas. En estos casos, las distancias de inter´es son lo suficientemente peque˜nas, y los tiempos de excitaci´on lo suficientemente cortos, como para que el estudio de modelado de la transmisi´on de calor deba contemplar la ecuaci´on hiperb´olica. As´ı, el objetivo de esta comunicaci´on es obtener la soluci´on anal´ıtica del problema de ablaci´on en la cornea bajo el punto de vista del modelo hiperb´olico y compararla con la obtenida con el modelo parab´olico. El inter´es de este trabajo se basa en que los modelos matem´aticos propuestos hasta la fecha para la predicci´on de la temperatura corneal durante el calentamiento con radiofrecuencia no contemplan la transferencia de calor con velocidad finita. El modelo que presentamos permitir´a una mejor estimaci´on de la temperatura, y por lo tanto, una mejora en la comprensi´on de los fen´omenos biof´ısicos implicados en estas terapias. 2. Planteamiento del problema Con objeto de poder abordar el problema de forma anal´ıtica hemos considerado un modelo que supone una simplificaci´on del escenario real. En concreto hemos considerado un electrodo activo esf´erico de radio 45 micras, incrustado en el tejido corneal que se supone de extensi´on infinita. En un punto en el infinito se situar´ıa el electrodo dispersivo. Es natural en el problema que vamos a estudiar el empleo de coordenadas esf´ericas en la ecuaci´on (1), siendo todas las variables independientes de los ´angulos esf´ericos. La fuente de calor viene dada por la expresi´on Q(r, t) = Pr0 4πr4H(t), donde H(t) es la funci´on de Heaviside, Pla energ´ıa aplicada y r0el radio del electrodo. As´ı, la ecuaci´on de gobierno del problema es −αµ∂2T ∂r2(r, t) + 2 r ∂T ∂r (r, t)¶+∂T ∂t (r, t) + τ∂2T ∂t2(r, t) = P α r0 4π k r4³H(t) + τδ(t)´,(2) siendo kla conductividad t´ermica ³α=k ρc ´yδ(t) la distribuci´on de Dirac. 2 Ecuaci´on hiperb´olica del calor en la ablaci´on corneal Las condiciones de iniciales y de contorno son T(r, 0) = T0,∂T ∂t (r, 0) = 0 ∀r > r0 l´ım r→∞ T(r, t) = T0,τ ρ0c r0 3kµ1 τ ∂T ∂t (r0, t) + ∂2T ∂t2(r0, t)¶=∂T ∂r (r0, t)∀t > 0.(3) Para obtener la ´ultima condici´on de (3), junto con la expresi´on del flujo hiperb´olico (ver [3]), se ha hecho una simplificaci´on adicional consistente en suponer que la temperatura del electrodo es constante en todos sus puntos, con lo que el calor que entra en ´el por unidad de tiempo se puede calcular mediante la f´ormula elemental ρ0c0 4πr3 0 3 ∂T ∂t (r0, t), siendo c0yρ0el calor espec´ıfico y la densidad del electrodo. Adimensionalizamos el problema utilizando las siguientes variables ρ:= r r0 ;ξ=α t r2 0 ;λ=α τ r2 0 ;V(ρ, ξ) = 4π k r0 PµTµr0ρ, r2 0ξ α¶−T0¶; donde T0es la temperatura ambiente. La formulaci´on del problema adimensional es −µ∂2V ∂ρ2+2 ρ ∂V ∂ρ ¶+∂V ∂ξ +λ∂2V ∂ξ2=1 ρ4³H(ξ) + λ δ(ξ)´(4) V(ρ, 0) = 0 ,∂V ∂ξ (ρ, 0) = 0 ∀ρ > 1 (5) l´ım ρ→∞ V(ρ, ξ) = 0 ,∂V ∂ξ (1, ξ) + λ∂2V ∂ξ2(1, ξ) = 3 m ∂V ∂ρ (1, ξ)∀ξ > 0,(6) siendo m=ρ0c0 ρc .Tomando transformadas de Laplace L(ρ, s) := Lξ[V(ρ, ξ](ρ, s) respecto aξobtenemos −µ∂2L ∂ρ2+2 ρ ∂L ∂ρ ¶+ (s+λ s2)L=1 ρ4µ1 s+λ¶(7) l´ım ρ→∞ L(ρ, s) = 0 ,(s+λs2)L(1, s) = 3 m ∂L ∂ρ (1, s).(8) Utilizando la nueva funci´on z(ρ, s) := ρL(ρ, s) y el m´etodo de variaci´on de las constantes llegamos a la soluci´on general de (7) L=Ã−1 2pA(s)µ1 s+λ¶Zρ 1 e−√A(s)u u3du +M1(s)!e√A(s)ρ ρ +Ã1 2pA(s)µ1 s+λ¶Zρ 1 e√A(s)u u3du +M2(s)!e−√A(s)ρ ρ,(9) 3 E. J. Berjano, J. A. L´opez Molina, M. J. Rivera, M. Trujillo donde A(s) = s+λ s2yM1(s) y M2(s) son funciones que verifican las condiciones de contorno (8) M1(s) = 1 2pA(s)µ1 s+λ¶Z∞ 1 e√A(s)u u3du M2(s) = −e2√A(s) 2pA(s)µ1 s+λ¶mA(s)−3pA(s)+3 mA(s)+3pA(s)+3 Z∞ 1 e−√A(s)u u3du. Introduciendo el valor de estas funciones en (9), la temperatura adimensional ser´a la transformada inversa de Laplace V(ρ, ξ) : = L−1 s"1 2ρµ1 s+λ¶"e√A(s)ρ pA(s)Z∞ ρ e−√A(s)u u3du +e−√A(s)ρ pA(s) ×ÃZρ 1 e√A(s)u u3du −e2√A(s)mA(s)−3pA(s) + 3 mA(s)+3pA(s) + 3 Z∞ 1 e−√A(s)u u3du!# =L−1[F1] + L−1[F2]−L−1[F3].(10) As´ı, tenemos que calcular la transformada inversa de F1,F2yF3. En primer lugar L−1[F1] = L−1"e√A(s)ρ 2ρpA(s)µ1 s+λ¶Z∞ ρ e−√A(s)u u3du# =Z∞ ρÃL−1"1 s e√A(s)(ρ−u) pA(s)#+λL−1"e√A(s)(ρ−u) pA(s)#! du 2ρu3= y por el teorema de traslaci´on de las transformadas de Laplace y la f´ormula (36) de la secci´on 5.6 del cap´ıtulo 5 de [2] =Z∞ ρ H¡ξ−√λ(u−ρ)¢µλe−ξ 2λI0µ1 2λpξ2−λ(u−ρ)2¶ +Zξ √λ(u−ρ) e−v 2λI0µ1 2λpv2−λ(u−ρ)2dv¶!du 2√λ ρu3= =Zρ+ξ √λ ρµλe−ξ 2λI0µ1 2λpξ2−λ(u−ρ)2¶ +Zξ √λ(u−ρ) e−v 2λI0µ1 2λpv2−λ(u−ρ)2dv¶!du 2√λ ρu3.(11) Procediendo de forma an´aloga al c´alculo de la inversa de F1obtenemos L−1[F2] = Zρ 1 H³ξ−√λ(ρ−u)´µλ e−ξ 2λI0µ1 2λpξ2−λ(ρ−u)2¶ +Zξ √λ(ρ−u) e−v 2λI0µ1 2λpv2−λ(ρ−u)2¶dv!du 2√λ ρ u3= 4 Ecuaci´on hiperb´olica del calor en la ablaci´on corneal =Hµρ−ξ √λ−1¶Zρ ρ−ξ √λµλ e−ξ 2λI0µ1 2λpξ2−λ(ρ−u)2¶ +Zξ √λ(ρ−u) e−v 2λI0µ1 2λpv2−λ(ρ−u)2¶dv!du 2√λ ρ u3 +Hµξ √λ−ρ+ 1¶Zρ 1µλ e−ξ 2λI0µ1 2λpξ2−λ(ρ−u)2¶ +Zξ √λ(ρ−u) e−v 2λI0µ1 2λpv2−λ(ρ−u)2¶dv!du 2√λ ρ u3(12) donde I0(z) denota la funci´on modificada de Bessel de orden 0. El c´alculo de L−1[F3] resulta ser el m´as complicado y sutil del art´ıculo. Para ello se utiliza el teorema de convoluci´on en la forma L−1[F3] = Z∞ 1 L−1"µ1 s+λ¶e−√A(s)(u+ρ−2) pA(s) mA(s)−3pA(s)+3 mA(s) + 3pA(s)+3#du 2ρ u3 =Z∞ 1ÃL−1"e−√A(s)(u+ρ−2) pA(s)#∗L−1"µ1 s+λ¶Ã1−6pA(s) mA(s)+3pA(s)+3!#! du 2ρu3. La primera inversa que aparece en la expresi´on anterior es G(ρ, ξ, u) : = L−1"e−√A(s)(u+ρ−2) pA(s)# =H³ξ−√λ¡u+ρ−2¢´e−ξ 2λ √λI0µ1 2λpξ2−λ(u+ρ−2)2¶ =H³ξ−√λ¡u+ρ−2¢´G1(ρ, ξ, u).(13) Para calcular la segunda inversa L−1"µ1 s+λ¶Ã1−6pA(s) mA(s) + 3pA(s)+3!# =¡1 + λδ(ξ)¢−L−1"µ1 s+λ¶6pA(s) mA(s) + 3pA(s)+3# utilizamos la f´ormula de Bromwich, mediante el circuito habitual en el caso de tener puntos de ramificaci´on. En nuestro caso tales puntos son s= 0 y s=−1/λ y adem´as existen dos polos simples s1=1 2λm µ−m−qm2+ 2λ¡9−6m+√81 −108m¢¶ s2=1 2λm µ−m−qm2+ 2λ¡9−6m−√81 −108m¢¶. 5 E. J. Berjano, J. A. L´opez Molina, M. J. Rivera, M. Trujillo interiores al contorno de Bromwich. Obs´ervese que s1ys2son n´umeros reales si 0 ≤m < 3 4. Como las situaciones pr´acticas m´as frecuentes est´an cubiertas con este caso, a partir de ahora supondremos que el valor de mest´a en este intervalo. Llamando r1=1 √2mq9−6m+√81 −108m;r2=1 √2mq9−6m−√81 −108m y utilizando el teorema de los residuos se obtiene R(ξ) : = L−1"µ1 s+λ¶6pA(s) mA(s)+3pA(s) + 3#=−12 2 X k=1 eskξ(1 + λ sk)r2 k sk(1 + 2λ sk)(2rkm−3) −6 πZ0 −1 λ es ξ 1 + λ s s¡m(s+λs2)+3¢√−s−λs2 (m(s+λs2) + 3)2−(s+λs2)ds . De (13) se tiene que G(ρ, ξ, u)∗¡1 + λ δ(ξ)−R(ξ)¢=Zξ 0 G(ρ, ξ −v, u)(1 + λδ(v)−R(v)¢dv =λG(ρ, ξ, u) + Zξ 0 G(ρ, ξ −v, u)¡1−R(v)¢dv y por tanto L−1[F3] = Z∞ 1µλG(ρ, ξ, u) + Zξ 0 G(ρ, ξ −v, u)¡1−R(v)¢dv¶du 2ρu3 =Hµξ √λ−ρ+ 1¶Zξ √λ−ρ+2 1 (λ G1(ρ, ξ, u) +Zξ−√λ(u+ρ−2) 0 G1(ρ, ξ, u)¡1−R(v)¢dv!du 2ρu3.(14) Concluyendo, la temperatura adimensional que quer´ıamos calcular V(ρ, ξ, λ, m) se obtiene a partir de (10), (11), (12) y (14). 3. Comparaci´on con el modelo parab´olico Para comparar los resultados anteriores con los obtenidos con el modelo cl´asico de Fourier calculamos la temperatura parab´olica TF(r, t, m) para el mismo problema. Utilizando las mismas variables adimensionales ρyξ, la temperatura adimensional de Fourier VF(ρ, ξ, m) es la soluci´on del problema −µ∂2V ∂ρ2+2 ρ ∂V ∂ρ ¶+∂V ∂ξ =1 ρ4 V(ρ, 0) = 0 ,∂V ∂ξ (ρ, 0) = 0 ∀ρ > 1 l´ım ρ→∞ V(ρ, ξ) = 0 ,∂V ∂ξ (1, ξ) = 3 m ∂V ∂ρ (1, ξ)∀t > 0. 6 Ecuaci´on hiperb´olica del calor en la ablaci´on corneal Utilizando los mismos m´etodos que en el caso hiperb´olico obtenemos (si 0 ≤m < 3 4) VF(ρ, t, m) = Z∞ 1 1 2ρu3 Zξ 0 e−(u−ρ)2 4v √πv dv du −Z∞ 1 1 2ρu3 Zξ 0Ã1−6 πZ∞ 0 (3 −my)e−y(ξ−v) √y(m2y2+ (9 −6m)y+ 9) dy!e−(u+ρ−2)2 4v √πv  dv . 4. Resultados y conclusiones En las simulaciones empleamos la caracter´ısticas el´ectricas y t´ecnicas de la c´ornea de un trabajo previo [1] y el tiempo de relajaci´on t´ermica de 0,1s. Suponemos un electrodo de platino [4]. La potencia aplicada fue de 30 mW durante 600 ms y la temperatura inicial 35◦C. La figura 1 muestra que al comienzo del calentamiento con la ecuaci´on parab´olica se obtienen valores m´as bajos de temperatura que con la hiperb´olica, lo cual es altamente relevante en las aplicaciones quir´urgicas. Desde el punto de vista matem´atico, es significativa la presencia de singularidades de la soluci´on localizadas a lo largo de la recta (en variables adimensionales) ξ=√λ(ρ−1),que refleja la naturaleza hiperb´olica de la ecuaci´on de gobierno (2). 0.00006 0.00008 0.0001 0.00012 0.00014 0.00016 0.00018 r 40 50 60 70 80 90 T t=0.05 t=0.03 t=0.01 t=0.05 t=0.03 t=0.01 Figura 1: Temperaturas de la c´ornea (hiperb´olica trazo continuo, Fourier discontinuo) desde r= 45 10−6mhasta r= 185 10−6men diferentes tiempos. Agradecimientos Este trabajo ha sido financiado parcialmente por el MEC y la FEDER, Proyecto MTM2004-02262 y la red de investigaci´on MTM2006-26627-E, el ”Plan Nacional de Investigaci´on Cient´ıfica, Desarrollo e Innovaci´on Tecnol´ogica del Ministerio de Educaci´on y Ciencia”de Espa˜na (TEC 2005-04199/TCM) y el Programa de Apoyo a la Investigaci´on y Desarrollo (PAID-04-07) de la UPV. 7 E. J. Berjano, J. A. L´opez Molina, M. J. Rivera, M. Trujillo Referencias [1] E. J. Berjano, J. L. Alio & J. Saiz, Modeling for radio-frequency conductive keratoplasty: implications for the maximum temperature reached in the cornea, Physiol. Meas. 26, (2005), 157-172. [2] A. Erd´elyi, M. F. Oberhettinger & F. G. Tricomi, Tables of Integral Transforms. Vol. I and II, McGraw-Hill, New York, 1954. [3] M. N. ¨ Ozi¸sik & D. T. Tzou, On the wave theory in heat conduction, ASME Journal of Heat Transfer, 116, (1994) 526-535. [4] D. Panescu, J. G. Whayne, S. D. Fleishman, M. S. Mirotznik, D. K. Swanson & J. G. Webster Threedimensional finite analysis of current density and temperature distributions during radio-frequency ablation, IEEE Trans. Biomed. Eng., 42, (1995), 879-890. 8