Resolución numérica de algunos sistemas parabólico–elípticos no lineales

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–7) Resoluci´on num´erica de algunos sistemas parab´olico–el´ıpticos no lineales Jos´ e Manuel D´ ıaz Moreno1, Mar´ ıa Teresa Gonz´ alez Montesinos2, Francisco Orteg´ on Gallego1 1Dpto. de Matem´aticas, Universidad de C´adiz, CASEM, Campus del R´ıo San Pedro, 11510 Puerto Real, C´adiz. E-mails: josemanuel.di[email protected], [email protected]. 2Dpto. de Matem´aticas, Universidad de C´adiz, Escuela Superior de Ingenier´ıa de C´adiz, 11002 C´adiz. E-mail: mariateresa.gonzal[email protected]. Palabras clave: Problema del termistor, silistor, ecuaciones el´ıpticas no lineales, ecuaciones parab´olicas no lineales, soluciones de capacidad, soluciones d´ebiles, resoluci´on num´erica de problemas no lineales, m´etodo de los elementos finitos. Resumen En este trabajo se aborda la resoluci´on num´erica del problema del termistor en dimensi´on dos. La mayor dificultad de dicho an´alisis reside en el hecho que se supone conducci´on met´alica y adem´as se satisface la ley de Wiedemann–Franz. Para la simulaci´on num´erica se ha considerado el caso de un termistor PTC de medida o silistor. En tal caso, los par´ametros f´ısicos reales conducen un problema adimensional con difusi´on muy peque˜na, lo que puede dar lugar a la aparici´on de una capa l´ımite a lo largo de la frontera del silistor. 1. Introducci´on El calor producido por una corriente el´ectrica que atraviesa un semiconductor est´a descrito por el llamado problema del termistor, consistente en un sistema de dos ecuaciones acopladas, una parab´olica no lineal y otra el´ıptica, cuyas inc´ognitas son la temperatura, u, y el potencial el´ectrico, ϕ. Las leyes de Ohm y Fourier vienen dadas por J=σ(u)EyQ=−a(u)∇u, respectivamente, donde Jes la intensidad de corriente el´ectrica, Qel flujo de calor, E=−∇ϕel campo el´ectrico, y σ(u) y a(u) son sendas conductividades el´ectrica y t´ermica. El problema 1 J. M. D´ıaz Moreno, M. T. Gonz´alez Montesinos, F. Orteg´on Gallego del termistor se deduce a partir de las leyes de conservaci´on de la corriente y la energ´ıa, a saber, ∇ · J = 0, ρc∂u ∂t +∇ · Q =J · E, siendo ρla densidad del semiconductor, y csu calor espec´ıfico. As´ı, obtenemos el siguiente problema: ρc∂u ∂t − ∇ · (a(u)∇u) = σ(u)|∇ϕ|2en ΩT= Ω ×(0, T), ∇ · (σ(u)∇ϕ) = 0 en ΩT, u= 0 sobre ΓT=∂Ω×(0, T), ϕ=ϕ0sobre ΓT, u(·,0) = u0en Ω,              (1) donde Ω ⊂RN, dominio ocupado por el dispositivo el´ectrico, es un abierto acotado y regular, N≥2 y T > 0. Son escasos los trabajos en los que se aborda la resoluci´on num´erica de (1), pero entre ellos cabe destacar [1] y [13], en los que se estudia el problema del termistor con conductividad t´ermica constante y condiciones de contorno mixtas. Por otro lado, en [11] se analiza un problema similar considerando el llamado modelo de entalp´ıa y, de nuevo, a es una funci´on constante. Resumiendo, en estos art´ıculos y en tantos otros, la funci´on a siempre se supone constante. Ahora bien, cuando la conductividad t´ermica es de tipo Wiedemann–Franz, es decir, a(s) = Lsσ(s), siendo L > 0 la constante de Lorenz, y se produce conducci´on met´alica, esto es, σ(s) = O(s−1) para |s| → +∞, el estudio del problema (1) se complica sobremanera debido al car´acter degenerado de la ecuaci´on parab´olica y al no uniformemente el´ıptico de la ecuaci´on el´ıptica. Actualmente, bajo estas hip´otesis sobre las conductividades, la existencia de soluciones d´ebiles de (1) constituye un problema abierto; no obstante, en [8] los autores han demostrado la existencia de un cierto tipo de soluci´on que se adapta al marco funcional de (1): la soluci´on de capacidad. La finalidad de este trabajo estriba en mostrar algunos de los resultados obtenidos en la resoluci´on num´erica del sistema (1) en el caso bidimensional, suponiendo que la conductividad t´ermica satisface la ley de Wiedemann–Franz y adem´as se produce conducci´on met´alica, lo cual se corresponde con los casos f´ısicamente importantes y m´as complicados desde el punto de vista te´orico. 2. Soluci´on de capacidad Supongamos las siguientes hip´otesis sobre los datos del sistema (1): (H.1) u0∈L2(Ω) es tal que u0≥0 casi por doquier en Ω. (H.2) ϕ0∈L2(0, T;H1(Ω)) ∩L∞(ΩT). (H.3) σ∈C(R) y 0 < σ(s)≤σ0, para cualquier s∈R. (H.4) a∈C(R) y 0 < a(s)≤a0, para cualquier s6= 0, a(0) = 0. (H.5) Para cada δ > 0 existe una constante aδ>0 tal que ´ınf ese|s|>δa(s)≥aδ. 2 Resoluci´on num´erica del problema del termistor Este mismo problema se estudia en [4], donde se demuestra la existencia de soluciones d´ebiles del sistema (1), pero ahora estamos suponiendo que la funci´on σno est´a acotada inferiormente por una constante positiva. Este cambio en la hip´otesis (H.3) sobre la conductividad el´ectrica puede parecer insignificante pero es crucial, ya que el an´alisis del problema es mucho m´as complejo: las ecuaciones parab´olica y el´ıptica de (1) ahora van a ser degenerada y no uniformemente el´ıptica, respectivamente. De este modo, la existencia de soluciones d´ebiles no est´a garantizada y hemos de tratar con otro tipo de soluciones, a saber, las soluciones de capacidad (v´eanse [7, 8, 12]). N´otese adem´as que (H.3)–(H.5) incluyen el caso de conducci´on met´alica para σ, mientras que apuede definirse haciendo uso la ley Wiedemann-Franz. Sea A(s) = Rs 0a(τ) dτ. Claramente A(0) = 0, A∈C1(R), Aes estrictamente creciente y globalmente lipschitziana. Definici´on 1 Se dice que la terna (u, ϕ, Φ) es soluci´on de capacidad de (1) si se cumplen las siguientes condiciones: (C.1) u∈L∞(0, T;L1(Ω)),du dt∈L2(0, T;H−1(Ω)),A(u)∈L2(0, T;H1 0(Ω)) ∩Lq(ΩT), para cualquier q < 2 + 2/N,ϕ∈L∞(ΩT)yΦ∈L2(ΩT)N. (C.2) (u, ϕ, Φ) verifica el sistema de ecuaciones diferenciales ∂u ∂t −∆A(u) = ∇ · (ϕΦ) en L2(0, T;H−1(Ω)), ∇ · Φ = 0 en L2(0, T ;H−1(Ω)). (C.3) Para cada S∈W1,∞(R)tal que sop Ses compacto, se tiene que S(A(u))ϕ−S(0)ϕ0∈L2(0, T;H1 0(Ω)), S(A(u))Φ = σ(u) [∇(S(A(u))ϕ)−ϕ∇S(A(u))] . (C.4) u(·,0) = u0. En [8] se prueba el resultado de existencia que se enuncia a continuaci´on. Teorema 1 Bajo las hip´otesis (H.1)–(H.5), el sistema (1) admite soluci´on de capacidad u en el sentido de la definici´on 1. M´as a´un, u≥0casi por doquier en ΩT, el gradiente de uest´a definido casi por doquier en ΩTy es tal que ∇uχ{u>δ}∈L2(ΩT)para cualquier δ > 0. Finalmente, si S∈L1 loc(R)es tal que S′∈L∞(R)ysop S⊂R\(−δ0, δ0)para alg´un δ0>0, entonces S(u)∈L2(0, T;H1 0(Ω)) y∇S(u) = S′(u)∇uen ΩT. 3. Resoluci´on num´erica del problema Un termistor PTC (positive temperature coefficient) es un dispositivo el´ectrico semiconductor cuya resistencia aumenta con la temperatura. El modelo que presentamos en este trabajo se corresponde con el de un termistor PTC de medida, basado en silicio dopado, tambi´en conocido como silistor. 3 J. M. D´ıaz Moreno, M. T. Gonz´alez Montesinos, F. Orteg´on Gallego Nuestro objetivo es llevar a cabo un an´alisis num´erico del comportamiento de un silistor t´ıpico, cuyo rango de temperatura oscila aproximadamente entre −60oC = 213oK y 300oC = 573oK, y cuya resistencia el´ectrica viene dada R(u) = R0r(u), siendo r(u) = 1 + α(u−u0) + β(u−u0)2,(2) donde u0= 25oC = 298oK, R0=R(u0) y las constantes αyβdependen de las propiedades del termistor. Las conductividades el´ectrica y t´ermica vienen dadas pues por σ(u) = H SR(u), a(u) = Luσ(u) = LHu SR(u),(3) donde Ses el ´area de la secci´on del dispositivo, H×H×H/2 son las dimensiones del mismo, se satisface la ley de Wiedemann-Franz y se produce conducci´on met´alica. Teniendo en cuenta las hip´otesis (H.1)–(H.5), en lugar del problema (1), vamos a considerar este otro: ρc∂u ∂t − ∇ · (a(u)∇u) = σ(u)|∇ϕ|2en ΩT, ∇ · (σ(u)∇ϕ) = 0 en ΩT, a(u)∂u ∂n +h(x, y)(u−u0) = 0 sobre ΓT, ϕ=±ϕ0sobre Γ± S×(0, T), ∂ϕ ∂n = 0 sobre ΓN×(0, T), u(·,0) = u0en Ω,                          (4) donde Ω = (0, H)×(0, H), Γ = ∂Ω = ΓS∪ΓN, con ΓS∩ΓN=∅, Γ± Sson los lados del dispositivo correspondientes a la entrada y salida de corriente el´ectrica, y hes el coeficiente de transferencia de calor, que depende del silistor. Concretamente, estudiamos el comportamiento de un silistor durante una hora de funcionamiento, integrado en un circuito que recibe un voltaje inicial V0y cuya temperatura m´axima es uM. Sus datos f´ısicos reales son los que siguen: ρ= 2.3·103kg m−3, c = 701 J kg−1oK−1, H = 0.01 m, R0= 500 Ω, uM= 600oK, α = 7.874 ·10−3, β = 1.874 ·10−5, V0= 250 V, L= 2.44 ·10−8WΩoK−2, u0= 298oK, ϕ0(t) = 1 2cos π 544 t+ 3, T = 3600 s. (5) En cuanto a h, un valor est´andar suele ser 4 ·102WoK−1m−2; sin embargo, en vista de que los contactos del silistor son met´alicos, podemos tomar h= 102WoK−1m−2en los contactos, ΓS, y h= 10 WoK−1m−2en el resto de la frontera, ΓN. 3.1. Modelo adimensional El paso previo a la resoluci´on num´erica del sistema (4) pasa por considerar un modelo adimensional del mismo. Por ello, se introducen las siguientes variables adimensionales (designadas con una raya) u−u0=uM¯u, ϕ =V0 2¯ϕ, x =H¯x, y =H¯y, t =T¯ t, R(u) = R0¯ R(¯u), σ(u) = σ0¯σ(¯u), a(u) = a0¯a(¯u), (6) 4 Resoluci´on num´erica del problema del termistor donde σ0=σ(u0) y a0=a(u0). Por otro lado, gracias a (2) y (3), ¯ R(¯u) = r(u) = 1 + αuM¯u+βu2 M¯u2,¯σ(¯u) = 1 ¯ R(¯u),¯a(¯u) = u0+uM¯u u0¯ R(¯u).(7) En vista de (6) y (7), y prescindiendo de las barras usadas hasta ahora, obtenemos el problema ∂u ∂t −γ∇ · (a(u)∇u) = ησ(u)|∇ϕ|2en ΩT, ∇ · (σ(u)∇ϕ) = 0 en ΩT, a(u)∂u ∂n +κu = 0 sobre ΓT, ϕ=±ϕ0sobre Γ± S×(0, T), ∂ϕ ∂n = 0 sobre ΓN×(0, T), u(·,0) = 0 en Ω,                          (8) donde Ω = (0,1) ×(0,1), T= 1, γ=Lu0T ρcSHR0,η=V2 0T 4ρcSHLuMR0yκ=SR0h Lu0. Para la resoluci´on num´erica de (8), consideramos un esquema de Euler impl´ıcito en tiempo, dividiendo el intervalo [0,1] en Nsubintervalos de longitud τ, y escribiremos fn(x, y) = f(x, y, nτ), n≥0, para cualquier funci´on fdefinida en ΩT. De este modo, para n≥0, se considera la siguiente sucesi´on de problemas el´ıpticos: dados unyϕn, resolver un+1 −un τ−γ∇ · (a(un)∇un+1) = ησ(un)|∇ϕn|2en Ω, ∇ · (σ(un+1)∇ϕn+1) = 0 en Ω, a(un)∂un+1 ∂n +κun+1 = 0 sobre Γ, ϕn+1 =±ϕn+1 0sobre Γ± S, ∂ϕn+1 ∂n = 0 sobre ΓN,                      (9) donde u0= 0 y ϕ0∈H1(Ω) es la soluci´on de ∇ · (σ(u0))∇ϕ0= 0 en Ω, ϕ0=ϕ0 0sobre ∂Ω.(10) La resoluci´on de los problemas (9)–(10) se realiza mediante el m´etodo de los elementos finitos. 3.2. Resultados num´ericos Teniendo en cuenta (5), el coeficiente γde la ecuaci´on para ues muy peque˜no y pueden surgir capas l´ımites a lo largo de ∂Ω. Por esa raz´on, se construye una malla m´as densa en nodos cerca de la frontera que en el resto del dominio; as´ı se consigue evaluar mejor la funci´on a(u) y resolver adecuadamente el problema para ude (9). Para que la malla se adapte mejor a la soluci´on num´erica, se hace uso de un algoritmo adaptativo ([10]) que actualiza la malla cada 100 ciclos de tiempo (v´ease figura 1). Las figuras 2 y 3 muestran diversos gr´aficos de las temperaturas obtenidas con el algoritmo (9)- (10) durante una hora de funcionamiento del silistor. Se observa la pronunciada pendiente de la capa l´ımite a lo largo de toda la frontera. 5 J. M. D´ıaz Moreno, M. T. Gonz´alez Montesinos, F. Orteg´on Gallego (a) 1 iteraci´on (b) 500 iteraciones (c) 1000 iteraciones Figura 1: Mallas seg´un dististas iteraciones en t. (a) 1 iteraci´on (b) 500 iteraciones (c) 1000 iteraciones Figura 2: Funci´on useg´un distintas iteraciones en t. Para estudiar con m´as detalle el comportamiento de la funci´on ucerca de Ω se han realizado tres cortes en superficie integral (figura 3). 0 0.1 0.2 0.3 0.4 0.5 0 0.02 0.04 0.06 0.08 0.1 u x Corte 40 "cortea40.dat" (a) x∈(0,0.1) e y= 0.5 0 0.1 0.2 0.3 0.4 0.5 0 0.02 0.04 0.06 0.08 0.1 u y Corte 40 "corteb40.dat" (b) x= 0.5 e y∈(0,0.1) 0 0.1 0.2 0.3 0.4 0.5 0 0.02 0.04 0.06 0.08 0.1 u x=y Corte 40 "cortec40.dat" (c) x=y∈(0,0.1) Figura 3: Cortes de la gr´afica de la funci´on upara t= 0.5. Agradecimientos Este trabajo ha sido parcialmente financiado por el Ministerio de Educaci´on y Ciencia, proyecto I+D+I MTM2006-04436, con la participaci´on del FEDER, y la Junta de Andaluc´ıa, III Plan Andaluz de Investigaci´on, grupo FQM-315. 6 Resoluci´on num´erica del problema del termistor Referencias [1] A. C. Fowler, I. Frigaard, S. D. Howinson, Temperature surges in current–limiting circuit devices, SIAM J. Appl. Math., vol. 52, 4 (1992), 998–1011. [2] P. J. Frey, MEDIT, An interactive mesh visualization software, INRIA, 0253, Th`eme 4 (2001). [3] M. T. Gonz´alez Montesinos, Estudio matem´atico de algunos problemas no lineales del electromagnetismo relacionados con el problema del termistor, Tesis, Universidad de C´adiz, Espa˜na, 2002. [4] M. T. Gonz´alez Montesinos, F. Orteg´on Gallego, The evolution thermistor problem with degenerate thermal conductivity, Comm. Pure Appl. Anal., vol. 1, 3 (2002), 313–325. [5] M. T. Gonz´alez Montesinos, F. Orteg´on Gallego, On certain non–uniformly and singular non–uniformly elliptic systems, Nonlinear Anal., 54 (2003), 1193–1204. [6] M. T. Gonz´alez Montesinos and F. Orteg´on Gallego, Algunos resultados sobre el problema del termistor, Bol. Soc. Esp. Mat. Apl., 31 (2005), 109–138. [7] M. T. Gonz´alez Montesinos, F. Orteg´on Gallego, Existence of a capacity solution to a coupled nonlinear parabolic–elliptic system. Comm. Pure Appl. Anal., vol. 6, 1 (2007), 23–42. [8] M. T. Gonz´alez Montesinos, F. Orteg´on Gallego, The thermistor problem with degenerate thermal conductivity and metallic conduction. (Aparecer´a en Discrete Contin. Dyn. Syst.–Suplement, 2007). [9] R. Glowinski. Numerical Methods for Nonlinear Variational Problems, Springer-Verlag, New York, 1984. [10] F. Hecht, O. Pironneau, A. Le Hyaric, K. Ohtsuka. FreeFem++, Version 2.11, 2006. [11] P. Shi, M. Shillor, X. Xu, Existence of a solution to the Stefan problem with Joule’s heating, J. Differential Equations, 105 (1993), 239–263. [12] X. Xu, A strongly degenerate system involving an equation for parabolic type and an equation of elliptic type, Comm. Partial Differential Equations, vol. 18, 1&2 (1993), 199–213. [13] S. Zhou, D. R. Westbrook, Numerical solutions of the thermistor equations, J. Comput. Appl. Math., 79 (1997), 101–118. 7