CARACTERIZACIÓN DE AGUJEROS CIRCULARES
EN PLACAS PLANAS INFINITAS A PARTIR DE LA
TEORÍA CLÁSICA DE ELASTICIDAD
Proyecto de grado de Ingeniería Mecánica
Felipe Silva Ardila
200210869
Alejandro Marañon
Asesor
Facultad de Ingeniería
Departamento de Ingeniería Mecánica
Universidad de Los Andes
Bogotá, D.C. Junio 16 del 2008
Tabla de Contenidos
1. Introducción 1
2. Objetivos 3
2.1. Definición del problema: . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3. Teoría Clásica de Elasticidad 5
3.1. Ley de Hook Generalizada . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2. Ecuaciones de Equilibrio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.3. Función de esfuerzos de Airy . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.4. Solución del problema 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.5. Sistema de coordenadas polares - Matriz de rotación . . . . . . . . . . . . . 10
3.6. Solución del Problema Planteado . . . . . . . . . . . . . . . . . . . . . . . . 11
4. Elementos Finitos - Simulaciones 15
4.1. El procedimiento de Elementos Finitos . . . . . . . . . . . . . . . . . . . . . 16
4.2. Simulaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.2.1. Elementos a utilizar y suposiciones . . . . . . . . . . . . . . . . . . 16
4.3. Resultados de las simulaciones . . . . . . . . . . . . . . . . . . . . . . . . . 20
5. Algoritmos Genéticos 22
I
5.1. ¿Qué son los Algoritmos Genéticos? . . . . . . . . . . . . . . . . . . . . . . 22
5.2. Implementación del Algoritmo Genético para la solución del Problema Plantea-
do . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.2.1. Momentos Geométricos . . . . . . . . . . . . . . . . . . . . . . . . 23
5.2.2. Componentes del Algoritmo Genético . . . . . . . . . . . . . . . . . 25
5.2.3. Etapas del algoritmo genético . . . . . . . . . . . . . . . . . . . . . 27
6. Resultados 30
6.1. Resultados del Algoritmo Genético . . . . . . . . . . . . . . . . . . . . . . . 31
7. Conclusiones 36
8. Apéndice A 38
8.1. Introducción y descripción de propiedades mecánicas . . . . . . . . . . . . . 38
8.2. El análisis de esfuerzos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
8.3. El análisis de deformación . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
9. Apéndice B 47
9.1. Otros resultados de las simulaciones . . . . . . . . . . . . . . . . . . . . . . 47
10. Apéndice C 49
10.1. Ventana experimental . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
11. Apéndice D 51
11.1. Algoritmo Genético . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
12. Bibliografía 58
II
Lista de Figuras
2.1. Definición del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.1. Fuerzas actuando en un cuerpo arbitrario. . . . . . . . . . . . . . . . . . . . 6
3.2. Placa infinita con agujero sometida a esfuerzos uniformes. . . . . . . . . . . 12
4.1. Elemento PLANE2 de Ansys utilizado en las simulaciones. . . . . . . . . . . 17
4.2. Diagrama Esfuerzo vs Deformación unitaria. Ensayo de tensión Universal . . 18
4.3. Análisis estadístico. ANOVA . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.4. Simulación placa completa. Estado de deformación enx (arriba), Estado de
deformación eny (abajo) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.5. Mallado final. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.6. Estado de deformaciones en la dirección x (izq) y en la dirección y (der). . . . 20
5.1. Resultados de población inicial, padres y nueva generación en una etapa pri-
maria. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5.2. Etapas del Algoritmo Genético. . . . . . . . . . . . . . . . . . . . . . . . . . 29
6.1. Comportamiento convergente de la función desempeño. . . . . . . . . . . . . 31
III
6.2. Comportamiento convergente de la función desempeño. Las lineas rojas rep-
resetan los valores reales utilizados en las simulaciones para obtener las de-
formaciones experimentales. Para ver otros resultados remitirse al apéndice
B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.3. Promedios de todos los parametros despues de 15 simulaciones. . . . . . . . 33
6.4. Desviacion estándar de todos los parametros despúes de 15 simulaciones. . . 33
6.5. Coeficiente de variacion de todos los parametros despúes de 15 simulaciones. 34
6.6. Resultados del algoritmo para 15 simulaciones con variaciones sobre el parámetro
k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
6.7. Resultados del algoritmo para 20 simulaciones. Parametrok = 0,55 . . . . . 35
6.8. Resultados del algoritmo para 20 simulaciones. Parametrok = 0,65 . . . . . 35
8.1. Diagrama de esfuerzos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
8.2. Fuerzas actuando en un elemento tetrahedríco. . . . . . . . . . . . . . . . . . 43
8.3. Deformaciones axiales y cortantes. . . . . . . . . . . . . . . . . . . . . . . . 45
9.4. Comportamiento convergente de la función desempeño. Las lineas rojas rep-
resetan los valores reales utilizados en las simulaciones para obtener las de-
formaciones experimentales. . . . . . . . . . . . . . . . . . . . . . . . . . . 47
9.5. Comportamiento convergente de la función desempeño. Las lineas rojas rep-
resetan los valores reales utilizados en las simulaciones para obtener las de-
formaciones experimentales. . . . . . . . . . . . . . . . . . . . . . . . . . . 48
IV
Capítulo 1
Introducción
Dentro de las técnicas de mantenimiento es importante poder caracterizar agujeros pre-
sentes en elementos mecánicos para evitar situaciones desfavorables; en la actualidad no
siempre es posible caracterizar dichos agujeros sin pretender destruir la pieza a tratar. Es así
como este proyecto de grado pretende caracterizar agujeros en placas a partir de los resultados
de pruebas no destructivas o simulaciones junto con la teoría clásica de elasticidad.
Para modelar cualquier estructura (ej. placas), se puede emplear un método muy bien
conocido y desarrollado a lo largo de los últimos 50 años, el de elementos finitos, que es capaz
de mostrar el estado de deformaciones de una pieza, sujeta a un estado de cargas determinado.
Si en un modelo de este tipo se introduce una pieza con fallas en este caso un agujero, éste
es capaz de cálcular las deformaciones y modelar la respuesta de una manera confiable. El
objetivo de esta tesis es solucionar el problema a la inversa, es decir, si se conoce cual es
el estado de deformaciones de una pieza, su geometría externa y además las cargas, ¿Sería
posible determinar la presencia de un agujero, así como su ubicación y tamaño?.
Si la respuesta a la pregunta fuera afirmativa, y si el modelo matemático y computacional
entregaran respuestas confiables, tendríamos un método de caracterización de agujeros que a
1
pesar de lo sofisticado, sería poco costoso, entregando una opción a la industria para deter-
minación de falla con una prueba no destructiva y confiable.
Para dicho fin se pretende entonces optimizar una función objetivo, restringida a las condi-
ciones de frontera y que como parámetro de entrada reciba la respuesta o estado de deforma-
ciones, ya sea simulada o experimental. Para caracterizar el agujero se pretende emplear
algoritmos genéticos ya que éste método permite encontrar óptimos globales de la función
objetivo sin mayores complicaciones.
2
Capítulo 2
Objetivos
2.1. Definición del problema:
Si se supone la existencia de una placa infinita con un agujero a la cual se aplica un esfuer-
zo axial constante. Al asumir que ésta tiene un comportamiento elástico, lineal e isotrópico,
es posible plantear el problema a examinar en este Proyecto de Grado de la siguiente forma:
Una vez dado el campo de desplazamientos o deformaciones medidos por medio de técnicas
de interferometría o simulaciones (según sea el caso), ¿es posible encontrar la ubicación y el
diámetro del agujero presente en la placa? (Fig. 2.1).
Para dicho fin, se busca plantear un modelo analítico matemático en un material elástico
e isotrópico usando la teoría clásica de elasticidad. Así mismo, usar la teoría de problemas
inversos para el planteamiento de una función objetivo a optimizar, por medio de algoritmos
genéticos.
3
Fig. 2.1: Definición del problema
Objetivos Específicos:
Estudiar y comprender la teoría clásica de elasticidad para estructuras en forma de
placas infinitas con agujeros.
Estudiar la definición y solución de problemas inversos basados en la teoría de elasti-
cidad para poder plantear un modelo matemático.
Obtener los desplazamientos por medio de experimentos de interferometría o simula-
ciones computacionales.
Plantar un algoritmo genético que permita solucionar el problema inverso especificado.
Realizar comparaciones entre teoría, experimentos y simulaciones computacionales.
4
Capítulo 3
Teoría Clásica de Elasticidad
En este capítulo, se presentarán definiciones y resultados importantes de propiedades
mecánicas necesarias para la solución del problema planteado anteriormente [1], [2] y [3].
Para una mayor introducción de las propiedades mecánicas favor remitirse al apéndice A.
3.1. Ley de Hook Generalizada
Cuando se tiene un estado de esfuerzos como el mostrado en la figura 8.1, el esfuerzoσx
genera una elongación en la direcciónx, y la deformacióne′x en la direcciónx es:
e′x =σx
E. (3.1)
Asi mismo,σy genera una contracción en la direcciónx lo cual genera una deformacióne′′x
con magnitud
e′′x = −ν σy
E, (3.2)
de la misma formaσz genera
e′′′x = −ν σz
E. (3.3)
5
Si estas tres deformaciones se superponen y si se usa el argumento anterior para las deforma-
ciones eny y z se tiene la ley de Hook Generelazida así:
ex =1
E[σx − ν (σy + σz)] ,
ey =1
E[σy − ν (σx + σz)] ,
ez =1
E[σz − ν (σx + σy)] . (3.4)
3.2. Ecuaciones de Equilibrio
Consideremos un cuerpo en equilibrio con volumenV y superficieS. Supongamos que la
región tiene una distribución de fuerzasF y tracciones superficialesT n tal y como lo muestra
la figura 3.1. Para el quilibrio estático es necesario que se conserve el momentum lineal, y
Fig. 3.1: Fuerzas actuando en un cuerpo arbitrario.
esto implica que las fuerzas que actúan sobre el cuerpo sean balanceadas, por ende la fuerza
resultante debe ser cero. Lo anterior se puede escribir matemáticamente así:∫∫s
T ni dS +
∫∫∫V
FidV = 0. (3.5)
6
Usando las ecuaciones (8.9) para el vector de tracciones, se puede rescribir la ecuación (3.5)
en terminos de esfuerzos: ∫∫s
σjinjdS +
∫∫∫V
FidV = 0. (3.6)
Aplicando el teorema de la divergencia visto en cálculo vectorial [8] sobre la integral de
superficie, tenemos que: ∫∫∫s
(σji,j + Fi) dV = 0. (3.7)
Puesto que la regiónV es arbitraria y el integrando (3.7) es continuo se tiene, por el teorema
del valor cero [1],[8]
σji,j + Fi = 0. (3.8)
Este resultado representa tres relaciones escalares llamadas lasecuaciones de equilibrio. Es-
critas en notación escalar son:
∂σx
∂x+∂τyx
∂y+∂τzx
∂z+ Fx = 0,
∂τxy
∂x+∂σy
∂y+∂τzy
∂z+ Fy = 0,
∂τxz
∂x+∂τyz
∂y+∂σz
∂z+ Fz = 0. (3.9)
3.3. Función de esfuerzos de Airy
Muchas soluciones para problemas de deformación y/o de esfuerzos en el plano pueden
ser determinadas por medio de una función de esfuerzos particular. Este método usa lafunción
de esfuerzos de Airyy su función principal es reducir la formulación general a una ecuación
que depende de sólo una incognita. La ecuación resultante puede ser resuelta por diferentes
métodos de las matemáticas aplicadas. Supongamos que las fuerzas pueden ser derivables de
una funciónpotencialV tal que
7
Fx = −∂V∂x
, Fy = −∂V∂y
. (3.10)
Esta suposición no es muy absurda puesto que muchas fuerzas encontradas en aplicaciones
como la fuerza de gravedad, caen dentro de esta categoria. Bajo las ecuaciones (3.10) las
condiciones de equilibrio 3.9 se pueden escribir asi:
∂(σx − V )
∂x+∂τxy
∂y= 0,
∂τxy
∂x+∂(σy − V )
∂y= 0. (3.11)
Se puede observar que estas ecuaciones se satisfacen si se escoge la siguiente represetación
para los esfuerzos
σx =∂2φ
∂y2+ V,
σy =∂2φ
∂x2+ V,
τxy = − ∂2φ
∂x∂y. (3.12)
Dondeφ = φ(x, y) es una forma arbitraria llamadaFunción de esfuerzos de Airy.Sustituyen-
do (3.12) en las ecueciones de equilibrio (3.9) en términos de desplazamientos [1] tenemos:
∂4φ
∂x4+ 2
∂4φ
∂x2∂y2+∂4φ
∂y4= ∇4φ = 0. (3.13)
Aquí el operador∇4 = ∇2∇2 es llamado el operadorbiarmónico. La ecuación biarmónica
en términos de coordenadas polares es ([1] pag 134):
∇4 =
(∂2
∂r2+
1
r
∂
∂r+
1
r2
∂2
∂θ2
) (∂2
∂r2+
1
r
∂
∂r+
1
r2
∂2
∂θ2
)φ = 0. (3.14)
8
3.4. Solución del problema 2D
Usando la función biarmónica de Airy debemos encontrar una solución a la siguiente
ecuación:
∇4 =
(∂2
∂r2+
1
r
∂
∂r+
1
r2
∂2
∂θ2
)2
φ = 0. (3.15)
Primero se debe encontrar una solución general suponiendo queφ(r, θ) es separable de la
forma φ(r, θ) = f(r)ebθ, dondeb es un parametro a determinar. Sustituyendo en (3.15) y
cancelando los terminosebθ tenemos
f ′′′′ +2
rf ′′′ − 1− 2b2
r2f ′′ +
1− 2b2
r3f ′ +
b2(4 + b2)
r4f = 0. (3.16)
Para resolver esta ecuación diferencial, se puede hacer el siguiente cambio de variabler = eξ,
lo cual transforma la ecuacion (3.16) en una ecuación diferencial con coeficientes constantes
f ′′′′ − 4f ′′′ + (4 + 2b2)f ′′ − 4b2f ′ + b2(4 + b2)f = 0. (3.17)
donde las primas denotand/dξ. La ecuación característica para la sustituciónf = eaξ de la
ecuación diferencial es
(a2 + b2)(a2 − 4a+ 4 + b2) = 0. (3.18)
Las raíces de esta ecuación son
a = ±ib, a = 2± ib
or (3.19)
b = ±ia, b = ±i(a− 2).
Debemos considerar solo soluciones periódicas enθ, y estas se obtienen escogiendob = in,
donden es un entero. Para un valor particular den, ocurren raíces repetidas, y esto implica
9
considereaciones especiales a la hora de desarrollar la solución final. Detalles de la solu-
ción completa son acreditados a Michell (1899) [1] (pag 158). La forma final de la solución
(Solución de Michell) se puede escribir asi
φ = a0 + a1 log r + a2r2 + a3r
2 log r
+ (a4 + a5 log r + a6r2 + a7r
2 log r)θ
+ (a11r + a12r log r +a13
r+ a14r
3 + a15rθ + a16rθ log r) cos θ
+ (b11r + b12r log r +b13r
+ b14r3 + b15rθ + b16rθ log r) sin θ (3.20)
+∞∑
n=2
(an1rn + an2r
2+n + an3r−n + an4r
2−n) cosnθ
+∞∑
n=2
(bn1rn + bn2r
2+n + bn3r−n + bn4r
2−n) sinnθ.
3.5. Sistema de coordenadas polares - Matriz de rotación
Ya que el problema planteado tiene una geometria circular debido al agujero presente
en la placa, es más comodo trabajar en un sistemas de coordenadas polares. Para esto es
necesario realizar un pequeño cambio de variables de la siguiente manera:
x = r cos(θ), y = r sin(θ), z = z. (3.21)
Debido a este cambio de coordenadas el campo de deformaciones característico del sistema
en terminos de las nuevas variablesr y θ no es el mismo que el campo de deformaciones para
x y y. Sin embargo, por medio de una rotación se puede obtener una relación este las de-
formaciones normales y cortanteer, eθ y erθ y las correspondientes componentes cartesianas
10
así
er =ex + ey
2+ex − ey
2cos(2θ) + exy sin(2θ),
eθ =ex + ey
2− ex − ey
2cos(2θ)− exy sin(2θ),
erθ =ey − ex
2sin(2θ) + exy cos(2θ). (3.22)
De igual manera se puede obtener una relacion para los esfuerzos así:
σr =σx + σy
2+σx − σy
2cos(2θ) + τxy sin(2θ),
σθ =σx + σy
2− σx − σy
2cos(2θ)− τxy sin(2θ),
τrθ =σy − σx
2sin(2θ) + τxy cos(2θ). (3.23)
Como el problema se resolvió para las deformaciones en coordenadas polares y las deforma-
ciones experimentales o simulaciones van a estar en coordedanas cartesianas, es necesario
encontrar una matriz de rotacion para trabajar en un solo sistema de coordenadas. De las
ecuaciones (3.22) se puede obtener la siguiente matriz de rotación:
ex
ey
exy
=
(1 + cos(2θ))
2
(1− cos(2θ))
2− sin(2θ)
(1− cos(2θ))
2
(1 + cos(2θ))
2sin(2θ)
sin(2θ)
2−sin(2θ)
2cos(2θ)
er
eθ
erθ
. (3.24)
3.6. Solución del Problema Planteado
Consideremos una placa infinita con la presencia de un agujero circular y sometida a un
campo de esfuerzo uniforme, tal y como se ve en la figura 3.2. Dada esta situación, el interés
es poder obtener el campo de deformaciones como función del radio del agujeroa y de la
11
posición enr y θ sobre la placa, para así poder realizar una comparación entre el campo de
deformaciones teórico y experimental.
Fig. 3.2: Placa infinita con agujero sometida a esfuerzos uniformes.
Inicialmente consideremos el estado de esfuerzos suponiendo que no existe agujero, para
este caso se tiene que
σx = T σy = τxy = 0. (3.25)
Ahora usando función de esfuerzos de Airy se tiene que:
σx =∂2φ
∂y2= T ⇒ Integrando respecto ay dos veces (3.26)
φ =1
2Ty2 =
1
2Tr2 sin2(2θ) =
1
2Tr2(1− cos(2θ)). (3.27)
Si se considera la presencia del agujero, está afecta la uniformidad del campo y se espera que
cuando nos ubiquemos a una distancia lejos del agujero, dicha perturbación caiga a cero y el
campo se mantenga uniforme. Teniendo en cuenta lo anterior, se supone una solución para la
función de Airy (Michell) que incluye los términos de asimetría ycos(2θ) así
φ = ao + a1 log r + a2r2 + a3r
2 log r
+ (a21r2 + a22r
4 + a23r−2 + a24) cos(2θ). (3.28)
12
Ahora usando las ecuaciones (3.27) y (3.28) se tiene que:
σr = a3(1 + 2 log r) + 2a2 +a1
r2− (2a21 +
6a23
r4+
4a24
r2) cos(2θ),
σθ = a3(3 + 2 log r) + 2a2 −a1
r2+ (2a21 + 12a22r
4 +6a23
r4) cos(2θ),
τrθ = (2a21 + 6a22r2 − 6a23
r4− 2a24
r2) sin(2θ). (3.29)
Para obtener las condiciones iníciales en términos der y θ basta con sustituir las condiciones
3.25 en las ecuaciones 3.23 para tener:
σr(a, θ) = τrθ(a, θ) = 0,
σr(∞, θ) =T
2(1 + cos(2θ),
σθ(∞, θ) =T
2(1− cos(2θ),
τrθ(∞, θ) = −T2
sin(2θ). (3.30)
Analizando las ecuaciones 3.29 y teniendo en cuenta que el esfuezo debe ser finito cuando
r 7→ ∞ se tiene quea3 = a22 = 0. Por último, utilizando las condiciones de frontera (3.30)
se tienen cinco ecuaciones y cinco incógnitas, de donde se obtiene que:
2a1 +a2
a2= 0,
2a21 +6a23
a4+
4a24
a2= 0,
2a21 −6a23
a4− 2a24
a2= 0,
2a21 +T
2= 0,
2a2 −T
2= 0. (3.31)
Solucionando para las constantesaij se tiene que
a1 = −a2T
2, a2 =
T
4, a21 = −T
4, a23 = −a
4T
4, y a24 =
a2T
2. (3.32)
13
Finalmente, sustituyendo los valores de las constantes en las ecuaciones 3.29 se tiene que:
σr =T
2
(1− a2
r2
)+T
2
(1 +
3a4
r4− 4a2
r2
)cos(2θ), (3.33)
σrθ =T
2
(1 +
a2
r2
)− T
2
(1 +
3a4
r4
)cos(2θ), (3.34)
τrθ = −T2
(1− 3a4
r4+
2a2
r2
)sin(2θ). (3.35)
Las ecuaciones 3.35 representan el estado de esfuerzo para el problema planteado, como de
las mediciones experimentales o simulaciones se van a obtener los estados de deformación
debemos usar la ley de Hooke generalizada (3.4) para pasar del estado teórico de esfuerzos
al estado de deformaciones. A continuación se presentan dichas relaciones en coordenadas
polares
er =1
E(σr − νσθ), (3.36)
eθ =1
E(σθ − νσr), (3.37)
erθ =1 + ν
Eτrθ. (3.38)
14
Capítulo 4
Elementos Finitos - Simulaciones
El análisis de estructuras complejas y otros sistemas en una formulación matricial es ac-
tualmente inimaginable sin el método de elementos finitos. Se cree que el origen de esta
aplicación tan rica y valiosa, no puede ser atribuida específicamente a tan solo una persona o
grupo investigativo sino mas bien a un conglomerado de descubrimientos científicos realiza-
do por un vasto grupo de personas.
La noción de divisiones geométricas se marca tiempo atrás con el filósofo Archimides,
quien para calcular el área de una superficie compleja, la dividió en pequeños triángulos y
cuadriláteros cuyas áreas podrían ser fácilmente calculadas; y elensamblajede estas áreas
individuales daban como resultado el área deseada.
Inicialmente este método era conocido esencialmente por la mecánica estructural y la in-
dustria aeroespacial. Actualmente, éste ha logrado tan alta acogida y maduréz que es pensado
como el método para solucionar problemas generales espaciales en todas las áreas de la inge-
niería y la ciencia de la física. Una definición moderna del Método de Elementos finitos dice
15
que es simplemente una técnica numérica para obtener soluciones aproximadas a un grupo
de ecuaciones diferenciales parciales.
4.1. El procedimiento de Elementos Finitos
Este proceso involucra tres fases generales: Una fase de pre-procesamiento, un fase de
análisis, y una fase de post-procesamiento. Posiblemente, la fase más compleja y la que más
demanda tiempo es la fase del pre-procesamiento. Dado un problema físico de interés, el
matemático junto con conceptos de ingeniería debe usar un juicio para transformar el prob-
lema del mundo real a un modelo matemático (discreto) que contenga toda la física esencial
y el cual sea posible solucionar numéricamente. Después de hacer suposiciones razonables,
dando así con un modelo matemático, donde normalmente se involucran ecuaciones difer-
enciales junto con condiciones de frontera o iníciales, se debe pasar a decidir qué tipo de
análisis se va a realizar: lineal, no-lineal, transciente, estado estable, etc.
Posteriormente se debe escoger el tipo de material para poder establecer las propiedades
mecánicas mencionadas en secciones anteriores así como establecer el tipo de elemento fini-
to, y sobre todo delimitar una malla refinada en regiones donde se esperan largos gradientes
del campo a tratar. Afortunadamente, gracias a la tecnología la fase de pre-proceso se ha
convertido en una fase más amigable por medio de software especializado.
4.2. Simulaciones
4.2.1. Elementos a utilizar y suposiciones
El programa que se va a utilizar para el desarrollo de las simulaciones esANSYS 10.0.
Puesto que se va a modelar una placa en 2D con un pequeño espesor de1,2 mm (Aluminio
16
obtenido en el mercado) se decidió utilizar una simulación de esfuerzo en el plano con es-
pesor. El tipo de elemento seleccionado para realizar el mallado fue PLANE2, el cual es un
triángulo con 6 nodos tal y como lo muestra la figura 4.1. El elemento tiene comportamiento
de desplazamiento cuadrático y es muy bueno para simulaciones con mallas irregulares por
lo cual es adecuado para el tipo de problema a analizar.
Fig. 4.1: Elemento PLANE2 de Ansys utilizado en las simulaciones.
Por otro lado la simulación se escogió como estructural, lineal, elástica e isotrópica con
los valores deE y ν del material obtenido. Para poder obtener estas propiedades, se realizó
un prueba de tensión estándar en la máquina INSTRON 5586 bajo la normaASTM B577-
02a. El ensayo se realizó para una muestra de 7 probetas, a una velocidad de5 mm/min, con
una temperatura de23oC y con una humedad del50 %. Los resultados de estos ensayos se
presenta en la figura 4.2.
Con esta información se realizaron tres aproximaciones del modulo de Young: La primera
tomando los dos primeros datos experimentales y hallando la pendiente; la segunda, tomando
los 3 primeros datos y tomar su regresión lineal y la última haciendo una regresión con los
3 primeros datos y obligándola a pasar por el cero. De estos datos se realizó un análisis de
varianza (ver Fig 4.3) para verificar la validéz de los datos. Del análisis estadístico se definió
un módulo de Young de66378 MPa
Inicialmente, se realizó una simulación con una plata completa de200 mm por200 mm
17
Fig. 4.2: Diagrama Esfuerzo vs Deformación unitaria. Ensayo de tensión Universal
Fig. 4.3: Análisis estadístico. ANOVA
y con un agujero de15 mm de radio para verificar la simetría tal y como lo muestra la figura
4.4. Debido a esto, se decidió simular solamente el primer cuadrante puesto que se obtiene
la misma información que si se simulara toda la placa pero a un menor costo de tiempo y
memoria computacional.
Para las simulaciones definitivas se tomó un cuarto de placa de100 mm por100 mm con
radios de agujero variando entre5 y 25 mm con intervalos de5 mm cada uno. Para poder
obtener una simulación adecuada, el tamaño de la malla no debe ser muy grande y se debe
18
Fig. 4.4: Simulación placa completa. Estado de deformación enx (arriba), Estado de deformación eny (abajo)
refinar en los sitios de interés o en los cambios bruscos de geometría, en este caso cerca al
agujero. Un mallado final se puede observar en la figura 4.5.
Fig. 4.5: Mallado final.
Lo último que se debe definir para la simulación antes de pasar a solucionar, son las condi-
ciones de frontera. Se desea obtener las deformaciones causadas por un esfuerzo uniforme,
para esto se definieron dos restricciones de movimiento: Los elementos de la línea vertical
19
izquierda estaban restringidos para moverse en la dirección horizontal y los elementos de
la línea horizontal inferior estaban restringidos para moverse verticalmente. Así mismo, se
definió un esfuerzo axial constante en la dirección horizontal con una magnitud de75 MPa.
4.3. Resultados de las simulaciones
Los datos que se querían obtener de las simulaciones, eran el estado de deformaciones
para una ubicación determinada sobre la placa; es decir, una ventana de deformaciones. Para
lograr esto, es necesario crear unos caminos (path en el lenguaje de programación de AN-
SYS), los cuales permiten obtener para una posición enx y en y dada, la información del
estado de esfuerzos y deformaciones, entre muchas otras variables de salida. Se definieron 10
caminos igualmente espaciados en la dirección vertical, donde para cada camino se definieron
las mismas 10 posiciones enx igualmente espaciadas, para lograr así una ventana cuadrada
de10 por 10 datos. El espesor de la ventana como su espaciamiento pueden ser modificado
simplemente. Para efectos de este documento se utilizó un tamaño de ventana cuadrada de
10 mm. El código en APDL que se utilizó para obtener estos datos en un formato de texto se
presenta en el apéndice C, al final de este documento. A continuación, se presentan imágenes
de los estados de deformación obtenidos en las simulaciones.
Fig. 4.6: Estado de deformaciones en la dirección x (izq) y en la dirección y (der).
20
De estos estados, se puede ver que se tiene una deformación máxima enx en la parte su-
perior del agujero, mirando los valores de los esfuerzos en este punto se observo un esfuerzo
máximo de75,45 MPa, lo cual indica que el agujero actúa como un concentrador de esfuer-
zos de3 veces el esfuerzo aplicado. Por el lado de las deformaciones eny se observa que se
generan más contraciones que elongaciones, teniendo un máximo ubicado aparentemente a
45o del eje horizontal.
21
Capítulo 5
Algoritmos Genéticos
5.1. ¿Qué son los Algoritmos Genéticos?
Los Algoritmos Genéticos son rutinas numéricas de optimización inspiradas en la selec-
ción natural y la genética, procesos que se derivan de la evolución biológica. Fueron desar-
rollados por primera vez en la década de los sesentas por John Holland. El método es muy
general y puede ser aplicado en un amplio rango de problemas. Este algoritmo es muy sencil-
lo de entender y el código computacional es fácil de implementar. En pocas palabras, lo que
hace el algoritmo genético es modificar repetidamente una población de posibles soluciones,
en cada paso el éste selecciona unos padres de forma aleatoria y los usa para generar una nue-
va población. Después de varias generaciones, la poblaciónevolucionaa una solución óptima.
Por lo general un Algoritmo genético consta de los siguientes componentes:
Una población inicial que se encuentre dentro del conjunto de posibles soluciones al
problema.
Una forma de calcular qué tan bueno o malo es el desempeño de cada individuo dentro
22
de la población.
Un método para cruzar las mejores soluciones (padres) y así conformar nuevos indi-
viduos que podrían llegar a tener mejor desempeño.
Un operador de mutación para aumentar la diversidad en los individuos y así mejorar
la búsqueda.
El algoritmo genético usa tres operaciones principales en cada iteración para la construcción
de la siguiente generación:
1. Selección:Permite escoger individuos llamados padres, los cuales van a contribuir sus
propiedades para la construcción de la siguiente generación.
2. Cruce:Combina dos padres para formar un hijo de la siguiente generación.
3. Mutación: Aplica cambios aleatorios a los padres para formar una diversidad de hijos.
Más adelante en el documento se van a explicar el tipo de operadores que se utilizarán.
5.2. Implementación del Algoritmo Genético para la solu-
ción del Problema Planteado
5.2.1. Momentos Geométricos
Antes de implementar el algoritmo, es necesario definir las variables que se van a utilizar
dentro de algoritmo en términos de las variables de nuestro problema. Para lograr esto, se
realizarán una serie de definiciones a continuación. Primero se debe definir la función a op-
timizar, la cual en este caso sería el error entre los datos experimentales o simulados y los
teóricos. Sin enmbargo, como se tiene toda unaventanade datos, lo recomendado es tratar
23
de reducir la cantidad de datos contenidos en la ventanta y al mismo tiempo mantener su
significado físico, esto se logra cálculando los momentos geométricos de la función de defor-
macionex(x, y). Físicamente, los momentos son una proyección de la función de densidad
ex(x, y) sobre un conjunto de funciones base. El propósito de esta proyección es capturar
la información global sobre la ventana de deformaciones y reemplezar el espacio contínuo
de ex(x, y) por un espacio discreto de momentos [4]. Los momentos geométricos (GM) se
definen matemáticamente en términos de la integral de Riemann de la siguiente manera:
mpq =
∫ ∞
−∞
∫ ∞
−∞xpyqex(x, y)dxdy (5.1)
Dondep+ q es el orden del momento. La integral de la ecuación (5.1) puede ser aproximada
discretamente por medio de la siguiente doble sumatoria
mpq =N∑
i=1
M∑j=1
xpyqex(x, y)∆x∆y (5.2)
donde∆x y ∆y son los intervalos de los caminos definidos para la obtención de la ventana
de deformaciones yMxN es el área de dicha ventanta. En muchas ocaciones el valor de los
momentos puede ser muy pequeño dependiendo del valor de la función densidad como por
ejemplo, en este caso las deformaciones unitarias son del orden de10−3. Para solucionar este
problema se realiza la siguiente ponderación
mpq =mpq
m00
(5.3)
Con este concepto de momentos geométricos introducido, se puede pasar a definir la función
objetivo del problema inverso de optimizacion de la siguiente manera
Función objetivo (Error):
mınx∈(a,x,y)
ψ(x) = mınx∈(a,x,y)
√ ∑p+q≤1
(mpqexp(x)− mpqteo(x)
)2
Sujeto a: a2 ≤ x2 + y2 0 ≤ x ≤ xmax 0 ≤ y ≤ ymax (5.4)
24
5.2.2. Componentes del Algoritmo Genético
Función de Desempeño:
La función de desempeño es aquella que se busca optimizar en el algoritmo. Para algoritmos
estándar de optimización, ésta se conoce como función objetivo. Sin embargo, en nuestro
problema lo que se quiere es maximizar úna función de desempeño que implique una mini-
mización de la función objetivo (5.4). Al definir la funcion desempeño de la siguiente manera,
se logra este objetivo
F (x) =1
1 + ψ(x)(5.5)
de esta definición es claro queF (x) logra su máximo cuandoψ(x) logra su minimo, y el
valor del máximo en este caso sería uno.
Individuos: Un individuo es cualquier punto al cual se le puede aplicar la función de desem-
peño. El valor de la función para dicho individuo se denominarádesempeño. En esta caso el
individuo es una 3-tupla(a, r, θ) tal quex pertenezca al rango de la ventana de medicion. Al
individuo x se le dicecromosomay a los componentes de la 3-tupla,a, x y y se les llama
comunmentegenes.
Población: Una población es un arreglo de individuos. Por ejemplo, si el tamaño de la
población es50 y el numero de cromosomas en la función de desempeño es3 , la población
va a estar representada por una matriz de50x3.
Operadores:
Selección:Existen muchas formas de hacer la selección, entre ellas están: asignación
en fila, selección por ruleta, muestreo estocástico universal, selección local, selección
truncada, selección por torneo, entre otros. En este caso se decidió trabajar con la se-
lección por ruleta, por su eficiencia y simplicidad.
25
La selección por ruleta funciona de la siguiente manera: A cada individuo se le asigna
un valor de la función desempeño, un valor entre cero y uno. A continuación, se orga-
nizan los valores de desempeño de mayor a menor. Luego se suman todos estos valores
y a este valor le llamaremossuma. Se genera aleatoriamente un valor entre cero y suma
el cual llamaremoss. Por último, se empieza un proceso iterativo empezando con el
mayor valor de desempeño, si este valor es mayor ques se escoge el individuo al que
corresponde dicho valor de desempeño, de lo contrario se suma el valor asignado al
siguiente individuo y se realiza la misma comparación, hasta que un individuo sea se-
leccionado. Para seleccionar el siguiente individuo, se realiza el mismo proceso pero
sin incluir al previamente seleccionado.
Este método es llamado ruleta porque la probabilidad de selección es igual a la de una
ruleta en la que las rebanadas son proporcionales al desempeño del individuo, siendo
mas probable de ser escogido si se tiene una rebanada más grande.
Cruce: Se van a utiliar dos tipos de cruce: Intermedio y Lineal
1. Intermedio: Supongamos que tenemosp1 y p2 como los padres, entonces los
nuevos hijosφ1 y φ2 serán:
φid = pi
1 + vid
(pi
2 − pi1
)(5.6)
Coni =1, 2 y 3 puesto que cada padre consta de 3 componentes,d = 1, 2 y donde
vid son valores aleatorios. De la ecuación (5.6), se puede ver que los hijos serán
escogidos en alguna vecindad de los padres.
2. Lineal: Similar al cruce intermedio, los hijosφ1 y φ2 se escogerán alrededor de
los padres de la siguiente manera
φd = p2 + vd (p1 − p2) (5.7)
26
Donded = 1, 2 y vd valores aleatorios.
Mutación: Esta función cambiará el gen de uncromosomapi, por un nuevo valor. Se
escoge al azar un gen del cromosoma a mutar y al valor de este gen se le suma o resta
un porcentaje de su valor actual. La escogencia del gen tiene una probabilidad uniforme
sobre todos los genes del cromosoma. El valor porcentual a sumar o restar varía entre
2 y un 15 %.
Criterios de detención:Existen varios criterios de dentención
Número de generaciones:El algoritmo para cuando se llegue a un número de genera-
ciones dado.
Tiempo límite: El algoritmo para cuando se sobrepase el tiempo límite de ejecución.
Límite de desempeño:El algoritmo para cuando el valor de la función de desempeño
para el mejor individuo de la población sobre pase un valor límite dado.
Estancamiento:El algorimto para cuando la función de desempeño no mejore despues
de un número dado para el estancamiento.
Para la solución del problema se utilizarán como métodos de detención los criterios de
Número de generaciones y Límite de desempeño.
5.2.3. Etapas del algoritmo genético
Población Inicial: El algoritmo empieza creando una población aleatoria de invidivi-
dos que estan dentro del rango de la ventana. Esto se va a lograr por medio de una
distribución uniforme en los paremetrosx, y y paraa se distribuye uniformemente
entre0 y√x2 + y2.
27
Siguiente Generación:El algoritmo selecciona individuos de la población actual, es
decir los padres, por medio del método ruleta, para crear sus respectivos hijos. A la
siguiente generación siempre van a pasar los primero dos padres, y luego por medio
de itereaciones entre los operacionescruce lineal, cruce intermedio,y mutaciónse
generan hijos hasta completar el tamaño de la poblacion inicial.
De lo anterior se puede ver que el tamaño de cada poblacion no cambia, sin embargo
algunos algoritmos tienen el tamaño de la población como un parametro más, comen-
zando desde un tamaño más grande, para luego ir disminuyendo su tamaño a medida
que la funcion de desempeño tiende a cero. Esto se hace para mejorar los tiempos de
convergencia.
La figura 5.1 muestra la ubicación en las coordenasx y y para una poblacion inicial,
sus respectivos padres y una nueva generación.
Fig. 5.1: Resultados de población inicial, padres y nueva generación en una etapa primaria.
28
Criterio de detención: Finalmente se mira si se cumplen los requisitos de detención y
de no ser asi, se genera una nueva población a analizar.
A continuación, se presenta un esquema que resume las etapas del algoritmo genético previ-
amente explicadas (ver Fig 5.2) :
Fig. 5.2: Etapas del Algoritmo Genético.
29
Capítulo 6
Resultados
Como se mencionó anteriormente este trabajo pretende caracterizar agujeros en placas
planas. Se realizaron simulaciones para radios de agujeros que varían desde 5 hasta 25 mm.
De estos resultados se obtuvieron lasdeformaciones experimentales, asi mismo, por medio de
la teórica clásica de elasticidad se obtuvieron los valores de las deformaciones como función
de los parámetrosa, x y y. Finalmente, se desea optimizar el error de dichas deformaciones
para los parámetros descritos por medio de algoritmos genéticos y así lograr caracterizar el
agujero.
Toda la implementación del código de optimización por medio de algoritmos genéticos
con las funciones descritas en los capítulos anteriores, se realizo en el programa MATLAB
7.0.1.(Apéndice D) Esta es una herramienta computacional para un manejo práctico y efi-
ciente de matrices, por ende es un entorno adecuado para solucionar el problema planteado
ya que las variables de interés se pueden representar matricialmente.
30
6.1. Resultados del Algoritmo Genético
Lo primero que se quiere verificar antes de ver que tan buenos son los resultados, es ver
si la función objetivo es minimizada. Debido a la definición de la función de desempeño de
nuestro algoritmo genético, si ésta logra su valor máximo es decir igual a uno, implicaría que
la función objetivo o de error si se está minimizando. En la figura 6.1 se presenta una imagen
del comportamiento típico de una función de desempeño convergente.
Fig. 6.1: Comportamiento convergente de la función desempeño.
Esta figura muestra que efectivamente se está logrando optimizar la función objetivo de-
seada. Es por esto que las gráficas obtenidas por el modelo de algoritmo genético se deben
comportar así. Como el interés es observar el comportamiento de las variablesa, x y y, la
salida del algoritmo tanto gráfica como numérica, son estas variables cuando la entrada es el
campo de deformación junto con las propiedades mecánicas del material simulado. La figura
6.2 muestra las gráficas finales generadas por el algoritmo genético después de realizar 10
simulaciones simultaneas.
31
Fig. 6.2: Comportamiento convergente de la función desempeño. Las lineas rojas represetan los valores reales utilizados en las simulaciones
para obtener las deformaciones experimentales. Para ver otros resultados remitirse al apéndice B
De la figura 6.2 se puede concluir que sí se está logrando optimizar el error planteado
dentro del modelo y que la función de desempeño si converge a uno en todos los escenarios.
Además, las posiciónes enx y en y se logran estimar correctamente ya que convergen a la
línea roja como es lo deseado. Sin embargo, se puede ver que el valor del tamaño del agujero
es inestable y no converge a su estado deseable. Antes de concluir sobre los datos se considera
necesario presentar los valores numéricos obtenidos con sus respectivos errores.
32
Fig. 6.3: Promedios de todos los parametros despues de 15 simulaciones.
Fig. 6.4: Desviacion estándar de todos los parametros despúes de 15 simulaciones.
Debido a que las escalas de los parámetros no son iguales, se decidió calcular el coefi-
ciente de variación para que todos los datos sean comparables dentro de escalas diferentes ya
que este parámetro es independiente de la escala en la que se está trabajando si se asegura
que todos los datos a utilizar sean positivos.
Los datos presentados en las figuras (6.3), (6.4) y (6.5) muestran la excelente conver-
gencia de los parámetrosx y y. Sin embargo se ve un error significativo en los resultados
referentes al tamaño del agujero. Para tratar de solucionar este problema se realizó una mod-
ificación a la función objetivo, manteniendo sus propiedades y criterios de convergencia. Lo
que se hizo fue replantear su definición de la siguiente manera
33
Fig. 6.5: Coeficiente de variacion de todos los parametros despúes de 15 simulaciones.
F (x) =k
k + φ(x)(6.1)
Se puede ver de la definición (6.1) que la función de desempeño mantiene sus propiedades.
Si la funciónF 7−→ 1 entoncesφ 7−→ 0 lo cual significa que se esta minimizando la función
objetivo. Los resultados del algoritmo para diferentes valores del parámetrok se presentan
en la figura 6.6
Fig. 6.6: Resultados del algoritmo para 15 simulaciones con variaciones sobre el parámetrok .
34
Como se puede ver en la figura anterior (ver Fig. 6.6) el valor dek que mejora la conver-
gencia del parámetroa esk = 0,6. Con este valor, se paso a evaluar el modelo para valores
de k cercanos a0,6. Los mejores resultados fueron para valores dek = 0,55 y k = 0,65.
Las figuras 6.7 y 6.8 muestran los valores obtenidos después de 20 corridas del algoritmo
genético.
Fig. 6.7: Resultados del algoritmo para 20 simulaciones. Parametrok = 0,55
Fig. 6.8: Resultados del algoritmo para 20 simulaciones. Parametrok = 0,65
Finalmente se logró mejorar la convergencia del parámetroa tal y como lo muestran las
figuras 6.7 y 6.8. Sin embargo, vale la pena observar que se tiene una desviación estándar
bastante grande y por consiguiente se debe tener cuidado a la hora de usar los valores del
tamaño de agujero dado por el modelo al momento de analizar estructuras.1
1Nota: En los anexos se encuentran más resultados de las simulaciones y los códigos en archivos.m del
algoritmo genético.
35
Capítulo 7
Conclusiones
Se logró una técnica eficiente, económica y sencilla de caracterización de agujeros
en placas planas. Esta técnica podría ser utilizada como método de mantenimiento
preventivo en la industria Colombiana, disminuyendo los costos y logrando la exactitud
que dicho proceso requiere.
Se pudo ver de los resultados, que independiente de la ubicación de la ventana y del
tamaño del agujero se logro estimar correctamente la ubicación exacta de la falla, tanto
su posición enx, como su posición eny, permitiendo que el modelo sea aplicable a
distintas geometrías.
En cuanto al tamaño del agujero, este parámetro no ha convergido perfectamente, se
debe modificar al algoritmo para que se estabilice en los valores deseados. Se cree
que una posible forma de lograr esto, es hacer que las restricciones del parametroa se
modifiquen en cada nueva generación, tratando de seguir una distribución normal, con
paramétros el promedio y la desviación de la poblacion actual.
La teórica clásica de elasticidad modela correctamente el comportamiento de una placa
36
sometida a un esfuerzo uniforme. A pesar de que las ecuaciones que rigen dicho modelo
son un poco abstractas, es posible entender su derivación.
El método de los momentos geométricos para reducir el tamaño de la información sin
perder el significado presente en los datos fue satisfactorio. Además, si se tiene infor-
mación de las deformaciones de una placa por medio de métodos de interferometría
seria indispensable utilizar esta herramienta pues de lo contrario los cálculos serían
exhaustivos, inclusive para un computador.
Se comprobó que los métodos de optimización por medio de algoritmos genéticos, son
fáciles de implementar y permiten una conexión entre la ingeniería y las ciencias de la
naturaleza para la solución de problemas aplicados, los cuales serían muy difíciles (por
no decir casi imposibles) de solucionar por otros métodos más sofisticados.
37
Apéndice A
8.1. Introducción y descripción de propiedades mecánicas
Esfuerzo: Consideremos una sección de área subdividida en áreas pequeñas llamadas
∆A. Si queremos reducir cada vez más el tamaño de∆A, debemos tener en cuenta las sigu-
ientes consideraciones. Se trata de un material continuo, en donde la materia esta distribui-
da uniformemente en vez de estar compuesto por un número finito de átomos o moléculas.
Además, todas las partículas están conectadas entre si, es decir se trata de un material co-
hesivo, en donde no se presentan grietas, huecos o separaciones. Ahora consideremos que
una fuerza esta actuando sobre esta área, como es sabido esta fuerza estará actuando en una
única dirección, pero por razones prácticas se hablará de sus tres componentes,∆Fx ,∆Fy y
∆Fz. Mientras∆A se acerca a cero, también lo hace∆F y sus componentes; sin embargo, el
cociente de dicha fuerza en el área se acercará a un número finito. A este cociente se le llama
ESFUERZO, y describe la intensidad de la fuerza interna en un área especifica.
Esfuerzo normal: La intensidad de la fuerza por unidad de área, actuando normal a∆A,
se define como el esfuerzo normalσ
σz = lım∆A→0
∆Fz
∆A. (8.1)
Esfuerzo cortante:La intensidad de fuerza por unidad de área, actuando tangente a∆A,
38
se define como el esfuerzo cortante
τzx = lım∆A→0
∆Fx
∆Aτzy = lım
∆A→0
∆Fy
∆A. (8.2)
Estado general de esfuerzos:Si analizamos planos en las direcciones x-z y y-z se puede
obtener una pequeño elemento cúbico del material el cúal representa el estado de esfuerzo
actuando alrededor de un punto específico del cuerpo.
Fig. 8.1: Diagrama de esfuerzos.
Elongación:La elongación o contracción de un segmento de línea por unidad de longitud
se define como elongaciónε donde:
ε =lo − lflf
= ∆S. (8.3)
Ley de Hooke: Si se hace un diagrama de esfuerzo - deformación unitaria se puede
observar que los materiales de ingeniería muestran una zona en la que se tiene una relación
lineal (zona a la cual se le llama región elástica), como consecuencia un aumento en esfuerzo
genera un aumento proporcional en deformación. Este fenómeno fue descubierto por Robert
39
Hooke en 1676. Matemáticamente esta ley se puede escribir de la siguiente forma:
σ = Eε. (8.4)
En dondeE representa la constante de proporcionalidad de cada material, la cual es llamada
el módulo de elasticidad o el módulo de Young.
Módulo de Poisson:Cuando un cuerpo es sometido a fuerzas de tensión, el cuerpo no
solo se deforma expandiéndose sino que también se contrae lateralmente. Análogamente,
cuando un cuerpo se contrae sus paredes laterales tienden a expandirse. Si suponemos que la
elongación esδ y que la deformación radial fueδ′ tenemos que las elongaciones longitudinal
y radial son
εradial =δ′
rεlongitudinal =
δ
L. (8.5)
En los primeros años del sigloXIX, el científico francés S.D. Poisson descubrió que dentro
del rango elástico (en donde se cumple la ley de hooke) el cociente entre los valores de
deformación se mantiene constante, puesto que las deformaciones son proporcionales. A esta
constante de proporcionalidad se le llama módulo de Poisson y es un valor numérico único
para cada material que es homogéneo e isotrópico. Matemáticamente el módulo de Poisson
se define asi:
ν = − εradial
εlongitudinal
. (8.6)
El signo negativo es debido a que uno de las dos deformaciones tiene un valor negativo, sea
por una elongación longitudinal (positivo) la cual causa una contracción lateral (negativo)
o vice versa. La constante de Poisson es adimensional y para la mayoría de materiales no
porosos su valor está entre14
y 13.
40
Módulo de rigidez: Cuando se tienen esfuerzos cortantes se miden los ángulos de defor-
mación de las caras del cuadrado y para estas condiciones la ley de Hooke se puede escribir
de la siguiente manera:
τ = Gγ. (8.7)
DondeG es la constante de rigidez de un material homogéneo e isotrópico.
8.2. El análisis de esfuerzos
Para especificar el estado de esfuerzos de un puntoP en el continuo, consideremos un
elemento∆S de una superficie que pasa a través deP . Sean un vector unitario normal
a dicha superficie en el puntoP , llendo del lado 1 al lado 2. Las fuerzas ejercidas por el
material en el lado 2 de∆S sobre el material del lado 1 se pueden representar por una fuerza
∆T enP y de un momento∆G
lım∆S→0
∆T
∆S= Tn, lım
∆S→0
∆G
∆S= 0, (8.8)
dondeT n es un vector llamadotraccióno vector de esfuerzosa través de la superficie en el
puntoP . En generalT n no tendrá la misma dirección quen, pero se puede descomponer en
cada una de las tres direcciones de coordenadas, así:
Tn = σnxi + σnyj + σnzk. (8.9)
Si tenemos por ejemplon = x, entonces la tracción en la superficie con vector normali es:
Ti = σxxi + σxyj + σxzk. (8.10)
Los valoresσxx, σyy y σzz son llamados esfuerzosnormales; y los esfuerzosσxy, σxz, σyz son
llamados esfuerzoscortantes.
41
Se ha encontrado que las fuerzas actuando sobre un volumen∆V se pueden dividir entre
fuerzas de superficie causadas por las fuerzas de tracción o esfuerzos descritos anteriormente
y por fuerzas de cuerpo rígido. El ejemplo más claro de lo anterior es la gravedad, la cual
ejerce una fuerza hacia abajo con magnitud(ρg∆V ) en un elemento de volumen∆V y den-
sidad de masaρ. Si se asume que las fuerzas actuando en un pequeño elemento de volumen
alrededor de un puntoP son equivalentes a∆R enP y a un momento∆M , entonces
lım∆V→0
∆R
∆V= F, lım
∆V→0
∆M
∆V= 0. (8.11)
El vectorF es llamado fuerza por unidad de volumen. Si se asume queF es una funcion
continua de las coordenadas(x, y, z) se puede escribir
F = Fxi + Fyj + Fzz. (8.12)
Se puede demostrar facilmente [3] que los esfuerzos cortantes son simétricos es decir:σij =
σji parai 6= j.
Ahora consideremos un pequeño elemento rectangular de un tetrahedroPABC con tres
lados parelelos a los planos coordenados y un cuarto lado con vector normaln = (l,m, n),
tal como se ve en (Fig. 8.2).
Consideremos las fuerzas actuando sobre el elemento en la direcciónx. Si el áreaABC es
∆S, y la alturaPN esh, entonces las áreasPBC,PCA y PAB seránl∆S,m∆S, n∆S, y el
volumen del tetrahedro será∆V = h∆S/3. Las condiciones de equilibrio sobre el elemeneto
en la direcciónx generan la siguinete ecuación:
∆Rx + ∆S(σnx − lσxx −mσyx − nσzx) = 0. (8.13)
Los esfuerzos de superficie, cuandoh y ∆S son suficientemente pequeños, se pueden
42
Fig. 8.2: Fuerzas actuando en un elemento tetrahedríco.
aproximar por los valores en el puntoP . Dividiendo la ecuación (8.13) por∆S y consideran-
do límites tenemos:
lım∆V→0
∆Rx
∆S= lım
∆V→0
∆Rx
∆V· ∆V
∆S= lım
∆V→0
∆Rx
∆V· h3
= 0. (8.14)
De donde obtenemos la primera de las tres ecuaciones
43
σnx = lσxx +mσyx + nσzx,
σny = lσxy +mσyy + nσzy,
σnz = lσxz +mσyz + nσzz.
(8.15)
Muchas veces es conveniente usar sufijos numéricos1, 2, 3, en vez dex, y, z. En está no-
tación, las componenesτi de los esfuerzos o tracciones a través de una superficion con vector
normalν1, ν2, ν3 dados por (8.15) es:
τi =3∑
j=1
σijνj. (8.16)
También es común el uso de la convención de doble sufijo de Einstein y reescribimos lo
anterior así:
τi = σijνj. (8.17)
Esta ecuación nos permite mostrar que las variablesσij son las componentes de untensorde
orden 2. Si cambiamos de eje de coordenadasi1, i2, i3 por otros ejesi′1, i′2, i
′3, tenemos con la
notación de doble sufijo,
ii = liji′i, (8.18)
entonces los esfuerzosσ′ij en el nuevo eje de coordenadas son
σij = likljlσ′kl. (8.19)
8.3. El análisis de deformación
Consideremos un rectangulo en el planoxy con ladosdx y dy (8.3). Bajo la aplicación
de una carga axial el elemento se extiende en la direcionx una cantidaddux. La deformación
44
por unidad de área se define como:
εxx =ux(x+ dx, y, , z)− ux(x, y, z)
dx=∂ux
∂x. (8.20)
De la misma forma una carga puede ser aplicada en la direccióny causando una deformación
εyy y una carga puede ser aplicada en la direcciónz causando una deformaciónεzz.
εxx =∂ux
∂x, εyy =
∂uy
∂y, εzz =
∂uz
∂z. (8.21)
Fig. 8.3: Deformaciones axiales y cortantes.
Ahora supongamos que el lado vertical del rectangulo se desplaza en la direcciónx por
dux y el lado horizontal se mueve en el planoy porduy. Como resultado, el lado vertifal rota
en la dirección contrareloj por un anguloθ2 mientras el desplazamiento verticalduy induce
una rotacion dirección relojθ1. Debido a estas rotaciones el angulo originalπ/3 se convierte
enφ y el cambio de ángulo es:
45
γ =π
2− φ = θ1 + θ2. (8.22)
En el contexto de la teoria de pequeña deformaciones tenemos:
tan θ1 =∂uy
∂x≈ θ1, tan θ2 =
∂ux
∂y≈ θ2, (8.23)
Ahora podemos definir como deformaciones cortantes las siguientes cantidades:
εxy =
(∂ux
∂y+∂uy
∂x
),
εxz =
(∂ux
∂z+∂uz
∂x
),
εyz =
(∂uy
∂z+∂uz
∂y
)(8.24)
Las dos deformaciones cortantesεxx, εyy son obtenidas considerando desplazamientos
angulares en los planosxz y yz, respectivamente. Luego tenemos que las cantidades
ε11 = εxx, ε12 = εxy, ε13 = εxz, ..., , ε33 = εzz, (8.25)
forman las componentes de untensor. Diciendo esto, se quiere decir que (al igual que con los
esfuerzos), si cambiamos de eje de coordenadasi1, i2, i3 por otros ejesi′1, i′2, i
′3, tenemos con
la notación de doble sufijo,
ε′ij = likljlεkl. (8.26)
Además, este tensor de deformaciones también es simétrico, es decir
εij = εji. (8.27)
También es importante notar que podemos escribir las ecuaciones (8.24), en la siguiente
forma:
εij =1
2(ui,j + uj,i) . (8.28)
46
Apéndice B
9.4. Otros resultados de las simulaciones
Para un agujero de diametro 10mm
Fig. 9.4: Comportamiento convergente de la función desempeño. Las lineas rojas represetan los valores reales utilizados en las simulaciones
para obtener las deformaciones experimentales.
47
Para un agujero de diametro 5mm
Fig. 9.5: Comportamiento convergente de la función desempeño. Las lineas rojas represetan los valores reales utilizados en las simulaciones
para obtener las deformaciones experimentales.
48
49
Apéndice C 10.1. Ventana experimental
Este es el código que se utilizo para obtener los campos de deformación
experimental, para utilizarlo basta que al final de realizar una simulación den ANSYS se
escriba el siguiente comando:
"\input, ventana (nombre el archivo), txt (la extensión)"
Código:
/POST1
PP1X=12 !X punto inf. izq de la ventana MODIFICABLE
PP1Y=20 !Y punto inf. izq. de la ventana MODIFICABLE
LV=10 !Largo de la ventana en dir x MODIFICABLE
AV=10 !Ancho en dir Y MODIFICABLE
P2X=PP1X+LV
P2Y=PP1Y
DIV=AV/9 ! NUMERO DE DIVISIONES EN VENTANA
50
*DO,K,1,10
PATH,CAMINO%K%,2,7,15 !DEFINE EL CAMINO
(2PUNTOS,6VARIABLES,10DIVISIONES)
PPATH,1,,PP1X,PP1Y,0 !PUNTO 1 DEL CAMINO
PPATH,2,,P2X,P2Y,0 !PUNTO 2 DEL CAMINO
PDEF,DESX,EPEL,X !DEFINIR VARIABLES A INTERPOLAR
DESPLAZAMIENTO EN X EPEL = deformaciones S= Stress
PDEF,DESY,EPEL,Y
PDEF,DESXY,EPEL,XY !DEFINIR VARIABLES A INTERPOLAR
DESPLAZAMIENTO EN Y
PAGET,RES%K%_,TABLE !COPIA VARIABLES EN TABLA
*CFOPEN,placa25x50y30,txt,,APPEND !COPIA RESULTADOS EN ARCHIVO DE TEXTO
*VWRITE,RES%K%_(1,1),RES%K%_(1,2),RES%K%_(1,5),RES%K%_(1,6),RES%K%_(1,7)
(F20.15,',',F20.15,',',E30.18,',',E30.18,',',E30.18)
*CFCLOSE
PP1Y=PP1Y+DIV !PASA A SIGUIENTE CAMINO
P2Y=P2Y+DIV
*del,,prm_ !borra los parámetros viejos
*ENDDO
51
Apéndice D 11.1. Algoritmo Genético.
El algoritmo genético que se presenta a continuación es una archivo.m de MatLab.
Para correr esta función se utiliza el siguiente comando en la ventana principal de MatLab:
“x=algo_genxy(Deformaciones, # de iteraciones, Carga aplicada, E, Poisson, Xmax, Ymax,
Tamaño de la población)”
Donde se debe manejar el mismo sistema de unidades y la variable Deformaciones
son las deformaciones que se obtienen a partir de la ventana descrita anteriormente. Por
otro lado, la función va generando las graficas con los resultados correspondientes y cuando
se cumple su criterio de detención devuelve la información del elemento de la última
población que obtuve mejor valor de la función desempeño.
Algoritmo Genético – Código en MatLab
function x=algo_genxy(def,iter,T,E,mu,xmax,ymax,tampoblacion) %Los valores deseados se deben ir cambiando segun sea el caso adeseado=15; xdeseado=20; ydeseado=20; msim=momentosxy(def); [h,j]=size(def); ancho=round(abs(def(1,1)-def(h,1))); alto=round(abs(def(1,2)-def(h,2)));
52
nven=round(sqrt(h)); po=poblacion_inicial(xmax,ymax,tampoblacion); figure(1); subplot(2,2,1),ylabel('Funcíon de desempeño'),xlabel('Poblacion No.') subplot(2,2,2),ylabel('Tamaño Agujero'),xlabel('Poblacion No.') subplot(2,2,3),ylabel('Posicion x'),xlabel('Poblacion No.n') subplot(2,2,4),ylabel('Posicion y'),xlabel('Poblacion No.'), pause(1e-6); i=1; while (i<iter) fo=fun_objxy(po,msim,T,E,mu,alto,ancho,nven); rul=ruleta(fo); %funcion(i)=max(fo(:,1)); funcion(i)=rul(1,1); funx(i)=rul(1,3); funy(i)=rul(1,4); funa(i)=rul(1,2); po=nueva_poblacion(rul(:,2:4),xmax,ymax,tampoblacion); figure(1); subplot(2,2,1), hold on, plot(funcion) subplot(2,2,2), hold on, plot(funa); subplot(2,2,3), hold on, plot(funx); subplot(2,2,4), hold on, plot(funy) , pause(1e-6); hold on if (funcion(i)>=0.9999) i=iter; end i=i+1; disp(rul(1,:)) end hold on figure(2) subplot(2,2,1), hold on, plot(funcion,'k'), ylabel('Función de Desempeño','FontSize',14),xlabel('Poblacion No.','FontSize',14),title('Parámetro 0.55','FontSize',14),grid on subplot(2,2,2), plot(funa,'k') ,ylabel('Tamaño Agujero','FontSize',14),xlabel('Poblacion No.','FontSize',14),title('Parámetro 0.55','FontSize',14),grid on subplot(2,2,3), plot(funx,'k') ,ylabel('Posicion x','FontSize',14),xlabel('Poblacion No.','FontSize',14),title('Parámetro 0.55','FontSize',14),grid on subplot(2,2,4), plot(funy,'k') ,ylabel('Posicion y','FontSize',14),xlabel('Poblacion No.','FontSize',14),title('Parámetro 0.55','FontSize',14),grid on ,pause(1e-6); hold on figure(1) subplot(2,2,2), hold on ,plot([1,iter],[adeseado,adeseado],'r'); subplot(2,2,3), hold on ,plot([1,iter],[xdeseado,xdeseado],'r'); subplot(2,2,4), hold on ,plot([1,iter],[ydeseado,ydeseado],'r'); x=rul(1,:); end %Momentos experimentales: %Los datos deben venir en este orden x,y,ex,ey,exy function B=momentosxy(datos)
53
x=datos(:,1); y=datos(:,2); ex=datos(:,3); ey=datos(:,4); exy=datos(:,5); deltax=abs(x(1)-x(2)); deltay=deltax; momento_00=sum(ex*deltax*deltay.*ey); momento_10=sum((x.*ex)*deltax*deltay.*ey); momento_01=sum((y.*ex)*deltax*deltay.*ey); momento_11=sum(((x.*ex).*y)*deltax*deltay.*ey); %momento_21=sum((x.*(x.*ex).*y)*deltax*deltay.*exy.*ey); %momento_12=sum((y.*(x.*ex).*y)*deltax*deltay.*exy.*ey); B(1)=momento_00/momento_00; B(2)=momento_10/momento_00; B(3)=momento_01/momento_00; %B(4)=momento_11/momento_00; %B(5)=momento_21; %B(6)=momento_12; B=B'; end function [A]=poblacion_inicial(x,y,tampo) for i=1:tampo A(i,2)=x*rand(); A(i,3)=y*rand(); amax=sqrt(A(i,2)^2+A(i,3)^2); %para asegurar que el el hueco no este fuera de los valores de x y y A(i,1)=amax*rand(); end end
function [FO]=fun_objxy(po,msim,T,E,mu,alto,ancho,nven) [n,m]=size(po); FO=zeros(n,1); for i=1:n A=ventana(po(i,1),po(i,2),po(i,3),alto,ancho,nven,T,E,mu); mteo=momentosxy(A); error=(mteo-msim); cerror=error.*error; sumerror=sqrt(sum(cerror)); FO(i,1)=0.55/(0.55+sumerror); end FO=[FO,po]; end % la idea de esta funcion es generar la ventana teorica para cada punto function D=ventana(a,x0,y0,alto,ancho,n,T,E,mu) deltax=ancho/(n-1); deltay=alto/(n-1); C=zeros(n*n,3); Z=zeros(n*n,3); %A contiene las posiciones en x,y for i=1:n for j=1:n A(j+n*(i-1),1)=x0+(j-1)*deltax;
54
A(j+n*(i-1),2)=y0; end y0=y0+deltay; end B=erreteta(A(:,1),A(:,2)); %B contiene los r, teta para los x,y de la matriz A long=n*n; for i=1:long sr=sigmar(T,a,B(i,1),B(i,2)); steta=sigmateta(T,a,B(i,1),B(i,2)); srteta=taorteta(T,a,B(i,1),B(i,2)); C(i,1)=er(sr,steta,E,mu); C(i,2)=eteta(sr,steta,E,mu); C(i,3)=erteta(srteta,E,mu); end %C contiene las deformaciones para los r,teta de la matriz B for i=1:long Z(i,:)=rotacion(C(i,:),B(i,2)); end D=[A,Z]; function x=sigmar(T,a,r,teta) x=((T/2)*(1-(a/r)^2))+((T/2)*(1+(3*(a/r)^4)-(4*(a/r)^2))*cos(2*teta)); end function x=sigmateta(T,a,r,teta) x=(T/2)*(1+(a/r)^2)-(T/2)*(1+3*(a/r)^4)*cos(2*teta); end function x=taorteta(T,a,r,teta) x=-(T/2)*(1-3*(a/r)^4+2*(a/r)^2)*sin(2*teta); end function x=er(sr,steta,E,mu) x=(sr-(mu*steta))/E; end function x=eteta(sr,steta,E,mu) x=(steta-(mu*sr))/E; end function x=erteta(taorteta,E,mu) x=taorteta*(1+mu)/E; end function x=er(sr,steta,E,mu) x=(sr-(mu*steta))/E; end % Para pasar a coordenadas cartesianas a partir de las cilindricas % entra el vector en forma de fila y sale en forma de fila, teta es el % angulo de rotacion function [y]=rotacion(z,teta)
55
Q=[(1+cos(2*teta))/2 (1-cos(2*teta))/2 sin(2*teta);(1-cos(2*teta))/2 (1+cos(2*teta))/2 -sin(2*teta); -sin(2*teta)/2 sin(2*teta)/2 cos(2*teta)]; y=inv(Q)*(z'); y=y'; end % aqui enta la poblacion con la funcion objetivo en su primera componente function A=ruleta(poblacion) [n,m]=size(poblacion); if n<=4 disp('¡La población es muy pequeña para realizar ruleta!'); else organizados=sortrows(poblacion,[-1]); copia=organizados; A(1,:)=copia(1,:); A(2,:)=copia(2,:); organizados=[copia(3:n,:)]; suma=sum(organizados(:,1)); suma1=suma*rand(); i=1; copia2=organizados; temp=0; while(i<=n-2) temp=temp+organizados(i,1); if(temp>=suma1) p1=organizados(i,:); organizados=zeros(n-3,m); organizados=[copia2(1:(i-1),:);copia2((i+1):n-2,:)]; i=n-2; end i=i+1; end suma=sum(organizados(:,1)); suma2=suma*rand(); i=1; temp=0; while(i<=n-3) temp=temp+organizados(i,1); if(temp>=suma2) p2=organizados(i,:); i=n-3; end i=i+1; end A(3,:)=p1(1,:); A(4,:)=p2(1,:); end end %esta funcion crea una nueva poblacion dados los padres, usando los cruces %y mutuaciones function A=nueva_poblacion(padres,xmax,ymax,n) p1=padres(1,:); p2=padres(2,:);
56
p3=padres(3,:); p4=padres(4,:); %primeros hijos, cruzando los mejores padres A1=primeros_hijos(p1,p2); [m,p]=size(A1); n=n-m; %Ahora mutando el primer padre v1=round(rand()); if(v1==0) singo=1; else singo=-1; end p1max=sqrt(p1(2)^2+p1(3)^2); mp1=p1; mp1(1)=mp1(1)*(1+singo*0.1); if (mp1(1)>p1max || mp1(1)<=0) mp1(1)=4+p1max*rand(); end %cruces fuertes A=[A1;mp1]; A=cruces(A,p3,p2,n-6,xmax,ymax); end function A=primeros_hijos(p1,p2) c1=p1; c1(2)=p2(2); amax1=sqrt(c1(2)^2+c1(3)^2); if(c1(1)>amax1) c1(1)=amax1*rand(); end c2=p1; c2(3)=p2(3); amax2=sqrt(c2(2)^2+c2(3)^2); if(c2(1)>amax2) c2(1)=amax2*rand(); end c3=p2; c3(2)=p1(2); amax3=sqrt(c3(2)^2+c3(3)^2); if(c3(1)>amax3) c3(1)=amax3*rand(); end c4=p2; c4(3)=p1(3); amax4=sqrt(c4(2)^2+c4(3)^2); if(c4(1)>amax4) c4(1)=amax4*rand(); end A=[p1;c1;c2;c3;c4]; end function [z]=cruce_inter(p3,p4,xmax,ymax) x=p3; y=p4; v1=rand(1,3); z=zeros(1,3);
57
for i=1:3 z(i)=x(i)+v1(i)*(x(i)-y(i)); end if(z(1)>(sqrt(z(2)^2+z(3)^2)) ||z(1)<0) z(1)=v1(1)*sqrt(z(2)^2+z(3)^2)/2; end if (z(2)>xmax ||z(2)<0) z(2)=xmax*v1(2); end if(z(3)>ymax || z(3)<0) z(3)=ymax*v1(3); end end function [z]=cruce_lineal(p3,p4,xmax,ymax) x=p3; y=p4; v1=rand(); z=x+(v1*(x-y)); if(z(1)>(sqrt(z(2)^2+z(3)^2)) ||z(1)<0) z(1)=v1*sqrt(z(2)^2+z(3)^2)/2; end if (z(2)>xmax ||z(2)<0) z(2)=xmax*v1; end if(z(3)>ymax || z(3)<0) z(3)=ymax*v1; end end function z=mutacion(p,xmax,ymax,por) z=p; posicion=round(2*rand+1); v1=round(rand()); if(v1==0) singo=1; else singo=-1; end if(posicion==1) z(1)=z(1)*(1+singo*por); end if(posicion==2) z(2)=z(2)*(1+singo*por); end if(posicion==3) z(3)=z(3)*(1+singo*por); end %prueba que estén dentro de los rangos if(z(1)>(sqrt(z(2)^2+z(3)^2)) ||z(1)<0) z(1)=rand*sqrt(z(2)^2+z(3)^2); end if (z(2)>xmax ||z(2)<0) z(2)=xmax*rand; end if(z(3)>ymax || z(3)<0) z(3)=ymax*rand; end end
Bibliografía
[1] SADD, Martin H. ELASTICITY, Theory, Applications, and Numerics.ELSEVIER, But-
terworth Heinemann 2005.
[2] R.C. HibbelerMechanics of materials. 6th Edition, Prentice Hall, New Jersey 2004.
[3] LAZARUS, Teneketzis Tenek y ARGYRIS, John.Finite Element Analysis for Composite
Structures. 1st Edition, Kluwer Academic Publishers, Dordrecht, Boston y London, 1998.
[4] MARAÑON, A. RUIZ, P.D. NURSE, A.D. HUNTLEY, J.M. RIVERA, L. ZHOU, G.
Identification of subsurface delaminations in composite laminatesWolfson School of
Mechanical and Manufacturing Engineering, Lougborough, October 18, 2006.
[5] SILVA, Felipe. Estructuras, formas y materiales topológicamente óptimos.Tesis de
Matemáticas, Departamento de la Universidad de los Andes, 2008.
[6] GOSZ, Michael R.Finite Element Method. Applications in Solids, Structures, and Heat
Transfer. 1st Edition, Taylor and Francis, Boca Raton, London and New York
[7] ANSYS Help.Release 10.0 Documentation for ANSYS. 2005 SAS IP, Inc.
[8] STEWART, James.Calculus EARLY TRANSCENDENTALS. Brooks/Cole 3rd. Edition,
1995.
58
[9] COLEY David A. An Introduction to Genetic Algorithms for Scientist and Enginners.
Singapore, River Edge NJ. World Scientific. 1999.
[10] HURTADO, Carlos Arturo.Caracterización de Delaminación en un Material Com-
puesto Empleando Modelaje por Vigas de Euler y Algoritmos Genéticos.Proyecto de
grado de Ingeniería Mecánica. 2006
[11] GEATbx: Genetic and Evolutionary Algorithm Toolbox for use with MATLAB Docu-
mentation. Version 3.70 (released November 2005) Autor: Hartmut Pohlheim.
[12] Genetic Algorithm and Direct Search ToolboxFor Use with MATLAB. User´s Guide,
Version 1 The MathWorks. June 2004.
59
Top Related