Implementación de elementos finitos singulares especiales para grietas en interfaces adhesivas

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Ingeniería Aeroespacial Implementación de elementos finitos singulares especiales para grietas en interfaces adhesivas en modo III Autor: Manuel Romero Laborda Tutor: Vladislav Mantic Lescisin, Mar Muñoz-Reja Moreno Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Grado Ingeniería Aeroespacial Implementación de elementos finitos singulares especiales para grietas en interfaces adhesivas en modo III Autor: Manuel Romero Laborda Tutor: Vladislav Mantic Lescisin, Mar Muñoz-Reja Moreno Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Grado: Implementación de elementos finitos singulares especiales para grietas en interfaces adhesivas en modo III Autor: Manuel Romero Laborda Tutor: Vladislav Mantic Lescisin, Mar Muñoz-Reja Moreno El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Resumen Las soluciones clásicas de la mecánica de la fractura lineal (MFEL) para los campos de tensiones y desplazamientos cercanos a la punta de la grieta son bien conocidos. Para el campo de tensiones estudiado a lo largo de una línea radial r , existe una singularidad en la punta de la grieta proporcional a Kr−1/2 , donde dicha K se corresponde con el factor de intensificación de tensiones diferente según qué modo de fractura sea el empleado (KI,KII,KIII) Figura 1 Modos de carga MFEL. Será por tanto una misión fundamental de la mecánica de la fractura el cálculo numérico o analítico del factor de intensificación de tensiones para saber cuando se producirá propagación de la grieta. Debido a las numerosas hipótesis usadas en el MEF, definidas mediante diferentes funciones polinómicas en los elementos, este método no es capaz de obtener una representación exacta del comportamiento cercano a la región de la singularidad, este error se podría disminuir mediante un mayor refinamiento cerca de dicho punto, sin embargo esto es costoso numéricamente y computacionalmente, por tanto se decide no por un refinado tan detallado sino por la generación de un mallado correcto con diversas funciones de formas y elementos capaces de captar el comportamiento de la grieta. En este caso el principal problema de estudio corresponderá con el análisis del problema a fractura a modo III o problema antiplano, cuya solución corresponde directamente con la resolución de la ecuación de Laplace adaptada a la elasticidad. Como se demostrará posteriormente este problema estará adaptado a superficies con interfaces adherentes. Esta interfaz adherente será modelada mediante una distribución continua de muelles cuya condición de contorno es denominada como condición de contorno de Robin. El análisis de la misma sumada a las diferentes de contorno de Dirichlet (en desplazamientos) o Neumann (en tensiones) permitirá la resolución de la ecuación de Laplace para numerosas geometrías y problemas hasta el análisis del problema final, consistente en la aplicación de ciertas funciones de forma singulares, capaces de captar el comportamiento de estas tensiones en la punta de la grieta. I Índice Resumen I 1 Generación y obtención funciones de forma desplazamientos 1 1.1 Obtención de funciones de forma singulares 4 1.1.1 Uso de elementos definidos por funciones de forma integradas numéricamente 4 1.1.2 Generación funciones interpoladoras mediante polinomios de Lagrange 5 Set inicial conocido, inclusión N3 y posterior iteración 7 Generación del set inicial mediante polinomios de Lagrange, inclusión N3 y posterior iteración 7 Funciones de forma No singulares 9 Funciones de forma singulares 14 2 Origen y definición del problema 17 2.1 Ecuación Laplace y problema antiplano 17 2.1.1 Similitud con la ecuación de calor 18 Solución exacta 19 2.1.2 Formulación variacional para un problema antiplano con condiciones de contorno Dirichlet, Neumann y Robin 20 2.2 Programas base 21 2.2.1 Método Alberty/SH 21 Problema inicial 21 2.2.2 Evolución del programa 23 Placa rectangular con diferentes condiciones Neumann y Dirichlet 23 Inclusión de condición de Robin para placa rectangular 27 2.2.3 Solución analítica 31 (D-R) 90o32 (N-R) 90o34 (D-R) 45o36 (N-R) 45o38 (D-R) 135o39 (N-R) 135o41 2.2.4 Tensiones con elementos normales 42 Tensiones problema D-R 90º 45 Tensiones problema N-R 90o46 3 Análisis grieta mediante funciones de forma no singulares y singulares 47 3.1 Definición del problema 47 3.2 Elementos no singulares 48 3.2.1 Solución en desplazamientos 48 3.2.2 Solución en tensiones 51 3.3 Análisis mediante funciones de forma singulares 59 3.3.1 Cambios de variable y Jacobiano 59 III 6Capítulo 1. Generación y obtención funciones de forma desplazamientos La idea principal del algoritmo será la modificación por tanto del conjunto de funciones (n-m) que no verifiquen las propiedades, y tras ello producir una actualización a los nodos iniciales para mantener las condiciones activas, es decir, el seguimiento lleva a cabo 3 pasos: •Paso 1: Nm+1(r)←Nm+1(r)−∑m a=1Nm+1(ra)Na(r) Nm+1(rm+1)−∑m a=1Nm+1(ra)Na(rm+1)(1.30) •Paso 2: Na(r)←Na(r)−Na(rm+1)Nm+1(r),a=1,2,...,m(1.31) •Paso 3: Si m+1<nreemplazar m por m+1 y repetir pasos 1-3,si por el contrario m+1=nparar. Este algoritmo posee por tanto una serie de consecuencias: • Tras completar los pasos 1 y 2 las funciones de forma verifican la propiedad de interpolación en los primeros m+1 nodos, para completarlo en todos basta con mantener el número de iteraciones. • Aunque en numerosos casos como este el algoritmo puede ser calculado a mano, este posee una fácil implementación en una subrutina de creación de funciones de forma. •Las derivadas de dichas funciones pueden ser construidas de manera similar Antes de realizar la explicación de la obtención de las funciones de forma de este trabajo se decide poner un ejemplo para su mejor comprensión: Ejemplo: Considerando espacio unidimensional con 3 nodos posicionados en r1=0 , r2=1/2 , r3=1 , se desea construir el conjunto capaz de representar 1,r,rα, para ello se parte del set de funciones de forma: N1(r) = 1−2r;N2(r) = 2r;N3(r))rα(1.32) Como se puede observar Na(rb) = δab si 1≤a,b≥2 de donde se puede deducir por tanto en términos del algoritmo anterior m=2,n=3, es decir, se verifican las características con la existencia de dos nodos únicamente pero la presencia de un tercer nodo con esa forma genera la alteración de todas las funciones de forma mediante el algoritmo anterior, quedando de tal manera que: N3(r)←N3(r)−N3(0)N1(r)−N3(1/2)N2(r) N3(1)−N3(0)N1(1)−N3(1/2)N2(1)=rα−2(1/2)αr 1−2(1/2)α(1.33) N1(r)←N1(r)−N1(1)N3(r) = 1−2r+rα−2(1/2)αr 1−2(1/2)α(1.34) N2(r)←N2(r)−N2(1)N3(r) = 2r−2rα−2(1/2)αr 1−2(1/2)α(1.35) Este caso encuentra como se ha comentado un orden rα , pero la combinación de diversos nodos y diversas funciones de forma pueden conseguir la obtención de diferentes bases como (1,r,rα,rβ) con β variante con α=1,2,...,β=1,2,... . Esto de aquí genera un alto número de posibilidades de generación de las cuales se ha encontrado una solución adecuada. Dicho orden de posibilidades corresponde: •Variación en el número de nodos:2,3... En nuestro caso el número de nodos a uso será 3 • Variación en la posición de los nodos intermedios del elemento:mitad de sección, posicionamiento a un cuarto del elemento,...,etc. • Partición de una base lineal conocida sumado a la adición de una última función de forma con la forma rlnr . Todo esto iría acompañado por una posterior iteración hasta la verificación completa de las propiedades. • Obtención del set inicial de funciones de forma mediante polinomios de Lagrange y posterior iteración hasta verificación completa de las propiedades. A continuación se muestran las dos opciones con resultados coherentes y finitos obtenidas para 3 nodos. 1.1 Obtención de funciones de forma singulares 7 Set inicial conocido, inclusión N3 y posterior iteración Para la obtención y formación de este conjunto de funciones de forma se partirá del posicionamiento de los nodos en las posiciones r= [0,1 2,1], además de un set inicial de funciones de forma conocidas: N1(r) = 1−2r;N2(r) = 2r;N3(r) = rln(r)(1.36) Como se puede comprobar con facilidad m=2,n=3, debido a la presencia de 3 nodos(n) y al funcionamiento de Na(rb) = δab si 1≤a,b≥2.Esto genera un conjunto final de funciones de forma: N3(r)←N3(r)−N3(0)N1(r)−N3(1/2)N2(r) N3(1)−N3(0)N1(1)−N3(1/2)N2(1)=r[ln(1/2)−ln(r)] ln(1/2)(1.37) N1(r)←N1(r)−N1(1)N3(r) = 1−2r+r[ln(1/2)−ln(r)] ln(1/2)(1.38) N2(r)←N2(r)−N2(1)N3(r) = 2r−2r[ln(1/2)−ln(r)] ln(1/2)(1.39) y como se puede comprobar: 3 ∑ i=1 Ni(r) = 1(1.40) además de la condiciones: N1(0) = 1; N2(0) = 1; N3(0) = 0(1.41) N1(1/2) = 0; N2(1/2) = 1; N3(1/2) = 0(1.42) N1(1) = 0; N2(1) = 0; N3(1) = 1(1.43) Generación del set inicial mediante polinomios de Lagrange, inclusión N3 y posterior iteración En este apartado se utilizarán el posicionamiento de nodos r= [0,1 4,1] , la principal diferencia aquí será la obtención del set inicial de funciones. Para ello se propone generación de dos funciones de forma mediante polinomios de Lagrange, y la inclusión de Ni(r) = rln(r) como última función de forma para el nodo restante, a todo ello habría que incluir la posterior iteración. Las numerosas opciones de generación del set inicial a una previa iteración comprende: N1=l1(r);N2=l2(r);N3=rln(r)(1.44) N1=l1(r);N2=rln(r);N3=l3(r)(1.45) N1=rln(r);N2=l2(r);N3=lr(3)(1.46) Donde cabe destacar la(f) corresponde con los polinomios de Lagrange en el punto a de la función f explicado con anterioridad. Desafortunadamente ninguna de estos conjuntos de funciones iniciales van producir conjuntos finales de funciones de forma correctos tras iteración, ya sea por presencia de funciones con tendencia infinita antes o después de la iteración. A partir de aquí, se decide probar con la obtención del set inicial mediante polinomios de Lagrange con una excepción que el anterior, la función sobrante no se generará mediante Ni=rln(r) sino mediante Ni=la(rln(r)) , aplicando los polinomios de Lagrange a la función f=rln(r), obteniendo así varias posibilidades: N1=l1(r);N2=l2(r);N3=l3(rln(r)) (1.47) N1=l1(r);N2=l2(rln(r));N3=l3(r)(1.48) N1=l1(rln(r));N2=l2(r);N3=lr(3)(1.49) 8Capítulo 1. Generación y obtención funciones de forma desplazamientos Con todo ello se obtiene una solución viable correspondiente con: N1=l1(r);N2=l2(rln(r));N3=l3(r)(1.50) quedando como resultado: N1(r) = l1(r) = (r−1 4)(r−1) 1 4 (1.51) N3(r) = l3(r) = r(r−1 4) 3 4 (1.52) N2(r) = l2(rln(r)) = rln(r)(rln(r)−1) (1 4ln(1 4))2(1.53) No obstante a pesar del cumplimiento de las propiedades: 3 ∑ i=1 Ni(r) = 1(1.54) N1(0) = 1; N2(0) = 1; N3(0) = 0(1.55) N1(1/4) = 0; N2(1/4) = 1; N3(1/4) = 0(1.56) N1(1) = 0; N2(1) = 0; N3(1) = 1(1.57) Se ha desechado esta opción debido a la necesidad de un cumplimiento matemático tal y como es: limr→0(rln(r)) = 0(1.58) A partir de aquí podrían surgir ciertos problemas innecesarios los cuales se evitan usando el primer conjunto de funciones de forma obtenidas. Una vez obtenidas dichas funciones de forma en el espacio unidimensional falta adaptarla al espacio de referencia bidimensional ξ−η , lo cual se realiza mediante la generación de un elemento de referencia con 6 nodos lineal en ηy con las funciones de forma específicas adaptadas a ξ: N3(r) = r[ln(1/2)−ln(r)] ln(1/2)(1.59) N1(r) = 1−2r+r[ln(1/2)−ln(r)] ln(1/2)(1.60) N2(r) = 2r−2r[ln(1/2)−ln(r)] ln(1/2)(1.61) A continuación una vez obtenida la evolución de cada función de forma en el eje ξ , bastaría con adaptarlo al espacio bidimensional, siendo así necesaria realizar una evolución lineal de dos nodos en el eje η como se observa en las siguientes funciones de forma. 1.1 Obtención de funciones de forma singulares 9 Quedando por tanto el siguiente conjunto de funciones de forma: N1(ξ,η) = (1−η)1−ξ−ξln(ξ) ln(1/2) N2(ξ,η) = (1−η)2ξln(ξ) ln(1/2) N3(ξ,η) = (1−η)ξ−ξln(ξ) ln(1/2) N4(ξ,η) = ηξ−ξln(ξ) ln(1/2) N5(ξ,η) = η2ξln(ξ) ln(1/2) N6(ξ,η) = η1−ξ−ξln(ξ) ln(1/2) (1.62) Junto a sus características mostradas en la tabla a continuación (1.1). Tabla 1.1 Funciones de forma en desplazamientos. Niu ∂Ni ∂ξ ∂Ni ∂η i (1−η)1−ξ−ξln(ξ) ln(1/2)(1−η)−1−1 ln(1/2)−ln(ξ) ln(1/2)−1−ξ−ξln(ξ) ln(1/2)1 (1−η)2ξln(ξ) ln(1/2)2(1−η)ln(ξ) ln(1/2)+1 ln(1/2)−2ξln(ξ) ln(1/2)2 (1−η)ξ−ξln(ξ) ln(1/2)(1−η)1−ln(ξ) ln(1/2)−1 ln(1/2)−ξ−ξln(ξ) ln(1/2)3 ηξ−ξln(ξ) ln(1/2)η1−ln(ξ) ln(1/2)−1 ln(1/2) ξ−ξln(ξ) ln(1/2)4 η2ξln(ξ) ln(1/2)2ηln(ξ) ln(1/2)+1 ln(1/2) 2ξln(ξ) ln(1/2)5 η1−ξ−ξln(ξ) ln(1/2)η−1−1 ln(1/2)−ln(ξ) ln(1/2)η1−ξ−ξln(ξ) ln(1/2)6 Como se mostrará en los apartados más adelante, estas funciones de forma serán denominadas en desplazamientos ya que tratan de modelar tanto desplazamientos y tensiones, serán fundamentales para la obtención de dichos en la matriz de rigidez, mientras que para el modelado de la matriz Jacobiana y la geometría se usarán funciones de forma geométricas más sencillas. Debido al cambio de espacio entre espacio real x−y y espacio de referencia ξ−η , se considera oportuno por tanto la evolución de dichas funciones de forma tanto en el espacio de referencia como en el espacio real. Para ello se ha decidido representar funciones de forma cuadráticas en ambos espacios y posteriormente nuestras funciones de forma con motivo de observar la diferencia de comportamiento. Cabe destacar que el espacio real x-y dependerá del mallado en cada punto del espacio x-y del problema correspondiente por tanto se representará sobre un triangulo del espacio x-y generalizado, por tanto para analizar cómo evolucionan estas funciones en los diferentes espacios real y referencia observado en las figuras 1.6, 1.7, como se ha comentado, se tratará primeramente unas funciones de formas típicas y posteriormente las obtenidas por el procedimiento anterior. Funciones de forma No singulares Para este apartado, se utilizarán las ideas del artículo [ 10 ], con lo que se analizará un elemento de referencia rectangular de 8 nodos y su evolución al elemento real triangular general de 8 nodos también de los cuales 3 estarán colapsados, por tanto los espacios a tratar serán los mostrados en la figura 1.6, para ello al ser un elemento de 8 nodos tendrá 8 funciones de forma que se encuentran generalizadas para la siguiente fórmula obtenida en [10]: Ni(ξ,η) = [(1+ξξi)(1+ηηi)−(1−ξ2)(1+ηηi)−(1−η2)(1+ξξi)]ξ2 iη2 i/4+(1.63) (1−ξ2)(1+ηηi)(1−ξ2 i)η2 i/2+(1−η2)(1+ξξi)(1−η2 i)ξ2 i/2 10 Capítulo 1. Generación y obtención funciones de forma desplazamientos Figura 1.6 Espacios real y referencia. Figura 1.7 Espacios real y referencia para elemento rectangular de 8 nodos. Para distinguir entre las funciones de forma se tiene: N1(ξ,η)→ξi=−1;ηi=−1(1.64) N2(ξ,η)→ξi=0;ηi=−1(1.65) N3(ξ,η)→ξi=1;ηi=−1(1.66) N4(ξ,η)→ξi=1;ηi=0(1.67) N5(ξ,η)→ξi=1;ηi=1(1.68) N6(ξ,η)→ξi=0;ηi=1(1.69) 1.1 Obtención de funciones de forma singulares 11 N7(ξ,η)→ξi=−1;ηi=1(1.70) N8(ξ,η)→ξi=−1;ηi=0(1.71) A partir de aquí se observará la evolución de cada una de las funciones de forma en el espacio de referencia y en el espacio real. Es importante destacar la naturaleza de la punta del triángulo en el espacio real, cuyo resultado es un colapso de las funciones de forma 7,8,1 tal que: (N1(ξ,η)+N8(ξ,η)+ N7(ξ,η))espacio ξ−η=N1espacio real(x,y)(1.72) Tal y como se puede ver a continuación por la imagen del colapso además de los resultados numéricos: Figura 1.8 Colapso funciones de forma para formar espacio triangular,sacada de [10]. Figura 1.9 N2yN2en espacio triangular. 12 Capítulo 1. Generación y obtención funciones de forma desplazamientos Figura 1.10 N3yN3en espacio triangular. Figura 1.11 N4yN4en espacio triangular. Figura 1.12 N5yN5en espacio triangular. 1.1 Obtención de funciones de forma singulares 13 Figura 1.13 N6yN6en espacio triangular. Figura 1.14 N7yN8. Figura 1.15 N1yN7+N8+N1=N1tEn espacio triangular. 14 Capítulo 1. Generación y obtención funciones de forma desplazamientos Funciones de forma singulares En esta sección se realizará el mismo análisis para las funciones obtenidas en la sección anterior. En este caso el espacio de referencia desde donde parte la generación de funciones será distinto al anterior, y la representación del espacio x-y (que en cada problema dependerá del mallado en cuestión) será otro espacio elegido generalizado, quedando: Figura 1.16 Espacios real y referencia. Estas funciones de forma representadas como evolución en ξquedarían tal que: Figura 1.17 Evolución en ξ. Se observa además como la función N2 supera la unidad en 2 ocasiones, esto no afectará nuestro proyecto ya que el conjunto entero de funciones de forma cumplen las diversas características mencionadas anteriormente, e individualmente, sigue siendo valor unidad en el mismo nodo. A partir de aquí se realiza un análisis semejante al anterior pero con el nuevo set de funciones de forma. Es importante observar en el nuevo conjunto de gráficas la evolución logarítmica en el eje ξ mientras que mantiene una evolución lineal en el eje η: 1.1 Obtención de funciones de forma singulares 15 Figura 1.18 N2yN2En espacio triangular. Figura 1.19 N3yN3En espacio triangular. Figura 1.20 N4yN4En espacio triangular. 22 Capítulo 2. Origen y definición del problema Para su resolución por tanto se realiza el mallado mostrado en la anterior imagen. Para un análisis más detallado se ha realizado diferentes análisis según los términos de la ecuación los cuales se corresponden con: •Término de fuerzas de volumen ( f) • Término de fuerzas en el contorno ( g ) aplicado en este caso a los bordes formado por los siguientes nodos : –[1→3] –[5→7] Ambos bordes se denominarán de ahora en adelante como bordes con condiciones de contorno de Neumann (ΓN) •Contorno con desplazamiento nulo o contorno de Dirichlet (ΓD) Solución para f=1,g=0 :Para este primer caso se estudia como evoluciona el dominio frente a fuerzas de volumen en el mismo y fuerzas nulas en ambos bordes de Neumann. Además se compara el resultado de Alberty y Sanz-Herrera por lo que se puede comprobar su similitud. Observando ambas imágenes se confirma como el desplazamiento para los bordes Dirichlet es nulo frente al desplazamiento de resto de nodos producido por fuerzas en el volumen, además de tener el máximo desplazamiento en el centro al estar lo más alejado posible de dichos bordes. Figura 2.5 Solución de Alberty y Solución de SH. Se puede observar así perfectamente la solución para el problema planteado con anterioridad. Como ya se conoce en todos los bordes con condición de contorno de Dirichlet se produce un desplazamiento nulo, mientras que por el contrario todos los nodos sometidos a fuerzas de volumen sufrirán un desplazamiento en dirección de estas fuerzas de volumen, en este caso se aprecia como los nodos centrales al estar sometido a la acción de las fuerzas de volumen y estar separados lo máximo posible del centro sufren los máximos desplazamientos uz(x,y) . Además de esto los dos bordes comentados anteriormente no poseen condiciones Dirichlet sino que están sometidos a condiciones Neumann produciendo así el desplazamiento de los mismos. A partir de ahora se tratarán por tanto las soluciones de Alberty por comodidad. 2.2 Programas base 23 Solución para f=0,g=1 :En este caso no existirán fuerzas en el volumen y solo en el contorno de Neumann quedando tal que: Figura 2.6 Solución Alberty para f=0,g=1. Solución para f=1,g=1 :En este caso no existirán fuerzas en el volumen y solo en el contorno de Neumann quedando tal que: Figura 2.7 Solución Alberty para f=1,g=1. En este último caso se observa como además de las fuerzas de volumen en todo el dominio existe un desplazamiento mayor para los bordes de Neumann al estar sometidos además al efecto de las propias fuerzas en el contorno, además de aumentar también el desplazamiento global. Cualquier combinación posible y/o similar tendrá la misma distribución de desplazamientos pero con mayor o menor módulo de desplazamientos, por ejemplo una distribución f=8,g=0, la cual aumentará en 8 los desplazamientos de la figura 2.6 2.2.2 Evolución del programa A lo largo del trabajo, el programa base de Alberty se ha ido desarrollando para diferentes recintos, diferentes mallados y distintos tipos de problemas que se analizarán a continuación. Placa rectangular con diferentes condiciones Neumann y Dirichlet Para comprobación del correcto funcionamiento del programa se ha dispuesto diferentes configuraciones en una placa rectangular. Dichas configuraciones corresponden diferentes fuerzas (volumen o contorno), así 24 Capítulo 2. Origen y definición del problema como diferentes configuraciones en las condiciones de contorno. Placa rectangular N-D: Donde N-D representan condiciones Neumann-Dirichlet en un contorno únicamente por condición como se muestra en las siguientes imágenes, dentro de las cuales se tratarán un par de casos: •Fuerzas aplicadas en el contorno cuya solución conocida es lineal(g=0) •Fuerzas de volumen aplicadas en el dominio( f=1) f=0,g=1: Con solución lineal conocida: Figura 2.8 Descripción del problema. Figura 2.9 Solución del problema. 2.2 Programas base 25 f=1,g=0: En el que actúan las fuerzas de volumen Xpor lo que la solución no será lineal: Figura 2.10 Descripción del problema. Figura 2.11 Solución del problema. Placa rectangular N-3-D Denominada así debido a la presencia de una condición de Neumann y 3 bordes del contorno con condición Dirichlet. En este caso distinguiremos también entre la aplicación de fuerzas de volumen o fuerzas sobre el contorno: 26 Capítulo 2. Origen y definición del problema f=0,g=1: Aplicación sobre el contorno Neumann únicamente tal que: Figura 2.12 Descripción del problema. Figura 2.13 Solución del problema. 2.2 Programas base 27 f=1,g=0: Aplicación de fuerzas de volumen sobre el dominio tal que: Figura 2.14 Descripción del problema. Figura 2.15 Solución del problema. Inclusión de condición de Robin para placa rectangular A partir de aquí se incluirá la condición de Robin en el borde ΓR para observar su evolución y comportamiento en el programa. Demostración término de Robin: Partiendo de la formulación variacional de la ecuación 2.32 y utilizando únicamente el término referido a la matriz de rigidez aportada por la parte de Robin: Kmn =ZΓr κuvdΓ(2.43) 28 Capítulo 2. Origen y definición del problema Donde se encuentran representados los desplazamientos reales y virtuales u,v . Como consecuencia y demostrado anteriormente, la solución en el problema presentado en la ecuación 2.29 y 2.32 el siguiente término derivará en: Kmn =ZΓr κψnψmvdΓ(2.44) Nuestro conocimiento comprende de: ψm=ψm(ξ,η)(2.45) Por tanto esta integral en la longitud del dominio real pasará por tanto al espacio de referencia, sumado a esto añadimos el factor del término de Robin aplicado únicamente al lado Robin el cual denominaremos lado ξ, todo esto deriva en la siguiente formulación: Zl 0 f dl =Z1 0 fdl dξdξ(2.46) Mediante a la definición de dl: dl =dxcos(α)+dysin(α)(2.47) Teniendo en cuenta el valor de ambos tal que: cos(α) = x2−x1 longitud (2.48) sin(α) = y2−y1 longitud (2.49) Por tanto se queda tal que: Z1 0 fdx dξcos(α)+ dy dξsin(α)dξ=Z1 0 fdx dξ(x2−x1) longitud +dy dξ(y2−y1) longitud dξ(2.50) Conociendo por tanto: x=∑Nigxi(2.51) y=∑Nigyi(2.52) Y aprovechando las funciones de forma obtenidas en el artículo [4] se obtiene: dx dξ=x2−x1(2.53) dy dξ=y2−y1(2.54) Por tanto el resultado final queda de la forma: Z1 0 fdx dξ(x2−x1) longitud +dy dξ(y2−y1) longitud dξ=Z1 0 f(x2−x1)2 longitud +(y2−y1)2 longitud dξ(2.55) Siendo: (x2−x1)2+(y2−y1)2= (longitud)2(2.56) Y quedando finalmente: Zl 0 f dl =longitud ·Z1 0 f dξ(2.57) Además de la representación de Robin en una placa rectangular se ha visto oportuno el uso de la condición de Robin en la placa del problema inicial de Alberty/SH para verificar su veracidad, para ello se ha aplicado al lado: [1→2,2→3](2.58) Se puede observar por tanto como la presencia de muelles en dicho lado incrementa el matriz de rigidez en 2.2 Programas base 29 Figura 2.16 Problema original κ=0y Lado Robin con κ=20. Figura 2.17 Lado Robin con κ=100. dichos nodos, produciendo un endurecimiento de la zona y por tanto una menor cantidad de desplazamiento. El resto del dominio sin embargo se encuentra sometido a la misma fuerza de volumen. En la siguiente imagen 2.17 se resuelve el mismo problema con una κ=100 , para mostrar con más detenimiento el lado Robin se representa además en el siguiente problema en una placa rectangular como la vista con anterioridad con las siguientes condiciones de contorno: Figura 2.18 Condiciones contorno problema. 30 Capítulo 2. Origen y definición del problema Figura 2.19 Problema original κ=0y Lado Robin con κ=10. Figura 2.20 Problema original κ=0y Lado Robin con κ=10. Se ha decidido además representar para un número de alto de nodos para poder apreciar de manera más exacta el comportamiento del borde Robin, para ello se ha representado 3600 en la imagen de la izquierda en una representación normal y unos 10000 nodos aproximados en la segunda imagen representada mediante una interpolación en el mapa de colores de Matlab para no saturar la imagen de nodos: Figura 2.21 Problema original κ=0y Lado Robin con κ=10. Debido al interés se ha decidido representar los tiempos de computación según el número global de nodos, variable de interés debido al aumento de nodos, se producen soluciones más aproximadas, por tanto el aumento del número de nodos excesivamente puede no ser eficiente produciendo una mejora muy pequeña en los resultados pero un aumento exponencial del tiempo de computación. Para ello se ha hecho un ajuste polinómico con los diferentes tiempos de computación, se han aproximado por un polinomio lineal (n=1) y 2.2 Programas base 31 otro cuadrático (n=2), y se ha representado tanto normal como en escala logarítmica para el eje x. Se aprecia como un aumento de nodos produce un aumento del tiempo de computación en segundos a partir de ciertos valores. Figura 2.22 Tiempo de computación (s) para n=1. Figura 2.23 Tiempo de computación (s) para n=2. 2.2.3 Solución analítica A partir del artículo [ 8 ] en el cual se desarrolla un nuevo procedimiento analítico para la deducción de la serie asintótica de las soluciones singulares en desplazamientos y tensiones cerca del vértice de las esquinas isótropas elásticas lineales proponiendo solución a dos tipos de problemas: •Dirichlet – Robin (D-R) •Neumann – Robin (N-R) Bajo el supuesto de carga de corte antiplano, el problema elástico correspondiente se reduce a la ecuación de Laplace para el desplazamiento fuera del plano. En la deducción de tal solución singular, la variable compleja se utiliza para proponer un función armónica en forma de serie asintótica incluyendo tanto el poder como los términos logarítmicos. En ambos problemas se encontrará presente la condición de Robin la cual puede ser descrita mediante: σθz(r,θ2)+Ku(r,θ2) = 0(2.59) σrz(r,θ) = G∂u(r,θ) ∂r(2.60) σrz(r,θ) = G r ∂u(r,θ) ∂θ (2.61) 38 Capítulo 2. Origen y definición del problema (N-R) 45o La solución analítica se corresponde con: uj(r,θ) = u(0) j(r,θ)+u(1) j(r,θ)+u(2) j(r,θ)(2.82) u(0) j(r,θ) = r4j−4cos((4j−4)θ)(2.83) u(1) j(r,θ) = b(0) j,1r4j−3cos((4j−3)θ)(2.84) u(2) j(r,θ) = b(0) j,2r4j−2cos((4j−2)θ)(2.85) Con: b(0) j,1=γcos((4j−4)π) (4j−3)sin((4j−3)π/4)(2.86) b(0) j,2=b(0) j,1γcos((4j−3)π/4) (4j−2)sin((4j−2)π/2)(2.87) La comparación entre ambas soluciones queda por tanto: Aquí será también de interés la representación de Figura 2.38 Solución analítica (N-R 45 o ) u(r,θ) con γ=1,j=1 y Solución numérica (N-R 45 o ) u(r,θ) con γ=1,j=1. los lados θ=0o,θ=45opara observar la diferencia de ambas, quedando como resultado: Figura 2.39 u(r,θ=45o)con γ=1,j=1yu(r,θ=0o)con γ=1,j=1. 2.2 Programas base 39 Figura 2.40 Error total (N-R 45o)u(r,θ)con γ=1,j=1y Error total (N-R 45o)u(r,θ)con γ=1,j=1. Figura 2.41 Tiempo computación (s) n=1 y Tiempo computación (s) n=3. (D-R) 135o Debido al rango de ω∈[0o,360o] es importante analizar además soluciones mayores de 90 o , para ello se ha decidido analizar ambos problemas de 135o, cuya solución analítica para γ=1,j=1se corresponde con: uj(r,θ) = u(0) j(r,θ)+u(1) j(r,θ)+u(2) j(r,θ)(2.88) u(0) j(r,θ) = r2/3sin((2/3)θ)(2.89) u(1) j(r,θ) = (3/5)21/2r5/3sin((5/3)θ)(2.90) u(2) j(r,θ) = (9/40)r8/3sin((8/3)θ)(2.91) Gráficamente esta solución se corresponde con: Además de la representación de la línea θ=135o 40 Capítulo 2. Origen y definición del problema Figura 2.42 Solución analítica (D-R 135 o ) u(r,θ) con γ=1,j=1 y Solución numérica (D-R 135 o ) u(r,θ) con γ=1,j=1. Figura 2.43 Representación uz(r,θ=135o)Solución analítica frente numérica. En cuanto al error y tiempo de computación obtenidos: Figura 2.44 Error total (D-R 135 o ) u(r,θ) con γ=1,j=1 y Error total (D-R 135 o ) u(r,θ) con γ=1,j=1 . 2.2 Programas base 41 Figura 2.45 Tiempo computación (s) n=1 y Tiempo computación (s) n=3. (N-R) 135o Como última comprobación de solución analítica se tiene el problema de Neumann-Robin en para 135 o para γ=1,j=1en el que la solución se corresponde con: uj(r,θ) = u(0) j(r,θ)+u(1) j(r,θ)+u(2) j(r,θ)(2.92) u(0) j(r,θ) = 1(2.93) u(1) j(r,θ) = 21/2rcos(θ)(2.94) u(2) j(r,θ) = 0.5r2cos(2θ)(2.95) Figura 2.46 Solución analítica (N-R 135 o ) u(r,θ) con γ=1,j=1 y Solución numérica (N-R 135 o ) u(r,θ) con γ=1,j=1. En este problema nos vuelve a interesar los dos extremos que se corresponden con θ=0o,θ=135o quedando finalmente: 42 Capítulo 2. Origen y definición del problema Figura 2.47 u(r,θ=45o)con γ=1,j=1yu(r,θ=0o)con γ=1,j=1. Figura 2.48 Error total (N-R 135 o ) u(r,θ) con γ=1,j=1 y Error total (N-R 135 o ) u(r,θ) con γ=1,j=1 . Figura 2.49 Tiempo computación (s) n=1 y Tiempo computación (s) n=3. 2.2.4 Tensiones con elementos normales Una vez obtenidos los desplazamientos, una de las variables con alta importancia, es el cálculo de las tensiones. Para la obtención de las mismas: σxz σyz=G ∂u ∂x ∂u ∂y!=G[J−1] ∂u ∂ξ ∂u ∂η !(2.96) 2.2 Programas base 43 Donde: [J] = ∂x ∂ξ ∂y ∂ξ ∂x ∂η ∂y ∂η !(2.97) Además se conoce: ∂u ∂ξ = 3 ∑ i=1∂Ni(ξ,η) ∂ξ ui(2.98) ∂u ∂η = 3 ∑ i=1∂Ni(ξ,η) ∂η ui(2.99) Las funciones de forma empleadas hasta aquí por tanto se corresponden a las funciones de forma empleadas por Alberty y SH en el montaje de sus programas y definidas además en el artículo [4] donde se conoce: x=xj+(xj 2−xj 1)ξ+(xj 3−xj 1)η(2.100) y=yj+(yj 2−yj 1)ξ+(yj 3−yj 1)η(2.101) Quedando como resultado: [J] = (x2−x1) (y2−y1) (x3−x1) (y3−y1)(2.102) [J−1] = (y3−y1) (y1−y2) (x1−x3) (x2−x1)1 det[J](2.103) Además de: ∂u ∂ξ =u2−u1(2.104) ∂u ∂η =u3−u1(2.105) Cabe destacar como las tensiones dentro de cada uno de los elementos no singulares será por tanto constante debido a la forma lineal de cada una de sus funciones de forma del elemento. A continuación se muestran resultados para las diferentes tensiones en diferentes problemas: Figura 2.50 Problema lineal solución en desplazamientos. 44 Capítulo 2. Origen y definición del problema Y como cabe de esperar en dicho problema lineal la solución en tensiones será: Figura 2.51 σxz yσyz. Como cabe de esperar en los problemas lineales la solución para tensiones será constante. En este problema al existir únicamente tensión σxz esta será constante de valor unidad mientras que la tensión σyz será nula. Como siguiente problema se tendrá: Figura 2.52 Problema 2 solución en desplazamientos. Y como cabe de esperar en dicho problema lineal la solución en tensiones será: Figura 2.53 σxz yσyz. 2.2 Programas base 45 Por último se analizará: Figura 2.54 Problema 3 solución en desplazamientos. Figura 2.55 σxz yσyz. A partir de aquí será interesante también conocer las tensiones obtenidas en los problemas relacionados con la solución analítica, para ello se analizarán los dos casos correspondientes a 90º: Tensiones problema D-R 90º Es importante recalcar la existencia de 4 tipos de tensiones: σxz,σyz,σrz,σθz Tensiones σxz,σyz: Figura 2.56 σxz yσyz. 46 Capítulo 2. Origen y definición del problema Tensiones σrz,σθz: Para estas tensiones se ha recurrido a la solución analítica obtenida anteriormente, por lo que lo mostrado a continuación es la distribución de tensiones en (r,θ) analíticamente: Figura 2.57 σrz yσθz. Tensiones problema N-R 90o Tensiones σxz,σyz : Tensiones σrz,σθz: Para estas tensiones se ha recurrido a la solución analítica obtenida Figura 2.58 σxz yσyz. anteriormente, por lo que lo mostrado a continuación es la distribución de tensiones en ( r,θ ) analíticamente: Figura 2.59 σrz yσθz. 3 Análisis grieta mediante funciones de forma no singulares y singulares 3.1 Definición del problema Para el estudio del problema se analizará por tanto una grieta en interfaz adhesiva modelada como distribución continua de muelles. Como se ha ido viendo a lo largo de las diferentes secciones del trabajo, dicha interfaz adhesiva necesitará de su modelo teórico, en este caso, este comprenderá del borde con condiciones de contorno tipo Robin, deducida y analizada con mayor detalle en apartados anteriores. Este modelo por tanto, comprenderá de una geometría rectangular, con una distribución de tensiones σyz procedentes del modo III de fractura, y una distribución de muelles, procedentes del modelo de interfaz adhesiva tal como la siguiente: Figura 3.1 Geometría inicial. •c: Altura de la placa •b: Anchura de la placa •a: Longitud de la grieta respecto al ancho de la placa •d: Longitud de la interfaz adhesiva 47 54 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.14 σyzpara k=0.1 y σxz para k=1. Figura 3.15 σyzpara k=10 y σxz para k=100. Una de las misiones más importantes es la capacidad de seguir la condición de: σyz +κ·uz(x,y) = 0in ΓR(3.1) A continuación, por tanto, se evaluarán junto a los desplazamientos uz(x,y): 3.2 Elementos no singulares 55 Figura 3.16 σyzpara k=0.1 y σxz para k=1 frente κ·uz(x,y=0). Figura 3.17 σyzpara k=10 y σxz para k=100 frente κ·uz(x,y=0). Para una mejor apreciación se reproducen las mismas gráficas para un número menor de elementos, lo cual permitirá visualizar los límites de cada elemento marcados con marcadores rojos, además de con más resolución: Figura 3.18 σyzpara k=0.1 y σxz para k=1 frente κ·uz(x,y=0). 56 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.19 σyzpara k=10 y σxz para k=100 frente κ·uz(x,y=0). Se observa, que la máxima tensión no es alcanzada en el elemento final anterior a la grieta, en este caso posicionado en x=1, sino que el pico de tensión se encuentra varios elementos desplazados a la derecha. Como se ha comentado con anterioridad esto se debe a la presencia de elementos no singulares con tensión constante interna dentro de cada elemento. La misión del elemento singular será por tanto encontrar la discontinuidad y el pico de tensiones el en elemento situado en x=1, es decir, aquel situado en la punta de la grieta. También se observa el cumplimiento de la ecuación 3.1 a gran escala. Como se ha comentado con anterioridad, y término de gran relevancia, la tensión en la punta de la grieta no se corresponde con la máxima tensión, siendo esta la presente en elementos diferentes del posicionado para x=1. Este error será corregido con los elementos singulares, donde la máxima tensión sí estará en este elemento que posee variación interna creciente debido a su naturaleza. Esto además mostrará la discontinuidad de tensiones presente en la punta de la grieta. Por ello las gráficas a continuación representan el valor de tensión máximo que se alcanza en el dominio, distinto de aquel presente en la punta de la grieta. En ambas gráficas se muestra la evolución de la tensión máxima del recinto en función del número de nodos Figura 3.20 σyzMAX para k=0.1 frente a número global de nodos,escalas lineal y semilogarítmica. para escalas lineal y logarítmica y k=0.1. Se puede observar una evolución lenta a pesar del gran incremento de nodos que satura completamente el tiempo de computación llegando hasta un total de N=25000 nodos globales. Se puede observar una evolución similar para el resto de valores de rigidez de distribución de muelles: 3.2 Elementos no singulares 57 Figura 3.21 σyzMAX para k=1 frente a número global de nodos,escalas lineal y semilogarítmica. Figura 3.22 σyzMAX para k=10 frente a número global de nodos,escalas lineal y semilogarítmica. Figura 3.23 σyzMAX para k=0.1 frente a número global de nodos,escalas lineal y semilogarítmica. No obstante, estos análisis están realizados frente al número de nodos globales, teniendo la misma distribución de nodos tanto en x como en y. Por tanto, se decide hacer una comparativa junto al número global de nodos, teniendo en cuenta una distribución de 50 nodos en el eje y(eje de menor importancia) y una distribución incremental de nodos en el eje x de mayor importancia por la presencia de la grieta. Siguiendo esta distribución mencionada con anterioridad se observa una mayor tendencia a un número constante en menor espacio de nodos globales tal que: 58 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.24 σyzMAX para k=100 frente a número global de nodos, nodos en y=50. Se observa por tanto con este mallado constante en y e incremental en x una tendencia hacia un valor constante, sin embargo y como se ha visto, el aumento del número de nodos en y mejora la solución pero no de manera muy convergente, por ello se muestra la evolución con el número de nodos para diferentes k y 2 mallados diferentes en y: Figura 3.25 σyzMAX para varias k, frente a número global de nodos, nodos en y=10,20. 3.3 Análisis mediante funciones de forma singulares 59 3.3 Análisis mediante funciones de forma singulares Por último se procede a resolver el anterior problema mediante el uso de las funciones de forma definidas en la sección [1], con el objetivo de observar su comportamiento en una grieta tal y como se ha ido comentando en los diferentes apartados de este estudio. Para ello se parte de las funciones de forma presentes en la tabla (1.1) junto a sus respectivas derivadas parciales respecto ξyη. Es importante destacar que todos los procedimientos consultados sobre obtención de matrices de rigidez emplean funciones de forma no singulares, sin embargo al introducir estos con el objetivo de captar el comportamiento de los desplazamientos va a ser necesario distinguir entre dos tipos de funciones de forma: 1. Funciones de forma "para desplazamientos": Denominadas como Niu debido a ir asociada en la ecuación de Laplace solo a los términos que procedan de: •Desplazamientos u •Desplazamientos virtuales v 2. Funciones de forma "geométricas" asociadas a la geometría: Denominadas Nig y relacionadas tanto con el Jacobiano como con los términos x,y Al estar integrando en el espacio real (x,y) y tener conocimiento de las funciones de forma en el espacio de referencia (ξ,η) se necesitará por tanto recalcular las integrales en este último y obtener así la solución analítica para los diferentes términos de la ecuación de Laplace en los que está tanto la matriz de rigidez como el término de Robin procedente de la formulación variacional. 3.3.1 Cambios de variable y Jacobiano Para el posterior desarrollo de la matriz de rigidez será necesario definir tanto los cambios de variable como la matriz Jacobiana y su evolución: ∂u ∂ξ ∂u ∂η != [J] ∂u ∂x ∂u ∂y!(3.2) Por tanto y también aplicado a la geometría: dx dy= [J]dξ dη(3.3) La integral de área en el espacio real se ajusta al espacio de referencia mediante: dxdx =|J|dξdη(3.4) donde el Jacobiano viene definido mediante: [J] = "∂x ∂ξ ∂y ∂ξ ∂x ∂η ∂x ∂η #(3.5) Para ello se definirán x,y tal que: x= N ∑ i=1 Nigxi y= N ∑ i=1 Nigyi (3.6) Donde cabe recalcar que dichas funciones de forma se corresponden únicamente con funciones de forma geométricas. Estas pueden estar asociadas tanto a un espacio de referencia de 6 nodos como uno de 4 nodos, por simplicidad se utilizará este último quedando las funciones de forma: Esto, generado mediante un elemento de forma de referencia rectangular de 4 nodos: x=x1(1−ξ)(1−η)+x2ξ(1−η)+x3ξη +x4(1−ξ)η(3.7) 60 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Tabla 3.2 Funciones de forma geométricas. Nig ∂Nig ∂ξ ∂Nig ∂η i (1−ξ)(1−η)−(1−η)−(1−ξ)1 ξ(1−η) (1−η)−ξ2 ξη η ξ 3 (1−ξ)η−η(1−ξ)4 y=y1(1−ξ)(1−η)+y2ξ(1−η)+y3ξη +y4(1−ξ)η(3.8) Teniendo en cuenta que el elemento de referencia de 4 nodos se transforma en un elemento triangular real de 4 nodos donde los nodos 1 y 4 se colapsan quedando: x4=x1(3.9) y4=y1(3.10) Esto deja como resultado: J=x2(1−η)+x3η−x1y2(1−η)+y3η−y1 ξ(x3−x2)ξ(y3−y2)(3.11) 3.3.2 Evolución del determinante Jacobiano Al aparecer en las diferentes integrales de la ecuación va a ser importante la forma del determinante Jacobiano en los espacios de referencia y real, para ello se muestra en las imágenes a continuación su evolución en el espacio de referencia obtenida mediante Mathematica y Matlab, además de un espacio real triangular de genérico. Es importante tener en cuenta para el cálculo del determinante Jacobiano, que este se realiza de manera lineal por tanto se realiza para un elemento rectangular de 4 nodos. No obstante y a pesar de esto el elemento usado para la matriz de rigidez comprende de un elemento rectangular de 6 nodos por tanto es Figura 3.26 Apreciación espacios formación matriz rigidez y formación Jacobiano . 3.3 Análisis mediante funciones de forma singulares 61 importante como la posterior matriz de rigidez analítica obtenida se realiza mediante el elemento rectangular de 6 nodos ya que la deducción en la primera parte del proyecto de las funciones de forma fue mediante este elemento de referencia, mientras que por otro lado para el Jacobiano no será necesario añadir esta evolución cuadrática en ξya que aportaría alta complejidad a los resultados finales. Por tanto las posiciones en el elemento rectangular de 4 nodos se corresponden con las siguientes posiciones en el elemento rectangular de 6 nodos quedando: x1x2x3x4 y1y2y3y4=x1x3x4x6 y1y3y4y6(3.12) Teniendo en cuenta igualmente que antes que la transformación al espacio real triangular conlleva consigo: x1=x6(3.13) y1=y6(3.14) Quedando un Jacobiano lineal en función de los nodos del elemento rectangular de 6 nodos tal que: J=x3(1−η)+x4η−x1y3(1−η)+y4η−y1 ξ(x4−x3)ξ(y4−y3)(3.15) Es decir para implementación posterior total de la matriz de rigidez se ha decidido usar xi con i=1,...,6 eliminando los nodos medios solo para realizar un Jacobiano lineal. det[J] = (x3(1−η)+x4η−x1)(ξ(y4−y3))−(y3(1−η)+ y4η−y1)(ξ(x4−x3)) (3.16) Partiendo del siguiente mallado: Figura 3.27 Mallado en roseta. Cogiendo el siguiente elemento en el espacio real se obtiene una evolución del determinante Jacobiano tal que: 62 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.28 Espacio real y det[J(ξ,η)]. 3.3.3 Montaje Matriz rigidez Partiendo de la ecuación de Laplace, mediante el término representativo de la matriz de rigidez con los diferentes desplazamientos (u) y desplazamientos virtuales (v): GZΩf∂v ∂x ∂v ∂y"∂u ∂x ∂u ∂y#dxdy (3.17) Donde las posiciones de los nodos en el espacio real se corresponden con los diferentes xi,yi . Utilizando el conjunto de ecuaciones formado por las ecuaciones [3.2-3.6] queda como resultado final: GZ1 0Z1 0∂v ∂ξ ∂v ∂η ([J]−1)T[J]−1"∂u ∂ξ ∂u ∂η #|J|dξdη(3.18) Aquí entra el objetivo de las funciones de forma en desplazamientos para la evaluación de los mismos en los diferentes elementos tal que: u= N ∑ i=1 Niu(ξ,η)ui v= N ∑ i=1 Niu(ξ,η)vi (3.19) 3.3 Análisis mediante funciones de forma singulares 63 Por tanto sus derivadas en el espacio de referencia: ∂u ∂ξ = N ∑ i=1∂Niu(ξ,η) ∂ξ ui ∂v ∂ξ = N ∑ i=1∂Niu(ξ,η) ∂ξ vi ∂u ∂η = N ∑ i=1∂Niu(ξ,η) ∂η ui ∂v ∂η = N ∑ i=1∂Niu(ξ,η) ∂η vi (3.20) Teniendo en cuenta el colapso entre los nodos 1 y 6 en el elemento triangular del espacio real, para cuadrar esto en la ecuación final se generan las matrices tal que: v1... v5G      Z1 0Z1 0        ∂(N1+N6) ∂ξ ∂(N1+N6) ∂η ∂N2 ∂ξ ∂N2 ∂η . . .. . . ∂N5 ∂ξ ∂N5 ∂η        ([J]−1)T[J]−1 ∂(N1+N6) ∂ξ ··· ∂N5 ∂ξ ∂(N1+N6) ∂η ··· ∂N5 ∂η !|J|dξdη          u1 . . . u5    Donde cabe destacar como al ser un espacio de 6 nodos y estar colapsados 2 nodos en el punto 1 del espacio real (x,y) se produce como resultado la suma de ambas funciones de forma como se puede ver en el artículo [10] resultando la matriz final como: Ki j5x5=G      Z1 0Z1 0        ∂(N1+N6) ∂ξ ∂(N1+N6) ∂η ∂N2 ∂ξ ∂N2 ∂η . . .. . . ∂N5 ∂ξ ∂N5 ∂η        ([J]−1)T[J]−1 ∂(N1+N6) ∂ξ ··· ∂N5 ∂ξ ∂(N1+N6) ∂η ··· ∂N5 ∂η !|J|dξdη      5x5 (3.21) Resultado teórico Matriz Rigidez Al resultar como resultado la integral término a término de una matriz 5x5 se ha procedido a realizar los cálculos mediante el Mathematica. Para ello se han definido varias constantes tal que: x1x2x3x4x5x6=a b c d e f (3.22) donde cabe recalcar: x1=x6→a=f(3.23) similarmente: y1y2y3y4y5y6=g h i j k l(3.24) por tanto: y1=y6→g=l(3.25) Todo ello debido al colapso entre el nodo 1 y 6 generando un único nodo en el espacio real formando una única esquina dentro del triángulo. El resultado final será por tanto la integral doble de cada uno de los términos de la integral mostrada en la ecuación 3.21,gracias al código empleado en Mathematica se ha obtenido por tanto una solución analítica exacta integrada de tal manera que: 70 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.31 Mallados con nθ=12,24. Como se ha ido desarrollando en la sección 3, cada uno de los elementos que conforman la sección circular en la roseta de punta de la grieta serán elementos singulares. Como se desarrolla en dicho apartado, cada uno de estos elementos estarán conformados por 6 nodos,además de poseer una tensión variable en el interior del elemento, lo que nos permitirá analizar de manera más específica la tensión de punta de la grieta. El resto de elementos se irá incrementando para observar la convergencia del elemento, y cada uno de ellos estará conformado por 3 nodos, mostrando por tanto una tensión constante en cada uno de ellos. El primer problema a mostrar será el mismo estudiado en la sección 2.2.2, correspondiente con la placa rectangular N-D, problema ya conocido con solución lineal. La sección geométrica de este problema será: Tabla 3.3 Geometría inicial del problema. c(Alto) b(Ancho) a(longitud de grieta) d(longitud interfaz adhesiva) σyz 22 1 1 1 Por tanto a continuación se muestran numerosos resultados para los diferentes mallados, y al comprobar tanto su linealidad como su similitud, tanto entre ellos como con las soluciones anteriores ya existentes en la sección 2.2.2 y conocidas por la elasticidad, se puede verificar la validez de esta nueva matriz de rigidez obtenida para las nuevas funciones de forma analizadas en la parte [1] del trabajo presentes en la tabla (1.1). Figura 3.32 Solución lineal para nθ=24. En esta primera imagen se muestra el mallado completo para un número global de nodos considerable (approx ∼ 400), pero debido al entorpecimiento y dificultad gráfica que esto supone, todas las imágenes mostradas a continuación poseen la orden de Matlab shading interp quedando de tal manera: La similitud entre los diferentes mallados con la solución real da pie a la interpretación correcta del cálculo 3.4 Comprobación matriz de rigidez 71 Figura 3.33 uz(x,y)para mallados con nθ=2,4. Figura 3.34 uz(x,y)para mallados con nθ=12,24. de la matriz de rigidez. No obstante para verificarlo de manera más completa se decide también por el cálculo de un segundo problema ya estudiado con anterioridad en la sección 2.2.2 definido por: Dirichlet condition, x=0; Neumann condition, y=2; (3.52) Figura 3.35 Solución conocida. En la imagen anterior se muestra por tanto la solución conocida, obtenida mediante un mallado completo únicamente con funciones de forma no singulares y elementos no singulares con 3 elementos, a partir de aquí la comparación con el resto de mallados con elementos singulares queda tal como: 72 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.36 uz(x,y)para mallados con nθ=2,4. Figura 3.37 uz(x,y)para mallados con nθ=12,24. No solo los desplazamientos tendrán un papel importante, sino que también la comprobación de las tensiones, debido a su naturaleza variable dentro del elemento, tendrán una alta importancia en el análisis, por ello se muestran además las tensiones de los problemas vistos anteriormente, todo ello teniendo en cuenta: • Las tensiones siguientes se mostrarán para los valores nθ=2,24 , ya que representan los valores más bajos y más altos tomados dentro de este mallado. • Debido a tener tensiones constantes en el interior de todos los elementos menos los elementos pertenecientes a la roseta circular, se realiza la representación con pocos elementos en el mallado global para no dificultar la visualización. • Debido a la naturaleza de los elementos singulares( f(ln) ), tendrán un pequeño pico de tensión en el punto x=1,y=0 donde se encuentran estos elementos. Al resolver problemas cuya tensión es constante en todo el dominio estos darán una tensión constante con un fallo numérico asociado a la precisión del programa de orden O(10−15). • La representación individual de cada uno de los elementos singulares se realiza mediante el comando mesh de Matlab, debido a la naturaleza de su obtención. Esto producirá un pequeño contraste visual el cual no afectará al resultado. 3.4 Comprobación matriz de rigidez 73 El primer resultado corresponde con las tensiones σxz,σyz con nθ=2 , para el problema lineal quedando de resultado: Figura 3.38 σxz,σyz para problema lineal con nθ=2. Como se ha comentado con anterioridad al existir un fallo de la precisión de Matlab con orden O(10−15)∼ O(10−12) el mapa de colores que Matlab asocia a cada elemento es ligeramente diferente, no influyendo en el resultado. Figura 3.39 Muestra del error asociado a la precisión de Matlab en σyz. Se comprueba así, una solución con tensiones constantes en todo el dominio producidas por la generación de un desplazamiento lineal,producido por una tensión σxz unitaria, por tanto los resultados concuerdan con la solución esperada de tensión σxz =1,σyz =0. 74 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Para el valor de mallado de roseta nθ=24 Figura 3.40 σyz,σxz para problema lineal con nθ=24. Posteriormente se muestran las soluciones asociadas al segundo problema analizado no lineal para ambas nθ : Figura 3.41 σyz,σxz para problema lineal con nθ=2. Figura 3.42 σyz,σxz para problema lineal con nθ=24. 3.4 Comprobación matriz de rigidez 75 Además se representa una segunda vista de la tensión σyz , de este segundo problema con nθ=2 tal que: Figura 3.43 σyz. Para comparar la solución anterior se presenta con un mayor mallado sin dificultar el análisis gráfico: Figura 3.44 σyz. 76 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares 3.4.1 Solución para nθ=2 Para este primer mallado explicado con detenimiento en la sección 3.4 se muestran las soluciones en desplazamientos para los siguientes valores de la rigidez de la distribución de muelles: κ=0.1,1,100 Figura 3.45 uz(x,y)para κ=0.1,1. En la siguiente imagen se tiene el mismo problema con κ=100 para 2 mallados diferentes (el menor mallado permite una mejor visualización gráfica). Figura 3.46 uz(x,y)para κ=100 con mayor y menor mallado. A continuación se muestran las tensiones en los elementos singulares. En este caso se trata solo de 2 elementos singulares ( nθ=2 ), por tanto se decide a mostrar su evolución no solo en el espacio real del dominio x-y sino que además se muestra en el espacio ξ−η . Es importante recalcar, como para la representación y obtención de la variación interior de las tensiones en los elementos no singulares, dependerán de un mallado independiente de cada uno de estos elementos. Para mostrar la diferencia entre ambos tipos de mallado se muestran ambos espacios para un mallado: NUM =20 yNUM =200 (3.53) 3.4 Comprobación matriz de rigidez 77 Figura 3.47 σyz(ξ,η)elemento izquierdo y derecho en espacio de referencia ξ,ηcon NUM=20. Figura 3.48 σyz(ξ,η)elemento izquierdo y derecho en espacio de referencia ξ,ηcon NUM=200. Como se puede observar, el mallado interior de dichos elementos es de alta importancia también ya que un mallado interior garantizará mejores resultados, no obstante incrementará también el tiempo de obtención de resultados. A la hora de representarlos en el espacio real x-y triangular conocido se obtiene: Figura 3.49 σyz(x,y)elemento izquierdo y derecho en espacio real x,ycon NUM=200. 78 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Esto, visto en el espacio general junto al resto de tensiones: Figura 3.50 σyz(ξ,η) formado únicamente por elementos no singulares e inclusión de elementos singulares. Sin embargo al tener un mallado con solo dos elementos singulares, no posee un buen comportamiento individual en el interior del elemento como se ve en las anteriores imágenes, por ello se realiza para diferentes valores de nθ 3.4.2 Solución para nθ=24 Numerosos programas capaces de resolver problemas similares o capaces del empleo del método FEM(finite element method), recomiendan el uso de un mallado de punta de grieta mediante una roseta de un número de elementos de roseta mayor de 2 elementos como el estudiado con anterioridad. Esto tiene su base en la capacidad de capturar de manera correcta los desplazamientos y la evolución de las tensiones en dirección radial, por ello se decide realizar el análisis final con un mallado del tipo nθ=24 explicado en la sección 3.4. Figura 3.51 Ejemplos de mallado punta de grieta. Con este mallado en la punta de la grieta se solucionará el problema estudiado y presentado en la sección 3. Al igual que antes lo primero será el análisis de resultados de los desplazamientos, en este caso más exactos que el anterior gracias a este mallado. 3.4 Comprobación matriz de rigidez 79 Figura 3.52 uz(x,y)para κ=1 y 100. Ahora, gracias a la existencia de las funciones de forma singulares, al igual que antes se tendrán tensiones variables dentro de dichos elementos. En las dos imágenes a continuación se representan la tensión global únicamente de los elementos no singulares, para poder observar posteriormente la presencia de estos: Figura 3.53 σyz,σxz para κ=1sin elementos singulares. Analizando individualmente las tensiones en la roseta se tendrán: Figura 3.54 σyz,σxz para κ=1en elementos singulares. 86 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.69 σyz máxima frente a número de nodos global con nθ=24. 3.5 Comprobación final Es importante siempre asegurar los resultados obtenidos por un programa numérico montado y calculado manualmente con valores conocidos deducidos de manera analítica. Por todo ello se decide una última comprobación del funcionamiento del propio programa evolucionado del programa de Alberty, hasta el programa final. Dicho problema poseerá cierta solución para los desplazamientos que serán implementados como condición de contorno en la solución analítica, es decir, las condiciones de contorno impuestas en desplazamientos dependerán del valor UD .Para su comprobación, se representa, primeramente la solución analítica conocida para el siguiente problema y su misma solución con la opción shading interp como comando de Matlab: Figura 3.70 uz(x,y)solución analítica. Esta solución analítica comprende de: Condición Neumann con σyz =0en x∈[0,1],y=0(3.55) Condición Robin en x∈[1,2],y=0(3.56) y condiciones en desplazamientos tales como: 3.5 Comprobación final 87 Figura 3.71 UDborde izquierdo. Figura 3.72 UDborde derecho. Figura 3.73 UDborde superior. La solución numérica será, al igual que antes, con una representación recurrente al comando de Matlab shading interp: 88 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares Figura 3.74 uz(x,y)solución numérica para 10100 nodos y nθ=12. Se observa por tanto una alta similitud entre la solución proporcionada por la formulación analítica para dichas condiciones y la solución numérica obtenida. Es importante además calcular el error global en el sistema y en el propio dominio a partir de la diferencia en el dominio entre ambas soluciones, por tanto el error se definirá como: error =(solución analítica-solución numérica) máximo valor de uz(x,y)(3.57) Como ya se vio en la sección 2 del trabajo, en los problemas lineales como N-R 90 º el error será del orden de la precisión mínima de Matlab(O( ∼10−15 )), no obstante a medida que aumenta la no linealidad aumentará este error numérico, en este caso un error bajo del orden (O(∼10−4)/O(∼10−3)). Figura 3.75 Error en el dominio para 10100 nodos y nθ=12. Como se aprecia en la anterior figura, se obtienen errores del orden O(10−4) en la representación de la solución numérica, aproximando de manera adecuada los desplazamientos y el cumplimiento de las soluciones Dirichlet. Gracias al conocimiento de la solución analítica se sabe que los desplazamientos para x,y=0 para dicho problema se corresponden con: Como resultado, se conoce el desplazamiento para la punta de la grieta, en este caso situada en (x=1,y=0) tal y como:uz(1,0) = 1. 3.5 Comprobación final 89 Figura 3.76 Representación linea y=0 solución analítica del problema. A partir de aquí se analizan dos mallados globales con un mallado refinado formados por 10100,90300 nodos respectivamente, y 3 mallados menos refinados formados por 110,930,2550 nodos, todo ello para varios valores de elementos singulares en la roseta quedando: Tabla 3.4 Valor desplazamiento de punta de grieta para varios mallados. nθ110 nodos 930 nodos 2550 nodos 10100 nodos 90300 nodos 2 0.82300 0.90796 0.934599 0.96021 0.98301 4 0.95384 0.98128 0.98756 0.99284 0.99712 8 0.95354 0.98119 0.98751 0.99282 0.99708 24 0.95349 0.98116 0.98749 0.99278 0.99698 Cabe destacar, como la evolución con el mallado no posee un comportamiento trivial, ya que se comporta de manera creciente en función del número de elementos global frente a precisión sin embargo para ciertos tipos de mallados, el comportamiento con el aumento de número de elementos de roseta es ligeramente peor con nθ creciente a partir de nθ=4 , posiblemente por el comportamiento radial, afecta a la aplicación o rigidez que la condición Robin aporta sobre dicho lado del elemento singular , es por dicho motivo, por el que diferentes programas encargados del cálculo FEM tales como ANSYS o ABAQUS, dedican gran parte de su tiempo al estudio de la influencia de diferentes mallados en los problemas singulares para su mejora. Esto marcaría el final de este estudio, dejando así una posible continuación mencionada en el apartado final: Conclusiones y desarrollos futuros. En la última imagen por tanto se observa la representación analítica y numérica juntas, viendo un comportamiento casi similar excepto por el error de O(10−4) prácticamente inapreciable: Figura 3.77 Representación linea y=0 solución analítica vs numérica del problema. 90 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares 3.6 Conclusiones y desarrollos futuros Desde el punto de vista autocrítico y con el conocimiento de la evolución del desarrollo académico, es importante realizar, no solo una evaluación general del trabajo, sino dar pie para inspirar la continuación del mismo hacia un estudio más detallado, de cara a la total comprensión de dicho elemento y a la posible mejora de los resultados obtenidos en este. Como primer punto de mejora, se conoce las funciones de forma analizadas en el mismo. Para ello, se ha analizado un método iterativo de obtención de funciones de forma sobre un espacio unidimensional de 3 nodos, con la finalidad de adaptación a un elemento rectangular final en el espacio de referencia ξ,η de 6 nodos, lineal en y, y con comportamiento con base 1,r,ln(r) en x. A partir de ello se deduce un amplio y extenso abanico de posibilidades de obtención de sets de funciones de forma. Dicho abanico de posibilidades se extendería mediante el análisis de colocación de 4 o más nodos en el eje x frente a los 3 nodos únicamente estudiados en este trabajo, sumado a las diferentes posibilidades de colocación de dichos nodos dentro de la misma línea unidimensional en comparación con la situación fija de nuestro set de funciones establecido situados en los extremos y punto medio de la línea unidimensional. A este alto número de posibilidades abiertas de obtención de set de funciones de forma en un solo eje habrá que añadir la complicación del segundo eje o eje y. En este análisis dicho eje y ha tenido una influencia lineal de dos nodos, sin embargo un análisis altamente complejo y más detallado permitiría generar el set se funciones de forma con una evolución diferente a la lineal en dicho eje, esta evolución pasaría por cuadrática, cúbica o incluso cualquier base en el espacio de funciones como por ejemplo la base 1,r,rα,rβ. Todas estas posibilidades de intercambio de posiciones dentro de los elementos, sumados al extenso conjunto de funciones de forma de posible análisis a través del cual realizar la primera iteración da lugar a una posible rama de investigación de sets de funciones de forma para elementos de X nodos con base f,f1,f2... donde f representa las funciones base. A partir de la obtención de cada uno de los elementos más complejos cada vez será necesario un análisis completo de cada uno de ellos para observar su relevancia dentro de la obtención de desplazamientos y tensiones. Esto implicaría el cálculo para cualquier set de funciones de forma de la propia matriz de rigidez cuya deducción no sería para nada trivial y digno de un extenso y complejo análisis. Por último, una vez llegado al final del análisis de cualquier tipo de set o pack de funciones de forma a partir de la cual se haya montado los elementos singulares, cabría destacar el papel o la influencia que aportarían diferentes tipos de mallado además de su convergencia exacta. Es por esto último por lo que diferentes programas FEM tales como ANSYS o ABAQUS, dedican gran parte de su tiempo al estudio de la influencia de diferentes mallados en los problemas singulares para su mejora, además de la obtención del mallado óptimo mediante numerosos algoritmos de mallado dentro de cada tipo de problema diferente, siendo el entendimiento de este comportamiento fundamental para la aproximación correcta de la solución final. 3.7 Anexos 91 3.7 Anexos En este apartado de anexos se incluiran los principales códigos empleados tanto en Matlab como en Mathematica. 3.7.1 Código principal de Alberty, funciones de forma no singulares Como ya se comentó con anterioridad, el programa principal ha sido generado a partir de la resolución de la ecuación de Laplace propuesta por Alberty. A partir de aquí se han ido desarrollando diferentes programas capaces de solucionar las diferentes geometrías analizadas, para los elementos no singulares. El siguiente código comprende el programa base capaz de generar la matriz de rigidez A y el vector de fuerzas b a partir de geometrías introducidas junto a sus mallados (elements3.dat,elements4.dat) además de las condiciones nodales de contorno (neumann.dat,dirichlet.dat). Código 3.1 Código principal de Alberty. % FEM2D two-dimensional finite element method for Laplacian. % Initialisation clear all; clc; load coordinates.dat; coordinates(:,1)=[]; eval('load elements3.dat; elements3(:,1)=[];','elements3=[];'); eval('load elements4.dat; elements4(:,1)=[];','elements4=[];'); eval('load neumann.dat; neumann(:,1) = [];','neumann=[];'); load dirichlet.dat; dirichlet(:,1) = []; FreeNodes=setdiff(1:size(coordinates,1),unique(dirichlet)); %ZFD = sparse(size(coordinates,1),size(coordinates,1)); A = sparse(size(coordinates,1),size(coordinates,1)); b = sparse(size(coordinates,1),1); %--------------------- % x=[1 2] % ZDF=f(x) %-------------------- %% % Assembly %%%%elementos triangulares for j = 1:size(elements3,1) A(elements3(j,:),elements3(j,:)) = A(elements3(j,:), ... elements3(j,:))+stima3(coordinates(elements3(j,:),:)); %stima3 devuelve M, que es igual a R de S.H, es decir, la matriz de %rigidez. S.H. multiplica por 0 su matriz de masa. Ver la comparación de las funciones stima3 y local end % %%%%%elementos cuadrangulares % for j = 1:size(elements4,1) % A(elements4(j,:),elements4(j,:)) = A(elements4(j,:), ... % elements4(j,:)) + stima4(coordinates(elements4(j,:),:)); % end A1=A; %% % Volume Forces %%%%elementos triangulares 92 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares %Vamos a sumar fuerzas de volumen en todos los nodos PARA CADA TRIÁNGULO for j = 1:size(elements3,1) b(elements3(j,:)) = b(elements3(j,:)) + ... det([1,1,1; coordinates(elements3(j,:),:)']) * ... %%%jacobiano de la transformación f(sum(coordinates(elements3(j,:),:))/3)/6; %%%f(baricentro)/6*[1;1;1] %SDFA=f(sum(coordinates(elements3(j,:),:))/3)/6 %%%b la monta igual que Sanz Herrera, la única diferencia es que la %%%fuente en Alberty es 1 y en S.H es cero (ver la función f de cada %%%uno) end b %%%%%elementos cuadrangulares % for j = 1:size(elements4,1) % b(elements4(j,:)) = b(elements4(j,:)) + ... % det([1,1,1; coordinates(elements4(j,1:3),:)']) * ... % f(sum(coordinates(elements4(j,:),:))/4)/4; % end %% % Neumann conditions %Aqui vamos a sumar solo los nodos libres for j = 1 : size(neumann,1) b(neumann(j,:))=b(neumann(j,:)) + ... norm(coordinates(neumann(j,1),:)-coordinates(neumann(j,2),:))*...% longitud del lado g(sum(coordinates(neumann(j,:),:))/2)/2; end A3=A %% % Dirichlet conditions u=sparse(size(coordinates,1),1); u(unique(dirichlet))=u_d(coordinates(unique(dirichlet),:)); b=b-A*u; A4=A %% % Computation of the solution u(FreeNodes)=A(FreeNodes,FreeNodes)\b(FreeNodes); %% % graphic representation show(elements3,coordinates,full(u)); % % graphic representation % show(elements4,coordinates,full(u)); El programa es capaz de resolver por tanto el sistema de la forma Ax =b mediante el comando ya existente en Matlab A\by su posterior representación mediante la función @show. Código 3.2 Fuerzas de volumen. function VolumeForce = f(x); %F Volume force in considered domain. % Y = F(X) returns values of forces at N discrete points in the considered 3.7 Anexos 93 % domain. This input data has to be chosen by the user. X has dimension N % x 2 and Y has dimension N x 1. % % % See also FEM2D, U_D, and G. % J. Alberty, C. Carstensen and S. A. Funken 02-11-99 % File <f.m> in $(HOME)/acf/fem2d/ % This volume force is used to compute Fig. 3 in % "Remarks around 50 lines of Matlab: Short finite element % implementation". VolumeForce = ones(size(x,1),1); % VolumeForce = zeros(size(x,1),1); Código 3.3 Valor de condición de contorno Neumann. function Stress = g(x) %G Data on the Neumann boundary % Y = G(X) returns values of the normal-derivative at N discrete points on % the Neumann boundary. This input data has to be choosen by the user. X % has dimension N x 2 and Y has dimension N x 1. % % % See also FEM2D, F, and U_D. % % J. Alberty, C. Carstensen and S. A. Funken 02-11-99 % File <g.m> in $(HOME)/acf/fem2d/ % This Neumann boundary data is used to compute Fig. 3 in % "Remarks around 50 lines of Matlab: Short finite element % implementation" Stress = ones(size(x,1),1); Código 3.4 Solución Alberty Matriz rigidez. function [R,M,jacob]=local(vertices) jacob=abs(det([1 1 1;vertices'])); BKt=vertices(2:3,:)-[vertices(1,:);vertices(1,:)]; CK=inv(BKt)'*inv(BKt); Kxx=(1/2)*[1 -1 0; -1 1 0; 0 0 0]; Kyy=(1/2)*[1 0 -1;0 0 0;-1 0 1]; Kxy=(1/2)*[1 0 -1;-1 0 1;0 0 0]; K0=(1/24)*[2 1 1;1 2 1;1 1 2]; R=jacob*(CK(1,1)*Kxx+CK(2,2)*Kyy+CK(1,2)*(Kxy+Kxy')); M=jacob*K0; % FORMA ALTERNATIVA de Sanz Herrera: % G=[1 1 1;vertices']\[0 0;1 0; 0 1]; % jacob=det([1 1 1;vertices']); % R=jacob*G*G'/2; % M=jacob/24*[2 1 1;1 2 1;1 1 2]; 94 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares % del artículo de Alberty: % d=size(vertices,2); % G=[ones(1,d+1);vertices']\[zeros(1,d);eye(d)]; % M=det([ones(1,d+1);vertices'])*G*G'/prod(1:d); %nota: %d=2, entonces ones(1,d+1)=1 1 1 y [zeros(1,d);eye(d)]=[0 0;1 0; 0 1] %la G de Alberty y de Sanz Herrera son la misma %la matriz M de Alberty es la misma que la R de Sanz Herrera Por último y como se ha deducido en la parte [2] de esta trabajo, el código para las tensiones generadas por elementos no singulares capaces de generar tensiones constantes dentro del elemento: Código 3.5 Tensiones para funciones de forma no singulares. function [sigma_xz,sigma_yz]=sigma_fun(u,x,y) det=(x(2)-x(1))*(y(3)-y(1))-(x(3)-x(1))*(y(2)-y(1)); J=[y(3)-y(1),y(1)-y(2);x(1)-x(3),x(2)-x(1)]; J=J/det; [u(2)-u(1);u(3)-u(1)]; sigma=J*[u(2)-u(1);u(3)-u(1)]; sigma_xz=sigma(1); sigma_yz=sigma(2); end En los dos siguientes códigos se muestra la inclusión final de la condición de Robin para las funciones de forma no singulares: Código 3.6 Inclusión condición de Robin. %% Matriz Robin busca componentes triángulo(det incluye x1 x2 x3) for j=1:numLadosRob nodos_forman_lado=robin(j,:); %esta funcion busca el triangulo al que pertenece T=buscatriangulo(nodos_forman_lado,numTriangulos,triangulos); %pertenece al triangulo T %formado por : nodos p.e(2 4 5) %nos dan las coordenadas de (x2 y2;x4 y4;x5 y5) coorNodos=coordenadas(triangulos(T,:),:); longitudLado=norm(coordenadas(2,:)-coordenadas(1,:)); K(nodos_forman_lado,nodos_forman_lado)=K(nodos_forman_lado,nodos_forman_lado) +... robin_cambio_espacio1(coorNodos,k_robin,triangulos(T,:),nodos_forman_lado); end K=K+KROBIN; Con la función: Código 3.7 Condición de Robin. function R=robin_cambio_espacio(vertices,k_R,Tri,nodos_lado) for i=1:2 busqueda=nodos_lado(i); p(i)=find(Tri==busqueda); end longitud=norm(vertices(p(2),:)-vertices(p(1),:)); 3.7 Anexos 95 NUM=20; x=linspace(vertices(1,1),vertices(2,1),NUM); y=linspace(vertices(1,2),vertices(2,2),NUM); h=linspace(0,longitud,NUM); h1=1/NUM; R(1,1)=k_R*longitud*sum(phi(1,NUM).*phi(1,NUM).*h1); R(1,2)=k_R*longitud*sum(phi(1,NUM).*phi(2,NUM).*h1); R(2,1)=k_R*longitud*sum(phi(2,NUM).*phi(1,NUM).*h1); R(2,2)=k_R*longitud*sum(phi(2,NUM).*phi(2,NUM).*h1); end Código 3.8 Función phi. function phi_sol=phi(t,NUM) %% eleccion funcion de forma if t==1 phi_t=@(xi,eta)((1-eta).*(1-xi)); phi_sol=phi_t(linspace(0,1,NUM),0); elseif t==2 phi_t=@(xi,eta)(xi.*(1-eta)); phi_sol=phi_t(linspace(0,1,NUM),0); end end 3.7.2 Obtención de los valores singulares mediante Mathematica Para la resolución de la matriz de rigidez de elementos singulares, además de las tensiones, se obtienen los siguientes códigos mediante el Mathematica: Código 3.9 Obtención matriz rigidez mediante integración en Mathematica. M6[xi_, eta_] = (1 - xi - xi*Log[xi]/Log[1/2]) M5[xi_, eta_] = 2 xi*(Log[xi]/Log[1/2]) M4[xi_, eta_] = (xi - xi*Log[xi]/Log[1/2]) M3[xi_, eta_] = -(xi - xi*Log[xi]/Log[1/2]) M2[xi_, eta_] = -2 (xi*Log[xi]/Log[1/2]) M1[xi_, eta_] = -(1 - xi - xi*Log[xi]/Log[1/2]) N6[xi_, eta_] = (eta)*(-1 - 1/(Log[1/2]) - Log[xi]/Log[1/2]) N5[xi_, eta_] = (2 eta)*(1/(Log[1/2]) + Log[xi]/Log[1/2]) N4[xi_, eta_] = (eta)*(1 - 1/(Log[1/2]) - Log[xi]/Log[1/2]) N3[xi_, eta_] = (1 - eta)*(1 - 1/(Log[1/2]) - Log[xi]/Log[1/2]) N2[xi_, eta_] = 2*(1 - eta)*(1/(Log[1/2]) + Log[xi]/Log[1/2]) N1[xi_, eta_] = (1 - eta)*(-1 - 1/(Log[1/2]) - Log[xi]/Log[1/2]) J[xi_, eta_] = {{c (1 - eta) + d*eta - a, i (1 - eta) + j*eta - g}, {xi (d - c), xi (j - i)}} determinanteJacobiano[xi_, eta_] = Det[J[xi, eta]] inversaJ[xi_, eta_] = Inverse[J[xi, eta]] inversaJTrans[xi_, eta_] = Transpose[inversaJ[xi, eta]] A[xi_, eta_] = inversaJTrans[xi, eta].inversaJ[xi, eta] CC[xi_, eta_] = {{N1[xi, eta] + N6[xi, eta], M1[xi, eta] + M6[xi, eta]}, {N2[xi, eta], M2[xi, eta]}, {N3[xi, eta], M3[xi, eta]}, {N4[xi, eta], M4[xi, eta]}, {N5[xi, eta], M5[xi, eta]}} 102 Capítulo 3. Análisis grieta mediante funciones de forma no singulares y singulares a=x(1); % b=x(2); c=x(3); d=x(4); % e=x(5); % f=x(1); g=y(1); % h=y(2); i=y(3); j=y(4); % k=y(5); % l=y(6); sigma_yz=@(xi,eta)((2*(-a + c*(1 - eta) + d*eta)*xi*log(xi))/((-c*g*xi +... d*g*xi + a*i*xi - d*i*xi - a*j*xi + c*j*xi)*log(2))-... (2*(-c + d)*(1 - eta)*xi*(-(1/log(2)) - log(xi)/log(2)))/(-c*g*xi +... d*g*xi + a*i*xi - d*i*xi - a*j*xi + c*j*xi))*desplazamientos(2)+... (-((2*(-a + c*(1 - eta) + d*eta)*xi*log(xi))/((-c*g*xi +... d*g*xi + a*i*xi - d*i*xi - a*j*xi + c*j*xi)*log(2)))-... ((2*(-c + d)*eta*xi*(-(1/log(2)) - log(xi)/log(2)))/(-c*g*xi +... d*g*xi +a*i*xi - d*i*xi - a*j*xi + c*j*xi)))*desplazamientos(5)... -desplazamientos(1)*(((-c + d)*xi*((1 - eta)*(-1 + 1/log(2) + log(xi)/log(2)) +... eta*(-1 + 1/log(2) + log(xi)/log(2))))/(-c*g*xi + d*g*xi + a*i*xi -... d*i*xi - a*j*xi + c*j*xi))+... (-(((-c + d)*(1 - eta)*xi*(1 + 1/log(2) + log(xi)/log(2)))/(-c*g*xi +... d*g*xi + a*i*xi - d*i*xi - a*j*xi + c*j*xi))+... ((-a + c*(1 - eta) + d*eta)*(-xi -(xi*log(xi))/log(2)))/(-c*g*xi +... d*g*xi + a*i*xi - d*i*xi - a*j*xi + c*j*xi))*desplazamientos(3)+... (-(((-c + d)*eta*xi*(1 + 1/log(2) + log(xi)/log(2)))/(-c*g*xi + d*g*xi +... a*i*xi - d*i*xi - a*j*xi + c*j*xi))+... ((-a + c*(1 - eta) + d*eta)*(xi + (xi*log(xi))/log(2)))/(-c*g*xi + ... d*g*xi + a*i*xi - d*i*xi - a*j*xi + c*j*xi))*desplazamientos(4); % sigma_xz=@(xi,eta)(((-i + j).*xi.*((1 - eta).*(-1 + 1./log(2) + log(xi)./log (2)) +... % eta.*(-1 + 1./log(2) + log(xi)./log(2))))./(-c*g.*xi + d*g.*xi + a*i.*xi -... % d*i.*xi - a*j.*xi + c*j.*xi)).*desplazamientos(1)+... % ((2.*(g - (1 - eta).*i - eta.*j).*xi.*log(xi))./((-c.*g.*xi +... % d.*g.*xi + a.*i.*xi - d.*i.*xi - a.*j.*xi + c.*j.*xi).*log(2)) +... % (2.*(1 - eta).*(-i + j).*xi.*(-(1./log(2)) - log(xi)./log(2)))./(-c.*g.* xi +... % d.*g.*xi + a.*i.*xi - d.*i.*xi - a.*j.*xi + c.*j.*xi)).*desplazamientos (2)+... % (((1 - eta).*(-i + j).*xi.*(1 + 1./log(2) + log(xi)./log(2)))/(-c.*g.*xi +... % d.*g.*xi + a.*i.*xi - d.*i.*xi - a.*j.*xi +... % c.*j.*xi) + ((g - (1 - eta).*i - eta.*j).*(-xi - (xi.*log(xi))./log(2))) ./(-c*g.*xi +... % d.*g.*xi + a.*i.*xi - d.*i.*xi - a.*j.*xi + c.*j.*xi)).*desplazamientos(3) +... % ((eta.*(-i + j).*xi.*(1 + 1./log(2) + log(xi)./log(2)))./(-c.*g.*xi +... % d.*g.*xi + a.*i.*xi - d.*i.*xi - a.*j.*xi +... % c*j.*xi) + ((g - (1 - eta).*i - eta.*j).*(xi + (xi.*log(xi))./ log(2)))./(- c*g.*xi +... % d*g.*xi + a*i.*xi - d*i.*xi - a*j.*xi + c*j*xi)).*desplazamientos(4)+... 3.7 Anexos 103 % (-((2.*(g - (1 - eta).*i - eta.*j).*xi.*log(xi)))./((-c.*g.*xi +... % d*g.*xi + a*i.*xi - d*i.*xi - a*j.*xi + c*j*xi).*log(2))+... % (2*eta.*(-i + j).*xi.*(-(1./log(2))- log(xi)./log(2)))./(-c.*g.*xi + d.*g .*xi +... % a.*i.*xi - d.*i.*xi - a.*j.*xi + c.*j.*xi)).*desplazamientos(5); Código 3.14 Introducción matriz Robin singular en Matlab. RobMat(1,1)=(4+log(8)*(-5+log(64)))/(54*(log(2))^2); RobMat(1,2)=-4/(27*(log(2))^2)+1/log(4)-2/log(512); RobMat(1,3)=1/6+2/(27*(log(2)^2))-1/log(16); RobMat(2,1)=RobMat(1,2); RobMat(2,2)=8/(27*(log(2))^2); RobMat(2,3)=(2*(-2+log(8)))/(27*(log(2))^2); RobMat(3,1)=RobMat(1,3); RobMat(3,2)=RobMat(2,3); RobMat(3,3)=(2+(-2+log(8))*log(8))/(27*(log(2))^2); Índice de Figuras 1 Modos de carga MFEL I 1.1 Distorsión de un elemento con 8 nodos obtenida de [9] 1 1.2 Elemento a estudiar unidimensional obtenida de [9] 2 1.3 Evolución dξ/dr frente a r 3 1.4 r(ξ)yξ(r)3 1.5 R2(s,t)con p=1/2 y u(s,t)obtenidas de [12] 5 1.6 Espacios real y referencia 10 1.7 Espacios real y referencia para elemento rectangular de 8 nodos 10 1.8 Colapso funciones de forma para formar espacio triangular,sacada de [10] 11 1.9 N2yN2en espacio triangular 11 1.10 N3yN3en espacio triangular 12 1.11 N4yN4en espacio triangular 12 1.12 N5yN5en espacio triangular 12 1.13 N6yN6en espacio triangular 13 1.14 N7yN813 1.15 N1yN7+N8+N1=N1tEn espacio triangular 13 1.16 Espacios real y referencia 14 1.17 Evolución en ξ14 1.18 N2yN2En espacio triangular 15 1.19 N3yN3En espacio triangular 15 1.20 N4yN4En espacio triangular 15 1.21 N5yN5En espacio triangular 16 1.22 N1yN616 1.23 N1en espacio triangular 16 2.1 Grieta en interfaz adhesiva modelada como distribución continua de muelles 17 2.2 Diferentes modos de fractura 17 2.3 Balance de calor en elemento diferencial obtenida de [5] 18 2.4 Problema inicial 21 2.5 Solución de Alberty y Solución de SH 22 2.6 Solución Alberty para f=0,g=123 2.7 Solución Alberty para f=1,g=123 2.8 Descripción del problema 24 2.9 Solución del problema 24 2.10 Descripción del problema 25 2.11 Solución del problema 25 2.12 Descripción del problema 26 2.13 Solución del problema 26 2.14 Descripción del problema 27 2.15 Solución del problema 27 2.16 Problema original κ=0y Lado Robin con κ=20 29 105 106 Índice de Figuras 2.17 Lado Robin con κ=100 29 2.18 Condiciones contorno problema 29 2.19 Problema original κ=0y Lado Robin con κ=10 30 2.20 Problema original κ=0y Lado Robin con κ=10 30 2.21 Problema original κ=0y Lado Robin con κ=10 30 2.22 Tiempo de computación (s) para n=131 2.23 Tiempo de computación (s) para n=231 2.24 Sección circular (D-R) y Sección circular (N-R) 32 2.25 Solución analítica (D-R 90o)u(r,θ)con γ=1,j=1y Solución numérica (D-R 90o)u(r,θ)con γ=1,j=133 2.26 Solución analítica (D-R 90o)u(r,θ)con γ=1,j=1,θ=90oy Zoom solución numérica (D-R 90º) u(r,θ)con γ=1,j=1,θ=90o33 2.27 Error obtenido u(r,θ)con γ=1,j=1,θ=90oy Error obtenido (D-R 90o)u(r,θ)con γ=1,j= 1,θ=90o34 2.28 Tiempo computación (s) n=2 y Tiempo computación n=3 34 2.29 Error global problema N-R 90o34 2.30 Solución analítica (N-R 90o)u(r,θ)con γ=1,j=1y Solución numérica (N-R 90o)u(r,θ)con γ=1,j=135 2.31 u(r,θ=90o)con γ=1,j=1yu(r,θ=0o)con γ=1,j=135 2.32 Tiempo computación (s) n=1 y Tiempo computación (s) n=3 35 2.33 Solución analítica (D-R 45o)u(r,θ)con γ=1,j=1y Solución numérica (D-R 45o)u(r,θ)con γ=1,j=136 2.34 u(r,θ=45o)Solución analítica frente solución numérica 36 2.35 Error global en el dominio 37 2.36 Error total (D-R 45o)u(r,θ)con γ=1,j=1y Error total (D-R 45o)u(r,θ)con γ=1,j=137 2.37 Tiempo computación (s) n=1 y Tiempo computación (s) n=3 37 2.38 Solución analítica (N-R 45o)u(r,θ)con γ=1,j=1y Solución numérica (N-R 45o)u(r,θ)con γ=1,j=138 2.39 u(r,θ=45o)con γ=1,j=1yu(r,θ=0o)con γ=1,j=138 2.40 Error total (N-R 45o)u(r,θ)con γ=1,j=1y Error total (N-R 45o)u(r,θ)con γ=1,j=139 2.41 Tiempo computación (s) n=1 y Tiempo computación (s) n=3 39 2.42 Solución analítica (D-R 135o)u(r,θ)con γ=1,j=1y Solución numérica (D-R 135o)u(r,θ)con γ=1,j=140 2.43 Representación uz(r,θ=135o)Solución analítica frente numérica 40 2.44 Error total (D-R 135o)u(r,θ)con γ=1,j=1y Error total (D-R 135o)u(r,θ)con γ=1,j=140 2.45 Tiempo computación (s) n=1 y Tiempo computación (s) n=3 41 2.46 Solución analítica (N-R 135o)u(r,θ)con γ=1,j=1y Solución numérica (N-R 135o)u(r,θ)con γ=1,j=141 2.47 u(r,θ=45o)con γ=1,j=1yu(r,θ=0o)con γ=1,j=142 2.48 Error total (N-R 135o)u(r,θ)con γ=1,j=1y Error total (N-R 135o)u(r,θ)con γ=1,j=142 2.49 Tiempo computación (s) n=1 y Tiempo computación (s) n=3 42 2.50 Problema lineal solución en desplazamientos 43 2.51 σxz yσyz 44 2.52 Problema 2 solución en desplazamientos 44 2.53 σxz yσyz 44 2.54 Problema 3 solución en desplazamientos 45 2.55 σxz yσyz 45 2.56 σxz yσyz 45 2.57 σrz yσθz46 2.58 σxz yσyz 46 2.59 σrz yσθz46 3.1 Geometría inicial 47 3.2 uz(x,y)con κ=0.148 3.3 uz(x,y)con κ=149 3.4 uz(x,y)con κ=100 49 3.5 uz(x,y=0)para diferentes κ50 Índice de Figuras 107 3.6 uz(x,y=0)para diferentes κ50 3.7 σyz yσxz para k=0.1 51 3.8 σyz yσxz para k=1 51 3.9 σyz yσxz para k=10 51 3.10 σyz yσxz para k=100 52 3.11 Evolución de uz(x,y)frente longitud de grieta (a) 52 3.12 Evolución de σyz 53 3.13 Evolución de σyz frente longitud de grieta (a) 53 3.14 σyzpara k=0.1 y σxz para k=1 54 3.15 σyzpara k=10 y σxz para k=100 54 3.16 σyzpara k=0.1 y σxz para k=1 frente κ·uz(x,y=0)55 3.17 σyzpara k=10 y σxz para k=100 frente κ·uz(x,y=0)55 3.18 σyzpara k=0.1 y σxz para k=1 frente κ·uz(x,y=0)55 3.19 σyzpara k=10 y σxz para k=100 frente κ·uz(x,y=0)56 3.20 σyzMAX para k=0.1 frente a número global de nodos,escalas lineal y semilogarítmica 56 3.21 σyzMAX para k=1 frente a número global de nodos,escalas lineal y semilogarítmica 57 3.22 σyzMAX para k=10 frente a número global de nodos,escalas lineal y semilogarítmica 57 3.23 σyzMAX para k=0.1 frente a número global de nodos,escalas lineal y semilogarítmica 57 3.24 σyzMAX para k=100 frente a número global de nodos, nodos en y=50 58 3.25 σyzMAX para varias k, frente a número global de nodos, nodos en y=10,20 58 3.26 Apreciación espacios formación matriz rigidez y formación Jacobiano 60 3.27 Mallado en roseta 61 3.28 Espacio real y det[J(ξ,η)] 62 3.29 Mallados con nθ=2,4 69 3.30 Mallados con nθ=6,8 69 3.31 Mallados con nθ=12,24 70 3.32 Solución lineal para nθ=24 70 3.33 uz(x,y)para mallados con nθ=2,4 71 3.34 uz(x,y)para mallados con nθ=12,24 71 3.35 Solución conocida 71 3.36 uz(x,y)para mallados con nθ=2,4 72 3.37 uz(x,y)para mallados con nθ=12,24 72 3.38 σxz,σyz para problema lineal con nθ=2 73 3.39 Muestra del error asociado a la precisión de Matlab en σyz 73 3.40 σyz,σxz para problema lineal con nθ=24 74 3.41 σyz,σxz para problema lineal con nθ=2 74 3.42 σyz,σxz para problema lineal con nθ=24 74 3.43 σyz 75 3.44 σyz 75 3.45 uz(x,y)para κ=0.1,176 3.46 uz(x,y)para κ=100 con mayor y menor mallado 76 3.47 σyz(ξ,η)elemento izquierdo y derecho en espacio de referencia ξ,ηcon NUM=20 77 3.48 σyz(ξ,η)elemento izquierdo y derecho en espacio de referencia ξ,ηcon NUM=200 77 3.49 σyz(x,y)elemento izquierdo y derecho en espacio real x,ycon NUM=200 77 3.50 σyz(ξ,η)formado únicamente por elementos no singulares e inclusión de elementos singulares 78 3.51 Ejemplos de mallado punta de grieta 78 3.52 uz(x,y)para κ=1 y 100 79 3.53 σyz,σxz para κ=1sin elementos singulares 79 3.54 σyz,σxz para κ=1en elementos singulares 79 3.55 σyz,σxz para κ=1en todo el dominio 80 3.56 σyz en elemento singular x,y=0 para 2 tipos de mallado 80 3.57 σyz en elemento singular x,y=0 para κ=0.1,181 3.58 σyz en x∈[1,2],y=0para κ=0.1,181 3.59 σyz en x∈[0,2],y=0para κ=100 82 3.60 σyz en x∈[1,2],y=0para κ=1( izquierda ) y κ=100 (derecha) sin elemento singular 82 3.61 σyz en x∈[1,2],y=0para κ=1(izquierda) y κ=100 (derecha) con elemento singular 83 3.62 σyz máxima frente a número de nodos global con nθ=683 108 Índice de Figuras 3.63 σyz máxima frente a número de nodos global con nθ=884 3.64 σyz máxima frente a número de nodos global con nθ=12 84 3.65 σyz máxima frente a número de nodos global con nθ=24 84 3.66 σyz máxima frente a número de nodos global con nθ=685 3.67 σyz máxima frente a número de nodos global con nθ=885 3.68 σyz máxima frente a número de nodos global con nθ=12 85 3.69 σyz máxima frente a número de nodos global con nθ=24 86 3.70 uz(x,y)solución analítica 86 3.71 UDborde izquierdo 87 3.72 UDborde derecho 87 3.73 UDborde superior 87 3.74 uz(x,y)solución numérica para 10100 nodos y nθ=12 88 3.75 Error en el dominio para 10100 nodos y nθ=12 88 3.76 Representación linea y=0 solución analítica del problema 89 3.77 Representación linea y=0 solución analítica vs numérica del problema 89 Índice de Tablas 1.1 Funciones de forma en desplazamientos 9 3.1 Geometría inicial del problema 48 3.2 Funciones de forma geométricas 60 3.3 Geometría inicial del problema 70 3.4 Valor desplazamiento de punta de grieta para varios mallados 89 109 Índice de Códigos 3.1 Código principal de Alberty 91 3.2 Fuerzas de volumen 92 3.3 Valor de condición de contorno Neumann 93 3.4 Solución Alberty Matriz rigidez 93 3.5 Tensiones para funciones de forma no singulares 94 3.6 Inclusión condición de Robin 94 3.7 Condición de Robin 94 3.8 Función phi 95 3.9 Obtención matriz rigidez mediante integración en Mathematica 95 3.10 Obtención tensiones mediante integración en Mathematica 98 3.11 Obtención condición Robin mediante integración en Mathematica 99 3.12 Introducción matriz rigidez singular en Matlab 100 3.13 Introducción tensiones singular en Matlab 101 3.14 Introducción matriz Robin singular en Matlab 103 111