Análisis del flujo de calor en materiales graduados funcionalmente en contacto mediante el método de los elementos de contorno

MATE RI ALES COM P UESTOS 01 2 49 ANÁLISIS DEL FL U JO DE CALOR EN MATERIALES GRADUADOS FL"N C IO NALMEN TE EN CONTACTO ME DIANTE EL MÉTODO DE LOS ELEMENTOS DE CONTORNO J. P. Núñcz*.#, V. Ma nt ic *, F. Par ís *, J. Be r gc r *'' *Escuela Superi or de Ingenieros. Univ. de Se vi ll a. Camino de l os Descubt imicntos s/n, 41092 Sevilla **Colorado Sch oo l of Mines. Gol den, 1500 lllinois St reet, CO 8040 l, USA e-ma il: manti c@cs i.us.es, paris@es i. us.es, jberger@ mines.edu Palabras Clave: Materiales Graduados Funcionalmente (FGM5), Materiales Anisótropos, Conducción de calor, Revestimientos térmicos, Método de los Elementos de Contorno. Res ume n: Se desarrolla una aplicación del Método de los El ementos de Contorno para los problemas del flujo de calor estacionario en los sólidos anisótropos en contacto. Las conductividades térmicas en cada sólido se cons id eran constantes o con una variación exponencial en una dirección arbitraria pero fija. Los resultados numéricos presentados demuestran la utilidad del c ód igo computacional desarrollado para el estudio de las aplicaciones de los materiales graduados .f únciona/menle (FCMs) como las barreras !érmicas de los álabes en los motores térmicos some!idos a afias temperaluras. Abstrac t.· A Boundary Element ana~ )lsis of steady-state heat lransfer problems in anisotropic so/ids in contact is developed. Therma/ conductivities in each so/id are considered to be constan! or with an exponential variation in a .fixed but arbi/rQiy direction. N um er ical results presented demonstrate use.fulness o/ !he computational code developed for an analysis o.f app!ications of Functiona/ly Graded Materials (FCMs) as thermal barrier coatings in high-lemperalure turbine engine app/icalions. 1.- INT RODUCCIÓN La potencia que generan los motores térmicos (por ejemplo 1 as turbinas de gas o motores de aviones) de p en de fun da menta lmente de la temperatura máxima del ciclo que desc ri be el fluido que circula por el moto r. Sin emba rgo, esta t em peratura máxima está limitada por el material que conforma los álabes de la tur bi na el e tal forma que cualquier descubrimiento o avance en la mejora de las propiedad es de l os materiales y su resistencia a altas temper at uras in fl uye de ma nera im por tante en el rend im iento de los motores térmicos. "Ahora en GEC I, 5-7 rue Fran<;:ois Truffaut, 91000 Evry, France 25 0 MA TE RIAL ES COMPUES TO S 01 Para proteg er las partes metá li cas de un calo r excesivo en las aplicaciones de alta te mp eratura se han ut il izado l os revestimientos convencional es hom ogéneos de ma t er iales cerámicos. La dificultad en su utilización consiste en que, debido a las altas tension es residuales y tém1icas que apa r ece n en la interfase unidas a la baja r es istencia al desp egue de la misma así co mo a la ba ja tenacidad del revestimiento, estos materiales son propicios al ag ri etamiento y d es conchado, lo que co nduce a una pérdida de la protección tén nica y finalme nt e al fallo del co mpon ent e estru ctu ral '. En los nu evos di seños el e los revestimientos térmicos se apli ca n las capas o las entre -ca pas c on una c ompo sición graduada del metal y de la cerámica de tal forma que las pr op iedad es termomecánicas del r eves timiento c er ca del sub str ato metálico y las de l propio su bs trato metálico fueran simila re s. Estos nu ev os materiales se pueden considerar c om o materiales c ompu estos no ho mo gé neos de partículas en los cuales la fracción vol umé tri ca de Jos constituye nt es y posiblemente la micro -estructura del material va ría a tr avés del espesor. La aplicación de estos ma teriales, llamados material es graduad os funcionalmente (Functionally Gr ad ed Mat er ials - F GM s)2, permite reducir las tension es r es idual es e increm en tar la resistencia al desp eg ue y ademá s, debido al a um ento del contenido del metal en el revestimiento, se aume nt a su tenac id ad3. Su proceso de fabri cac ión se basa en un control exhaus ti vo de la porosidad, que se hará may or ce rca del substr at o y menor en la superficie el e tal forma qu e al difundir una aleación metálica por el ma terial se obtenga el ef ecto deseado. El pr esente trabajo está enmarcado dentro el e un trabajo má s amplio que es el desarro ll o de una herramienta de cá lc ul o para ana li zar proble ma s térmicos y term oe lásti co s planos que incluyen l os FGMs anisótropos, y en particular es tudiar el efecto el e la pr esencia de g ri etas y co ntacto en los especimenes sometidos a un fluj o de calo r. Debido a las bien conoc id as ventajas del Método de los Elementos de Co ntorno (M EC) para el análisis de l as g ri etas y el co nt act o4, sobre todo a su inherente precisión en la evaluación de los parámetros de los estados tensional es sing ul ar es , se ha optado por este método en el desarrollo de dicho có digo. En este ar tí c ul o se presenta la pr imera parte del trabajo d escr ito, que es el desarro ll o de un có digo para el análisis del flujo de calor estacionario en problemas con FGMs. El MEC es un méto do numé ri co apli ca do a la so lu ción de las Ec u ac iones Integrales de C ont orno ( EI C) cor respondient es a un problema dete rmin ado. La fornmlación ge neral de las EIC para el problema del fl ujo_ de ca lor est ac ionario en medios anisótropos homog én eos a trozos desarro ll ad a previamente' se ha implementado6 en un có digo del M EC más ge neral ya existente en el Grupo de Elas ti cidad y Re sistencia de Mat er ial es de la ESI de Sevi ll a. El énfasis de este artículo es tá en la fo rmulación de la EIC, y su implementación en dicho código, para los FGM anisótropos co n una va ri ación exponencial de las conductivi da d es térmicas. El en foque aquí pr esentado se podría f ác ilmente extender a ot ros ti pos de variación de las co nductivi da des 7, co mo po r ejemplo la cua dr ática. Finalmente se pr esentan soluciones numéri ca s de dos problemas. El primero co n una solución ana lí tica co n oc ida, para demostr ar un cor recto fu nci ona mie nt o del código desarrollado. El segundo repr ese nta una simulación senci ll a del problema de un sustrato homogéneo y un revestimiento de ti po FGM con un de sp eg ue asimilable a g ri e ta situada en la in te rf ase. M ATE RIALES CO MP UEST OS 01 25 1 2 .- EL PROBLEMA DEL FL U JO DE CALOR EN FGM EXPONEN CI ALES Se co nsidera el flujo de calor es taciona ri o en un sólido bidimensional D e 912 en el plano (x¡, x2) y con un conto rn o oD suave a trozos. El material del só li do se supone en general anisótropo con una variación de las conductividad es de f om 1 a: ku( x )= K u ex p{2 ~ · x } x E D. i,J=1,2 , donde el vect or constante ~ = (/] 1, /] 2) representa la dirección y mag nitud de la va 1' ia ció n. ~ se cons id era real o co mplejo imagi nar io. La mat ri z K es simét ri ca y positi va definida. Para mate ri ales isótropos K ij = K bij . El fl ujo de ca lor se calcula como h¡(x) =-ku (x )T_ 1(x) do nde T(x) es la temperatura del sólido. El flujo n orm al asoci ado a un vector normal unita ri o n(x) se eva lúa como q(x) = n¡( x)h¡(x ). La L ey de Fourier h¡ ¡(x) = Q( x ), d ond e Q(x) es la intensidad del flu jo de calor, se puede escribir en términos de temperatura como: - (K ii ~ ¡¡ (x) + 2/]¡ K ¡¡ r, 1 (x) ) e xp { 2 ~ · x} = Q( x) . (1) En cada parte del conto rn o oD = oDr uo Dq uo DR uoD cp u8 D ci definimos un tipo de condiciones de contorno o de contacto térmi co con otros sólidos. Así, tenemos: don de T y q son la temperatura y el fl uj o normal prescritos, To es la temperatura ambiental y Rr es el co eficiente de r es istencia térm ica (puede ser una var iable de punto). Adicionalmente: A B A B A "D q + q =O, T -T = Rrq en o ci donde o Dcp y oDci representan r es pectivamente las zon as de contacto perfecto e imperfecto en tre el sólido aquí ana li za do, seña lad o co mo A, y el otro señalado como B. 3.- DESARROLLO DEL MEC PARA LOS FGM EXPON EN Ci ALES 3.1.- Solución fundamental para los FGM exponenciales La solución fundamental es la solución en temperatura en todo el plan o originada p or una fuente de ca lor puntual4, Q(x) = 6 (x - y) , y E 912 , donde 6( x - y) es la función delta. La solución fundamental de la ecuación ( 1) tiene la expresión 8: 252 MATERIALES COMPUESTOS 01 K (KR)e - Jl ·(X+ ~ ) G(x,y)=- 0 2Tr det K donde y es el punto de la fuente de calor. x es el punto donde se evalúa la temperatura, K= pT Kp , R = r TK - 1r siendo r = x -y, Ko eslafuncióndc Bessel modifícadade orden cero. Para los materiales isótropos esta exp resión coincide con la obtenida previamente9. Se puede observar que C(x,y) = C(y, x ). El fl ujo normal correspondiente a C(x, y) se calcula como su derivada cononnal K: donde n = n(x) y K1 es la función de Be ssc l modifícada de primer orden. 3.2.- Ec uaci ón Int eg ral de Conto rn o para los FGM exponenciales Consideremos el esta do real en el teorema de Grccn de reciproci dad para el só li do D como aquél defínido en el contorno iJD por la temperatura T(x) y el flujo normal q(x) asociado a la normal exterior a 8D , y con Q(x) = O. Entonces si el estado auxiliar es el definido por C(x.y) y iJ G( x ,y) ovr para y E D uo D , se puede obtener con un procedimie nto estándar4·5·6 la EIC de temperatura para los FGM exponencial es : • cK (y)T(y) + f (oC (x ,y)T( x )-C( x, y)q( x)J ds x = O, (2) avr (D . MATERIALES COMPUESTOS 01 253 o K 11 =5, K 22 = l. K12= K 21 = O ~ , =- 1.5, ~ 2= 0 q0 -l 11 o o o . Fig. 1. Configuración de placa cuadrada de FG\11. colocación. Se ha utilizado un enfoque emianalítico para el cálculo de las integrales sobre los elementos que incluyen el punto de colocación. Los términos singular es obtenidos de los desarrollos asintóticos de las funcion es de Bessel se han integrado analíticamente y l os términos regulares restantes se han integrado numéricamente. 4.- EJEMP LOS 4.1.- Plac a cuadr a da de FGM ani sótro po ~e consi de ra un problema para una placa se9.ún la f¡gura 1 en el cual se conoce la solución ana lí tica~ : T(x 1 .x 2) = 100~ -e -2 P1 x1 ~ - e2 P1 J1 La solución en q es constante en toda la placa. El exponente fJ 1=-l.5 corresponde a una reducción de las co nductividades térmicas entre las caras x1=0 y x1=1 de valor 0.0498. El contorno se ha di sc retizado con 120 elementos uniformemente distribuido . 100 fO. E+OO "' -o 75 fo -3.E - 05 -o "' .~ 50 "' -6 . E05 E 25 o ·9.E-05 z o g LlJ -l.E -04 o 0.5 XI o 0.2 0.4 0.6 0.8 X¡ donde e K (y)= 1 si y E D, e K (y)= 0.5 si y es un punto en una parte suave de rD y Fig. 2. a) Evolución de T. b) Er ror nonnalizado de T. c , (y)= 1-1 so n{n(l )/\ n< 2l )arcco/ n<I ). Kn < 2 )) 1<. • 2 2Tr o ~ ~ .1<. .1\ ' si y es el vértice de una esquina de ( D con los vectores normales exteriores n (1) y n <2> . Resolviendo EIC (2) para y E cD de spués de aplicar las condiciones de contorno (y en su caso las de contacto) se obtienen lo s valor es de T y q en iJ D desconocidos. Una vez co nocidos T y q en t odo iJ D, se eva lúa T en D utilizando (2) para y E D. 3.3.- Impl eme nt a ción del có di go del MEC El algoritmo ele resolución numérica de la EIC (2) se ha i mp lementado 6 de acuerdo al procedimiento general del MEC de colocaciones de sc rito en París y Cañas 4. El contorno del só lido y las variables se han discrctiza clo mediante el eme ntos lineales continuos. La integrales se evalúan numéricamente sobre los elementos que no incluyen el punto de 100.- ------ ---- -------------- ~ 50 25 1 -o-Ana l itic~ --MEC 1 0+----r ---- r----.---,--~ o 0.2 0.4 ~ 0.6 0.8 1 ::r l.E+OO .--------------- ---, ~ 0 l.E-01 -o - ~ I.E-02 '" E l.E03 o z l.E· 04 l -o-Cara derecha J -- Cara i7quierdaJ g l. E05 +-----r----.---~----.--- ~ U.J o 0.1 0.2 0.3 0.4 0.5 Fig. 3. a) Evolución de q en la cara de recha. b) Va l or absoluto del error normalizado de q. Se presenta una comparación de la solución analítica en T y en q con la obtenida por el MEC. En la figura 2a se observa una excelente coincidencia entre la solución numérica y la analítica en T. Para conocer la precisión de la solución obtenida se presentan en la figura 2b los errores de T normalizados al valor máximo de T en la placa. Aunque en la figura 3a tampoco es apreciable la diferencia entre la so lución analíti ca y numérica en q, los errores 254 MATERIALES COMPUESTOS O 1 E E o a= 4mm o 11 v K= 2.09Wm ·2K·1 a .!:l 11 o K=25.5 1 Wm· 2K·1 T0= 295K L= l6mm Fig. 4. Con fi guración de sustrato con revestimiento. 0. 01 16 0.0087 0.0058 0.0029 500 T( K) + () Cru a izq uie r da -eCara derech a 1000 en q son ma y or es que l os de T. Se pueden obse rvar pequeñas oscilaciones debido a la proximidad de las esquinas. La m ag nitud de estas oscilaciones es aproximadam ente de un 1% del valor de q ana lí ti co , como se puede ver en la figura 3b. 4.2.- Sustrato isótropo homogéneo con revestimiento FGM isótropo Se cons id era una configuración de su strato con revestimiento sometida a una elev ada Ten la cara superior del revestimiento y con un despegue en la int er fase modelado co mo una grieta, figura 4. La conductividad térmica del revestimi ento coincide con la del sustrato en la interf ase. Las resistencias térmicas en las caras exteri ore s supe ri or e inf erior son R -¡== 5x L0 '4W.1 nlK y R -r== L04W.1m2K 1500 r espect ivamente. Aunque el códi go apli cado permita considera r una resistencia térmica finita y variable a lo largo de la grieta (que en un modelo má s realista podría ser función de la presión ele contacto en los labios de la grieta o de la separaci ón de los mismos ), l os resultados aquí presentados corresponden a una resistencia infinita, es decir la grieta actúa como un aislante perfecto. Para discre ti zar el sustr at o y Fig 5. Evolución de Ten las caras laterales. r eve stimiento se han utilizado 77 y 93 elementos de contorno respectivamente, discretizando la interfase c on una ma ll a confo rm e. 1600 0Cara in fcrí01 1400 --8C1ra su pe ri or ---AInter fase revest im iento 1200 - * - l n tc r t:1sc~us t rato 7 - 1000 - 800*- --~---~---*--~~ "" Xx -x.x - x-~ ¡= 600 4 00 .... .., ... . .., .. .. -<)-----<) -- --..,-----<:>- --- -<:>-- ---<:>----"'· -- -<> 200 0.004 O .OO R x ,( m) 0. 01 2 0.016 Fig. 6. Evolución de T en las cara s horizontales y en la interfase. Los resultados de T ob t en i dos en l as cara s laterales se pr es entan en la fi gura 5. En la ca ra derecha la g ri eta provoca un salto en T con pendient es nulas en los labios de la grieta de acuerdo con la condición de co ntorno, un aislamiento perfecto. La evolución encont r ada de T aproximadame nte lin ea l en el su strato y exponenc ial en el rev es timiento en la cara izquier da es tá de acuerdo con l as soluciones mono -dimensionales en cada tipo de material y se corresponde con el hecho de que la perturb ación del flujo de ca lor debido a la grieta se ha ya amorti guado a esta di st ancia. MATERIALES COMPUE ST OS 0 1 255 5.E+06 4 E+ 06 3 . E~06 2.E 106 -<:;- (;¡ ra i nfer i or -aCara ~ u pcrior --h;-ln tl!rf ase r rvcs1i mi emo -* [nt cr fase s u ~tratc · ·C: i. E ~06 k--_,;. &-~fr-&'~ 6="""'~~ ~ ~ Ob OO l -l= ~ >f'H3-Eif&-EH!H3-fB.¡;¡..{?Xl,tli~¡;Bg$!!illlllD "'-i.E+06 ¡.,_ j x : -2.E1-06 ~ ! -3.E ... 06 • .¡[<1)6 ~ -S E+06 ;------~--- -- -r-----'i-- -- ~ o 0.00 -l 0.008 x1 (m ) 0.012 0.016 fig. 7. Evo lución de q en las caras horizontales y en la in t er fase. i.E+09 H ' -+q E 3 R ecta de '5' regresi ón ~ 1.L+ 08 ~ ,,J q 584 1 4r-tJ ~ ~~ l. L-OS I. F. -07 I E06 1.E-05 r( m) Fig. 8. Evaluación del término s in g ul ar de q. En la figura 6 se puede observ ar el efecto de la g ri eta, primero en un au men to considerable de la T, alrededor ele 260 °C, en la ca ra superior del reves ti miento encima de la grieta, y segundo en la ev olución de l salto de T a lo largo de la grieta ele acuerdo con el desarro ll o asintó ti co 10 cuyo pr i mer término es proporcional a la fu nción raíz cuadrada de la distancia al vértice de la grieta. En la evolución horizontal de q, figura 7, la presencia de la grieta implica una disminuci ón del flujo en la ca ra supe r io r del revestimiento situada enci ma de la grieta, y en una singula ri dad de q en la interfase en el entorno del vértice de la grieta. Esta singularidad se analiza con más detalle en la figura 8, donde se puede ver un acu erdo excelente del orden de singularidad calculado (med iante una r egr es ión li neal de la evolución de q en función de la distancia al v ért ice de la grieta en escala logarítmica) con el predicho analíticamente 10 (-0.5). La exce le nt e con cor dancia de los resultados numéricos del MEC con l as características ele la solución ana lí tica conocid as apri or i permite supo ner que est os resultados numé ricos a pr oxi ma n adecuadamente T y q del problema pl anteado. 5.- CONCLUSIONES Se ha pr es ent ado la formulaci ón teór ica (l os núcl eos integrales y la Ecuación Int egr al de Conto mo c or re sp on dientes) y la implementaci ón de un cód igo del Método de l os Elementos el e Co ntorno (M EC) con l os element os li neales continu os pa ra los problemas de flujo de cal or es tacionar io en vari os só lidos anisótro pos hom~g~neo s o F GM ex p o~e?c1ales en contacto. Lo s resultados num éricos demu es tran que este cod1 go puede ser mu y u t!l en el estudio de los revestimientos térmi cos de tipo FGM, en particular en presencia de f allos de tipo despegue - gr ieta, q ue or iginan un cam bio br usco de resis t e~c i a térmi ca en la i nterf ase y por lo tanto una fuerte perturbación del flujo de calor. El trabaJO pr ese ntado representa la primera parte del desarro ll o en curso de un có digo del MEC para los problema s ter moe l ás ticos en las con ~ i g ur ac iones d esc ritas. 256 MATERIAL ES COMPUESTOS 01 6.- AGRADECIMIENTOS . Se agradece la ayuda recib id a de la Comisi ón de Intercambio Cultural, Educativo y Científico entre España y los EE.UU. de América (Proyecto del Programa Fulbright N. 99271) y de la DGESfC del Ministerio de Educaci ón y Cultura (Proyecto N. P B98-Il 18). 7.- REFERENCIAS . l. Miller T. C. ( 1998), Fi nite Element Analysis of a The rm ally Loaded Interface Crack in a Ceramic Coating. Eng in ecring Fracture Mechanics, Vol. 59, N° 2, pp 203-2 1 4. 2. Proceed in gs of the 6th In terna ti onal Symposium on Functionally Graded Materials (FGM2000), Estes Park, Colorado, USA, Sept. 2000 (in press). 3. Lee Y -D., Erdogan F. ( 199 8), Interface Crack in g of FGM Coatings under Steady-State Heat Flow. Enginecring Fracture Mechanics, Vol. 59, N° 2, pp 361-380. 4. París F., Ca i'í as J. ( 1997), Boundary Elemcnt Method. F un damentals and App li cations, Oxford Un ivcrs it y Press. 5. Mantic V., París F. (1995), On Free Terms and Singul ar Integrals in Tsotropic and A ni sotropic Potenti al Thcory. Computa ti onal Mechanics 95 , Atluri S. N., Yagawa G. , C ru sc T. A. , Springer, pp. 280628 11. 6. Núi'iez J. P. (2000), Aná li sis del Problema Tér mi co con Contacto para Medios Anisótropos con Variación Gradual de las Propiedades mediante el Método de los Elementos de Co ntorno, Proyecto Fin el e Carrera, ESJ , Universidad de Sevilla. 7. Clements D. L. , Budhi W. S. ( 1999), A Boundary Element Methocl of th e Solution of a Class ofSteady-State Problems for An isotropic Media. J. Heat Transfer, Vol. 121, pp. 462465. 8. Berger, J. R. , Martin P. A., Mantic V., Gray L. J. , Fundamental Solutions for SteadyStatc Heat Transfer in an Ex ponentia ll y-Graded A ni sotropic Material (submitted for publication). 9. Gray L. J. , Kaplan T., Richardson J. D. , Pa ulino G. H .. Green's Functions and Boundary Integral Equations for F un c ti ona ll y Graded Materials: Heat Conduction, J. Applied Mechanics (in press). 1 O. Mantic V., Pa rí s F ., Bcrger J., S in g ul ar ities in 2D Anisotropic Po tential Proble ms in Multimaterial Co rn ers. Real Variable Approach (submitted for publication).