Estudio de la estabilidad dinámica de pares eje–cojinete en problemas evolutivos de lubricación

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Estudio de la estabilidad din´amica de pares eje–cojinete en problemas evolutivos de lubricaci´on J. Durany1, J. Pereira-P´ erez1, F. Varas1 1Dpto. Matem´atica Aplicada II, Universidad de Vigo, Campus Marcosende, E-36313 Vigo. E-mails: [email protected], [email protected]. Palabras clave: Termohidrodin´amica, lubricaci´on, elementos finitos, vol´umenes finitos, elementos de contorno, estabilidad. Resumen En este trabajo se presenta un modelo matem´atico y una combinaci´on de m´etodos num´ericos para la resoluci´on de problemas termohidrodin´amicos evolutivos en lubricaci´on de pares eje-cojinete. En concreto, un m´etodo de elementos finitos para el problema hidrodin´amico, junto a un algoritmo de dualidad para la parte no lineal de frontera m´ovil (cavitaci´on). Adem´as, un esquema de vol´umenes finitos, de tipo cell– vertex de orden dos, para el c´alculo de la temperatura del fluido, y un m´etodo de elementos de contorno para la t´ermica del cojinete. Una de las aplicaciones importantes, por su inter´es industrial, se refiere al an´alisis de la estabilidad del dispositivo mec´anico. En este sentido, la resoluci´on num´erica que se presenta permite comparar el comportamiento din´amico del par eje–cojinete con la curva de estabilidad neutra que se obtiene al utilizar modelos anal´ıticos simplificados. 1. Modelo termohidrodin´amico de lubricaci´on El modelo termohidrodin´amico que se describe es cl´asico en el an´alisis de la lubricaci´on de pares eje-cojinete (ver [4], entre otros). Primero, se formula la ecuaci´on de Reynolds con un modelo de cavitaci´on de Elrod-Adams para el c´alculo de la presi´on del fluido lubricante. A continuaci´on, se plantea el problema t´ermico en el dispositivo mediante la ecuaci´on de la energ´ıa para la temperatura del fluido y unas ecuaciones de conducci´on t´ermica en el eje y en el cojinete. El acoplamiento de las ecuaciones en derivadas parciales viene dada, por una parte, por una ley de variaci´on de la viscosidad en funci´on de la temperatura en la ecuaci´on del flujo y, por otra parte, por medio de la influencia del campo de velocidades en la ecuaci´on de la energ´ıa. Adem´as, el acoplamiento t´ermico del fluido con el dispositivo viene dado por las condiciones de contorno en las paredes de contacto con eje y cojinete. 1 J. Durany, J. Pereira-P´erez, F. Varas Para simplificar los c´alculos se emplean variables adimensionales y el dominio del fluido se transforma en un paralelep´ıpedo. Las ecuaciones para este cambio de variable son: x=Raθ, y =h¯y, z =L¯z, ;u=ωRa¯u, v =ωC¯v, w =ωRa¯w, t =ωτ , (1) donde los valores Ra, L, C y h representan el radio del eje, la longitud del par, el juego radial y la altura de la pel´ıcula de fluido, respectivamente, y (u, v, w) las componentes de la velocidad. La velocidad angular de rotaci´on del eje est´a denotada por ω, y el tiempo real por la variable τ. Se definen tambi´en los siguientes cambios de variable que afectan a otras magnitudes: h=C¯ h, η =Ra/L, T =T0¯ T, µ =µ0¯µ, p =µ0ωRa C2 ¯p. (2) En estas ecuaciones T0yµ0representan la temperatura y la viscosidad del fluido en la ranura de alimentaci´on, que se utilizar´an de referencia, mientras que pyTrepresentan la presi´on y temperatura de la pel´ıcula de lubricante. Realizando el cambio de variable en el espacio comprendido entre el eje y el cojinete, el dominio adimensional se convierte en un volumen comprendido entre la superficie rectangular Ω1= [0,2π]×[0,1] para (θ, ¯z) y la altura: ¯ h= 1 + εcosθ, (3) con εrepresentando el coeficiente de excentricidad del par eje–cojinete. El modelo matem´atico del problema de lubricaci´on con cavitaci´on se realiza mediante el modelo de Elrod-Adams (ver [3], por ejemplo), introduciendo una inc´ognita adicional que representa la concentraci´on o ´ındice de cavitaci´on: ϑ. Esta funci´on ϑvale uno en la zona activa (Ω+ 1) y toma valores en el intervalo (0,1) en la zona cavitada (Ω0 1). De este modo, el problema se formula en los siguientes t´erminos: Hallar (p, ϑ) tales que se verifiquen las condiciones: ∂ ∂θ ¯ h3¯ G∂¯p ∂θ +∂ ∂˜z ¯ h3¯ G∂¯p ∂z =∂ ∂θ ¯ h1− ¯ I2 ¯ J2+∂(ϑ¯ h) ∂t ,(4) ¯p > 0y ϑ = 1 en Ω+ 1×[0, tF] ∂ ∂θ ϑ¯ h1−¯ I2 ¯ J2+∂(ϑ¯ h) ∂t = 0,¯p= 0,0≤ϑ≤1en Ω0 1×[0, tF],(5) ¯ h3¯ G∂¯p ∂~n = (1 −ϑ)¯ h1− ¯ I2 ¯ J2cos(~n,~ i),¯p= 0 en Σ×[0, tF],(6) ¯p= 0 en ∂Ω1×[0, tF], ϑ =ϑ0en Γ0×[0, tF],(7) ¯p(0, θ, ˜z) = ¯p0(θ, ˜z) , ϑ(0, θ,˜z) = ϑ0(θ,˜z) ,(8) donde se ha realizado un nuevo cambio de variable, ˜z = ¯z η, para simplificar las ecuaciones, transformando el dominio en Ω1= [0,2π]×[0,1 η]. Adem´as, Σ representa la frontera libre entre la zona activa Ω+ 1y la zona cavitada Ω0 1,~n es el vector normal a la frontera Σ, ~ ies 2 Estabilidad din´amica de pares eje–cojinete el vector unitario en la direcci´on θ, Γ0la frontera de alimentaci´on de fluido (θ= 0) y tF el tiempo final. Los coeficientes viscosos que intervienen en el problema est´an dados por: ¯ G=Z1 0 ¯y ¯µ¯y− ¯ I2 ¯ J2d¯y , ¯ I2=Z1 0 ¯y ¯µd¯y , ¯ J2=Z1 0 d¯y ¯µ,(9) y la viscosidad del fluido, ¯µ, se considera variable con la temperatura seg´un la ley: ¯µ=e−βT0(¯ T−1) (10) donde βes el coeficiente termoviscoso y T0la temperatura de referencia. Cuando las condiciones de funcionamiento de un par eje–cojinete son normales, los gradientes axiales de temperatura son despreciables y la ecuaci´on de la energ´ıa para la temperatura del fluido lubricante ¯ Tfse puede plantear en un plano medio en coordenadas de longitud y espesor (θ, ¯y) (ver [4], entre otros): ρcp C2 k ∂¯ Tf ∂t +Pe¯u∂¯ Tf ∂θ +1 ¯ h¯v−¯u¯y∂¯ h ∂θ ∂¯ Tf ∂¯y−1 ¯ h2 ∂2¯ Tf ∂¯y2= =Nd ¯µ ¯ h2∂¯u ∂¯y2 ,en Ωf×[0, tF],(11) donde ρ,cpykdenotan densidad, calor espec´ıfico y conductividad, Pees el n´umero de Peclet, Ndun par´ametro adimensional asociado a la disipaci´on viscosa, Ωf= [0,2π]×[0,1] el dominio adimensional de una seci´on longitudinal y las componentes de la velocidad, ~v = (¯u, ¯v), est´an definidas a partir de la presi´on del fluido, ¯p, por: ¯u=¯ h2∂¯p ∂θ ¯ I− ¯ I2 ¯ J2 ¯ J+¯ J ¯ J2 ; ¯v=−¯ hZ¯y 0∂¯u ∂θ∗dξ , (12) con ¯ I2,¯ J2dadas en (9) y las nuevas integrales, ¯ I,¯ J, definidas por: ¯ I=Z¯y 0 ξ ¯µdξ ;¯ J=Z¯y 0 dξ ¯µ.(13) La ecuaci´on de conservaci´on de la energ´ıa (11) debe reformularse sobre la regi´on cavitada para tener en cuenta la presencia de la mezcla. As´ı, con las adimensionalizaciones correspondientes, el problema t´ermico se formula en Ωfmediante la misma expresi´on (11), pero ahora con ρ,cp,k,Pe yNd dependiendo de la concentraci´on del fluido: ρcp=ρfcfϑ+ρaca(1 −ϑ) ; k=kfϑ+ka(1 −ϑ); Pe =ωC2[ρfcfϑ+ρaca(1 −ϑ)] [kfϑ+ka(1 −ϑ)] ;Nd =µ0ω2Ra2 [kfϑ+ka(1 −ϑ)]T0 ,(14) siendo ρf, cf, kf, ρa, ca, kala densidad, el calor espec´ıfico y la conductividad del fluido y el gas, respectivamente. Esta ecuaci´on (11) se completa con las condiciones Dirichlet en la frontera de entrada θ= 0 y en las fronteras de contacto con el eje (¯y= 1) y el cojinete (¯y= 0). 3 J. Durany, J. Pereira-P´erez, F. Varas La ecuaci´on de la energ´ıa en el cojinete se reduce a un problema evolutivo de difusi´on junto a las condiciones de contorno que la acoplan, por una parte, con la temperatura del fluido y, por otra parte, con las condiciones ambientales: ρccpc ∂Tc ∂t −div(kc∇Tc) = 0 en Ωc×[0, tF],(15) −kc ∂Tc ∂~n =hc(Tc−Ta)en Γe c×[0, tF],(16) i)Tc=T0¯ Tfyii)−kc ∂Tc ∂~n =kT0 C¯ h ∂¯ Tf ∂¯yen Γi c×[0, tF],(17) Tc(¯x, 0) = Tc0(¯x)en Ωc,(18) donde Tces la funci´on temperatura en el cojinete, Ωcel dominio anular del cojinete con frontera interior Γi cy exterior Γe c,kcyklos coeficientes de difusi´on en el cojinete y en el fluido, ρcla densidad, cpcel calor espec´ıfico, ~n un vector, unitario, normal y exterior a las fronteras, Tala temperatura ambiente, hcel coeficiente de convecci´on al exterior y Tc0la temperatura inicial. Finalmente, en el eje se considera un modelo de temperatura uniforme en cada instante de tiempo debido a las altas velocidades de rotaci´on. As´ı, el calentamiento o enfriamiento se corresponde con el balance entre el calor recibido del fluido y el flujo t´ermico con el exterior. Es decir: ZΓe kT0 C¯ h d¯ Tf d¯ydl =he(Te−Ta) + Ae ∂Te ∂t (19) donde Tees la temperatura en el eje, Γela frontera con el fluido, heel coeficiente de convecci´on al aire, Aeun coeficiente de capacidad calor´ıfica y la derivada normal en la frontera se corresponde con la vertical. 2. Resoluci´on num´erica del problema acoplado La resoluci´on num´erica del problema acoplado termohidrodin´amico requiere la soluci´on de cada uno de los sub–modelos que se han descrito en la secci´on anterior. En particular, la soluci´on del problema hidrodin´amico (4)–(8) en el fluido ha sido el objetivo de algunos trabajos previos de los autores. En esta secci´on se presenta un esquema del m´etodo num´erico desarrollado en Durany–Pereira–Varas [2], [3], para problemas estacionarios y transitorios, respectivamente. En concreto, el esquema pasa por la construcci´on de una formulaci´on d´ebil en t´erminos del operador mult´ıvoco de Heaviside, H, para tratar el problema de frontera m´ovil de la cavitaci´on (ver Durany–Garc´ıa-V´azquez [1], para los detalles). A continuaci´on, se realiza una semidiscretizaci´on en tiempo de la derivada material utilizando el m´etodo de las caracter´ısticas, junto a una aproximaci´on espacial con elementos finitos lineales a trozos. Finalmente, se aplica un algoritmo de dualidad para resolver la no linealidad debida al operador de Heaviside. La ecuaci´on de la energ´ıa en el fluido lubricante se resuelve mediante un esquema de vol´umenes finitos de tipo cell-vertex de orden dos, siguiendo el trabajo de Morton-Stynes- S¨uli [5]. Este m´etodo combina unas buenas propiedades de monoton´ıa en la recuperaci´on de capas l´ımite con una convergencia de orden dos bajo ciertas hip´otesis de regularidad 4 Estabilidad din´amica de pares eje–cojinete (ver Durany–Pereira–Varas [2]). Adem´as, conduce a un ensamblado muy sencillo, parecido al del m´etodo de elementos finitos, que se simplifica si las mallas son muy estructuradas. La resoluci´on del problema t´ermico en el cojinete se aborda empleando un m´etodo de reciprocidad dual combinado con una discretizaci´on de la frontera, tanto interna como externa, mediante elementos de contorno P1. Los nodos de la discretizaci´on son los v´ertices de los arcos que, en la frontera interna de contacto con el fluido, coinciden con los nodos de la malla utilizada en el m´etodo de vol´umenes finitos. Se imponen las condiciones de continuidad del flujo en la frontera interior y la igualdad del flujo a la convecci´on al ambiente en la frontera exterior. Este procedimiento tiene la ventaja de no tener que mallar todo el cojinete y, junto con la simetr´ıa del problema, permite reducir en gran medida el coste computacional de esta etapa (ver Durany–Pereira–Varas [3]). Finalmente, el problema t´ermico (19) en el eje puede resolverse expl´ıcitamente de manera muy simple, ya que se trata de un problema unidimensional de balance de energ´ıa. 3. An´alisis de la estabilidad din´amica del eje Se puede determinar la fuerza ejercida por el fluido sobre el eje y relacionar el desplazamiento del centro del eje y la velocidad de ese desplazamiento con la variaci´on de la fuerza (ver Frene et al. [4], para los detalles). Para ello, se presenta en primer lugar un m´etodo cl´asico basado en la obtenci´on de un modelo anal´ıtico de la presi´on y una aproximaci´on lineal en el entorno del punto de equilibrio. Esta t´ecnica simplificada tiene el inconveniente de prefijar la posici´on de la zona cavitada y, en consecuencia, es importante comparar los resultados que se obtienen con los correspondientes a utilizar el modelo hidrodin´amico de Elrod–Adams para la cavitaci´on descrito previamente. La posici´on del eje queda definida por el vector ~ OcOacuyas coordenadas polares son el ´angulo de desplazamiento del eje,φ, y la excentricidad dimensional del eje respecto del cojinete, e=Cε ( es habitual considerar dos sistemas de coordenadas: el sistema fijo (x, y) y el sistema (r, t) que se desplaza con el centro del eje). Para resolver anal´ıticamente la ecuaci´on de Reynolds, se puede introducir una simplificaci´on conocida como modelo de cojinete infinitamente largo. Esta aproximaci´on se basa en considerar que ∂p/∂θ es mucho mayor que ∂p/∂z, con lo que la ecuaci´on para la presi´on queda reducida a un caso unidimensional: 1 R2 a ∂ ∂θ h3 µ ∂p ∂θ = 6 (ω−2˙ φ)∂h ∂θ + 2 ˙ecos θ,(20) donde Raes el radio del eje, h=C(1 + εcosθ) la altura de la pel´ıcula de fluido (ver ecuaci´on (3)), µla viscosidad del fluido, pla presi´on, ωla velocidad de rotaci´on del eje, y˙ φ, ˙e, son las derivadas temporales de las coordenadas polares de ~ OcOa. Esta ecuaci´on (20) se corresponde con tomar el origen de ´angulos, θ, sobre la l´ınea de centros ¯ OcOaen el instante de tiempo considerado y, en consecuencia, es v´alida en el sistema de coordenadas m´ovil (r, t). Si adem´as se realiza la imposici´on de situar la zona activa (no cavitada) en 0≤θ≤π, se puede resolver (20) anal´ıticamente y obtener una expresi´on expl´ıcita de la 5 J. Durany, J. Pereira-P´erez, F. Varas presi´on: p= 6µRa C2"(ω−2˙ φ)εsenθ(2 + εcos θ) (2 + ε2)(1 + εcos θ)2+˙ε ε1 (1 + εcos θ)2−1 (1 + ε)2 −8 ˙ε π(1 −ε2)2(ε2+ 2) (ε2+ 2) 2arccos ε+cosθ 1 + εcosθ  −2εp1−ε2senθ 1 + εcosθ +ε2√1−ε2 2 senθ(ε+cosθ) (1 + εcosθ)2)#.(21) Una vez conocida la presi´on, p, se pueden determinar las componentes de la fuerza ejercida por el fluido sobre el eje, mediante la f´ormula: ~ F=Fr Ft=cos φ sen φ −sen φ cos φ·Fx Fy=ZL 0Zθ2 θ1 pcos θ senθ Radθdz, (22) siendo θ1, θ2los l´ımites de la zona activa, Lla longitud del par eje–cojinete y se ha introducido el cambio de variable correspondiente. A partir de las expresiones de Fr, Ft, en funci´on de ε, φ, ˙ε y ˙ φse realiza una aproximaci´on en serie de Taylor de la fuerza en el entorno del punto considerado: Fx Fy=−Axx Axy Ayx Ayy ·¯x ¯y−Bxx Bxy Byx Byy ·˙ ¯x ˙ ¯y,(23) que relaciona el desplazamiento del centro del eje y la velocidad de ese desplazamiento con la variaci´on de la fuerza. En el caso particular de disponer de la presi´on anal´ıtica (21), las componentes de la fuerza en coordendas (r, t) son: Fr=−12µLRaRa C2"ε2(ω−2˙ φ) (2 + ε2)(1 −ε2)+˙ε (1 −ε2)3/2π 2−8 π(2 + ε2)#,(24) Ft= 12µLRaRa C2"πε(ω−2˙ φ) 2√1−ε2(2 + ε2)+2ε˙ε (2 + ε2)(1 −ε2)#,(25) y los coeficientes matriciales de (23): Axx =2[8ε2(2 + ε4) + π2(1 −ε2)(2 −ε2+ 2ε4)] (2 + ε2)(1 −ε2)[4ε2+π2(1 −ε2)]3/2,(26) Axy =π[8ε2(2 + ε4) + π2(1 −ε2)(2 −ε2+ 2ε4)] ε(2 + ε2)√1−ε2[4ε2+π2(1 −ε2)]3/2,(27) Ayx =π[4ε4−π2(1 −ε2)2] ε√1−ε2[4ε2+π2(1 −ε2)]3/2;Ayy =2[4ε2+π2(2 −ε2)] [4ε2+π2(1 −ε2)]3/2,(28) Bxx =2[2ε2(π2(2 + ε2)−16) + π2(1 −ε2)(π2(1 −ε2) + 8ε2)] πε√1−ε2(4ε2+π2(1 −ε2))3/2,(29) 6 Estabilidad din´amica de pares eje–cojinete Bxy =Byx =2(π2−8)(ε2+ 2) (4ε2+π2(1 −ε2))3/2;Byy =π√1−ε2(π2−8)(ε2+ 2) ε(4ε2+π2(1 −ε2))3/2.(30) Cuando se buscan soluciones de tipo peri´odico para el sistema (23), se obtiene una condici´on de existencia para las mismas basada en el criterio de Routh-Hurwitz (ver [4]). Esta condici´on establece que para la ecuaci´on de la velocidad de rotaci´on adimensional del cojinete: ˜ω=ω0pMC/W0,(31) con Mla masa del eje y W0la carga soportada por el eje cuando se encuentra en equilibrio, existe un valor l´ımite para el par´ametro que interviene, ω0(velocidad cr´ıtica), dado por: ˜ω2 0=BxxByy −BxyByx (Bxx+Byy )(Axx Ayy−Axy Ayx) BxxAyy +ByyAxx −Bxy Ayx−Byx Axy −Bxx Axx+Byy Ayy +Bxy Ayx+ByxAxy Bxx+Byy .(32) De manera que para valores ˜ω < ˜ω0el movimiento es estable: ante una perturbaci´on respecto a la posici´on de equilibrio el sistema retorna a dicha posici´on. Para valores ˜ω > ˜ω0 las perturbaciones separan al eje de la posici´on de equilibrio estacionario y lo llevan a ciclos l´ımite, causando la inestabilidad. Sustituyendo las f´ormulas (26)– (30) en (32), se determina la velocidad cr´ıtica para distintos valores de la excentricidad. Con los datos: Ra= 0,0254, L= 0,0254, C= 0,000127, µ= 0,0254, W0= 150 N, y Mvariable, en la Figura 1 se observa esta velocidad cr´ıtica, de manera que las relaciones (ε, ˜ω0) situadas por debajo de la curva presentan situaciones de estabilidad, mientras que las que se sit´uan por encima son casos inestables. Adem´as, la velocidad cr´ıtica tiene una as´ıntota vertical y es un n´umero imaginario para valores altos de ε. En consecuencia, no se podr´ıa garantizar el comportamiento del par para valores de la excentricidad mayores que dicha as´ıntota. Cuando el estudio din´amico se realiza a partir del modelo hidrodin´amico de cavitaci´on de Elrod–Adams (4)–(8), se incluye el c´alculo de la fuerza ~ F, dada por (22), en cada iteraci´on temporal del m´etodo num´erico descrito en la Secci´on 2. Por otra parte, considerando el eje como un s´olido r´ıgido, se plantea la resoluci´on num´erica, mediante un m´etodo cl´asico de Newmark, del sistema de ecuaciones diferenciales ordinarias que rige la aceleraci´on del centro de gravedad del eje, (xg, yg), en funci´on de dicha fuerza: M¨xg=ZΓe p nxdΓ ; M¨yg=ZΓe p nydΓ,(33) siendo (nx, ny) las componentes del vector unitario, normal e interior a la frontera del eje. Los resultados obtenidos para los valores cr´ıticos de estabilidad se muestran en la Figura 1. Por ejemplo, para el valor de la excentricidad ε= 0,5, en la Figura 2 se puede comprobar el movimiento del centro del eje en tres comportamientos diferentes: caso estable, caso cr´ıtico y caso inestable, para los valores de ˜ω0= 1,10,2,5 y 3,66, respectivamente. Estos valores deber´ıan ser m´as realistas que los del caso anal´ıtico, precisamente por no imponer la posici´on de la zona cavitada. 7 J. Durany, J. Pereira-P´erez, F. Varas 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Figura 1: Velocidad cr´ıtica, ˜ω0, frente a la excentricidad ε(anal´ıtica: puntos; num´erica: c´ırculos) −0.5 0 0.5 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −0.5 0 0.5 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −0.5 0 0.5 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Figura 2: Movimiento del centro del eje (´angulo, excentricidad), para excentricidad inicial ε= 0,5. Caso estable, cr´ıtico e inestable, de izquierda a derecha. Agradecimientos Este trabajo ha sido parcialmente financiado por el Proyecto MTM2004-05796-C02-2 del MEC de Espa˜na, incluyendo fondos FEDER (UE), y el Proyecto PGIDIT05PXIC32202PN de la Xunta Galicia. Referencias [1] J. Durany, G. Garc´ıa, C. V´azquez. Numerical simulation of a lubricated hertzian contact problem under imposed load. Finite Elem. in Anal. and Design. Vol. 38 (2002), 645–658. [2] J. Durany, J. Pereira, F. Varas. A cell-vertex finite volume method for thermohydrodynamic problems in lubrication theory . Comp. Meth. Appl. Mech. Engrg., Vol. 195 (44-47) (2006), 5949-5961. [3] J. Durany, J. Pereira, F. Varas, Resoluci´on num´erica de problemas evolutivos de lubricaci´on termohidrodin´amica y estabilidad de pares eje-cojinete. En: Proceedings CMNE/CILAMCE 2007, APMTAC, Portugal, 2007. [4] J. Frene, D. Nicolas, B. Degueurce, D. Berthe, M. Godet. Lubrification hydrodynamique, Eyrolles, Paris, 1990. [5] K.W. Morton, M. Stynes, E. Suli. Analysis of a cell-vertex finite volume method for convectiondiffusion problems. Math. of Comput. Vol. 66 (220) (1997), 1389–1406, . 8