Tratamiento de las integrales de dominio en el análisis de la ecuación de poisson por el método de los elementos de contorno

ANALES DE INGENIERIA MECANICA. AÑO l. N° 1, 1983 TRATAMIENTO DE LAS INTEGRALES DE DOM INI O EN EL ANALIS IS DE LA ECUACION DE POISSON POR EL METODO DE LOS ELEMENTOS DE CONTORNO Federico París 1 Segundino de Le ón Escuela T écnica Superior de Ingenieros Industriales de Sevilla Resumen.- El Método de los elementos de contorno basa su princ i pal ventaja, como alternativa numérica de cálculo, en que basta discretizar el contorno del dominio en estudio. Cuando la naturaleza del problema anal izado incluye una sol ic i tación de dominio, como es el caso de la ecuación de Poisson, esta ventaja desapa rece. Para ciertas expresiones particulares de la s~ 1 icitación de dominio, es posible su transformación al contorno mediante la aplicación del - teorema de la divergencia. Para una so li citación general, se ensaya en este trabajo, una --- aproximación numérica mediante integrales de contorno, que son idénticas a las que aparecen en la ecuación original, por lo que la aproximación de la inte g ral de dominio no representa un excesivo cálculo adiciona l. l. INTRODUCCION. Existen en la actualidad diferentes alternativas numéricas para el tratamiento aproximadode la mayor parte de los prob lemas de la FísicaMatemática, aunque ninguno de ellos presente ven tajas universales. El grupo de soluciones de do~ minio, fundamentalmente el Método de los Elementos Finitos, basados en una formulación débil -- del problema en estudio conducen a un sistema de ecuac i ones simétrico y para el caso del M.E.F. - en banda, lo que, pese al usualmente gran tamaño del s istema de ecuaciones, simplifica notablemen te la res o lución. El grupo de las alternativas~ de contorno, basada en un teorema de reciproci-- dad, conduce a un sistema de ecuaciones 1 leno y no simétrico pero de tamaño más reducido. Exis-- ten un gran número de publicaciones en que se - - compara la efic i encia computac i onal de ambos métodos de cálculo, 111, 121, 131 pueden serv ir de ejemplo. Sin embargo, independientemente de las ventajas de una u otra alternativa en l os aspectos numéricos como tamaño y tiempos de genera--- ción y resolución del sistema de ecuaciones que dependen fundamentalmente del ejemplo concreto - sobre el que se establezca la comparación, hay - una diferencia objetiva a favor de la alternativa de contorno, que cada día tiene un peso mayor en la explotación comercial de los programas que recogen ambas posibilidades. El tiempo de preparación de datos se disminuye drásticamente si só lo es preciso discretizar el contorno del domi-~ nio en estudio. Sin embargo, existen casos en -- que esta ventaja desaparece debido a la existencia de soli c itaciones de dominio, fuerzas en volumen en problemas elastostáticos, carga trans-- 26 versal en placas o la presencia de un segundo - - miembro en la ecuación de Pois son . En estos ca-- sos es preciso evaluar una integral de dominio, lo que requiere la discretización del mismo. El hech o de que esta discretización no necesite especiales requerimientos de finura debido a que - no incluye ninguna incógnita del problema y notiene por tanto trascendencia en el sistema deecuaciones no quita para que una de las princip~ les ventajas de la alte rnati va de contorno haya desaparecido. La transformación de la integral - de dominio en integral de contorno mediante la - aplica c ión del teorema de la divergencia sólo es posible para expresiones muy particulares de la sol icitació n y desde luego en ningún caso cuando la información de ésta es en forma discreta. La alternativa que se real iza en es te trabajo está en la 1 ínea de las propuestas real izadas por Nardini y Brebbia 141 y Garcia Suarez ISI · - El primero de ellos la aplica a problemas dinám~ cos con par ticularidades que en el desarrollo-- del presente artículo se pondrán de manifiesto Y el segundo utiliza funciones de la misma natura· leza pero formalmente diferentes para introducir la aproximación. 2. PLANTEAMIENTO DEL PROBLEMA. Sobre un dominio genérico Vcon un contorno él V, está definida una función <P gobernada por la ecuación de Poisson (1). él V1 U él V2 él V en él V1 <P '$ en él V2 q él<P añ= Fig. 1. Configuración del dominio bajo estudio. = -p ( 1 ) q Sobre dicho contorno se conoce bien el po-- tencial <Po el flujo q. La segunda fórmula de Green establece una - relación integral entre ~ y una f -~ 1 iar 1/J '+' unc 1on auxi--- 2 2 r¡ 1/J - 1/J 'V <P)dv ~~ - 1/J ~)ds (2) Si elegimos 1/J como la soluc ión fundamental de - la ecuación de L aplace, (2) queda en la forma <P(x) + J<P(y) ~ ~ (x ,y ) ds(y).= él V = J ljJ (x,y) ~~ (y) ds(y) + él V Jp(z) V ljl (x, z) . dA(z) (3) con x,z € V ; y € él V. Las incognita s del problema en (3) son ---- <P(x)' ~ (y) en él V1 Y <P(y) en él V2. Para obten er incognitas exclusivamente 1 en e conto rn o hace-- mos tender x a ~ v. En tal o caso obt enem os c(x) <P(x) + J<P(x) ~~ (x,y) ds(y) él V donde: = Jw(x,y) ~~ (y) ds(y) + av + Jp(z) 1/J(x,z) dA(z) V z € V. (4) x,y € a v; c(x) es una constante que depende de la forma real del contorno en x. La ecu · ~ (4) do de 1 ac,on constituye la base del mé to incluyeo~ e!em~ntos de contorno puesto que sól~ av él<P ncognJtas sobre el contorno, <P(y) en - _2 y él x (y) en éJVl. Una vez calculadas las incognitas mediante la resolución numérica de 27 (4), que se describe en el apartado se puede util i zar (3) para calcular puntos del dominio que se desee. siguiente, <P en los -- 3. FORMULACION NUMERICA. APROXIMACION DE LA INTEGRAL DE DOMINIO. La resolución de (4) implica una aproxima-- ción de <P(y) Y ~~ (y) en él V. As imi smo implica - una evaluación de la integral de dominio lo __ ~ue const!tuye el núcleo central de este' trabaJO. En P:~mer l~gar se describe brevemente la - formulac1on numerica de las integrales de con-- torno, es decir, e l caso en que al ser p = o _ la ec~~ción de campo (1), se transforma en 1~ ecuac1on de Laplace. Esta formulac ión pu ede encontrarse en detalle en las referencias 13 1 y - l ~l· Poster!ormente en 3.1 se anal iza la evalua c~o~ de la 1nte gra l de volumen para un valor ge n~rlco ~e ~ Y en 3.2 para el caso particular de discontinuidades en p. De cual quier forma tal como pu ede desprenderse de (4) 1 • t' y 1 de d · - , a 1n egra O~ lnl o, en el problema que estamos anal izan do, solo representa una constante a añadir - da · ~ d 1 • a ca ecuac1on e Sistema que se generará a contT nuación. Se va a establecer una - ~ él <P aprox1macion sobre <P y añ en av en la for m a: (S) donde Ni. son func ione s de pegueño.soporte 1 inea les asociadas al punto i y ~ ~ y q1 r - 1 '+' epresentan va ores de potenc ial y flujo en el punto i d 1 conto rno. Dicho contorno se aproxima también e __ por el~ment~s rectilíneos, representandose unpunto 1nter1or de un elemento cualquiera por: (6) donde: xj(c) d d ~ coo: ena_a cartesiana j de ·un punto Interior definido por la coorde denada natural ~ - - Na( ~ ) función de forma asociada al extre mo a del elemento. x~ coordenada cartesiana j del extremo a . La Fi g. 2 ilustra estas aprox 1 "mac· 1ones. La ecuación (4), sin tener en c uenta la · _ tegr~l ~e dominio, se puede escribir en la f~~ ma s1gu1ente, teniendo en cuenta la sustitu~ión del contorno real en uno discreto de N elemen-- tos. c(.Q.) <P(2) + N ¿ k=l (7) Si introducimos la aproximación (S) sobre cada - - elemento, (7) puede representarse en 1 a forma: (a) Definición de un (b) funciones a propunto del contorximantes no aproximado. F i g. 2. N t'(k) } c( 9, ) <j¡(Q.) + l: (A 1 A2) k=l ~<Pk (k+1) N r,(k) l í: (B1 B2) qk (k+1) (8) k=1 donde los subíndices de l as variables <P y q están asociados al elemento sob re el que se integra y los paréntesis a l nodo al que están aso-- ciados. A. y B. representan las integraciones - 1 1 sobre l os elementos una vez introducida (5). - - Así, A = f N et(f,;) d\jl dsk Ct a-n av k (9) B = f N \ji dsk Ct Ct av Si la ecuac1on (9) se escribe para todos los n~ dos del contorno , obtenemos un sistema de ecuaciones en la forma H <j¡ _ G q (NxN) (Nxl)-(NxN) (Nxl) (1 O) donde N es el número de nodos del contorno. Reordenando datos e incógnitas (10) se transforma en: K • X (NxN) (Nx l) F (:lxll ( 11 ) La integral de dominio que vamos a evaluar a -- co ntinu ac •on se rá preciso in cl ui rla en F, para cada nodo ( 9. ) del contorno. 28 3.1. Caso general. Aproximamos la distribución p(z) por una s~ rie de funciones f(z) definidas sobre todo el - dominio, en la forma: p(z ) M ¿ Ct . f.(J,z) J J M ¡: a . f.(r(J,z)) J J j=l j=l ( 1 2) donde f. (J,z) representa el valor de la función J asociada a un punto J del dominio o del contorno en un punto z. a. representa el coeficiente ' J asociado a la función anterior. En definitiva, f. juega el papel de un sistema de coordenadas J general izadas de dimensión M en el que pretend~ mos representar p(z) y Ct. son los componentes - J sobre dicho sistema que se ajustarán de acuerdo a algún criterio. Puesto que p(z) está definido en todo el dominio y contorno, los puntos J pu~ den elegirse sobre uno u otro. Por razones que aparecerán evidentes, más adelante se incluirán siempre los N nodos del contorno y los M-N puntos del dominio que se consideren necesarios p~ ra conseguir una buena aproximación de p(z) en V. Con la aproximac1on (12) la integral de dominio que aparec e en (4) se ha transformado en M inte grales de la forma: a . fr . (J,z) \j! (x,z) dA(z) J J V ( 13 ) Para que cada una de estas inte grales pueda ser expresada a partir de inte grales de contorno, - basta s i las asimilamos a la inte gral de domi - nio d~ (4) que f. represente a p, lo que equiv~ J . le a decir q ue debe existir una función ~ - asociada a cada f. que cump la la relación J J f. ( 14) J Puesto que no hay definidas condiciones de co nt orno sobre ~- o sus derivadas, basta con -- J que ~ . sea una solución particular de (14) , que J puede encontrarse para cada expresión de f. ( r (J, z)). J Asumiendo catorce y tomando como referencia (4), cada una de las integrales (13) puede ser expresada a part ir de integrales de contorno - en la forma: -J f j (J ,z) ljl (x ,z) dA(z) V C(x) ~ .(x) + J f f d~ .(J,y) + ~ j ( J , y) ~ (x, y) d s (y) - ----"'J..,.<l_nav av l)l(x ,y) ds(y) (15) Int rod uc iend o (12) y (15) en (4): e( X) <j¡ (X) + f <j¡ (y) ~~ (X , y) d S (y) av = J l)l (x,y) ~~(y) ds(y) - aj /c(x) ct>j(x) + av + Jct>j(J,y) ~~ (x,y) ds(y) -J av av . l)l (x,y) ds (y) 1 <lct> .(J,y) J ( 16) (16) incluye exclusivamente integrales de con-- torno de la misma naturaleza y cuya discretización ya se ha anal izado posteriormente. Puesto que los núcleos de la integración son los mis-- mos (integraciones de \ji y~~), (16) puede escri birse introduciendo la aproximación (5) para -- <lct> . ct> j y () ~ en la siguiente forma: H <P G q 1 H cp (NxN) (Nxl) (NxN) (Nx1 l (NxN) (NxM) q cp 1 G a ( 1 7) (NxN) (NxM) (Mx1l donde N es el número de puntos del contorno uti 1 izad o en la dlscretización del mismo y M el . nT de puntos utilizados para la aproximación de p y que incluye los N del contorno y M-N de dominio. La formación del sistema de ecuaciones para obtener (11) y su resolución no tienen ninguna particularidad, con respecto al caso de la ecua ción de Laplace. Asimismo la evaluación de laintegral para el cálculo de <P en los puntos inte rn os se real Iza en la misma forma que para el contorno. El va l or de ~ se determina ajustando p y a . fj mediante colocación en los M puntos. J p (Mxll F a (Mx M) 11"X1 ( 18) De (18) podemos obtener el valor de a para in-- e 1 u i r en ( 1 7) . 3.2. Caso de función discontinua. Vamos a anal izar el caso que se presenta -- con frecue · d . . nc1a e que p sea constante en subdomlnlos, pero de diferente valor. En la Flg. 3 - se presenta esta situac ión. Es preciso buscar funciones ~que verifi--- quen: ( 19) La integral de d om •'n¡ 'o d se escompone en: 29 av 1 av e p P¡ en V1 p p2 en v2 P¡ f. p2 Fig. 3. Función p discontinua. J p(z) l)l (x,z) dA(z) V JP ¡ (z) l)l (x,z) dA(z) + v, + Jp 2( z) l)l (x,z) dA(z) v2 (20) Y teniendo en cuenta (19) las integraciones sobre los subdominios v1 y v2 se pueden poner en notación abreviada, al igual que se hizo en ( 15) : f Pa \ji dA V a - le ci> + f ~ ~\ji ds a a a on av +av acp a e f \ji <l na dsj (21) av +av a e (21) exige la int egrac ión sobre la línea de dis continuidad, pero puesto que ~ son soluc ionesa particulares de ( 19), podemos ajustarlos de manera que: cp l ~2 act> 1 act> 2 a;;-= an en av e (22) con lo cual la s Integrale s de (21) se reducen - al contorno del dominio que ya ha tenido que -- ser discretizado. El tratamiento para este caso particular de fu nci ón discontinua es diferente en cuanto que no se utiliza (12), pero es similar en cuantoque la sust i tuc ión de la integral de dominio se hace igual que en el caso general según (15). 4. EJEMPLOS. Tres ejemplos van a ser considerados, compa rando en el los el valor de la Integral calcula~ do numéricamente a partir de integraciones numé ricas sobre el domin io discretizado con el va-~ lor que se obtiene a partir de las lnte ·g racio-- nes de contorno que se proponen, o bien los resultados finales de la resolución del problema . Para los dos primeros ejemplos que cor re sponden al caso general se han utiliz ado funciones 11-- neales f. (J ,z) J 1 - r(J,z) (23) Las funciones <!l . asociadas se calculan en el -- J apend ice 1 . Para el tercer caso que corresponde a una - función sigular, las funciones <!la se calculan - en el apendice 2. Ej emp 1 o 1. Se trata de resolver la ecuac1on IJ q· <!> = P = cte. que representa, por ejemplo, la ec~~ción de la torsión de un prisma en la formulac1on de Prandtl. Para el contorno, según se indica en la Fig. 4, se han utilizado discretizacion:~ de 16 32 elementos y para aproximar la func1on en - ~~ interior se han uti 1 izado 1, 5 y 9 puntos in ternos. La discretización sobre los lados es re gu 1 a r . 2 1 3 (25) 1" • {1, 9 ( 17) B JL- ---- "' 5 ( 9) v A 17 (3 3) F. 4 D1'scretización de 16 y (32) el emen1 g. . tos. A continuación , en la tabla 1 se obtienen - l os resultados de la integración desde lo s puntos A y B. Los resultados de la integración numé ri ca son los siguientes. Inte g ra c ión de s de A Inte grac ión desde B Tabla 1. Valor de la Inte g ral de dominio. N° de ele me n tos en av 16 32 7126 6888 1 7510 6028 7259 6073 5 9 7632 6288 6349 7383 6330 6392 En la tabla no se in c luyen discretizaciones con 8 el e mentos po rque la s func i on es <ll .,al ser f. 1 i nea l , son cúbicas, y debe h aber Jal menos J 4 puntos para su definición sobre un lad~. Tampoco se incluyen mayor número de puntos Inter-- nos ni de elementos de contorno, porque los resultados no mejoran sensiblemente. Ejemplo 2. En este ejemplo el dominio sigue sien~o un cuadrado que representa una placa. La var1able de campo cp representa la flecha w de la placa. La solicitación de dominio_definida por (24) r~ presenta el valor de M/D = (Mx+My)/(1+v)D en -- cualquier punto de la placa para el caso de car ga transversal uniforme q0• ce ce 2 2 a 1 mn p(x,y) -~ ¡: ¡:; (m +n ) 2 2 2 TI D n=1 m=1 (~) a 2 m Ti x m TI y (24) TI sen sen a a h3 16q E o con D 12 ( 1v3) a -2mn TI mn 2 6 Kg/cm 2 q = 5000 Kg/m E = 2' 1.1 o o a = 2 m. h 1 o cm. ' V = 0,3. La Fig. 5 recoge la discretización y lo s p unto s que sirven para la comparación de resultados. 30 y .... ~ 13 9 a 1 -- B e D [ X v 17 A 5 Fig. 5. Discretización de la placa. (l<j> La Tabla 2 recoge el valor de an en A, co~ pa rando la so luci ón analítica con los resulta-- dos del M.E.C. con 16 elementos de conto~n~ Y 1, 5 y 9 puntos para definir p en el dominiO· CJ <)> +4 Tab la 2. Valor de añ . 10 en A. Analítica 5 9 -+ Puntos 1 1 interno s -2,80 -1 ,52 -2,43 -2,72 lo s va l ores de <)> en el do La Ta bl a 3 re cog e min io con 16 e le mentos de con t orno y 5 ó 9 puntos internos, así como la solución analítica. Tabla 3. Valor de cp. 10+2 en puntos internos . e D E ANALITICA 0,87 1 ,48 · 1 ,69 16 av + 5 v 0,76 1 ,28 1 ,48 16 av + 9 v 0,84 1 ,42 1, 63 La elección de este ejemplo está motivada - por la importancia de esta alternativa en el es tudio de las placas, como continuación de los~ triilbajos iniciados en 17 1. Los resultados quese representan en las Tablas 2 y 3 reflejan una rápida convergencia en l os resultados para un - número de puntos interno s reducido, si tenemos en cuenta la complejidad de la función definida en (24), y para una discretización de contorno que no necesita ninguna modificación respecto a la que se utilizaría para una evaluación numéri ca de la integral de dominio. Ejemplo 3. Es te ejemplo corresponde a una aplicación_del cas o de ca rg a discontinua y se describe en la Fig . 6. B A p = O en V, v2 2 v1 --e 0 ,6 7 1 • 3 3 Fig. 6. Caso de carga discontinua. La s funciones <!l 1 y <!l 2 para este caso se cal cul an en el apendice 2. La Tabla 4 recoge los - resultados de la integral desde los puntos ABC. analítica y numéricamente para 16 y 32 elemen-- tos. Tabla 4. Va lor de la integral de dominio. A B e 452 6051 5683 ANAL 1 T 1 CA 269 6017 5865 405 6042 57 1 1 16 av} k NUMERICA 32 él Vk _L os re s ultado s son excel ente s pa ra discreti zaclones r educ ida s , salvo en el caso de inte-- ~ 31 grar desde la discontinuidad, lo que requiere - una discretización más fina. 5. CONCLUSIONES. El cálculo numérico de la integral de dominio que aparece en la aplicación del M.E.C. a ~ la ecuación de Poisson, por medio de integrales de contorno que se ha propuesto, evita la dis-- cretización del dominio para la evaluación de - dicha integral, lo que permite mantener la gran ventaja del M.E.C. consistente en la discretiza ción del contorno exclusivamente. SÍ es impres~ cindible dar información de los valores de la - f unción en algunos puntos del do minio, pero esto es inherente a la existencia de una función conocida sobre él. A cambio de esta clara ventaja, es preciso real izar algunas integraciones más, tantas como puntos del dominio se utilizan para aproximarla función. A veces, en estos puntos se está in teresado en calcular el valor de la incógnita Ó su/s derivada/s, de la ecuación de campo y -- las integraciones, habría entonces que real izar las de cualquier forma . Lo que sí es preciso es invertir una matriz MxM. Para el caso de funciones discontinuas seha propuesto un método elegante y rápido de eva luar la integra l de dominio. En el caso en que - haya más de dos subdominios o uno de el los sea cerrado , se hace preciso inte grar sobre líneas interiores de los subdomi nios , lo que en el ca - so presentado en el artículo no era necesar io. De cualquier forma, la inte g ral de dominio puede evitarse. L os ejemplos que se muestran, indican que - un número muy reducido de puntos en el dominio es suf iciente para conseguir una precisión equivalente en los resultados a la evolución nu mérica sobre el dominio. Se hace preciso, no -~ obstante, un mayor estudio sobre la naturaleza de las funciones f. a utilizar para una mayorJ precisión y ·para g arantizar la existencia de la matriz inver sa que permite definir lo s coef i- -- cientes a . 6. REFERENCIAS. 1. - Lachat, J.C. y Watson, J.O. "Appli cat ion de la methode des equations integrales au calcul des structures". C.E.T.I.M., n° 25, 1976. 2.- Brebbia, C. y Domínguez J. "Boundary Ele--- ment Methods Versus Finite Elements". Applied Numeri ca l Modell ing. Pentech Pre ss , 1978. 3. - Pari s , F. "El método de los Elementos de -- contorno en la Teoría del Potenc ia l y la -- Elasticidad". Tesis E.T.S. 1.1. de Madrid, 197 9. 4.- Nardini, P. y Brebbia, C.A. "A new Approach to Free Vibration Analysis Using BoundaryElements". Proceedings of the Fourth lnternational Conference on Boundary Element. -- Methods in Engineering. Southampton, 1982. 5.- García Suarez, c.· "Trata m iento de cie rta s - integrales de volumen que aparecen en el Mé todo de los Elementos de Contorno". Anal esde ingeniería, n° 2, 1982. 6.- Ala rcón , E., Martín, A. y Paris, F. "Bounda ry Elements in Potential and Elasticity Theory". Computers and Structures, Vol. 10, pp. 351-362, 1979. 7.- P aris, F. y León, S. "El método de los elementos de contorno aplicado al cálculo de placas apoyadas". Anales de Ingeniería Mecá nica, n° 2, 1982. APENDICE 1. Soluciones part i culares del caso g~ neral 1/ <P . f .. J J Tomando f. = 1 J -r, es preciso resolver () 2 <P d <P 1 --+ - crr -r d r2 r Haciendo el cambio 2¡cjl S obtenemos ar ~+_!_ S = 1 - r dr r cuya solución es de la forma: S= para -1 n r e -r que en 1 f ( 1r) l n r 1 e d r + k r 22 r: 1 r=O, <P 1 00 ' k o. 2 S r r = 2 T y una solución pa r ticu l ar de cjl es por tant o , r2 r3 <P = T 9 APENDICE 11. Solución part icul ar para p discontinua. La Fig. 7 indica la configuración del dominio que se considera. Para que las funcione s <P x cumplan han de s er de l a fo rm a : en V a 32 y D 1 -3D e D 2 p p Fig. 7. Definición de p en el dominio. cjl a 2 a 2 a a a al X + a2 y + a3 xy + a4 X + a5 y + a a + a6 y para que se cumplan las ecuaciones de campo - en los subdominios imponiendo la cont inui dad de <P y su derivada - a con re s pecto a x en 3V siguiente s relaciones:c (x=x1) se obtienen las 1 2 a2 a2 1 2 a4 a4 1 p 1p2 2 al -- 2-+ al 1 (P 2-pl) a3 xl + 2 a3 1 2 pl -p2 2 a5 xl --2-+ a5 Las expresiones de <1> 1 y cjl 2 quedan: pl-p2 1 2 P2 2 2 <1>1 (- 2-+ al) X + (- - a ) y + 2 1 2 + a2 xy + h (p2-pl) + 2 a31 X + 2 2 pl-p2 2 + a4 y + x1 - 2- + a5 2 2 p2 2 2 2 <1> 2 al xl + (- al ) y + a2 xy + 2 2 2 2 + a3 X + a4 y + as una s olución particular de esta s <P a es: p1 -p2 2 2 2 2 p, -p2 <P l + Ly (p2 -p l )x + xl - - -2X 2 + xl 2 cjl 2 P2 2 T Y