scieee AI-readable full text Open interactive document viewer

Análisis del problema de contacto sin fricción con el método de los elementos de contorno

París Carballo, Federico; Garrido, José Antonio

Abstract

Se formula la aplicación del método de los elementos de contorno al problema de contacto, en ausencia de fricción, de dominios bidimensionales - en el campo elástico. Los resultados obtenidos se comparan con otras alterna tivas numéricas de contorno y dominio.

Full text

ANALES DE INGENIERIA MECANICA -AÑO 2 - N.0 1 - 1984 ANALISIS DEL PROBLEMA DE CONTACTO SIN FRICCION CON EL METODO DE LOS ELEMENTOS DE CONTORNO F. Paris* y J.A. Garrido** * E.T.S.I.I. de Sevilla / ** E.T.S.I.I. de Las Palmas Resumen. Se formula la aplicación del método de los elementos de contorno al problema de contacto, en ausencia de fricción, de dominios bidimensionales - en el campo elástico . Los resultados obtenidos se comparan con otras alterna tivas numéricas de contorno y dominio. 1. INTRODUCCION Existe un gran número de problemas, tanto dentro de la Ingeniería Mecánica - como Civil, en los que se producen ióter a cc i one s directas entre los contornos~ de dos o más dominios en estado de defor mación. El estudio de uniones atornilladas, problemas de zunchado, cojinetes, así co mo los pilotes de fricción, la interac-= ción terreno estructura, los problemas - relacionados con el impacto y la modelización de los planos de deslizamiento en la teoría de la plasticidad,son ejemplos que resaltan la importancia del estudio del contacto entre sólidos deformables. No es de extrañar, por tanto, que -- desde que Hertz expusiera los primeros - trabajos 111 hayan aparecido numerosaspublicaciones sobre este tema. 121 repr~ senta una revisión actual del estado de conocimientos . Sin embargo, y dado que - las soluciones analíticas sólo son posibles en casos de geometría simples, para condiciones ideales, teniendo en cual--- quier caso un desarrollo matemático complejo, el conocimiento del problema ha - ido parejo al uso de métodos experimenta les y más recientemente numéricos. En - = concreto el Método de los Elementos Fini tos ha sido usado en la última década en forma profusa para abordar este proble-- ma. Referencias 131 a 1111 representanuna selecta colección de artículos que - usan el Método de los Elementos Finitos. 22 En general el tratamiento del problema de contacto admite · dos enfoques basados en modelar independientemete ambos dorü-- nios utilizando luego unas ecuaciones - de restricción en la zona de contacto, ó bien en usar elementos de enlace en - la zona de contacto cuyas propiedades - reflejen el comportamiento de las super ficies. - El Método de los Elementos de Con-- torno ha sido más recientemente utiliza do ¡¡>or Andersson 112j , 113.1 , 1141 y -= 1151 para abordar el cas o plano con y si :- ! fricción . En el esquema de un método de contorno sólo la primera de las altern~ tivas usadas en el M.E.F. cobra sentido. Ciñendonos al problema de contacto sin fricción, para el caso de que la s~ perficie de contacto sea constante se - trata de un problema directo, es decir las condiciones de contorno pueden esta blecerse de forma definitiva a priori y la carga total puede aplicarse de una - vez. En cambio si la superficie de contacto puede variar con la aplicación de las cargas, las condiciones de contorno no pueden establecerse previamente y e~ ben dos posibilidades. Aplicar directamente toda la carga con una superficie de contacto supuesta y ajustar poste- - - riormente esta superficie en sucesivas aproximaciones hasta que los resultados coincidan cualitativamente con-la supe! ficie de contacto supuesta. Esta . es la alternativa de Andersson ¡¡>ara problemas sin friccion 1121 y 1151. La otra al- ANALES DE INGENIERIA MECANICA AÑO 2 N.0 1 1984 ternativa consiste en dar incrementos de carga que correspondan con variaciones - discretas de la superf i cie de contacto - hasta que se termina de aplicar la car-- ga. Esta técnica es la utilizada por --- ejemplo en lsl. En el presente trabajo se ha seguido e~ta seg~nda alternativa que no requiere n1nguna 1dea previa sobre la solución -- del problema y que permitirá enlazar con la forma de resolver el problema con --- fricción que necesariamente ha de ser in cremen tal. En el apartado 2 se define el probl~ ma. En el apartado 3 se establece el --- planteamiento integral para el problema anteriormente presentado. En el apartado 4 se pr e senta la aproximación numérica con elementos discontinuos. En el 5 se - establece el proceso de resolución, presentando los resultados en el apartado 6 para el problema de Punch elástico sobre cimentación elástica. 2. DEFINICION DEL PROBLEMA Consideremos la situación indicada - en la Fig. 1 en la que dos dominios D y D8 con contornos íl D y íl D se encuent~an A B sometidos a la acción de acciones exte-- riores sobre dicho contorno tA -A tB -B i' ui Y i' ui respectivamente. Bajo la acción de di chas acciones ambos dominios alcanzan -- una situación deformada con una superficie de contacto. Fig. 1: Definición del problema. Ambos contornos pueden particionarse en la siguiente forma: K = A,B donde íl DK L representa la parte del contor no de K que queda fuera de la zona decontacto que asumimos. La h . ó · 1p tes1s de pequeñas deformacio- 23 nes, permite poner: Aunque la Fig. 1 representa la inte racción entre dos dominios planos y laformulación que sigue recoge este caso la extrapolación a más de dos contorno~ o más de una zona de contacto es inme-- diata. Se supone que las tensiones que se alcanzan mantienen al material en su -- comportamiento elástico y se admite la posibilidad de que ambos dominios e s tén caracterizados por propiedades elásti-- cas diferentes EK, v K. Para poder establecer las condiciones ~e contorno del problema, se define un s1sterna local de coordenadas tal corno se indica en la Fig. 2. Fig. 2: Definición del sistema lo-- cal de coordenadas . La dirección 1 es normal a la super ficie y la dirección 2 es su perpendicÜ lar en sentido de las agujas del reloj~ En la zona libre se darán algunas - de estas tres condiciones: - t. t. 1 1 - u. u. 1 1 -r u . 1 i;zoj i,j=1,2 ti t . 1 En la zona de contacto, existen 4 - variables asociadas a un mismo punto M en cada uno de los dominios, K K K K u1 (M), u2 (M), t1 (M), t2 (M) K=A,B pero entre estas ocho variables existen las siguientes relaciones, supuesto que el contacto se realiza sin fricción ANALES DE INGENIERIAMECANICAAÑO 2- N.o l - 1984 A u1 (M) + B u1 (M) = o A B t2 = t2 = o ( 1 ) t~(M) t~ (M) = o con t1 (M) < o Por lo que en un punto de la zona - de contacto existirán cuatro incógnitas. A A A B u1 (M), t1 (M) , u2 (M), u2 (M) 3. FORMULACION DEL PROBLEMA EN ECUACIONES INTEGRALES La formulación del problema elástico en ecuaciones integrales se obtiene a través de la fórmula de Somigliana -- que relaciona tensiones y desplazamientos en el contorno. Asf, para un domi-- nio (Fig. 3) D, con contorno ao dicha - fórmula toma la siguiente expresión, con fuerzas de volumen nulas: C .. (x) u. (x) + JT .. (X , y) u1 . (y) ds (y) l.J J l.J - a o Juij(x,y) tj (y) ds(y) a o (2) X Fig. 3: Configuración del dominio. donde: x,y E ao uJ. (y), tJ. (y) son desplazamientos y tensiones en ao U .. (x,y), l.J e .. (x) l.J T .. (x,y) representan lal.J solución de Kelvin, -- desplazamientos y tensiones en y en la di-- rección j al aplicar - carga unidad en x en - la dirección i. es el término libre cu yo valor depende de la forma local del contar no en x. Si X t ao y la norm a l es c o ntinua en 1 x. c . . (x) = -2 6 . . (x). l.J . l. J 24 Una más detallada obtención de (2) pue de encontrarse en 1161. Como es bien= conocido (2) representa la ecu~ción bá sica para resolver un problema elásti= co mediante el método de los elementos de contorno. Por ello veamos la forma que adopta (2) al aplicarla a los doro~ nios A y B del problema definido en la Fig. 1. f <j (x,y)K u~(y) ds(y) + f T:j (x,y)K u~(y) ds(y) aoL aoc = J ~j(x,y) t~(y) ds(y) + aoL + f Tf..(x,y) t~(y) ds(y) K=A,B (3) l.J J a oc donde T~. representa el núcleo de inte l.J - gración de (2) incluyendo el efecto -- del término libre. Si introducimos en (3) las condi-- ciones de contorno (1) deJ problema, - las dos ecuaciones adoptan la siguiente forma u~ ds J u~ ds J A u1 ds t A. ds + f UA J i 1 a oc A u1 ds t A. ds + J UA J i 1 a oc (4) + f Tiz u~ ds = a oc ( 5) (4) y (5) representan un sistema de -- dos ecuaciones integrales que pueden - aplicarse en puntos del contorno de A y B respectivamente. Por tanto, para - cada punto de A y B que no pert~nec~n a la superficie de contacto y que tl.enen dos incógnitas se dispone de dos - ecuaciones integrales, y para cada pu~ to de la superficie de contacto se di~ pone de cuatro ecuaciones integrales ~ (dos de (4) y dos de (5)) para determ~ nar las cuatro incógnitas que anterio! mente se habfan especificado en el --- apartado 2. 4. APLICACION DEL METODO DE LOS ELEME~ TOS DE CONTORNO. DISCRETIZACION CON ELEMENTOS DISCONTINUOS El Método de lo s E lement os de Co nt or no se ba sa , como es bi en conocid o , ANALES DE INGENIERlAMECANICA AÑO 2 N.0 1 1984 en realizar una aproximación sobre las variables u. y t. que aparecen en la - l. l. ecuación integral (2). Esta aproxima-- ción realizada en la forma: u N .• a. l. l. t = N .• b. l. l. (6) donde ai y bi representan los desplaz~ mientas y tensiones en puntos del contorno fijados como referencia para establecer la variación, permite que la ecuación (2), al aplicarla en diferentes puntos x del contorno,de origen a un sistema lineal de ecuaciones (7) cu yas incógnitas son los valores descono cidos de u y t en el contorno. NE 1 r T~ . (x ,y) ¡: u . (y)> al\ d~ - K=1 ) l. J J aok -J uij(x,y) t. (y)> aD dsk 1 o aok J k (7) donde NE representa el número de elemen tos (J DK utilizados para representar el - contorno a D. La aproximación que se va a seguir en este estudio es lineal no conforme. El contorno aD se aproxima por elemen- - tos lineales aDK, sobre los cuales se - supone una variación lineal de la fun-~ ción a partir de los valores de dicha - función en dos puntos interiores del -- elemento, lo que no garantiza la continuidad interelemental, al igual que sucede con los elementos constantes. La - formulación detallada de este tipo de - e lementos puede encontrarse en 1171. Es ta aproximación permite expresar l as va riables de (7) u. y t. en aDK' a partir l. l. de los valores de las tensiones y des-- plazamientos en dos puntos interiores - del elemento en un sistema local de --- co ordenadas (Fig. 4). ~) aD ~( ~ ) k !! k ( 8) k . 2 (a,b) ~nL k (8) !) Cl D ~( ~ ) k !!k( 8 ) tk . 2 (a,b) k -nL ,- = _, ·= =1 Fig . 4: Defi ni c i ón de l e l ement o uti l i zad o e n l a ap ro xim ac ión. - 25 La introducción de (8) en (7) permi te obtener un sistema lineal de ecuacio nes que puede representarse en la for= m a H u (4NEx4NE) (4NEx1) G (4NEx4NE) t (4NEx1) 1 (9) donde ~ y Q son matrices cuyos coefi--- cientes dependen de los núcleos de inte gración de la Ecuación de Somigliana -= T .. y U .. , y de las matrices N, Q y L - l.J l.J - - - introducidas en (8). u y t representan los valores de los desplazamientos y -- las tensiones en el sistema local de -- coordenadas definido para cada elemento. Para el problema de contacto que se está tratando se puede obtener un siste roa de ecuaciones similar al (9) para ca da uno de los dominios implicados en el problema, estando estos sistemas rela-- cionados a través de (1). Si tomamos como referencia las ecua cienes integrales (4) y (5) que ya - contemplaban las condiciones ( 1 ), e intro ducimos en ellas la aproximación, no re sulta dificil establecer una configura= ción genérica del sistemas de ecuacio-- nes que aparece. En orden a la claridad de dicha representación, y sin que se - establezca restricción alguna sobre los problemas a tratar, supongamos que los contornos exteriores a la zona de con-- tacto de ambos dominios (J DL se pueden - particionar en una zona a Du donde están prescritos los desplazamientos y otra - (J Dt donde están prescritas las tensio-- nes. a D~ + a D~ a DK + a DK a K u t + De K= A,B (10) Llamamos NK el número de e l ementos de cada dom ini o que estará particiona- K K K . K K do en N , Nt y Nc asoc1.ados a ao , aD K u u t y aoc respectivamente, cumpliéndose obviamente: NK NK u + NK t + NK e NA e NB e ( 11 ) Con estas premisas, ( 3) para cada -- contorno K darfa lugar a las siguientes NK ecuaciones lineales: K ~L - K ~L K · ~1C K ~ 2 C ANALES DE INGENIERIA MECANICAAÑO 2 - N.0 1- 1984 ~~~ 1 GK -u K ~1 1 ~21 -K !_L tK -L K !c1 K !c2 ( 1 2) El significado de cada subrnatriz está - en función de las definiciones (10) y - (11), siendo sus dimensiones: HK + (4NK, 4NK) -t t HK + -u (4NK, 4NK) u K É!c1 .... (4NK, 2NK) e K É!c2 .... (4NK, 2NK) e Finalmente las ecuaciones integra-- les (4) y (5) en las que se han aplicado las condiciones de contacto producirán , siguiendo la representación (12) el siguiente sistema: HA -t o -GA -u o A É!c1 o HB o -GB A -t -u -!!e, )-HA -A A -u ~L + ~t l-HB . -B GB -u ~L + -t A É!c2 o o B É!c2 -A ) !.L -B !_L A -~1 B -~1 A ~L B ~L A ~1C A ~2C B ~2C A .i1c ( 1 3) (13) representa un sistema de 4(NA+NB) - ecuaciones con 4NA + 4NB + 4N~ + 4N~ + A A Bu Au + 2Nc + 2Nc + 2Nc + 2Nc incógnitas por lo que teniendo en cuenta (1) se tratade un sistema determinado cuya resolu--- ción permitirá, en la forma que se expli ca en el apartado siguiente, calcular -= las incógnitas (tensiones y/ó desplaza-- mient os ) en los contornos de ambos dominios. 5. EL PROCESO DE RESOLUCION En el caso en que la superficie de - 26 contacto final del problema sea conocida, el sistema de ecuaciones (13) puede resolverse directamente aplicando toda la solicitación exterior tK y uK ya que durante la actuación de estas acciones exteriores, las condiciones de contorno de ambos dominios no varfan. En el caso en que la superficie de contacto final sea desconocida puesto - que dependa de la carga aplicada (con-- tacto entre dos esferas, por ejemplo), el problema no puede resolverse de una sola vez, existiendo dos técnicas para abordarlo, al tratarse sólo de una nolinealidad geométrica, es decir, cambio de la longitud total de la zona de contacto y no cambio cualitativo de las -- condiciones dentro de dicha zona. Pues- to que en el problema con fricción además de la primera existe también esta - segunda no-linealidad, parece razonable de las dos alternativas, elegir aquella que permite también el tratamiento de - los problemas de contacto con fricción, es decir, la técnica incremental. Por - ello, en el caso de que la superficie - de contacto sea susceptible de variar - al evolucionar la solicitación, todas - las variables del problema, tanto las - incógnitas u y t corno los datos u y t - que aparecen en-(13) deben interpretarse en sentido incremental. El valor del incremento está rela-- cionado con la variación de solicita--- ción que provoca el contacto entre dos puntos. Puesto que en la zona prevista de contacto siempre mayor que la real, se ha modelado igual A y B (misma longi tud de los elementos y misma situaciónde los nodos en ellas) los pares de nodos candidatos a contactar están claramente definidos. Ahora bien, aunque cualquier pequeño incremento de carga permite conocer las direcciones de los desplazamientos en dos nodos candidatos, pA, pB, ello - no implica la definición inmediata del incremento para que se produzca el contacto ya que ambos desplazamientos no - necesariamente provocan que los dos nodos lleguen para el mismo valor de la - carga al punto p. (Fig. 5). Por ello es preciso establecer un - criterio de contacto entre los dos puntos en función de las proyecciones de - los vectores desplazamientos u(pA) y -- u(pB) sobre una dirección úni~a, tornanaose un valor intermedio entre las normales a los elementos a que pertenecen A B p y p . La abertura in ici al A entre dos pun- ANALES DE INGENIERIA MECANICA AÑO 2 N.o 1 1984 Fig. 5: Definición de la abertura - entre dos puntos candidatos al contacto. tos que van a contactar se establece -- pues corno: 1:. = r AB ( p) • V AB ( p) ( 1 5) El incremento de la carga, viene pues condicionado porque las proyeccio~ nes de los desplazamientos totales -- ( B) · A AB u p y u(p ) sobre la direcci~n v -- sea igual a 1:. . Si en los nodos candidatos a entrar en la zona de contacto se utiliza como sistema local de coordenadas, el defin~ do por v AB y su correspondiente normal, lo que puede obtenerse a través de la - matriz L( 8 ) en (8), la solución del sis terna (13) permite directamente estable= cer la comparación para determinar el - menor valor del incremento. Asimismo - este sistema local de coordenadas co~ún a 717rnentos en contacto, es el que se - ut1l1za para establecer las condiciones de contorno definidas en (1). El uso de e lementos discontinuos exige un cuidado extremo al considerar ~uevas zonas de contacto, puesto que la lntroducción solamente de los dos prime ros nodos del elemento en contacto pue= de_conducir a irncornpatibilidades geornétrlcas corno se indica en la Fig. 6. Por ello es preciso recurrir a un - Valor promediado de los incrementos de carga necesarios para unir los nodos p ~ los q_ de A y B, o bien, es preciso in 1 roduc1r restricciones cinernáticas en= dos desplazamientos de los cuatro nodos e lo s elementos qu e van a contactar pa ra que dicho contacto se produzca de -= una sola ve z. 27 o • Fig. 6: Incompatibilidades geornétri cas en elementos disconti-= nuos. 6. RESULTADOS Se presentan los resultados obtenidos en la compresión de un dominio rectangular sobre una cimentación elástica también rectangular, corno se indica en la Fig. 7. 8 100 ' 2 1,2 Nw/rrm EA ~ Nw/mm 2 VA VB 0,35 llllllllllllllf o A t 100 J r B 160 UJD Fig. 7: Configuración del ejemplo. La simetría del problema permite es tudiar solamente la configuración repre sentada en la Fig. 8 con la discretiza= ción que en ella se indica. Los resultados se presentan en las Figs. 9 y 10. La primera representa el valor de la tensión normal en la zona de contacto referido a la carga unitaria. Se compara con los resultados obtenidos por Andersson 1141 con elementos parabó licos. La Fig. 10 representa los desli= zarnientos relativos que se producen en la zona de contacto, relativos al desli zarniento máximo. Los resultados son coro parados co n los obtenidos por Fr e driks= son 131 con el Método de los Ele m entos Finitos, utilizando elementos triangul~ ANALES DE INGENIERIA MECANICAAÑO 2 - N. 0 1- 1984 Po tf 20 elerrentos J3 30 elementos Nc 5 elerrentos Fig. 8: Discretización utilizada. res con variación lineal de desplaza--- miento con 496 elementos, 712 grados de libertad y 27 nodos de contacto. 1/"'""\ '· . ~ elementos discontinuos ~ e lement os parabólicos (l\ndersson) dD e Fig. 9: Tensiones normales en la zo na de contacto. Fig. 10 : VE) (M.E.F) = 1.04 . !1'aX -2 . 10 ntr -2 (M.E .C .) = 1,00.10 mn a Dc Deslizamientos relativos - en la zona de contacto. 7 . CONCLUSIONES Se ha presentado una f o rmulación en ecuac i ones integrales que permite la -- aplicación del Método de los Elementos de Contorno para el problema de contac - to de dominios bidimensionales sin frie ción. 28 La formulación no incluye las varia bles intermedias utilizadas por Anders= son !121, no resultando necesario establecer a priori la ubicación de la zona común de contacto. La comparación con Elementos Fini-- tos resulta para este problema,ventajosa para la formulación presentada. Por un lado, conceptualmente se trata de un problema en que todas las variables que interesan están en el contorno. En el - M.E.F. es preciso utilizar · elementos su periores o necesariamente se producir~n errores al aplicar las condiciones de - contorno en tensiones (Fredriksson 131). Se trata de un problema que hay que resolver incrementalmente, una parte que puede ser considerada costosa computa-- cionalmente en el M. E.C, como es el cál culo de desplazamientos y tensiones enpuntos internos, si es que se está inte resada en muchos puntos, no tiene paraeste problema una importancia destaca-- ble en el tiempo total de cálculo. Los elementos discontinuos presenta dos en este trabajo para el problema de contacto han dado excelentes resultados lo que unido a la simplicidad para apli car condiciones de contorno, para evi-= tar el problema de las esquinas, hace - pensar en su ventajosa aplicación para problemas tridimensionales. Los mayores esfuerzos para dismi--- nuir los tiempos de computación deben - encaminarse a mejorar la resolución del sistema de ecuaciones siguiendo lineas de actuación como las propuestas por -- !!,.a:rks 111 11 y Bettess 1181 . 8. REFERENCIAS 1. Hertz. "Gesamelte Werke, Bdo. 1". Leipzig, 1895. 2. Gladwell, G.M.L. "Contact· Problems in the classical Theory of Elasticity". Sijthoff and Noordhoff (1980). 3. Fredriksson, B. "Finite Element solution of surface nonlinearities in structural mechanics with special emphasis to contact and fracture mechanics problems". Computers and Structures, Vol. 6, pp. 281-290, 1976. 4. Gaertner, R. "Investigation of plane elastic contact allowing for friction". Computer and Struc ture, Vol. 7, pp. 59-63, 1977. - 5. Okamoto, N., y Nabazawa, M. "Finite element incremental cont act ana lysi s with various frictional conditions". Int . J. Num. M. Eng., Vol. 14, pp . 337 - 357 , 1979. ANALES DE INGENIERIAMECANICA AÑO 2 N.o 1 1984 6. Herrmann, L.R. "Finite Elements analysis of contact problems", ASCE, pp. 1043-1057, octubre, 1978. 7. Nguyen Dang Hung y Gery de Saxce "Frictionless contact of elastic bodies by Finite Element Method and Mathematical programming technique" Computers and Structure, Vol. 11, pp. 55-67, 1980. 8. Mahmoud, F.F., Salamon, N.J., Marks, W.R. "A Direct automated procedure for frictionless contact problems". Int. J. Num. M. Eng., Vol. 8, pp. 245- 257 (1982). 9. Haber, R.B., Abel, J.F. "Contact-Slip Analysis using Mixed Displacement". Journal of Eng. Mech. ASCE. Vol. 109, NQ 2. Abril, 1983. 10. Torstenfelt, B. "Contact Problems with friction in general purpase Finite Element computer Programs". Computers and Structures, Vol. 16, NQ 1-4, pp. 487-493 (1983). 11. Marks, W.R. y Salamon, N.J. "A modified Conjugate Gradient Method for frictionless Contact Problems". ASME. Vol. 105. April, 1983. 12. Andersson, T., Fredriksson, B., Allan Persson, B.G. · "The Boundary Elcment Me thod applied to two dimensional contact Problems". New Developments in Boundary Element Methods Ed. by C.A. Brebbia CML P ublications. 1980. 13. Andersson, T. "The .Boundary Elment Method applied to two-dimensio- nal contact Problem with Friction". Boundary Element Method. Ed. by C.A. Brebbia. Springer-Verlag, 1981. 14. Andersson, T. "Th e Second generation Boundary Element Contact Program". Boundary Element Methods in Engineering Springer-Verlag, 1982. 15. Andersson, T. y Allan Persson B G "The Boundary Element Method applied .. to two dimensional Contact problems" . Progress in Boundary Element Method. Ed. by C.A. Brebbia. Vol. 2. Pentech Press. 1983 . 16. Paris, F . "El Método de los Element~s de Contorno en la Teorfa del Potenc1al Y la Elasticidad" Tesis Doctoral. Madrid: 1979. 17. Paris F G ·d "El , :Y arr1 o, J.A. ementos d1scontinuos lineales en el a~álisis elástico bidimensional con el ~ntodo de los Elementos de Contorno". 0 ales de Ingenier ía M ecánica . Año 1, N_ l, pp. 20 - 25. 1983. 29 18. Bettess, J.A. "Economical solution technique for boundary integral matrices". Int. J. Num. Mat. Eng. Vol. 19, pp. 1073-1077. 1983.