Aplicación de la ecuación integral de contorno axisimétrica para elasticidad y termoelasticidad.

Yomar A. González¹, Adrian P. Cisilino², Miguel E. Cerrolaza¹

¹ Instituto Nacional de Bioingeniería, Universidad Central de Venezuela, Caracas, Venezuela, yomar.gonzalez@inabio.edu.ve, http://www.ucv.ve/cebio.htm.

² División de Soldadura y Fractomecánica, Instituto Tecnológico de Materiales - INTEMA, Universidad Nacional de Mar del Plata, Mar del Plata, Argentina, cisilino@fi.mdp.edu.ar, http://intema.fi.mdp.edu.ar/.

Resumen

Para el tratamiento de problemas con simetría axial, la ecuación integral de contorno para elasticidad y/o termoelasticidad en 3D sufre una transformación de coordenadas y un proceso de integración adicional para generar las soluciones axisimétricas requeridas.

Este trabajo presenta la implementación de las soluciones axisimétricas para el estudio de problemas definidos en un medio elástico-lineal, isótropo, homogéneo y en régimen estacionario, en un código de Elementos de Contorno multi zona. Las expresiones analíticas de las soluciones y el desarrollo asintótico de las mismas fueron extraídas de la literatura. Para mantener un alto grado de exactitud en la solución usando elementos lineales isoparamétricos, algunos términos singulares contenidos en la ecuación integral de contorno fueron integrados analíticamente. La nueva herramienta  es aplicada a diferentes ejemplos con solución analítica conocida, cotejando los resultados y comprobando la precisión y eficacia del código.

Palabras claves: elementos de contorno, simetría axial, elasticidad, soluciones fundamentales, termoelasticidad, conductividad térmica.

Bem axisymmetric formulation for elasticity and thermoelasticity.

ABSTRACT

Dealing with axisymmetric problems involves a co-ordinate transformation system and additional integration process of the well-known classic boundary integral equation for 3D elasticity and thermoelasticity theories.

An axisymmetric boundary element formulation is implemented in a steady-state framework for multi-domain elastic, isotropic and homogeneous problems. The fundamental solutions and their singular behaviour were taken from available literature. In order to obtain a high accuracy in the solution, employing continuous linear elements, some singular terms appearing in the boundary integral equation have been integrated analytically. Several numerical examples are provided to illustrate de versatility of this approach for practical engineering problems.

Key words: Boundary element, axisymmetry, elasticity, fundamental solutions, thermoelasticity, thermal conductivity.

Recibido: 30/04/07  Revisado: 07/08/07 Aceptado: 20/08/07

1. IntroducCión

Este trabajo presenta un código computacional basado en un método de frontera capaz de simular y caracterizar el estado de tensión-deformación de un sólido. En general, los algoritmos basados en el método de elementos de contorno se encuentran disponibles para aplicaciones a problemas planos o tridimensionales. Sin embargo para  los problemas con simetría axial, la formulación 3D para elasticidad y/o termoelasticidad sufre una transformación de coordenadas y un proceso de integración adicional para generar las soluciones axisimétricas requeridas, las cuales suelen expresarse en términos de las integrales elípticas completas de primer y segundo orden. El comportamiento asintótico de éstas últimas es determinante para establecer el procedimiento de integración más adecuado en los casos singulares.

A corto y mediano plazo, se pretende utilizar el código como una herramienta numérica capaz de simular procesos biológicos complejos (Ej. consolidación ósea), los cuales se caracterizan por ser no lineales y altamente heterogéneos, siguiendo las teorías de comportamiento comúnmente utilizadas en la literatura especializada, con el subsecuente beneficio del ahorro en tiempo computacional asociado al método.

Afortunadamente, el método de elementos de contorno cuenta con las soluciones fundamentales cuya implementación permitiría a mediano y a largo plazo,  simular este tipo de fenómenos complejos, no lineales y altamente heterogéneos. Una primera aproximación conservadora sería considerar modelos axialmente simétricos en el campo elástico-lineal, para lo cual se hace necesario la incorporación de la formulación axisimétrica presentada por Balas, Slâdek y Slâdek (1989), siguiendo la metodología propuesta por Graciani et al. (2005) en un código multi zona basado en el Método de los Elementos de Contorno desarrollado por Beer (2001). Más adelante, se pretende complementar el código, incorporando nuevas formulaciones y metodologías para la evaluación e integración de variables adicionales más complejas características del fenómeno biológico dentro de la ecuación de contorno.

La nueva herramienta es aplicada a diferentes ejemplos con solución analítica conocida, cotejando los resultados y comprobando la precisión y eficacia de los códigos.

2. FORMULACIÓN AXISIMÉTRICA EN ELEMENTOS DE CONTORNO

En esta sección, se presenta en forma general la ecuación integral de contorno axisimétrica para los desplazamientos correspondiente a un punto de colocación ubicado en el contorno Г de la geometría.

