Estudio numérico de grietas en sólidos 3-D transversalmente isótropos

Estudio numérico de grietas en sólidos 3-D transversalmente isótropos M.P. Ariza, M. Solís, J. Domínguez Escuela Superior de Ingenieros, Universidad de Sevilla, Camino de los Descubrimientos sin, 41092 Sevilla. Tlf: 954 487293 e-mail: [email protected]. Resumen En este artículo se presentan resultados de problemas de Mecánica de la Fractura dinámica armónica en materiales transversalmente isótropos. Estos resultados han sido obtenidos mediante la aplicación de la formulación mixta del Método de los Elementos de Contorno desarrollada por los autores. Dicha formulación permite que tan solo sea necesario discretizar la superficie de la grieta y el contorno externo del problema. La formulación se valida mediante la comparación de algunos resultados con los obtenidos por otros autores. Posteriormente se resuelven otros problemas de gran interés pero más complicados para los que no existe solución previa conocida por los autores Palabras Clave: transversalmente isótropo, fractura, elementos de contorno, dinámica. Abstract Results for severa! problems of time harmonic Fracture Mechanics in transversely isotropic solids are presented in this paper. A Mixed Boundary Element Formulation developed by the authors has been used in order to obtain these results. Only the crack surface and the externa! boundary need to be discretized. Results obtained for sorne simple problems are compared with those obtained by other authors. Then, other interesting and more complex problems, for which no previous results are known by the authors, are solved. Keywords: transversely isotropic, fracture, boundary elements, dynamics l. Introducción El comportamiento transversalmente isótropo se puede considerar como un caso particular de anisotropía, en el que el material presenta comportamiento isótropo en todos los planos perpendiculares a una dirección fija, conocida como eje de simetría del material. Esta ley de comportamiento se define por cinco constantes elásticas independientes, y se puede escribir como sigue (siendo la dirección 3 el eje de simetría del material) 0"11 e,, C,2 C13 oo o 0"22 c,2 C11 C13 oo o 0"33 C13 C13 C33 oo o 0"23 oo o C44 o o CT3¡ oo oo C44 o CT12 oo oo o (C11-C12)/2 2221 Mecánica de la Fractura León, 2004 CNIM XV f Esta relación se puede reducir al caso isótropo haciendo C 11 =C 33 =)..+2µ, C 44 = µ y C,2 = C13 = A (siendo A yµ las constantes elásticas de Larné). Diversos materiales presentan este tipo de comportamiento (hielo, materiales de estructura cristalina de tipo hexagonal, algunos terrenos arcillosos, materiales compuestos o cornposites fo rm ados por fibras orientadas en una dirección y embebidas en una matriz polirnérica, etc.). Existen por tanto un buen número de aplicaciones en las que resulta de interés para los ingenieros conocer el comportamiento de estos materiales. Esto incluye el estudio de su comportamiento ante la presencia de grietas o fisuras, de cara a predecir su posible propagación o estudiar corno puede afectar su presencia al comportamiento de la pieza oel medio que se trate. Este comportamiento se puede estudiar bajo el efecto de cargas estáticas o dinámicas. Entre estas últimas se incluiría el caso de la difracción de ondas elásticas debido a la presencia de una grieta, lo cual resulta de especial interés para la detección de grietas en este tipo de materiales mediante técnicas no destructivas, corno por ejemplo ultrasonidos. En este artículo se presentan algunos resultados para problemas de interés en la Mecánica de la Fractura de sólidos transversalmente isótropos bajo cargas armónicas. Para obtener estos resultados se ha aplicado el Método de los Elementos de Contorno (MEC), que es una herramienta numérica especialmente adecuada para resolver problemas de este tipo, debido a su capacidad para reducir el problema al contorno. Esto hace que no sea necesario discretizar el dominio del problema, con la consecuente simplificación del proceso de mallado, sobre todo en problemas tridimensionales. Tampoco se hace necesario por tanto la discretización de la zona cercana al frente de grieta, donde tiene lugar un comportamiento singular del campo de tensiones. Por otro lado, se pueden resolver problemas en medios infinitos de una fo rm a natural ' satisfaciendo automáticamente la condición de radiación y sin necesidad de recurrir a contornos ficticios que cierren el problema. Estas ventajas del método quedan del todo patentes utilizando una Fo rm ulación Mixta del MEC, en la que se hace uso de la formulación hipersingular a la hora de plantear las 2222 XV f CN Mecánica de la Fractura IM león, 2004 __________________ .......:...;_;_ __ ecuaciones en la superficie de la grieta, evitando así la técnica de subdominios necesaria en la fo rm ulación clásica del método. 2. Formulación La Ecuación Integral de Contorno (EIC) en desplazamientos, para un punto y de un cuerpo Q, de contorno r, bajo carga a rm ónica se puede escribir como c 1k (y) · u 1 (y,@ ) + f p,:(x,y,m) · u k (x,@ ) · dr -fu 1 : (x,y,@) · A(x,m) · dr = O r r (1) para l,k- 1,2,3, donde Uk y Pk son la componente k del desplazamiento y las tracciones, respectivamente, ro es la frecuencia, y u* 1k , p*1k son los campos de desplazamientos y tracciones de la solución fundamental armónica, elástica y tridimensional para materiales transversalmente isótropos. El coeficiente Ctk es un coeficiente que depende de la superficie del contorno en el punto y, siendo su valor Y2 para un punto suave del contorno y 1 para un punto interno. Los detalles sobre la formulación elastodinárnica del MEC se pueden encontrar en [ 1 J. Los términos de la solución fundamental se pueden descomponer en una parte estática y otra parte dinámica, que tiende a cero con la frecuencia. Los té rm inos de la parte estática fueron obtenidos por Pan y Chou [2], mientras que la parte dinámica se toma de la solución fundamental dinámica a rm ónica obtenida por Wang y Achenbach [3] Cuando existe una grieta, la coincidencia geométrica entre las superficies de ambos labios de la grieta impide que se pueda formar a partir de la EIC en desplazamientos un número suficiente de ecuaciones linealmente independientes para resolver el problema. Esto obliga a recurrir a la división del dominio en subdominios cuyo contorno contenga a la superficie de la grieta, y después aplicar condiciones de compatibilidad entre los subdominios. Este inconveniente se puede solventar mediante la utilización de la EIC en tracciones, también conocida como fo rm ulación hipersingular del MEC. Esta ecuación se obtiene a partir de la EIC en desplazamientos, derivándola en el punto de colocación, combinando estas ecuaciones derivadas según la ley de comportamiento del material y proyectándolas según la no rm al unitaria exterior en dicho punto de colocación (N(y)), 2223 Mecánica de la Fractura León, 2004 CNIM XV ( de forma que en el término libre de la ecuación aparezcan las tracciones en este punto (entiéndanse por tracciones la proyección del tensor de tensiones). La EIC en tracciones se puede escribir como donde s*imk y d*imk son combinaciones lineales de p*1k y u*1k respectivamente, y sus expresiones se pueden encontrar en [ 4] y [5] para materiales isótropos y transversalmente isótropos respectivamente. La formulación hipersingular no plantea ninguna complicación desde el punto de vista conceptual, pero necesita un tratamiento especial desde un punto de vista numérico. Una revisión de los trabajos más relevantes se puede encontrar en [ 6] y [7]. Las dificultades numéricas en la implementación son debidas a que los núcleos s*imk y d*tmk son hipersingulares y fuertemente singulares, respectivamente. Es decir, presentan singularidades de orden r- 3 y r· 2, cuya integración numérica plantea dificultades. Para solventar esta dificultad, Ariza y Domínguez han desarrollado un proceso de regularización para materiales isótropos y transversalmente isótropos [ 4, 5], con el que se reduxen las integrales hipersingulares y fuertemente singulares a otras integrales de superficie y de línea que son a lo sumo débilmente singulares y pueden ser evaluadas numéricamente sin dificultad. La EIC en tracciones, una vez realizado el proceso de regularización se pede escribir como c1 k · p¡(y,m) + J-s, , kNm { [u k(x,m)- uk(y,m) - uk,h(y,m)(xh -yh) ]- d1 ,kNnJA (x,m) -pk(y,m) ] }df + [u k(y,m)l 1 k + uk.h(y,m)J 1hk + A (y,m)K1k ]+ J. {s 1 : k (x,y, m) · N"' (y).uk (x,m)- d1 :,1 (x,y, m) · N 1 11 (y)· A (x,m) }ir (3) donde los términos ltk , J¡hk y K1k recogen las nuevos términos integrales, y sus expresiones se pueden encontrar en los trabajos de Ariza y Domínguez [ 4,5]. Se puede observar en la ecuación (3) que el proceso de regularización se aplica tan solo a la parte estática de la solución fundamental, ya que la parte dinámica es regular y no plantea ninguna dificultad en cuanto a su integración numérica. Además, por hacer más simple la implementación, la regularización se aplica tan solo al elemento que contiene al punto de colocación. 2224 --- XVI Mecánica de la Fractura CNIM León,2004 ----------------------- Por otro lado, la derivación de la EIC en desplazamientos en el punto de colocación requiere que se cumpla continuidad C 1 en dicho punto. Este requerimiento se cumple mediante el Método de Multicolocación propuesto por Gallego y Domínguez en 2-D [8] y extendido al caso tridimensional por Ariza y Domínguez [5]. Esta técnica consiste en ubicar los puntos de colocación en posiciones cercanas a cada nodo pero dentro de un elemento, de forma que se cumpla el requerimiento de continuidad. Esto hace que en el elemento contiguo se tengan que evaluar integrales que son cuasi-singulares de orden r· 3. Para evaluar estas integrales de manera eficiente se aplica la transformación de coordenadas propuesta por Telles [9], cuyo efecto quedó reflejado en los resultados presentados en [10]. La formulación mixta del MEC desarrollada por Ariza y Domínguez en materiales isótropos y transversalmente isótropos consiste en formular la ecuación integral en desplazamientos en los puntos de colocación pertenecientes al contorno externo, cuando éste existe, mientras que en la superficie de la grieta caben dos opciones: discretizar tan solo una de las superficies y formular en ella la EIC en tracciones, o discretizar ambas superficies y formular en una de ellas la EIC en tracciones y en la otra la EIC en desplazamientos. Con la primera opción se tienen como variables en los nodos de la grieta los Desplazamientos de Apertura de Grieta (DAG), que son suficientes para calcular los FIT, que es normalmente lo que interesa obtener en este tipo de problemas. La segunda opción permite conocer los desplazamientos en uno y otro labio de la grieta. Los FIT se pueden obtener a partir de los DAG de los nodos a 14 de los elementos del mismo tipo que se ubican a lo largo del frente de grieta. Para el modo I, que es el que tiene lugar en los ejemplos numéricos que se presentan, se utiliza la siguiente expresión K = ¡-;¡ ful /Ji l u = _3_ _ _3_ l+1ni 1+ donde u z es el DAG en el nodo a Vi, L el tamaño del elemento a Vi en sentido normal al frente de grieta, y 1 , m 1 y m 2 son factores que dependen de las constantes elásticas del material. 2225 Mecánica de la Fractura 3. Resultados numéricos León, 2004 CNIM XVI A continuación se presenta la solución de dos problemas en un medio infinito en presencia de una grieta circular en el primer ejemplo y de forma elíptica en el segundo, siendo el plano de la grieta perpendicular al eje de simetría del material. En ambos casos, la grieta se encuentra bajo la acción de una carga armónica. 3.1 Grieta circular en medio infinito Se consideran dos tipos de materiales, uno isótropo (G=106 Pa, v=0.3, p=lOO kg/m 3 ) y un material compuesto grafito-epoxy (C 11= 13.92 GPa, C 33 =160.7 GPa, C 12= 6.92 GPa, C1 3= 6.44 GPa, C44=7.07 GPa, p= 1578 kglm\ Se presenta la amplitud del DAG a lo largo del radio de la grieta para tres frecuencias distintas, y se comparan los resultados con los obtenidos por Kundu y Bostrom [11]. En la figura 1 se representan las dicretizaciones empleadas en la resolución. La malla más refinada se ha empleado para resolver el problema con las dos frecuencias más altas. (a) (b) Figura l. Discretizaciones de grieta circular en medio infinito En las figuras 2 a, b y c se representan los resultados obtenidos para cada una de las frecuencias en el material isótropo y en las figuras 3 a, b y c en el caso transversalmente isótropo. Las frecuencias consideradas son aquellas que hacen que el producto k r *a adopte el valor 2, 10 y 20, siendo ky = ro/(p/C44) 112 el número de onda asociado a cada frecuencia y a el radio de la grieta. Los resultados se encuentran adimensionalizados con respecto al DAG en el centro de la grieta en el caso de carga estática de la misma amplitud. En las figuras se puede comprobar como la superficie de la grieta se deforma de forma más compleja conforme aumenta la frecuencia. 2226 XVI CNIM León, 2004 Mecánica de la Fractura Se puede observar que los resultados están de acuerdo con los obtenidos por Kundu y Bostrom, salvo una ligera discrepancia en el caso de alta frecuencia. Estos mismos autores advierten sin embargo en [ 11] que para valores tan altos de la frecuencia la aproximación que utilizan necesitaría un mayor número de términos. 0.8 0.6 0.4 02 -"--- .. ····· KlXldu & Bostrom ',, 1 """· ', 1 ··---Ariza et al. '\: \1 0.4 . · 0.35 J 0.3 i º 0.15 0.1 0.05 O+-----------------' 0+--- --- ------- 0.2 0.4 0.6 0.8 0 0.1 02 0.3 0.4 rli.5 0.6 0.7 0.8 0.9 ''ª º·3 r ...... __ m_ ........... _. .. ·-- -·-··-- Material isotropo Kt*a=20 m ............. -.................... ,- .. __. .. ¡ 0.2S ' \ \· ··· Ariza et al. ( 2 \ -Kundu & Bostrom l a 0.1S \ 1 ' 0 0.1 \._ ,; ;..: '-.. ' .. ·/""' .-;,.-,,.,. °', \¿).:_;;;,/ -- , -J ·- ::::>" ,/ ''\ O.O, 0.1 02 0.3 0.4 O.S 06 0.7 O.?I O'i ,,. Figura 2. DAG para grieta circular en medio infinito isótropo bajo ondas incidentes en dirección normal al plano de la grieta y distintas frecuencias 0.45 0.4 0.35 1 0.3 o 0.25 0.2 0.15 0.1 o.os · ···· Ariza et at -k"1Clu & Bos1rom O U U M U M 1 "' OA ,..,__ __ 1.6 ! 1.4 J 1, o 0.8 0.6 0.4 º·' 02 0.35 Material graflto-epoxy Kt"a=20 0.3 -K"""" & B0&trom 025 0.2 0.15 0.1 0.05 0.2 0.4 0.6 08 -Arizaetal. 0.4 ... 0.8 ,,. Figura 3. DAG para grieta circular en medio infinito transversalmente isótropo bajo ondas incidentes en dirección normal al plano de la grieta y distintas frecuencias 2227 Mecánica de la Fractura 3.2 Grieta elíptica en medio infinito León, 2004 CNIM XV 1 Los resultados de este caso son para una grieta elíptica, con una relación entre semiejes a/b=2 inmersa en un medio isótropo y otro transversalmente isótropo, con las mismas propiedades que en el ejemplo anterior. La discretización de la grieta se muestra en la figura 4. A la hora de representar los resultados se han considerado los dos puntos más significativos: 0=0º y 0=90º, siendo 0 el ángulo definido en la figura 4a. Se representa la evolución del FIT en estos puntos para un rango de frecuencia O kr a lü y O kr b lO, según el punto de que se trate. El FIT se representa normalizado con respecto al caso estático. '" Figura 4. Discretización grieta elíptica Se puede observar como la evolución es cualitativamente similar al caso isótropo, y como para 0=90º el valor del FIT alcanza un valor más alto y presenta un mayor grado de oscilación con la frecuencia. 4. Conclusiones Se han mostrado algunos resultados de problemas de interés en la Mecánica de la Fractura de materiales transversalmente isótropos, obtenidos con una formulación numérica precisa y eficiente que permite tratar frecuencias relativamente altas. Esta formulación permitirá resolver problemas más complejos para los que no hay soluciones previas. Estos nuevos ejemplos se presentarán en nuevo artículo en el que se tratarán más en detalle algunos aspectos relativos a la implementación de la formulación empleada. -grafüo-epoxy n2 -1sotropo " 1_0 ,e o.e º' º' -graf1to-epoxy -isótropo º·º o � ��-�-��----'10 º·ºo������-�_, K,a K.,a Figura 5. FIT modo I en grieta elíptica. a) 0=0°, b) 0=90º 2228 XVI CNIM León, 2004 Mecánica de la Fractura 5. Referencias l. J. Domínguez, Boundary Elements in Dynamics, Computational Mechanics Publications, Southampton y Elsevier Applied Science, London, 1993. 2. YC Pan. and TW. Chou, Point force solution for an infinite transversely isotropic solid,J. Appl. Mech., 43, (1976) 608--612. 3. CY. Wang, and JD Achenbach, 3-D Time-harmonic elastodynamic green's functions for anisotropic solids, Proc. Roy. Soc. London, A449, (1995), 441--458. 4. J. Dominguez, M.P. Ariza, and R. Gallego, Flux and traction boundary elements without hypersingular or strongly-singular integrals, lnt. J. Numer. Meth. Engng. 48, (2000), 111--135. 5. MP. Ariza y J. Dominguez, Boundary element formulation for 3-D transversely isotropic cracked bodies, lnt. J. Numer. Meth. Engng. (In press). 6. MH. Aliabadi, Boundary element formulations in fracture mechanics, Appl. Mech. Review ASME, 50, ( 1997), 83--96. 7. DE. Beskos, Boundary element methods in dynamic analysis:Part II (1986-1996), Appl. Mech. Rev., 50, (1997), 149--197. 8. R. Gallego and J. Dominguez, Hypersingular BEM for transient elastodynamic, lnt. J. Numer. Meth. Eng., 39, (1996), 1681--1705. 9. J.C.F. Telles, A self-adaptative coordinate transformation for efficient numerical evaluation of general boundary element integrals, lnt. J. Numer. Meth. Engng., 24, (1987) 959-973 10. M. Solís, M. P. Ariza, J. Domínguez, Resolución de problemas de grietas 3-D en materiales transversalmente isótropos mediante el método de los elementos de contorno, Anales de Mecánica de la Fractura, XXI, 2004, 126-131. 11. T. Kundu and A. Bostrom, Axisymmetric scattering of aplane longitudinal wave by a circular crack in a transversely isotropic solid, ASME J. Appl. Mech., 58 (1991) 695-- 702. 5. Agradecimientos Este trabajo ha sido posible gracias a la financiación recibida por parte del Ministerio de Ciencia y Tecnología (DPI2000-12l 7-C02-0l y DPI2001-2377-C02-0l). 2229