Análisis de problemas elásticos 2d utilizando la técnica de descomposición de dominio paralelo-iterativa y el método de los elementos de contorno

Brizeida Gámez ^1,4, Eduardo Divo ², Alain Kassab ³, Miguel Cerrolaza ⁴

¹ Departamento de Diseño Mecánico y Automatización, Facultad de Ingeniería, Universidad de Carabobo, Naguanagua, Edo. Carabobo, Venezuela, Código Postal 2005, Teléfono: +58-416-3312551, Fax: 0241-8672814, brizeida.gamez@gmail.com

² Department of Engineering Technology, University of Central Florida, Orlando, FL, USA

³ Department of Mechanical, Materials and Aerospace Engineering, University of Central Florida, Orlando, FL, USA

⁴ Instituto Nacional de Bioingeniería, Universidad Central de Venezuela, Caracas, Venezuela

RESUMEN

Para resolver problemas de elasticidad con el Método de Elementos de Contorno (MEC) sólo se requiere discretizar la superficie; sin embargo las matrices resultantes se caracterizan por no ser diagonalmente dominantes y estar completamente llenas. Esto representa un reto para resolver problemas de gran escala, dado los requerimientos de almacenamiento de datos y la solución de sistemas de gran número de ecuaciones. En este artículo se desarrolla una técnica de descomposición de dominio o multiregiones eficiente y efectiva, con algoritmos iterativos región por región, ensamblados especialmente para paralelismo computacional. La aproximación de descomposición de dominio efectivamente reduce el número de condiciones de los sistemas de ecuaciones algebraicas resultantes, al mismo tiempo que incrementa la eficiencia de los procesos de solución, convergiendo rápidamente a una solución estable. El paralelismo computacional satisface perfectamente la técnica iterativa de la descomposición de dominio. Los resultados demuestran la calidad de la aproximación al comparar con las soluciones de problemas en una región y las respectivas soluciones analíticas.

Palabras clave: descomposición de dominio, elasticidad, método de elementos de contorno, paralelismo computacional.

Analysis of 2d elastic problems by using a parallel-iterative domain decomposition technique and the boundary element method

ABSTRACT

The Boundary Element Method (BEM) requires only a surface mesh to solve elasticity problems; however, the resulting matrices are fully-populated and non-diagonally dominant. This poses serious challenges for large-scale problems due to storage requirements and the solution of large sets of non-symmetric systems of equations. In this article, an effective and efficient domain decomposition, or artificial sub-sectioning technique, along with a region-by-region iteration algorithm particularly tailored for parallel computation to address these issues is developed. The domain decomposition approach effectively reduces the condition numbers of the resulting algebraic systems, while increasing efficiency of the solution process and decreasing memory requirements. The iterative process converges very efficiently while offering substantial savings in memory. The iterative domain decomposition technique is ideally suited for parallel computation. Results demonstrate the validity of the approach by providing solutions that compare closely to single-region BEM solutions and benchmark analytical solutions.

Keywords: domain decomposition, elasticity, boundary element method, parallel computation.

Recibido: diciembre de 2006 Revisado: abril de 2007

INTRODUCCIÓN

El Método de Elementos de Contorno (MEC) ha demostrado ser una técnica numérica muy poderosa en muchas disciplinas de ingeniería. El atractivo de este método se atribuye a la reducción de la dimensionalidad del problema, ya que para análisis en dos dimensiones sólo necesita ser discretizada la línea del contorno del dominio y para problemas en tres dimensiones únicamente se requiere la discretización de la superficie del dominio. Esto significa que comparándolo con métodos de análisis de dominio, el mismo garantiza una sustancial reducción en la preparación de datos y en los sistemas de ecuaciones que deben ser resueltos (Xiangqiao et al., 2006); sin embargo las matrices resultantes se caracterizan por no ser diagonalmente dominantes y estar completamente llenas lo cual representa un reto para resolver problemas de gran escala, dado losrequerimientos de almacenamiento de datos y la solución de un gran número de sistemas de ecuaciones.

Dado lo anterior, en el presente artículo se desarrolla una técnica de descomposición de dominio o multiregiones eficiente y efectiva, con algoritmos iterativos región por región, particularmente ensamblados por paralelismo computacional. La aproximación de descomposición dedominio efectivamente reduce el número de condiciones de los sistemas de ecuaciones algebraicos resultantes mientras incrementa la eficiencia de los procesos de solución y disminuye los requerimientos de memoria. El procesoiterativo converge eficientemente para todas las situaciones reportadas a la vez que proporciona ahorros significativos en la memoria del computador para soluciones eficientes de problemas de elasticidad en 2D. El paralelismo computacional satisface perfectamente la técnica iterativa de la descomposición de dominio. Se obtienen resultados para geometrías con diferentes tipos de condiciones de contorno impuestas. Los resultados demuestran la validez de la aproximación comparando con las soluciones de problemas en una región y las respectivas soluciones analíticas.