Se asume que el plano B en dos dimensiones (figura 1a) rota 360° con respecto al eje z. Esto permite formar una geometría axisimétrica bajo una carga también axisimétrica, donde r y φ son las direcciones radial y angular respectivamente. Se puede afirmar que tanto los desplazamientos como las tensiones calculadas son constantes para cualquier valor de φ, por lo que esta clase de problemas pueden ser analizados eficientemente, con tan solo considerar al cuerpo B en vez de todo el dominio tridimensional. Como resultado de esta simetría axial, las direcciones r y z son suficientes para definir el problema.

Figura 1. (a) Definición del contorno G y de la región B en un problema tridimensional, (b) Definición del sistema de coordenadas cilíndricas (r,φ,z), (c)  Representación 3D de un cuerpo (Ω+Г), (d)  Región genérica (Ω* + Г*) donde la carga unitaria es aplicada, con las mismas propiedades elásticas que (Ω+Г). También se evidencia la relación existente entre el punto fuente “P” y el punto de integración “Q”.   (Balas et al., 1989)

La formulación para la elasticidad lineal axisimétrica ha sido estudiada y desarrollada, entre otros, por Kermandis (1975), Cruse et al. (1977), y Rizzo y Shippy (1979). Actualmente De Lacerda et al. (2000) ha trabajado en la aplicación de ésta formulación en problemas de contacto, obteniendo excelentes resultados (Graciani et al., 2005). Por otra parte, Cruse (1975) y Rizzo y Shippy (1977), fueron los primeros en implementar con ciertas limitaciones el método de la ecuación integral de contorno para problemas termoelásticos con simetría axial (usualmente no acoplado), en los cuales la distribución de temperatura incide directamente sobre el campo tensión-deformación del problema. Seguidamente, Bakr y Fenner (1983), derivaron las ecuaciones explícitas axisimétricas formuladas para que cada nodo del contorno posea tres grados de libertad y así combinar tanto el análisis térmico como el elástico.

Se puede encontrar en la literatura (Bakr y Fenner, 1983; Brebbia et al., 1984; Dargush y Banerjee, 1992; Brebbia et al., 2003; Banerjee, 1994 y Graciani et al., 2005) un análisis detallado acerca de la derivación de la ecuación integral de Somigliana tomando en cuenta la simetría axial. Consecuentemente, tan solo se hará mención en el texto a las consideraciones más relevantes en el desarrollo matemático antes de presentar las soluciones fundamentales y su respectivo análisis en casos de integración singular.

La deducción de la ecuación integral de los desplazamientos (campo elástico-lineal) para un punto de colocación en el contorno axisimétrico G, comienza rescribiendo la expresión tridimensional del núcleo de desplazamientos (1), deducida considerando un problema genérico como el mostrado en la figura 1.c, en términos del nuevo sistema en coordenadas cilíndricas (2). La nueva ecuación  ahora puede ser integrada a lo largo de una trayectoria circular (anillo de carga unitaria)  con respecto al eje de simetría considerado (4) (Balas, Slâdek y Slâdek, 1989; Brebbia et al., 2003). De manera análoga se procede con la ecuación que describe el campo de tracciones.

, Siendo        (1)

    (2)

      (3)

T es la matriz que provee la transformación de coordenadas e Y relaciona las direcciones de las cargas unitarias en el punto P, con las correspondientes direcciones en el nuevo sistema de coordenadas. El nuevo campo de desplazamientos viene dado entonces por

      (4)

2.1 Soluciones fundamentales para elasticidad

La ecuación integral de contorno axisimétrica que describe el problema general elástico es de la forma

     (5)

Donde P es el punto del anillo donde la carga unitaria es aplicada, Q es el punto de integración, C_αη (P) son constantes las cuales dependen de la geometría en P , r (Q) es la distancia radial desde el eje de simetría hasta el punto de integración, u_r y u_z son los desplazamientos radiales y axiales respectivamente, t_r y t_z, son las tracciones radiales y axiales respectivamente,   y  y son los núcleos axisimétricos para los  

desplazamientos y tracciones respectivamente. (De Lacerda et al., 2000).

Haciendo uso de integrales particulares (Balas, Slâdek y Slâdek, 1989) y tomando en cuenta el término 2πr (Q) e introduciendo notaciones especiales (6), se pueden escribir las expresiones explícitas (7-10) para los desplazamientos elásticos en términos de r, φ y z (coordenadas cilíndricas) considerando el caso general donde R>0.

    (6)

En las expresiones (6), R y r representan al radio del punto de colocación e integración respectivamente, k es el argumento para la integral elíptica de primer orden “K” y de segundo orden “E”, Z y z, son las coordenadas axiales de los puntos de colocación e integración respectivamente.

(*) Comparando esta expresión (Balas, Slâdek y Slâdek, 1989) con la presentada por Brebbia et al. (1984), se hace necesario que el argumento k quede elevado a una potencia cuadrada.

   (7)

   (8)

    (9)

  (10)

(**) La notación empleada aquí es la transpuesta a la presentada por Balas, Slâdek y Slâdek (1989).

