Aplicación de criterios de orientación de grieta en problemas de crecimiento de grieta por fatiga bajo carga no proporcional

UNIVERSIDAD POLITÉCNICA DE VALENCIA Departamento de Ingeniería Mecánica y de Materiales TESIS DE MÁSTER APLICACIÓN DE CRITERIOS DE ORIENTACIÓN DE GRIETA EN PROBLEMAS DE CRECIMIENTO DE GRIETA POR FATIGA BAJO CARGA NO PROPORCIONAL Presentada por D. Francisco Gelardo Rodríguez Dirigida por Dr. D. Eugenio Giner Maravilla Valencia, Diciembre de 2014 - 2 - - 3 - A mis padres, hermanos, compañeros de piso y amigos. A Laura. - 4 - - 5 - - 6 - - 7 - RESUMEN Este trabajo se realiza como parte del Máster de Ingeniería Mecánica y de Materiales impartido en la Universidad Politécnica de Valencia. En este trabajo se hace un estudio de los parámetros que pueden afectar a la determinación de la orientación de la grieta en un problema de fretting fatiga bajo carga no proporcional con contacto completo. Los parámetros investigados han sido las cargas aplicadas, constante y variable, la relación de tensiones R que evalúa el nivel de tensión media y la diferencia de rigidez entre la probeta y el indentador. El criterio que se emplea para determinar la dirección de crecimiento de grieta es el criterio del mínimo incremento de la tensión tangencial ∆τ mín . A tal objeto se ha empleado el software de elementos finitos ABAQUS® para realizar el modelo, mallado y cálculo de resultados. Mediante el uso de distintas rutinas de Matlab® creadas por el Departamento de Ingeniería Mecánica, se han llevado a cabo los análisis para determinar la influencia de cada parámetro sobre la orientación de la grieta. Este trabajo se realiza mediante el MEF standard, ya que se trata de un análisis cualitativo del efecto de los distintos parámetros. Finalmente se realiza un análisis en X-FEM para confirmar los resultados y dejando abierta la posibilidad de realizar estudios posteriores en este campo. Palabras clave: Mecánica de la fractura, orientación de grieta, fretting fatiga, criterio del mínimo incremento de la tensión tangencial, MEF. - 8 - - 9 - ABSTRACT This work is done as part of the Master of Mechanical Engineering and Materials given by the Universidad Politécnica de Valencia. In this work a study of the parameters that can affect the orientation of the crack problem of fretting fatigue under nonproportional loading with full contact is done. The parameters investigated are the constant and variable load applied, the stress ratio R that assesses the medium stress level and the difference in stiffness between the specimen and the indenter. The criterion used to calculate the direction of crack growth is the criterion of minimum increase of shear stress ∆τ min . For this purpose the software ABAQUS® is used to model the geometry model, mesh and solve the problem. By using different Matlab® routines created by the Department of Mechanical Engineering, different analyses have been run to assess the influence of each parameter on the crack orientation. This work is done by using standard FEM, since this is a qualitative analysis of the effect of different parameters. Finally, an X - FEM analysis is performed to verify the results and provides the opportunity of further studies in this field. Keywords: Fracture mechanics, crack orientation, fretting fatigue, criterion of minimum range of shear stress, FEM. - 16 - Figura 4.12: ∆τ en cada substep del step 4. Caso 1. ________________________ - 67 - Figura 4.13: ∆τ en cada substep del step 4. Caso 3. ________________________ - 68 - Figura 4.14: ∆τ en cada substep del step 4. Caso 5. ________________________ - 69 - Figura 4.15: ∆τ en cada substep del step 4. Caso 7. ________________________ - 70 - Figura 4.16: ∆τ en cada substep del step 4. Caso 13. _______________________ - 71 - Figura 4.17: ∆τ en cada substep del step 4. Caso 13. _______________________ - 72 - Figura 4.18: Modelo mallado. ________________________________________ - 73 - Figura 4.19: Detalle de la malla en las inmediaciones de la zona donde se sitúa la grieta. Propagación 1. ______________________________________________ - 74 - Figura 4.20: ∆τ en cada substep del step 6. Propagación 1. _________________ - 75 - Figura 4.21: ∆σ en cada substep del step 6. Propagación 1. _________________ - 76 - Figura 4.22: Detalle de la malla en las inmediaciones de la zona donde se sitúa la grieta. Propagación 2. ______________________________________________ - 77 - Figura 4.23: Propagación de la grieta 1-2-3._____________________________ - 77 - Figura 4.24: Propagación de la grieta 4-5. ______________________________ - 78 - - 17 - 1 Introducción al problema 1.1 Presentación del problema El problema de la predicción y control de grietas es un tema de gran interés en la actualidad. Los costes de las materias primas y de la energía hacen indispensable el tratar de diseñar los componentes con el coeficiente de seguridad más bajo posible, asegurando la integridad estructural de los elementos que se estén diseñando. Es por ello que se hace necesario disponer de unas herramientas altamente fiables que permitan calcular con precisión la vida a fatiga de un determinado componente, para lo que es necesaria la correcta estimación y predicción de la dirección de propagación de las grietas que pueden aparecer bajo determinados estados tensionales. En este proyecto se llevará a cabo el estudio de orientación de grieta en problemas de fretting fatiga. El problema de fretting fatiga presenta un estado multiaxial de tensiones, en muchos casos con variación no proporcional, y requiere el empleo de criterios de fatiga multiaxial, adaptados a las particularidades del fretting. Los problemas en los que aparece fatiga bajo carga no proporcional están caracterizados por la aplicación de dos cargas las cuales no sufren la misma variación en el tiempo. 1.2 Objeto En problemas de fatiga 2D bajo carga no proporcional, la variación de los factores de intensidad de tensiones, K, en modo I y modo II no guarda la misma relación a lo largo de todo el ciclo, por lo que no es posible predecir una única dirección de propagación a lo largo del ciclo. Por tanto, no son de aplicación criterios muy utilizados, como el criterio de la máxima tensión circunferencial, MTS. Es necesario utilizar criterios que tengan en cuenta las variaciones de las magnitudes - 18 - relevantes (factores de intensidad de tensiones, tensión normal al plano de grieta, tensión tangencial) a lo largo de todo el ciclo. Además, para relaciones de tensiones de fatiga que impliquen cargas de compresión (p.ej. R=-1), aparece contacto con fricción entre las caras de grieta, cuyo efecto debe ser considerado y puede condicionar la dirección de propagación. El objeto de este trabajo es comprobar cómo afecta la variación de parámetros en la orientación de grieta una vez originada. Se realizarán simulaciones variando los parámetros de carga constante, carga alternante, R y rigidez de los materiales para un mismo valor de carga alternante, la carga constante. También se realiza para un mismo valor de carga constante, varios valores de R. 1.3 El problema de fretting fatiga En el ámbito de la ingeniería, fatiga es un término que se utiliza para definir la reducción de la resistencia mecánica al someter un material o componente a esfuerzos cíclicos que acaban por producir una reducción de la vida del componente en comparación con las propiedades intrínsecas del material. El problema de fretting fatiga se caracteriza porque las tensiones que originan y hacen crecer inicialmente las grietas son debidas al contacto entre dos componentes mecánicos. Habitualmente estos componentes en contacto, experimentan desplazamientos relativos de pequeña amplitud, lo que ocasiona un desgaste superficial conocido como fretting wear. El efecto de tensiones de contacto es análogo al de los concentradores de tensión: las altas tensiones cerca de la superficie hacen que la grieta se inicie. Además, en las primeras fases de crecimiento de grieta, estas mismas tensiones, unidas a las globales, provocan un crecimiento más acelerado. En la Fig.1.1 se muestra esquemáticamente la disposición de las fuerzas en una situación típica de fretting fatiga. - 19 - Figura 1.1: Esquema de las cargas que aparecen en fretting fatiga, donde A s es el área de la sección donde se aplica la tensión σ. La fuerza P mantiene en contacto los dos sólidos. La fuerza tangencial variable, Q, induce el deslizamiento entre las dos piezas. Generalmente, también existe una tensión global, σ, variable o no, aplicada a uno de los sólidos llamada “bulk”. En principio el fretting puede aparecer en cualquier máquina donde haya componentes en contacto como es el caso de los álabes y ejes acanalados de las turbinas y transmisiones, y en las uniones cubo-eje. Normalmente, es fácil de identificar a posteriori porque se observan marcas en las zonas que han estado en contacto. Éstas pueden tener un polvo característico o tener un aspecto compacto. El problema de fretting se aborda desde dos perspectivas distintas: • Fretting wear: Estudia el desgaste de las superficies en contacto entre dos cuerpos sometidos a cargas oscilantes y los efectos asociados. • Fretting fatiga: Estudia la iniciación y crecimiento de grietas por fatiga, donde además de la tensión aplicada sobre el componente se superponen las tensiones debidas al contacto. - 20 - Bryggman y Söderberg (1986) [3] distinguieron dos regímenes de deslizamiento: deslizamiento parcial y deslizamiento global. Concluyeron que, esencialmente, el fretting fatiga ocurre en las condiciones de deslizamiento parcial y el fretting wear en deslizamiento global. Para obtener información sobre la naturaleza del contacto, es útil clasificar los tipos de contacto: 1. Contacto completo: El área de contacto es independiente de la carga aplicada. 2. Contacto incompleto: El área de contacto es dependiente de la carga aplicada. Habitualmente, de forma experimental, se suelen estudiar dos configuraciones de contacto incompleto: cilíndrico (2D) y esférico (3D). A continuación se hará sólo un breve resumen del contacto completo y del contacto incompleto cilíndrico en dos dimensiones. 1.3.1 Contacto completo Hay muchos dispositivos mecánicos que involucran contacto completo, donde el fretting fatiga es el principal mecanismo del fallo. Un modelo de contacto completo se muestra en la Fig. 1.2 constituido por un indentador con un determinado ángulo en contacto con otro cuerpo de superficie plana. El indentador está sometido a una carga de compresión constante P, y al mismo tiempo a una carga tangencial cíclica Q paralela a la superficie del contacto y/o carga cíclica σ aplicada a la probeta. - 21 - Figura 1.2: Esquema de un problema de contacto completo. Para la configuración de la Fig. 1.2 de un indentador con esquinas a 90º sobre un semiplano infinito, la distribución de presión normal debida a una carga P por unidad de espesor, y en ausencia de fricción viene dada por la Ec. 1.1. 22 )( xa P xp − = π (1.1) Donde se observa que la ecuación anterior es singular en los extremos de la zona de contacto x=±a. En este esquema, 2a, define el ancho de la zona de contacto. - 22 - Figura 1.3: Distribución de tensiones normal a lo largo de la zona de contacto 2ª. Dicha distribución se muestra con línea discontinua y de forma aproximada en la Fig. 1.3 donde se observa que la tensión tiende a infinito al aproximarse a los bordes del indentador (Hills et al., 1993; Hills y Novell, 1994) [4, 5]. Muchos autores han estudiado el problema analíticamente, suponiendo distintas hipótesis simplificativas como por ejemplo tratar el problema con un indentador rígido, semiplano incompresible, indentador con deslizamiento o adhesión total, o indentador con esquinas redondeadas, etc. (Sackfield et al., 2001, 2002; Mugadu y Hills, 2002) [6, 7, 8]. Otros investigadores han propuesto resolver el problema utilizando un enfoque análogo al de la mecánica de la fractura, debido a la similitud del campo de tensiones alrededor de los bordes del modelo de dos grietas laterales en una placa. Este enfoque se denomina “crack like analogue” (Giannakopoulus et al., 2002; Conner et al., 2004) [9, 10]. En fretting es importante conocer las tensiones internas y en la superficie. Existen métodos de la Teoría de la Elasticidad que permiten el cálculo de las tensiones y los desplazamientos (Muskhelishvili, 1953) [11]. Hills y Nowell (1994) [5] desarrollaron estos métodos aplicándolos al problema de contacto completo. - 23 - 1.3.2 Contacto incompleto Como su nombre indica, el contacto incompleto cilíndrico es el que se produce entre dos cilindros. En los ensayos de fretting, en particular, uno de los cilindros tiene un radio infinito, es decir, se produce el contacto entre un cilindro y un plano, como muestra la Fig.1.4. En un ensayo de fretting se aplica una carga normal constante P que mantiene en contacto los sólidos y posteriormente se aplica la carga variable Q y la tensión σ. Debido a ellas, se distinguen dos tipos de zonas en la zona de contacto, una en el interior, de tamaño 2c, donde las superficies se mantienen adheridas, y otra en ambos extremos de la zona de contacto, donde se produce un deslizamiento parcial. Si la tensión axial σ es nula, la zona de adhesión estará centrada respecto a la zona de contacto. Si dicha carga es distinta de cero, la zona de adhesión se desplaza hacia un lateral una distancia e (excentricidad) (Hills y Nowell, 1994) [5]. Figura 1.4: Contacto cilíndrico con las fuerzas aplicadas; el semiancho de la zona del contacto es a, el de la zona de adhesión es c y la excentricidad es e. En esta geometría, a diferencia de las anteriores, la distribución de tensiones bajo la zona de contacto no presenta singularidades. Esta - 24 - geometría es una de las más usadas en ensayos de fretting si se desea tener una buena aproximación de las tensiones producidas en el contacto y al mismo tiempo tener la posibilidad de determinarlas analíticamente. Por el contrario, tiene la desventaja de que existen pocos casos reales de fatiga por fretting en los que aparezca esta geometría en el contacto. La expresión analítica para el campo de tensiones bajo el contacto en función de las cargas aplicadas en la superficie se puede encontrar en Johnson (1985) [12]. En función de los valores de Q y de σ, pueden aparecer dos casos muy distintos, uno en el que el deslizamiento se produce en el mismo sentido en todo el contacto y otro en el que los deslizamientos se producen en sentidos contrarios en las dos zonas de deslizamiento existentes (deslizamiento reverso) (Hills y Nowell, 1994; Tur et al., 2002) [5, 13]. 1.4 Hipótesis aplicadas Durante la realización de este proyecto se asumirá que el contacto se produce en la totalidad de la superficie entre ambos componentes, por lo que se tendrá un problema de contacto completo. Para modelar el problema se va a considerar una carga P constante y una tensión σ Bulk variable. Es habitual distinguir dos etapas claramente diferenciadas: nucleación de grieta y su posterior propagación. Debido a las fuertes tensiones en la zona de contacto es frecuente que los procesos de nucleación en problemas de contacto completo ocurran velozmente, consumiéndose la mayor parte de la vida en la fase de propagación. En la realización de este proyecto no se considerará la fase de nucleación o iniciación y se asumirá que la grieta está plenamente formada, con un tamaño suficiente para considerar su entorno como un medio continuo. Por tanto, para analizar la etapa de la propagación de grieta en estas condiciones es absolutamente necesario tener en cuenta la interacción contacto-grieta para la estimación de los FIT, lo que frecuentemente hace necesario el modelado numérico de esta interacción, por ejemplo, mediante el método de los elementos finitos extendido (X-FEM). - 25 - 2 Revisión de fundamentos 2.1 Introducción En este apartado se hará una revisión de las soluciones propuestas desde el inicio del estudio del problema de fretting-fatiga y se comentará la evolución de los diferentes métodos que han sido utilizados para su resolución. El Método de los Elementos Finitos se ha consolidado durante las últimas cuatro décadas como el método numérico más versátil para el análisis de problemas de la mecánica del sólido. Tras el establecimiento de las bases del método, muy pronto surgieron aplicaciones directas a la Mecánica de la Fractura (Watwood, 1969; Dixon y Pool, 1969) [14, 15]. Desde entonces el número de referencias en la literatura acerca de la aplicación del MEF y sus variantes como el X-FEM (Moës et al., 1999) [16] a la Mecánica de la Fractura ha crecido de forma imparable. A lo largo de los últimos 40 años han aparecido periódicamente revisiones de los métodos que permiten aplicar el MEF a la mecánica de la Fractura. La década de los 70 fue especialmente fructífera y pronto surgió la necesidad de revisar y ordenar la multitud de trabajos aparecidos y que establecieron la mayoría de los métodos disponibles hoy en día. Así se pueden destacar los trabajos de Rice y Tracey (1973) [17] y de Gallagher (1978) [18]. Posteriormente apareció el libro de Owen y Fawkes (1983) [19], de carácter introductorio y que incluye detalles acerca de la implementación de la mayoría de los métodos, y la detallada revisión editada por Atluri (1986) [20]. Raju y Newman (1984) [21] que presentaron un completo resumen acerca de la aplicación de métodos numéricos para el análisis de grietas 3D, incluyendo comparaciones en base a ejemplos. - 32 - Una tensión infinita no puede existir en un material real. Si la carga aplicada no es demasiado elevada, el material puede acomodar la existencia de una grieta inicial ideal de forma que la tensión teóricamente infinita se reduzca a un valor finito. En materiales dúctiles, como es el caso de muchos metales, aparecen grandes deformaciones plásticas en las inmediaciones del extremo de grieta. La región en la que el material fluye se denomina zona plástica. La deformación en el extremo de grieta da lugar a un extremo de grieta con un radio de curvatura pequeño (pero no infinitamente pequeño), de forma que el aspecto del extremo de grieta es romo. De este modo, la tensión no tiende a infinito, y la grieta se abre en su extremo una cantidad δ , denominada desplazamiento de apertura de extremo de grieta (CTOD). En todos los casos, el extremo de grieta experimenta una gran deformación y se desarrolla una separación finita en el extremo de grieta, se redistribuyen en una zona mayor. En el extremo de grieta se alcanza un valor finito de la tensión que puede ser resistido por el material, aunque a partir de una cierta distancia del extremo de la grieta, las tensiones son superiores a las correspondientes a la grieta ideal, de forma que se verifique el equilibrio global de cargas. En cualquier caso, en MFEL la zona plastificada es muy pequeña y queda englobada por los campos elásticos dominados por el FIT (hipótesis de small scale yielding). De todo lo anterior se deduce que en el Planteamiento Local de la MFEL es muy importante encontrar expresiones explícitas de K. Este a su vez depende de la configuración y geometría del problema, incluyendo la propia longitud de grieta a, como se puede ver en la siguiente ecuación: aCK nom ⋅= πσ Donde C es el llamado factor geométrico, siendo un parámetro dependiente del modo de apertura de la grieta, el tipo de carga aplicada y obviamente, de la geometría del componente analizado. La tensión nominal σ nom también depende del problema considerado y de las solicitaciones (flexión, torsión,…). En esta ecuación, a es el tamaño de (2.6) - 33 - grieta. Es importante remarcar que el FIT tiene unidades [ mMPa ] en el SI. Las ecuaciones 2.1, 2.2 y 2.3 son básicas en MFEL y como se observa están descritas en función de los Factores de Intensidad de Tensiones (FIT) como únicos parámetros caracterizantes. Cuando a finales de los años 60 el MEF comenzó a ser aplicado a problemas de la MFEL, pronto surgió la necesidad de representar correctamente el estado tensional dado por las ecuaciones 2.1, 2.2 y 2.3, ya que la formulación tradicional del MEF no está especialmente indicada para el modelado del comportamiento singular. Durante los últimos 50 años se han desarrollado numerosos métodos para la determinación del factor de intensidad de tensiones. Desde que Irwin estableció que el valor de los FIT caracteriza de forma unívoca el estado tensional en el entorno del extremo de una grieta en MFEL, su evaluación ha sido un objetivo prioritario en la aplicación de la Mecánica de la Fractura, lo que ha dado lugar a una gran diversidad de técnicas disponibles. Muchos de los planteamientos iniciales, de carácter analítico, han sido en la actualidad superados por la versatilidad que ofrecen los métodos numéricos. En este capítulo, la atención se centra fundamentalmente en este último tipo de métodos, y en particular, en aquellos relacionados con el empleo del Método de los Elementos Finitos (MEF) y del Método de los Elementos Finitos Extendido (X-FEM). Una visión global de las bases para la determinación de los FIT se puede encontrar en la colección de trabajos editada por Sih (1973) [28], muchos de ellos de carácter analítico, y sobre todo, en el libro de Aliabadi y Rooke (1991) [29] (con una parte sustancial orientada a la aplicación del método de los elementos de contorno a la Mecánica de la Fractura) o el completo resumen de Rooke (1994). 2.2.4 Planteamiento global de la MFEL El primer enfoque utilizado en el estudio de la propagación de una grieta en un cuerpo cargado con comportamiento elástico lineal fue el llamado “planteamiento global” o equivalentemente “planteamiento - 34 - energético”. En 1921, A.A. Griffith [22] publicó sus trabajos en los que utiliza el concepto clave de tasa de liberación de energía (strain energy release rate) denotada con el símbolo G en su honor. A veces, G es denominada también velocidad de relajación de energía. Figura 2.4: Cuerpo cargado con grieta de longitud a y superficie de grieta a·B. Para comprender el planteamiento energético, vamos a considerar el cuerpo de la Fig.2.4 bajo la acción de unas cargas exteriores que lo deforman, este almacenará una energía potencial. Si se supone que contiene una grieta de longitud a, en el momento en el que esta grieta avance un cierta cantidad ∆ a, cambiará su geometría, la distribución de sus tensiones, los puntos de aplicación de las cargas, etc. En general, ese cambio supondrá una variación en la energía disponible: • Si la variación de energía disponible es igual o mayor que la necesaria para romper la cohesión del material existente en el extremo de grieta, esa grieta progresará. Puede ser que el crecimiento sea de forma inestable, propagándose rápidamente y ocasionando en último término la rotura total de la pieza. • Si la variación de la energía disponible es menor que la necesaria para romper el material, la grieta no progresará. - 35 - Griffith expresó esta idea del siguiente modo: “El crecimiento de grieta sólo puede ocurrir si la energía requerida para formar nuevas superficies de grieta dA puede ser suministrada por el sistema”. Se debe entender ”sistema” como el conjunto del propio sólido que contiene la grieta como las cargas exteriores que actúen sobre el cuerpo. Esas “fuentes” de energía necesaria pueden tener el siguiente origen: 1. El sólido es capaz de proporcionar energía liberando parte de su energía de deformación elástica. 2. Las cargas exteriores son capaces de proporcionar energía a partir del trabajo que desarrollan cuando desplazan su punto de aplicación. Para formalizar una expresión matemática de G, se considera el caso más general tridimensional en el que un cuerpo (elástico o no) está sometido a unas cargas y presenta una grieta con un área A. Como se ha comentado, si las cargas cambian con el tiempo es posible que la grieta avance; por la ley de conservación de la energía es necesario que se cumpla para un cuerpo en equilibrio: Γ+++= & &&&& KUUT pe Donde T es el trabajo realizado por las fuerzas exteriores (supuestas constantes); U e es la energía potencial de deformación elástica; U p es el trabajo realizado en el caso que exista deformación plástica; K es la energía cinética del cuerpo; Γ es la energía consumida en generar el área de grieta, por rotura de la estructura del material en el extremo de la grieta. Como todos los cambios con respecto al tiempo son debidos a un cambio en el área de grieta, se puede escribir: dA d A dt d = (2.7) (2.8) - 36 - Por tanto es equivalente hablar de variación con respecto al tiempo que variación con respecto al área de grieta. La ley de conservación de la energía queda: dA d dA dK dA dU dA dU dA dT pe Γ +++= En el supuesto de cargas constantes con el tiempo y si la grieta se supone que crece lentamente, se considera el problema como cuasi estático y se puede despreciar el término que tiene en cuenta la energía cinética K. La ecuación se puede reordenar como: dA d dA dU dA UTd pe Γ += − )( En elasticidad se define el término de la energía potencial total ( Π ) de un sistema a la diferencia: TU e −=Π Y por tanto: dA d dA dU dA d p Γ += Π − El sentido físico implícito en esta ecuación es que la variación (decrecimiento) en el valor de Π cuando crece una grieta es igual a la variación de la energía consumida en deformación plástica y generación de nuevas superficies de grieta. En el caso elástico, el trabajo consumido en deformación plástica se puede considerar despreciable y la conservación de la energía queda dA d dA d Γ = Π − (2.9) (2.10) (2.11) (2.12) (2.13) - 37 - El primer término de la ecuación anterior es la definición formal de G, es decir la tasa de liberación de energía por unidad de área de grieta: dA d G Π −= Por otro lado, el término del lado derecho tiene que ver con la formación de nuevas superficies de grieta. A veces se escribe de las siguientes formas equivalentes: G=R o bien G=Gc Donde R se denomina tenacidad a fractura. Es decir, G representa la energía disponible para el crecimiento de grieta y R (Gc) representa la resistencia del material que debe ser vencida para que la grieta progrese. Es lógico pensar que R es una propiedad de cada material. De lo anterior se deduce el siguiente criterio de fallo en Mecánica de la Fractura Elástico Lineal: • Si G < R la grieta no llega a progresar. • Si G ≥ R la grieta crece. Su crecimiento puede ser estable o inestable. 2.3 Método de los Elementos Finitos (FEM) Para caracterizar las ecuaciones 2.1, 2.2 y 2.3 a través de los Factores de Intensidad de Tensiones (FIT) es necesaria la implementación de un método que permita obtener el resultado de los FIT teniendo en cuenta la singularidad producida por la propia grieta. Para la aplicación del Método de los Elementos Finitos aplicado a la simulación del crecimiento de grieta se aplica el siguiente esquema: (2.14) (2.15) - 38 - • Modelado de la geometría del problema. • Generación de la malla, teniendo presente la grieta. Por lo tanto se tendrá que realizar una malla que llegue hasta una de las caras de la grieta, la bordee y continúe por la otra cara, generando la discontinuidad como puede apreciarse en la Fig. 2.5. Hay que tener en cuenta que un factor determinante será el refinamiento en torno a la singularidad, con el correspondiente coste computacional asociado. Figura 2.5: Detalle de malla según el planteamiento FEM • Aplicación del FEM para la resolución de las ecuaciones del apartado anterior. Se obtienen sendos valores para el Factor de Intensidad de Tensiones a través, por ejemplo, de técnicas de extrapolación de tensiones, extrapolación de desplazamientos o aplicación de integrales de dominio. • Resolución del criterio correspondiente (en función de los FIT) para la obtención del ángulo que seguirá el siguiente incremento de grieta. • Generación de la nueva geometría añadiendo el nuevo incremento de grieta con su correspondiente ángulo calculado en el paso - 39 - anterior. Notar que se ha de desechar la geometría anterior con su malla asociada. • Mallado de la nueva geometría teniendo en cuenta los 2 incrementos de grieta. • Y así sucesivamente hasta alcanzar la longitud de grieta deseada o provocar la rotura de la pieza. Como puede observarse, el Método de los Elementos Finitos aplicado a la MFEL requiere de un gran consumo de tiempo, tanto de interacción con el usuario como computacional ya que para cada nuevo incremento de grieta se ha de generar una nueva malla que tenga en cuenta el nuevo tramo de grieta que se ha calculado para posteriormente calcular los FIT y aplicar el criterio correspondiente de orientación de grieta. 2.4 Método de los Elementos Finitos Extendido (X-FEM) En la formulación convencional de elementos finitos, la existencia de una grieta se modela explícitamente mediante la frontera de los elementos. En contraste, en el método X-FEM los lados de los elementos no tienen por que coincidir con la posición de la grieta, lo que proporciona una gran versatilidad. Los métodos que no requieren mallas, como el “Element Free Galerkin Method”, o métodos sin malla (Belytschko et al., 1994) [30], se empezaron a utilizar en la mecánica computacional para dar solución a problemas que representaban dificultades con el MEF. No obstante, estas técnicas requieren cálculos computacionales más complejos, relacionados con la generación de funciones representativas y cuadraturas adicionales. Además presentan dificultades para satisfacer las condiciones de contorno de Dirichlet. El método X-FEM se basa en el enriquecimiento del modelo de elementos finitos con grados de libertad adicionales en los elementos geométricamente intersectados por la grieta. De esta forma la discontinuidad se incorpora sin modificar la discretización de la malla, que es generada sin considerar la posición de la grieta. Obviamente, en la - 40 - implementación del X-FEM, es necesario conocer topológicamente la posición de la grieta respecto a la malla. Con este fin, se utiliza la técnica LS (Level Set Method) [31] para caracterizar los elementos y nodos afectados por la grieta (denominados nodos y elementos enriquecidos). En la Fig. 2.6 se muestra una porción de la malla utilizada en este trabajo mostrando con círculos los nodos enriquecidos con 2gdl adicionales (total 4 gdl) y con cuadrados los nodos enriquecidos con 8 gdl adicionales (total 10gdl). Los elementos enriquecidos son aquellos que contienen al menos un nodo enriquecido. Figura 2.6: Detalle de una malla X-FEM Los nodos con 2 gdl adicionales (uno para cada dirección del plano) tienen definidas funciones de forma que incluyen la función de Heaviside H(x) (módulo unitario y cambio de signo en la cara de la grieta). Físicamente, esta función introduce la discontinuidad entre caras de grieta. Los nodos con 8gdl adicionales son enriquecidos en las dos direcciones del plano con 4 funciones F j (x) que reproducen el comportamiento singular de la MFEL en tensiones. De esta forma, en el caso bidimensional, la interpolación de elementos finitos, considerando un punto de coordenadas x, resulta: ∑ ∑ = =      ++= malla nn i j i jjiiief bxFaxHuxNxu 1 4 1 )()()()( (2.16) - 41 - Donde nn malla es el número total de nodos de la malla, N i (x), u i son las funciones de forma y gdl convencionales de cada nodo i y a i , j i b los gdl de libertad adicionales asociados a las funciones de Heaviside H(x) y de extremo de grieta F j (x). Es importante indicar que en la ecuación los gdl adicionales a i y j i b sólo se añaden para aquellos nodos que son enriquecidos, según la topología grieta-malla. Como sucede en el MEF, es necesario realizar integraciones numéricas en el dominio del elemento para el cálculo de la matriz de rigidez. Sin embargo, el hecho de que exista la discontinuidad debida a la grieta, exige dividir previamente los elementos intersectados por ella en subdominios en los que la grieta sea uno de sus lados. - 48 - Las tensiones máximas principales tienen dirección horizontal en los extremos de la barra determinada por la tensión aplicada como se muestra en la Fig.3.5; excepto en la zona de la grieta; donde ésta se abre como se puede apreciar en la Fig.3.4 creándose una discontinuidad y un concentrador de tensiones que provoca que las líneas de fuerza se concentren en el frente de grieta. Figura 3.5: Barra sometida a tracción. Orientación y magnitud de las tensiones máximas principales. En la Fig.3.6 se observa el valor de las tensiones normales en relación con el ángulo en el que se toman en un elemento situado en el frente de grieta en el último step. También se representa ∆τ 12 . El criterio para establecer el origen de los ángulos se puede ver en la Fig.3.7. En éste gráfico queda patente que el caso de carga es proporcional, porque la máxima y la mínima tensión se producen en la misma dirección y en la gráfica se puede trazar una vertical que una los máximos y mínimos de cada substep. - 49 - -100 -80 -60 -40 -20 0 20 40 60 80 100 -20 -10 0 10 20 30 40 50 60 70 1 2 3 4 5 6 7 8 σ 2 max σ 2 min efec σ 2 min no efec ∆ σ 2 efec ∆ σ y no efec ∆ τ 12 Figura 3.6: Tensiones normales en el fondo de grieta y su variación con el ángulo de estudio θ (º) En la Fig.3.6, los valores mostrados corresponden a: - Los valores de las tensiones normales σ 2 para cada uno de los 8 substeps, numerados a la derecha de la gráfica, que componen el último step en el que se aplica una carga de tracción de 20MPa. - La tensión σ 2 máxima que coincide con la tensión del substep 8. - La tensión σ 2 mínima efectiva, al ser los valores de la tensión mínima negativos para todos los ángulos en el substep 1 por estar aplicando una carga de compresión, la tensión efectiva es 0. - La tensión σ 2 mínima no efectiva que coincide con la tensión del substep 1 cuando la σ b es todavía de compresión y está cerrando la grieta. - 50 - - El incremento de la tensión efectiva ∆σ 2 efec; es la diferencia entre la tensión máxima y mínima efectiva, en este caso al ser σ 2 min=0, ∆σ 2 efec= σ 2 max. - El incremento de la tensión no efectiva ∆σ y no efec= σ 2 maxσ 2 min; en este caso es diferencia entre las tensión en el substep 8 menos el substep 1. - El incremento de la tensión tangencial ∆τ 12 . Figura 3.7: Criterio de orientación de grieta. El criterio para establecer la orientación del ángulo es como se indica en la Fig.3.7, valor entre 0 y 90º para los ángulos situados en el 3 er cuadrante y entre 0 y -90º para los ángulos situados en el 4º cuadrante. El origen se establece en el fondo de grieta. La Fig.3.8 muestra las tensiones máximas en azul y las tensiones tangenciales en rojo que se producen en un elemento situado en el extremo de grieta, en el último step, para cada substep. El último step es un ciclo que comienza con la pieza sometida a compresión, -20MPa, y acaba en tracción a 20MPa. De ahí que durante los 3 primeros substeps las tensiones en el eje x sean de compresión. En el 4º substep, la tensión a la que se somete a la pieza es 0, con lo que las tensiones son nulas. Es a - 51 - partir del 5º substep cuando la pieza es sometida a tracción y las tensiones principales se hacen positivas y producen la apertura de la grieta. Figura 3.8: Valor de las tensiones normales y tangenciales frente al ángulo del elemento diferencial en el extremo de grieta. El ángulo θ de la Fig.3.8 es la orientación de un elemento diferencial en el extremo de la grieta que muestra el ángulo en el que se estudian las tensiones al ir variando el ángulo θ , de acuerdo al criterio de orientación visto en la Fig.3.7. - 52 - Figura 3.9: Representación de las tensiones en un elemento diferencial. La representación del ∆τ 12 para las orientaciones comprendidas entre 90 y -90º en el frente de grieta en el último step de carga se puede ver en la Fig. 3.10. En éste gráfico queda muy patente si el caso de carga es proporcional, como sucede en este caso, donde la máxima y la mínima tensión siempre se producen en la misma dirección y en la gráfica se puede trazar una vertical que una los máximos y mínimos de cada substep; o si el caso es no proporcional como pasa en el caso de tener una carga constante y una carga variable; es el caso de las simulaciones de fatiga fretting que se expondrán en éste trabajo. De acuerdo al criterio ∆τ min , la orientación de grieta es a 0, 90º, - 90º. Se descarta -90º porque en esa dirección no hay material por el cual pueda progresar la grieta. Y entre 0 y 90º, el criterio establece que el ángulo en el que la tensión principal máxima es mayor determina la orientación, que según la Fig.3.10, es a 90º. - 53 - Como conclusión la grieta progresa a 90º, totalmente vertical, dominada por el modo I de apertura de, tal como se puede comprobar experimentalmente. -100 -80 -60 -40 -20 0 20 40 60 80 100 -15 -10 -5 0 5 10 15 20 25 θ (º) [MPa] 12345678 τ12 max τ12 min ∆ τ12 Figura 3.10: Valor de las tensiones tangenciales frente al ángulo del elemento diferencial en el extremo de grieta. - 54 - 4 Modelado y resultados. 4.1 Modelo numérico MEF El modelo de fretting fatiga analizado en este proyecto corresponde al problema en condiciones de contacto completo de 2D esquematizado en la Fig. 4.1. Figura 4.1. Problema fretting-fatiga con contacto completo. Las dimensiones del modelo son h=2c=2B=10mm, 2L=40mm. El coeficiente de fricción que se ha tomado para modelar el contacto entre el indentador y la probeta es µ int =0,8 [39] [40]. Con objeto de minimizar el número de elementos empleados, teniendo en cuenta que la geometría de la pieza es simétrica tanto en el eje X como respecto del eje Y se modela sólo un cuarto de pieza y se aplican condiciones de contorno de simetría en la línea vertical izquierda (U1=0), y en la línea horizontal inferior (U2=0). Esta simplificación es - 55 - válida según los diferentes ensayos correlacionados con las simulaciones realizadas [34] [35]. La Fig.4.2 muestra una máquina de ensayos de tracción donde se realizan los ensayos de fretting-fatiga de contacto completo. Se pueden ver los utillajes empleados para garantizar la carga constante P aplicada en los indentadores en dirección perpendicular a la probeta. Figura 4.2: Máquina de ensayos de contacto completo, mostrando los elementos en contacto. El material de la probeta en forma de cruz con las dimensiones indicadas, está fabricada en aluminio EN AW-7075-T6 según norma EN485-2. Las características mecánicas de este aluminio para el espesor dado son: Estado de tratamiento Espesor (mm) Resistencia a la tracción R m (MPA) Límite elástico R p0,2 (MPA) Alargamiento mín. A50 (%) Dureza HBW T6 5 545 475 8 163 Tabla 1. Propiedades del aluminio 7075 T6. - 56 - El módulo de Young considerado para este material es de 72GPa. Las dimensiones de la probeta permiten que el concentrador de tensiones no afecte a los extremos de la pieza, donde se verifican las condiciones de contorno aplicadas. La dimensión de la grieta inicial es de 0.3mm. Una vez definido el modelo se realiza el mallado empleando las herramientas de ABAQUS. Figura 4.3: Modelo mallado. En la siguiente figura se puede observar un zoom de la malla anterior en la zona refinada donde se sitúa la grieta, que se puede ver justo en la esquina. - 57 - Figura 4.4. Detalle de la malla en las inmediaciones de la zona donde se sitúa la grieta. La malla se ha hecho de tal modo que en el contorno de la grieta, se tiene un mallado estructurado de elementos cuadriláteros cuadráticos CPE8R de tamaño a 0 /20, que en este problema es de 0,015mm. En las zonas colindantes se hace un mallado en barrido hasta alcanzar un tamaño de elemento de 0,5mm. En los extremos se vuelve a realizar un mallado estructurado de tamaño 0,5mm. Esto permite limitar el número de elementos del modelo sin perjuicio de su representatividad y acortar el tiempo empleado en ejecutar cada cálculo. Al estar la malla formada por elementos más pequeños en la zona donde se produce la singularidad y que se quiere estudiar y de mayor tamaño en las zonas alejadas que no son objeto de estudio como en el caso de un refinamiento h-adaptativo. El modelo completo contiene un total de 61985 nodos y 20444 elementos. 4.2 Aplicación del ciclo de carga. En las simulaciones realizadas en este proyecto, se han considerado 4 pasos de carga para llevar a cabo cada análisis. En el primer paso se aplica la carga constante σ P hasta alcanzar su valor máximo; esta carga permanece constante durante el resto de pasos (steps) siguientes. En los siguientes pasos se aplica la carga cíclica en la probeta - 64 - Figura 4.9: Caso 15. Orientación de las tensiones máximas principales. Los casos 16 al 21 hacen un barrido aplicando la misma carga constante y alternante σ P y σ B, pero incrementando la rigidez del indentador de 20000 MPa a 600000 MPa para ver cómo la el ángulo de la grieta va disminuyendo por debajo de 90º hasta los 77º donde se estabiliza al asignar la rigidez del aluminio tanto al indentador como a la probeta. En los casos 22, 23 y 24 lo que se hace es mantener muy alta la rigidez del indentador, 1000000 MPa y variar la carga alternante que disminuye hasta los 10 MPa para constatar que pese a tener una rigidez muy alta y aplicar una carga constante también muy alta frente a la carga alternante horizontal, esto no repercute en la orientación final de la grieta. En la Fig.4.10 se puede observar el caso 21 donde aplicando en un indentador de E=600000 MPa, con una carga constante de 200 MPa y disminuyendo la carga alternante a 50 MPa, el ángulo se sigue manteniendo en 79º. - 65 - -100 -80 -60 -40 -20 0 20 40 60 80 100 -400 -300 -200 -100 0 100 200 300 400 θ (º) [MPa] 1 2 3 4 5 6 7 8 τ12 max τ12 min ∆ τ12 Figura 4.10: ∆τ en cada substep del step 4. Caso 21. 4.3.3 Proporcionalidad del caso de carga Una de las características que define el problema de fretting fatiga modelado es que se trata de un problema de cargas no proporcionales. Esto se debe a que la carga que se aplica en el indentador es constante, mientras que la carga aplicada en la probeta de ensayo es alternante. - 66 - En el caso que se ha mostrado de ejemplo de comprobación de la metodología, la barra sometida a tracción-compresión, ésta está sometida únicamente a una carga alternante, por lo que la tensión tangencial máxima se da durante todos los substeps en la misma orientación del elemento diferencial a 45º haciéndose el ∆τ =0 a 0 y 90º, donde se produce la máxima tensión principal. -100 -80 -60 -40 -20 0 20 40 60 80 100 -15 -10 -5 0 5 10 15 20 25 θ (º) [MPa] 12345678 τ12 max τ12 min ∆ τ12 Figura 4.11: ∆τ en cada substep del step 4. Caso barra sometida a tracción-compresión. Conforme en los casos aumenta la tensión constante σ P , se observa que el desfase entre la máxima y la mínima tensión tangencial τ aumenta. Resultando este desfase en una medida de la no - 67 - proporcionalidad del caso de carga. Esto se aprecia en los ejemplos siguientes. Para el Caso 1, Fig.4.12; siendo σ P (MPa) =1e-6 / σ B (MPa) = 200 / R = -1 se mide el ángulo en el que se produce el máximo de la tensión tangencial en el substep 1, es 32 grados, la mínima en el substep 8 se produce a 34 grados. -100 -80 -60 -40 -20 0 20 40 60 80 100 -200 -150 -100 -50 0 50 100 150 200 250 θ (º) [MPa] 1 2 3 4 5 6 7 8 τ12 max τ12 min ∆ τ12 Figura 4.12: ∆τ en cada substep del step 4. Caso 1. - 68 - En el Caso 3, Fig.4.13; siendo σ P (MPa) =50 / σ B (MPa) = 200 / R = -1 se mide el ángulo en el que se produce el máximo de la tensión tangencial en el substep 1, es 15 grados, la mínima en el substep 8 se produce a 42 grados. -100 -80 -60 -40 -20 0 20 40 60 80 100 -200 -150 -100 -50 0 50 100 150 200 250 θ (º) [MPa] 1 2 3 4 5 6 7 8 τ12 max τ12 min ∆ τ12 Figura 4.13: ∆τ en cada substep del step 4. Caso 3. - 69 - En el Caso 5, Fig.4.14; siendo σ P (MPa) =100 / σ B (MPa) = 200 / R = -1 se mide el ángulo en el que se produce el máximo de la tensión tangencial en el substep 1, es 5 grados, la mínima en el substep 8 se produce a 49 grados. -100 -80 -60 -40 -20 0 20 40 60 80 100 -200 -150 -100 -50 0 50 100 150 200 250 θ (º) [MPa] 1 2 3 4 5 6 7 8 τ12 max τ12 min ∆ τ12 Figura 4.14: ∆τ en cada substep del step 4. Caso 5. - 70 - En el Caso 7, Fig.4.15; siendo σ P (MPa) =200 / σ B (MPa) = 200 / R = -1 se mide el ángulo en el que se produce el máximo de la tensión tangenciales en el substep 1, es -3 grados, la mínima en el substep 8 se produce a 59 grados. -100 -80 -60 -40 -20 0 20 40 60 80 100 -250 -200 -150 -100 -50 0 50 100 150 200 250 θ (º) [MPa] 1 2 3 4 5 6 7 8 τ12 max τ12 min ∆ τ12 Figura 4.15: ∆τ en cada substep del step 4. Caso 7. - 71 - En el Caso 13, Fig.4.16; siendo σ P (MPa) =200 / σ B (MPa) = 10 / R = -1 se mide el ángulo en el que se produce el máximo de la tensión tangencial en el substep 1, es -15 grados, la mínima en el substep 8 se produce a 74 grados. -100 -80 -60 -40 -20 0 20 40 60 80 100 -200 -150 -100 -50 0 50 100 150 200 θ (º) [MPa] 1 234 567 8 τ12 max τ12 min ∆ τ12 Figura 4.16: ∆τ en cada substep del step 4. Caso 13. En este último caso, se observa que al ser un caso donde la carga σ B << σ P , es prácticamente un caso estático donde no existe fatiga como se puede apreciar por el incremento de la tensiones tangenciales ∆τ≈ 0. - 72 - 4.4 Orientación de grieta aplicando XFEM. Para ilustrar el estado actual en el que se encuentra el cálculo y predicción de la orientación de grieta en la ingeniería mecánica, se ejecuta una simulación utilizando el método de elementos finitos extendido, abreviado XFEM. Para ello se modela una probeta que tiene las dimensiones indicadas en la Fig.4.17: Figura 4.17: ∆τ en cada substep del step 4. Caso 13. Las dimensiones del modelo son h=16mm; c=2B=10mm; L=50mm y con un espesor t de 1mm. El coeficiente de fricción que se ha tomado para modelar el contacto entre el indentador y la probeta es µ int =0,8 como se ha tomado anteriormente en los modelos. Se crea una grieta de tamaño a 0 =0,1mm e inclinada ésta 60º sobre la horizontal. El material empleado en esta simulación es aluminio EN AW7075-T6 según norma EN-485-2 cuyas características se pueden consultar en la Tabla 1, y cuyo módulo de elasticidad E es 72GPa. - 73 - Con objeto de minimizar el número de elementos empleados, teniendo en cuenta que la geometría de la pieza es simétrica en el eje X como respecto del eje Y se modela sólo la mitad superior de la pieza y se aplican condiciones de contorno de simetría en la línea horizontal inferior (U2=0). Dado que la zona de estudio es el contorno de la grieta, se restringen los movimientos de la línea vertical izquierda (U1=0) para que la pieza sea isostática. En el caso que nos ocupa, en la Fig.4.18, se puede ver una figura con la malla creada y un detalle en la Fig.4.19 de esta malla alrededor de la grieta. En esta Fig.4.19 se puede ver cómo se modela una grieta sin modificar la malla tal como se explicaba en el capítulo 2.4; en el extremo de grieta se ven los 4 nodos enriquecidos con 10gdl adicionales señalados en rojo donde se evalúa su propagación, en los alrededores de las caras de la grieta se ven los 8 elementos enriquecido con 4gdl adicionales marcados en azul; y ya más alejados se pueden ver los elementos enriquecidos con 2gdl en naranja. El modelo completo tiene un total de 34944 nodos y 34488 elementos. Figura 4.18: Modelo mallado. - 80 - 79º. En ese caso, el modo de apertura de la grieta está más afectado por el modo II de apertura. - El aumento del módulo de Young E del indentador produce que las líneas de fuerza de tracción-compresión originalmente en la dirección de la carga alternante σ B “fuguen” hacia el indentador, provocando el cambio de dirección de las tensiones alrededor de la grieta y provocando la disminución en el ángulo de propagación. - Se confirma la no proporcionalidad de los casos de carga a través de las simulaciones realizadas, lo que indica la necesidad de emplear un criterio válido bajo estas condiciones como es el que se ha empleado en este proyecto fin de master: el criterio de la mínima variación de la tensión tangencial ∆τ min . 5.2 Trabajos futuros En este apartado se plantean una serie de posibles estudios que pueden llevarse a cabo tomando como base los resultados y conclusiones obtenidos en este trabajo de investigación. - Evaluar el efecto de cambiar el material del indentador a otros diferentes con rigideces distintas de la de la probeta, como se ha realizado y documentado en el capítulo 4, utilizando XFEM y realizando ensayos experimentales para confirmar y correlacionar con las técnicas numéricas más avanzadas los resultados anticipados en este trabajo. - Aplicar a diferentes tipos de ensayos, como pudieran ser ensayos de flexión, el cambio de rigidez del indentador para predecir la orientación de la grieta en estos casos y generalizar o particularizar las conclusiones de este trabajo. - 81 - 6 Bibliografía [1] ABAQUS, Inc. ABAQUS v10.6, 2010. [2] Mathworks, Inc. MATLAB R3013b, 2013. [3] U. Bryggman and S. Söderberg. Contact conditions in fretting. Wear, 110:1-17, 1986. [4] D.A. Hills, D. Nowell and A. Sackfield. Mechanics of Elastic Contacts. Butterworth-Heinemann, Oxford, 1993. [5] D.A. Hills and D. Nowell. Mechanics of Fretting Fatigue, Solid mechanics and its applications, volume 30. Kluwer Academic Press, 1994. [6] A. Sackfield, A. Mugadu and D.A. Hills. The influence of an edge radius on the local stress field at edge of a complete fretting contact. International Journal of Solid and Structures, 39, 2002. [7] A. Sackfield, C.E. Truman, and D. A. Hills. The tilted punch under normal and shear load (with application to fretting test). International Journal of Mechanical Sciences, 43, 2001. [8] A. Mugadu and D.A. Hills. A generalised stress intensity approach to characterising the process zone in complete fretting contacts. International journal of Solid and Structures, 39, 2002. [9] A.E. Giannakopoulos, T.C. Lindley, and S. Suresh. Overview no. 129 – Aspects of equivalence between contact mechanics and fracture mechanics: Theoretical connections and a life-prediction methodology for fretting-fatigue. Acta Materialia, 46(9):2955-2968, 2002. [10] B.P. Conner, S. Suresh, and T.C. Lindley. Application of fracture mechanics based life prediction method for contact fatigue. International Journal of Fatigue, 26:511-520, 2004. [11] N.I Muskhelishvili. Some basic problems of the mathematical theory of elasticity. Noordhoff, Groningen, 1953. - 82 - [12] K.L. Johnson. Contact Mechanics. Cambridge University Press, 1985. [13] M. Tur, F.J. Fuenmayor, and J.J. Ródenas. Influence of bulk stress on contact conditions and stresses during fretting fatigue. Journal of Strain Analysis for Engineering Design, 37(6):479-492, 2002. [14] V.B. Watwood, Jr. The finite element method for prediction of crack behaviour. Nuclear Engineering and Design, 11:323-332, 1969. [15] J.R. Dixon and L.P. Pook. Stress intensity factors calculated generally by the finite element technique. Nature, 224:166-167, 1969. [16] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46(1):131-150, 1999. [17] J.R. Rice and D.M. Tracey. Computational fracture mechanics. In S.J. Fenves, N. Perrone, A.R. Robinson, and W.C. Schnobrich, editors, Numerical and Computer Methods in Structural Mechanics, pages 585623, New York, 1973. Academic Press. [18] R.H. Gallagher. A review of finite element techniques in fracture mechanics. In D.R.J. Owen and A.R. Luxmoore, editors, Numerical Methods in Fracture Mechanics, Proceedings 1st Conference, pages 1-25, Swansea, 1978. Pineridge Press. [19] D.R.J. Owen and A.J. Fawkes. Engineering Fracture Mechanics: Numerical Methods and Applications. Pineridge Press, Ltd. Swansea, UK., 1983. [20] S.N. Atluri. Computational Methods in the mechanics of Fracture, volume 2 of Computational Methods in Mechanics. North Holland (Elsevier Science), Amsterdam, 1986. [21] I.S. Raju and J.C. Newman, Jr. Methods for analysis of cracks in three-dimensional solids. Journal of Aero. Soc. Of India, 36(3):153-172, 1984. - 83 - [22] A.A. Griffith. The phenomena of rupture and flow in solids, volume 221:163-198 of Sereis A. Philosophical Transactions of the Royal Society of London, St. Louis, Missouri, 1921. [23] G.R. Irwin. Fracture Dynamics. Am. Soc. Metals, Cleveland, 1948. [24] H.M. Westergaard. Bearing Pressure and Cracks, volume 6. Journal of Applied Mechanics, 1937. [25] M.L. Williams. Stress singularities resulting from Various Boundary conditions in Angular Corners of Plates in Extension. Journal of Applied Mechanics 19:526-528, 1952. [26] E.E. Gdoutos. Fracture Mechanics: an Introduction. Solid Mechanics and its Applications. Kluwer Academic Publishers, Dordrecht, Holanda, 1993. [27] T.L. Anderson. Fracture Mechanics: Fundamentals and Applications. CRC Press, Boca Ratón, Florida, 2 nd edition, 1995. [28] G.C. Sih. Methods of Analysis and Solutions of Crack Problems, volume 1 of Mechanics of Fracture. Noordhoff International Publishing, Leyden, Netherland, 1973. [29] M.H. Aliabadi and D.P. Rooke. Numerical Fracture Mechanics. Computational Mechanics Publications: Solid Mechanics and its Applications. Kluwer Academic Publishers, United Kingdom, 1991. [30] T. Belytschko, Y.Y. Lu, and L.Gu. Element-free Galerkin Methods. International Journal for Numerical Methods in Engineering, 37, 1994. [31] E. Giner, A. Vercher, J.E. Tarancón, O.A. González and F.J. Fuenmayor. Análisis mediante X-FEM de la orientación de grieta en un problema de fretting-fatiga con contacto completo. In Anales de mecánica de la fractura, volume 23, pages 141-146. Albarracín, 2006. [32] P.J.E Forsyth. A two stage process of fatigue crack growth. In: Proc crack propagation symposium, the College of Aeronautics, vol. 1. Cranfield; 1961. p. 76-94. - 84 - [33] M.C. Dubourg and V. Lamacq. Stage II crack propagation direction determination under fretting fatigue loading: a new approach in accordance with experimental observations. In: Hoeppner D.W., et al., editors. Fretting fatigue: current technology and practices, ASTM STP 1367, West Conshohocken; 2000. p. 436-50. [34] M. Sabsabi. Modelado de grieta y estimación de vida en Fretting Fatiga mediante el Método de los Elementos Finitos Extendido X-FEM. PhD Thesis, Universidad Politécnica de Valencia, 2010. [35] E.Giner, M. Sabsabi, J.J. Ródenas, F.J. Fuenmayor. Direction of crack propagation in a complete contact fretting-fatigue problem. International journal of Fatigue 58:172-180. 2014. [36] B. Cotterell, J.R. Rice. Slightly curved or kinked cracks. Int J Fract 1980;16:155–69. [37] R. Ribeaucourt, M.C Baietto-Dubourg, A. Gravouil. A new fatigue frictional contact crack propagation model with the coupled XFEM/LATIN method. Comput Methods Appl Mech Eng 2007;196:3230– 47. [38] Giner E, Sabsabi M, Fuenmayor FJ. Calculation of KII in crack face contacts using X-FEM. Application to fretting fatigue. Eng Fract Mech 2011;78(2):428–45. [39] Giner E., Sukumar N., Denia F.D., and Fuenmayor F.J. Extended Finite Element method for fretting fatigue crack propagation. International Journal of Solid and Structures, 45:5675-5687, 2008. [40] Mutoh Y., Xu J.Q., and Kondoh K. Observation and analysis of fretting fatigue crack initiation and propagation. In S.E. Kinyon, D.W. Hoeppner, and Y.Mutoh editors, Advanced in Basic Understanding and Applications, pages 61-75, West Conshohocken, 2003. American Society for Testing and Materials ASTM STP 1425.