MÉTODO DE ELEMENTOS DE CONTORNO (MEC) EN ELASTICIDAD

Con el MEC es posible determinar las tracciones, desplazamientos y esfuerzos en el contorno y en los puntos internos de un dominio (Brebbia et al., 1989 y Cheng et al., 2001). Para establecer las ecuaciones integrales de contorno en elasticidad, es necesario indicar previamente la ecuación de equilibrio (1) y la constitutiva (2) que rigen dicho fenómeno, además de las relaciones de deformación (3) y las tracciones en el contorno (4), las cuales pueden ser expresadas como:

σ_ij representa el campo de esfuerzos, b_i son las fuerzas decampo, δ_ij es el delta de Kronecker, e_ij es la deformación, ui es el desplazamiento, t_i es la tracción y n_j es el vector normal unitario apuntando hacia afuera del contorno Γ del dominio Ω. Las propiedades μ y õ son el módulo de rigidez y la relación de Poisson, respectivamente. Combinando las ecuaciones (1) y (2) y sustituyendo la deformación de la ecuación (3), se obtiene la ecuación de Navier, la cual viene dada por:

Un planteamiento adecuado de un problema con valor en el contorno, requiere que en cada porción de la frontera sea prescrita cualquiera de las dos condiciones, tracción ti = ti‾ sobre el contorno Γ_t o desplazamiento u_i = u_i sobre el contorno Γ_u ; donde Γ_t U Γ_u = Γ es el contorno de solución del dominio Ω [Cheng et al., 2001]. Finalmente en una formulación del método de elementos de contorno basada en la ecuación integral de Somigliana (Wei et al., 2002 y Gaul et al., 2003), queda planteada una relación integral entre los desplazamientos en un punto de colocación «p» los desplazamientos u_i y las tracciones t_i en todo el contorno Γ . Adicionalmente, las fuerzas de campo b_i, permanecen relacionadas a través de una integral de dominio , como sigue:

donde:

H_ij y G_ij son las soluciones fundamentales en términos de tracción y desplazamiento, c^p_ij es una constante geométrica, tal que c^p_ij = 1 si p ∈ Ω , c^p_ij= 0 si p ∉ Ω y c^p_ij= ½ en el caso de un contorno suave (Brebbia et al., 1989). En el caso de dos dimensiones las soluciones fundamentales son: 

la distancia r está definida como p r = │x_i − x^p_i│

DESCOMPOSICIÓN DE DOMINIO ITERATIVA (DDI)

Para determinar la eficiencia de la técnica DDI es necesario, en principio, indicar los requerimientos de memoria para la resolución de un problema de dominio Ω en una región simple (figura 1), con la discretización característica del método de elementos de contorno.

Figura 1. Esquema de la discretización del contorno de una región simple.

Si se genera una discretización de un Número de Elementos (NE) con un Número de Nodos (NN) independientes por cada elemento, en una región simple, el sistema algebraico resultante tiene dimensiones N, donde N = d × NE × NN y d es el número de dimensiones espaciales (2 ó 3). De forma tal que:

donde:

el vector {x} representa los valores desconocidos de las tracciones y de los desplazamientos en el contorno. Para este caso el número de operaciones de puntos flotantes requeridos para generar el sistema algebraico es proporcional a N² , asimismo la asignación en la memoria es proporcional a N² .

La solución del sistema algebraico puede llevarse a cabo usando un método de solución directa tal como la descomposición Lower-Upper (LU) requiriendo operaciones de puntos flotantes proporcionales a N³ o usando un método indirecto tal como el de Gradiente Bi-Conjugado (BCG) o Residuos Mínimos Generalizados (GMRES), los cuales en general requieren operaciones de puntos flotantes proporcionales a N² para alcanzar la convergencia (Divo et al., 2005).