Las ecuaciones arriba descritas están deducidas a partir de su expresión general en términos de las funciones de Legendre de orden cero y de sus derivadas, las cuales, a su vez, suelen ser escritas en función de las integrales elípticas completas de primer y segundo orden  y .

Sustituyendo las expresiones (7-10) en las ecuaciones de las tensiones expresadas en coordenadas cilíndricas y tomando en cuenta el vector normal unitario en el punto de integración Q, se obtienen las tracciones (11-14) sobre el contorno G .

En las ecuaciones (7-14), G y υ representan a los módulos de corte y de Poisson respectivamente y sus valores varían según las regiones asumidas y, n_r y n_z, son las componentes radial y axial del vector normal unitario definido en el punto de integración Q.

   (11)

   (12)

  (13)

       (14)

De esta forma y      corresponden respectivamente, a las componentes axial y radial de los 

desplazamientos y del vector tensión  asociado a un plano cuya normal es n, provocados en un punto Q, cuando se aplican un anillo de carga unitaria, según la dirección axial y radial, en el punto P. (Graciani, 2006).

Cuando la coordenada radial del punto de colocación es cero (R = 0), las expresiones para los desplazamientos y las tracciones se reducen a las ecuaciones (15) y (16) respectivamente.

   (15)

        (16)

(***) La notación empleada aquí es la transpuesta a la presentada por Balas, Slâdek y Slâdek (1989) como consecuencia del cambio observado en (8) y (9).

2.2 Soluciones fundamentales para termo-elasticidad

Así como en el caso elástico, la formulación termoelástica tridimensional definida usualmente en un sistema cartesiano (Banerjee, 1994), se transforma a un sistema de coordenadas cilíndrico, dando lugar a los ya conocidos núcleos axisimétricos elásticos y a las nuevas soluciones para el campo de potencial, términos de acople y sus derivadas, quedando la ecuación integral (5) redefinida según (17), donde  P es el punto del anillo donde la carga unitaria es aplicada, Q es el punto de integración, C_αη (P) son constantes las cuales dependen de la geometría en P, _γ (Q) es la distancia radial desde el eje de simetría hasta el punto de integración. u_r , u_z y θ son los desplazamientos en dirección radial y axial y la temperatura asociada al punto de interés. t_r , t_z y q son las tracciones y el flujo de temperatura respectivamente.

   (17)

Al igual que en el desarrollo previo, el termino 2πr (Q) , el cual aparece después del proceso de integración en la dirección φ, es absorbido por los s respectivos. Los subíndices  a y h toman los valores de r y z, según sea el caso. Las expresiones  y de manera análoga, los términos correspondientes a las 

tracciones en la ecuación (17),  son idénticos a aquellos definidos para elasticidad axisimétrica descritos en el apartado 2.1 y presentados e implementados en González et al. (2006).   y  son los núcleos 

axisimétricos relacionados con el flujo de potencial (Banerjee, 1994; Bakr y Fenner, 1983).  

son los términos de acople (Banerjee, 1994; Bakr y Fenner, 1983). Estas soluciones fundamentales fueron derivadas considerando un medio homogéneo en un problema termoelástico en  régimen estacionario.

Manteniendo las notaciones especiales (6) e introduciendo otras nuevas (18), se pueden escribir las expresiones explícitas (19,20) para los términos de acople y flujo de potencial en términos de r, φ y z considerando el caso general donde R > 0 (Banerjee, 1994). Las notaciones (6) y (18) son relaciones matemáticas que buscan simplificar convenientemente la formulación a la hora de su implementación.

   (18)

R, Z, r y z representan las coordenadas del punto de colocación e integración respectivamente, λ y μ son las constantes de Lamé y β es un parámetro termoelástico definido en función del coeficiente de expansión térmica α. Por conveniencia, definamos ahora la distancia entre P y Q como ρ y el coeficiente de las integrales elípticas como m, manteniendo su expresión original, dejando a k como la notación para la conductividad térmica lineal.

  (19)

  (20)

n_r y n_z representan las componentes del vector normal unitario definido en el punto Q .

3. INTEGRACIÓN DE LAS SOLUCIONES FUNDAMENTALES

Sí el punto P está fuera del elemento “t” de integración, los integrando en (5) son regulares y corresponden a la evaluación de las expresiones (7-14), por lo que las integrales pueden ser computadas numéricamente mediante la implementación de la cuadratura de Gauss-Legendre. Sí por el contrario, el punto P coincide con alguno de los nodos del elemento en integración, algunos de los integrando son singulares y deben ser evaluados en base a las expansiones singulares correspondientes. Las integrales numéricas relacionadas a estos términos deben ser calculadas con una cuadratura especial o método correspondiente, por ejemplo, en el caso de la integración de los términos   y   , aunque su comportamiento es regular, ameritan de un 

método de integración numérica más refinado y/o adaptativo, como por ejemplo el Método Adaptativo de Romberg. Para el caso de   dada la función logarítmica que la acompaña, se emplea la cuadratura 

especial de Gauss-Laguerre. (Beer, 2001).

