Jhon Freddy Ochoa Avendano~ › 56030 › 7 › 1020729401.2016.pdfJhon Freddy Ochoa Avendano~...
Transcript of Jhon Freddy Ochoa Avendano~ › 56030 › 7 › 1020729401.2016.pdfJhon Freddy Ochoa Avendano~...
Uso del metodo de los elementos finitos con un esquema deexplosion de malla, para la determinacion de la trayectoria de
grietas en materiales microestructurados sometidos a cargamonotonica
Jhon Freddy Ochoa Avendano
Universidad Nacional de Colombia
Facultad de Ingenierıa
Departamento de Ingenierıa Mecanica y Mecatronica
Bogota D.C., Colombia
2016
Uso del metodo de los elementos finitos con un esquema deexplosion de malla, para la determinacion de la trayectoria de
grietas en materiales microestructurados sometidos a cargamonotonica
Jhon Freddy Ochoa Avendano
Tesis presentada como requisito parcial para optar al tıtulo de:
Magıster en Ingenierıa de Materiales y Procesos
Director(a):
Ph.D., MSc., I.M Diego Alexander Garzon Alvarado
Lınea de Investigacion:
Metodos Numericos en Ingenierıa
Grupo de Investigacion:
Grupo de Modelado y Metodos Numericos en Ingenierıa - GNUM
Universidad Nacional de Colombia
Facultad de Ingenierıa
Departamento de Ingenierıa Mecanica y Mecatronica
Bogota D.C., Colombia
2016
Dedicatoria,
A mama, por su incansable apoyo,
amor y confianza.
A Daniel, quien mas que mi her-
mano ha sido un ejemplo para mı.
A mis sobrinos Pipe y Alejo, quie-
nes le aportan tanta felicidad a mi vida.
A mi abuela Marıa, por su infinito carino.
Especialmente a la memoria de papa, a
quien no he dejado de recordar y extranar
desde el momento de su partida.
Jhon Ochoa
”. . . la sensacion de la vida, de la conciencia de
sı mismo, se duplicaba casi en esos instantes,
que no duraban mas que un relampago. La
mente y el corazon se veıan iluminados con
una desusada luz; todas sus emociones, todas
sus dudas, todas las inquietudes parecıan pacifi-
carse de pronto, resolviendose en una suprema
serenidad, plena jubilosa, clara y armonica es-
peranza, plena de razon y de causas definitivas.”
Fragmento de El idiota - Fedor M. Dos-
toevskij
Agradecimientos
Este trabajo es el resultado de un proceso largo que involucro a varias personas. Por tal
motivo quiero expresar mi gratitud hacia estas.
Especial agradecimiento a mi director de tesis PhD Diego Alexader Garzon Alvarado, quien,
con su paciencia, conocimiento y preocupacion, me ayudo tanto en mi crecimiento cientıfico
y personal. Siempre estare agradecido.
Al profesor Dorian Luis Linero Segrera, PhD., quien puso a disposicion gran parte de su
tiempo para asesorarme y guiarme durante el desarrollo del artıculo producto de este traba-
jo.
A mis companeros y amigos de la maestrıa con quienes innumerables veces dialogamos acerca
de nuestras respectivas tesis. Especial agradecimiento a Juan David G., Manuel P. y Hector
C. por su constante apoyo academico y personal.
A toda mi familia, que siempre estuvo pendiente de mi evolucion durante esta etapa y con
quienes deje de compartir varios momentos para poder cumplir con mis responsabilidades.
A mi Alma mater, la Universidad Nacional de Colombia que me ha acogido y brindado
tantas cosas en estos anos.
xi
Resumen
Esta tesis presenta la formulacion implementacion y validacion de un modelo cualitativo
para la prediccion de trayectorias de agrietamiento de solidos bajo condicion de esfuerzo
plano y deformaciones infinitesimales. Se basa en el metodo de los elementos finitos tradicio-
nales con una modificacion en la discretizacion tradicional, donde se introducen elementos
no lineales tipo vınculo rıgido entre las caras de los elementos lineales planos triangulares.
Se logran obtener trayectorias de fisuracion para tres problemas diferentes, con resultados
satisfactorios en comparacion con modelos experimentales. En este metodo se evita la com-
plejidad de un analisis no lineal y cualquier proceso de regeneracion o adaptacion de la malla.
Palabras clave:Mecanica de fractura Computacional; Patrones de agrietamiento de vi-
gas; Analisis no lineal mediante elementos finitos.
Abstract
This thesis presents the formulation, implementation, and validation of a simplified qua-
litative model to determinate the crack path of solids considering static loads, infinitesimal
strain and plane stress condition. This model is based on finite element method with a spe-
cial meshing technique, where non-linear link elements are included between the faces of
the linear triangular elements. The stiffness loss of some link elements represent the crack
opening. Three experimental tests of bending beams are simulated, where the cracking pat-
tern calculated with the proposed numerical model is similar to experimental result. The
simplified formulation of this model avoids a complex non-linear analysis and a regenerative
or adaptive meshes.
Keywords:Computational Fracture Mechanics; Cracking pattern of beams; Non-linear
analysis with finite element method.
Contenido
Agradecimientos IX
Resumen XI
Lista de figuras XIV
Lista de tablas XVII
1 Introduccion 1
1.1 Motivacion y Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Antecedentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Marco teorico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.1 Factor de intensidad de esfuerzos . . . . . . . . . . . . . . . . . . . . 7
1.3.2 Problema elastico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Problema Elastico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Esquema de la tesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Implementacion del Metodo de los elementos finitos con explosion de malla. 13
2.1 Formulacion del problema elastico bidimensional . . . . . . . . . . . . . . . . 13
2.2 Formulacion del elemento vınculo rıgido . . . . . . . . . . . . . . . . . . . . . 17
2.3 Algoritmo de propagacion implementado . . . . . . . . . . . . . . . . . . . . 19
2.4 Modelo implementado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3 Explosion de malla 25
3.1 Discretizacion inicial del problema. . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Definicion de explosion de malla. . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 Funcionamiento del programa de explosion de malla. . . . . . . . . . . . . . 26
3.3.1 Archivos de entrada. . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Archivos de salida. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.3 Variables y arreglos utilizados. . . . . . . . . . . . . . . . . . . . . . . 27
3.3.4 Descripcion del programa. . . . . . . . . . . . . . . . . . . . . . . . . 30
4 Ejemplos Numericos de aplicacion 43
4.1 Ejemplo 1: Viga de tres puntos con entalla excentrica. . . . . . . . . . . . . . 43
4.1.1 Ejemplo 2: Viga de tres puntos con entalla y agujeros excentricos. . . 46
xiv Contenido
4.1.2 Ejemplo 3: Viga con entalla y dos cargas. . . . . . . . . . . . . . . . . 49
5 Conclusiones 51
5.1 Trabajos a partir de esta tesis . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.2 Productos de la tesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Bibliografıa 53
Lista de Figuras
1-1. Pasos para la simulacion del crecimiento de grietas (tomado y modificado de
[1]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1-2. Representacion del comportamiento del campo de desplazamientos y su co-
rrespondiente capo de deformaciones para las tres clases de modelos. Columna
izquierda: perfil del campo desplazamientos, columna derecha: perfil campo
de deformaciones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1-3. Notacion para la definicion del campo de esfuerzos y desplazamientos en el
frente de grieta (tomado y modificado de [2]). . . . . . . . . . . . . . . . . . 7
1-4. Deformaciones longitudinales en un elemento diferencial sometido a esfuerzos
normales.(a) direccion x, (b)direccion y, (c)direccion z. . . . . . . . . . . . . 9
2-1. Union de caras comunes de elementos triangulares lineales por medio de ele-
mentos vınculo rıgido. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2-2. Rotacion elemento link. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2-3. Comparacion de algoritmo general para la prediccion de trayectorias de grietas
y algoritmo implementado en este trabajo. . . . . . . . . . . . . . . . . . . . 19
2-4. Ejemplo de explosion de malla. . . . . . . . . . . . . . . . . . . . . . . . . . 21
2-5. Calculo del esfuerzo nodal promedio. . . . . . . . . . . . . . . . . . . . . . . 22
2-6. Elementos barra asociados al nodo con el maximo esfuerzo principal. . . . . 23
2-7. (a) Angulo de inclinacion del esfuerzo principal maximos θσmax Y angulos
de inclinacion de los elementos barra asociados al nodo p θilink. (b) Error de
angulo en el crecimiento de grieta. . . . . . . . . . . . . . . . . . . . . . . . . 24
3-1. Discretizacion inicial del problema. . . . . . . . . . . . . . . . . . . . . . . . 26
3-2. Malla explotada que genera el codigo. . . . . . . . . . . . . . . . . . . . . . . 26
3-3. Malla convencional con numeracion. . . . . . . . . . . . . . . . . . . . . . . . 31
3-4. Nueva numeracion de los nodos. . . . . . . . . . . . . . . . . . . . . . . . . . 33
3-5. Nuevas conectividades. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3-6. Malla numerada despues de la explosion. . . . . . . . . . . . . . . . . . . . . 38
3-7. a) Baricentro y bisectrices de los elementos triangulares. b) Nueva posicion
de los vertices del triangulo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3-8. a) Malla original. b) Malla despues de todo el proceso. . . . . . . . . . . . . 40
4-1. Geometrıa de la viga de tres puntos con entalla excentrica. . . . . . . . . . . 43
xvi Lista de Figuras
4-2. Resultados ejemplo 1 - Caso I. Simulacion con mallas diferentes. . . . . . . . 44
4-3. Resultados ejemplo 1 - Caso II. Simulacion con malla Fina. . . . . . . . . . . 45
4-4. Geometrıa viga de tres puntos con entalla y agujeros excentricos. . . . . . . . 46
4-5. Resultados ejemplo 2 - Caso I. Simulacion con mallas diferentes. . . . . . . . 47
4-6. Resultados ejemplo 2 - Case I. Evolucion de la simulacion de crecimiento de
fisura. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4-7. Resultados ejemplo 2 - Caso II. Simulacion con mallas diferentes. . . . . . . 49
4-8. Geometrıa viga con entalla y dos cargas. . . . . . . . . . . . . . . . . . . . . 50
4-9. Comparison between presented model and benchmark of Arrea et al. in [3] . 50
Lista de Tablas
3-1. Tabla de variables y arreglos definidos en el programa de explosion de malla. 28
3-2. Tabla de variable y arreglos definidos en el programa. Continuacion. . . . . . 29
3-3. Tabla de variable y arreglos definidos en el programa. Continuacion. . . . . . 30
4-1. Variacion de la distancia a y b en la figura 4-1. . . . . . . . . . . . . . . . . 44
4-2. Mallas utilizadas en ejemplo 1 - Caso I. . . . . . . . . . . . . . . . . . . . . . 44
4-3. Variacion de la distancia a y b en la figura 4-4. . . . . . . . . . . . . . . . . 46
4-4. Mallas utilizadas en el Ejemplo 2 - Caso I. . . . . . . . . . . . . . . . . . . . 47
4-5. Mallas utilizadas en ejemplo 2 - Caso II. . . . . . . . . . . . . . . . . . . . . 48
1 Introduccion
La mecanica de fractura permite analizar el proceso de fisuracion progresiva de un solido
durante la aplicacion de una carga externa. Cuando se utilizan herramientas computaciona-
les para tal fin, se conoce como mecanica de fractura computacional (MFC) y dentro de sus
principales objetivos esta el estudio de la formacion y crecimiento de grietas en los solidos
[1]. La necesidad de predecir trayectorias de fisuracion y propiedades de materiales mediante
herramientas computacionales es planteada por Coterrell [4], pues de esta manera se obtie-
nen resultados que permiten estudiar nuevos materiales reduciendo el costo y tiempo que
conlleva el estudio mediante tecnicas experimentales.
Existen dos formas de predecir trayectorias de fisuracion. La primera se basa en un en-
foque netamente analıtico, donde es necesario conocer los valores del factor de intensidad
de esfuerzos K para el estado de crecimiento de la grieta. Este calculo presenta gran difi-
cultad especialmente en problemas de geometrıa compleja, pues el valor de K entre otras
cosas, depende de la longitud de la misma fisura [5], lo que plantea un nuevo problema a
medida que la grieta avanza en su crecimiento. La segunda alternativa se basa en enfoques
numericos, donde se crea un factor de intensidad de esfuerzos generalizado para geometrıas
y condiciones de frontera arbitrarias [6].
Mediante enfoques numericos es posible abordar una gama mas amplia de problemas. Sin
embargo, como consecuencia de que la fisura representa una singularidad en el campo de
esfuerzos y en el campo de deformaciones, se generan dos situaciones que ameritan conside-
racion especial cuando se utilizan herramientas computacionales [2]. La primera es que se
hace necesario definir una longitud de crecimiento de grieta para cada paso de avance de la
fisura. La segunda es el alto costo computacional que representa obtener el estado tensional
de los solidos fisurados, pues para cada avance de la grieta se hace necesario generar una
nueva discretizacion espacial, lo que implica resolver un nuevo problema despues de incre-
mento de tamano en la fisura [7].
En esta tesis se presenta un enfoque numerico para predecir la trayectoria de fisuras en
solidos fragiles, bajo una condicion de esfuerzo plano y deformaciones infinitesimales. Se
trata de un enfoque discreto y cualitativo, basado en el metodo de los elementos finitos
convencionales bajo una explosion de malla, donde las caras comunes de los elementos trian-
gulares son duplicadas y conectadas con elementos tipo barra. La metodologıa presentada,
2 1 Introduccion
permite predecir la trayectoria de fisuracion del solido sin modificar la discretizacion inicial
a pesar del crecimiento de la grieta. Adicionalmente el enfoque planteado evita calcular los
valores del factores de intensidad de esfuerzo, lo que permite disminuir el costo computacio-
nal durante la solucion del problema.
1.1. Motivacion y Objetivos
Los inicios de la mecanica de fractura se ubican entre los anos 1921 a 1924, cuando el ingenie-
ro A.A. Griffith [8] explico la fractura de materiales fragiles a partir de criterios energeticos.
Logro una excelente aproximacion a los resultados experimentales obtenidos para materiales
casi perfectamente elasticos1. Sin embargo, el aporte mas importante de su trabajo fue la
concluir que la resistencia de los materiales es controlada por las discontinuidades presentes
en el mismo. Ya hacia 1956, se desarrolla el concepto de tasa de liberacion de energıa por
Irwin [9], basado en la teorıa de Griffith. Poco tiempo despues Irwin fijo su atencion en
el trabajo hecho por Westergaard [10] donde se hacıa un analisis complejo de esfuerzos y
desplazamientos en el frente de grieta. Basado en dicho trabajo Irwin [11], logro mostrar que
los esfuerzos y desplazamientos en la punta de la grieta podıan describirse por una constante
que guardaba relacion con la tasa de liberacion de energıa. A esta constante se le dio del
nombre de factor de intensidad de esfuerzos K.
Junto con el desarrollo de estas teorıas y su posible aplicacion a problemas practicos, nacio la
necesidad de buscar alternativas que fueran aplicables a problemas reales. Esto debido a que
rara vez la geometrıa de un cuerpo y el material que lo compone permiten una formulacion
matematica simple. De manera que los metodos numericos comenzaron a ser relacionados
con la mecanica de fractura, dando sus inicios a lo que hoy se conoce como mecanica de
fractura computacional [12].
Segun Ingraffea [1], la simulacion del crecimiento de grietas en un proceso incremental, donde
una serie de pasos son repetidos para lograr la progresion del modelo. Cada paso requiere de
una informacion calculada en un paso anterior, tal y como se muestra en la figura 1-1.
1Materiales que producen muy poca deformacion plastica antes de la fractura, por ejemplo los materiales
ceramicos.
1.1 Motivacion y Objetivos 3
Figura 1-1: Pasos para la simulacion del crecimiento de grietas (tomado y modificado de
[1]).
En la figura 1-1, Ri corresponde al estado inicial del solido. Ri se conoce como estado de
representacion y posee informacion del modelo geometrico del dominio (con o sin grietas),
condiciones de frontera del problema y las propiedades del material (agrietado o no agrieta-
do). Ri se convierte en el estado de analisis Ai a traves de un proceso de discretizacion D
que a su vez contiene una funcion de enmallado M . Ver (1-1).
D(Ri,M(Ri))⇒ Ai (1-1)
El estado de analisis Ai, contiene una descripcion del problema adecuada para ser ingresa-
da en el proceso de solucion S. Este ultimo frecuentemente consiste en un programa para
el calculo de esfuerzos mediante elementos finitos, elementos de contorno o cualquier otro
metodo numerico que permita solucionar el sistema de ecuaciones arrojado por Ai. El objeti-
vo de S es convertir los datos de Ai a un estado de equilibrio Ei el cual consiste en variables
de campo tales como desplazamientos y esfuerzos. Adicionalmente el estado Ei contiene in-
formacion del material tales como factores de intensidad de esfuerzos o parametros de la
fuerza conductora Fi para cada uno de los frentes de grieta. Ver (1-2).
S(Ai)⇒ Ei, Fi (1-2)
El siguiente paso es una funcion de actualizacion U , cuyo objetivo es generar un nuevo estado
representacion Ri+1 a partir de la informacion inicial Ri y las variables contenidas en Ei.
U incluye a su vez funcion de representacion del incremento del tamano de la grieta C.
Este ultimo por su lado determina la direccion y la longitud del crecimiento de la grieta,
generalmente a partir del calculo del factor de intensidad de esfuerzos K. Ver 1-3.
U(Ei, Ri, C(Fi))⇒ Ri+1 (1-3)
El proceso antes mencionado es iterativo y se repite hasta que se logra una condicion deter-
minada por el usuario. Presenta principalmente dos componentes que resultan en alto costo
computacional cuando se intenta solucionar el problema mediante elementos finitos [13]. El
primero es el hecho de que el proceso de discretizacion D deba llevarse a cabo para cada
4 1 Introduccion
paso de crecimiento de la grieta, pues la fisura no puede atravesar el elemento finito porque
ya no se cumplirıa el principio del continuo [14]. En la implementacion computacional esto
implica un procedimiento de reenmallado y aplicacion de condiciones de contorno para cada
iteracion. El segundo es que de igual manera, el calculo de K en la funcion de actualizacion
U se debe hacer para cada avance de la grieta. Como se muestra el la seccion 1.3, el calculo
del factor de intensidad de esfuerzos se da a partir de la solucion del estado de tensiones del
dominio. Dado que la fisura representa una singularidad del campo de esfuerzos y deforma-
ciones en el dominio del problema, en la mayorıa de casos se requiere realizar un refinamiento
de malla en los alrededores del frente de grieta, lo que aumenta significativamente la cantidad
de operaciones que debe realizar el computador.
Surgio entonces la motivacion en predecir la trayectoria de fisuras a partir de un enfoque
de elementos finitos, que evitara cualquier proceso de regeneracion o adaptacion de malla,
buscando simplificar o suprimir el calculo del factor de intensidad de tensiones K, con el
objetivo general de reducir el costo computacional. De manera que en esta tesis se planteo
proponer un algoritmo para la prediccion de la trayectoria de grietas en materiales sometidos
a carga monotonica, haciendo uso del metodo de los elementos finitos bajo un esquema de
explosion de malla, a traves de los siguientes objetivos especıficos:
Implementar un algoritmo de solucion basado en el metodo de los elementos finitos
bajo un esquema de explosion de malla, para lograr una correccion en la trayectoria
de grietas en zonas crıticas de la microestructura del material.
Realizar simulaciones de la fractura de diferentes especımenes comparando los resulta-
dos con los presentados por otros autores.
1.2. Antecedentes
El comportamiento de materiales quasifragiles sometidos a solicitaciones mecanicas se carac-
teriza por un proceso de microagrietamiento en zonas caracterısticas, que posteriormente se
traduce en la aparicion de fisuras macroscopicas que se propagan para liberar esfuerzo [15]. A
pesar de los avances en la modelacion teorica y computacional logrados en los ultimos anos,
el problema de localizacion y evolucion se grietas sigue siendo complejo. De manera general,
el proceso de fractura en un solido puede ser representado de acuerdo a la descripcion del
comportamiento del campo de desplazamientos y deformaciones, figura 1-2. Dependiendo
del tipo de comportamiento de los campos de desplazamientos y deformaciones, los modelos
pueden ser divididos en tres [16]:
1.2 Antecedentes 5
(a) Discontinuidad Fuerte
(b) Banda de deformacion separada por dos discontinui-
dades debiles
(c) Perfıl continuo de deformaciones localizadas
Figura 1-2: Representacion del comportamiento del campo de desplazamientos y su co-
rrespondiente capo de deformaciones para las tres clases de modelos. Columna
izquierda: perfil del campo desplazamientos, columna derecha: perfil campo de
deformaciones.
1. Modelos con propagacion cohesiva de discontinuidades 2. El campo de desplazamientos
y el campo de deformaciones presenta un salto en un punto del dominio (en el caso
de 1-D), lo que funciona como una discontinuidad fuerte, figura 1.2(a). La hipotesis es
que el material se mantiene lineal elastico durante la formacion de la zona de proceso
de fractura (FPZ3). Se definen fuerzas de cohesion entre las caras de la discontinui-
dad mediante una ley de separacion; tales fuerzas desaparecen cuando las superficies
libres de la fisura superan una distancia determinada. De esta manera la singularidad
es removida y se da una formacion de FPZ con una longitud caracterıstica; hay un
crecimiento discreto de la fisura [15].
2Traduccion de Models with propagating cohesive discontinuities3Por sus siglas en ingles ”Fracture Process Zone”
6 1 Introduccion
La implementacion computacional en esta clase de modelos se ha dado con el metodo
de los elementos finitos(FEM) mediante la utilizacion de elementos especiales. Tam-
bien varias formas del metodo de los elementos de contorno (BEM) [17]. Para tratar la
discontinuidad de los desplazamientos, se requiere de reenmallado del solido para cada
avance de la fisura, por lo que varios autores han presentado tecnicas de reenmallado
[18, 19, 20, 21]. Una alternativa al reenmallado es utilizar discontinuidades embebidas
en los elementos finitos. Ası la interpolacion es enriquecida con terminos que repro-
ducen un salto en el campo de desplazamientos [22, 23, 24]. Sin embargo el metodo
se restringe a ciertos tamanos y modos de elementos para garantizar la unicidad de
la respuesta de los mismos [16]. Esto impide la utilizacion del enfoque en problemas
de tres dimensiones. Para evitar esta situacion algunos autores propusieron incorporar
el concepto de particion de unidad4 en el enriquecimiento de la interpolacion [25, 26].
En el caso de elementos finitos estandar, las funciones de forma son enriquecidas con
funciones seleccionadas por el usuario. Al seleccionar funciones de Heaviside por ejem-
plo, corresponde al Metodo de los Elementos Finitos Extendidos (X-FEM), metodo
ampliamente utilizado en propagacion de fisuras en tres dimensiones [27, 28].
2. Modelos de continuidad suave con regularizacion parcial.5. Son modelos que represen-
tan la FPZ mediante la localizacion una banda finita de deformacion, figura 1.2(b).
Aunque la descripcion del campo de desplazamientos es continua, hay presencia de
discontinuidades debiles entre las fronteras de la banda de material de la FPZ y el
material circundante [15]. El material de la banda de deformacion puede ser considera-
do con propiedades independientes en una simulacion mediante elementos finitos. En
principio Bazant fijo el tamano del elemento tan grande como el ancho de la FPZ [29],
lo que ocasiona una dependencia entre la malla y la banda de ablandamiento6; poste-
riormente la tecnica fue mejorada al ajustar el modulo de ablandamiento acorde con
el tamano del elemento. Lo anterior corresponde al enfoque de energıa de fractura, ya
que se intenta reducir la disipacion de energıa en la banda de ablandamiento [30, 31].
3. Modelos de continuidad suave con regularizacion completa 7. En este tipo de modelos
hay continuidad tanto en el campo de desplazamientos como en el de deformaciones.
La FPZ se representa por una banda de material con ablandamiento en la que la de-
formacion aumenta gradualmente; desde un valor mınimo en la frontera de la banda
hasta un valor maximo en el centro de la misma [15].
Los modelos de regularizacion completa generalmente son costosos desde el punto de
vista computacional, pero proporcionan varias mejoras con respecto a los modelos ba-
4Partition-of-unity5Traduccion de Softening continuum models with partial regularization6Traduccion de Softening Band7Traduccion de Regularizated Softenig Continua
1.3 Marco teorico 7
sados en el concepto de energıa de fractura. Si la discretizacion es suficientemente fina,
la influencia de la malla sobre el proceso de fractura es sustancialmente evitada. Este
tipo de modelos funcionan muy bien con esquemas de malla adaptativa [32]. El ajuste
adaptativo puede ser por refinamiento h donde el tamano de los elementos es cambiado
manteniendo el orden de los elementos constantes, o por refinamiento p en el que se
cambia el orden de las funciones de forma del elemento.
1.3. Marco teorico
1.3.1. Factor de intensidad de esfuerzos
La presencia de una fisura en un solido representa una singularidad en los campos de de-
formaciones y esfuerzos. La solucion analıtica del estado tensorial en el frente de grieta es
ya bien conocida y se pueden encontrar en los libros de mecanica de fractura como [33]. La
figura 1-3 representa el analisis del campo de esfuerzos y desplazamientos en dos dimensiones
para un solido agrietado.
Figura 1-3: Notacion para la definicion del campo de esfuerzos y desplazamientos en el
frente de grieta (tomado y modificado de [2]).
La expresiones 1-4,1-5 y 1-6, representan los esfuerzos de un elemento diferencial en el frente
de grieta para un material lineal elastico e isotropico:
σ11 =KI√2πr
cosθ
2
(1− sin
θ
2sin
3θ
2
)− KII√
2πrsin
θ
2
(2− cos
θ
2cos
3θ
2
)(1-4)
σ22 =KI√2πr
cosθ
2
(1 + sin
θ
2sin
3θ
2
)+
KII√2πr
sinθ
2cos
θ
2cos
3θ
2(1-5)
8 1 Introduccion
σ12 =KI√2πr
sinθ
2cos
θ
2cos
3θ
2+
KI√2πr
cosθ
2
(1− sin
θ
2sin
3θ
2
)(1-6)
De igual manera, el campo de desplazamientos para el mismo elemento diferencial de la
figura 1-3, es de la forma:
u1 =KI
2µ
√r
2πcos
θ
2(κ− cos θ) +
KII
2µ
√r
2πsin
θ
2(2 + κ+ cos θ) (1-7)
u2 =KI
2µ
√r
2πsin
θ
2(κ− cos θ) +
KII
2µ
√r
2πcos
θ
2(2− κ− cos θ) (1-8)
En 1-7 y 1-8, µ corresponde al modulo de rigidez del material, mientras que κ recibe el
nombre de constante de Kolosov y su valor depende del estado tensional del solido:
µ =
3− 4ν Deformacion plana3−ν1+ν
Esfuerzo plano(1-9)
En 1-9, ν corresponde al coeficiente de Poisson. Adicionalmente, las anteriores expresiones
estan en funcion de los factores de intensidad de esfuerzo K, por lo tanto la mayorıa de
trabajos se han enfocado en encontrar soluciones de K para diferentes tipos de problemas.
1.3.2. Problema elastico
1.4. Problema Elastico.
La respuesta elastica de un material isotropico se puede definir mediante el modulo de Young
E y el coeficiente de Poisson ν [34]. Con la ley de Hooke se establece una relacion lineal entre
el esfuerzo normal y la deformacion lineal en la misma direccion [35], donde precisamente la
constante de proporcionalidad de dicha relacion es E. Ambas constantes se determinan por
medio de ensayos mecanicos donde un especimen es cargado en la direccion axial mientras
se mide su elongacion.
Al ser el modulo de Young la pendiente de la curva (σ, ε), la deformacion del especimen
podrıa tomar la forma:
ε =σ
E(1-10)
Por su parte el coeficiente de Poisson relaciona la deformacion transversal y la deformacion
axial por medio de un cociente:
ν =εtransversalεaxial
(1-11)
Tomando un elemento diferencial del material al que se aplican cargas externas, se generan
esfuerzos en las direcciones σxx, σyy y σzz como se muestra en la figura 1-4. Como conse-
cuencia de las cargas aplicadas, se presentan deformaciones longitudinales εxx = σxxE
en la
1.4 Problema Elastico. 9
direccion x, εyy = σxxE
en la direccion y y εzz = σzzE
en la direccion z. Adicionalmente se
presentan las deformaciones laterales en la direccion x : εyy = εzz = −ν σxxE
, en la direccion
y: εxx = εzz = −ν σyyE
y en la direccion z: εxx = εyy = −ν σzzE
.
Figura 1-4: Deformaciones longitudinales en un elemento diferencial sometido a esfuerzos
normales.(a) direccion x, (b)direccion y, (c)direccion z.
Si se considera que los esfuerzos cortantes solo producen deformaciones angulares, las defor-
maciones longitudinales debidas a los esfuerzos normales estarıa determinadas por:
εxx =1
E(σxx − νσyy − νσzz)
εyy =1
E(−νσxx + σyy − νσzz) (1-12)
εzz =1
E(−νσxx − νσyy + σzz)
Al despejar las componentes de esfuerzo normal de las expresiones (1-12) se obtendrıa:
σxx =E
(1 + ν)(1− 2ν)[(1− ν)εxx + νεyy + νεzz]
σyy =E
(1 + ν)(1− 2ν)[(νεxx + (1− ν)εyy + νεzz] (1-13)
σzz =E
(1 + ν)(1− 2ν)[(νεxx + νεyy + [(1− ν)εzz]
Por otra parte, al formular la ley de Hooke para esfuerzo cortante, se establece una relacion
entre las deformaciones angulares y el esfuerzo cortante en la pieza, donde la constante de
proporcionalidad es G y se conoce como modulo de elasticidad a cortante. Este se relaciona
10 1 Introduccion
con el modulo de Young como se muestra en (1-14)
G =E
2(1 + ν)(1-14)
La deformacion angular del elemento diferencial de la figura 1-4, se expresarıa en el plano
xy como γxy, γxz en el plano xz y γyz en el plano yz de la siguiente manera:
γxy =σxyG
γxz =σxzG
(1-15)
γyz =σyzG
(1-16)
De igual manera, si en la anterior expresion se despejan las componentes de los esfuerzos
cortantes se tendrıa:
σxy =E
2(1 + ν)γxy
σxz =E
2(1 + ν)γxz (1-17)
σyz =E
2(1 + ν)γyz
Las componentes de esfuerzos normales mostradas en (1-13) y las componentes de esfuerzo
cortante mostradas en (1-17), se pueden escribir en forma matricial de la siguiente manera:
σxxσyyσzzσxyσxzσyz
=
E
(1 + ν)(1− 2ν)
(1− ν) ν ν 0 0 0
(1− ν) ν 0 0 0
(1− ν) 0 0 012(1− 2ν) 0 0
12(1− 2ν) 0
sim 12(1− 2ν)
εxxεyyεzzγxyγxzγyz
De acuerdo a lo anterior se define la matriz constitutiva de un material elastico e isotropico
D como:
D =E
(1 + ν)(1− 2ν)
(1− ν) ν ν 0 0 0
(1− ν) ν 0 0 0
(1− ν) 0 0 012(1− 2ν) 0 0
12(1− 2ν) 0
sim 12(1− 2ν)
La ecuacion constitutiva que relaciona los esfuerzos del material con sus deformaciones to-
marıa la forma:
σ = Dε (1-18)
1.5 Esquema de la tesis 11
1.5. Esquema de la tesis
La tesis esta organizada de la siguiente manera: En el capıtulo 2 se presenta la formulacion
matematica del problema elastico mediante elementos finitos con explosion de malla. Ini-
cialmente se realiza el analisis de las implicaciones de la explosion de la malla en el metodo
de los elementos finitos, posteriormente se plantean el problema en forma matricial para los
elementos planos y los elementos tipo barra o vınculo rıgido generados como consecuencia
de la explosion de la malla. En el capıtulo 3, se expone el algoritmo desarrollado para llevar
a cabo la explosion de la malla. Posteriormente, en el capıtulo 4, se presentan los ejem-
plos numericos desarrollados en la investigacion. Al final, en el capitulo 5 se presentan las
conclusiones de la tesis desarrollada.
2 Implementacion del Metodo de los
elementos finitos con explosion de
malla.
El objetivo de este capıtulo es exponer la metodologıa propuesta para la prediccion de
trayectorias de fisuracion. Inicialmente se presentan de manera individual la formulacion
tanto de elementos planos como de elementos barra. Posterior a esto, se describe de manera
general la metodologıa desarrollada, mostrando cada uno de los pasos necesarios para realizar
la simulacion de crecimiento de grieta.
2.1. Formulacion del problema elastico bidimensional
La condicion de equilibrio de un solido de dominio Ω, rodeado por contorno Γ se expresa
mediante:
∇ · σ + b = 0, ∀x ε Ω (2-1)
Donde:
σ · n = t∗, ∀x ε Γ (2-2)
El tensor σ corresponde al tensor de esfuerzos. Por su parte b hace referencia a las fuerzas
de volumen aplicadas al solido y t∗ a las fuerzas por unidad de area aplicadas sobre una
superficie de normal n.
De acuerdo con le principio de los trabajo virtuales [35], dado un vector de desplazamiento
virtual δυ, se presenta (2-1) de forma debil ası:∫Ω
δυ · (∇ · σ + b) dΩ = 0 (2-3)
Donde υ es una funcion de peso que correponde a un tensor de orden 1. Aplicando el principio
de la divergencia, la integral anterior se puede escribir de la forma:
−∫Ω
∇δυ : σ dΩ +
∫Γ
δυ · (σ · n) dΓ +
∫Ω
δυ · b dΩ = 0 (2-4)
14 2 Implementacion del Metodo de los elementos finitos con explosion de malla.
Al ser υ un tensor de orden 1, el termino ∇υ, puede descomponerse en una parte simetrica
y una parte antisimetrica como se muestra en 2-5.
∇υ =1
2(∇υ +∇υT ) +
1
2(∇υ −∇υT ) (2-5)
Se definen la parte simetrica y la parte antisimetrica como S y AS respectivamente:
S =1
2(∇δυ +∇δυT )
(2-6)
AS =1
2(∇δυ −∇δυT )
El termino ∇υ : σ de la ecuacion (2-4) se puede escribir de la siguiente forma:
σ : ∇υ = σ : (S + AS) = σ : S + σ : AS (2-7)
El producto σ : AS = 0, mientras la parte simetrica del gradiente del desplazamiento virtual
se define como la deformacion virtual δε ası:
S =1
2(∇δυ +∇δυT ) = δε (2-8)
Por lo tanto la ecuacion (2-4) en forma incremental quedarıa de la forma:
−∫Ω
δε : σ dΩ +
∫Γ
˙δυ · (σ · n) dΓ +
∫Ω
˙δυ · b dΩ = 0 (2-9)
El modelo propuesto en este trabajo supone que los subdominios fuera de la zona de fisura-
cion conservan un comportamiento lineal elastico durante carga y descarga, mientras que en
la zona de fisura se pierde totalmente la cohesion entre sus caras.
En general, la tasa del esfuerzo puede expresarse como σ = D : ε, donde D es un tensor
de cuarto orden que corresponde al tensor constitutivo tangente del material. Adicionalmen-
te (2-2) en forma incremental se escribe como σ · n = t∗, donde t
∗es la tasa de las fuerzas
de traccion de superficie. Reemplazando estos valores:
∫Ω
δε : (D : ε) dΩ−∫Γ
˙δυ · t∗ dΓ−∫Ω
˙δυ · b dΩ = 0 (2-10)
Para elementos finitos bidimensionales, los tensores ε, dε y ˙δυ de la ecuacion (2-10) se
presentan en forma matricial ası:
ε =
∂∂x
0
0 ∂∂y
∂∂y
∂∂x
[UxUy
](2-11)
2.1 Formulacion del problema elastico bidimensional 15
δε =
∂∂x
0
0 ∂∂y
∂∂y
∂∂x
[CxCy
](2-12)
˙δυ =
[CxCy
](2-13)
En las anteriores expresiones, el vector [Ux Uy]T es el campo de la tasa de desplazamientos
reales en dos dimensiones, mientras que el vector [Cx Cy]T es campo de la tasa de desplaza-
mientos virtuales.
Este modelo supone que la fisura se presenta entre los lados de los elementos bidimen-
sionales. Por lo tanto, el dominio es discretizado con elementos triangulares lineales que
representan las zonas excluidas de fisura. Para este tipo de elementos el vector de la tasa de
desplazamientos reales y virtuales son de la forma:
[UxUy
]=
[N1 0 N2 0 N3 0
0 N1 0 N2 0 N3
]
u1x
u1y
u2x
u2y
u3x
u3y
(2-14)
[CxCy
]=
[N1 0 N2 0 N3 0
0 N1 0 N2 0 N3
]
c1x
c1y
c2x
c2y
c3x
c3y
(2-15)
Siendo U = [u1x, u
1y, · · · , u3
y] y C = [c1x, c
1y, · · · , c3
y] los vectores de valores nodales de des-
plazamiento real y virtual del elemento. Los terminos Ni son las funciones de forma estandar
del elemento triangular lineal en el nodo i.
Las tres componentes de deformacion en el plano se presentan en un vector de la forma:
ε =
∂∂x
0
0 ∂∂y
∂∂y
∂∂x
[N1 0 N2 0 N3 0
0 N1 0 N2 0 N3
]
u1x
u1y
u2x
u2y
u3x
u3y
(2-16)
16 2 Implementacion del Metodo de los elementos finitos con explosion de malla.
O lo que es equivalente a:
ε =
∂N1
∂x0 ∂N2
∂x0 ∂N3
∂x0
0 ∂N1
∂y0 ∂N2
∂y0 ∂N3
∂y∂N1
∂y∂N1
∂x∂N2
∂y∂N2
∂x∂N3
∂y∂N3
∂x
u1x
u1y
u2x
u2y
u3x
u3y
(2-17)
Se define B como la matriz de operadores diferenciales actuando sobre las funciones de
forma ası:
[B] =
∂N1
∂x0 ∂N2
∂x0 ∂N3
∂x0
0 ∂N1
∂y0 ∂N2
∂y0 ∂N3
∂y∂N1
∂y∂N1
∂x∂N2
∂y∂N2
∂x∂N3
∂y∂N3
∂x
(2-18)
El termino δε : (D : ε) de la ecuacion 2-10 se pude escribir como:
[[c1x c1
y c2x c2
y c3x c3
y]]
∂N1
∂x0 ∂N1
∂y
0 ∂N1
∂y∂N1
∂x∂N2
∂x0 ∂N2
∂y
0 ∂N2
∂y∂N2
∂x
D
∂N1
∂x0 ∂N2
∂x0 ∂N3
∂x0
0 ∂N1
∂y0 ∂N2
∂y0 ∂N3
∂y∂N1
∂y∂N1
∂x∂N2
∂y∂N2
∂x∂N3
∂y∂N3
∂x
u1x
u1y
u2x
u2y
u3x
u3y
(2-19)
En forma tensorial la expresion anterior toma la forma:
C BT D B U (2-20)
De manera que la ecuacion (2-10) puede escribirse de la siguiente manera:
C
∫Ω
BT D B dΩ
U − C∫
Γ
NT t∗dΓ
− C∫
Ω
NT b dΩ
= 0 (2-21)
En la anterior expresion, el vector C tiene un valor arbitrario, constante y diferente de cero
asociado al desplazamiento virtual, el cual se cancela de la ecuacion. A su vez, la expresion
representa un problema algebraico de la forma:
[K]U − Q − P = 0 (2-22)
donde:
K =
∫Ω
BT D B dΩ (2-23)
2.2 Formulacion del elemento vınculo rıgido 17
Q =
∫Γ
NT t∗dΓ (2-24)
P =
∫Ω
NT b dΩ (2-25)
2.2. Formulacion del elemento vınculo rıgido
En este modelo se incorporan elementos vınculo rıgido que unen los nodos de las caras
originalmente comunes entre dos elementos triangulares lineales. La longitud de este tipo de
elementos es sustancialmente menor que la dimension promedio de los elementos triangulares
del resto del dominio, como lo indica sin escala la figura 2-1.
Figura 2-1: Union de caras comunes de elementos triangulares lineales por medio de ele-
mentos vınculo rıgido.
En los elementos vınculo rıgido se establece un desplazamiento axial en el sistema de coor-
denadas local [X ′, Y ′] en cada uno de los nodos extremos [u′1x , u′1y , u
′2x , u
′2y ], que descompone
vectorialmente al sistema global de coordenadas [X, Y ] contenido en el plano [u1x, u
1y, u
2x, u
2y],
como lo indica la figura 2-2.
18 2 Implementacion del Metodo de los elementos finitos con explosion de malla.
Figura 2-2: Rotacion elemento link.
La relacion de los desplazamientos entre los sistemas global y local se representa como:u′1xu′1yu′2xu′2y
=
cosθ sinθ 0 0
−sinθ cosθ 0 0
0 0 cosθ sinθ
0 0 −sinθ cosθ
u1x
u1y
u2x
u2y
Se define la matriz de transformacion [T ] como:
[T ] =
cosθ sinθ 0 0
−sinθ cosθ 0 0
0 0 cosθ sinθ
0 0 −sinθ cosθ
Las fuerzas y desplazamientos locales se transforman al sistema global por medio de f ′ =
[T ]f y u′ = [T ]u respectivamente.
El comportamiento del elemento vınculo rıgido se considera perfectamente elastico en la
metodologıa aca utilizada. Entonces la relacion entre las cargas y desplazamientos nodales
en el sistema local toman la forma: f ′ =[K ′]u′. En el sistema global es de la forma:
f =[K]u. Donde k es la rigidez del elemento vınculo rıgido y [K ′] la forma matricial
de k. La rigidez se define ası:
k =
→∞ si σlink < σu0 si σlink > σu
; [K ′] = k
1 0 −1 0
0 0 0 0
−1 0 1 0
0 0 0 0
(2-26)
2.3 Algoritmo de propagacion implementado 19
Dado que [T−1] = [T T ], se obtiene una expresion que relaciona las cargas globales f con
los desplazamientos locales u′:
[T T ][K ′][T ]u′ = f (2-27)
2.3. Algoritmo de propagacion implementado
En la seccion 1.1, se presento la figura 1-1 la cual expone el algoritmo tıpico para predecir
la trayectoria de una grieta. Se hablo de que para cada avance de la grieta, era necesario
hacer una iteracion del algoritmo presentado por Ingraffea [1]. Se expusieron dos situaciones
del algoritmo que generan costo computacional. La primera, el hecho de realizar para cada
iteracion el proceso de discretizacion D junto con su funcion de enmallado M . La segunda,
la ejecucion de la funcion U para cada iteracion, que generalmente corresponde al calculo
del factor de intensidad de esfuerzos K, ver figura 2.3(a).
(a) Algoritmo propuesto por Ingraffea
[1]
(b) Algorıtmo propuesto en la
tesis
Figura 2-3: Comparacion de algoritmo general para la prediccion de trayectorias de grietas
y algoritmo implementado en este trabajo.
En la figura 2.3(b) se presenta el algoritmo implementado en esta tesis. Inicia con un estado
representacion del solido Ri equivalente al presentado por el algoritmo de la figura 2.3(a),
pues se incluyen las condiciones de frontera y propiedades del material.Siguiendo la trayec-
toria punteada se llega al estado de analisis Ai, pasando antes a traves de un proceso de
discretizacion convencional D que contiene la funcion de enmallado M y el proceso de ex-
plosion de malla EX. Las trayectorias punteadas de Ri hasta D, de D hasta EX y de EX
hasta Ai indican que solo se recorren para el primer avance de grieta. De manera analoga
la elipse punteada que encierra a D y a M , indica que el proceso de discretizacion solo se
20 2 Implementacion del Metodo de los elementos finitos con explosion de malla.
realiza para el primer avance de crecimiento de la fisura.
En el estado Ai se tiene la descripcion del problema, que se resuelve por medio de la funcion
de solucion S similar a lo mencionado en la seccion 1.1. En el caso del modelo aca presentado,
en S se obtiene la solucion del campo de esfuerzos y deformaciones mediante elementos fini-
tos. Dichos valores se llevan al estado de equilibrio Ei donde se realiza el postprocesamiento
de los datos obtenidos en S. El postprocesamiento consiste basicamente en encontrar el valor
y direccion del primer esfuerzo principal maximo y en ubicar el nodo al que le corresponden
dichos valores.
La funcion C identifica la direccion de propagacion de la fisura. Para tal fin, la informa-
cion obtenida en Ei es analizada, junto con informacion proveniente de la malla explotada
en Ei. En este punto es importante mencionar que en C no es necesario calcular el factor de
intensidad de esfuerzos, lo que reduce de manera sustancial la cantidad de operaciones que
se realizan en la ejecucion del algoritmo.
El resultado de C, es un estado inicial Ri+1 que puede ser llevado nuevamente al proce-
so de analisis Ai sin pasar la discretizacion D. De tal forma se evita que cualquier proceso
de regeneracion o adaptacion de malla.
El objetivo de disminuir la cantidad y complejidad de operaciones que realiza el algorit-
mo, junto con la anulacion de procesos de discretizacion para cada avance en el crecimiento
de la grieta, podrıa significar la disminucion del costo computacional.
2.4. Modelo implementado
Dado un solido discretizado con elementos triangulares lineales convencionales, se genera una
nueva malla en la cual las caras compartidas de los elementos triangulares son separadas una
pequena distancia. Los nodos de caras compartidas por triangulos son duplicadas y luego
unidas con elementos tipo vınculo rıgido, ver figura 2-4. En esta tesis se le da el nombre de
explosion de malla a dicho procedimiento y se explica a profundidad en el capıtulo 3.
2.4 Modelo implementado 21
Figura 2-4: Ejemplo de explosion de malla.
Posterior a la explosion de la malla, se realiza el calculo de esfuerzos, se identifica la lo-
calizacion del maximo esfuerzo principal y su direccion. La grieta crece perpendicular a la
direccion del esfuerzo principal y avanza por entre los elementos triangulares planos cuando
la rigidez del los elementos tipo vınculo de hace cero, ver seccion 2.2.
Para realizar un simulacion de trayectoria de fisuracion, se realizan los siguientes pasos:
i Discretizacion inicial del problema.
ii Explosion de la malla.
iii Incremento de la accion aplicada.
iv Calculo de esfuerzos nodales.
v Calculo de esfuerzos principales nodales.
vi Seleccion del elemento vınculo asociado a la trayectoria de fisura.
vii Proceso incremental de carga.
Los pasos i y ii, se explican de manera detallada en el capitulo 3. A continuacion se presenta
el desarrollo de los pasos subsiguientes.
Incremento de la accion aplicada.
La carga es aplicada al solido de forma discreta. Se producen deformaciones que son calcula-
das a nivel elemental con (2-28). De igual manera los esfuerzos en el interior de cada elemento
son calculados mediante (2-29). En particular, el elemento triangular lineal ofrece el mismo
valor de esfuerzo en cualquier punto de su interior. Los calculos realizados se realizan en la
22 2 Implementacion del Metodo de los elementos finitos con explosion de malla.
malla original, figura 3-1.
ε(e) = [B]u(e) (2-28)
σ(e) = [D]ε(e) (2-29)
Calculo de esfuerzos nodales.
Recordando que el campo de esfuerzos y deformaciones es continuo en el intermedio de los
elementos pero discontinuo en los nodos, se calcula el promedio de la misma componente de
esfuerzo obtenida en todos los elementos que comparten un mismo nodo. El resultado sera
una componente de esfuerzo en cada nodo de la malla (figura 2-5).
Figura 2-5: Calculo del esfuerzo nodal promedio.
Calculo de esfuerzos principales nodales.
A partir del esfuerzo promedio nodal de cada componente σx−n, σy−n y σxy−n, se realiza
el calculo del primer y segundo esfuerzo principal, σmax−n y σmin−n respectivamente para
cada nodo, Ademas se calcula el angulo θσmax−n que indica la direccion del esfuerzo principal
maximo. Para lo anterior se utilizan las siguientes expresiones.
σmax−n =
√(σx−n − σy−n
2)2 + σxy−n +
(σx−n + σy−n
2
)
σmin−n =
√(σx−n − σy−n
2)2 + σxy−n −
(σx−n + σy−n
2
)(2-30)
θσmax−n =1
2arctan
(2σxy−n
σx−n − σy−n
)Una vez obtenido el esfuerzo principal maximo en cada nodo, se busca el valor mas alto en
la malla σmax−p, siendo p el nodo donde el esfuerzo principal maximo es el mayor de los
nodos de la malla. Este valor es comparado con la resistencia a tension del material σu, de
tal manera que si se cumple la condicion σmax−p > σut, ocurre fractura del material como lo
establece el criterio de falla de Rankine.
2.4 Modelo implementado 23
Seleccion del elemento vınculo asociado a la trayectoria de fisura.
El nodo con σmax−p es por donde se propaga la fisura. Se identifican los elementos vınculo
que surgieron en la explosion de la malla para este nodo. Ver figura 2-6.
Figura 2-6: Elementos barra asociados al nodo con el maximo esfuerzo principal.
El angulo de inclinacion de cada vınculo asociado al nodo de esfuerzo maximo se calcula
mediante la siguiente expresion:
δx = ‖xi − xj‖ δy = ‖yi − yj‖ θlink = arctan
(δyδx
)(2-31)
Los angulos de inclinacion de los vınculos asociados θlink se comparan con el angulo de incli-
nacion del esfuerzo maximo θσmax−p calculado en (2-30). Esta comparacion se realiza con el
objetivo anular la rigidez del elemento vınculo cuya inclinacion es mas cercana a la direccion
del esfuerzo maximo. Este modelo supone que la fisura tiende a propagarse en direccion
perpendicular al esfuerzo principal maximo.
En el caso de la figura 2.7(a), el vector que representa la direccion del esfuerzo principal
maximo del nodo p σmax−p, muestra que su inclinacion es mas cercana a la del elemento
vınculo conformado por los nodos 15 y 11. De acuerdo a lo anterior, el elemento entre for-
mado por los nodos 15 y 11 soporta mayor carga a traccion y en consecuencia su rigidez sera
anulada para representar la fisura.
24 2 Implementacion del Metodo de los elementos finitos con explosion de malla.
(a) Link elegido (b) Error de angulo
Figura 2-7: (a) Angulo de inclinacion del esfuerzo principal maximos θσmax Y angulos de
inclinacion de los elementos barra asociados al nodo p θilink. (b) Error de angulo
en el crecimiento de grieta.
Error del angulo de propagacion de la fisura.
Dado que la inclinacion del elemento vınculo desactivado difiere del valor exacto de pro-
pagacion θσmax−p , la direccion de propagacion de la fisura lleva un error implıcito θerr. Se
define dicho error como la diferencia entre θσmax−p y θlink del vınculo seleccionado, ver figura
2.7(b).
θerr = θσmax−p − θlink (2-32)
Proceso incremental de carga.
El procedimiento anteriormente indicado se repiten en cada paso de carga. En particular, la
rigidez del vınculo asociado a la direccion de la fisura se anula para el siguiente paso de carga.
El angulo calculado en el siguiente paso de carga con (2−31) se corrige sumando la magnitud
θi−1err , la cual corresponde al error angular obtenido en el paso anterior:
θiσmax = θi−1σmax + θi−1
err (2-33)
3 Explosion de malla
El metodo de los elementos ha sido muy bien adaptado a la mecanica de fractura permi-
tiendo la resolucion de problemas complejos. Sin embargo, la modelacion de crecimiento de
fisuras a traves de una malla de elementos finitos es de gran dificultad, pues cada avance en
el crecimiento de la grieta representa una modificacion topologica de la malla [?]. De manera
general, la representacion numerica de la fisura se hace mediante una discontinuidad en el
campo de desplazamientos del dominio. A partir de dicha interpretacion, se han generado
diferentes modelos numericos, cuya revision fue presentada en la seccion 1.2.
En esta investigacion, la singularidad en el campo de desplazamiento y la propagacion de
la fisura se hace mediante una explosion de malla. La simulacion de crecimiento de fisura
inicia con la discretizacion inicial del problema a traves de un enmallado convencional 2D
de elementos finitos. La malla inicial es sometida a un procedimiento que genera una nueva
discretizacion, al que se llamo explosion de malla en esta tesis. El objetivo de este capıtulo
es exponer la metodologıa de explosion de malla, para esto se define el procedimiento y a su
vez se expone el algoritmo desarrollado para generar la nueva discretizacion.
3.1. Discretizacion inicial del problema.
Se genera la geometrıa del especimen, se establecen sus condiciones de borde y se enmalla con
elementos triangulares convencionales de 3 nodos. Si las acciones aplicadas sobre el solido
generan una distribucion uniforme del estado de esfuerzos, cualquier punto material podrıa
ser el comienzo de una fisura. Con el fin de inducir la formacion de una fisura, se establece
una entalla en el lugar del solido donde se quiere iniciar el proceso de fractura. Se obtiene
una geometrıa como la que se muestra en la figura 3-1.
26 3 Explosion de malla
Figura 3-1: Discretizacion inicial del problema.
3.2. Definicion de explosion de malla.
Dado un solido discretizado con elementos triangulares lineales convencionales1, se genera
una nueva malla en la cual las caras compartidas de los elementos triangulares son separadas
una pequena distancia. Los nodos de caras compartidas por los elementos triangulares son
duplicados y luego unidos con elementos tipo vınculo rıgido. Para este procedimiento se
adopto el termino explosion de malla2, ver figura 3-2.
Figura 3-2: Malla explotada que genera el codigo.
3.3. Funcionamiento del programa de explosion de malla.
El codigo presentado en este documento realiza explosion de mallas en 2D generadas median-
te elementos triangulares convencionales planos. El codigo parte de mallas como la mostrada
1Se denominara en adelante como malla original2La discretizacion generada con este procedimiento se denomina malla explotada
3.3 Funcionamiento del programa de explosion de malla. 27
en la figura 3-2 (a) y las convierte en mallas como la presentada en la figura 3-2 (b).
3.3.1. Archivos de entrada.
Conect.txt: Contiene los parametros iniciales para dimensionar el problema. Allı se al-
macenan datos iniciales importantes como el numero de nodos, el numero de elementos, la
cantidad de nodos por elemento y la definicion de la frontera del dominio. Todos estos datos
corresponden a la malla original.
Nodos.inp: Contiene las coordenadas de los nodos de la malla original.
Conectividades.inp: Contiene las conectividades de la malla original. Esto es la nume-
racion de los nodos que conforman cada uno de los elementos de la malla original.
Analisis.inp: Contiene los grupos de nodos que conforman cada una de las fronteras del
dominio. Adicionalmente, contiene la informacion de las propiedades mecanicas del solido.
3.3.2. Archivos de salida.
NodosNew.inp: Almacena las coordenadas de los nodos de la explotada.
ConectividadesNew.inp: Contiene las conectividades de la malla explotada.
s1.inp, s2.inp, s3.inp: Se almacenan los arreglos S1, S2 y S3 respectivamente. Son arreglos
generados por el programa que contienen informacion tanto de la malla original como de la
malla explotada. En la tabla 3-1 la informacion de estos arreglos es ampliada.
contorno.inp: Este archivo contiene los grupos de nodos (o elementos) que resultan en
el contorno de la nueva malla.
3.3.3. Variables y arreglos utilizados.
En la tabla 3-1 se presenta el nombre de la variables y arreglos utilizados por el programa.
Tambien para cada una de ellos se presenta su tamano, tipo y una breve descripcion de la
funcion que cumple dentro del programa de explosion de malla.
28 3 Explosion de malla
Nombre de la Tamano Tipo Descripcion
Variable/Arreglo
NUMNODE - I Numero de nodos
de la malla original.
NELEMS - I Numero de elementos
de la malla original.
nnod - I Numero de nodos
por elemento.
Searstr - L Variable logica para
determinar si se
realiza o no la
lectura de datos.
dim - I Dimension del espacio
del dominio (1D, 2D, 3D).
conectividades (NELEMS,nnod+1) I Conectividades de
cada elemento
de la malla original.
nodes (NUMNODE,dim) R Coordenadas de cada nodo
de la malla original.
S1 (NUMNODE) I Cantidad de elementos
triangulares que confluyen
en cada nodo de la malla.
original.
S2 (NUMNODE,12) I Numero de cada uno
de los elementos
triangulares que confluyen
en un nodo de la malla
original.
Tabla 3-1: Tabla de variables y arreglos definidos en el programa de explosion de malla.
3.3 Funcionamiento del programa de explosion de malla. 29
Nombre de la Tamano Tipo Descripcion
Variable/Arreglo
S3 (NUMNODE,12) I Nueva numeracion de
los nodos.
Cada fila contiene
el numero de cada
nuevo nodo que se
formo a partir
del nodo explotado
en la malla original.
S4 (NUMNODE) I Cantidad de elementos
link nuevos que
surgen de cada uno
de los nodos explotados.
S5 (NUMNODE,12) I Numeracion de cada link
que surgio de cada nodo
despues de la explosion.
sumaNodes - I Suma total de nodos de
la nueva malla.
nodesNew (sumaNodes,dim) R Coordenadas de los
nuevos nodos creados
en la malla explotada.
conectividadesNEw2D (NELEMS,nnod+1) I Conectividades de
la malla explotada.
link - I Numero de elemento barrar a adicionar
en la malla explotada.
conectividades1Dlink (20*NUMNODE,2) I Conectividades de los links
de la malla explotada.
h - R Diametro equivalente
de cada elemento.
Tabla 3-2: Tabla de variable y arreglos definidos en el programa. Continuacion.
30 3 Explosion de malla
Nombre de la Tamano Tipo Descripcion
Variable/Arreglo
filasContorno(x) - I Numero de filas
del arreglo de nodos,
grupos o elementos.
(Archivo .inp)
Contorno(x) (filasContorno(x),6) I Arreglo que contiene grupo,
nodos o elementos.
Tabla 3-3: Tabla de variable y arreglos definidos en el programa. Continuacion.
3.3.4. Descripcion del programa.
Lectura de archivos
En esta parte del programa se realiza la lectura de los archivos de datos de entrada (.inp).
Se utiliza la funcion logica Searstr que hace la busqueda de palabras claves del archivo de
entrada. A partir de esa palabra clave se realiza la lectura de un arreglo.
A manera de ejemplo, en el caso de las primeras lıneas, se abre el archivo (con estatus
“old”) llamado conectividades.inp, luego se busca la cadena de caracteres clave: *ELE-
MENT,TYPE=U1,ELSET=UEL., en este caso dicha cadena de caracteres corresponde a
las conectividades de la malla original. De manera que se leen las conectividades de la malla
original, siendo NELEMS el numero de filas y nnod+1 en numero de columnas.
La parte del programa antes mencionada se muestra a continuacion, donde se realiza la
lectura de conectividades y de coordenadas de los nodos de la malla original:
! Lectura del archivo de conectividades
open(UNIT=1,file=’conectividades.inp’, status=’old’)
if (Searstr (1,’*ELEMENT,TYPE=U1,ELSET=UEL’)) then
READ(1,*)((conectividades(i,j),j=1,nnod + 1),i=1,NELEMS)
else
stop ’###..Error en lectura’
end if
close(1)
!
! Lectura del archivo de nodos y sus coordenadas
open(UNIT=2,file=’nodos.inp’, status=’old’)
if (Searstr (2,’*NODE,NSET=N2’)) then
3.3 Funcionamiento del programa de explosion de malla. 31
READ(2,*) (i,(nodes(i,j),j=1,2),i=1,NUMNODE)
else
stop ’###..Error en lectura’
end if
close(2)
Determinacion de elementos confluyentes en un nodo.
En esta parte del programa, se determina la cantidad y numeracion de cada uno de los ele-
mentos triangulares que llegan a un nodo en particular. Para esto, se recorre la matriz de
conectividades elemento por elemento. El objetivo es saber cuantas veces aparece un nodo
especıfico en dicha matriz. Este numero de veces que aparece el nodo i-esimo en la matriz
de conectividades es equivalente al numero de elementos triangulares que llegan a ese mismo
nodo. El dato obtenido es almacenado en el arreglo S1.
A la vez que se realiza conteo de las veces que aparece un nodo cualquiera en la matriz
de conectividades, tambien se identifica el numero de elemento que llega a un nodo especıfi-
co. Este dato de los elementos que confluyen en cada uno de los nodos se almacena en el
arreglo S2.
A manera de ejemplo, se define una malla triangular como la de la figura 3-3. Esta es
la malla original.
Figura 3-3: Malla convencional con numeracion.
De esta manera se tendrıa un arreglo (3-1) llamado nodes donde se almacenan las coorde-
nadas de cada uno de los nodos de la malla original. La primera columna corresponde a la
numeracion del nodo y las siguientes segunda y tercera, corresponden a la coordenada en X
e Y respectivamente.
32 3 Explosion de malla
nodes =
1 x1 y1
2 x2 y2
3 x3 y3
4 x4 y4
5 x5 y5
6 x6 y6
(3-1)
El arreglo (3-2) llamado conectividades contiene la informacion de los nodos que conforman
cada uno de los elementos de la malla original. La primera columna corresponde al numero de
elemento y las siguientes a cada uno de los nodos que conforman el elemento correspondiente.
conectividades =
1 1 6 5
2 6 4 5
3 2 6 1
4 3 4 6
5 2 3 6
(3-2)
El arreglo S1 para el ejemplo tratado se muestran en (3-3). La primera fila del mismo co-
rresponde a la cantidad de elementos triangulares que confluyen en el nodo 1. La segunda
fila a la cantidad de nodos que confluyen en el nodo 2 y ası sucesivamente.
S1 =[2 2 2 2 2 5
](3-3)
El arreglo S2, mostrado en (3-4), muestra en su primera fila los elementos que confluyen en
el nodo 1. En la segunda fila a los nodos que confluyen en el nodo 2 y ası sucesivamente.
S2 =
1 3 0 0 0
3 5 0 0 0
5 4 0 0 0
4 2 0 0 0
2 1 0 0 0
1 2 3 4 5
(3-4)
Conteo del total de nodos y nueva numeracion de los mismos.
Como se mostro en la figura 3-2, la explosion de la malla busca que cada uno de los elemen-
tos triangulares tenga nodos independientes, es decir, que cada triangulo tenga sus propios
lados y que esos lados no se compartan con ningun otro triangulo. Para lograr esto, primero
es necesario conocer la cantidad total de nodos que debera tener en la malla explotada. Este
valor total de nodos se puede encontrar sumando la cantidad de elementos que confluyen
3.3 Funcionamiento del programa de explosion de malla. 33
en cada nodo, pues es sabido que si un elemento comparte un lado con otro elemento, se
deben generar dos nuevos nodos. Por ejemplo, en la malla de la figura 3-3, se ve que los
elementos 3 y 5 comparten un lado del triangulo, o sea, comparten los nodos 2 y 6 . De
manera que tanto el nodo 2 , como el nodo 6 , quedaran convertidos en 4 nuevos nodos.
Es decir, cada nodo de la malla original explotara en un numero de nodos igual a la cantidad
de elementos que confluyan en el.
A continuacion se muestra la parte del programa en donde se realiza dicho conteo para
determinar el nuevo total de nodos de la malla explotada:
! Se totaliza el numero de nodos que incluyen los antiguos y los nuevos
sumaNodes = 0
do i=1,NUMNODE
sumaNodes = sumaNodes + s1(i)
enddo
!
El nuevo total de nodos se almacena en sumaNodes. Este valor sirve para realizar la nueva
numeracion de nodos, figura 3-4, mediante las siguientes lıneas de codigo.
! Numeracion de nuevos nodos en la nueva base
allocate(s3(NUMNODE,20))
s3 = 0
k = 1
do i=1,NUMNODE
do j=1,s1(i)
s3(i,j) = k
k = k + 1
enddo
enddo
!
Figura 3-4: Nueva numeracion de los nodos.
34 3 Explosion de malla
La nueva numeracion de los nodos se almacena en el arreglo S3. En este arreglo, cada fila se
contiene los nodos generados a partir del nodo antiguo. Es decir, en la primera fila quedan
almacenados los nodos generados a partir del nodo 1 de la malla original, en la segunda fila
los generados a partir del nodo 2 de la malla original y ası sucesivamente para todos los
nodos de la malla original. El arreglo S3 para la malla mostrada en la figura 3-4 se muestra
en (3-5).
S3 =
1 2 0 0 0
3 4 0 0 0
5 6 0 0 0
7 8 0 0 0
9 10 0 0 0
11 12 13 14 15
(3-5)
Asignacion de coordenadas a los nodos de la malla explotada.
Posterior a realizar la nueva numeracion de los nodos, se asignan las coordenadas a cada no-
do de la malla explotada. En este punto del codigo, las coordenadas asignadas a los nuevos
nodos, son iguales a las del nodo que se exploto, aunque en figura 3-4 se muestren los nodos
de manera separada. Por ejemplo el nodo 6 de la malla original, figura 3-3, exploto en 5
nuevos nodos, a saber, 11 , 12 , 13 , 14 y 15 . En la figura 3-4 se ven separados una
pequena distancia, sin embargo, estos 5 nuevos nodos tienen las mismas coordenadas entre
sı e iguales a las coordenadas del nodo 6 de la malla original. La separacion de los nodos,
se hace posteriormente en el codigo.
Nuevas conectividades.
Obtenida la numeracion de los nuevos nodos surgidos en la explosion de la malla, es necesario
realizar las conectividades entre ellos con el fin definir los elementos triangulares. Aunque
las conectividades varıen, es importante resaltar que la geometrıa de la malla original no
cambia, lo unico que ha cambiado hasta el momento es la numeracion de la misma.
En las siguientes lıneas del codigo se realizan las nuevas conectividades.
! Construccion de las conectividades: Nota: Inicialmente los triangulos no cambian!
do i=1,NELEMS
conectividadesnew2D(i,:)=i
enddo
3.3 Funcionamiento del programa de explosion de malla. 35
!
p = 1
do i=1,NUMNODE
do j=1,s1(i)
do k=1,nnod
if(conectividades(s2(i,j),k+1).eq.i)then
conectividadesnew2D(s2(i,j),k+1) = p
endif
enddo
p = p + 1
enddo
enddo
Hasta el momento la malla original de la figura 3-3 se ha convertido en la malla que se
muestra en la figura 3-5. Las nuevas conectividades de los elementos triangulares se muestra
en el arreglo conectividadesNew2D que se muestra en (3-6).
Figura 3-5: Nuevas conectividades.
conectividadesNew2D =
1 1 11 10
2 12 8 9
3 13 15 2
4 6 7 13
5 4 5 14
(3-6)
Construccion de los elementos vınculo rıgido.
En esta parte del programa se crean los elementos link que unen los nuevos nodos. Para
describir este procedimiento se toma como ejemplo nuevamente la malla de la figura 3-3.
36 3 Explosion de malla
La operacion incluye los siguientes pasos:
1. Buscar en S2 cuales son los elementos que confluyen en un nodo i-esimo. Como ejemplo
tomaremos la fila 4 del arreglo S2, la cual corresponde al nodo 4 de la malla original.
De aca se obtiene que los elementos que confluyen en el nodo 4 son los elementos 4
y 2 .
2. Con los datos del paso anterior, el programa va al arreglo de conectividades de la malla
original y se ubica en las filas correspondientes a los elementos encontrados en S2. En
este caso el programa ubica las filas 4 y 2 del arreglo de conectividades de la malla
original. Se obtienen los nodos correspondientes a dichos elementos, es decir, para el
elemento 4 se tiene conectividades(4, :) =[3 4 6
]y para el elemento 2 se tiene
conectividades(2, :) =[6 4 5
].
3. Se toma cada nodo del primer elemento y se compara con cada uno de los nodos
del otro elemento. Para el ejemplo, el primer nodo del elemento 4 es el nodo 3 ,
es decir conectividades(4, 2) = 3. Este ultimo nodo se compara con los nodos 6 ,
4 y 5 porque conectividades(2, :) =[6 4 5
]. Posteriormente se toma el segundo
nodo del elemento 4 conectividades(4, 2) = 4 con los nodos 6 , 4 y 5 porque
conectividades(2, :) =[6 4 5
]. Se realiza la misma operacion para el tercer nodo del
elemento 4 . El programa va contando las veces en que esos nodos comparados son
iguales. Cuando este numero es 2 (r = 2 en el programa), se sabe que los elementos
4 y 2 esta compartiendo 2 nodos, esto quiere decir que los elementos comparten uno
de sus lados. En el caso del ejemplo tratado, el lado compartido por los elementos 4
y 2 son los nodos 6 y 4 de la malla original.
4. Conocidos lados que se comparten entre elementos triangulares y los nodos que corres-
ponden a ese lado compartido, el programa va al arreglo S3 y determina en que nodos
de la malla explotada se convirtieron los nodos del paso anterior. En el caso puntual
del ejemplo en discusion, se puede ver que el nodo 6 se convirtio en los nodos 11 ,
12 , 13 , 14 y 5 porque S3(6, :) =[11 12 13 14 15
]y el nodo 4 se convirtio
en los nodos 7 y 8 porque s3(4, :) =[7 8
].
5. Se crea un elemento vınculo rıgido entre esos 2 nodos del lado compartido. En el caso
del ejemplo, se crearıan 2 vınculos. Uno entre el nodo 7 y 8 de la malla explotada
y otro entre los nodos 12 y 13 . Es importante notar que aunque el nodo 6 de la
malla original exploto en 5 nodos diferentes, el programa solo crea el vınculo entre los
nodos 12 y 13 . Esto se debe a que el programa tiene en cuenta los lados compartidos
por los elementos nada mas. De esta manera se evita que hayan elementos link que se
intercepten con otros, es decir, el programa evita por ejemplo que se cree un link entre
los nodos 14 y 12 o entre los nodos 15 y 13 de la malla explotada.
3.3 Funcionamiento del programa de explosion de malla. 37
6. El programa hace una sumatoria de elementos link creados y crea un arreglo de conec-
tividades llamado conectividades1Dlink donde se almacenan las conectividades de los
nodos link.
A continuacion se muestra la parte del programa encargada de crear los elementos vınculo
rıgido.
! Se construyen los elementos link
! Se adicionan elementos link, aumenta el numero de elementos
link = 0
do i=1,NUMNODE
nlinknod = 0
do j = 1,s1(i)-1
do k = j+1,s1(i)
r = 0
do l=1,3
do m=1,3
if(conectividades(s2(i,j),l+1).eq.conectividades(s2(i,k),m+1))then
r = r+1
endif
enddo
enddo
if(r.eq.2)then
link = link + 1
nlinknod = nlinknod + 1
conectividades1Dlink(link,1) = s3(i,j)
conectividades1Dlink(link,2) = s3(i,k)
s5(i,nlinknod) = link + NELEMS
endif
enddo
enddo
s4(i) = nlinknod
enddo
pause
! Numero total de links
NtotalLinks = link
!
En la figura 3-6 se muestran los vınculo rıgido creados, junto con su numeracion.
38 3 Explosion de malla
Figura 3-6: Malla numerada despues de la explosion.
Asignacion de la longitud a los elementos vınculo rıgido.
En la figura 3-6 separacion existente entre las caras de los elementos triangulares que antes
de las explosion eran compartidas. Para darle una longitud a los elementos vınculo, es nece-
sario generar un espacio entre los elementos triangulares. El espacio generado garantiza que
el dominio del problema no sufra cambios en su geometrıa ni en sus dimensiones. La forma
de generar el espacio entre elementos triangulares garantizando que el dominio no varıe, es
haciendo que los elementos triangulares disminuyan su tamano pero no su forma.
El programa realiza la disminucion de tamano de los triangulos a partir del baricentro y
de las bisectrices de los mismos. De manera que los vertices del triangulo se mueven sobre
las bisectrices y en direccion al baricentro. A continuacion se muestran las lıneas del progra-
ma que se encargan de este proceso:
! Explosion de los elementos para obtener elementos link con longitud
do i=1,NELEMS
xbar = 0.d0; ybar = 0.d0; xx = 0.d0; dir = 0.d0; x1 = 0.d0; x2 = 0.d0
y1 = 0.d0; y2 = 0.d0; area = 0.d0; h = 0.d0
! se calculara el baricentro
xbar = (1.0/3.d0)*(nodes(conectividades(i,2),1) + &
nodes(conectividades(i,3),1) + &
nodes(conectividades(i,4),1))
ybar = (1.0/3.d0)*(nodes(conectividades(i,2),2) + &
nodes(conectividades(i,3),2) + &
nodes(conectividades(i,4),2))
! Se calcula la direccion de cada nodo hacia el baricentro
xx(1,1) = xbar - nodes(conectividades(i,2),1)
3.3 Funcionamiento del programa de explosion de malla. 39
xx(2,1) = xbar - nodes(conectividades(i,3),1)
xx(3,1) = xbar - nodes(conectividades(i,4),1)
xx(1,2) = ybar - nodes(conectividades(i,2),2)
xx(2,2) = ybar - nodes(conectividades(i,3),2)
xx(3,2) = ybar - nodes(conectividades(i,4),2)
dir(1,:) = xx(1,:)/sqrt(dot_product(xx(1,:),xx(1,:)))
dir(2,:) = xx(2,:)/sqrt(dot_product(xx(2,:),xx(2,:)))
dir(3,:) = xx(3,:)/sqrt(dot_product(xx(3,:),xx(3,:)))
!
! Calculo de lados para hallar el area
x1 = xx(2,1) - xx(1,1)
x2 = xx(3,1) - xx(1,1)
y1 = xx(2,2) - xx(1,2)
y2 = xx(3,2) - xx(1,2)
!
! Area del elemento (para calcular el tamanho caracteristico)
area = 0.5d0*abs(x1*y2-y1*x2)
! tamano caracteristico
h = sqrt(area/3.141592654)
!
! Se reubican los nodos de cada elemento
NodesNew(conectividadesnew2D(i,2),1) = NodesNew(conectividadesnew2D(i,2),1) +...
+ 0.1*h*dir(1,1)
NodesNew(conectividadesnew2D(i,2),2) = NodesNew(conectividadesnew2D(i,2),2) +...
+ 0.1*h*dir(1,2)
NodesNew(conectividadesnew2D(i,3),1) = NodesNew(conectividadesnew2D(i,3),1) +...
+ 0.1*h*dir(2,1)
NodesNew(conectividadesnew2D(i,3),2) = NodesNew(conectividadesnew2D(i,3),2) +...
+ 0.1*h*dir(2,2)
NodesNew(conectividadesnew2D(i,4),1) = NodesNew(conectividadesnew2D(i,4),1) +...
+0.1*h*dir(3,1)
NodesNew(conectividadesnew2D(i,4),2) = NodesNew(conectividadesnew2D(i,4),2) +...
+ 0.1*h*dir(3,2)
enddo
En la figura 3-7 se muestra el baricentro y las bisectrices para cada triangulo y las nuevas
ubicaciones de los nodos de los triangulos.
40 3 Explosion de malla
Figura 3-7: a) Baricentro y bisectrices de los elementos triangulares. b) Nueva posicion de
los vertices del triangulo.
Posterior a esta reubicacion de los nodos, se trazan los nuevos elementos triangulares y
se ubican los vınculos rıgidos correspondientes. En la figura 3-8 se la malla despues de la
explosion.
Figura 3-8: a) Malla original. b) Malla despues de todo el proceso.
3.3 Funcionamiento del programa de explosion de malla. 41
Conteo total de elementos.
En las siguientes lıneas del programa, se realiza la sumatoria total tanto de los elementos
planos como de los elemento lineales.
do i=1,link
print*,i,conectividades1Dlink(i,1),conectividades1Dlink(i,2)
enddo
pause
!
! Numero total de elementos
NtotalElems = NELEMS + NtotalLinks
print*,NtotalElems
pause
4 Ejemplos Numericos de aplicacion
Para validar el metodo utilizado en esta tesis se resolvieron tres ejemplos. Los dos primeros
corresponden a vigas con entallas iniciales sometidas a una carga que genera un estado de
esfuerzos combinados en el frente de la entalla; en cada ejemplo se desarrollan dos casos en
los que se varıa la geometrıa de la entalla inicial. El tercer ejemplo es un benchmark conocido
como Viga de Iosipescu. Los resultados obtenidos fueron comparados con trabajos numericos
y experimentales de otros autores.
4.1. Ejemplo 1: Viga de tres puntos con entalla
excentrica.
La figura 4-1 muestra el problema para el primer ejemplo. Se trata de una viga de (PMMA)
Polimetilmetacrilato simplemente apoyada y cargada en la parte superior central. Tiene una
entalla que no esta alineada con la carga aplicada, de manera que se genera un modo mixto
de cargas en la punta de la entalla, lo que produce una trayectoria de fisuracion curvilınea
[36] (figura 4-1). Se realizaron simulaciones para dos casos donde los valores de a y b se
variaron como se muestra en el cuadro 4-1; esto permite observar trayectorias de fisuracion
diferentes. Los parametros de los casos I y II son dados por la referencia [36]: Modulo de
Young E = 3100MPa, Relacion de Poisson ν = 0,4, Resistencia a la traccion σu = 76MPa
y un valor de carga aplicada P que garantiza la fractura.
Figura 4-1: Geometrıa de la viga de tres puntos con entalla excentrica.
44 4 Ejemplos Numericos de aplicacion
Caso a [cm] b [cm]
I 6.0 1.0
II 5.0 1.5
Tabla 4-1: Variacion de la distancia a y b en la figura 4-1.
Ejemplo 1 - Caso I.
Se obtuvieron trayectorias de fisuracion para tres mallas diferentes definidas en el cuadro
4-2. En la figura 4-2, se muestran las trayectoras de fisuracion obtenidas con cada tipo de
malla. Los resultados son comparados con la trayectoria obtenida experimentalmente por
Ingraffea y otros [36]. La malla 1 muestra una fisura cercana a la experimental; se observa
que los valores mas alejados se presentan al principio y al final de la trayectoria. Dada la
baja cantidad de elementos de la malla 1, se ve una lınea de grieta con trazos mas rectos
que la curva experimental.
Malla No. Elementos No. Nodos
1 1413 765
2 4412 2310
3 9432 4869
Tabla 4-2: Mallas utilizadas en ejemplo 1 - Caso I.
(a) Malla 1 (b) Malla 2 (c) Malla 3
Figura 4-2: Resultados ejemplo 1 - Caso I. Simulacion con mallas diferentes.
4.1 Ejemplo 1: Viga de tres puntos con entalla excentrica. 45
El resultado obtenido con la malla 2 muestra valores alejados al inicio de la fisura, sin em-
bargo en este caso la trayectoria se corrige en la etapa final de la propagacion. Al considerar
la lınea de tendencia de este resultado, se observa una buena aproximacion a la trayecto-
ria de fisuracion experimental donde se corrobora una trayectoria curvilınea. La malla mas
fina muestra una trayectoria de fisura muy cercana al resultado experimental. Los valores
obtenidos permanecen cercanos tanto al inicio como al final de la fisura. Adicionalmente se
observa que la lınea de tendencia obtenida con la malla 3 presenta bastante similitud con
la trayectoria experimental; en especial para la primera mitad del recorrido las lıneas estan
sobrepuestas.
Ejemplo 1- Caso II.
El caso II del ejemplo 1 considera un cambio de ubicacion y profundidad de la entalla inicial
como lo indica el cuadro 4-1. En la figura 4.3(a) se muestran los resultados obtenidos para
la malla mas fina. Se observa que la lınea de tendencia de la trayectoria corresponde en gran
medida a la lınea de trayectoria experimental obtenida en [36]. En la figura 4.3(b) se ilustra
la trayectoria de fisura obtenida mediante la simulacion en la malla de 9636 elementos.
(a) Comparacion (b) Simulacion
Figura 4-3: Resultados ejemplo 1 - Caso II. Simulacion con malla Fina.
46 4 Ejemplos Numericos de aplicacion
4.1.1. Ejemplo 2: Viga de tres puntos con entalla y agujeros
excentricos.
El segundo ejemplo realizado se presenta en la figura 4-4. Es similar al ejemplo anterior,
pero esta vez el especimen tiene tres agujeros internos. La presencia de agujeros en la viga
distorsionan los campos de esfuerzo y de deformacion, lo que proporciona trayectorias de fi-
suracion curvilıneas bastante caracterısticas [6]. De igual manera se realizan dos casos donde
se varıan los valores de a y b de la figura 4-4, (cuadro 4-3). Este es un problema que ha
sido estudiado experimentalmente por Ingraffea y colaboradores [36] y numericamente en las
referencias [6], [37], [38], [39], [40], entre otros.
Figura 4-4: Geometrıa viga de tres puntos con entalla y agujeros excentricos.
Caso a [cm] b [cm]
I 6.0 1.0
II 5.0 1.5
Tabla 4-3: Variacion de la distancia a y b en la figura 4-4.
Ejemplo 2 - Caso I.
Se obtuvieron trayectorias de fisuracion para tres mallas diferentes, (cuadro 4-4). Los resulta-
dos se presentan en la figura 4-5(a), 4-5(b) y4-5(c) para las mallas 1, 2 y 3 respectivamente.
En la figura 4-5(d) se presenta la imagen de los resultados experimentales obtenidos por
Ingraffea y colaboradores [36].
Para las mallas 1 y 2 se obtuvieron trayectorias con tramos rectilıneos. Los resultados presen-
tan una tendencia similar a la fisura experimental, sin embargo los valores obtenidos al inicio
4.1 Ejemplo 1: Viga de tres puntos con entalla excentrica. 47
y al final de la propagacion, difieren de los resultados experimentales. La lınea de tendencia
esta bastante distanciada de la trayectoria real.
Para la malla 3, los resultados presentan una fuerte tendencia a seguir la trayectoria de
la fisura real. En este caso los valores obtenidos siempre estan cercanos a la trayectoria ex-
perimental. La lınea de tendencia se sobrepone a la trayectoria experimental. En la figura
4-6 se muestra la secuencia de propagacion de grieta para la malla mencionada.
Malla No. Elementos No. Nodos
1 2630 1398
2 10520 5428
3 23200 11846
Tabla 4-4: Mallas utilizadas en el Ejemplo 2 - Caso I.
(a) Malla 1 (b) Malla 2 (c) Malla 3
Figura 4-5: Resultados ejemplo 2 - Caso I. Simulacion con mallas diferentes.
48 4 Ejemplos Numericos de aplicacion
(a) Paso Inicial (b) Paso 10 (c) Paso 30 (d) Paso Final
Figura 4-6: Resultados ejemplo 2 - Case I. Evolucion de la simulacion de crecimiento de
fisura.
Ejemplo 2 - Caso II.
Para el caso II del segundo ejemplo, la ubicacion y geometrıa de la entalla inicial se modifican
como se menciono en la tabla 4-3. Se hace con el proposito de verificar trayecotrias de
fisuracion diferentes con respecto al caso anterior. Se obtuvieron trayectorias de fisuracion
para tres mallas con las caracterısticas presentadas en la tabla 4-5.
Malla No. Elementos No. Nodos
1 2191 1172
2 10284 5311
3 16296 8360
Tabla 4-5: Mallas utilizadas en ejemplo 2 - Caso II.
En la figura 4.7(a) se muestran la trayectoria de fisuracion obtenida con la malla 1. Se
observa una diferencia sustancial con la trayectorias de la obtenida experimentalmente en
[36]. Utilizando una malla de mas elementos, se obtiene una trayectoria mas cercana a la
experimental, sin embargo esta solo se acerca a la trayectoria experimental en la primera
mitad como se ve en la figura 4.7(b). Para la malla 3, figura 4.7(c), la fisura obtenida es
mas cercana a la trayectoria experimental reportada por Ingraffea y otros en [36]. Se observa
similitud a lo largo de toda su longitud y el punto de llegada al agujero del especimen es
4.1 Ejemplo 1: Viga de tres puntos con entalla excentrica. 49
muy cercano. En la figura 4.7(c) se muestra en lınea punteada la imagen de la trayectoria
experimental presentada en [36].
(a) Malla1 (b) Malla2 (c) Malla3
Figura 4-7: Resultados ejemplo 2 - Caso II. Simulacion con mallas diferentes.
4.1.2. Ejemplo 3: Viga con entalla y dos cargas.
La figura 4-8 muestra un ejemplo de aplicacion. Esta prueba fue presentada por Arrea y otros
en [3]. Se trata de una viga de concreto con entalla inicial sometida a dos cargas puntuales y
apoyos simples. Este tipo de geometrıa y de aplicacion de las cargas proporciona un esfuerzo
cortante constante y un momento de casi 0 [41]. El modulo de Young es de E = 24800Mpa,
relacion de Poisson de ν = 0,18 y resistencia a la traccion de σt = 4Mpa. El espesor de la
viga es de t = 152mm, las demas medidas indicadas en la figura 4-8 estan en mm.
50 4 Ejemplos Numericos de aplicacion
Figura 4-8: Geometrıa viga con entalla y dos cargas.
Se obtuvo una trayectoria simulada para una malla inicial de 8417 nodos y 16284 elementos.
La comparacion entre los resultados obtenidos con el modelo presentado, los reportados por
Arrea en [3] y la lınea de tendencia de la trayectoria obtenida se observan en la figura 4.9(a).
Se observa que cualitativamente las trayectorias son similares. Sin embargo, la trayectoria
obtenida con el modelo propuesto se aleja de la curva experimental en su primer tramo,
para luego acercarse al punto final de la fisura. Adicional a lo anterior, en la figura 4.9(b),
se muestra una imagen de la simulacion realizada.
(a) Comparacion resultados
ejemplo 3.
(b) Resultado final ejemplo 3.
Figura 4-9: Comparison between presented model and benchmark of Arrea et al. in [3]
5 Conclusiones
La formulacion y validacion de un modelo discreto para prediccion de trayectorias de fisura-
cion de solidos bajo condicion de esfuerzo plano fue presentado en esta tesis. Es un enfoque
cualitativo basado en el FEM bajo una extension de malla1, que permite determinar la tra-
yectoria de agrietamiento de solidos con comportamiento idealmente elastico. Funciona para
problemas en 2D y utiliza una discretizacion con elementos triangulares; se asume que la
fisura se propaga en direccion perpendicular al maximo esfuerzo principal. De los resultados
obtenidos y presentados previamente en el capitulo 4 se obtienen las siguientes conclusiones.
1. El concepto de explosion de malla presentado en el capitulo 3 sugiere una novedad en
la forma como se aborda la discontinuidad geometrica que representa la fisura en el
solido. Su principal ventaja es que la malla de elementos finitos no necesita ajustarse a
geometrıa de la fisura, incluso durante su crecimiento. Esto permite que las simulaciones
de trayectoria de fisuracion puedan hacerse con una unica malla inicial.
2. La metodologıa presentada permite predecir trayectorias de fisuracion sin la necesidad
de hacer calculos de factores de intensidad de esfuerzo, ni realizar algun tipo regene-
racion o adaptacion de malla. Esto permite reducir la cantidad de calculos durante
la simulacion de la trayectoria de fisuracion, permitiendo a su vez reducir el costo
computacional del procedimiento.
3. Las trayectorias de agrietamiento obtenidas, presentan buenas aproximaciones respecto
a los resultados experimentales y numericos reportados por otros autores. Sin embargo,
la desactivacion de elementos vınculo para el avance de la fisura hace que la trayec-
toria obtenida dependa de la topologıa de la malla inicial; el tamano de la malla y
la orientacion de sus elementos afectan de alguna manera la trayectoria de fisuracion
obtenida. Esto se manifiesta principalmente en trayectorias de fisuracion con curvas
pronunciadas.
4. Los resultados podrıan ser mas precisos si la malla tiene elementos mas pequenos,
especialmente en las zonas del solido cercanas a agujeros o con variaciones geometricas
sustanciales. Esto permitirıa obtener trayectorias de fisuracion mas suaves, pues al ser
un metodo que presenta una alta sensibilidad a la malla, el paso de avance de la fisura
queda determinado por el tamano del elemento finito plano.
1Explosion de malla
52 5 Conclusiones
5.1. Trabajos a partir de esta tesis
El trabajo subsiguiente a esta tesis es la programacion de un codigo, que permita obtener
la respuesta (σ − ε) durante el proceso de agrietamiento de los ejemplos tratados. A pesar
de que el codigo realizado determina una buena aproximacion de la trayectoria de grieta,
contar con la respuesta del sistema permitirıa realizar analisis mas profundos en cuanto a la
utilidad del esquema de explosion de malla propuesto.
Por otra parte, podrıan implementarse bechmarks con los codigos generados en esta te-
sis, con el objetivo de comparar los tiempos de solucion y validar de manera cuantitativa la
reduccion del costo computacional durante la simulacion de crecimiento de grietas.
5.2. Productos de la tesis
Artıculo: J. Ochoa-Avendano, D.A. Garzon-Alvarado, Dorian L. Linero. Simplified
qualitative discrete numerical model to determinate cracking pattern in
brittle materials by means of finite element method. Sometido: Junio de 2016
a: Computers and Concrete, An International Journal.
Bibliografıa
[1] A. Ingraffea. Encyclopedia of Computational Mechanics, volume 2 - Solids and Struc-
tures. John Wiley and Sons, 2007.
[2] M. Sabsabi. Modelado de grieta y estimacion de vida en Fretting Fatiga mediante el
Metodo de los Elementos Finitos Extendido X-FEM. PhD thesis, 2010.
[3] Manrique Arrea. Mixed-mode crack propagation in mortar and concrete. Cornell Uni-
versity, Jan., 1982.
[4] B. Coterrell. The past, present and future of fracture mechanics. Engineering of Fracture
Mechanics Journal, 60:533–553, 2002.
[5] Pietro Cornetti, Nicola Pugno, Alberto Carpinteri, and David Taylor. Finite fracture
mechanics: A coupled stress and energy failure criterion. Engineering Fracture Mecha-
nics, 73(14):2021–2033, 2006.
[6] T.N. Bittencourt, P.A. Wawrzynek, A. Ingraffea, and J.L. Sousa. Quasi automatic
simulation of crack propagation for 2d lefm problems. Engineering Fracture Mechanics
Journal, 55:321–324, 1996.
[7] L. P. Pook. Five decades of crack path research. Engineering Fracture Mechanics,
77(11):1619–1630, 2010.
[8] A. A. Griffith. The phenomena of rupture and flow in solids. Philosophical Transactions
of the Royal Society of London A: Mathematical, Physical and Engineering Sciences,
221(582-593):163–198, 1921.
[9] Irwin G.R. Onset of fast crack propagation in high strength steel and aluminum alloys.
Sagamore Research Conference Proceedings, 2:289–305, (1956).
[10] HM Westergaard. Bearing pressures and cracks. SPIE MILESTONE SERIES MS,
137:18–22, 1997.
[11] George R Irwin. Analysis of stresses and strains near the end of a crack traversing a
plate. Spie Milestone series MS, 137(167-170):16, 1997.
54 Bibliografıa
[12] G Anlas, MH Santare, and J Lambros. Numerical calculation of stress intensity factors
in functionally graded materials. International Journal of Fracture, 104(2):131–143,
2000.
[13] Timon Rabczuk. Computational methods for fracture in brittle and quasi-brittle solids:
state-of-the-art review and future perspectives. ISRN Applied Mathematics, 2013, 2013.
[14] J Oliver, AE Huespe, MDG Pulido, and E Chaves. From continuum mechanics to
fracture mechanics: the strong discontinuity approach. Engineering fracture mechanics,
69(2):113–136, 2002.
[15] Milan Jirasek and Borek Patzak. Models for quasibrittle failure: Theoretical and
computational aspects. In Second European Conference on Computational Mecha-
nics,(ECCM), Cracow, Poland, 2001.
[16] M. Jirasek. Conditions of uniqueness for finite elements with embedded cracks. In
Proceedings of the Sixth International Conference on Computational Plasticity.
[17] R.D. Hensehell and K.G. Shaw. Crack tip finite elements are unnecessary. International
Journal for Numerical Methods in Engineering, 9:495–507, 1975.
[18] P.A. Wawrzynek and A. Ingraffea. An interactive approach to local remeshing around
propagating crack. Finite Elements in Analysis and Design, 5:87–96, 1989.
[19] L. Martha, P.A Wawrzynek, and A. Ingraffea. Arbitrary crack propagation using solid
modeling. Engineering with Computers, 9:63–82, 1993.
[20] J. Cervenka. Discrete crack modeling in concrete structures. PhD thesis, 1994.
[21] B.J. Carter, C.S. Chen, A. Ingraffea, and P.A. Wawrzynek. A topology-based system
for modeling 3d crack growth in solid and shell structures. In Proceedings of the Ninth
international Congress on Fracture ICF9, pages 1923–1934.
[22] E.N. Dnorvik, A.M. Cuitino, and G. Giogia. Finite elements with displacement inter-
polated embedded localization lines insensitive to mesh size and distortions. Computer
Methods in Applied Mechanics and Engineering, 90:829–844, 1990.
[23] M. Klisinski, K. Runesson, and S. Sture. Finite element with inner softening band.
Journal of Engineering Mechanics, ASCE, 117:575–587, 1991.
[24] J.C. Simo and J. Oliver. A new approach to the analysis and simulation of strain
softening in solids. In Fracture and Damage in Quasibrittle Structures, pages 25–39.
[25] J.M. Melenk and I. Babuska. The partition of unity finite element method - basic
theory and applications. Computer Methods in Applied Mechanics and Engineering,
139:289–314, 1996.
Bibliografıa 55
[26] C.A Duarte and J.T. Oden. H-p clouds - an h-p meshless method. Numerical Methods
for Partial Differential Equations, 12:673–706, 1996.
[27] N. Sukumar, N. MoA((s, B. Moran, and T. Belytschko. Extended finite element method
for three-dimensional crack modeling. International Journal for Numerical Methods in
Engineering, 48:1549–1570, 2000.
[28] C. Daux, N. Moes, J Dolbow, N. Sukumar, and T. Belystschko. Arbitrary branched
and intersecting cracks with the extended finite element method. International Journal
for Numerical Methods in Engineering, 48:1741–1760, 2000.
[29] Zdenek P Bazant. Instability, ductility, and size effect in strain-softening concrete.
Journal of the engineering mechanics division, 102(2):331–344, 1976.
[30] ST Pietruszczak and Z Mroz. Finite element analysis of deformation of strain-softening
materials. International Journal for Numerical Methods in Engineering, 17(3):327–334,
1981.
[31] Z.P. Bazant and B.H. Oh. Crack band theory for fracture of concrete. Materials ans
Structures, 16:155–177, 1983.
[32] A. Huerta, P. Diez, A. Rodriguez-Ferran, and G. Pijaudier-Cabot. Error estimation and
finite element analysis of softening solids. In Ladev. and J.T. Oden, editors, Advances
in Adaptive Computational Methods in Mechanics, pages 333–348. Elsevier.
[33] Ted L Anderson and TL Anderson. Fracture mechanics: fundamentals and applications.
CRC press, 2005.
[34] Egor Paul Popov. Engineering mechanics of solids. Prentice Hall, 1990.
[35] E. Onate. Structural Analysis with the Finite Element Method. Linear Statics, volume
Volume 1 : The Basis and Solids. Springer, 2010.
[36] A. Ingraffea and M. Grigoriu. Probabilistic fracture mechanics a validation of predictive
capability, 1990.
[37] H Nguyen-Xuan, GR Liu, Stephane Bordas, S Natarajan, and T Rabczuk. An adaptive
singular es-fem for mechanics problems with singular field of arbitrary order. Computer
Methods in Applied Mechanics and Engineering, 253:252–273, 2013.
[38] Samuel Geniaut and Erwan Galenne. A simple method for crack growth in mixed mode
with x-fem. International Journal of Solids and Structures, 49(15):2094–2106, 2012.
[39] SM Hausler, K Lindhorst, and P Horst. Combination of the material force concept
and the extended finite element method for mixed mode crack growth simulations.
International Journal for Numerical Methods in Engineering, 85(12):1522–1542, 2011.
56 Bibliografıa
[40] Jean-Charles Passieux, Julien Rethore, Anthony Gravouil, and Marie-Christine Baietto.
Local/global non-intrusive crack propagation simulation using a multigrid x-fem solver.
Computational Mechanics, 52(6):1381–1393, 2013.
[41] M Prasad and C Krishnamoorthy. Computational model for discrete crack growth in
plain and reinforced concrete. Computer methods in applied mechanics and engineering,
191(25):2699–2725, 2002.