Por otra parte, si ahora se adopta una aproximación basada en una descomposición de dominio, el dominio original Ω es dividido en K sub-dominios Ω _l ∴ l = 1....k separados por interfases creadas artificialmente, entonces cada subdominio puede ser discretizado con elementos de contorno de manera independiente (figura 2), para un caso donde K=4. Sucesivamente, la solución en cada sub-dominio puede ser obtenida de forma independiente a través de un proceso estándar, siempre y cuando las condiciones de borde en las interfases artificiales entre los sub-dominios sean conocidas. Por ejemplo, el primer sub-dominio Ω₁ en la figura 2 es analizado independientemente. La aplicación del método de elementos de contorno en este sub-dominio genera un sistema algebraico de la siguiente forma:

Figura 2. Esquema de una descomposición de dominio y discretización típica.

donde la dimensión n=dxnexNN es obviamente una fracción de la dimensión N del sistema algebraico original. Si la resolución de la discretización en las interfases es similar a la del resto del contorno, la nueva dimensión n es aproximadamente equivalente a:

La composición de este sistema algebraico requiere operaciones de puntos flotantes proporcionales a n² . La asignación de memoria es también proporcional a n² lo que significa que el ahorro en memoria (Divo et al., 2003) es proporcional a:

En el caso de K=4 esto equivale al 16,6% de la memoria necesaria para una región simple. Los requerimientos de asignación de memoria directa para posteriores manipulaciones algebraicas descienden a una proporción n² mientras que los sistemas de matrices pueden ser almacenados fácilmente en disco para su uso posterior, luego que las soluciones en los sub-dominios restantes hayan sido resueltas efectivamente. En el caso de la solución de los sistemas algebraicos utilizando métodos directos, la proporcionalidad en operaciones de puntos flotantes es n³, lo que conduce a un ahorro de:

En el caso de K=4 la solución directa de cada uno de los sistemas de los sub-dominios consumiría un 6,4% del tiempo que tomaría la solución del problema de región simple. Evidentemente este ahorro no es directo ya que deben efectuarse un total de K soluciones. Sin embargo, este proceso es perfectamente implementable en paralelo y el ahorro puede ser alcanzado en su totalidad.

Naturalmente, las condiciones de borde en las interfases artificiales entre los sub-dominios son originalmente desconocidas así que debe concebirse un esquema que garantice la continuidad y equilibrio de las variables de campo entre los sub-dominios, esto implica que: 

donde:

los superíndices a y b indican cada lado de la interfase en cuestión. Para poder garantizar estas condiciones, es preciso en principio partir de una suposición de estos valores que sea suficientemente adecuada. Para esto se lleva a cabo un pre-proceso en donde estos valores se aproximan utilizando una discretización preliminar de baja resolución en una región simple. En consecuencia, un proceso iterativo es empleado para progresar desde la suposición inicial hasta una solución final que satisfaga una norma o nivel de convergencia adaptada a las condiciones de continuidad y equilibrio.

Particularmente, la progresión del proceso iterativo que se plantea en este articulo involucra dos etapas: en la primera etapa, cada interfase es impuesta individualmente con condiciones de primer tipo (prescribiendo desplazamientos u_i ). Si se prescriben los desplazamientos u_i en las interfases, se obtienen tracciones ti como consecuencia de la aplicación del método de elementos de contorno estándar en cada sub-dominio. Estas tracciones obtenidas no coinciden en ambos lados de cada interfase ya que la suposición original no es generalmente correcta. De forma que es necesario alterar estas tracciones para obligarlas a satisfacer las condiciones de equilibrio utilizando las siguientes relaciones:

De esta manera se garantiza que las tracciones actualizadas a ambos lados de la interfase tengan la misma magnitud pero signos opuestos satisfaciendo la condición fundamental de equilibrio. Una vez que estas tracciones sean actualizadas, se prosigue a la segunda etapa del proceso iterativo en donde las interfases son impuestas con condiciones de segundo tipo (prescribiendo tracciones t_i). Prescribiendo las tracciones obtenidas en la actualización anterior y prosiguiendo con soluciones estándar por el método de elementos de contorno en cada sub-dominio, un nuevo campo de desplazamientos es obtenido en las interfases. De nuevo, estos desplazamientos no coinciden en ambos lados de cada interfase y es necesario actualizarlas para satisfacer las condiciones de continuidad tal que:

Garantizando igual magnitud y signo en los desplazamientos a ambos lados de cada interfase. En este punto, puede ser formulada una norma que mida la distancia entre el desplazamiento actualizado y el anterior para condicionar elproceso iterativo a detenerse cuando satisfaga una norma mínima impuesta. Esto es:

donde: NNI es el número de nodos en las interfases. El proceso iterativo debe detenerse en el caso que la norma esea menor a una norma de convergencia impuesta δ, de locontrario, el proceso iterativo debe continuar desde laprimera etapa utilizando los desplazamientos actualizados como condiciones de borde.