Para las soluciones termo-elásticas, cuando los términos, que acompañan a K (m) se vuelven singulares, se pueden expandir los núcleos alrededor del punto P, obteniendo así las expresiones a 

ser evaluadas e integradas para estos casos. Otra alternativa menos precisa, la constituye la implementación del Método Adaptativo de Romberg.

A pesar de que la integración numérica de los núcleos elásticos sobre el elemento singular amerita 

especial atención, existe un método alternativo para evitar este paso basado en la combinación de dos técnicas: Rigid Body Motion en dirección z e Inflation Mode para axisimetria en dirección r, obteniendo así las integrales impropias más la influencia de los coeficientes C_an con gran exactitud (Banerjee, 1994; De Lacerda et al., 2000). En el caso de que R = 0 y cuando , las técnicas mencionadas anteriormente no se aplicaron, debido a las restricciones en los desplazamientos para puntos de colocación en el eje de simetría, en su lugar, los valores correspondientes en la diagonal de la matriz H del sistema inicial se obtuvieron 

mediante la integración por cuadratura de Gauss de los resultados producto de las expansiones, más la adición del único término libre distinto de “cero”, C_zz (Graciani et al., 2005). Para ello, se consideró la singularidad como débil. De esta manera el método sería aplicable a cualquier problema elástico axialmente simétrico.

En el caso de los núcleos , los coeficientes desconocidos h_i,j (j= 3,6,9,…) para cada bloque elemental 3x3  

de la matriz H del sistema   se obtuvieron planteando un problema con solución analítica 

conocida: cilindro sólido libre de tensiones sometido a un campo de temperatura constante, donde  ,  

y , siendo  x y  h  las coordenadas radiales y axiales de los nodos del problema. Las tracciones/flujos

 relacionadas con este campo de desplazamientos y temperatura son nulos.

Basados en las expresiones para todas las relaciones elásticas deducidas anteriormente, el método analítico también se presenta como una excelente opción para integrar con gran precisión los términos singulares de las expansiones, en aquellos casos donde se empleen elementos lineales continuos (Graciani et al., 2005).

4. IMPLEMENTACIÓN NUMÉRICA

Las expresiones mencionadas anteriormente, fueron implementadas siguiendo el algoritmo propuesto por Beer (2001), según la metodología propuesta por Graciani et al. (2005). El código original (Beer, 2001) está estructurado en una secuencia de subrutinas que permiten la obtención de los campos de desplazamiento/temperatura y tensión/flujo en el rango elástico-lineal según sea el caso. Para la implementación de la nueva formulación, se realizaron cambios en las subrutinas más importantes (aquellas relacionadas con la evaluación e integración de los núcleos para el ensamblaje del sistema de ecuaciones algebraicas), manteniendo la estructura inicial del algoritmo. Los núcleos fueron ajustados en función de la metodología utilizada por Beer (2001) para el proceso de integración y ensamblaje de la matriz del sistema. El cálculo se realiza por región, integrando cada elemento τ de la región n con respecto a cada uno de los puntos de colocación i, tratando por separado los casos regulares y singulares, tanto para  como para . Los resultados de la integración son guardados en arreglos elementales que posteriormente, y según las condiciones de contorno impuestas, serán ensamblados directamente en la matriz del sistema A (Ax=B). Después de que cada matriz elemental sea ensamblada, se adiciona el termino libre C_αη correspondiente, calculados a través de las metodologías mencionadas anteriormente (Graciani et al., 2005). Para resolver el sistema de ecuaciones por región se utilizó un método de aproximación iterativo. Este proceso se repite tantas veces como regiones existan en el modelo, en función de las condiciones de contornos particulares y de la imposición de las relaciones de compatibilidad (que permitirían incluir el efecto de una región sobre otra) sobre las superficies de contacto (de ser el caso). Las modificaciones fueron implementadas y verificadas en González et al. (2006).

La adición de los nuevos términos a la ecuación integral de contorno (núcleos termo-elásticos), no reportó mayores inconvenientes, manteniendo la estructura del código descrita anteriormente. Tan solo fue necesario expandir las dimensiones de los arreglos involucrados en los procesos de evaluación, integración y ensamblaje de todo el conjunto de expresiones que acopladas, reproducen el comportamiento termoelástico del material.

5. RESULTADOS NUMÉRICOS PARA ELASTICIDAD CONSIDERANDO UN MEDIO ISÓTROPO Y HOMOGÉNEO

5.1 Esfera hueca

Consideremos una esfera hueca con radio interno R1 y externo R2. La superficie interna se mantiene a una presión uniformemente distribuida P1 y  a la superficie externa se le aplica una presión P2. El modelo consistió de 41 elementos isoparamétricos con aproximación lineal con el objeto de representar con mayor precisión la geometría curva del problema y consecuentemente el vector normal asociado a los elementos a los cuales se les impuso la condición de carga. Los valores numéricos empleados son (Bakr et al., 1983)

Figura 2. (a) esquema representativo del problema simétrico en 3D: problema elástico, (b) discretización para BIE. Se puede observar el sentido y la numeración utilizada para los nodos.

