DocsTec-6988

Vista previa del material en texto

2 
INSTITUTO T E C N O L Ó G I C O Y DE ESTUDIOS SUPERIORES DE M O N T E R R E Y 
C A M P U S ESTADO DE M É X I C O 
Solución Numérica de Problemas de Mecánica de Sólidos para Cuerpos 
con Grietas y/o Inclusiones 
TESIS QUE PARA OPTAR EL GRADO DE 
DOCTOR EN CIENCIAS E INGENIERÍA DE MATERIALES 
PRESENTA 
Víctor Manuel Romero Medina 
Asesor: Dr. SERGUEI KANAOUN MIRONOV 
Comité de tesis: Dr. SERGUEI KANAOUN MIRONOV 
DR. VALERIY LEVIN 
DR. OLEKSANDR TKACHENKO 
DR. SADEGH BABAII KOCHEKSARAII 
Presidente 
Secretario 
Vocal 
Vocal 
Jurado: DR. VALERIY LEVIN 
DR. OLEKSANDR TKACHENKO 
DR. SADEGH BABAII KOCHEKSARAII 
Dr. SERGUEI KANAOUN MIRONOV 
Atizapán de Zaragoza, Edo. de Méx., 15 de Enero de 2009 
© [2009] por Víctor Manuel Romero Medina 
Todos los Derechos Reservados 
3 
4 
Agradecimientos 
Agradezco fervientemente al Dr. Serguei Kanaoun por todo el apoyo que me brindo desde que 
llegó al ITESM, CEM, en 1992, por invitarme a participar con él en sus proyectos de investigación, 
de donde surgió la posibilidad de ser mi asesor durante los últimos años dando por resultado este 
trabajo. 
También agradezco al comité de tesis, Dr. Oleksandr Tkachenko, Dr. Valeriy Levin y Dr. Sadegh 
Babaii, por su interés en el presente trabajo y de quienes recibí comentarios muy valiosos para el 
desarrollo del mismo. 
A mis amigos y compañeros de siempre en el ITESM, CEM, con quienes he compartido momen­
tos importantes: Dr. Armando Bravo, Dr. Ulises Figueroa, M. C. Armando Gómez, Técnico An­
tonio García Mecedo, Dra. Olimpia Salas y Dr. Joaquín Oseguera, de quienes siempre he recibido 
consejos y comentarios que me han servido durante mi desarrollo profesional y personal y con 
quienes he sembrado una gran amistad. 
Agradezco de forma especial al Dr. Pedro Grasa Soler, por su apoyo incondicional para que 
pudiera llevar a buen término esta investigación. 
Asimismo, agradezco a mis compañeros y amigos de la Universidad del Caribe, al Rector M. A. 
Fernando Espinoza de los Reyes Aguirre, la Dra. Estela Cerezo Acevedo, y principalmente al Ing. 
Hilario López Garachana, quien me apoyo grandemente para que realizara la última etapa del 
desarrollo de esta tesis. 
5 
"...y así, del poco dormir y del mucho leer se le secó el cerebro, de manera que vino a perder el 
juicio. Llenósele la fantasía de todo aquello que leía en los libros, así de encantamientos como de 
pendencias, batallas, desafíos, heridas, requiebros, amores, tormentas y disparates imposibles; 
y asentósele de tal modo en la imaginación que era verdad toda aquella máquina de aquellas 
soñadas invenciones que leía, que para él no había otra historia más cierta en el mundo." 
El Ingenioso Hidalgo Don Quijote de la Mancha 
Miguel de Cervantes Saavedra 
6 
Resumen 
El trabajo de investigación está enfocado en la aplicación de un método numérico novedoso, de¬ 
nominado Método de Puntos de Frontera, a la solución de las ecuaciones integrales en la frontera 
que se presentan en el estudio de mecánica de cuerpos deformables, lo que permite reducir en 
una dimensión el análisis en cuestión. Principalmente lo aplicamos a la solución de los siguientes 
problemas: 
1. El segundo problema de elasticidad estático de cuerpos homogéneos en dos dimensiones, 
2. El segundo problema de elasticidad estático de cuerpos homogéneos con grieta en dos di¬ 
mensiones, 
3. El segundo problema de plasticidad estático de cuerpos homogéneos en dos dimensiones, y 
4. El segundo problema de elasticidad dinámico de cuerpos homogéneos con grieta en dos 
dimensiones. 
En el método se utilizan funciones de aproximación tipo Gauss, lo que representa una gran ven¬ 
taja de este modelo, ya que no se requiere discretizar la frontera en pequeños elementos como en 
el BEM, sino que nos permite definir una cantidad f i n i t a de nodos o puntos sobre la frontera de 
los que sólo se requiere conocer características geométricas tales como sus coordenadas y la ori¬ 
entación de vectores unitarios normales a la frontera en dichos puntos. Esto nos lleva a definir 
los componentes de la matriz de coeic ientes mediante el cálculo analítico de integrales estandar 
cuyos valores se pueden guardar fácilmente en la memoria de la computadora para ser utilizados 
posteriormente en el cálculo de esfuerzos en el interior del dominio de solución, alcanzando una 
gran precisión. Los resultados obtenidos para los diferentes problemas son comparados con solu¬ 
ciones exactas existentes en los casos (1), (2) y (4), con resultados experimentales, como en el 
caso del factor de intensidad de esfuerzos del problema (2), y con modelos numéricos de elemen¬ 
tos f i n i t o s en el caso del problema (3), mostrando gran precisión en la solución. 
7 
Abstract 
This research work is focused in the application of a relatively new numerical method defined 
as Boundary Point Method to the solution of the boundary integral ecuations that are present in 
the continuum mechanics, reducing the dimension of the problem by one. Mainly we apply the 
method to the solution of the following problems: 
1. The second boundary value problem of static elasticity for 2D homogeneous bodies, 
2. The second boundary value problem of static elasticity for 2D homogeneous bodies with 
crack, 
3. The second boundary value problem of static el astolplasticity for 2D homogeneous bod-
ies,and 
4. The second boundary value problem of dynamic elasticity for 2D homogeneous bodieswith 
crack. 
The method requires the use of Gaussian approximation functions which represents an advantage 
because it is no neccesary to discretize the boundary with small elements as in BEM, so it is only 
necessary to define certain f i n i t e quantity of nodes or points on the boundary of the body from 
which only is necessary to know geometrical properties as coordinates and direction of normal 
unit vectors at the boundary at such points. This allows us to define the elements of the coeficient 
of main matrix of the linear sistem of equations by analitical calculus of standar integrals, the 
resulting data can be saved in the computer memory to be used afterwards in the calculus of 
stresses inside the domain of solution, with a very high pressition. For problems (1), (2) and (4), 
results obtained are compared with exact solutions found in the bibligraphy, for intensity stress 
factor of problem (2) results are compared with experimental data, and for problem (3) results are 
compared with another numerical model as F E M using ANSYS. All the results show a very high 
presicion in the solution. 
8 
Lista de figuras 11 
1 Introducción 14 
1.1 Métodos numéricos aplicados en la mecánica del medio contínuo 14 
1.1.1 Método de diferencias finitas (FDM) 15 
1.1.2 Método de elementos finitos (FEM) 15 
1.1.3 Método de elementos de frontera (BEM) 16 
1.2 Métodos sin malla 17 
1.2.1 Necesidad de los métodos sin malla 19 
1.2.2 La idea de los métodos sin malla 20 
1.2.3 Método de Puntos de Frontera (BPM) 21 
1.3 Meta y objetivos 22 
2 Funciones de aproximación 24 
3 Problema de elasticidad bidimensional 28 
3.1 Formulación del problema plano en lenguaje de ecuaciones integrales 28 
3.2 Acción de los operadores del problema sobre las funciones de Gauss 32 
3.3 Discretización de las ecuaciones integrales del segundo problema de elasticidad. 
Caso plano 35 
3.4 Ejemplos de solución numérica del problema 42 
3.5 Solución del problema para cuerpos con grieta 46 
3.6 Conclusiones 56 
4 Problema elasto-plástico para cuerpos en 2D 57 
Contenido 
9 
4.1.2 Teoría de la plasticidad 59 
4.1.3 Especificaciones para el endurecimiento 63 
4.1.3.1 Endurecimiento isotrópico multilineal y endurecimiento 
isotropico bilineal 63 
4.1.3.2 Especificaciones para endurecimiento isotrópico no lineal 63 
4.1.3.3 Especificaciones para el endurecimiento cinemático bilineal 64 
4.2 Formulación de las ecuaciones elasto-plásticas en forma de ecuaciones integrales . . . . 654.3 Acción de operadores integrales del problema elasto-plástico sobre las funciones 
de Gauss 68 
4.4 Ejemplo de aplicación. Deformación elastoplástica de una placa rectangular con 
un corte lateral 70 
4.4.1 La ley de plasticidad 70 
4.4.2 Algoritmo de la solución numérica 72 
4.4.3 Resultados numéricos 74 
4.4.4 Comparación de resultados del modelo utilizando ANSYS 76 
4.5 Conclusiones 78 
5 Problema de elasticidad dinámica para cuerpos con grieta 79 
5.1 Formulación de las ecuaciones del problema de elasticidad dinámica en forma de 
ecuaciones integrales 79 
5.2 Solución numérica del problema 2D 80 
5.3 Aplicación al problema de difracción de ondas P planas sobre una grieta 84 
5.3.1 Conclusiones 86 
6 Conclusiones 88 
Bibliografía 90 
4.1 Diferentes formas de la teoría de plasticidad 57 
4.1.1 Plasticidad independiente del t iempo o de la razón (rate independent plasticity) 58 
10 
A Tensores de cuarto orden especiales 93 
A.1 Las bases tensoriales de cuarto orden 93 
A.1.1 Base E 93 
A.1.2 Base P 95 
A.1.3 Base θ 97 
B Algunas integrales sobre la esfera unitaria 98 
C Métodos Numéricos 100 
C.1 Métodos de solución de sistemas de ecuaciones lineales 100 
C.1.1 Conjuntos de ecuaciones no singulares versus conjuntos singulares 100 
C.1.2 Representación matricial 102 
C.1.3 Eliminación de Gauss-Jordan 102 
C.1.4 Eliminación Gaussiana con substitución hacia atrás 103 
C.1.5 Descomposición L U y sus aplicaciones 104 
C.1.6 Efectuando la descomposición L U 106 
C.2 Interpolación en dos dimensiones 108 
C.2.1 Orden superior para mejorar la precisión 109 
C.2.2 Orden superior para suavizar la función mediante interpolación bicúbica 110 
D Programas computacionales en lenguaje C 112 
11 
Lista de figuras 
Fig. 2.1 Gaussianas 24 
Fig. 2.2 Gráfica de f (x) 24 
Fig. 2.3 Gráfica de f (x) aumentada 25 
Fig. 2.4 f1/2(x) y términos individuales 26 
Fig. 2.5 f 4 ( x ) y términos individuales 26 
Fig. 3.1 R (x) de la aproximación Gaussiana de la ec.(3.21) 33 
Fig. 3.2 R (x) de la aprox. Gaussiana de la ec.(3.21). D = 2 34 
Fig. 3.3 R (x) de la aprox. Gaussiana de la ec.(3.22). D = 2 34 
Fig. 3.4 Bases global ( e 1 , e 2 ) y local ( s
( i ) , n ( i ) ) de un área plana Ω con frontera T 36 
Fig. 3.5 Orificio de una placa plana infinita 43 
Fig. 3.6 Solución exacta. Componente b s 44 
Fig. 3.7 Solución exacta. Componente bn 45 
Fig. 3.8 Esfuerzos en la dirección y, σ11 45 
Fig. 3.9 Disco circular con R = 1, y cargas concentradas opuestas 46 
Fig. 3.10 Distribución de σ11 (x) en el disco sujeto a dos fuerzas de compresión 
concentradas 47 
Fig. 3.11 Distribución de σ11(y) en el disco sujeto a dos fuerzas de compresión 
concentradas 47 
Fig. 3.12 Error R ( x ) de la sol. numérica de una grieta recta en un plano infinito 48 
Fig. 3.13 Sistema de coordenadas local en la punta de la grieta 48 
Fig. 3.14 Área rectangular con grieta central 51 
Fig. 3.15 Distribución del esfuerzo σ11(x, 0) a lo largo de la grieta de un área rectangular. 54 
12 
Fig. 3.16 Distribución del esfuerzo σ22(x, 0) a lo largo de la grieta de un área rectangular. 54 
Fig. 3.17 Factores de intensidad de esfuerzos para una grieta central en un área plana 
rectangular 55 
Fig. 4.1 Diferentes superficies de cedencia 60 
Fig. 4.2 Curva esfuerzo deformación para el análisis de edurecimiento no lineal en 
ANSYS 64 
Fig. 4.3 Representación en 2D de un cuerpo con región de deformación plástica 65 
Fig. 4.4 Región de deformación plástica en una placa rectangular en 2D 71 
Fig. 4.5 Lineas A y B en la región de deformación para la comparación de σ and ε 75 
Fig. 4.6 Comparación de σ sobre la linea A para mallas con 2 1 , 4 1 , 61 y 81 nodos 
por lado 75 
Fig. 4.7 Comparación de ε sobre la linea A para mallas con 2 1 , 4 1 , 61 y 81 nodos 
por lado 76 
Fig. 4.8 Comparación de σ sobre la linea B para mallas con 2 1 , 41 , 61 y 81 nodos 
por lado 76 
Fig. 4.9 Comparación de ε sobre la linea B para mallas con 21 , 4 1 , 61 y 81 nodos 
por lado 77 
Fig. 4.10 Esfuerzo equivalente de plasticidad en el eje 77 
Fig. 4.11 Deformación equivalente de plasticidad en el eje 78 
Fig. 5.1 Distribución de esfuerzos en un disco sujeto a dos fuerzas 84 
Fig. 5.2 Apertura de la grieta y factores de intensidad de esfuerzos para el problema 
de difracción de ondas P 86 
Fig. 6.1 Malla cuadrada para interpolación bicúbica 108 
13 
1 Tensor identidad 
A , Aij Matriz de coeficientes del sistema de ecuaciones lineales 
b Vector del potencial de doble capa 
C , C i j k i Tensor módulo elástico 
D Parámetro adimensional de la función tipo Gauss, D = 2 
E, E i j k i Base tensorial de cuarto orden 
F (‧), (‧) Transformada de Fourier 
G , Gij Función de Green 
h Distancia entre puntos de frontera 
I, I ijkl Tensor identidad de cuarto orden 
KI, KII Factores de intensidad de esfuerzos 
n,ni Vector normal exterior a la frontera 
P, Pijkl Base tensorial de cuarto orden 
S, Sijkl, T, Tijkl Funciones tensoriales de cuarto orden 
U, Ui Función potencial de doble capa 
X Vector solución del sistema de ecuaciones lineales 
a = b/a Razón de dimensiones de placa para placa con grieta 
β Angulo que forma el vector normal con el eje y en el punto frontera 
٦ Recta tangente a la frontera en el punto frontera 
δij Delta de Kronecker 
δ (x) Función delta de Dirac 
εpl, εpl Tensores deformación plástica, deformación plástica equivalente 
η Coordenadas adimensionales 
λ,μ Coe ic ien tes de Lamé 
θ, θijkl Base tensorial de cuarto orden 
Γ Frontera del dominio de solución 
V Razón de Poisson 
p Densidad del cuerpo 
σ,σ, σpl Tensores de esfuerzos, esfuerzos equivalentes, esfuerzos de plasticidad 
ω Frecuencia 
Ω Dominio de solución 
Nomenclatura 
14 
1 Introducción 
Actualmente sabemos que si no existieran las técnicas numéricas, sería casi imposible resolver 
analíticamente problemas de ingeniería prácticos con un grado razonable de precisión. La mayoría 
de las técnicas numéricas en mecánica del medio continuo están basados en el principio de que es 
posible derivar algunas ecuaciones y relaciones que describan con precisión el comportamiento 
de una pequeña parte diferencial de un cuerpo. Dividiendo el cuerpo entero en un gran número 
de estas partes pequeñas y utilizando las ecuaciones correspondientes para conectarlas, es posible 
obtener una predicción de los valores de las variables con precisión razonable tales como esfuerzos 
y desplazamientos en el cuerpo. Conforme el tamaño de estas partes se hacen más pequeñas, la 
solución numérica se hace más precisa, pero el costo del t iempo de cómputo puede ser prohibitivo. 
N o exsite un substituto para la experiencia en la aplicación de las técnicas numéricas a problemas 
de ingeniería prácticos porque no hay una respuesta clara a la pregunta ¿Qué tan pequeñas deben 
ser las partes componentes para obtener la precisión óptima? 
Por lo anterior es importante revisar y comparar las diferentes técnicas numéricas que se han 
estado desarrollando para la solución de problemas de la mecánica del medio continuo. 
A continuación daremos una brede descripción de cada uno de ellos. 
1.1 Métodos numéricos aplicados en la mecánica del medio continuo 
Actualmete es posible clasificar los métodos numéricos aplicados en la mecánica del medio con-
tínuo como métodos con malla y métodos sin malla. Para los métodos con malla podemos distin¬ 
guir los más comunes como son el método de diferencias f i n i t a s (FDM, Finite Diference Method), 
el método de elementos f i n i t o s (FEM, Finite Element Method) y el método de elementos de fron¬ 
tera (BEM, Boundary Element Method) ; para los métodos sin malla, aunque se han desarrollado 
una gran variedad de ellos, nos enfocaremos al método que es de nuestro interés particular, el 
método de puntos de frontera (BPM, Boundary Point Method) . 
15 
1.1.1 Método de diferencias finitas (FDM) 
En este método, las derivadas en las ecuaciones diferenciales parciales que gobiernan el fenó¬meno físico son discretizadas en forma de ecuaciones en diferencias. Por lo tanto, el cuerpo o 
dominio de solución es discretizado o dividido como una malla de celdas y la aproximación en 
diferencias es aplicada en cada celda interior. Esto resulta en un sistema de ecuaciones algebrái-
cas lineales (con una matriz de coeic ientes en banda) que proporciona una solución única ya que 
las condiciones de frontera del problema real son satisfechas. 
La aproximación en diferencias f i n i t a s es la más simple y relativamente la más fácil de programar. 
Su limitación principal al aplicarlo a problemas de ingeniería prácticos es que no es adecuado para 
problemas con geometrías con fuertes irregularidades. Por otro lado, debido a que es dificil variar 
el tamaño de las celdas en regiones particulares, no es adecuado para problemas donde los valores 
de las variables cambian rápidamente, tal como en problemas con concentración de esfuerzos. En 
la actualidad, los métodos de diferencias f i n i t a s son populares para problemas de transferencia de 
calor y mecánica de fluidos, más que para problemas de análisis de esfuerzos. 
1.1.2 Método de elementos finitos (FEM) 
En la actualidad, los métodos de elementos finitos son ampliamente usados en análisis en inge-
niería. Estos métodos son usados extensamente en el análisis de sólidos y estructuras, en trans-
ferencia de calor y dinámica de fluidos; de hecho, los métodos de elementos finitos son útiles 
prácticamente en cada campo de análisis de ingeniería. 
En este método, el dominio de solución entero es dividido en pequeños elementos (elementos 
finitos). Sobre cada elemento, el comportamiento es descrito por las ecuaciones diferenciales que 
gobiernan el fenómeno físico. Todos estos pequeños elementos se mantienen conectados y los 
requerimientos de continuidad y equilibrio son satisfechos entre los elementos vecinos. Dado que 
las condiciones de frontera del problema real son satisfechas, se puede obtener una solución única 
con el sistema de ecuaciones lineales obtenido (en este caso, la matriz de coeic ientes está poblada 
de forma dispersa). 
16 
El F E M es muy adecuado para problemas prácticos de ingeniería de geometrías complejas. Para 
obtener buena precisión en regiones donde los valores de las variables cambian rápidamente se 
debe utilizar un gran número de elementos f i n o s . 
1.1.3 Método de elementos de frontera (BEM) 
Los ingenieros que han trabajado con programas desarrollados bajo el F E M se pueden preguntar 
por qué es necesario crear otras técnicas de cómputo. La respuesta es que el F E M ha demostrado 
ser inadecuado o ineficaz en muchas aplicaciones de la ingeniería. El análisis con elementos f i n i t o s 
sigue siendo un proceso comparativamente lento debido a la necesidad de definir o de redefinir 
mallas en la pieza o el dominio bajo estudio. 
El método de elementos de frontera (BEM) [1] surgió como una alternativa potencial a elementos 
f in i to s , particularmente para casos donde se requiere una mayor exactitud debido a los problemas 
tales como concentración de esfuerzos o donde el dominio se extiende al infinito. Sin embargo, 
la característica más importantes del método de elementos de frontera es que solamente requiere 
la discretización de la frontera más que del dominio de solución. Por lo tanto, los códigos de el¬ 
ementos de frontera son más fáciles de utilizar con los modeladores de sólidos y los generadores 
de mallas existentes. Esta ventaja es particularmente importante para el diseño, ya que el proceso 
implica generalmente una serie de modificaciones que son más difíciles de realizar usando ele¬ 
mentos f in i tos . Las mallas pueden ser generadas fácilmente y los cambios de diseño no requieren 
mallado completo. 
En este método, las ecuaciones diferenciales que gobiernan el fenómeno físico son transformadas 
en ecuaciones integrales que son aplicables sobre la frontera del dominio de solución. Estas ecua¬ 
ciones son integradas numéricamente sobre la frontera que se divide en pequeños segmentos de 
frontera (o elementos de frontera). Como en los métodos anteriores, dado que las condiciones de 
frontera son satisfechas, se obtiene un sistema de ecuaciones lineales para el que se obtiene una 
solución única. 
Con este método se pueden acomodar fronteras geométricamente complejas. Además, ya que 
todas las aproximaciones están restringidas a la frontera, se pueden modelar regiones con variables 
que cambian su valor rápidamente con mejor precisión que el FEM. 
17 
1.2 Métodos sin malla 
Los métodos sin malla, se utilizan para establecer un sistema de ecuaciones algebraicas para el 
dominio entero del problema sin el uso de una malla predefinida. Los métodos sin malla utilizan 
un sistema de nodos dispersos dentro del dominio del problema, así como sistemas de nodos 
dispersos en las fronteras del dominio para representar (no discretizar) el dominio del problema 
y sus fronteras. Estos sistemas de nodos dispersos no forman una malla, lo que significa que no se 
requiere ninguna información sobre la relación entre los nodos, por lo menos para la interpolación 
de las variables del campo. 
Hay un número de métodos sin malla, tales como el método de Galerkin libre de elementos (EFG) 
[2], el método local sin malla de Petrov-Galerkin (MLPG) [3], el método de la interpolación 
puntual (PIM) [4], el método del conjunto de puntos (PAM) [5], el método del punto finito [6], 
el método de diferencias finitas con mallas irregulares arbitrarias [7, 8], y otros más. Todos estos 
métodos comparten la misma característica que no requieren mallas predefinidas, por lo menos 
para la interpolación de las variables del campo. Los nombres para los varios métodos sin malla 
todavía se están discutiendo. Debido a que la metodología todavía está en etapa de desarrollo, 
constantemente se proponen nuevos nombres de métodos. Puede tardar cierto t iempo antes que 
todos los métodos se categorizen y se unifiquen adecuadamente para evitar confusiones en la 
comunidad. 
En contraste con el FEM, se prefiere el término método libre de elementos, y en contraste con 
FDM, se prefiere el término método de diferencias finitas usando mallas arbitrarias o irregulares. 
El requisito mínimo para considerar un método como método sin malla es: 
• Que una malla predefinida no es necesaria, por lo menos en la interpolación de las variables 
del campo. 
El requisito ideal para un método sin malla es: 
• Que no se requiere malla alguna durante el proceso de solución del problema de una geometría 
arbitraria dada gobernada por un sistema de ecuaciones diferenciales parciales sujetas a cualquier 
tipo de condiciones de frontera. 
18 
La realidad es que los métodos sin malla desarrollados hasta ahora no son realmente ideales, y 
fallan en una de las categorías siguientes: 
• Métodos que requieren celdas de respaldo para la integración de las matrices del sistema 
derivadas de la forma débil sobre el dominio del problema. Estos métodos no están realmente 
libres de malla. Los métodos E F G pertenecen a esta categoría. Estos métodos son prácticos en 
gran medida, siempre y cuando sea generalmente factible la creación de una malla de fondo, 
y siempre puede ser automatizada usando una malla triangular para dominios bidimensionales 
(2D) y mallas tetrahédricas para dominios tridimensionales (3D). 
• Métodos que requieren celdas de respaldo locales para la integración de las matrices del sistema 
sobre el dominio del problema. Los métodos M L P G pertenecen a esta categoría. Se dice que 
estos métodos son esencialmente libres de malla porque la creación de una malla local es más 
fácil que crear una malla para el dominio entero del problema; es una tarea más simple que se 
puede realizar automáticamente sin ninguna predefinición para la malla local. 
• Métodos que no requieren una malla, pero que son menos estables y menos precisos. Los méto¬ 
dos de colocacióny los métodos de diferencias f i n i t a s que utilizan mallas irregulares pertenecen 
a esta categoría. La selección de nodos basados en el tipo de un problema físico es todavía im¬ 
portante para obtener resultados estables y precisos [9, 10, 11]. La automatización de la selec¬ 
ción de los nodos y el mejoramiento en la estabilidad de la solución son algunos de los desafíos 
en estas clases de métodos. Este tipo de métodos tienen una ventaja muy significativa: son muy 
fáciles de implementar, porque no se requiere ninguna integración. 
• Métodos de partículas que requieren una predefinición de estas para sus volúmenes o masas. 
El algoritmo entonces realizará los análisis incluso si el dominio del problema experimenta 
la deformación y separación extremadamente grandes. Los métodos SPH pertenecen a esta 
categoría. Este tipo de métodos sufren de problemas en la imposición de las condiciones de 
frontera. Además, la predefinición de las partículas todavía requiere técnicamente una cierta 
clase de malla. 
19 
1.2.1 Necesidad de los métodos sin malla 
El F E M es robusto y se ha desarrollado a fondo para el análisis de esfuerzos estáticos y dinámi¬ 
cos, lineales o no lineales de sólidos y estructuras, así como de flujos de fluidos. La mayoría de los 
problemas prácticos de la ingeniería relacionados con sólidos y estructuras se solucionan actual¬ 
mente usando una gran cantidad de paquetes bien desarrollados del F E M que están disponibles en 
el mercado. Sin embargo, las siguientes limitaciones del F E M se están poniendo cada vez más de 
manifiesto: 
1. La creación de una malla para el dominio del problema es un requisito previo al usar los 
paquetes del FEM. El analista pasa generalmente la mayor parte de su t iempo creando la 
malla, y se convierte en un componente importante del costo de un proyecto de simulación 
porque el costo del t iempo de la CPU (unidad central de procesamiento) está disminuyendo 
drásticamente. La preocupación es más por el t iempo de la mano de obra, y menos por 
el t iempo de computadora. Por lo tanto, idealmente el proceso de mallado sería realizado 
completamente por la computadora sin la intervención humana. 
2. Durante los cálculos, los esfuerzos obtenidos son discontinuos y menos exactos usando los 
paquetes de FEM. 
3. Al manejar deformaciones grandes, la precisión se pierde considerable debido a las distor¬ 
siones de los elementos. 
4. Es muy difícil simular crecimiento de grieta con las trayectorias arbitrarias y complejas y las 
transformaciones de la fase debido a las discontinuidades que no coinciden con las líneas 
nodales originales. 
5. Es muy difícil simular la fractura del material en una gran cantidad de fragmentos pues el 
F E M esencialmente se basa en la mecánica del medio continuo, en la que los elementos for¬ 
mulados no pueden estar divididos. Los elementos pueden ser, o totalmente fragmentados o 
permanecer como una pieza completa. Esto lleva generalmente a una mala representación de 
la trayectoria de la fractura. Un error serio puede ocurrir porque la naturaleza del problema 
es no lineal; y por lo tanto, los resultados son altamente dependientes de la trayectoria. 
20 
6. Se han propuesto aproximaciones de re-mallado para manejar estos tipos de problemas en el 
FEM. En la aproximación del re-mallado, el dominio del problema es re-mallado en cada 
paso durante el proceso de la simulación para prevenir la distorsión severa de elementos y 
para permitir que las líneas nodales sigan siendo coincidentes con las fronteras de la discon¬ 
tinuidad. Con este i n , se han desarrollado procesadores complejos, robustos, y generadores 
de mallas adaptivas. Sin embargo, estos procesadores son solamente realizables para proble¬ 
mas en 2D. N o hay procesadores confiables disponibles para crear mallas hexahédricas para 
problemas 3D debido a la dificultad técnica. 
7. Los procesadores adaptantes requieren el mapeo de las variables del campo entre las mallas 
en etapas sucesivas durante la solución del problema. Este proceso de mapeo lleva a menudo 
a realizar cálculos adicionales, así como una degradación de la exactitud. Además, para los 
problemas grandes en 3D, el costo de cómputo de re-mallar en cada paso llega a ser muy 
alto, incluso si un esquema adaptivo está disponible. 
8. Los F D M trabajan muy bien para una gran cantidad de problemas, especialmente para solu¬ 
cionar problemas de la dinámica de l u idos . Pero adolecen de una desventaja importante, en 
que dependen de nodos regularmente distribuidos. Por lo tanto, se han efectuado estudios du¬ 
rante mucho t iempo para desarrollar métodos que utilicen mallas irregulares. Actualmente, 
se siguen realizando esfuerzos en esta dirección. 
1.2.2 La idea de los métodos sin malla 
Una evaluación enfocada a estas dificultades asociadas al F E M revela la raíz del problema: la 
necesidad de utilizar los elementos, que son la base de construcción del FEM. Se requiere una 
malla con conectividad predefinida para formar los elementos. Mientras se requiera utilizar ele¬ 
mentos, los problemas mencionados anteriormente no serán fáciles de solucionar. Por lo tanto, la 
idea de eliminar los elementos, y por lo tanto las mallas, ha evolucionado naturalmente. se han 
propuesto los conceptos de métodos libres de elementos o libres de malla, en los que el dominio 
del problema está representado por un sistema de nodos arbitrariamente distribuidos. 
21 
Los métodos sin malla tienen gran potencial para solucionar los problemas difíciles menciona¬ 
dos arriba. Los esquemas adaptivos pueden ser desarrollados fácilmente, pues no hay malla, y por 
lo tanto ningun concepto de conectividad esta implicado. Así, no hay necesidad de proporcionar 
a priori ninguna información sobre la relación entre los nodos. Esto proporciona f lexibi l idad en 
la adición o la supresión de nodos o puntos siempre y donde sean necesarios. Por ejemplo, para 
el análisis de esfuerzos de un dominio sólido, existen a menudo áreas de concentración de es¬ 
fuerzos, incluso singularidades. Uno puede agregar, con relativa libertad, puntos en el área de 
concentración de esfuerzos sin la preocupación de su relación con los otros puntos existentes. En 
problemas de crecimiento de grieta, los puntos se pueden agregar fácilmente alrededor de la nariz 
o punta de la grieta para capturar la concentración de esfuerzos con la exactitud deseada. Este rei¬ 
namiento nodal se puede desplazar junto con la propagación de la grieta con un arreglo de puntos 
de apoyo asociados con la geometría global. Las mallas adaptivas de una gran variedad de prob¬ 
lemas, en 2D o 3D, incluyendo análisis de esfuerzos lineales y no lineales, estáticos y dinámicos, 
se pueden tratar con mucha eficacia en los métodos sin malla de una forma relativamente simple. 
Debido a que no hay necesidad de crear una malla, y los puntos se pueden crear por una com¬ 
putadora en una manera completamente automatizada, el t iempo que un ingeniero pasaría en la 
generación convencional de una malla puede ser disminuído. Esto se puede traducir en ahorros 
substanciales del costo y del t iempo en proyectos del modelado y de la simulación. 
1.2.3 Método de Puntos de Frontera (BPM) 
El método de puntos de frontera (BPM) es un método sin malla relativamente nuevo desarrolla¬ 
do para obtener la solución de una amplia gama de ecuaciones integrales de la físico matemática 
propuesto por Vladimir Maz ' ya [13, 14] y V. Maz 'ya y G. Shmidt [15]. El análisis de multires-
olución basado en estas funciones fue propuesto en los trabajos de V. Maz 'ya y G. Shmidt [16]. 
La utilización de estas funciones para la solución de las ecuaciones integrales de elasticidad tiene 
dos ventajas principales. La primera, es que la acción de los operadores integrales del problema 
sobre estas funciones son combinaciones de unas cuantas funciones estandar. Estas, simplemente 
pueden ser tabuladas, mantenidas en la memoria dela computadora y posteriormente ser usadas 
para la solución de cualquier problema de elasticidad. Como resultado, el t iempo para el cálculo 
22 
1.3 Meta y objetivos 
Por lo descrito anteriormente podemos establecer que esta tesis tiene como meta fundamental 
presentar una metodología en la que se aplica el Método de Puntos de Frontera (BPM, Boundary 
Points Method) a la solución de problemas de elasticidad, tanto estáticos como dinámicos, de 
cuerpos homogéneos con grietas y/o inclusiones y de plasticidad. Para alcanzar esta meta se tienen 
establecidos los siguientes objetivos generales: 
• Aplicar la metodología a la solución de la ecuación integral del segundo problema de elas¬ 
ticidad de cuerpos homogéneos sin y con grieta, comprobando su aplicabilidad con algunos 
problemas cuya solución exacta es conocida. 
• Aplicar la metodología a la solución de la ecuación integral del segundo problema de elasto-
plasticidad de cuerpos homogéneos, comprobando su aplicabilidad con otros métodos, como 
por ejemplo, FEM, mediante el uso del programa comercial ANSYS. 
de la matriz del sistema de ecuaciones lineales obtenido después de la discretización del problema 
es reducido sustancialmente comparado con el BEM. 
También es importante hacer notar que el B P M no requiere la construcción de una malla en la 
frontera. En lugar de eso, utiliza las coordenadas de un número finito de puntos en la frontera y la 
orientación de la frontera a partir de vectores unitarios normales exteriores a la frontera en cada 
punto. Esto representa menor información que la requerida en el BEM, donde la forma de los 
elementos de las fronteras deben ser definidas. Entonces los puntos de frontera en el B P M juegan 
el papel de los elementos de frontera en el BEM. La idea general es utilizar una función base 
localizada (tal como una distribución Gaussiana) centrada alrededor de cada punto de la frontera 
concentradas en los planos tangentes a la frontera en cada punto o nodo. La teoría de aproximación 
en la aplicación más simple del método, las integrales de frontera resultantes son aproximadas por 
integrales sobre la tangente a cada punto de la frontera, llevando una aproximación de segundo 
orden a la solución exacta. La simetría radial de las funciones base ayudan en la determinación de 
expresiones analíticas para las integrales de frontera requeridas. 
23 
• Aplicar la metodología a la solución de la ecuación integral del problema de elasticidad dinámi¬ 
ca de cuerpos homogéneos con grieta, comprobando su aplicabilidad con algunos problemas 
cuya solución exacta es conocida. 
2 Funciones de aproximación 
24 
Fig. 2 .1 : Gaussianas. 
Por supuesto, la función f (x) está acotada, es positiva y suave; además, f (x + 1) = f (x) , es 
decir, es periódica con período 1, por lo que esperamos que la gráfica de f debería verse como 
una curva ondulada periódica agradable. Sin embargo, es interesante encontrar que esta gráfica es 
una constante, como se muestra en la f i g u r a (2.2). 
Supongamos que se nos da la tarea de dibujar la gráfica de la función 
obtenida mediante la suma de las Gaussianas desplazadas, como se muestra en la figura (2.1) 
Fig. 2.2: Gráfica de f (x) . 
25 
Fig. 2.3: Gráfica de f (x) aumentada. 
Con diferentes valores del parámetro D > 0. Las figuras (2.4, 2.5) muestran las gráficas de fD 
para los parámetros D = 0 , 5 y D = 4 , respectivamente. La gráfica de la función f1/2(x) muestra 
el comportamiento oscilante, mientras que f4 parece ser una constante. En efecto, f4 también 
es oscilante alrededor de 3 , 5 4 4 9 0 7 7 0 1 8 1 1 0 3 2 0 5 ± 1 , 4 3 x 1 0 - 1 5 , que es muy difícil de graficar. 
Podemos concluir, que la función oscilante fD tiende a una constante si D se incrementa. 
Para explicar rigurosamente las peculiaridades de las gráficas, consideremos la serie de Fourier de 
la función 
En efecto, esta impresión superficial está equivocada. Si la escala del eje y se modifica como 
en la figura (2.3), entonces podemos ver que f (x) no es constante y oscila entre 2 , 5 0 6 6 2 8 2 6 y 
2 , 5 0 6 6 2 8 2 9 . 
Se obtiene la misma gráfica si el procedimiento se repite para la suma 
(2.1) 
Sus coeficientes se pueden calcular como sigue 
26 
Fig. 2.5: f4(x) y términos individuales. 
El orden de la sumatoria y la integración se puede intercambiar debido a la convergencia absoluta 
de la sumatoria infinita. Entonces obtenemos la serie de Fourier 
(2.2) 
Esta representación de la función 6 es un caso especial de la conocida sumatoria de Poisson 
Fig. 2.4: f1/2(x) y términos individuales. 
(2.3) 
27 
donde F u denota la transformada de Fourier de la función u. De la ecuación (2.2), tenemos 
Los coeficientes e - π 2 D v 2 , v = 1, 2 , . . . , pueden ser muy pequeños dependiendo de D, como se 
puede ver de la relación e - π 2 = 5.172 3 x 1 0 - 5 . En particular, si D ≥ 1, entonces para cualquier 
x, el módulo de la ecuación (2.4) es menor que 1,04 x 1 0 - 4 D . Nótese que para los casos D = 2 
y D = 4 la diferencia es comparable a las presiciones simple y doble, respectivamente, en la 
aritmética de las computadoras modernas; es decir, en estos casos donde la función θ (x , D) es 
numéricamente la función constante 1. Además, la diferencia |θ(x , D) — 1| se puede hacer menor 
que cualquier tolerancia positivo prescrita ε, escogiendo D suficientemente grande. Para ello es 
su ic iente considerar 
Nota. La función θ está fuertemente conectada con la función Teta de Jacobi V3, que está definida 
como [17] 
es decir, nuestra función θ (x , D) difiere de 1 por la serie infinita 
(2.4) 
D > π2 (log |ε| - l o g 2 ) . 
mediante la relación θ(x, D) = V3 ( Π X | i Π D ) . 
28 
Donde δ (x) es la función delta de Dirac y δij es el símbolo delta de Kronecker. 
3 Problema de elasticidad bidimensional 
3.1 Formulación del problema plano en lenguaje de ecuaciones inte¬ 
grales 
Consideremos un cuerpo elástico que ocupa una región en el espacio de dos o tres dimensiones 
con una frontera suave Γ, el material del cuerpo es homogéneo con el tensor módulo elástico 
C ( C i j k l son los componentes de este tensor). El vector desplazamiento u (x) de los puntos del 
cuerpo satisface la ecuación 
Aquí, q es el vector de fuerzas de cuerpo, xi (i = 1, 2, 3 para el caso tridimensional o i = 1, 2 
para el caso bidimensional) son las coordenadas cartesianas de cada punto x del cuerpo. En la 
notación consideramos que se realiza la suma con respecto a índices repetidos. El segundo prob-
lema de elasticidad de valores en la frontera es la solución de la ecuación (3.1) con las siguientes 
condiciones en la frontera Γ: 
Aquí, σ (x) es el tensor de esfuerzos con componentes σij, n ( x ) es un vector normal dirigido 
hacia el exterior de la frontera Γ en el punto x EΓ, f (x) es el vector de fuerzas aplicadas en la 
frontera del cuerpo, el punto significa operación producto escalar entre tensores y vectores; es 
decir, σ • n → σ i j n j . 
La función de Green, G (x) , es el operador diferencial en el lado izquierdo de la ecuación (3.1) 
que satisface la ecuación siguiente 
(3.1) 
(3.2) 
(3.3) 
σ (x) • n ( x ) ¦ Γ = f (x) , σ i j ( x ) = C i j k l 
u l (x). 
29 
La integral de frontera 
se llama el potential de doble capa con vector densidad b (x ) . Este potential satisface la ecuación 
(3.1) en todo el dominio excepto sobre la frontera Γ, en donde tiene una discontinuidad. El salto 
del potential U (x) sobre Γ es igual a su densidad b (x) [19, 20, 24]. 
En esta última ecuación, U ( x 0 + 0) es un valor límite del potential U (x) cuando x → x 0 E Γ 
desde el lado del vector normal externo n ( x 0 ) en el punto x 0 , y U ( x 0 — 0) es el mismo límite 
cuando x → x 0 E Γ desde el lado opuesto del vector normal n ( x 0 ) . 
El campo de esfuerzos Σ (x) que corresponde al potential U (x) tiene la forma siguiente: 
Debido a que la función U (x) tiene un salto sobre la frontera Γ, esta integral contiene una funcióndelta singular concentrada sobre Γ Introduzcamos la función tensorial de rango cuatro, S (x) , 
como sigue: 
Entonces, el potencial 
coincide con el potencial Σ (x) en todas partes excepto en la frontera Γ. Es posible mostrar [20, 
24], que Σ1 (x) no tiene una componente singular concentrada sobre Γ, y Σ1 (x) satisface la 
ecuación 
en todo el espacio. El salto del tensor Σ1 (x) sobre la frontera Γ tiene la forma [20, 24] 
(3.4) 
(3.5) 
(3.6) 
(3.7) 
(3.8) 
(3.9) 
[U (x0)] = U (x0 + 0) - U (x0 - 0) = b ( x 0 ) , x0 E Γ. 
Σ1ij(x) = Sijkl (x - x ' ) nk (x ') bl (x') dΓ. 
[Σ1ij (x0)] = - S i j k L (n0) n0kbL (x0) , x0 E Γ , 
30 
Aquí n 0 = n ( x 0 ) , S ( n 0 ) = S (k) ¦ k = n o y S (k) es la transformada de Fourier de la función S (x) 
definida en la ecuación (3.7). De la ecuación (3.7) se obtiene la siguiente expresión para S (k) : 
S i j k l ( k ) = C i j k l k k k m G l s ( k ) C m s p q - C i j p q , Gls ( k ) = [ k i k j C i l j s ]
- 1 , ( 3 . 1 0 ) 
donde G (k) es la transfromada de Fourier de la función de Green definida en la ecuación (3.3) 
Para cuerpos isotrópicos Gls (k) tiene la forma 
donde para el caso tridimensional, siendo λ0 y μ0 los coeficientes de Lamé del 
medio. Para el problema plano, N0 toma una de las formas siguientes, según sea el caso: 
para estado de deformación plana. 
para estado de esfuerzos planos. 
y v0 es la razón de Poisson. De acuerdo con las ecuaciones (3.9, 3.10) la componente normal de 
Σ1 (x) es continua sobre la frontera Γ; es decir, 
[n0 • Σ
1 ( x 0 ) ] = 0. (3.12) 
Las propiedades mencionadas de Σ1 (x) nos permite utilizar este potencial para la determinación 
del tensor de esfuerzos σ (x) en la región Ω con frontera Γ. Entonces, encontremos σ (x) en la 
forma 
(3.11) 
σij (x) = Σ
1
ij (x) = S i j k l (x - x ' ) nk (x ' ) bl (x ') dΓ'. 
La ecuación integral para la densidad desconocida b (x) en esta ecuación se obtiene a partir de la 
condición de frontera (3.2), y toma la forma 
T (x, x ' ) • b (x') dΓ = f ( x ) , x E Γ, 
Tij (x, x ' ) = nk (x) S k i j l (x - x ' ) nl ( x ) . 
(3.14) 
(3.15) 
El núcleo del operador integral en la ecuación (3.14) tiene la sigularidad fuerte 
T (x, x ' ) ̃ ¦x ̶ x ' ¦ - 3 cuando x ' → x 
31 
Aquí, Γ se considera como una frontera cerrada suave; la integral en el lado derecho, se entiende 
en el sentido de su valor principal de Cauchy (v.p.), b (x) es una función finita suave. La misma 
fórmula de regularización se puede aplicar a una frontera abierta si b (x) → 0 cuando x → Γ, 
donde Γ es el contorno frontera de Γ. 
En el caso de cuerpos con grietas, las fronteras de las grietas ΓC deben ser incluidas como parte de 
la frontera total del cuerpo. Las condiciones de frontera sobre Γc tienen la forma siguiente: 
n (x) ‧ σ (x) |Γ c = 0 (3.17) 
si los lados de las grietas están libres de esfuerzos. Aquí, n (x) es un vector normal arbitrario a Γc 
en el punto x E ΓC. 
Nota 1. Las ecuaciones (3.13, 3.14) también se pueden aplicar a la solución del problema de 
elasticidad para cuerpos con grietas cuyos lados no estén libres de esfuerzos; pero las fuerzas 
aplicadas a los diferentes lados de las fronteras de la grieta deben ser de la misma magnitud y 
de direcciones opuestas. en este caso, la condición de frontera (3.17) se cambia por la condición 
n ‧ σ|Γc = f, donde f es la fuerza aplicada en el lado positivo de la grieta (Cualquier lado de la 
grieta se puede escoger como positivo y n es un vector normal externo a este lado.) 
Nota 2. El potencial U (x) en la ecuación (3.4) representa los desplazamientos de los puntos del 
cuerpo que corresponde al tensor de esfuerzos en la ecuación (3.13). A partir de la propiedad 
(3.5) de U (x) , se sigue que el vector b (x) en las ecuaciones (3.13, 3.14) es el salto del vector de 
desplazamiento sobre la frontera de la grieta (apertura de la grieta). 
Nota 3. Usualmente, para la solución del segundo problema de elasticidad con valores en la fron¬ 
tera, el vector de desplazamiento se presenta en la forma del potencial de capa simple (denominado 
método indirecto.) Como resultado, el núcleo del operador en la ecuación integral del problema 
no tiene la singularidad fuerte como la del núcleo de la ecuación (3.14); pero el potencial de capa 
por lo que se debe definir un procedimiento de regularización para el cálculo de la integral en 
la ecuación (3.14) . Se demuestra en [20, 24], que la integral en la ecuación (3.14) se puede 
comprender en el sentido siguiente: 
T (x, x ' ) ‧ b (x') dΓ ' = v.p. T (x, x ' ) ‧ [b (x') - b (x)] dΓ'. (3.16) 
32 
simple no puede utilizarse para las simulaciones de grietas ya que es continuo sobre las fronteras 
de su distribución. 
Como se mencionó anteriormente, la ventaja de las ecuaciones (3.13, 3.14) es que las fronteras 
de grieta se pueden incluir como parte de la frontera total Γ del cuerpo. La solución única de la 
ecuación (3.14) existe si su lado derecho, f (x) , (el vector de las fuerzas de frontera) satisface las 
condiciones 
F0 = f (x) dΓ = 0, M0 = x x f (x) dΓ = 0. (3.18) 
(La fuerza total, F 0 , y el momento total, M 0 , de las fuerzas aplicadas f (x) deben ser iguales a 
cero. Las fronteras de grieta no participan en esta condición). 
3.2 Acción de los operadores del problema sobre las funciones de 
Gauss 
Para la solución numérica de la ecuación integral (3.14) utilizaremos la clase especial de funciones 
de aproximación descritas en 2. Consideremos una función finita suave u (x) dependiente de la 
coordenada x en la región — ∞ < x < ∞. Esta función se puede aproximar mediante las series 
siguientes 
|R (x)| ≥ (||u | | + | | u ' | | ) R 0 (D ) + | |u"|| (3.20) 
(3.19) 
Aquí, mh ( m = 0, ± 1 , ± 2 , . . . ) son las coordenadas de los nodos de aproximación, h es la distancia 
entre los nodos, D es un parámetro adimensional. Es posible demostrar, de acuerdo con [21, 14], 
que la siguiente evaluación para el error de aproximación se cumple, 
u (x) = uh (x) + R (x) , 
R0 (D ) = O ( e x p ( - π
2 D ) ) . 
32 
simple no puede utilizarse para las simulaciones de grietas ya que es continuo sobre las fronteras 
de su distribución. 
Como se mencionó anteriormente, la ventaja de las ecuaciones (3.13, 3.14) es que las fronteras 
de grieta se pueden incluir como parte de la frontera total Γ del cuerpo. La solución única de la 
ecuación (3.14) existe si su lado derecho, f (x) , (el vector de las fuerzas de frontera) satisface las 
condiciones 
F0 = f (x) dΓ = 0, M0 = x x f (x) dΓ = 0. (3.18) 
(La fuerza total, F 0 , y el momento total, M 0 , de las fuerzas aplicadas f (x) deben ser iguales a 
cero. Las fronteras de grieta no participan en esta condición). 
3.2 Acción de los operadores del problema sobre las funciones de 
Gauss 
Para la solución numérica de la ecuación integral (3.14) utilizaremos la clase especial de funciones 
de aproximación descritas en 2. Consideremos una función finita suave u (x) dependiente de la 
coordenada x en la región — ∞ < x < ∞. Esta función se puede aproximar mediante las series 
siguientes 
|R (x)| ≥ (||u | | + | | u ' | | ) R 0 (D ) + | |u"|| (3.20) 
(3.19) 
Aquí, mh ( m = 0, ± 1 , ± 2 , . . . ) son las coordenadas de los nodos de aproximación, h es la distancia 
entre los nodos, D es un parámetro adimensional. Es posible demostrar, de acuerdo con [21, 14], 
que la siguiente evaluación para el error de aproximación se cumple, 
u (x) = uh (x) + R (x) , 
R0 (D ) = O ( e x p ( - π
2 D ) ) . 
33 
Fig. 3 .1: R (x) de la aprox. Gaussiana de la ec.(3.21). D variable. 
Así, el error de aproximación R (x) , ecuación (3.20), depende de dos parámetros: la dispersión D 
y el paso de aproximación h. Si h es suficientemente pequeño y D = O (1) este error se puede 
considerar prácticamente despreciable. Consideremos la aproximación (3.19) en algunos casos 
particulares. 
P R I M E R CASO. Consideremos la función u (x) como el impulso unitario en el intervalo ̶1 < 
x < 1; es decir, 
(3.21) 
Las gráficas de loserrores de aproximación 
R ( x ) = u (x ) ̶ u h ( x ) , u h ( x ) = φ(x ̶ m h ) 
se presentan en la figura (3.1) para M = 50 y D = 0,2 (línea sólida), D = 2 (línea con triángulos), 
D = 20 (línea con puntos); y en la figura (3.2) para D = 2 y M = 20 (línea sólida), M = 40(línea 
con triángulos), M = 60 (línea con puntos). Aquí, 2 M + 1 es el número de puntos dentro del 
intervalo (̶1,1), h = 
En estas gráficas se puede observar que el error R (x) está concentrado cerca de los extremos 
del intervalo (̶1,1) y es mínimo para D D = 2. El máximo del error decrece con el paso de la 
aproximación h. 
34 
Fig. 3.2: R (x) de la aprox. Gaussiana de la ec.(3.21). D = 2. 
Error R (x) de la aproximación Gaussiana de la función impulso en la ecuación (3.21) 
S E G U N D O CASO. Consideremos la función u (x) definida por la ecuación 
u (x) = (3.22) 
Las gráficas de los errores de aproximación 
R (x) = u (x) — uh ( x ) , uh (x) = 
se presentan en la figura (3.3) para D = 2 y M = 20 (línea sólida), M = 40 (línea con triángulos), 
M = 60 (línea con puntos.) Como en el primer caso, el error es máximo en los extremos del 
intervalo (—1,1), pero debido a la continuidad de u (x ) , este error es menor que en el primer caso, 
y el mínimo del error se obtiene nuevamente para D = 2. 
El análisis detallado de la aproximación (3.19) se presenta en los trabajos de Maz 'ya [21, 14, 15]. 
Fig. 3.3: R ( x ) de la aprox. Gaussiana de la ec.(3.22). D = 2. 
(x — m h ) 
35 
donde h es la distancia entre nodos vecinos, con D = 2. Los vectores b ( i ) en la ecuación (3.24) 
son las incógnitas principales del problema. 
Si sustituímos la ecuación (3.24) en la ecuación (3.23), esta última se convierte en la suma de 
potenciales concentrados en las líneas tangentes γi en cada nodo de la frontera Γ. 
3.3 Discretización de las ecuaciones integrales del segundo problema 
de elasticidad. Caso plano 
(3.25) 
σ ( x, y) = [S (x — x', y — y') ‧ n (x', y')] ‧ b (x', y') Γ (x', y') dx'dy' ≈ 
I
( i ) ( x , y ) • b ( i ) , 
I
( i ) ( x , y ) = S (x - x', y - y') • n
( i ) φ i (x', y') dγ'. 
(3.26) 
(3.27) 
Consideremos un cuerpo isotrópico homogéneo que ocupa una área plana Ω con frontera Γ, como 
se muestra en la figura (3.4). La solución del segundo problema de elasticidad de valores en la 
frontera para esta área se puede determinar a partir de las ecuaciones (3.13, 3.14). La ecuación 
(3.13) para el tensor de esfuerzos σ (x, y) se puede reescribir de la forma siguiente: 
donde Γ (x, y) es la función delta concentrada sobre el contorno de frontera Γ, la integral se realiza 
sobre todo el espacio bidimensional. Para la construcción de la solución numérica de este prob-
lema, seleccionemos un conjunto de puntos ( x ( i ) , y(i)) sobre la frontera Γ, con distancias iguales 
entre ellos, y cambiemos el potential (3.23) concentrado sobre Γ por la suma de potenciales con-
centrados sobre las líneas tangentes γi en cada punto i, ver figura (3.4). Entonces, la distribución 
b (x, y) Γ (x, y) en la ecuación (3.23) es aproximada por la suma siguiente 
(3.23) 
(3.24) 
donde γi (x, y) es la función delta concentrada sobre la línea tangente γi. De acuerdo a la aproxi-
mación (3.19), la función φi (x) tiene la forma 
36 
Fig. 3.4: Bases global (e1 e 2 ) y local ( s
( i ) , n ( i ) ) de un área plana Ω con frontera Γ. 
Consideremos la integral que corresponde al nodo i-ésimo en el sistema de coordenadas local 
(s , z) de este nodo, como se muestra en la figura (3.4), 
I ( s , z ) = S ( s - s ' , z ) • n
( i ) φ (s') ds'. (3.28) 
Aquí, s es la coordenada a lo largo de la tangente del i-ésimo nodo, z es la coordenada a lo largo 
del eje paralelo al vector n ( i ) normal a Γ en el i-ésimo nodo. En este sistema, n(i)1 = 0 y n 2
( i ) = 1. 
La integral en la ecuación (3.28) se puede reescribir como la integral sobre todo el plano infinito 
(s , z) si se utiliza la función delta de Dirac, δ (z) , 
Iijl (s,z) = 
S i j k l (s̶s', z ̶ z') n k φ (s') δ (z') ds'dz'. (3.29) 
Debido a que esta es una integral del tipo convolución, se puede reescribir mediante transformadas 
de Fourier de las funciones del integrando en la forma 
Iijl (s,z) = 
S i j k l ( k 1 , k2) n k φ ( k 1 ) e
-i(k1s+k2z)dk1dk2. (3.30) 
Aquí, la transformada de Fourier S (k) del tensor S (x) esta definida en la ecuación (3.10); y para 
el caso de un medio isotrópico toma la forma (ver Apéndice A) 
S (k) = -4μ0N0 [ E 1 - 2 E 5 ( m ) + E 6 (m)] , k = ¦k¦. (3.31) m = 
37 
Los tensores E 1 , E 5 ( m ) y E 6 ( m ) en esta ecuación tienen los componentes siguientes (ver Apéndice 
A) 
E 1 i j k l = δi(kδj)l, E 5 i j k l ( m ) = m ( i δ i ) ( k m l ) , E 6 i j k l ( m ) = mimjmkml, (3.32) 
en donde los paréntesis en los subíndices significan simetrización. 
La transformada de Fourier (k 1 ) de la función φ (s) tiene la forma 
(3.33) 
Después de calcular la integral en la ecuación (3.29) con respecto a k 2 obtenemos 
(3.34) 
donde 
Aquí s ( i ) , n ( i ) son los vectores unitarios sobre los ejes s y z, respectivamente, y 
Aquí, i = Después de sustituir la ecuación (3.34) en la ecuación (3.30) e integrando con 
respecto a k1, obtenemos la siguiente expresión del tensor I 
I (s, z) = -4μ0N0 [J 1 1 (s , z) t11 + J12 (s , z) t 1 2 + J21 (s , z) t 2 1 + J22 (s , z) t 2 2 ] , (3.36) 
con 
y 
38 
Aquí, las cuatro funciones ji (ζ, η) tienen la forma de las integrales siguientes 
(3.37) 
Estas integrales dependen de las variables adimensionales (ζ, η) y simplemente pueden ser tabu­
ladas y mantenidas en la memoria de la computadora para cálculos posteriores. Para 
10 estas integrales pueden ser cambiadas por sus expresiones asimptóticas siguientes 
(3.38) 
i ) , n ( i ) ) son 
βi. 
Nótese que las ecuaciones siguientes conectan las funciones con la función de error 
especial Erf(z) 
Introduzcamos un sistema de coordenadas global con la base Cartesiana ( e 1 ; e 2 ) , y el sistema de 
coordenadas local, ( s ( i ) , n ( i ) ) , conectado con los nodos sobre la frontera Γ. ( s ( i ) , n ( i ) ) son vectores 
unitarios de la base local en el nodo x ( i ) , ver figura (3.4) 
En la base local ( s ( i ) , n ( i ) ) , el vector b ( i ) de la ecuación (3.24) tiene la forma 
(3.40) 
La conexión entre la base global y la base local en el i-ésimo nodo tiene la forma 
s(i) = cos(βi)e1 + sin(βi)e2 n(i) = - sin(βi)e1 + cos(βi)e2. (3.41) 
Aquí, βi es el ángulo entre los vectores unitarios e2 y n(i), figura (3.4).Las coordenadas de un punto arbitrario x(x, y) en la base local ( s (
x = x e 1 + y e 2 = r ( i ) s ( i ) + rs(i)n(i), (3.42) 
con 
rs(i) = (x - x(i)) cos(βi) + (y - y(i)) sin(βi), rn(i) = - (x - x(i)) sin(βi) + (y - y(i)) cos(βi). 
La conección entre las bases locales entre los nodos i-ésimo y j - é s i m o tiene la forma 
s(j) = cos( γji)s(i) + sin(γji)n(i), n(j) = - sin( γji)s(i) + cos( γji)n(i), γji = βj - 
(3.43) 
(3.39) 
39 
La expresión (3.26) para el tensor de esfuerzos en un punto arbitrario x del medio es aproximada 
por la ecuación siguiente 
σ ( x , y ) (x - x(i),y - y(i)) • b ( i ) = 
Usando las ecuaciones (3.43), el tensor de esfuerzos en la base global ( e 1 , e 2 ) tiene la forma 
siguiente 
donde, para k, l = 1, 2, 
El vector de fuerzas de frontera que actúa en el i-ésimo nodo tiene la forma siguiente 
f (x(i),y(i))= f (i); (3.46) 
donde 
En la base local del i-ésimo nodo, este vector tiene la forma 
f ( i ) = f s ( i ) s ( i ) + f n ( i ) n ( i ) , 
(3.45) 
(3.47) 
40 
f ( i ) = n ( i ) • σ ( x , y ) = n ( i ) • (S (x - x', y - y') • n (x', y')) • b (x', y') dΓ' 
• I (x - x(j),y - y(j)) • b ( j ) . 
y los componentes de A son 
donde los componentes de la fuerza de frontera son 
41 
(3.48) 
(3.49) 
r s ( i j ) = ( x ( i ) - x ( j ) ) cos (β j ) + ( y ( i ) - y(j)) sin (β j ) , 
r n ( i j ) = - ( x ( i ) - x(j)) s in (β j ) + ( y ( i ) - y(j)) cos ( β j ) . 
(3.50) 
El sistema de ecuaciones algebráicas líneales para los componentes de los vectores b (i) 
en la baselocal, se pueden obtener a partir de las condiciones de frontera de la ecuación (3.2) que 
serán satisfechas en todos los nodos (método de colocación). 
Introduzcamos un vector columna X para las incógnitas que están conectadas con los compo-
nentes por las relaciones 
X = ||Xj|| , j = 1, 2, 3 , . . . , 2N, 
X 2 i - 1 = bS(i), X2i = bn(i), i = 1, 2, 3,. . . ,M . 
(3.51) 
Aquí, M es el número total de puntos de frontera. El vector columna F define las fue rzas que 
actúan en los puntos de frontera 
F = ||Fj|| , j = 1, 2, 3 , . . . , 2M, 
F 2 i - 1 = fs(i), F2i = fn(i) i = 1, 2, 3, . . . ,M . 
(3.52) 
Aquí,(fs(i), fn(i)) son los valores de los componentes de las fuerzas (tangencial, fs, y normal, fn) 
aplicadas en los puntos de frontera, y que son conocidas a partir de las condiciones de frontera. 
Aquí, los argumentos γij- están definidos en la ecuación (3.43), y 
4 2 
Nótese que los componentes de la matriz B en las ecuaciones (3.54, 3.49) están expresadas me­
diante las funciones estandar (3.37), y el t iempo de computadora para sus cálculos es pequeño. 
Después de la solución de la ecuación (3.53), el tensor de esfuerzos en un punto arbitrario del 
cuerpo se puede determinar utilizando la ecuación (3.45). 
Aquí los componentes de la matriz B =||Bij|| están definidas por las ecuaciones siguientes 
(3.54) 
Las ecuaciones para la determinación del vector X tienen la forma 
(3.53) Bij Xj = Fi, j = 1, 2, 3 , . . . , 2 M . 
3.4 Ejemplos de solución numérica del problema 
Para evaluar la capacidad del método numérico propuesto, lo aplicaremos a los dos problemas 
siguientes que tienen solución exacta: 
1. Un plano elástico infinito con un orificio circular de radio unitario. 
2. Un disco elástico de radio unitario. 
P R I M E R CASO. Consideremos un plano elástico con un orificio circular de radio unitario (R = 
1), sujeto a un campo tensorial uniaxial constante, σ0, aplicado en el infinito a lo largo de la 
coordenada x, ver figura (3.5). 
σ0 = e1 x e1 (3.55) 
La condición de frontera en el borde Γ del orificio es 
n • σ|r = 0. (3.56) 
Nót
dian
Des
cue
Aqu
Las
4 2 
ese que los componentes de la matriz B en las ecuaciones (3.54, 3.49) están expresadas me­
te las funciones estandar (3.37), y el t iempo de computadora para sus cálculos es pequeño. 
pués de la solución de la ecuación (3.53), el tensor de esfuerzos en un punto arbitrario del 
rpo se puede determinar utilizando la ecuación (3.45). 
í los componentes de la matriz B =||Bij|| están definidas por las ecuaciones siguientes 
(3.54) 
 ecuaciones para la determinación del vector X tienen la forma 
(3.53) Bij Xj = Fi, j = 1, 2, 3 , . . . , 2 M . 
3.4 Ejemplos de solución numérica del problema 
Para evaluar la capacidad del método numérico propuesto, lo aplicaremos a los dos problemas 
siguientes que tienen solución exacta: 
1. Un plano elástico infinito con un orificio circular de radio unitario. 
2. Un disco elástico de radio unitario. 
P R I M E R CASO. Consideremos un plano elástico con un orificio circular de radio unitario (R = 
1), sujeto a un campo tensorial uniaxial constante, σ0, aplicado en el infinito a lo largo de la 
coordenada x, ver figura (3.5). 
σ0 = e1 x e1 (3.55) 
La condición de frontera en el borde Γ del orificio es 
n • σ|r = 0. (3.56) 
43 
El campo de esfuerzos en el medio se puede encontrar de forma similar a la ecuación (3.13) 
σ (x) = σ0 + [S (x - x ' ) • n (x')] • b (x') dΓ' . 
La ecuación integral para el vector b (x) se puede expresar a partir de las condiciones de frontera 
ecuación (3.2), y toma la forma 
(3.57) 
n (x) • [S(x - x ' ) • n (x')] • b (x ') dΓ' = - n (x) • σ0, x E Γ . (3.58) 
Nótese que en la base local (s, n ) , el lado derecho de esta ecuación tiene la forma 
Donde φ es el ángulo polar sobre el plano (x, y), ver figura (3.5). 
Fig. 3.5: Orificio de una placa plana infinita. 
Es posible mostrar que la solución exacta de la ecuación (3.58) tiene la forma 
b (φ) = bs (φ) s (φ) +bn (φ) n ( φ ) , (3.59) 
con componentes 
bs (φ) = 1/2 sin (2φ) , bn (φ) = ( c o s 2 φ - 1/4) . 
Entonces, la solución exacta del campo de esfuerzos en el medio está representado por la ecuación 
44 
Aquí, 1 es el tensor unitario de rango dos. 
En las figuras (3.6, 3.7), la solución exacta (3.59) (línea sólida) es comparada con la solución 
numérica de la ecuación (3.58), obtenida con el Método de Puntos de Frontera. La distribución de 
esfuerzos, σ11, a lo largo de la coordenada y (x = 0) se presenta en la figura (3.8). Nótese que la 
distribución exacta, σ11 (0, y) , que de acuerdo con la ecuación (3.60), tiene la forma 
y la línea sólida en la f i g u r a (3.8) corresponde a esta ecuación. 
Se puede observar que para M = 80, las soluciones numérica y exacta prácticamente coinciden. 
(3.61) 
Fig. 3.6: Solución exacta. Componente b s . 
S E G U N D O CASO. Un disco elástico de radio unitario (R = 1) sujeto a dos fuerzas concen¬ 
tradas aplicadas a lo largo de su diámetro, ver f i g u r a (3.9). La solución exacta de este problema, 
Fig. 3.7: Solución exacta. Componente b n. 
45 
Fig. 3.8: Esfuerzos en la dirección y, σ11. 
(3.62) 
donde los ángulos Φ1 y Φ2, y los intervalos r1 y r2 se muestran en la figura (3.9). 
La solución de este problema se encuentra a partir de la ecuación (3.23) que a continuación se 
repite 
σ (x) = [S (x - x ' ) • n (x')] • b (x ') dΓ' . 
considerando que | F | = 1, tiene la forma [25] 
y la ecuación para el vector b (x) se obtiene a partir de las condiciones de frontera, y tiene la 
forma integral siguiente 
n (x) • [S (x - x ' ) • n (x')] • b (x') dΓ' = - F δ (φ) + Fδ (φ - π) , x E Γ. (3.63) 
Para la solución numérica, las fuerzas concentradas aplicadas en los puntos de frontera para los 
ángulos φ = 0 y φ = π en la base local de estos puntos de frontera tienen la forma 
Fig. 3.9: Disco circular con R = 1, y cargas concentradas opuestas. 
46 
donde h es la distancia entre los puntos de frontera vecinos a lo largo de Γ Para los demás puntos 
de frontera, f ( i ) = 0. 
La distribución de esfuerzos σ11 (x, y) en diferentes secciones paralelas al eje del disco se presen­
tan en las figuras (3.10, 3.11). Las líneas sólidas representan a las soluciones exactas; las líneas 
interrumpidas corresponden a 60 puntos equidistantes sobre la frontera. Se puede observar que el 
error de la solución numérica es esencial solamente en una pequeña vecindad de los puntos de 
frontera sobre los que se encuentran aplicadas las fuerzas. 
(3.64) 
3.5 Solución del problema para cuerpos con grieta 
Consideremos la solución numérica del problema de elasticidad para una grieta recta que ocupa 
un intervalo (|x| < 1, y = 0) en un plano infinito. El plano está sujeto a un campo de esfuerzos 
y 
fo
P
á
d
d
L
ta
in
e
fr
la ecuación para el vector b (x) se obtiene a partir de las condiciones de frontera, y tiene la 
rma integral siguiente 
n (x) • [S (x - x ' ) • n (x')] • b (x') dΓ' = - F δ (φ) + Fδ (φ - π) , x E Γ. (3.63) 
ara la solución numérica, las fuerzas concentradas aplicadas en los puntos de frontera para los 
ngulos φ = 0 y φ = π en la base local de estos puntos de frontera tienen la forma 
Fig. 3.9: Disco circular con R = 1, y cargas concentradas opuestas. 
46 
onde h es la distancia entre los puntos de frontera vecinos a lo largo de Γ Para los demás puntos 
e frontera, f ( i ) = 0. 
a distribución de esfuerzos σ11 (x, y) en diferentes secciones paralelas al eje del disco se presen­
n en las figuras (3.10, 3.11). Las líneas sólidas representan a las soluciones exactas; las líneas 
terrumpidas corresponden a 60 puntos equidistantes sobre la frontera. Se puede observar que el 
rror de la solución numérica es esencial solamente en una pequeña vecindad de los puntos de 
ontera sobre los que se encuentran aplicadas las fuerzas. 
(3.64) 
3.5 Solución del problema para cuerpos con grieta 
Consideremos la solución numérica del problema de elasticidad para una grieta recta que ocupa 
unintervalo (|x| < 1, y = 0) en un plano infinito. El plano está sujeto a un campo de esfuerzos 
Fig. 3.10: Distribución de σ11(x) en el disco sujeto a dos fuerzas de compresión concentradas. 
Fig. 3.11: Distribución de σ11 (y) en el disco sujeto a dos fuerzas de compresión concentradas. 
σ0(x) en el infinito. Encontremos la solución de este problema de forma similar a la ecuación 
(3.57), y la ecuación integral de este problema toma la forma de la ecuación (3.58), donde Γ es la 
línea de la grieta. 
Para un campo de esfuerzos constante dirigido a lo largo del eje y 
μ0 = 1 , N0 = 1 , la solución de esta ecuación toma la forma 
y tomando 
b0 (x) = b0 (x) n , b0 (x) = (3.65) 
47 
donde n = e 2 es el vector normal a la línea de la grieta. 
El error R (x) = b n (x) — b0 (x) de la solución numérica para varias cantidades de puntos de 
frontera sobre la línea de la grieta se presenta en la f i g u r a (3.12).. En esta gráfica se puede ver 
que el máximo del error está concentrado en las vecindades de las puntas de la grieta. Pero las 
asimptóticas de la solución cerca de estas puntas nos proporcionan información importante para 
aplicaciones de ingeniería. Es bien sabido [20, 24] que el brinco exacto del vector desplazamiento, 
Fig. 3.12: Error R ( x ) de la sol. numérica de una grieta recta en un plano infinito. 
b ( x ) , sobre la grieta tiene la siguiente ecuación asimptótica cerca de la punta de la grieta 
Aquí r es la distancia entre el punto x E Γ y la punta de la grieta, como se ve en la figura (3.13). 
(3.66) 
48 
Fig. 3.13: Sistema de coordenadas local en la punta de la grieta. 
Los componentes ( β s , β n ) del vector β en la ecuación (3.66) están conectadas con los factores de 
intensidad de esfuerzos KI y KII mediante las ecuaciones [20, 24] 
y las asimptóticas de los componentes σ22 y σ 1 2 del tensor de esfuerzos cerca de la punta de la 
grieta tienen las siguientes formas conocidas 
(3.67) 
(3.68) 
Para calcular los coeficientes del vector β en la ecuación (3.66) y los factores de intensidad de 
esfuerzos con alta precisión, tenemos que m o d i i c a r el método desarrollado previamente. Esta 
modificación está basada en el t e o r e m a de conservatividad polinomial [24] que aplicaremos a 
continuación. Consideremos una grieta recta (|x| < l, y = 0) en un plano infinito. El teorema nos 
dice que si un campo externo σ0 (x, y) es polinomial sobre la línea de la grieta como sigue 
49 
(3.69) σ0(x, 0) = σ 1 + σ 2 x + σ 3 x 2 + • • • + σ m x m - 1 , |x | < l, 
( σ k = constante, k = 1, 2 , . . . , m ) , entonces el vector de la apertura de la grieta b (x) toma la 
forma siguiente 
B m (x) = a 1 + a 2 x + a 3 x 2 + • • • + a m x m - 1 , (3.70) 
donde los a k (k = 1, 2 , . . ., m ) son vectores constantes y B m (x) es una función polinomial de la 
misma potencia que el campo externo ao(x, 0). 
Este teorema nos permite encontrar la solución de la ecuación integral del problema de la grieta 
de forma similar a la ecuación (3.70). Para hallar la solución, introduzcamos un conjunto de nodos 
auxiliares sobre la línea de la grieta x a ( j ) ( j = 1, 2 , . . . . M a ) y utilicemos estos nodos para realizar 
la aproximación Gaussiana (3.19) de la función b (x) en la ecuación (3.70) 
Aquí, h a es la distancia entre los nodos auxiliares. Para una buena aproximación de b (x) , el 
número M a debe ser suficientemente grande ( M a > 50). Después de sustituir la ecuación (3.71) de 
b (x) en la ecuación (3.57) de los componentes del tensor de esfuerzos, obtenemos las ecuaciones 
siguientes para r, s = 1, 2: 
(3.72) 
(3.71) 
Donde (as(k),an(k)) son los componentes de los vectores a k de la ecuación (3.70) en la base local 
(s, n ) , es decir, 
ak = a S ( k ) s + a n ( k ) n . 
Los componentes de las fuerzas que actúan sobre la línea de la grieta tienen las ecuaciones 
y si los lados de la grieta están libres de esfuerzos, entonces tenemos 
fs (x) = fn (x) = 0, |x | ≥ l. 
Para obtener las 2 m incógnitas aS(k) y an(k) de la ecuación (3.71) tenemos que satisfacer las ecua­
ciones (3.73) y (3.74) en m puntos de frontera sobre la línea de la grieta. La cantidad de estos 
puntos de frontera coinciden con la potencia máxima de la función polinomial de la ecuación 
(3.70). En realidad, esta potencia no es muy grande y podemos escoger las coordenadas de estos 
puntos de frontera en la región interior de la grieta, donde el error de la aproximación Gaussiana 
es pequeño. El sistema de ecuaciones final para toma la forma 
(3.74) 
(3.75) 
donde las x ( 1 ) son las coordenadas de los puntos de frontera sobre la línea de la grieta. 
Nótese que las asimptóticas de la expresión (3.70) cerca de las puntas de la grieta, x = ± l , tienen 
la forma 
(3.76a) 
Consideremos la solución numérica del sistema (3.75) cuando el campo aplicado es de tensión 
constante a lo largo del eje y, es decir, 
σ011 = σ012= 0, σ022 = 1. 
50 
51 
Para m = 5, M a = 50, l = 1 y una distribución homogénea de los puntos de frontera, 
la solución del sistema (3.75) para k = 1 , . . . , 5, es 
Tomando en cuenta que la solución exacta de este problema tiene la forma de la ecuación (3.65), 
y por lo tanto 
podemos ver que el error de la solución numérica es menor que el 1 % en este caso. 
Ahora consideremos una área rectangular con una grieta recta central sometida a tensión constante 
en la dirección y, ver figura (3.14).La solución de este problema se puede encontrar en la forma 
Fig. 3.14: Área rectangular con grieta central. 
donde la primera integral se reliza sobre la línea de la grieta Γc y la segunda integral sobre el borde 
externo del área Γb. 
(3.77) 
as(k) = 0 , 
a n ( 1 ) = 0,9912, a n ( 2 ) = 0, a n ( 3 ) = 0,0037, a n ( 4 ) = 0, a n ( 5 ) = 0,0048 . 
as (k)= 0, 
a n ( 1 ) = 1, an(l) = 0, 
k = 1 2 3 . . . ; 
l = 2, 3 , . . . , 
a ( x, y) = [S (x - x ' , y - y') • n (x' , y')] • b c (x' , y') dΓ' + 
[S (x - x ' , y - y') • n (x ' ,y ' ) ] • b (x ' ,y ' ) dΓ', 
Las fuerzas 
f (x) = n (x) • σ (x) , x E Γ. 
sobre la frontera del cuerpo y las fronteras de la grieta son calculadas a partir de la ecuación (3.78). 
Los valores de estas fuerzas en los puntos de frontera tienen la forma siguiente en las bases locales 
están definidas en la ecuación (3.48). donde las funciones 
donde 
(3.78) 
Utilizando las ecuaciones (3.70, 3.71) para la representación del vector b c (x) hemos escogido 
un conjunto de nodos auxiliares x a ( j ) ( j = 1, 2 , . . . , M a ) sobre la línea de la grieta, así como la 
potencia m — 1 de la función polinomial B m (x) en la ecuación (3.70). Entonces, se deben definir 
las coordenadas de los m puntos de frontera sobre la línea de la grieta y las condiciones de frontera 
se deben satisfacer en estos puntos. Enumeremos los puntos de frontera sobre la línea de la grieta 
desde 1 hasta m. También tenemos que definir los puntos de frontera sobre la frontera exterior del 
cuerpo donde el vector b (x) es aproximado mediante las ecuaciones (3.24, 3.71). La numeración 
de estos puntos de frontera se haces desde m + 1 hasta M , donde M es el número total de puntos 
de frontera, o nodos principales. Como resultado, el tensor de esfuerzos a (x, y) en la base global 
( e 1 , e 2 ) es aproximado por las ecuaciones obtenidas a partir de las ecuaciones (3.72, 3.45, 3.48) 
52 
( s ( i ) , n ( i ) ) : 
f ( x ( i ) , y ( i ) ) = f s ( i ) s ( i ) + f n ( i ) n ( i ) , (3.79) 
53 
donde la función B0 (x) está definida en la ecuación (3.71) y las funciones A l k S S, AlkSn,AlknS,y Alknn 
están dadas en la ecuación (3.49). 
El sistema f i n a l de coeic ientes desconocidos que definen la solución b (x) sobre la 
línea de la grieta, y los coeic ientes en los puntos sobre la frontera exterior del área, se 
obtienen a partir de las condiciones de frontera y toman la forma 
BcklXl = F k , k = 1, 2 , . . . , 2 M . (3.80) 
aquí el vector X de incógnitas tiene los componentes siguientes: 
X 2 l - 1 = aS (k), X 2 1 = an(k) l =1 2 , . . . , m ; (3 81) 
X 2 l - 1 = bS(k), X21 = bn(k), l = m + 1 ,m + 2 , . . . , M . 
El lado derecho, F , de la ecuación (3.80) es: 
F 2 k - 1 = f s ( k ) , F2k = f n ( k ) , k = 1 , 2 , . . . , M . (3.82) 
54 
Los elementos de la matriz B c toman la forma 
(3.83) 
La solución de la ecuación (3.80) se obtuvo para varios tamaños del rectángulo y diferentes lon¬ 
gitudes de grieta. En las figuras (3.15, 3.16), las gráficas de las funciones σ 1 1(x,0) y σ22(x,0), 
respectivamente, son presentadas para a = 2, b = 3, y l = 1. El número de nodos auxiliares so­
bre la línea de la grieta se escogió como M a = 70, el número de puntos de frontera fue m = 5 y 
M = 27, 47 y 67. El incremento del número de nodos auxiliares y de puntos de frontera respecto 
a los indicados no cambió la solución. 
Fig. 3.15: Distribución del esfuerzo σ11(x, 0) a lo largo de la grieta de un área rectangular. 
Fig. 3.16: Distribución del esfuerzo σ22(x, 0) a lo largo de la grieta de un área rectangular. 
El resultado de los cálculos de los factores de intensidad de esfuerzos en las puntas de la grieta se 
presentan en la figura (3.17);en esta figura 
R c 2 l - 1 , 2 k - 1 = C l k s s , B c 2 l - 1 , 2 k = Clksn, 
B c 2 l , 2 k - 1 = Clkns, Bc2l,2k = Clknn l.k=1,2,...,m; 
Bc2l-1,2k-1= Alkss, Bc2l-1,2k= Alksn 
B c 2 l , 2 k - 1 = Alkns, Bc2l,2k= Alknn, l, k = m + 1, m + 2,...,M. 
55 
ω m Ma M 
0,5 5 70 105 
1,5 11 70 71 
Se puede ver que una distribución homogénea de los puntos de frontera sobre la línea de la gri¬ 
eta lleva a obtener errores relativamente grandes en el cálculo de los factores de intensidad de 
esfuerzos. Este error disminuye si la densidad de puntos de frontera es mayor en los extremos de 
la grieta. Los experimentos numéricos muestran que la posición óptima de los puntos de frontera 
x ( i ) sobre la línea de la grieta corresponde a las raices de los polinomios de Chebyshev; es decir, 
x ( i ) = l cos , i= 1 , 2 , . . . , m . 
Los resultados numéricos en la f i g u r a (3.17) fueron obtenidos para estas coordenadas de los puntos 
de frontera sobre la línea de la grieta. 
Fig. 3.17: Factores de intensidad de esfuerzos para una grieta central en un área plana rectangular. 
donde Bmn (x) es la componente normal del vector B m (x) definido en la ecuación (3.70), y 
son los factores de intensidad de esfuerzos KI y los coeficentes βn que coresponden a una 
placa infinita con grieta de longitud 1. En este caso, considerando μo = 1 y No = 1. Las 
líneas sólidas en la figura (3.17) son soluciones numéricas presentadas en [26] (la precisión es del 
1%), las líneas discontinuas son las soluciones numéricas obtenidas con el método desarrollado. 
Para los valores de ω indicados, se consideraron las siguientes cantidades de puntos de frontera y 
de nodos auxiliares: 
56 
3.6 Conclusiones 
Los métodos numéricos desarrollados son una herramienta efectiva para la solución del segundo 
problema con valores en la frontera de elasticidad. La precisión del método depende de la densidad 
de los puntos de frontera en la frontera del cuerpo. 
En esta versión desarrollada para este método, las distancias entre puntos de frontera vecinos se 
consideraron del mismo tamaño. Estas distancias deben ser de alrededor de 0,1R c para obtener 
una precisión aceptable de la solución ( ~ 1 %) . 
Aquí, R c es el radio mínimo de curvatura de la frontera del cuerpo u otra longitud característica. 
Nótese que en la vecindad de los puntos de fronteras anguladas o en los puntos de aplicación de 
las fuerzas concentradas, la función desconocida b (x) debe ser aproximada de acuerdo a la forma 
analítica de las asimptóticas de la solución exacta en estas regiones (ver sección 3.5). Esta modifi¬ 
cación incrementa esencialmente la precisión del método. El método también permite considerar 
distribuciones no homogéneas de puntos de frontera sobre la super ic ie del cuerpo considerado. 
57 
4.1 Diferentes formas de la teoría de plasticidad 
En esta sección describiremos la c lasi icación y descripción de las diferentes formas de la teoría 
de plasticidad que se presenta en ANSYS. 
El comportamiento no lineal de los materiales cuando están sujetos a diferentes tipos de cargas, se 
deben a la relación no lineal entre los esfuerzos y las deformaciones, es decir, el esfuerzo es una 
función no lineal de la deformación. Esta relación también es dependiente de la trayectoria, por lo 
que el esfuerzo depende tanto de la historia de la deformación como de la deformación misma. 
En la actualidad existen diferentes estudios sobre las no linealidades de materiales estructurales 
tales como plasticidad, elasticidad no lineal, hiperelasticidad, viscoelasticidad, y otros. 
Los métodos desarrollados para determinar las relaciones no lineales entre los esfuerzos y las 
deformaciones se pueden resumir como sigue: 
• Plasticidad independiente del tiempo. Está caracterizado por la deformación instantanea ir¬ 
reversible que ocurre en un material. 
• Plasticidad dependiente del tiempo. Ocurre cuando las deformaciones plásticas se desarrollan 
durante un intervalo de tiempo; también denominada viscoplasticidad. 
• Ondulamiento. Es una deformación irreversible que ocurre en un material y es dependiente 
del tiempo. El marco del t iempo para el ondulamiento generalmente es mayor que en el caso 
de plasticidad dependiente del t iempo. 
• Elasticidad no lineal. Se presenta en situaciones en donde la relación esfuerzo-deformación 
es no lineal. Además todas las deformaciones son reversibles. 
• Viscoelasticidad. Es la caracterización de propiedades del material que dependen del t iempo 
y que incluye una contribución viscosa a la deformación elástica. 
4 Problema elasto-plástico para cuerpos en 2D 
58 
‧ Concreto. Los materiales similares al concreto incluyen formación de grietas y desmoron¬ 
amiento. 
‧ Expansión. Aquí, los materiales se elongan debido a la presencia de flujo de neutrones. 
Definiciones de la deformación 
Para el caso de materiales no lineales, la definición de deformación elástica tiene la forma sigu¬ 
iente: 
εel = ε - εpl - εt - εo - εex (4.1) 
que nos dice que la deformación elástica se puede obtener a partir de la diferencia de la deforma¬ 
ción total con las demás deformaciones (pl=plástica, t=térmica, o=ondulamiento, ex=expansión). 
La teoría de plasticidad ofrece relaciones matemáticas que caracterizan la respuesta elastoplástica 
de los materiales. Para el presente trabajo, sólo consideraremos la relación entre las deformaciones 
total, elástica y plástica y se describirá el modelo que fue utilizado para comparar las simulaciones 
realizadas con el Método de Puntos de Frontera. 
4.1.1 Plasticidad independiente del tiempo o de la razón (rate independent 
plasticity) 
La plasticidad independiente del t iempo está caracterizada por la deformación irreversible que 
ocurre en un material una vez que cierto nivel de esfuerzos es alcanzado. Se considera que la 
deformación plástica se desarrolla instantáneamente, es decir, independiente del tiempo. La forma 
en que se presenta la deformación plástica se puede clasificar o caracterizar de acuerdo a los 
diferentes tipos de comportamiento del material como sigue: 
‧ Endurecimiento isotrópico no lineal 
‧ Endurecimiento isotrópico bilineal 
‧ Endurecimiento isotrópico multilineal 
Endurecimiento cinemático bilineal clásico 
59 
• Endurecimiento cinemático multilineal 
• Endurecimiento cinemático no lineal 
4.1.2 Teoría de la plasticidad 
La teoría de la plasticidad proporciona relaciones matemáticas que caracterizan la respuesta elasto-
plástica de los materiales. Hay tres ingredientes en la teoría de plasticidad independiente del tiem-
po: el criterio de cedencia, la regla de flujo y la regla de endurecimiento. 
Criterio de cedencia 
El criterio de cedencia determina el nivel de esfuerzos en el que se inicia la deformación plástica 
o cedencia. Para esfuerzos multicomponentes, este está representado