Adicionalmente, si se emplea un método de solución directa tal como la factorización LU para todos los sistemas algebraicos en cada uno de los sub-dominios, los factoresLU de las matrices de coeficientes para todos los subdominios se pueden calcular en el primer paso y almacenaren la memoria del disco para su uso posterior en el proceso iterativo. De esta manera, para llegar a la solución de cada sistema algebraico en cada paso del proceso iterativo serequeriría de un proceso de substitución simple y eficiente utilizando factores LU calculados previamente. Estacaracterística permite una reducción significativa en la cantidad de operaciones de punto flotante así como del tiempo computacional.

IMPLEMENTACIÓN PARALELA

El algoritmo propuesto ha sido codificado e implementado en un computador con procesador Intel Centrino Duo,velocidad de 1.66GHz y 1GB de memoria RAM, corriendo bajo ambiente Windows XP y compilado con Visual FORTRAN.

Con relación a la puesta en operación del código bajo MPI (Message Passing Interface), el modelo descompuesto en sub-dominios es asignado de acuerdo al número de regionesa cada procesador, de manera que cada uno de ellos correun componente (sub-dominio) del programa codificado enparalelo; esto con el fin de incrementar la eficiencia de losprocesos de solución y la utilización de la memoria del computador (Erhart et al., 2006.

Todo lo anterior demuestra que la formulación del métodode elementos de contorno para descomposición de dominio satisface perfectamente el paralelismo computacional.

El proceso anterior esta basado en las experiencias mostradas detalladamente en la referencia (Divo et al., 2005), en la cual se ha resuelto una pequeña muestra de un problema del MEC sobre todos los procesadores para identificar sucomportamiento relativo, seguido de una rutina de balancede carga para asignar dominios de manera óptima a cada procesador a través de minimización de una función objetivo.

VALIDACIÓN NUMÉRICA Y RESULTADOS

En esta sección se muestra la aplicación de la técnica de la descomposición de dominio iterativa en paralelo usando elmétodo de elementos de contorno, basada en dos ejemplos clásicos de la mecánica de sólidos con característicasbidimensionales. Para todos los modelos se obtienenresultados para una, dos y cuatro regiones, los cuales severifican con la solución exacta o con resultados producidospor otros autores. En todos los ejemplos, las propiedadesmecánicas empleadas son: ì=60GPa y õ=0,3. Es importante destacar que el tiempo total reportado para lograr los resultados comprende la resolución del campo elástico apartir de la obtención de los valores de los esfuerzos y desplazamientos en los puntos internos y en el contorno.

Viga en voladizo

Se considera una viga empotrada-libre, de longitud L=0.25m y altura h=0.05 m, sometida a una carga distribuida (P) de 1x10⁶Pa actuando sobre toda la longitud, como se indica en la figura 3. El modelo se discretiza con una malla compuesta por cien elementos cuadráticos. En las figuras 4 y 5 se distinguen los resultados de las líneas de deformaciones y esfuerzos de Von Mises, respectivamente. El modelo converge rápidamente como puede apreciarse en la figura 6, la cual muestra la relación entre el número de iteraciones y la norma, en una, dos y cuatro regiones, para la obtención de resultados debido a las condiciones elásticas. Por otra parte se indica el tiempo requerido para solucionar el modelo (figura 7), en donde es notable la eficiente disminución en el tiempo de resolución del problema, para múltiples regionescomparado con el modelo resuelto en una sola región, ya que se alcanza una reducción de 54,35% para el modelo de dos regiones y 57,60% para cuatro regiones. Finalmente cabe destacar que la diferencia entre los valores de desplazamientos obtenidos (u_DDI) para el centro del extremo libre de la viga, y la solución analítica (u_anal) presenta un error de 3,17%, como se muestra en la tabla 1.

Figura 3. Condiciones de contorno para viga en voladizo.

Figura 4. Diagrama de las líneas de desplazamientos para (a) una región, (b) dos regiones, (c) cuatro regiones.

Figura 5. Diagrama de las líneas de Esfuerzos de Von Mises para (a) una región, (b) dos regiones, (c) cuatro regiones.

Figura 6. Diagrama de convergencia para el modelo de viga en voladizo.

Figura 7. Tiempo requerido para la resolución del modelo viga en voladizo.

Tabla 1. Solución analítica vs. solución DDI

Viga a tracción

En la figura 8 se muestra una viga empotrada en un extremo y libre en el otro, de dimensiones L=0.25m. y h=0.005 m,sometida a una carga de tracción (P) de 1x10⁶ Pa, aplicada en la superficie libre de la viga, la misma se discretiza con cien elementos cuadráticos.

Figura 8. Condiciones de contorno para viga a tracción.

Las figuras 9 y 10 muestran las líneas de deformaciones y esfuerzos de Von Mises para el ejemplo considerado en una, dos y cuatro regiones. Al igual que en el caso anterior esnotable la eficiente convergencia del modelo, como puedeapreciarse en la figura 11. Adicionalmente se indica el tiemporequerido para obtener los resultados de la viga analizadapara una, dos y cuatro regiones (figura 12), en la misma se reporta una reducción del tiempo, en comparación con el modelo resuelto en una sola región, de 56,52% y 59,78% para dos y cuatro regiones, respectivamente.

Figura 9. Diagrama de las líneas de desplazamientos para (a) una región, (b) dos regiones, (c) cuatro regiones.

Figura 10. Diagrama de las líneas de Esfuerzos de Von Mises para (a) una región, (b) dos regiones, (c) cuatro regiones.

Figura 11. Diagrama de convergencia para el modelo de viga a tracción.

Figura 12. Tiempo requerido para la resolución del modelo viga a tracción.

CONCLUSIONES

Se ha presentado la aplicación de la técnica de descomposición de dominio paralela iterativa para la evaluación de problemas elásticos. La misma fue implementada empleando un programa basado en la formulación para problemas de elasticidad usando elMétodo de Elementos de Contorno (MEC), el cual es capazde calcular esfuerzos y desplazamientos en el contorno y en los puntos internos del dominio, utilizando las correspondientes ecuaciones integrales.

A partir de los análisis de una serie de ejemplos, con diferentesgeometrías y condiciones de contorno impuestas sobre superficies regulares, se demuestra la eficiencia de la técnica propuesta, lo que explica que el paralelismo computacional satisface perfectamente la técnica iterativa de la descomposición de dominio.

El proceso iterativo converge eficientemente para todas las situaciones reportadas a la vez que proporciona ahorros significativos en la memoria del computador y del tiempo computacional, lo cual conduce a soluciones eficientes de problemas de elasticidad en 2D.

AGRADECIMIENTOS

Los autores desean manifestar su agradecimiento a la Universidad de Carabobo, University of Central Florida-USA y al Fondo Nacional de Ciencia y Tecnología-Venezuela,por el apoyo académico y financiero otorgado para la realización de este trabajo.

REFERENCIAS

1. BREBBIA, C., DOMÍNGUEZ, J., (1989): «Boundary Element: An Introductory Course». Computational Mechanics Publ. Boston copublisher McGraw Hill. New York, p. 314.        [ Links ]

2. CHENG, A.H., CHEN, C.S., GOLBERG, M.A., RASHED, Y.F., (2001): «BEM for thermoelasticity and elasticity with body forcea revisit». Engineering Analysis with Boundary Elements, pp. 25, 377-387.        [ Links ]

3. DÍAZ, F., YATES, J.R., PATTERSON E.A., (2004): «Some improvements in the analysis of fatigue cracks using thermoelasticity». International Journal of Fatigue, pp. 26, 365-376.        [ Links ]

4. DIVO, E., KASSAB, A., (2005): «A meshless method for conjugate heat transfer problems». Engineering Analysis with Boundary Elements, 29, pp. 136-149.        [ Links ]

5. DIVO, E., KASSAB, A., RODRÍGUEZ, F., (2003): «Parallel domain decomposition approach for large-scale threedimensional boundary-element models in linear and nonlinear heat conduction». Numerical Heat Transfer, Part B. 44:417-437.        [ Links ]

6. ERHART, K., DIVO, E., KASSAB, A., (2006): «A parallel domain decomposition boundary element method approach for the solution of large scale transient heat conduction problems». Engineering Analysis with Boundary Elements 30, 553-563.        [ Links ]

7. GAUL, L., KÖEL, M., WAGNER, M., (2003): «Boundary Element Methods for Engineers and Scientists: An introductiory Course with Advanced Topics». Springer-Verlag. Berling-Heilderberg, p. 488.        [ Links ]

8. WEI GAO, X., DAVIES, T., (2002): «Boundary Element Programming in Mechanics». Cambridge University Press, United Kingdom, p. 249.        [ Links ]

9. XIANGQIAO, Y., (2006): «A boundary element modeling of fatigue crack growth in a plane elastic plate». Mechanics Research Communications 33: 470–481.        [ Links ]