Para este caso se cuenta con la solución analítica correspondiente (Timoshenko et al., 1951):

     (21)

      (22)

Tabla 1. Resultados numéricos para el problema de la esfera hueca usando una malla de  41 elementos (figura 2.b). u_R y σ_RR son los desplazamientos y  tensiones respectivamente en dirección R, radial al centro de la esfera, mientras que  σ_tt es la tensión en dirección tangencial  a las superficies de la esfera.

 

(*) Valor promedio en base a los resultados obtenidos en los nodos concernientes.

(°) Valor calculado en base a un sistema de coordenadas local al elemento (τ,φ,n)

(^)Valor calculado según el estado plano de tensiones en una partícula del sólido, rotado a un nuevo sistema de referencia (R,τ):

5.2 Esfera hueca encamisada

Consideremos ahora un problema similar al caso anterior, pero la esfera hueca (radio interno, R1) se encuentra ahora embebida dentro de un caso esférico de radio externo R2. Por lo que el problema puede ser analizado como multi-región. Llamemos Rm, al radio en la interfase donde no se consideran los efectos de interferencia. La superficie interna de la esfera hueca se encuentra bajo una presión uniformemente distribuida P1, mientras que la superficie externa del casco esférico se encuentra sometida a una presión P2. La esfera interna posee propiedades E1 y υ, mientras que las del casco externo, son E2 y  υ. La geometría se modeló con dos regiones interconectadas, con un total de 128 elementos isoparamétricos con aproximación lineal (figura 3.b). Los valores numéricos empleados son:

Figura 3. (a) esquema representativo del problema simétrico multiregión en 3D: problema elástico, (b) discretización para BIE. Se puede observar el sentido y la numeración utilizada para los nodos.

Las expresiones analíticas de este problema se presentan a continuación (Timoshenko et al., 1951), tomando a  pm= :1.364946

   (23)

  (24

   (25)

Tabla 2. Resultados numéricos para el problema de la esfera multi-región (figura 3.b). u_R y σ_RR  son los desplazamientos y  tensiones respectivamente en dirección radial, mientras que  σ_φφ  es la tensión en dirección circunferencial. Los puntos internos se localizan en ambas regiones del sistema.

 

(*) Valor promedio en base a los resultados obtenidos en los nodos concernientes.

(°) Valor calculado en base a un sistema de coordenadas local al elemento (τ,φ,n)

(^)Valor calculado según el estado plano de tensiones en una partícula del sólido, rotado a un nuevo sistema de referencia (R,τ):

6. RESULTADOS NUMÉRICOS PARA TERMOELASTICIDAD CONSIDERANDO UN MEDIO ISÓTROPO Y HOMOGÉNEO

6.1 Cilindro hueco

Consideremos un cilindro hueco con radio interno R1 y externo R2 y longitud H. La superficie interna se mantiene bajo una presión uniformemente distribuida P, mientras que la superficie externa se encuentra libre de tensiones. Las superficies se mantienen a temperatura constante θ. El desplazamiento axial no está restringido en los extremos del cilindro. El modelo consistió en 48 elementos isoparámetricos con interpolación lineal. Los valores numéricos empleados son (Bakr et al., 1983):

Figura 4. (a) esquema representativo del problema simétrico en 3D: problema termo-elástico, (b) discretización para BIE. Se puede observar el sentido y la numeración utilizada para los nodos.

Las expresiones analíticas de este problema se presentan a continuación (Timoshenko et al., 1951):

   (26)

Veamos un sencillo análisis de convergencia tomando como referencia los campos de desplazamiento (ur) y tensiones radiales y circunferenciales (S11 , S22 ). Para ello, se plantearon tres modelos discretizados de manera uniforme en todos sus lados con elementos lineales isoparámetricos. El modelo 1 consistió de 12 elementos, un segundo modelo fue discretizado con 24 elementos y el último modelo consistió de 48 elementos (figura 4.b).

Figura 5. Esquema escalonado de convergencia en el campo de tensiones radiales (S₁₁)                           y circunferenciales (S₂₂)  para el problema del cilindro hueco termoelástico: comparación de los resultados numéricos (BIE) obtenidos a partir de distintas discretizaciones con respecto a la solución analítica exacta.

Figura 6. Esquema de convergencia en el campo de desplazamientos en dirección radial para el problema del cilindro hueco termoelastico: comparación de los resultados numéricos (BIE) obtenidos a partir de distintas discretizaciones con respecto a la solución analítica exacta.

Los resultados numéricos presentados a continuación fueron obtenidos mediante el análisis del modelo 3, en vista de la convergencia presentada por el modelo con una desviación estándar de 0.2865, comparada con los valores de las desviaciones 0.3412 y 0.31 reportados por los modelos 1 y 2 respectivamente. A pesar de que el análisis de convergencia realizado no fue en extremo riguroso, ayudó a clarificar el grado de discretización requerido para obtener una solución aceptable al problema termoelástico planteado.

Tabla 3: Resultados numéricos en el contorno para el problema del cilindro hueco termoelástico usando una malla de 48 elementos (figura 4.b). u_r y σ_φφ son los desplazamientos en dirección radial y las  tensiones circunferenciales, mientras que  σ_eq es la tensión equivalente de Von Mises.

(‘) Valor promedio en base a los resultados obtenidos en los nodos encontrados en la superficie interna o externa del cilindro, según sea el caso.

6.2 Cilindro hueco axialmente infinito

Consideremos un cilindro hueco con radio interno R1 y externo R2 y longitud H. La superficie interna se mantiene a temperatura constante θ1 bajo una presión uniformemente distribuida P, mientras que la superficie externa se encuentra a una temperatura θ2 libre te tensiones. Se asume una condición adiabática en los extremos del cilindro, implicando conducción de calor radial únicamente. Los extremos también están restringidos de movimiento en la dirección axial, por lo que los resultados pueden ser cotejados con las soluciones analíticas para el caso de deformación plana. El modelo consistió de 48 elementos con el objeto de afinar la precisión sobre todo en el cálculo del campo de temperatura y del flujo de potencial, tomando en cuenta el comportamiento de las soluciones analíticas para este problema. Los valores numéricos empleados son (Bakr et al., 1983):

Figura 7.  (a) esquema representativo del problema simétrico en 3D: problema termoelástico, (b) discretización para BIE. Se puede observar el sentido y la numeración utilizada para los nodos.

Las expresiones analíticas de este problema se presentan a continuación (Bakr et al., 1983):

   (27)

Las constantes fueron deducidas a partir de las condiciones de contorno impuestas.

Tabla 4. Resultados numéricos en el contorno para el problema del cilindro hueco axialmente infinito usando una malla de 48 elementos (figura 7). u_r es el desplazamiento en dirección radial, σ_φφ representa el campo de tensiones en dirección circunferencial y σ_eq es la tensión equivalente  de Von Mises. θ_r es el campo de temperatura en dirección radial.

(*) Valor promedio calculado.

6.3 Esfera hueca encamisada con distribución de temperatura

Consideremos ahora el mismo problema del caso 5.2, pero la superficie interna de la esfera hueca se mantiene bajo una presión uniformemente distribuida P1 y se adiciona una distribución de temperatura constante θ1, mientras que la superficie externa del casco esférico se encuentra sometida a la misma presión P2 y a una temperatura θ2. La esfera interna  posee propiedades E1, υ, k, α1 mientras que las del casco externo, son E2, υ, k y α2. La geometría se modeló con dos regiones interconectadas con propiedades definidas, con un total de 128 elementos isoparamétricos con aproximación lineal. Para estudiar la convergencia de los valores reportados en puntos internos, se realizó un segundo análisis discretizando el modelo con 226 elementos. Los valores numéricos empleados son:

             

Las expresiones analíticas de este problema se presentan a continuación (Timoshenko et al., 1951):

      (28)

        (29)

   (30)

Las constantes fueron deducidas a partir de las condiciones de contorno impuestas.

Tabla 5. Resultados numéricos para el problema termoelástico de la esfera multi-región (128 elementos). u_RR y  σ_RR  son los desplazamientos y  tensiones respectivamente en dirección R, radial al centro de la esfera, mientras que  σ_tt es la tensión en dirección tangencial a las superficies de la esfera. θ_r es el campo de temperatura.

 

(*) Valor promedio en base a los resultados obtenidos en los nodos concernientes.

(°) Valor calculado en base a un sistema de coordenadas local al elemento.

Tabla 6. Resultados numéricos alternativos en puntos internospara el problema termoplástico de la esfera multi-región considerando ahora 226 elementos, haciendo énfasis en la superficie de interface.

7. DISCUSIÓN

Se comparan los desplazamientos en dirección radial, así como las tensiones en dirección radial, longitudinal y circunferencial y las campos de temperatura obtenidos de las expresiones analíticas tomadas de la literatura, con los resultados numéricos obtenidos mediante la aplicación de los nuevos códigos, para los puntos más representativos de los problemas analizados, simulando inicialmente un comportamiento puramente elástico y luego acoplando el campo de temperatura para simular el comportamiento termoelástico.

En los casos evaluados en el apartado 5 se presentan los resultados elásticos obtenidos. En todos los casos, las aproximaciones numéricas de las variables en estudio siguen el comportamiento impuesto por la solución analítica. En relación a los desplazamientos calculados en condiciones ideales de contacto e interferencia en el problema 5.2, los valores obtenidos no reportaron el grado de exactitud esperado, aunque fueron del mismo orden de magnitud. En elementos de contorno, el tratamiento de las interfaces amerita especial atención, al momento de imponer las relaciones de compatibilidad que permitirían incluir el efecto de una región sobre otra. Los autores decidieron respetar las relaciones del código original (Beer, 2001).

Por otra parte, el campo de desplazamientos en general del problema 5.2 muestra una rápida tendencia a disminuir en comparación con los valores reportados por el problema analizado en el ejemplo 5.1, cuya tendencia es menos pronunciada. La principal razón es el aumento en el espesor al considerar ahora dos regiones, lo cual opone mayor resistencia al movimiento en dirección radial. En la interface se presenta un cambio en el comportamiento de la curva desplazamiento (pendiente), debido a que las propiedades de los materiales involucrados fueron consideradas diferentes para este ejemplo.

En la figura 6, se evidencia una tendencia numérica hacia la curva de desplazamientos obtenida de la solución analítica. Esta aproximación mejora significativamente al incrementar el numero de elementos, sobretodo a nivel de la base, lo cual era de esperarse considerando que en el método de los elementos de contorno, al aumentar la cantidad de puntos de colocación (nodos) el error de aproximación tiende a disminuir (Beer, 2001). De igual manera ocurre con la evaluación del campo de desplazamientos en los puntos internos.

En cuanto al comportamiento del campo de tensiones, se presenta una tendencia similar a la reportada en el análisis anterior. Los resultados obtenidos fueron producto de un sencillo estudio de convergencia, con una desviación estándar mínima de 0.2865, lo que aseguró un grado de exactitud aceptable. El comportamiento mostrado por las otras variables también coincidió notablemente con el resultado analítico.

Considerando los resultados obtenidos en los problemas del apartado 5 y 6, se puede argumentar que en casos con múltiples regiones, el grado de discretización debe ser aún más fino, sobretodo a nivel de las interfaces. Sí tomamos el ejemplo 5.2 y acoplamos el fenómeno térmico (problema 6.3), los campos de temperatura y de tensiones en puntos internos, pierden precisión al ser calculados mediante el proceso de integración convencional a medida que éstos se acercan a las fronteras, aún cuando los resultados en el contorno fueron satisfactorios. Al realizar pruebas sucesivas con grados de discretización variables, la precisión en los cálculos mejoró considerablemente. Para ilustrar esta observación, se aumentó la cantidad de elementos utilizados en las superficies curvas de la esfera multi-región, haciendo énfasis en la interfase. El modelo inicial contó con 128 elementos, mientras que el más completo consistió de un total de 226 elementos, reportando mejoras en los valores de los campos ya mencionados.

8. CONCLUSIONES Y RECOMENDACIONES

Los problemas que implican los sólidos tridimensionales de simetría axial (o sólidos de revolución) sometidos a carga axial simétrica se reducen a  problemas bidimensionales, muy útiles a la hora simplificar el problema y disminuir los tiempos de cómputo.

El ya versátil Método de los Elementos de Contorno (BEM) también pudo incorporar esta ventaja, gracias a la forma general la ecuación integral de contorno axisimétrica deducida, entre otros, por Balas et al. (1989) y con correcciones y nuevas metodologías propuestas por Graciani (2006), mediante el uso de soluciones fundamentales cuasiestáticas definidas en el contorno que no requieren de integración en el dominio interno del problema. Las singularidades inherentes a los campos de desplazamientos y tensiones requirieron de un estudio asintótico más profundo dentro de la formulación, lo que condicionó la evaluación e integración de términos adicionales. Sin embargo, la incorporación de las nuevas soluciones fundamentales en un algoritmo clásico basado en la teoría de elasticidad para BEM (Beer, 2001) fue relativamente fácil, debido a la versatilidad presentada por la estructura de programación de dicho código.

 Los ejemplos clásicos seleccionados para la validar la implementación, permitieron verificar la eficiencia en la evaluación de la nueva formulación. El proceso de integración requirió especial cuidado debido a la presencia de las integrales elípticas completas de primer y segundo orden, cuyo comportamiento asintótico fue determinante para establecer el procedimiento de integración más adecuado en los casos singulares.

La incorporación a la ecuación integral de contorno de los términos de potencial y de acople, permitieron establecer con buena aproximación el campo termoelástico de tensiones y deformaciones para problemas que toman en cuenta la variación de temperatura debido al fenómeno de conducción de calor en estado estacionario, tal y como se demuestran en los ejemplos tratados. La presente formulación involucra solo variables de superficie mediante la implementación de soluciones fundamentales en medio infinito. Nuevamente, al igual que en el caso puramente elástico, se consideraron estados de carga y condiciones de contorno con simetría axial.

A pesar de que la evaluación e integración de la expansión singular de los núcleos, en aquellos casos cuando el radio de un punto de colocación es “cero” amerita de cierto cuidado, el hecho de contar con una cantidad de elementos importantes en la discretización, suele garantizar una buena aproximación en los resultados obtenidos, sobretodo, en las adyacencias al eje de simetría, donde en muchos casos se acostumbra a modificar la coordenada radial a conveniencia para disminuir la influencia de la singularidad geométrica presente, aún cuando formulaciones especiales sean implementadas para estos casos (expansiones singulares). Los resultados reportados demuestran una vez más, la eficacia de la formulación implementada.

9. AGRADECIMIENTOS

• Proyecto ALFA II-0357-FA-FCD-FI (ELBEnet).

• División de Soldadura y Fractomecánica, Instituto Tecnológico de Materiales (INTEMA), Universidad Nacional de Mar del Plata, Mar del Plata, Argentina.

• Instituto Nacional de Bioingeniería (INaBio), Universidad Central de Venezuela, Caracas, Venezuela.

• Grupo GEMM. Instituto de Investigación en Ingeniería de Aragón, Zaragoza, España.

• Fondo Nacional de Ciencia, Tecnología e Innovación (FONACIT), Caracas, Venezuela.

• Consejo de Desarrollo Científico y Humanístico (CDCH), Universidad Central de Venezuela, Caracas, Venezuela.

• Prof. Enrique Graciani, Escuela Superior de Ingenieros, Departamento de Mecánica de Medios Continuos, Universidad de Sevilla, España.

10. REFERENCIAS

1. Bailón-Plaza Alicia, Van der Meulen Marjolein C.H., “Beneficial effects of moderate, early loading and adverse effects of delayed or excessive loading on bone healing”, Journal of Biomechanics, 36, 1069-1077 (2003).        [ Links ]

2. Brebbia C.A., Dominguez J., Boundary element: an introductory course, 2th edition, Computational Mechanics Publications,U.K. (2003).        [ Links ]

3. Beer G., Programming the boundary Element Method: An introduction for Engineers, 1th edition, John Wiley & Sons, Ltd., U.K. (2001).        [ Links ]

4. Banerjee P.K., The Boundary Element Methods in Engineering, 2th edition, McGraw-Hill International, U.K. (1994).        [ Links ]

5. Brebbia C.A., Telles J.C.F., Wrobel L.C., Boundary Element Techniques: Theory and applications in Engineering, 1th edition, Springer-Verlag, Berlin Heidelberg New York Tokyo (1984).        [ Links ]

6. Balas J., Slâdek J., Slâdek V., Stress analysis by boundary element method, 1th edition, Amsterdam: Elsevier, (1989).        [ Links ]

7. Bakr A. A., Fenner R. T., “Boundary integral equation analysis of axisimmetric thermoelastic problems”, Journal of strain analysis, 18(4), 239-251 (1983).        [ Links ]

8. Chapetti Mirco D., Mecánica de Materiales, primera edición, Ediciones al Margen, La Plata, Buenos Aires (2005).        [ Links ]

9. Claes L. E., Heigele C. A., “Magnitudes of local stress and strain along bony surfaces predict the course and type of fracture healing”, Journal of Biomechanics, 32, 255-266 (1999).        [ Links ]

10. Cowin S.C., “Bone poroelasticity”, J. Biomech, 32, 218-238 (1999).        [ Links ]

11. Chopra M. B. and Dargush G.F., “Thermal stress analysis of axisymmetric bodies via the boundary element method”, Computer Methods in Applied Mechanics and Engineering, 108, 53-71 (1993).        [ Links ]

12. Cruse T.A., Snow D.W., Wilson R.B., “Numerical solutions in axisymmetric elasticity”, Compt. Strut., 7, 445-451  (1977).        [ Links ]

13. De Lacerda L.A., Wrobel L.C., “Frictional contact analysis of coated axisymmetric bodies using the boundary element method”, Journal of strain analysis, 35(5), 423-440 (2000).        [ Links ]

14. Duda G., Kirchner H., Wilke H., Claes L., “A method to determine the 3-D stiffness of fracture fixation devices and its application to predict inter-fragmentary movement”, Journal of Biomechanics, 31, 247-525 (1998).        [ Links ]

15. Dargush G. F., Banerjee P.K., “Time dependent axisymmetric thermoelastic boundary element analysis”, Int. J. Numer, Methods Eng., 33, 695-717 (1992).         [ Links ]

16. González Y., Cisilino A., Cerrolaza M.. “Aplicación del método de elementos de contorno a problemas con simetría axial: implementación de la formulación para elasticidad lineal”. XV Congreso sobre Métodos Numéricos y sus Aplicaciones ENIEF 2006, Santa Fe, Argentina (2006).        [ Links ]

17. Gómez-Benito M. J., García-Aznar J. M., Kuiper J. H., Doblaré M., “Influence of fracture gap size on the pattern of long bone healing: a computational study”, Journal of Theoretical Biology, 235, 105-119 (2005).        [ Links ]

18. Graciani E., Mantic V., Paris F., Blazquez A., “Weak formulation of axisymmetric frictionless contact problems with boundary elements Application to interface cracks”, J. Comp. Estruc., 83, 836-855 (2005).        [ Links ]

19. García J. M., Kuiper J. H., Doblaré M., Richardson J. B., Gómez M. J., “Simulación numérica del proceso de reparación de fracturas óseas”, V Congreso sobre Métodos Numéricos en Ingeniería (SEMNI), España, (2002).        [ Links ]

20. Rizzo F.J. and Shippy D.J., “An advanced boundary integral equation method for three-dimensional thermoelasticity”, Int. J. Numer. Methods Eng., 11, 1753-1768 (1977).        [ Links ]

21. Timoshenko S., Goodier J. N., Theory of elasticity, Second edition, McGrawh-hill book company Tokyo (1951).        [ Links ]